Update drude_openmm_run.py

transformato/bin/drude_openmm_run.py

@@ -94,10 +94,10 @@
     10 / picosecond,
     1 * kelvin,
     200 / picosecond,
-    0.0005 * picoseconds,
+    0.0001 * picoseconds,
 )
 
-# integrator.setMaxDrudeDistance(0.2 * angstroms)
+integrator.setMaxDrudeDistance(0.2 * angstroms)
 # if integrator.getMaxDrudeDistance() == 0:
 #     print("No Drude Hard Wall Contraint in use")
 # else:
@@ -153,8 +153,8 @@
 
 # Energy minimization
 if inputs.mini_nstep > 0:
-    print("\nEnergy minimization: %s steps" % inputs.mini_nstep)
-    simulation.minimizeEnergy(maxIterations=inputs.mini_nstep)
+    print("\nEnergy minimization:")
+    simulation.minimizeEnergy()
     print(simulation.context.getState(getEnergy=True).getPotentialEnergy())
 
 # Generate initial velocities
@@ -164,6 +164,29 @@
         simulation.context.setVelocitiesToTemperature(inputs.gen_temp, inputs.gen_seed)
     else:
         simulation.context.setVelocitiesToTemperature(inputs.gen_temp)
+
+## Do some additional pre-equilibration when using Drude particles
+print("Doing a first equilibration run")
+simulation.step(100_000)
+
+print("Doing a second equilibration run")
+simulation.integrator.setStepSize(0.0002 * picoseconds)
+simulation.context.reinitialize(preserveState=True)
+simulation.step(100_000)
+
+print("Doing a third equilibration run")
+simulation.integrator.setStepSize(0.0003 * picoseconds)
+simulation.context.reinitialize(preserveState=True)
+simulation.step(100_000)
+
+print("Doing a fourth equilibration run")
+simulation.integrator.setStepSize(0.0004 * picoseconds)
+simulation.context.reinitialize(preserveState=True)
+simulation.step(100_000)
+
+print("Starting the actual simulation")
+simulation.integrator.setStepSize(0.0005 * picoseconds)
+simulation.context.reinitialize(preserveState=True)
 
 # Production
 print("\nMD run: %s steps" % inputs.nstep)