Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Funnel ABF and Funnel Metadynamics #302

Draft
wants to merge 17 commits into
base: main
Choose a base branch
from
Draft

Funnel ABF and Funnel Metadynamics #302

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
1be24c7
Add funnel methods and function
gustavor101 Feb 14, 2024
72bcae4
Add funnel_cv
gustavor101 Feb 14, 2024
9681c97
Add funnel loggers
gustavor101 Feb 14, 2024
ad11e68
add funnel in init
gustavor101 Feb 14, 2024
40d228d
Add examples
gustavor101 Feb 14, 2024
8d1dd2a
Add files via upload
gustavor101 Feb 23, 2024
71c34b6
Add files via upload
gustavor101 Feb 23, 2024
29fcfd0
Add files via upload
gustavor101 Feb 23, 2024
db1f036
Add files via upload
gustavor101 Feb 23, 2024
926d7c4
Add files via upload
gustavor101 Feb 23, 2024
a8b859e
Update run.py
gustavor101 Mar 14, 2024
7c7f52f
Update funnel_function.py
gustavor101 Mar 14, 2024
02b63e8
Merge branch 'SSAGESLabs:main' into funnel_abf
gustavor101 Aug 20, 2024
1a974c7
Add files via upload
gustavor101 Aug 21, 2024
cbaddec
Add files via upload
gustavor101 Aug 21, 2024
89b035c
Update funnel_abf.py
gustavor101 Aug 23, 2024
fb6bfad
Merge branch 'SSAGESLabs:main' into funnel_abf
gustavor101 Oct 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file added examples/openmm/funnel_abf/complex-wat-prod.rst
View file
Open in desktop
Binary file not shown.
33,349 changes: 33,349 additions & 0 deletions examples/openmm/funnel_abf/complex-wat.prmtop
View file
Open in desktop

Large diffs are not rendered by default.

144 changes: 144 additions & 0 deletions examples/openmm/funnel_abf/referencesh.pdb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
205 changes: 205 additions & 0 deletions examples/openmm/funnel_abf/run.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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)

# %%
# %%
Binary file added examples/openmm/funnel_metad/complex-wat-prod.rst
View file
Open in desktop
Binary file not shown.
Loading
Loading