-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathcheck_refeps.py
executable file
·100 lines (93 loc) · 3.98 KB
/
check_refeps.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
90
91
92
93
94
95
96
97
98
99
100
#!/usr/bin/env python
from __future__ import division
from mdoutremd import HRemMdout
import numpy as np
from optparse import OptionParser
from remd import HRemLog
import sys
C = '\033[91m'; EC = '\033[0m'
usage='%prog [Options] | [mdout1] [mdout2] [mdout3] ... [mdoutN]'
parser = OptionParser(usage=usage)
parser.add_option('-l', '--log', dest='input_log', metavar='FILE', default=None,
help='Input log file to parse')
parser.add_option('-d', '--deltas', dest='deltas', default=False,
action='store_true', help='Print out deltas')
opt, arg = parser.parse_args()
if opt.input_log is None:
if len(arg) == 0:
parser.print_help()
quit(1)
remlog = HRemMdout(arg)
else:
remlog = HRemLog(opt.input_log)
running_left_fes = np.zeros(len(remlog.reps))
running_right_fes = np.zeros(len(remlog.reps))
left_fes = [np.zeros(remlog.numexchg) for i in range(len(remlog.reps))]
right_fes = [np.zeros(remlog.numexchg) for i in range(len(remlog.reps))]
left_runs = [np.zeros(remlog.numexchg) for i in range(len(remlog.reps))]
right_runs = [np.zeros(remlog.numexchg) for i in range(len(remlog.reps))]
left_nexch = [[0 for i in range(remlog.numexchg)] for i in range(len(remlog.reps))]
right_nexch = [[0 for i in range(remlog.numexchg)] for i in range(len(remlog.reps))]
n_r_exch = [0 for i in range(len(remlog.reps))]
n_l_exch = [0 for i in range(len(remlog.reps))]
#fls = [open('python.00%d' % d, 'w') for d in range(len(remlog.reps))]
#print ' PotEne 1 PotEne 2'
for i, rep in enumerate(remlog.reps):
for j in range(remlog.numexchg):
nrep = remlog.reps[rep.neighbor_index[j]]
# fls[i].write('='*80 + '\n')
if (j % 2 != i % 2):
n_l_exch[i] += 1
running_left_fes[i] += np.exp((rep.potene1[j] - nrep.potene2[j])
* 503.01 / rep.temp[j])
left_fes[i][j] = 1 / 503.01 * rep.temp[j] * np.log(running_left_fes[i]
/ n_l_exch[i])
left_runs[i][j] = running_left_fes[i]
left_nexch[i][j] = n_l_exch[i]
try:
right_fes[i][j] = right_fes[i][j-1]
except IndexError:
pass
else:
n_r_exch[i] += 1
running_right_fes[i] += np.exp((rep.potene1[j] - nrep.potene2[j])
* 503.01 / rep.temp[j])
right_fes[i][j] = 1 / 503.01 * rep.temp[j] * np.log(running_right_fes[i] /
n_r_exch[i])
right_runs[i][j] = running_right_fes[i]
right_nexch[i][j] = n_r_exch[i]
try:
left_fes[i][j] = left_fes[i][j-1]
except IndexError:
pass
for j in range(remlog.numexchg):
print 'Exchange %8d' % (j+1)
for i in range(len(remlog.reps)):
d1 = abs(left_fes[i][j] - remlog.reps[i].left_fe[j])
d2 = abs(right_fes[i][j] - remlog.reps[i].right_fe[j])
if opt.deltas:
if d1 > 0.02 or d2 > 0.02:
b = C; e = EC
else:
b = e = ''
print (b
+ 'Replica %4d: Delta Left FE = %10.4f Delta Right FE = %10.4f'
+ e ) % (i, d1, d2)
else:
if d1 > 0.02 or d2 > 0.02:
b = C; e = EC
else:
b = e = ''
print (b
+ 'Replica %4d: Left FE = %10.4f (%10.4f) Right FE = %10.4f (%10.4f)'
+ e) % (i, left_fes[i][j], remlog.reps[i].left_fe[j], right_fes[i][j],
remlog.reps[i].right_fe[j])
# print 'Group 1: num_right_exchg = %10d %10d' % (right_nexch[i][0],
# right_nexch[i][1])
# print 'Group 1: num_left_exchg = %10d %10d' % (left_nexch[i][0],
# left_nexch[i][1])
# print 'Group 1: total_right_fe = %s %s' % (right_runs[i][0],
# right_runs[i][1])
# print 'Group 1: total_left_fe = %s %s' % (left_runs[i][0],
# left_runs[i][1])
print '='*80