Add files via upload

examples/openmm/funnel_metad/run_metad.py

@@ -1,82 +1,60 @@
 #!/usr/bin/env python3
 
+import dill as pickle
+import matplotlib.pyplot as plt
+import numpy as np
+from parmed import load_file
+from parmed.openmm import StateDataReporter
 
+import pysages
 from pysages.approxfun import compute_mesh
 from pysages.backends import SamplingContext
 from pysages.colvars import Projection_on_Axis_mobile
 
 # %%
-from pysages.methods import (
-    CVRestraints,
-    Funnel_MetadLogger,
-    Funnel_Metadynamics,
-    get_funnel_force,
-)
+from pysages.methods import Funnel_MetadLogger, Funnel_Metadynamics, get_funnel_force
 from pysages.utils import try_import
 
 # %%
-# from IPython import get_ipython
 openmm = try_import("openmm", "simtk.openmm")
 unit = try_import("openmm.unit", "simtk.unit")
 app = try_import("openmm.app", "simtk.openmm.app")
-import math
-
-import dill as pickle
-import jax.numpy as np
-
-# import numpy as np
-import matplotlib.pyplot as plt
-import numpy
-from openmm import *
-from openmm.unit import *
-
-# ParmEd Imports
-from parmed import load_file
-from parmed import unit as u
-from parmed.openmm import NetCDFReporter, StateDataReporter
-
-import pysages
 
 # %%
-pi = numpy.pi
-kB = 0.008314462618
+pi = np.pi
+kB = 0.008314462618  # kJ/mol
 
 
 # %%
-def generate_simulation(T=298.15 * kelvin, dt=1.0 * femtoseconds):
+def generate_simulation(T=298.15 * unit.kelvin, dt=1.0 * unit.femtoseconds):
     print("Loading AMBER files...")
     ala2_solv = load_file("complex-wat.prmtop", "complex-wat-prod.rst")
-    pdb = app.PDBFile("last.pdb")
+    #    pdb = app.PDBFile("last.pdb")
     system = ala2_solv.createSystem(
         nonbondedMethod=app.PME,
         rigidWater=True,
-        switchDistance=1.0 * nanometer,
-        nonbondedCutoff=1.2 * nanometer,
+        switchDistance=1.0 * unit.nanometer,
+        nonbondedCutoff=1.2 * unit.nanometer,
         constraints=app.HBonds,
     )
     # Create the integrator to do Langevin dynamics
     integrator = openmm.LangevinIntegrator(
-        300 * kelvin,  # Temperature of heat bath
-        1.0 / picoseconds,  # Friction coefficient
-        1.0 * femtoseconds,  # Time step
+        300 * unit.kelvin,  # Temperature of heat bath
+        1.0 / unit.picoseconds,  # Friction coefficient
+        1.0 * unit.femtoseconds,  # Time step
     )
 
     # Define the platform to use; CUDA, OpenCL, CPU, or Reference. Or do not specify
-    platform = Platform.getPlatformByName("CPU")
+    platform = openmm.Platform.getPlatformByName("CPU")
     # Create the Simulation object
 
     sim = app.Simulation(ala2_solv.topology, system, integrator, platform)
 
     # Set the particle positions
-    #    sim.context.setPositions(ala2_solv.positions)
-    sim.context.setPositions(pdb.getPositions(frame=-1))
+    sim.context.setPositions(ala2_solv.positions)
     sim.reporters.append(app.PDBReporter("output.pdb", 200000))
     sim.reporters.append(app.DCDReporter("output.dcd", 200000))
     sim.reporters.append(StateDataReporter("data.txt", 20000, step=True, separator=" "))
-
-    # Minimize the energy
-    #    print('Minimizing energy')
-    #    sim.minimizeEnergy()
     return sim
 
 
@@ -109,7 +87,7 @@ def save_energy_forces(result, plot_grid):
     alpha = -(5.3 / (5.0))
     A = energy(mesh) * alpha
     surface = (A - A.min()).reshape(plot_grid.shape)
-    numpy.savetxt("FES.csv", numpy.column_stack([mesh, surface]))
+    np.savetxt("FES.csv", np.column_stack([mesh, surface]))
 
 
 # %%
@@ -123,7 +101,6 @@ def save_energy_forces(result, plot_grid):
 indices_sys = [ligand, host, [anchor]]
 A = [4.6879, 4.8335, 1.12045]
 B = [4.7829, 5.7927, 2.8729]
-# anchor = 89
 box = [5.5822, 5.5095, 5.4335]
 Z_0 = 0.610225
 Zcc = 1.2
@@ -147,13 +124,10 @@ def save_energy_forces(result, plot_grid):
             temp_pos.append(float(p) / 10.0)
         ref_loop1.append(temp_pos)
 coordinates.close()
-# print(ref_loop1)
 # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-# print(indices_sys)
 cvs = (
     Projection_on_Axis_mobile(
         indices_sys,
@@ -166,9 +140,6 @@ def save_energy_forces(result, plot_grid):
     ),
 )
 
-grid = pysages.Grid(lower=(cv_min,), upper=(cv_max,), shape=(32,), periodic=False)
-
-restraints = CVRestraints(lower=(-0.2,), upper=(2.0,), kl=(10000.0,), ku=(10000.0,))
 # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 funnel_force = get_funnel_force(
@@ -190,16 +161,19 @@ def save_energy_forces(result, plot_grid):
 funnel_file = "funnel.dat"
 callback = Funnel_MetadLogger(funnel_file, 10)
 sampling_context = SamplingContext(methodo, generate_simulation, callback)
-# la linea de abajo se comenta para cuando se carga el pickle
 state = pysages.run(sampling_context, timesteps)
-# estas lineas se descomentan para usar el pickle para continuar la simulacion
+# restart from pickle file
+
 # with open("restartb.pickle", "rb") as f:
-#        state = pickle.load(f)
+#    state = pickle.load(f)
 # state = pysages.run(state, generate_simulation, timesteps)
-# se le puede cambiar e nombre del restart para no sobreescribirlo
+
+# save pickle
 with open("restartb.pickle", "wb") as f:
     pickle.dump(state, f)
 result = pysages.analyze(state)
+# for plotting
+grid = pysages.Grid(lower=(cv_min,), upper=(cv_max,), shape=(32,), periodic=False)
 plot_energy(result, grid)
 save_energy_forces(result, grid)