Funnel ABF and Funnel Metadynamics

examples/openmm/funnel_abf/referencesh.pdb

@@ -0,0 +1,144 @@
+ATOM      1  O1  UNL     1      42.960  51.309  12.991  1.00  0.00           O  
+ATOM      2  C2  UNL     1      42.250  51.087  11.987  1.00  0.00           C  
+ATOM      3  N3  UNL     1      42.132  51.881  10.904  1.00  0.00           N  
+ATOM      4  C4  UNL     1      41.211  51.345   9.905  1.00  0.00           C  
+ATOM      5  C5  UNL     1      42.624  53.276  10.905  1.00  0.00           C  
+ATOM      6  N6  UNL     1      43.768  53.520  10.076  1.00  0.00           N  
+ATOM      7  C7  UNL     1      43.644  53.715   8.695  1.00  0.00           C  
+ATOM      8  C8  UNL     1      45.003  53.907  10.532  1.00  0.00           C  
+ATOM      9  O9  UNL     1      45.290  53.998  11.710  1.00  0.00           O  
+ATOM     10  N10 UNL     1      45.722  54.401   9.503  1.00  0.00           N  
+ATOM     11  C11 UNL     1      45.007  54.311   8.292  1.00  0.00           C  
+ATOM     12  C12 UNL     1      46.986  55.077   9.524  1.00  0.00           C  
+ATOM     13  N13 UNL     1      48.042  54.188   9.281  1.00  0.00           N  
+ATOM     14  C14 UNL     1      48.656  54.033   7.916  1.00  0.00           C  
+ATOM     15  C15 UNL     1      48.857  53.624  10.205  1.00  0.00           C  
+ATOM     16  O16 UNL     1      48.667  53.694  11.389  1.00  0.00           O  
+ATOM     17  N17 UNL     1      50.057  53.272   9.625  1.00  0.00           N  
+ATOM     18  C18 UNL     1      50.004  53.417   8.137  1.00  0.00           C  
+ATOM     19  C19 UNL     1      51.300  53.097  10.357  1.00  0.00           C  
+ATOM     20  N20 UNL     1      51.659  51.702  10.495  1.00  0.00           N  
+ATOM     21  C21 UNL     1      52.283  50.954   9.445  1.00  0.00           C  
+ATOM     22  C22 UNL     1      51.746  50.929  11.635  1.00  0.00           C  
+ATOM     23  O23 UNL     1      51.447  51.296  12.742  1.00  0.00           O  
+ATOM     24  N24 UNL     1      52.516  49.791  11.422  1.00  0.00           N  
+ATOM     25  C25 UNL     1      52.818  49.627  10.118  1.00  0.00           C  
+ATOM     26  C26 UNL     1      53.069  48.967  12.433  1.00  0.00           C  
+ATOM     27  N27 UNL     1      52.530  47.616  12.524  1.00  0.00           N  
+ATOM     28  C28 UNL     1      53.033  46.490  11.734  1.00  0.00           C  
+ATOM     29  C29 UNL     1      51.692  47.153  13.555  1.00  0.00           C  
+ATOM     30  O30 UNL     1      51.044  47.909  14.261  1.00  0.00           O  
+ATOM     31  N31 UNL     1      51.760  45.782  13.668  1.00  0.00           N  
+ATOM     32  C32 UNL     1      52.511  45.265  12.526  1.00  0.00           C  
+ATOM     33  C33 UNL     1      51.481  44.908  14.823  1.00  0.00           C  
+ATOM     34  N34 UNL     1      50.244  44.227  14.701  1.00  0.00           N  
+ATOM     35  C35 UNL     1      50.174  42.937  14.092  1.00  0.00           C  
+ATOM     36  C36 UNL     1      49.106  44.646  15.317  1.00  0.00           C  
+ATOM     37  O37 UNL     1      48.921  45.709  15.883  1.00  0.00           O  
+ATOM     38  N38 UNL     1      48.315  43.598  15.310  1.00  0.00           N  
+ATOM     39  C39 UNL     1      48.800  42.454  14.514  1.00  0.00           C  
+ATOM     40  C40 UNL     1      47.013  43.358  15.950  1.00  0.00           C  
+ATOM     41  N41 UNL     1      45.844  43.751  15.153  1.00  0.00           N  
+ATOM     42  C42 UNL     1      45.078  42.645  14.428  1.00  0.00           C  
+ATOM     43  C43 UNL     1      45.117  44.914  15.225  1.00  0.00           C  
+ATOM     44  O44 UNL     1      45.452  45.987  15.672  1.00  0.00           O  
+ATOM     45  N45 UNL     1      43.850  44.583  14.778  1.00  0.00           N  
+ATOM     46  C46 UNL     1      43.796  43.291  14.049  1.00  0.00           C  
+ATOM     47  C47 UNL     1      42.592  45.339  15.160  1.00  0.00           C  
+ATOM     48  N48 UNL     1      42.173  46.213  14.055  1.00  0.00           N  
+ATOM     49  C49 UNL     1      41.390  45.635  12.903  1.00  0.00           C  
+ATOM     50  C50 UNL     1      42.070  47.549  14.091  1.00  0.00           C  
+ATOM     51  O51 UNL     1      42.570  48.267  14.930  1.00  0.00           O  
+ATOM     52  N52 UNL     1      41.298  47.942  13.082  1.00  0.00           N  
+ATOM     53  C53 UNL     1      40.702  49.222  12.929  1.00  0.00           C  
+ATOM     54  C54 UNL     1      40.788  46.928  12.195  1.00  0.00           C  
+ATOM     55  N55 UNL     1      41.313  50.087  11.855  1.00  0.00           N  
+ATOM     56  C56 UNL     1      40.786  50.060  10.538  1.00  0.00           C  
+ATOM     57  N57 UNL     1      41.379  48.997   9.731  1.00  0.00           N  
+ATOM     58  C58 UNL     1      40.809  47.659   9.741  1.00  0.00           C  
+ATOM     59  C59 UNL     1      42.022  49.526   8.611  1.00  0.00           C  
+ATOM     60  O60 UNL     1      42.427  48.897   7.667  1.00  0.00           O  
+ATOM     61  N61 UNL     1      41.865  50.902   8.713  1.00  0.00           N  
+ATOM     62  C62 UNL     1      42.242  51.798   7.631  1.00  0.00           C  
+ATOM     63  N63 UNL     1      43.498  52.478   7.874  1.00  0.00           N  
+ATOM     64  C64 UNL     1      44.603  52.172   7.193  1.00  0.00           C  
+ATOM     65  O65 UNL     1      44.781  51.203   6.461  1.00  0.00           O  
+ATOM     66  N66 UNL     1      45.560  53.131   7.403  1.00  0.00           N  
+ATOM     67  C67 UNL     1      46.705  53.360   6.501  1.00  0.00           C  
+ATOM     68  N68 UNL     1      47.982  53.028   7.114  1.00  0.00           N  
+ATOM     69  C69 UNL     1      48.756  51.900   6.906  1.00  0.00           C  
+ATOM     70  O70 UNL     1      48.451  50.946   6.233  1.00  0.00           O  
+ATOM     71  N71 UNL     1      50.051  52.169   7.393  1.00  0.00           N  
+ATOM     72  C72 UNL     1      51.162  51.302   7.264  1.00  0.00           C  
+ATOM     73  N73 UNL     1      51.365  50.453   8.432  1.00  0.00           N  
+ATOM     74  C74 UNL     1      51.185  49.106   8.517  1.00  0.00           C  
+ATOM     75  O75 UNL     1      50.473  48.404   7.794  1.00  0.00           O  
+ATOM     76  N76 UNL     1      52.020  48.555   9.528  1.00  0.00           N  
+ATOM     77  C77 UNL     1      52.586  47.217   9.364  1.00  0.00           C  
+ATOM     78  N78 UNL     1      52.311  46.442  10.558  1.00  0.00           N  
+ATOM     79  C79 UNL     1      51.477  45.332  10.445  1.00  0.00           C  
+ATOM     80  O80 UNL     1      50.875  44.995   9.444  1.00  0.00           O  
+ATOM     81  N81 UNL     1      51.659  44.572  11.545  1.00  0.00           N  
+ATOM     82  C82 UNL     1      51.262  43.188  11.803  1.00  0.00           C  
+ATOM     83  N83 UNL     1      50.090  43.081  12.606  1.00  0.00           N  
+ATOM     84  C84 UNL     1      48.834  42.802  12.162  1.00  0.00           C  
+ATOM     85  O85 UNL     1      48.402  43.040  11.024  1.00  0.00           O  
+ATOM     86  N86 UNL     1      48.113  42.343  13.235  1.00  0.00           N  
+ATOM     87  C87 UNL     1      46.922  41.524  13.010  1.00  0.00           C  
+ATOM     88  N88 UNL     1      45.734  42.391  13.169  1.00  0.00           N  
+ATOM     89  C89 UNL     1      45.027  42.920  12.103  1.00  0.00           C  
+ATOM     90  O90 UNL     1      45.392  42.935  10.936  1.00  0.00           O  
+ATOM     91  N91 UNL     1      43.878  43.435  12.609  1.00  0.00           N  
+ATOM     92  C92 UNL     1      42.676  43.698  11.878  1.00  0.00           C  
+ATOM     93  N93 UNL     1      42.273  45.038  11.788  1.00  0.00           N  
+ATOM     94  C94 UNL     1      42.244  45.800  10.712  1.00  0.00           C  
+ATOM     95  N95 UNL     1      41.362  46.840  10.845  1.00  0.00           N  
+ATOM     96  O96 UNL     1      42.785  45.526   9.650  1.00  0.00           O  
+ATOM     97  H97 UNL     1      40.394  52.073   9.844  1.00  0.00           H  
+ATOM     98  H98 UNL     1      42.724  53.563  11.958  1.00  0.00           H  
+ATOM     99  H99 UNL     1      41.824  53.861  10.435  1.00  0.00           H  
+ATOM    100 H100 UNL     1      42.728  54.288   8.509  1.00  0.00           H  
+ATOM    101 H101 UNL     1      44.842  55.283   7.815  1.00  0.00           H  
+ATOM    102 H102 UNL     1      47.020  55.494  10.536  1.00  0.00           H  
+ATOM    103 H103 UNL     1      46.925  55.789   8.693  1.00  0.00           H  
+ATOM    104 H104 UNL     1      48.576  54.988   7.384  1.00  0.00           H  
+ATOM    105 H105 UNL     1      50.853  54.028   7.811  1.00  0.00           H  
+ATOM    106 H106 UNL     1      51.160  53.542  11.349  1.00  0.00           H  
+ATOM    107 H107 UNL     1      52.001  53.612   9.691  1.00  0.00           H  
+ATOM    108 H108 UNL     1      53.120  51.482   8.973  1.00  0.00           H  
+ATOM    109 H109 UNL     1      53.866  49.324  10.019  1.00  0.00           H  
+ATOM    110 H110 UNL     1      54.133  48.801  12.227  1.00  0.00           H  
+ATOM    111 H111 UNL     1      52.972  49.450  13.412  1.00  0.00           H  
+ATOM    112 H112 UNL     1      54.106  46.488  11.512  1.00  0.00           H  
+ATOM    113 H113 UNL     1      53.285  44.628  12.972  1.00  0.00           H  
+ATOM    114 H114 UNL     1      52.271  44.154  14.918  1.00  0.00           H  
+ATOM    115 H115 UNL     1      51.450  45.511  15.738  1.00  0.00           H  
+ATOM    116 H116 UNL     1      51.012  42.277  14.344  1.00  0.00           H  
+ATOM    117 H117 UNL     1      48.654  41.458  14.947  1.00  0.00           H  
+ATOM    118 H118 UNL     1      46.956  42.275  16.110  1.00  0.00           H  
+ATOM    119 H119 UNL     1      47.025  43.862  16.923  1.00  0.00           H  
+ATOM    120 H120 UNL     1      44.997  41.716  15.004  1.00  0.00           H  
+ATOM    121 H121 UNL     1      42.850  42.781  14.264  1.00  0.00           H  
+ATOM    122 H122 UNL     1      41.694  44.752  15.383  1.00  0.00           H  
+ATOM    123 H123 UNL     1      42.911  46.008  15.968  1.00  0.00           H  
+ATOM    124 H124 UNL     1      40.635  44.931  13.271  1.00  0.00           H  
+ATOM    125 H125 UNL     1      40.765  49.752  13.886  1.00  0.00           H  
+ATOM    126 H126 UNL     1      39.648  49.147  12.638  1.00  0.00           H  
+ATOM    127 H127 UNL     1      39.697  47.031  12.156  1.00  0.00           H  
+ATOM    128 H128 UNL     1      39.712  49.844  10.530  1.00  0.00           H  
+ATOM    129 H129 UNL     1      39.722  47.711   9.868  1.00  0.00           H  
+ATOM    130 H130 UNL     1      40.998  47.202   8.763  1.00  0.00           H  
+ATOM    131 H131 UNL     1      41.486  52.577   7.481  1.00  0.00           H  
+ATOM    132 H132 UNL     1      42.363  51.132   6.770  1.00  0.00           H  
+ATOM    133 H133 UNL     1      46.677  54.374   6.087  1.00  0.00           H  
+ATOM    134 H134 UNL     1      46.613  52.685   5.643  1.00  0.00           H  
+ATOM    135 H135 UNL     1      52.043  51.950   7.181  1.00  0.00           H  
+ATOM    136 H136 UNL     1      51.006  50.611   6.428  1.00  0.00           H  
+ATOM    137 H137 UNL     1      52.030  46.844   8.496  1.00  0.00           H  
+ATOM    138 H138 UNL     1      53.674  47.269   9.248  1.00  0.00           H  
+ATOM    139 H139 UNL     1      51.129  42.640  10.863  1.00  0.00           H  
+ATOM    140 H140 UNL     1      52.114  42.730  12.319  1.00  0.00           H  
+ATOM    141 H141 UNL     1      47.045  41.123  11.998  1.00  0.00           H  
+ATOM    142 H142 UNL     1      46.851  40.827  13.853  1.00  0.00           H  
+ATOM    143 H143 UNL     1      42.928  43.415  10.850  1.00  0.00           H  
+ATOM    144 H144 UNL     1      41.833  43.057  12.160  1.00  0.00           H  

