Skip to content

Commit c08d2cb

Browse files
kocimil1kocimil1
committed
Wrappers around computeFreeEnergy() function
1 parent 67988ea commit c08d2cb

5 files changed

Lines changed: 205 additions & 1 deletion

File tree

‎cpp/common/molecular/MolWorld_sp3_multi.h‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -208,6 +208,30 @@ class MolWorld_sp3_multi : public MolWorld_sp3, public MultiSolverInterface { pu
208208

209209
virtual int getMolWorldVersion() const override { return (int)MolWorldVersion::GPU; };
210210

211+
// ==================================
212+
// Free Energy Calculation
213+
// ==================================
214+
215+
double computeFreeEnergy(double lamda1, double lamda2, int n, int *dc, int nbStep = 100, int nMDsteps = 100000, int nEQsteps = 10000, double tdamp = 100.0, double T = 300, double dt = 0.5){
216+
printf("MolWorld_sp3_multi::computeFreeEnergy() - Starting\n");
217+
printf(" Parameters: lamda1=%g, lamda2=%g, n=%i, nbStep=%i, nMDsteps=%i, nEQsteps=%i, tdamp=%g, T=%g, dt=%g\n",
218+
lamda1, lamda2, n, nbStep, nMDsteps, nEQsteps, tdamp, T, dt);
219+
220+
// Save the loaded molecule to verify it was loaded correctly
221+
// Save each replica to a separate file to verify correct loading
222+
for (int iSys = 0; iSys < nSystems; ++iSys) {
223+
char verify_fname[256];
224+
sprintf(verify_fname, "loaded_DA_%i.xyz", iSys);
225+
printf(" Saving system %d to: %s\n", iSys, verify_fname);
226+
saveSysXYZ(iSys, verify_fname, "# Loaded DA.mol2 for free energy calculation", false, "w", Vec3i{1,1,1});
227+
printf(" System %d saved to %s\n", iSys, verify_fname);
228+
}
229+
230+
// TODO: Implement actual free energy calculation
231+
printf("MolWorld_sp3_multi::computeFreeEnergy() - Not yet fully implemented\n");
232+
return 0.0;
233+
}
234+
211235
// ==================================
212236
// Initialization
213237
// ==================================

‎cpp/common_resources/DA.mol2‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
@<TRIPOS>MOLECULE
2+
*****
3+
44 44 0 0 0
4+
SMALL
5+
GASTEIGER
6+
7+
@<TRIPOS>ATOM
8+
1 C -5.6868 -5.1707 0.8745 C.ar 1 UNL1 0.0000
9+
2 C -4.4174 -5.2554 1.4498 C.ar 1 UNL1 0.0000
10+
3 N -4.2571 -5.9048 2.6306 N.ar 1 UNL1 0.0000
11+
4 C -5.3167 -6.4737 3.2598 C.ar 1 UNL1 0.0000
12+
5 H -3.3072 -5.9677 3.0603 H 1 UNL1 0.0000
13+
6 C -6.5887 -6.3915 2.6877 C.ar 1 UNL1 0.0000
14+
7 N -6.7673 -5.7377 1.4935 N.ar 1 UNL1 0.0000
15+
8 H -5.1768 -6.9894 4.2016 H 1 UNL1 0.0000
16+
9 H -7.4072 -6.8563 3.2176 H 1 UNL1 0.0000
17+
10 H -5.8007 -4.6490 -0.0690 H 1 UNL1 0.0000
18+
11 H -3.5658 -4.8050 0.9556 H 1 UNL1 0.0000
19+
12 C -8.0756 -5.6018 0.8170 C.3 1 UNL1 0.0000
20+
13 C -9.2614 -6.2543 1.5493 C.3 1 UNL1 0.0000
21+
14 H -7.9992 -6.0594 -0.1947 H 1 UNL1 0.0000
22+
15 H -8.2984 -4.5183 0.6942 H 1 UNL1 0.0000
23+
16 C -10.5621 -6.0563 0.7702 C.3 1 UNL1 0.0000
24+
17 H -9.3752 -5.7920 2.5542 H 1 UNL1 0.0000
25+
18 H -9.0742 -7.3451 1.6576 H 1 UNL1 0.0000
26+
19 SI -12.0498 -6.8896 1.7148 Si 1 UNL1 0.0000
27+
20 H -10.4797 -6.5183 -0.2366 H 1 UNL1 0.0000
28+
21 H -10.7788 -4.9727 0.6561 H 1 UNL1 0.0000
29+
22 H -11.7760 -8.4075 1.8672 H 1 UNL1 0.0000
30+
23 H -13.3559 -6.6756 0.9084 H 1 UNL1 0.0000
31+
24 H -12.1934 -6.2454 3.1171 H 1 UNL1 0.0000
32+
25 C -1.6863 -7.5949 7.5698 C.ar 2 UNL2 0.0000
33+
26 C 0.1592 -6.6775 8.7442 C.ar 2 UNL2 0.0000
34+
27 C -2.1167 -6.3166 7.2064 C.ar 2 UNL2 0.0000
35+
28 C -0.2433 -5.3844 8.3988 C.ar 2 UNL2 0.0000
36+
29 C -1.3886 -5.2046 7.6250 C.ar 2 UNL2 0.0000
37+
30 H 0.3351 -4.5314 8.7332 H 2 UNL2 0.0000
38+
31 H -3.0095 -6.1974 6.6042 H 2 UNL2 0.0000
39+
32 N -0.5608 -7.7728 8.3301 N.ar 2 UNL2 0.0000
40+
33 H 1.0519 -6.7728 9.3446 H 2 UNL2 0.0000
41+
34 H -2.2643 -8.4493 7.2359 H 2 UNL2 0.0000
42+
35 C -0.1924 -9.1675 8.6585 C.3 2 UNL2 0.0000
43+
36 C 1.0785 -9.3267 9.5117 C.3 2 UNL2 0.0000
44+
37 H -0.0419 -9.7257 7.7075 H 2 UNL2 0.0000
45+
38 H -1.0391 -9.6361 9.2083 H 2 UNL2 0.0000
46+
39 H 1.9469 -8.8926 8.9702 H 2 UNL2 0.0000
47+
40 C 1.3699 -10.7998 9.7921 C.3 2 UNL2 0.0000
48+
41 H 0.9430 -8.8031 10.4830 H 2 UNL2 0.0000
49+
42 H 1.5232 -11.3501 8.8393 H 2 UNL2 0.0000
50+
43 H 0.5244 -11.2611 10.3455 H 2 UNL2 0.0000
51+
44 SI 3.1486 -10.9881 10.9822 Si 2 UNL2 0.0000
52+
45 O -1.7506 -4.0828 7.3167 O.2 2 UNL2 0.0000
53+
@<TRIPOS>BOND
54+
1 1 2 ar
55+
2 2 3 ar
56+
3 3 4 ar
57+
4 3 5 1
58+
5 4 6 ar
59+
6 6 7 ar
60+
7 1 7 ar
61+
8 4 8 1
62+
9 6 9 1
63+
10 1 10 1
64+
11 2 11 1
65+
12 7 12 1
66+
13 12 13 1
67+
14 12 14 1
68+
15 12 15 1
69+
16 13 16 1
70+
17 13 17 1
71+
18 13 18 1
72+
19 16 19 1
73+
20 16 20 1
74+
21 16 21 1
75+
22 19 22 1
76+
23 19 23 1
77+
24 19 24 1
78+
25 27 25 ar
79+
26 26 28 ar
80+
27 28 29 ar
81+
28 28 30 1
82+
29 27 29 ar
83+
30 27 31 1
84+
31 26 32 ar
85+
32 26 33 1
86+
33 32 25 ar
87+
34 25 34 1
88+
35 32 35 1
89+
36 35 36 1
90+
37 35 37 1
91+
38 35 38 1
92+
39 36 39 1
93+
40 36 40 1
94+
41 36 41 1
95+
42 40 42 1
96+
43 40 43 1
97+
44 40 44 1
98+
45 29 45 2

‎cpp/libs_OCL/MMFFmulti_lib.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -327,5 +327,8 @@ void setSwitches_multi( int CheckInvariants, int PBC, int NonBonded, int MMFF, i
327327
#undef _setbool
328328
}
329329

330+
double computeFreeEnergy(double lamda1, double lamda2, int n, int* dc, int nbStep, int nMDsteps, int nEQsteps, double tdamp, double T, double dt){
331+
return W.computeFreeEnergy(lamda1, lamda2, n, dc, nbStep, nMDsteps, nEQsteps, tdamp, T, dt);
332+
}
330333

331334
} // extern "C"

