Mercurial > repos > chemteam > biomd_extract_clusters

comparison NEQGamma.py @ 1:078dfd7fb26d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit 79589d149a8ff2791d4f71d28b155011672db827"
author chemteam
date Fri, 11 Sep 2020 21:55:13 +0000
parents
children e0ecaf2d05fb
comparison
equal deleted inserted replaced
0:ad49025ba90d 1:078dfd7fb26d
1 #!/usr/bin/env python
2
3 # coding: utf-8
4 # This script is a modified version of a script written
5 # by Steffen Wolf under the GPL v3.0.
6 # The original version can be accessed at
7 # https://github.com/moldyn/dcTMD/blob/master/NEQGamma.py
8
9 import argparse
10 import json
11 import sys
12
13 import numpy as np
14
15 import pandas as pd
16
17 import scipy
18 import scipy.integrate
19 from scipy.ndimage.filters import gaussian_filter
20
21
22 def get_file_names(list_file):
23 with open(list_file) as f:
24 return [line for line in f.read().split('\n') if line]
25
26
27 def NEQGamma(file_names, output, output_frict, vel, T, av, sigma):
28 N = len(file_names)
29 RT = 0.0083144598 * T
30
31 sys.stdout.write("reading data...\n")
32
33 # read in initial data to get length of necessary array
34 test_file = pd.read_csv(file_names[0], delim_whitespace=True,
35 header=None, skiprows=17, dtype=float)
36 length_data = len(test_file[0].values)
37 full_force_set = np.zeros((N, length_data))
38 x = np.zeros(length_data)
39 t = np.zeros(length_data)
40 t = test_file[0].values
41 x = test_file[0].values * vel
42
43 # read in data
44 for i in range(0, N):
45 current_file_name = file_names[i]
46 sys.stdout.write("reading file {}\n".format(current_file_name))
47 input_file_data = pd.read_csv(current_file_name, delim_whitespace=True,
48 header=None, skiprows=17, dtype=float)
49 full_force_set[i, :] = input_file_data[1].values
50
51 # preprocessing
52 # * force average: calculate $\left< f_c (t)\right>_N$.
53 # **Important:** this is an ensemble average over the trajectory ensemble
54 # $N$, not the time average over $t$
55 av_force = np.zeros(length_data)
56 av_forceintegral = np.zeros(length_data)
57 for i in range(length_data):
58 av_force[i] = np.mean(full_force_set[:, i])
59 av_forceintegral[1:] = scipy.integrate.cumtrapz(av_force, x)
60
61 # calculate $\delta f_c(t) = f_c(t) - \left< f_c (t) \right>_N$ for all $t$
62 sys.stdout.write("calculating fluctuations...\n")
63 delta_force_set = np.zeros((N, length_data))
64 for i in range(length_data):
65 delta_force_set[:, i] = full_force_set[:, i] - av_force[i]
66
67 # evaluation
68 # * optimized algorithm for numerical evaluation:
69 # * integrate: $\int_0^t dt' \delta f_c(t')$ for all $t'$
70 # * multiply by $\delta f_c(t)$ to yield
71 # $\int_0^t dt'\delta f_c(t) \delta f_c(t')$ for $t$
72 # with all $t' \leq t$ each
73 # * then calculate the ensemble average
74 # $\left< \int_0^t dt' \delta f_c(t) \delta f_c(t') \right>$
75 int_delta_force_set = np.zeros((N, length_data))
76 for n in range(N):
77 int_delta_force_set[n, 1:] = scipy.integrate.cumtrapz(
78 delta_force_set[n, :], t)
79
80 sys.stdout.write("averaging and integrating...\n")
81 intcorr = np.zeros((N, length_data))
82
83 for n in range(N):
84 for i in range(length_data):
85 intcorr[n, i] = delta_force_set[n, i] * int_delta_force_set[n, i]
86 if i % 1000 == 0:
87 sys.stdout.write("Trajectory {:2d} {:3.1f} % done\r".format(
88 n + 1, (i / length_data) * 100))
89
90 # shape of $\int_0^t dt' \delta f_c(t) \delta f_c(t')$:
91 sys.stdout.write("final average...\n")
92 av_intcorr = np.zeros(length_data)
93 for i in range(length_data):
94 av_intcorr[i] = np.mean(intcorr[:, i]) / RT
95
96 # autocorrelation function evaluation:
97 # * calculate $\left< \delta f_c(t) \delta f_c(t') \right>$
98 # for the last $t$
99
100 corr_set = np.zeros((N, length_data))
101 autocorr_set = np.zeros(length_data)
102
103 sys.stdout.write("calculating and processing ACF...\n")
104 for n in range(N):
105 for i in range(length_data):
106 corr_set[n, i] = delta_force_set[
107 n, i] * delta_force_set[n, length_data - 1]
108
109 for i in range(length_data):
110 autocorr_set[i] = np.mean(corr_set[:, i])
111
112 # * Gauss filter:
113 sys.stdout.write("applying Gauss filter...\n")
114 blurr = sigma
115 blurred = np.zeros(length_data)
116 blurred = gaussian_filter(av_intcorr, sigma=blurr)
117
118 # * sliding window average:
119 sys.stdout.write("applying sliding window average...\n")
120 window = av
121 runn_av = np.zeros(length_data)
122 runn_av = np.convolve(av_intcorr, np.ones((window,)) / window, mode='same')
123
124 # * $W_{diss}$ from integration:
125 wdiss = np.zeros(length_data)
126 wdiss[1:] = scipy.integrate.cumtrapz(av_intcorr, x) * vel
127
128 sys.stdout.write("writing output...\n")
129 dist = open(output, "w")
130 frict = open(output_frict, "w")
131
132 dist.write(
133 "#x force_integral frict_coeff wdiss corrected_force_integral\n")
134 for i in range(length_data):
135 dist.write("{:15.8f} {:20.8f} {:20.8f} {:20.8f} {:20.8f}\n".format(
136 x[i], av_forceintegral[i], av_intcorr[i], wdiss[i],
137 av_forceintegral[i] - wdiss[i]))
138
139 frict.write("""#x ACF frict_coeff """
140 """gauss_filtered_frict_coeff av_window_frict_coeff\n""")
141 for i in range(length_data):
142 frict.write("{:15.8f} {:20.8f} {:20.8f} {:20.8f} {:20.8f}\n".format(
143 x[i], autocorr_set[i], av_intcorr[i], blurred[i], runn_av[i]))
144
145 dist.close()
146 frict.close()
147
148 sys.stdout.write("Done!\n")
149
150
151 def main():
152 parser = argparse.ArgumentParser(description="""dcTMD friciton correction
153 (please cite: Wolf, S., Stock, G. Targeted Molecular Dynamics
154 Calculations of Free Energy Profiles Using a Nonequilibrium
155 Friction Correction. J. Chem. Theory Comput. 2018, 14(12), 6175-6182,
156 DOI: 10.1021/acs.jctc.8b00835). Integrates a constraint force file via
157 trapezoid rule, calculates the NEQ memory friction kernel and friction
158 factors, and performs a friction correction. First column: reaction
159 coordinate in nm calculated via t * vel. Second column: force integral,
160 i.e. the work profile. Third column: friction factors. Fourth column:
161 trapezoid integral (final value) of friction work along reaction
162 coordinate. Fourth column: friction corrected work profile. ATTENTION:
163 Use with python3 or higher!""")
164 parser.add_argument('-i', metavar='<xvg force file>', type=str,
165 help="""List of xvg constraint force files prefix
166 as given by Gromacs mdrun -pf option before running
167 number.""")
168 parser.add_argument('-o', metavar='<combined results>', type=str,
169 help="""file to write x, dG(x), friction coefficeint by
170 integration (time resolved), and the friction-corrected
171 dG(x).""")
172 parser.add_argument('-ofrict', metavar='<combined friction results>',
173 type=str,
174 help="""file to write x, ACF, friction coefficeint by
175 integration (time resolved), gauss filtered friction
176 coefficient, and slide window averaged friction.""")
177 parser.add_argument('-vel', metavar='<pull velocity>', type=float,
178 help="""pull velocity in nm/ps for converting simulation
179 time t into distance x""")
180 parser.add_argument('-T', metavar='<temperature>', type=float,
181 help='temperature in K')
182 parser.add_argument('-av', metavar='<average window>', type=int,
183 help="""size of averaging window for displaying
184 Gamma(x) (recommended: 4 to 20 per 100 data points)""")
185 parser.add_argument('-sigma', metavar='<gauss blurr>', type=int,
186 help="""sigma value for Gauss filter for displaying
187 Gamma(x) (recommended: 4 per 100 data points)""")
188 parser.add_argument('-json', metavar='<json>', type=str,
189 help='JSON file defining cluster membership')
190
191 args = parser.parse_args()
192
193 file_names = get_file_names(args.i)
194 if args.json:
195 with open(args.json) as f:
196 j = json.load(f)
197 file_names_dct = {n: [file_names[int(m)] for m in j[n]] for n in j}
198
199 for cluster in file_names_dct:
200 NEQGamma(file_names_dct[cluster],
201 'cluster{}_{}'.format(cluster, args.o),
202 'cluster{}_{}'.format(cluster, args.ofrict),
203 args.vel, args.T, args.av, args.sigma)
204 else:
205 NEQGamma(file_names, args.o, args.ofrict, args.vel,
206 args.T, args.av, args.sigma)
207
208
209 if __name__ == "__main__":
210 main()