Optionally write PQR files

bin/enrgmd

@@ -30,6 +30,9 @@ parser.add_argument('--backbone-scale', type=float, default=1.0,
 parser.add_argument('--enrg-md-steps', type=int, default=1e6,
                     help='Number of ENRG-MD steps')
 
+parser.add_argument('--write-pqr', action='store_true',
+                    help='Write PQR file in addition PDB?')
+
 parser.add_argument('input_file', type=str,
                     help="""Any of the following:
                     (1) a cadnano JSON file;
@@ -100,7 +103,7 @@ http://dx.doi.org/10.1093/nar/gkw155
             os.makedirs('output')
         except FileExistsError:
             ...
-        model.atomic_simulate( output_name = prefix, num_steps=args.enrg_md_steps, dry_run=True )
+        model.atomic_simulate( output_name = prefix, num_steps=args.enrg_md_steps, write_pqr = args.write_pqr, dry_run=True )
         
     except:
         raise

bin/mrdna

@@ -74,6 +74,9 @@ parser.add_argument('--enrg-md-steps', type=int, default=1e6,
 parser.add_argument('--debug', action='store_true',
                     help='Run through the python debugger?')
 
+parser.add_argument('--write-pqr', action='store_true',
+                    help='Write PQR file in addition PDB?')
+
 parser.add_argument('--draw-cylinders', action='store_true',
                     help='Whether or not to draw the cylinders')
 parser.add_argument('--draw-tubes', action='store_true',
@@ -154,7 +157,7 @@ def main():
         model._clear_beads()
         model._generate_atomic_model(scale=args.backbone_scale)            
         model.write_atomic_ENM( prefix )
-        model.atomic_simulate( output_name = prefix, dry_run=True )
+        model.atomic_simulate( output_name = prefix, write_pqr=args.write_pqr, dry_run=True )
     elif args.draw_tubes is True:
         output_name = prefix + ".tubes.tcl"
         model.vmd_tube_tcl( output_name )
@@ -166,7 +169,7 @@ def main():
         model._clear_beads()
         model._generate_atomic_model(scale=args.backbone_scale)            
         model.write_atomic_ENM( prefix )
-        model.atomic_simulate( output_name = prefix, dry_run=True )
+        model.atomic_simulate( output_name = prefix, write_pqr=args.write_pqr, dry_run=True )
     else:
         run_args = dict(
             model = model,
@@ -187,6 +190,7 @@ def main():
             oxdna_steps = args.oxdna_steps,
             oxdna_output_period = args.oxdna_output_period,
             enrg_md_steps = args.enrg_md_steps,
+            write_pqr = args.write_pqr,
             run_enrg_md = args.run_enrg_md
         )
 

mrdna/arbdmodel/__init__.py

@@ -665,6 +665,36 @@ class PdbModel(Transformable, Parent):
 
         return
 
+    def write_pqr(self, filename):
+        if self.cacheInvalid:
+            self._updateParticleOrder()
+        with open(filename,'w') as fh:
+            ## Write header
+            fh.write("CRYST1{:>9.3f}{:>9.3f}{:>9.3f}  90.00  90.00  90.00 P 1           1\n".format( *self.dimensions ))
+
+            ## Write coordinates
+            formatString = "ATOM {idx:>6.6s} {name:^4.4s} {resname:3.3s} {chain:1.1s} {resid:>5.5s}   {x:.6f} {y:.6f} {z:.6f} {charge} {radius}\n"
+            for p in self.particles:
+                data = p._get_psfpdb_dictionary()
+
+                idx = data['idx']
+                if np.log10(idx) >= 5:
+                    idx = " *****"
+                else:
+                    idx = "{:>6d}".format(idx)
+                data['idx'] = idx
+
+                x,y,z = p.collapsedPosition()
+                data['x'] = x
+                data['y'] = y
+                data['z'] = z
+                assert(data['resid'] < 1e5)
+                data['resid'] = "{:<4d}".format(data['resid'])
+                if 'radius' not in data:
+                    data['radius'] = 2 * (data['mass']/16)**0.333333
+                fh.write( formatString.format(**data) )
+        return
+
     def writePsf(self, filename):
         if self.cacheUpToDate == False:
             self._updateParticleOrder()
@@ -876,7 +906,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, 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, write_pqr=False, log_file=None, dry_run = False):
         assert(type(gpu) is int)
         num_steps = int(num_steps)
 
@@ -923,6 +953,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 )
             # os.sync()
@@ -1511,6 +1542,7 @@ run {num_steps:d}
         if output_directory == '': output_directory='.'
         self.writePdb( output_name + ".pdb" )
         self.writePdb( output_name + ".fixed.pdb", beta_from_fixed=True )
+        if write_pqr: self.write_pqr( output_name + ".pqr" )        
         self.writePsf( output_name + ".psf" )
         self.write_namd_configuration( output_name, output_directory = output_directory, minimization_steps=minimization_steps, num_steps=num_steps )
 

mrdna/simulate.py

@@ -88,6 +88,7 @@ def multiresolution_simulation( model, output_name,
                                 backbone_scale=1.0,
                                 oxdna_steps = None,
                                 oxdna_output_period = None,
+                                write_pqr = False,
                                 run_enrg_md = False,
                                 enrg_md_steps = 1e6,
                                 arbd = None
@@ -156,13 +157,15 @@ def multiresolution_simulation( model, output_name,
 
         """ Minimization """
         if minimization_steps > 0:
-            simargs = dict(timestep=200e-6, output_period=minimization_output_period, gpu=gpu, arbd=arbd)
+            simargs = dict(timestep=200e-6, output_period=minimization_output_period, write_pqr=write_pqr,
+                           gpu=gpu, arbd=arbd)
             run_step(num_steps=minimization_steps, simargs=simargs)
             model.clear_beads()
             model.generate_bead_model( **coarse_model_args )
 
         """ Coarse simulation """
-        simargs = dict(timestep=100e-6 if coarse_local_twist else 200e-6, output_period=coarse_output_period, gpu=gpu, arbd=arbd)
+        simargs = dict(timestep=100e-6 if coarse_local_twist else 200e-6, output_period=coarse_output_period,
+                       write_pqr=write_pqr, gpu=gpu, arbd=arbd)
         if bond_cutoff > 0 and minimization_steps <= 0:
             run_step(num_steps=0.25*coarse_steps, simargs=simargs)
             model.clear_beads()
@@ -172,7 +175,7 @@ def multiresolution_simulation( model, output_name,
             run_step(num_steps=coarse_steps, simargs=simargs)
 
         """ Fine simulation """ 
-        simargs = dict(timestep=40e-6, output_period=fine_output_period, gpu=gpu, arbd=arbd)
+        simargs = dict(timestep=40e-6, output_period=fine_output_period, write_pqr=write_pqr, gpu=gpu, arbd=arbd)
         if not fix_linking_number:
             model.clear_beads()
             model.generate_bead_model( 1, 1, local_twist=True, escapable_twist=True )
@@ -204,7 +207,7 @@ def multiresolution_simulation( model, output_name,
             except FileExistsError:
                 pass
 
-            model.atomic_simulate( output_name = output_prefix, num_steps=enrg_md_steps, dry_run = not run_enrg_md )
+            model.atomic_simulate( output_name = output_prefix, num_steps=enrg_md_steps, write_pqr=write_pqr, dry_run = not run_enrg_md )