‎examples/tFreeEnergy_multi/run.py‎

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
import sys
2+
import numpy as np
3+
import matplotlib.pyplot as plt
4+
5+
sys.path.append("../../")
6+
from pyBall import atomicUtils as au
7+
from pyBall import MMFF_multi as mmff
8+
9+
# Example usage of computeFreeEnergy function
10+
def main():
11+
print("=== Free Energy Calculation Example ===")
12+
13+
# Initialize the system
14+
nSys = 10
15+
xyz_name = "data/DA.mol2" # Load the DA molecule pair
16+
17+
print(f"Initializing MMFF_multi with {nSys} systems")
18+
print(f"Loading molecule from: {xyz_name}")
19+
20+
# Parameters for free energy calculation
21+
lamda1 = 0.0 # Initial lambda value
22+
lamda2 = 1.0 # Final lambda value
23+
dc = [0, 1, 2] # Distance constraint indices (example)
24+
nbStep = 100 # Number of steps
25+
nMDsteps = 100000 # Number of MD steps
26+
nEQsteps = 10000 # Number of equilibration steps
27+
dt = 0.5 # Time step
28+
tdamp = 100.0 # Temperature damping
29+
30+
# Initialize MMFF_multi with the DA.mol2 file
31+
mmff.init(
32+
nSys_=nSys,
33+
xyz_name=xyz_name,
34+
sElementTypes="../../common_resources/ElementTypes.dat",
35+
sAtomTypes="../../common_resources/AtomTypes.dat",
36+
sBondTypes="../../common_resources/BondTypes.dat",
37+
sAngleTypes="../../common_resources/AngleTypes.dat",
38+
bMMFF=True,
39+
bEpairs=False,
40+
gamma = 1 / (dt * tdamp), # Temperature damping
41+
T = 300.0 # Temperature in Kelvin
42+
43+
)
44+
print("Initialization complete")
45+
46+
47+
48+
# Call computeFreeEnergy
49+
print(f"Computing free energy with parameters:")
50+
print(f" lamda1 = {lamda1}")
51+
print(f" lamda2 = {lamda2}")
52+
print(f" dc = {dc}")
53+
print(f" nbStep = {nbStep}")
54+
print(f" nMDsteps = {nMDsteps}")
55+
print(f" nEQsteps = {nEQsteps}")
56+
print(f" dt = {dt}")
57+
58+
result = mmff.computeFreeEnergy(
59+
lamda1=lamda1,
60+
lamda2=lamda2,
61+
dc=dc,
62+
nbStep=nbStep,
63+
nMDsteps=nMDsteps,
64+
nEQsteps=nEQsteps,
65+
dt=dt
66+
)
67+
68+
print(f"\nFree energy result: {result}")
69+
70+
if __name__ == "__main__":
71+
main()

