Fixed psf writer for atomsitic simulations; add NAMD configuration writer; fixed bug with ENM

arbdmodel.py

@@ -545,7 +545,7 @@ class PdbModel(Transformable, Parent):
                     counter = 0
                 else:
                     fh.write(" ")
-            fh.write("\n")
+            fh.write("\n" if counter == 0 else "\n\n")
 
             ## Write out angles
             angles = self.get_angles()
@@ -559,7 +559,7 @@ class PdbModel(Transformable, Parent):
                     counter = 0
                 else:
                     fh.write(" ")
-            fh.write("\n")
+            fh.write("\n" if counter == 0 else "\n\n")
 
             ## Write out dihedrals
             dihedrals = self.get_dihedrals()
@@ -572,8 +572,8 @@ class PdbModel(Transformable, Parent):
                     fh.write("\n")
                     counter = 0
                 else:
-                    fh.write(" ")
-            fh.write("\n")
+                    fh.write(" ") 
+            fh.write("\n" if counter == 0 else "\n\n")
 
             ## Write out impropers
             impropers = self.get_impropers()
@@ -587,7 +587,17 @@ class PdbModel(Transformable, Parent):
                     counter = 0
                 else:
                     fh.write(" ")
-            fh.write("\n")
+            fh.write("\n" if counter == 0 else "\n\n")
+
+            fh.write("\n{:>8d} !NDON: donors\n\n\n".format(0))
+            fh.write("\n{:>8d} !NACC: acceptors\n\n\n".format(0))
+            fh.write("\n       0 !NNB\n\n")
+            natoms = len(self.particles)
+            for i in range(natoms//8):
+                fh.write("      0       0       0       0       0       0       0       0\n")
+            for i in range(natoms-8*(natoms//8)):
+                fh.write("      0")
+            fh.write("\n\n       1       0 !NGRP\n\n")
 
 
 class ArbdModel(PdbModel):
@@ -960,6 +970,137 @@ component "data" value 3
                 item = tuple(int(p.idx) for p in ex)
                 fh.write("EXCLUDE %d %d\n" % item)
 
+    def set_dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
+        ## TODO: cache coordinates using numpy arrays for quick min/max
+        raise(NotImplementedError)
+
+    def write_namd_configuration( self, output_prefix, num_steps = 1e6,
+                                  output_directory = 'output',
+                                  update_dimensions=True, extrabonds=True ):
+
+        format_data = self.__dict__.copy() # get parameters from System object
+
+        format_data['extrabonds'] = """extraBonds on
+extraBondsFile $prefix.exb
+""" if extrabonds else ""
+
+        if self.useTclForces:
+            format_data['margin'] = ""
+            format_data['tcl_forces'] = """tclForces on
+tclForcesScript $prefix.forces.tcl
+"""
+        else:
+            format_data['margin'] = """margin              30
+"""
+            format_data['tcl_forces'] = ""
+
+        if update_dimensions:
+            self.set_dimensions_from_structure()
+
+        for k,v in zip('XYZ', self.dimensions):
+            format_data['origin'+k] = -v*0.5
+            format_data['cell'+k] = v
+
+        format_data['prefix'] = output_prefix
+        format_data['num_steps'] = int(num_steps//12)*12
+        format_data['output_directory'] = output_directory
+        filename = '{}.namd'.format(output_prefix)
+
+        with open(filename,'w') as fh:
+            fh.write("""
+set prefix {prefix}
+set nLast 0;			# increment when continueing a simulation
+set n [expr $nLast+1]
+set out {output_directory}/$prefix-$n
+set temperature    300
+
+structure          $prefix.psf
+coordinates        $prefix.pdb
+
+outputName         $out
+XSTfile            $out.xst
+DCDfile            $out.dcd
+
+#############################################################
+## SIMULATION PARAMETERS                                   ##
+#############################################################
+
+# Input
+paraTypeCharmm	    on
+parameters          charmm36.nbfix/par_all36_na.prm
+parameters	    charmm36.nbfix/par_water_ions_na.prm
+
+wrapAll             off
+
+# Force-Field Parameters
+exclude             scaled1-4
+1-4scaling          1.0
+switching           on
+switchdist           8
+cutoff              10
+pairlistdist        12
+{margin}
+
+# Integrator Parameters
+timestep            2.0  ;# 2fs/step
+rigidBonds          all  ;# needed for 2fs steps
+nonbondedFreq       1
+fullElectFrequency  3
+stepspercycle       12
+
+# PME (for full-system periodic electrostatics)
+PME                 no
+PMEGridSpacing      1.2
+
+# Constant Temperature Control
+langevin            on    ;# do langevin dynamics
+# langevinDamping     1   ;# damping coefficient (gamma); used in original study
+langevinDamping     0.1   ;# less friction for faster relaxation
+langevinTemp        $temperature
+langevinHydrogen    off    ;# don't couple langevin bath to hydrogens
+
+# output
+useGroupPressure    yes
+xstFreq             4800
+outputEnergies      4800
+dcdfreq             4800
+restartfreq         48000
+
+#############################################################
+## EXTRA FORCES                                            ##
+#############################################################
+
+# ENM and intrahelical extrabonds
+{extrabonds}
+{tcl_forces}
+
+#############################################################
+## RUN                                                     ##
+#############################################################
+
+# Continuing a job from the restart files
+cellBasisVector1 {cellX} 0 0
+cellBasisVector2 0 {cellY} 0
+cellBasisVector3 0 0 {cellZ}
+
+if {{$nLast == 0}} {{
+    temperature 300
+    fixedAtoms on
+    fixedAtomsForces on
+    fixedAtomsFile $prefix.pdb
+    fixedAtomsCol B
+    minimize 2400
+
+    fixedAtoms off
+    minimize 2400
+}} else {{
+    bincoordinates  {output_directory}/$prefix-$nLast.restart.coor
+    binvelocities   {output_directory}/$prefix-$nLast.restart.vel
+}}
+
+run {num_steps:d}
+""".format(**format_data))
+
     def atomic_simulate(self, outputPrefix, outputDirectory='output'):
         if self.cacheUpToDate == False: # TODO: remove cache?
             self._countParticleTypes()
@@ -968,4 +1109,5 @@ component "data" value 3
         if outputDirectory == '': outputDirectory='.'
         self.writePdb( outputPrefix + ".pdb" )
         self.writePsf( outputPrefix + ".psf" )
+        self.write_namd_configuration( outputPrefix, output_directory = outputDirectory )
         os.sync()

segmentmodel.py

@@ -30,7 +30,7 @@ TODO:
  + remove performance bottlenecks
  - test for large systems
  + assign sequence
- - ENM
+ + ENM
  - rework Location class 
  - remove recursive calls
  - document
@@ -2244,7 +2244,7 @@ class SegmentModel(ArbdModel):
                         if nt2 is not None:
                             other.append((nt2,'stack'))
                             nt2 = nt2.basepair
-                            if nt2 is not None:
+                            if nt2 is not None and strand_piece.is_fwd:
                                 other.append((nt2,'cross'))
 
                         for nt2,key in other:
@@ -2392,6 +2392,18 @@ proc calcforces {} {
 }
 """)
 
+    def set_dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
+        positions = []
+        for s in self.segments:
+            positions.append(s.contour_to_position(0))
+            positions.append(s.contour_to_position(0.5))
+            positions.append(s.contour_to_position(1))
+        positions = np.array(positions)
+        dx,dy,dz = [(np.max(positions[:,i])-np.min(positions[:,i])+30)*padding_factor for i in range(3)]
+        if isotropic:
+            dx = dy = dz = max((dx,dy,dz))
+        self.dimensions = [dx,dy,dz]
+
     def _generate_atomic_model(self, scale=1):
         self.children = self.strands
         first_atomic_index = 0