box stabilisation

bin/openmmMD.py

@@ -1,5 +1,5 @@
 #!/usr/bin/env python3
-"""openmm-MD
+"""openmmMD
 
 Minimise the potential energy of the wildtype protein structure and mutated variants, according to AMBER14 potentials 
 Usage:
@@ -112,14 +112,14 @@ def setup_system(modeller, forcefield, solvmol: str, no_restraints: bool):
     return system
 
 def setup_simulation(modeller, system):
-    integrator = mm.LangevinMiddleIntegrator(10*kelvin, 1/picosecond, 0.001*picoseconds)
+    integrator = mm.LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.001*picoseconds)
     simulation = app.Simulation(modeller.topology, system, integrator)
     simulation.context.setPositions(modeller.positions)
     return simulation, integrator
 
 def energy_minimization(modeller):
     forcefield = setup_forcefield()
-    solvmol = 3000
+    solvmol = 500
     system = setup_system(modeller, forcefield, solvmol, arguments['--no_restraints'])
     simulation = setup_simulation(modeller, system)[0]
     integrator = setup_simulation(modeller, system)[1]
@@ -138,8 +138,8 @@ def md_nvt(simulation, csvname: str, totalsteps: int, reprate: int, pdbname, int
     init_pe = init_state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
     for temp in range(0, 300, 10):
         integrator.setTemperature(temp*kelvin)
-        simulation.step(1000)
-    simulation.context.setVelocitiesToTemperature(temp*kelvin)
+        simulation.step(100)
+    simulation.context.setVelocitiesToTemperature(300*kelvin)
     simulation.reporters.append(app.PDBReporter(pdbname, reprate))
     simulation.reporters.append(app.StateDataReporter(stdout, reprate, step=True, potentialEnergy=True, temperature=True, volume=True))
     prepdf = {'Step':[], 'Potential Energy (kJ/mole)':[], 'Temperature (K)':[], 'Box Volume (nm^3)':[]}
@@ -167,6 +167,7 @@ def md_nvt(simulation, csvname: str, totalsteps: int, reprate: int, pdbname, int
 
 def rmsf_analysis(pdb_traj: str, rmsfcsv: str):
     u = mda.Universe(pdb_traj)
+    u.trajectory.dt = 0.001
     average = align.AverageStructure(u, u, select='protein and name CA', ref_frame=0).run()
     ref = average.results.universe
     aligner = align.AlignTraj(u, ref, select='protein and name CA', in_memory=True).run()
@@ -179,7 +180,7 @@ def rmsf_analysis(pdb_traj: str, rmsfcsv: str):
     df.to_csv(rmsfcsv, index=False)
 
 def main():
-    arguments = docopt(__doc__, version='openmm-fb_MD.py')
+    arguments = docopt(__doc__, version='openmmMD.py')
     new_rep_cleaned = get_data(arguments['--incsv'], arguments['--from-col'])
     reformat_data(arguments['--incsv'], new_rep_cleaned)
     pdb = clean_wildtype(arguments['--in-pdb'], arguments['--pH'])
@@ -193,13 +194,13 @@ def main():
         integrator = wt_pdb[4]
         csvname = str("wt_traj" + strj + ".csv")
         pdbname = str("wt_traj" + strj + ".pdb")
-        sim_run = md_nvt(simulation, csvname, 200000, 10000, pdbname, integrator)
+        sim_run = md_nvt(simulation, csvname, 10000, 1000, pdbname, integrator)
         wt_ref = str("wt_reference" + strj + ".pdb")
         app.PDBFile.writeFile(sim_run[4], sim_run[5], open(wt_ref, "w"), keepIds=True)
         rmsfout = str("wt_rmsf" + strj + ".csv")
         rmsf_analysis(pdbname, rmsfout)
 
 if __name__ == '__main__':
-    arguments = docopt(__doc__, version='openmm-fb_MD.py')
+    arguments = docopt(__doc__, version='openmmMD.py')
     logging.getLogger().setLevel(logging.INFO)  
     main()