‎pyBall/MMFF_multi.py‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -631,7 +631,15 @@ def plot_selection(sel=None,ax1=0,ax2=1,ps=None, s=100):
631631
asel=ps[sel]
632632
plt.scatter( asel[:,ax1], asel[:,ax2], s=s, facecolors='none', edgecolors='r' )
633633

634-
634+
# ========= Free Energy Calculation
635+
636+
# double computeFreeEnergy(double lamda1, double lamda2, int n, int* dc, int nbStep, int nMDsteps, int nEQsteps, double tdamp, double T, double dt)
637+
lib.computeFreeEnergy.argtypes = [c_double, c_double, c_int, c_int_p, c_int, c_int, c_int, c_double, c_double, c_double]
638+
lib.computeFreeEnergy.restype = c_double
639+
def computeFreeEnergy(lamda1, lamda2, dc, nbStep=100, nMDsteps=100000, nEQsteps=10000, tdamp=100.0, T=300, dt=0.5):
640+
n = len(dc)
641+
dc = np.array(dc, dtype=np.int32)
642+
return lib.computeFreeEnergy(lamda1, lamda2, n, _np_as(dc,c_int_p), nbStep, nMDsteps, nEQsteps, tdamp, T, dt)
635643

636644
# ====================================
637645
# ========= Test Functions

0 commit comments

Comments
 (0)