diff --git a/bin/mrdna b/bin/mrdna
index 440e45073556469bcd61e04c43067e14f15afd5b..2430aadeb170b2e17ffeab2c646e60f2cc51cd99 100755
--- a/bin/mrdna
+++ b/bin/mrdna
@@ -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
diff --git a/mrdna/model/nbPot.py b/mrdna/model/nbPot.py
index 305d34bcc14cad0dae86659c83ff573097fbe84a..c26c3b07f116bc67121642c0a7dbf5c1d28df374 100644
--- a/mrdna/model/nbPot.py
+++ b/mrdna/model/nbPot.py
@@ -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")
diff --git a/mrdna/readers/__init__.py b/mrdna/readers/__init__.py
index 1d6a7893752fcb01ed1c664dd8d2e88939e8daa8..44e5a32cdebff722ac2bf604378220756147017e 100644
--- a/mrdna/readers/__init__.py
+++ b/mrdna/readers/__init__.py
@@ -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)
diff --git a/mrdna/readers/cadnano_segments.py b/mrdna/readers/cadnano_segments.py
index 689b138a88c81263da59501c711acabd805a5037..612b5f95406a1bf6466b730f9e3d9cfdb85d8b5c 100644
--- a/mrdna/readers/cadnano_segments.py
+++ b/mrdna/readers/cadnano_segments.py
@@ -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:
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index cfacf7a29d61fa729c1fa93ccb82da2d737f6fa6..ad49d1220334e18f4596415d68b0bd0116daee5b 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -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