-
Notifications
You must be signed in to change notification settings - Fork 3
/
freeEn_correlated_data.py
89 lines (67 loc) · 2.77 KB
/
freeEn_correlated_data.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
from alchemlyb.parsing.gomc import extract_dHdl, extract_u_nk
from alchemlyb.estimators import MBAR, BAR, TI
import pandas as pd
import numpy as np
import os
bd = os.path.dirname(os.path.realpath(__file__))
temprature = 298 #temperature (K)
k_b = 1.9872036E-3 #kcal/mol/K
k_b_T = temprature * k_b
def get_delta(est):
""" Return the change in free energy and standard deviation for TI and MBAR estimators.
"""
delta = est.delta_f_.iloc[0, -1] * k_b_T
d_delta = est.d_delta_f_.iloc[0, -1] * k_b_T
return delta, d_delta
def get_delta2(est):
""" Return the change in free energy and standard deviation for BAR estimator.
"""
ee = 0.0
for i in range(len(est.d_delta_f_) - 1):
ee += est.d_delta_f_.values[i][i+1]**2
delta = est.delta_f_.iloc[0, -1] * k_b_T
d_delta = k_b_T * ee**0.5
return delta, d_delta
##################################################
numFile = 22
compounds = ["F3O", "F5O", "F7O", "F9O", "F11O", "F13O", "F15O", "F17O"]
fname = "Free_Energy_BOX_0_PRODUCTION_"
ext = ".dat"
NREP = 5
print("%16s, %16s, %16s, %16s, %16s, %16s, %16s" % ("#Compound","TI(kcal/mol)","stdev","MBAR(kcal/mol)","stdev","BAR(kcal/mol)","stdev"))
#loop through all compounds
for com in compounds:
tis = []
mbars = []
bars = []
for nr in range(NREP):
#read the free energy files
data_loc = bd + "/" + com +"/TI_"+str(nr+1)+"/state_"
files = []
for i in range(numFile + 1):
freeEn_file = fname + str(i) + ext
file_path = data_loc + str(i) + "/" + freeEn_file
print("%4s: Reading File: %s " % (com, file_path))
files.append(file_path)
#for TI estimator
#print("Working on TI method ...")
dHdl = pd.concat([extract_dHdl(f, T=temprature) for f in files])
ti = TI().fit(dHdl)
sum_ti, sum_ds_ti = get_delta(ti)
tis.append(sum_ti)
#for MBAR estimator
#print("Working on MBAR method ...")
u_nk = pd.concat([extract_u_nk(f, T=temprature) for f in files])
mbar = MBAR().fit(u_nk)
sum_mbar, sum_ds_mbar = get_delta(mbar)
mbars.append(sum_mbar)
#for BAR estimator
#print("Working on BAR method ...")
bar = BAR().fit(u_nk)
sum_bar, sum_ds_bar = get_delta2(bar)
bars.append(sum_bar)
print("Replica-%d, %4s, %7.4f, %7.4f, %7.4f, %7.4f, %7.4f, %7.4f" % (nr, com, sum_ti, sum_ds_ti, sum_mbar, sum_ds_mbar, sum_bar, sum_ds_bar))
tis = np.array(tis)
mbars = np.array(mbars)
bars = np.array(bars)
print("Average, %4s, %7.4f, %7.4f, %7.4f, %7.4f, %7.4f, %7.4f" % (com, np.average(tis), np.std(tis), np.average(mbars), np.std(mbars), np.average(bars), np.std(bars)))