From ac887d189651b6841a1c56b44facd1abcf45e4be Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Tue, 21 Aug 2018 16:01:58 -0500
Subject: [PATCH] Fixed psf writer for atomsitic simulations; add NAMD
configuration writer; fixed bug with ENM
---
arbdmodel.py | 152 ++++++++++++++++++++++++++++++++++++++++++++++--
segmentmodel.py | 16 ++++-
2 files changed, 161 insertions(+), 7 deletions(-)
diff --git a/arbdmodel.py b/arbdmodel.py
index 4d0a60c..a270c5b 100644
--- a/arbdmodel.py
+++ b/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()
diff --git a/segmentmodel.py b/segmentmodel.py
index bd12bac..3a6a458 100644
--- a/segmentmodel.py
+++ b/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
--
GitLab