examples/openmm/funnel_abf/run.py

@@ -0,0 +1,205 @@
+#!/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.backends import SamplingContext
+from pysages.colvars import Projection_on_Axis_mobile
+
+# %%
+from pysages.methods import CVRestraints, Funnel_ABF, Funnel_Logger, get_funnel_force
+from pysages.utils import try_import
+
+# %%
+openmm = try_import("openmm", "simtk.openmm")
+unit = try_import("openmm.unit", "simtk.unit")
+app = try_import("openmm.app", "simtk.openmm.app")
+
+# %%
+pi = np.pi
+kB = 0.008314462618  # kJ/mol
+
+
+# %%
+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")
+    system = ala2_solv.createSystem(
+        nonbondedMethod=app.PME,
+        rigidWater=True,
+        switchDistance=1.0 * unit.nanometer,
+        nonbondedCutoff=1.2 * unit.nanometer,
+        constraints=app.HBonds,
+    )
+    # Create the integrator to do Langevin dynamics
+    integrator = openmm.LangevinIntegrator(
+        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 = 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.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=" "))
+
+    return sim
+
+
+# %%
+# functions for ploting and storing data
+def plot_energy(result):
+    fig, ax = plt.subplots()
+
+    ax.set_xlabel("CV")
+    ax.set_ylabel("Free energy $[\\epsilon]$")
+
+    free_energy = np.asarray(result["free_energy"])
+    x = np.asarray(result["mesh"])
+    ax.plot(x, free_energy, color="teal")
+
+    fig.savefig("energy.png")
+
+
+def plot_forces(result):
+    fig, ax = plt.subplots()
+
+    ax.set_xlabel("CV")
+    ax.set_ylabel("Forces $[\\epsilon]$")
+
+    forces = np.asarray(result["mean_force"])
+    x = np.asarray(result["mesh"])
+    ax.plot(x, forces, color="teal")
+
+    fig.savefig("forces.png")
+
+
+def plot_histogram(result):
+    fig, ax = plt.subplots()
+
+    ax.set_xlabel("CV")
+    ax.set_ylabel("Histogram $[\\epsilon]$")
+
+    hist = np.asarray(result["histogram"])
+    x = np.asarray(result["mesh"])
+    ax.plot(x, hist, color="teal")
+
+    fig.savefig("histogram2.png")
+
+
+def save_energy_forces(result):
+    Energy = np.asarray(result["free_energy"])
+    Forces = np.asarray(result["mean_force"])
+    Grid = np.asarray(result["mesh"])
+    hist = np.asarray(result["histogram"])
+    np.savetxt("FES.csv", np.column_stack([Grid, Energy]))
+    np.savetxt("Forces.csv", np.column_stack([Grid, Forces]))
+    np.savetxt("Histogram.csv", np.column_stack([Grid, hist]))
+
+
+#    numpy.savetxt("Forces.csv", numpy.column_stack([Grid, Forces]))
+
+# %%
+# %%
+# %%
+host = list(range(0, 144))
+ligand = list(range(144, 168))
+weights_host = np.ones(len(host)) / len(host)
+weights_lig = np.ones(len(ligand)) / len(ligand)
+anchor = 89
+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
+R_cyl = 0.6
+k_cone = 10000.0
+k_cv = 0.0
+cv_min = 0.0
+cv_max = 2.0
+cv_buffer = 0.05
+coordinates = open("referencesh.pdb", "r")
+ref_loop1 = []
+for line in coordinates:
+    lista = line.split()
+    id = lista[0]
+    if id == "ATOM":
+        atomid = int(lista[1])
+        residue = int(lista[4])
+        position = lista[5:8]
+        temp_pos = []
+        for p in position:
+            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,
+        references=ref_loop1,
+        weights_lig=None,
+        weights_prot=None,
+        A=A,
+        B=B,
+        box=box,
+    ),
+)
+
+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(
+    indices_sys, ref_loop1, A, B, Zcc, Z_0, R_cyl, k_cone, k_cv, cv_min, cv_max, cv_buffer, box
+)
+# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+height = 0.2
+sigma = 0.02
+stride = 500
+timesteps = 500
+ngauss = timesteps // stride + 1
+# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+methodo = Funnel_ABF(cvs, grid=grid, N=1000, ext_force=funnel_force, restraints=restraints)
+funnel_file = "funnel.dat"
+callback = Funnel_Logger(funnel_file, 10)
+sampling_context = SamplingContext(methodo, generate_simulation, callback)
+state = pysages.run(sampling_context, timesteps)
+# for restart simulation
+# with open("restart1.pickle", "rb") as f:
+#    state = pickle.load(f)
+# state = pysages.run(state, generate_simulation, timesteps)
+# save pickle file
+with open("restart1.pickle", "wb") as f:
+    pickle.dump(state, f)
+topology = (8, 8)
+result = pysages.analyze(state, topology=topology)
+plot_energy(result)
+plot_forces(result)
+plot_histogram(result)
+save_energy_forces(result)
+
+# %%
+# %%