Sync with arbdmodel 4784882; syntax changes have not yet been propogated to SegmentModel

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()

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):