Commit 47325593 authored by cmaffeo2's avatar cmaffeo2
Browse files

Added ability to adjust screening length

parent 80f9a1c5
......@@ -23,6 +23,9 @@ parser.add_argument('-d','--directory', type=str, default=None,
help='Directory for simulation; does not need to exist yet')
parser.add_argument('-g','--gpu', type=int, default=0,
help='GPU for simulation; check nvidia-smi for availability')
parser.add_argument('--debye-length', type=float, default=None,
help='Adjust the electrostatic repulsion matching osmotic pressure data from Rau and Parsegian under 25 mM MgCl2 condition with a Debye-Hueckel correction from the default 11.1 Angstroms.')
parser.add_argument('--output-period', type=float, default=1e4,
help='Simulation steps between DCD frames')
# parser.add_argument('--minimization-steps', type=float, default=0,
......@@ -39,6 +42,7 @@ parser.add_argument('--oxdna-steps', type=float, default=None,
parser.add_argument('--oxdna-output-period', type=float, default=1e4,
help='Simulation steps between oxDNA configuration and energy output')
parser.add_argument('--backbone-scale', type=float, default=1.0,
help='Factor to scale DNA backbone in atomic model; try 0.25 to avoid clashes for atomistic simulations')
parser.add_argument('--debug', action='store_true',
......@@ -78,7 +82,7 @@ def main():
else:
raise Exception("Unrecognized input file '{}'".format(infile))
model = read_model( str(infile) )
model = read_model( str(infile), debye_length=args.debye_length )
if args.output_prefix is not None:
prefix = args.output_prefix
......
......@@ -23,6 +23,8 @@ class AbstractNbDnaScheme(NonbondedScheme):
def __init__(self, *args,**kwargs):
NonbondedScheme.__init__(self,*args,**kwargs)
self.debye_length = self.debye_length0 # set by child class
x,y = self.load_pmf()
self.target_x = x
self.target_y = y
......@@ -73,6 +75,15 @@ class AbstractNbDnaScheme(NonbondedScheme):
bounds_error = False, fill_value="extrapolate",
assume_sorted=True)(x).T
u[x<13] = u[x<13] + 0.5 * (x[x<13]-13)**2
## Debye correction
if self.debye_length != self.debye_length0:
## units "e**2/(4 * pi * 80 epsilon0 AA)" kcal_mol
A = 4.1507964
du = (A/x) * (np.exp(-x/self.debye_length) - np.exp(-x/self.debye_length0))
du[x < 10] = du[x>=10][0]
u = u + du
return u-u[-1]
def get_bead_distributions(self, interhelical_distances):
......@@ -320,6 +331,11 @@ def _add_parameter_to_cache_dict( args ):
class nbDnaScheme100Mg(AbstractNbDnaScheme):
def __init__(self,*args,**kwargs):
## units "sqrt( 80 epsilon0 295 k K / ((4*100 mM + 200 mM) e**2/particle) )" AA
self.debye_length0 = 5.5771297
AbstractNbDnaScheme.__init__(self,*args,**kwargs)
def get_package_cache_file(self):
return get_resource_path("nbPot.rau_refined_100Mg.dat")
......@@ -369,6 +385,12 @@ class nbDnaScheme100Mg(AbstractNbDnaScheme):
class nbDnaScheme20Mg(AbstractNbDnaScheme):
def __init__(self,*args,**kwargs):
## units "sqrt( 80 epsilon0 295 k K / ((4*25 mM + 50 mM) e**2/particle) )" nm
## 25 mM MgCl2, sorry about name of this class/files
self.debye_length0 = 11.154259
AbstractNbDnaScheme.__init__(self,*args,**kwargs)
def get_package_cache_file(self):
return get_resource_path("nbPot.rau_refined_20Mg_200Na.dat")
......
......@@ -15,7 +15,7 @@ from .segmentmodel_from_pdb import SegmentModelFromPdb
def read_cadnano(json_file, **model_parameters):
data = read_json_file(json_file)
return model_from_cadnano_json(data)
return model_from_cadnano_json(data, **model_parameters)
def read_vhelix(maya_file, **model_parameters):
data = parse_maya_file(maya_file)
......
......@@ -180,9 +180,7 @@ def combineCompactRegionLists(loHi1,loHi2,intersect=False):
class cadnano_part(SegmentModel):
def __init__(self, part,
max_basepairs_per_bead = 7,
max_nucleotides_per_bead = 5,
local_twist = False
**kwargs
):
self.part = part
try:
......@@ -198,10 +196,7 @@ class cadnano_part(SegmentModel):
self._add_crossovers()
self._add_prime_ends()
SegmentModel.__init__(self, self.segments,
local_twist = local_twist,
max_basepairs_per_bead = max_basepairs_per_bead,
max_nucleotides_per_bead = max_nucleotides_per_bead,
dimensions=(5000,5000,5000))
**kwargs)
def _cadnano_part_to_segments(self,part):
try:
......@@ -625,13 +620,11 @@ def decode_cadnano_part(json_data):
def package_archive( name, directory ):
...
def read_model(json_data, sequence=None):
def read_model(json_data, sequence=None, **kwargs):
""" Read in data """
part = decode_cadnano_part(json_data)
model = cadnano_part(part,
max_basepairs_per_bead = 5,
max_nucleotides_per_bead = 5,
local_twist=False)
**kwargs)
# TODO
# try:
......
......@@ -1536,10 +1536,12 @@ class SegmentModel(ArbdModel):
def __init__(self, segments=[], local_twist=True, escapable_twist=True,
max_basepairs_per_bead=7,
max_nucleotides_per_bead=4,
dimensions=(1000,1000,1000), temperature=291,
dimensions=(5000,5000,5000), temperature=291,
timestep=50e-6, cutoff=50,
decompPeriod=10000, pairlistDistance=None,
nonbondedResolution=0,DEBUG=0):
nonbondedResolution=0,DEBUG=0,
debye_length = None
):
self.DEBUG = DEBUG
if DEBUG > 0: print("Building ARBD Model")
ArbdModel.__init__(self,segments,
......@@ -1559,6 +1561,9 @@ class SegmentModel(ArbdModel):
self.grid_potentials = []
self._generate_bead_model( max_basepairs_per_bead, max_nucleotides_per_bead, local_twist, escapable_twist)
if debye_length is not None:
nbDnaScheme.debye_length = debye_length
self.useNonbondedScheme( nbDnaScheme )
self.useTclForces = False
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment