From d3c116bf8aab3651e3c261a5243dcda021c60508 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Tue, 18 Feb 2020 18:03:40 -0600
Subject: [PATCH] Sync with arbdmodel 4784882; syntax changes have not yet been
propogated to SegmentModel
---
mrdna/arbdmodel/__init__.py | 80 +++++++++++++++++++++++++++------
mrdna/arbdmodel/interactions.py | 6 +--
2 files changed, 69 insertions(+), 17 deletions(-)
diff --git a/mrdna/arbdmodel/__init__.py b/mrdna/arbdmodel/__init__.py
index 1e396a0..1708211 100644
--- a/mrdna/arbdmodel/__init__.py
+++ b/mrdna/arbdmodel/__init__.py
@@ -674,11 +674,13 @@ class PdbModel(Transformable, Parent):
class ArbdModel(PdbModel):
- def __init__(self, children, dimensions=(1000,1000,1000), temperature=291,
- timestep=50e-6, particle_integrator = 'Brown',
- cutoff=50, decompPeriod=1000, pairlistDistance=None, nonbondedResolution=0.1,
+ 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=""):
+
PdbModel.__init__(self, children, dimensions, remove_duplicate_bonded_terms)
+ self.origin = origin
self.temperature = temperature
self.timestep = timestep
@@ -686,11 +688,11 @@ class ArbdModel(PdbModel):
self.particle_integrator = particle_integrator
- if pairlistDistance == None:
- pairlistDistance = cutoff+30
+ if pairlist_distance == None:
+ pairlist_distance = cutoff+30
- self.decompPeriod = decompPeriod
- self.pairlistDistance = pairlistDistance
+ self.decomp_period = decomp_period
+ self.pairlist_distance = pairlist_distance
self.extra_bd_file_lines = extra_bd_file_lines
@@ -752,9 +754,12 @@ class ArbdModel(PdbModel):
# self.initialCoords = np.array([p.initialPosition for p in self.particles])
def useNonbondedScheme(self, nbScheme, typeA=None, typeB=None):
- self.nbSchemes.append( (nbScheme, typeA, typeB) )
+ self.add_nonbonded_scheme(nbScheme, typeA, typeB)
+
+ def add_nonbonded_scheme(self, nonbonded_scheme, typeA=None, typeB=None):
+ self.nbSchemes.append( (nonbonded_scheme, typeA, typeB) )
if typeA != typeB:
- self.nbSchemes.append( (nbScheme, typeB, typeA) )
+ 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, log_file=None, dry_run = False):
assert(type(gpu) is int)
@@ -932,10 +937,16 @@ class ArbdModel(PdbModel):
params['outputPeriod'] = outputPeriod
for k,v in zip('XYZ', self.dimensions):
- params['origin'+k] = -v*0.5
params['dim'+k] = v
+
+ if self.origin is None:
+ for k,v in zip('XYZ', self.dimensions):
+ params['origin'+k] = -v*0.5
+ else:
+ for k,v in zip('XYZ', self.origin):
+ params['origin'+k] = v
- params['pairlistDistance'] -= params['cutoff']
+ params['pairlist_distance'] -= params['cutoff']
""" Find rigid groups """
rigid_groups = []
@@ -975,9 +986,9 @@ outputEnergyPeriod {outputPeriod}
outputFormat dcd
## Infrequent domain decomposition because this kernel is still very slow
-decompPeriod {decompPeriod}
+decompPeriod {decomp_period}
cutoff {cutoff}
-pairlistDistance {pairlistDistance}
+pairlistDistance {pairlist_distance}
origin {originX} {originY} {originZ}
systemSize {dimX} {dimY} {dimZ}
@@ -1313,7 +1324,7 @@ if {{$nLast == 0}} {{
run {num_steps:d}
""".format(**format_data))
- def atomic_simulate(self, output_name, output_directory='output'):
+ def atomic_simulate(self, output_name, output_directory='output', dry_run = False, namd2=None, log_file=None, num_procs=None, gpu=None):
if self.cacheUpToDate == False: # TODO: remove cache?
self._countParticleTypes()
self._updateParticleOrder()
@@ -1325,3 +1336,44 @@ run {num_steps:d}
self.write_namd_configuration( output_name, output_directory = output_directory )
os.sync()
+ if not dry_run:
+ if namd2 is None:
+ for path in os.environ["PATH"].split(os.pathsep):
+ path = path.strip('"')
+ fname = os.path.join(path, "namd2")
+ if os.path.isfile(fname) and os.access(fname, os.X_OK):
+ namd2 = fname
+ break
+
+ if namd2 is None: raise Exception("NAMD2 was not found")
+
+ if not os.path.exists(namd2):
+ raise Exception("NAMD2 was not found")
+ if not os.path.isfile(namd2):
+ raise Exception("NAMD2 was not found")
+ if not os.access(namd2, os.X_OK):
+ raise Exception("NAMD2 is not executable")
+
+ if not os.path.exists(output_directory):
+ os.makedirs(output_directory)
+ elif not os.path.isdir(output_directory):
+ raise Exception("output_directory '%s' is not a directory!" % output_directory)
+
+ if num_procs is None:
+ import multiprocessing
+ num_procs = max(1,multiprocessing.cpu_count()-1)
+
+ cmd = [namd2, '+p{}'.format(num_procs), "%s.namd" % output_name]
+ cmd = tuple(str(x) for x in cmd)
+
+ print("Running NAMD2 with: %s" % " ".join(cmd))
+ if log_file is None or (hasattr(log_file,'write') and callable(log_file.write)):
+ fd = sys.stdout if log_file is None else log_file
+ process = subprocess.Popen(cmd, stdout=subprocess.PIPE, universal_newlines=True)
+ for line in process.stdout:
+ fd.write(line)
+ fd.flush()
+ else:
+ with open(log_file,'w') as fd:
+ process = subprocess.Popen(cmd, stdout=log_file, universal_newlines=True)
+ process.communicate()
diff --git a/mrdna/arbdmodel/interactions.py b/mrdna/arbdmodel/interactions.py
index 8b7ae92..1c4b9b6 100644
--- a/mrdna/arbdmodel/interactions.py
+++ b/mrdna/arbdmodel/interactions.py
@@ -25,14 +25,14 @@ class NonbondedScheme():
class LennardJones(NonbondedScheme):
def potential(self, r, typeA, typeB):
- epsilon = sqrt( typeA.epsilon**2 + typeB.epsilon**2 )
+ epsilon = np.sqrt( typeA.epsilon**2 + typeB.epsilon**2 )
r0 = 0.5 * (typeA.radius + typeB.radius)
r6 = (r0/r)**6
r12 = r6**2
u = epsilon * (r12-2*r6)
u[0] = u[1] # Remove NaN
return u
-LennardJones = LennardJones()
+# LennardJones = LennardJones()
class HalfHarmonic(NonbondedScheme):
def potential(self, r, typeA, typeB):
@@ -41,7 +41,7 @@ class HalfHarmonic(NonbondedScheme):
u = 0.5 * k * (r-r0)**2
u[r > r0] = np.zeros( np.shape(u[r > r0]) )
return u
-HalfHarmonic = HalfHarmonic()
+# HalfHarmonic = HalfHarmonic()
class TabulatedPotential(NonbondedScheme):
def __init__(self, tableFile, typesA=None, typesB=None, resolution=0.1, rMin=0):
--
GitLab