From 1a6e08bcc4bebd8f3752e905d73c1c72bdbd812a Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Tue, 28 Mar 2023 14:47:54 -0500
Subject: [PATCH] Committing things...
---
arbdmodel/__init__.py | 73 ++++++++++++++++------
arbdmodel/abstract_polymer.py | 2 +-
arbdmodel/coords.py | 4 +-
arbdmodel/hps_polymer_model.py | 4 +-
arbdmodel/kh_polymer_model.py | 68 +++++++++++++++-----
arbdmodel/kh_polymer_model_pair_epsilon.py | 8 +++
arbdmodel/onck_polymer_model.py | 4 +-
7 files changed, 121 insertions(+), 42 deletions(-)
diff --git a/arbdmodel/__init__.py b/arbdmodel/__init__.py
index db37aac..2ddc874 100644
--- a/arbdmodel/__init__.py
+++ b/arbdmodel/__init__.py
@@ -388,10 +388,10 @@ class ParticleType():
excludedAttributes = ("idx","type_",
"position",
"children",
- "parent", "excludedAttributes","rigid"
+ "parent", "excludedAttributes",
)
- def __init__(self, name, charge=0, parent=None, **kargs):
+ def __init__(self, name, charge=0, parent=None, rigid=False, rigid_body_potential=None, **kargs):
""" Parent type is used to fall back on for nonbonded interactions if this type is not specifically referenced """
if parent is not None:
@@ -402,7 +402,9 @@ class ParticleType():
self.name = name
self.charge = charge
self.parent = parent
-
+ self.rigid = rigid
+ self.rigid_body_potential = rigid_body_potential
+
for key in ParticleType.excludedAttributes:
assert( key not in kargs )
@@ -810,7 +812,7 @@ class ArbdModel(PdbModel):
def __init__(self, children, origin=None, dimensions=(1000,1000,1000), temperature=291, timestep=50e-6,
particle_integrator = 'Brown',
cutoff=50, decomp_period=1000, pairlist_distance=None, nonbonded_resolution=0.1,
- remove_duplicate_bonded_terms=True, extra_bd_file_lines=""):
+ remove_duplicate_bonded_terms=True, dummy_types=tuple(), extra_bd_file_lines=""):
PdbModel.__init__(self, children, dimensions, remove_duplicate_bonded_terms)
self.origin = origin
@@ -832,7 +834,8 @@ class ArbdModel(PdbModel):
self.numParticles = 0
self.particles = []
self.type_counts = None
-
+ self.dummy_types = dummy_types
+
self.nbSchemes = []
self._nbParamFiles = [] # This could be made more robust
self.nbResolution = 0.1
@@ -869,7 +872,8 @@ class ArbdModel(PdbModel):
return s
elif typeA.is_same_type(A) and typeB.is_same_type(B):
return s
- raise Exception("No nonbonded scheme found for %s and %s" % (typeA.name, typeB.name))
+ # raise Exception("No nonbonded scheme found for %s and %s" % (typeA.name, typeB.name))
+ # print("WARNING: No nonbonded scheme found for %s and %s" % (typeA.name, typeB.name))
def _countParticleTypes(self):
## TODO: check for modifications to particle that require
@@ -881,8 +885,11 @@ class ArbdModel(PdbModel):
type_counts[t] += 1
else:
type_counts[t] = 1
+ for t in self.dummy_types:
+ if t not in type_counts:
+ type_counts[t] = 0
self.type_counts = type_counts
-
+
def _updateParticleOrder(self):
## Create ordered list
self.particles = [p for p in self]
@@ -906,7 +913,7 @@ class ArbdModel(PdbModel):
if typeA != typeB:
self.nbSchemes.append( (nonbonded_scheme, typeB, typeA) )
- def simulate(self, output_name, output_directory='output', num_steps=100000000, timestep=None, gpu=0, output_period=1e4, arbd=None, directory='.', restart_file=None, replicas=1, write_pqr=False, log_file=None, dry_run = False):
+ def simulate(self, output_name, output_directory='output', num_steps=100000000, timestep=None, gpu=0, output_period=1e4, arbd=None, directory='.', restart_file=None, replicas=1, random_seed=None, write_pqr=False, log_file=None, dry_run = False):
assert(type(gpu) is int)
num_steps = int(num_steps)
@@ -955,7 +962,7 @@ class ArbdModel(PdbModel):
self.writePdb( output_name + ".pdb" )
if write_pqr: self.write_pqr( output_name + ".pqr" )
self.writePsf( output_name + ".psf" )
- self.writeArbdFiles( output_name, numSteps=num_steps, outputPeriod=output_period, restart_file=restart_file )
+ self.writeArbdFiles( output_name, numSteps=num_steps, outputPeriod=output_period, restart_file=restart_file, random_seed=random_seed )
# os.sync()
## http://stackoverflow.com/questions/18421757/live-output-from-subprocess-command
@@ -994,7 +1001,7 @@ class ArbdModel(PdbModel):
# Methods for printing model #
# -------------------------- #
- def writeArbdFiles(self, prefix, numSteps=100000000, outputPeriod=10000, restart_file=None):
+ def writeArbdFiles(self, prefix, numSteps=100000000, outputPeriod=10000, restart_file=None, random_seed=None):
## TODO: save and reference directories and prefixes using member data
d = self.potential_directory = "potentials"
if not os.path.exists(d):
@@ -1019,7 +1026,7 @@ class ArbdModel(PdbModel):
self._writeArbdProductPotentialFile()
self._writeArbdGroupSitesFile()
self._writeArbdPotentialFiles( prefix, directory = d )
- self._writeArbdConf( prefix, numSteps=numSteps, outputPeriod=outputPeriod, restart_file=restart_file )
+ self._writeArbdConf( prefix, numSteps=numSteps, outputPeriod=outputPeriod, restart_file=restart_file, random_seed=random_seed )
# def _writeArbdCoordFile(self, filename):
# with open(filename,'w') as fh:
@@ -1070,17 +1077,17 @@ class ArbdModel(PdbModel):
fh.write(" ".join(str(p.idx) for p in c) + "\n")
- def _writeArbdConf(self, prefix, randomSeed=None, numSteps=100000000, outputPeriod=10000, restart_file=None):
+ def _writeArbdConf(self, prefix, random_seed=None, numSteps=100000000, outputPeriod=10000, restart_file=None):
## TODO: raise exception if _writeArbdPotentialFiles has not been called
filename = "%s.bd" % prefix
## Prepare a dictionary to fill in placeholders in the configuration file
params = self.__dict__.copy() # get parameters from System object
- if randomSeed is None:
+ if random_seed is None:
params['randomSeed'] = ""
else:
- params['randomSeed'] = "seed %s" % randomSeed
+ params['randomSeed'] = "seed %s" % random_seed
params['numSteps'] = int(numSteps)
# params['coordinateFile'] = "%s.coord.txt" % prefix
@@ -1174,10 +1181,15 @@ transDamping {g} {g} {g}
""".format(mass=pt.mass, g=gamma)
else:
raise ValueError("Unrecognized particle integrator '{}'".format(self.particle_integrator))
+ if pt.rigid_body_potential is not None:
+ particleParams['rigid_potential'] = "rigidBodyPotential {}".format(pt.rigid_body_potential)
+ else:
+ particleParams['rigid_potential'] = ''
fh.write("""
particle {name}
num {num}
{dynamics}
+{rigid_potential}
""".format(**particleParams))
if 'grid' in particleParams:
grids = []
@@ -1192,9 +1204,26 @@ num {num}
fh.write("gridFile {}\n".format(" ".join(grids)))
fh.write("gridFileScale {}\n".format(" ".join(scales)))
-
else:
fh.write("gridFile {}/null.dx\n".format(self.potential_directory))
+ if 'rb_grid' in particleParams:
+ grids = []
+ scales = []
+ for entry in pt.rb_grid:
+ try:
+ g, s = entry
+ except:
+ g = entry
+ s = 1
+ ## TODO, use Path.relative_to?
+ try:
+ grids.append(str( g.relative_to(os.getcwd()) ))
+ except:
+ grids.append(str(g))
+ scales.append(str(s))
+ fh.write("rigidBodyPotential {}\n".format(" ".join(grids)))
+ if any([s!=1 for s in scales]):
+ raise NotImplementedError
## Write coordinates and interactions
fh.write("""
@@ -1207,8 +1236,9 @@ tabulatedPotential 1
## The i@j@file syntax means particle type i will have NB interactions with particle type j using the potential in file
""".format(**params))
for pair,f in zip(self._particleTypePairIter(), self._nbParamFiles):
- i,j,t1,t2 = pair
- fh.write("tabulatedFile %d@%d@%s\n" % (i,j,f))
+ if f is not None:
+ i,j,t1,t2 = pair
+ fh.write("tabulatedFile %d@%d@%s\n" % (i,j,f))
## Bonded interactions
restraints = self.get_restraints()
@@ -1305,7 +1335,14 @@ component "data" value 3
scheme = self._getNbScheme(t1,t2)
scheme.write_file(f, t1, t2, rMax = self.cutoff)
self._nbParamFiles.append(f)
-
+ # try:
+ # scheme = self._getNbScheme(t1,t2)
+ # scheme.write_file(f, t1, t2, rMax = self.cutoff)
+ # self._nbParamFiles.append(f)
+ # except:
+ # self._nbParamFiles.append(None)
+
+
def _getNonbondedPotential(self,x,a,b):
return a*(np.exp(-x/b))
diff --git a/arbdmodel/abstract_polymer.py b/arbdmodel/abstract_polymer.py
index 732a956..9c80f61 100644
--- a/arbdmodel/abstract_polymer.py
+++ b/arbdmodel/abstract_polymer.py
@@ -195,7 +195,7 @@ class PolymerSection(ConnectableElement):
# self.sequence = None
def __repr__(self):
- return "<{} {}[{:d}]>".format( type(self), self.name, self.num_monomers )
+ return "<{} {}[{:d}]>".format( type(self), self.segname, self.num_monomers )
def set_splines(self, contours, coords):
tck, u = interpolate.splprep( coords.T, u=contours, s=0, k=1)
diff --git a/arbdmodel/coords.py b/arbdmodel/coords.py
index 6c9515f..52b8f89 100644
--- a/arbdmodel/coords.py
+++ b/arbdmodel/coords.py
@@ -29,9 +29,9 @@ def minimizeRmsd(coordsB, coordsA, weights=None, maxIter=100):
def minimizeRmsd(coordsB, coordsA, weights=None):
- q,comA,comB = _minimizeRmsd(coordsB, coordsA, weights)
+ q,comB,comA = _minimizeRmsd(coordsB, coordsA, weights)
assert( np.all(np.isreal( q )) )
- return quaternion_to_matrix(q),comA,comB
+ return quaternion_to_matrix(q),comB,comA
## http://onlinelibrary.wiley.com/doi/10.1002/jcc.21439/full
diff --git a/arbdmodel/hps_polymer_model.py b/arbdmodel/hps_polymer_model.py
index f56d43b..1d246b9 100644
--- a/arbdmodel/hps_polymer_model.py
+++ b/arbdmodel/hps_polymer_model.py
@@ -256,8 +256,8 @@ class HpsModel(ArbdModel):
kwargs['timestep'] = 10e-6
kwargs['cutoff'] = max(4*debye_length,20)
- if 'decompPeriod' not in kwargs:
- kwargs['decompPeriod'] = 1000
+ if 'decomp_period' not in kwargs:
+ kwargs['decomp_period'] = 1000
""" Assign sequences """
if sequences is None:
diff --git a/arbdmodel/kh_polymer_model.py b/arbdmodel/kh_polymer_model.py
index 4d732a0..3d56f2b 100644
--- a/arbdmodel/kh_polymer_model.py
+++ b/arbdmodel/kh_polymer_model.py
@@ -114,9 +114,15 @@ _types = dict(
sigma = 5.86,
)
)
-for k,t in _types.items():
- t.resname = t.name
+for k,t in list(_types.items()):
+ t.resname = t.name
+ t.is_idp = False
+
+ ## Add types for IDPs
+ _types[k+'IDP'] = ParticleType(t.name+'IDP', mass=t.mass, charge=t.charge, sigma=t.sigma, is_idp=True, resname=t.resname)
+
+
class KhNonbonded(NonbondedScheme):
def __init__(self, debye_length=10, resolution=0.1, rMin=0):
NonbondedScheme.__init__(self, typesA=None, typesB=None, resolution=resolution, rMin=rMin)
@@ -134,13 +140,24 @@ class KhNonbonded(NonbondedScheme):
u_elec = (A*q1*q2/D)*np.exp(-r/ld) / r
""" KH scale model """
- # alpha = 0.228
- # epsilon0 = -1.0
- # e_mj = epsilon_mj[(typeA.name,typeB.name)]
- # epsilon = alpha * np.abs( e_mj - epsilon0 )
- epsilon = epsilon_mj[(typeA.name,typeB.name)]
- lambda_ = -1 if epsilon > 0 else 1
- epsilon = np.abs(epsilon)
+ A_is_idp = B_is_idp = False
+ try:
+ A_is_idp = A.is_idp
+ except:
+ pass
+ try:
+ B_is_idp = B.is_idp
+ except:
+ pass
+
+ _idp_scale = (int(A_is_idp)+int(A_is_idp)) * 0.5
+ alpha = 0.159 + _idp_scale * (0.228 - 0.159)
+ epsilon0 = -1.0 + _idp_scale * (-1.36 + 1.0)
+
+ e_mj = epsilon_mj[(typeA.name,typeB.name)]
+ epsilon = alpha * np.abs( e_mj - epsilon0 )
+ # epsilon = epsilon_mj[(typeA.name,typeB.name)]
+ lambda_ = -1 if epsilon0 > e_mj else 1
sigma = 0.5 * (typeA.sigma + typeB.sigma)
@@ -193,6 +210,16 @@ class KhBeadsFromPolymer(Group):
if len(sequence) != polymer.num_monomers:
raise ValueError("Length of sequence does not match length of polymer")
+
+ try:
+ polymer.idp_array
+ if len(polymer.idp_array) != polymer.num_monomers:
+ raise Exception
+ self.idp_array = polymer.idp_array
+ except:
+ print("Warning: KhBeadsFromPolymer processing a polymer without 'idp_array' set; assuming IDP")
+ self.idp_array = np.ones(len(polymer.num_monomers),dtype=np.bool)
+
Group.__init__(self, **kwargs)
def _clear_beads(self):
@@ -205,6 +232,8 @@ class KhBeadsFromPolymer(Group):
c = self.polymer.monomer_index_to_contour(i)
r = self.polymer.contour_to_position(c)
s = self.sequence[i]
+ if self.idp_array[i]:
+ s = s + 'IDP'
bead = PointParticle(_types[s], r,
name = s,
@@ -232,6 +261,7 @@ class KhModel(ArbdModel):
sequences = None,
debye_length = 10,
damping_coefficient = 10,
+ idp_array = None,
DEBUG=False,
**kwargs):
@@ -242,8 +272,8 @@ class KhModel(ArbdModel):
kwargs['timestep'] = 10e-6
kwargs['cutoff'] = max(4*debye_length,20)
- if 'decompPeriod' not in kwargs:
- kwargs['decompPeriod'] = 1000
+ if 'decomp_period' not in kwargs:
+ kwargs['decomp_period'] = 1000
""" Assign sequences """
if sequences is None:
@@ -253,21 +283,25 @@ class KhModel(ArbdModel):
self.sequences = sequences
ArbdModel.__init__(self, [], **kwargs)
+
""" Update type diffusion coefficients """
self.types = all_types = [t for key,t in _types.items()]
self.set_damping_coefficient( damping_coefficient )
""" Set up nonbonded interactions """
- nonbonded = KhNonbonded(debye_length)
- for i in range(len(all_types)):
- t1 = all_types[i]
- for j in range(i,len(all_types)):
- t2 = all_types[j]
- self.useNonbondedScheme( nonbonded, typeA=t1, typeB=t2 )
+ self.kh_nonbonded = KhNonbonded(debye_length)
+ for t in all_types:
+ self._add_nonbonded_interaction(t)
""" Generate beads """
self.generate_beads()
+ def _add_nonbonded_interaction(self, type_):
+ i = self.types.index(type_) if type_ in self.types else 0
+ for j in range(i,len(self.types)):
+ t = self.types[j]
+ self.useNonbondedScheme( self.kh_nonbonded, typeA=type_, typeB=t )
+
def update_splines(self, coords):
i = 0
for p in self.polymer_group.polymers:
diff --git a/arbdmodel/kh_polymer_model_pair_epsilon.py b/arbdmodel/kh_polymer_model_pair_epsilon.py
index e696b2c..e84acd9 100644
--- a/arbdmodel/kh_polymer_model_pair_epsilon.py
+++ b/arbdmodel/kh_polymer_model_pair_epsilon.py
@@ -224,6 +224,14 @@ for line in _txtdata.split('\n'):
__add_value((n1,n2), v)
__add_value((n2,n1), v)
+
+ __add_value((n1+'IDP',n2), v)
+ __add_value((n2,n1+'IDP'), v)
+ __add_value((n2+'IDP',n1), v)
+ __add_value((n1,n2+'IDP'), v)
+ __add_value((n2+'IDP',n1+'IDP'), v)
+ __add_value((n1+'IDP',n2+'IDP'), v)
+
if __name__ == '__main__':
print(epsilon_mj)
diff --git a/arbdmodel/onck_polymer_model.py b/arbdmodel/onck_polymer_model.py
index 861f90d..a63429d 100644
--- a/arbdmodel/onck_polymer_model.py
+++ b/arbdmodel/onck_polymer_model.py
@@ -270,8 +270,8 @@ class OnckModel(ArbdModel):
kwargs['timestep'] = 20e-6
kwargs['cutoff'] = max(5*debye_length,25)
- if 'decompPeriod' not in kwargs:
- kwargs['decompPeriod'] = 1000
+ if 'decomp_period' not in kwargs:
+ kwargs['decomp_period'] = 1000
""" Assign sequences """
if sequences is None:
--
GitLab