-
Notifications
You must be signed in to change notification settings - Fork 1
/
bedamtempt_async_re.py
219 lines (189 loc) · 8.52 KB
/
bedamtempt_async_re.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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
import sys, time, random, math
from pj_async_re import async_re_job
from bedam_async_re import bedam_async_re_job
class bedamtempt_async_re_job(bedam_async_re_job):
def _checkInput(self):
async_re_job._checkInput(self)
#make sure BEDAM + TEMPERATURE is wanted
if self.keywords.get('RE_TYPE') != 'BEDAMTEMPT':
self._exit("RE_TYPE is not BEDAMTEMPT")
#BEDAM runs with IMPACT
if self.keywords.get('ENGINE') != 'IMPACT':
self._exit("ENGINE is not IMPACT")
#input files
self.extfiles = self.keywords.get('ENGINE_INPUT_EXTFILES')
if not (self.extfiles is None):
if self.extfiles != '':
self.extfiles = self.extfiles.split(',')
#list of lambdas
if self.keywords.get('LAMBDAS') is None:
self._exit("LAMBDAS needs to be specified")
lambdas = self.keywords.get('LAMBDAS').split(',')
#list of temperatures
if self.keywords.get('TEMPERATURES') is None:
self._exit("TEMPERATURES needs to be specified")
temperatures = self.keywords.get('TEMPERATURES').split(',')
#build parameters for the lambda/temperatures combined states
self.nreplicas = self._buildBEDAMStates(lambdas,temperatures)
def _buildBEDAMStates(self,lambdas,temperatures):
self.stateparams = []
for lambd in lambdas:
for tempt in temperatures:
st = {}
st['lambda'] = lambd
st['temperature'] = tempt
self.stateparams.append(st)
return len(self.stateparams)
def _buildInpFile(self, replica):
"""
Builds input file for a BEDAM replica based on template input file
BASENAME.inp for the specified replica at lambda=lambda[stateid] for the
specified cycle.
"""
basename = self.basename
stateid = self.status[replica]['stateid_current']
cycle = self.status[replica]['cycle_current']
template = "%s.inp" % basename
inpfile = "r%d/%s_%d.inp" % (replica, basename, cycle)
lambd = self.stateparams[stateid]['lambda']
temperature = self.stateparams[stateid]['temperature']
# read template buffer
tfile = self._openfile(template, "r")
tbuffer = tfile.read()
tfile.close()
# make modifications
tbuffer = tbuffer.replace("@n@",str(cycle))
tbuffer = tbuffer.replace("@nm1@",str(cycle-1))
tbuffer = tbuffer.replace("@lambda@",lambd)
tbuffer = tbuffer.replace("@temperature@",temperature)
# write out
ofile = self._openfile(inpfile, "w")
ofile.write(tbuffer)
ofile.close()
# update the history status file
ofile = self._openfile("r%d/state.history" % replica, "a")
ofile.write("%d %d %s %s\n" % (cycle, stateid, lambd, temperature))
ofile.close()
def _doExchange_pair(self,repl_a,repl_b):
"""
Performs exchange of lambdas for BEDAM replica exchange.
"""
kb = 0.0019872041
cycle_a = self.status[repl_a]['cycle_current']
sid_a = self.status[repl_a]['stateid_current']
lambda_a = float(self.stateparams[sid_a]['lambda'])
temperature_a = float(self.stateparams[sid_a]['temperature'])
# u_a: binding energy of replica a
# h_a: total energy of replica a (includes kinetic energy as we are not
# doing velocity rescaling here)
(u_a,h_a) = self._extractLast_BindingEnergy_TotalEnergy(repl_a,cycle_a)
(u_a,h_a) = (float(u_a),float(h_a))
cycle_b = self.status[repl_b]['cycle_current']
sid_b = self.status[repl_b]['stateid_current']
lambda_b = float(self.stateparams[sid_b]['lambda'])
temperature_b = float(self.stateparams[sid_b]['temperature'])
(u_b,h_b) = self._extractLast_BindingEnergy_TotalEnergy(repl_b,cycle_b)
(u_b,h_b) = (float(u_b),float(h_b))
# Acceptance criterion is based on exp(-Delta) where
# Delta = -(beta_b - beta_a)*[H_b-H_a] - (lmbd_b - lmbd_a)[beta_a*u_b-beta_b*u_a]
# To derive this start from the Boltzmann weight exp[-F(x|lambda,beta)] where
# F(x|lambda,beta) = beta*[H_0(x)+lambda*u(x)], set up the usual Metropolis
# exchange rules noticing that:
# H_b(x_a) = H_a(x_a) + (lmbd_b - lmbd_a)*u(x_a)
#
beta_a = 1./(kb*temperature_a)
beta_b = 1./(kb*temperature_b)
dl = lambda_b - lambda_a
du = beta_a*u_b - beta_b*u_a
db = beta_b - beta_a
dh = h_b - h_a
delta = -dl*du - db*dh
if self.keywords.get('VERBOSE') == "yes":
print "Pair Info"
print "%d %f %f %f %f" % (repl_a, lambda_a, u_a, beta_a, h_a)
print "%d %f %f %f %f" % (repl_b, lambda_b, u_b, beta_b, h_b)
print "dl = %f du = %f dh = %f delta = %f" % (dl,du,dh,delta)
csi = random.random()
if math.exp(-delta) > csi:
if self.keywords.get('VERBOSE') == "yes":
print "Accepted %f %f" % (math.exp(-delta),csi)
print (self.status[repl_a]['stateid_current'], self.status[repl_b]['stateid_current'])
self.status[repl_a]['stateid_current'] = sid_b
self.status[repl_b]['stateid_current'] = sid_a
if self.keywords.get('VERBOSE') == "yes":
print (self.status[repl_a]['stateid_current'], self.status[repl_b]['stateid_current'])
else:
if self.keywords.get('VERBOSE') == "yes":
print "Rejected %f %f" % (math.exp(-delta),csi)
def _extractLast_lambda_BindingEnergy_TotalEnergy(self,repl,cycle):
"""
Extracts binding energy from Impact output
"""
output_file = "r%s/%s_%d.out" % (repl,self.basename,cycle)
datai = self._getImpactData(output_file)
nf = len(datai[0])
nr = len(datai)
# [nr-1]: last record
# [nf-2]: lambda (next to last item)
# [nf-1]: binding energy (last item)
# [2]: total energy item (0 is step number and 1 is temperature)
#
# (lambda, binding energy, total energy)
return (datai[nr-1][nf-2],datai[nr-1][nf-1],datai[nr-1][2])
def print_status(self):
"""
Writes to BASENAME_stat.txt a text version of the status of the RE job
It's fun to follow the progress in real time by doing:
watch cat BASENAME_stat.txt
"""
logfile = "%s_stat.txt" % self.basename
ofile = self._openfile(logfile,"w")
log = "Replica State Lambda Temperature Status Cycle \n"
for k in range(self.nreplicas):
stateid = self.status[k]['stateid_current']
log += "%6d %5d %s %s %5s %5d \n" % (k, stateid, self.stateparams[stateid]['lambda'], self.stateparams[stateid]['temperature'], self.status[k]['running_status'], self.status[k]['cycle_current'])
log += "Running = %d\n" % self.running
log += "Waiting = %d\n" % self.waiting
ofile.write(log)
ofile.close()
def _getPot(self,repl,cycle):
(lmb, u, etot) = self._extractLast_lambda_BindingEnergy_TotalEnergy(repl,cycle)
# removes lambda*u from etot to get e0. Note that this is the lambda from the
# output file not the current lambda.
e0 = float(etot) - float(lmb)*float(u)
return (e0,float(u))
def _getPar(self,repl):
sid = self.status[repl]['stateid_current']
lmb = float(self.stateparams[sid]['lambda'])
tempt = float(self.stateparams[sid]['temperature'])
kb = 0.0019872041
beta = 1./(kb*tempt)
return (beta,lmb)
def _reduced_energy(self,par,pot):
# par: list of parameters
# pot: list of potentials
# This is for temperature/binding potential beta*(U0+lambda*u)
beta = par[0]
lmb = par[1]
e0 = pot[0]
u = pot[1]
return beta*(e0 + lmb*u)
if __name__ == '__main__':
# Parse arguments:
usage = "%prog <ConfigFile>"
if len(sys.argv) != 2:
print "Please specify ONE input file"
sys.exit(1)
commandFile = sys.argv[1]
print ""
print "===================================="
print "BEDAM Asynchronous Replica Exchange "
print "===================================="
print ""
print "Started at: " + str(time.asctime())
print "Input file:", commandFile
print ""
sys.stdout.flush()
rx = bedamtempt_async_re_job(commandFile, options=None)
rx.setupJob()
rx.scheduleJobs()