hdiff output

r27479/compute_rates.py 2016-03-16 18:33:46.587194776 +0000 r27478/compute_rates.py 2016-03-16 18:33:46.807197039 +0000
177:     177:     
178:     return generator.rate_constants, generator.Peq, generator.rate_norm178:     return generator.rate_constants, generator.Peq, generator.rate_norm
179: 179: 
180: 180: 
181: description="""Compute rates from a pathsample database.  The transition181: description="""Compute rates from a pathsample database.  The transition
182: states and minima data are read from min.data and ts.data.  The product182: states and minima data are read from min.data and ts.data.  The product
183: and reactant states are read from min.A and min.B.  The rates are computed183: and reactant states are read from min.A and min.B.  The rates are computed
184: by solving a set of linear equations using a sparse solver.184: by solving a set of linear equations using a sparse solver.
185: """185: """
186: 186: 
187: def run(directory, T, A, B, out_prefix, tstart, reverse=False):187: def run(directory, T, A, B, out_prefix, tstart):
188:     if not reverse: 
189:         source = "A" 
190:         destination = "B" 
191:         direction = "A->B" 
192:     else: 
193:         source = "B" 
194:         destination = "A" 
195:         direction = "B->A" 
196:  
197:     rate_constants, Peq, knorm = make_rates(directory, T)188:     rate_constants, Peq, knorm = make_rates(directory, T)
198:             189:             
199: 190: 
200:     if True:191:     if True:
201:         fname = "{}.rate_consts".format(out_prefix)192:         fname = "{}.rate_consts".format(out_prefix)
202:         print "saving rate constants to:", fname193:         print "saving rate constants to:", fname
203:         with open(fname, "w") as fout:194:         with open(fname, "w") as fout:
204:             fout.write("#starting_minimum ending_minimum rate_constant\n") 
205:             for (u,v), k in sorted(rate_constants.iteritems()):195:             for (u,v), k in sorted(rate_constants.iteritems()):
206:                 fout.write("{} {} {}\n".format(u, v, k/knorm))196:                 fout.write("%6d %6d %s\n" % (u,v,k/knorm))
207: 197: 
208:     print "checking and reducing the graph structure"198:     print "checking and reducing the graph structure"
209:     try:199:     try:
210:         rate_constants = reduce_rates(rate_constants, B, A=A)200:         rate_constants = reduce_rates(rate_constants, B, A=A)
211:     except Exception:201:     except Exception:
212:         analyze_graph_error(rate_constants, A, B)202:         analyze_graph_error(rate_constants, A, B)
213:         raise203:         raise
214:     204:     
215:     print "computing mean first passage times"205:     print "computing mean first passage times"
216:     calculator = TwoStateRates(rate_constants, A, B, weights=Peq, check_rates=False)206:     calculator = TwoStateRates(rate_constants, A, B, weights=Peq, check_rates=False)
217:     calculator.compute_rates(use_umfpack=True)207:     calculator.compute_rates(use_umfpack=True)
218:     print "computing rates"208:     print "computing rates"
219:     kAB = calculator.get_rate_AB()209:     kAB = calculator.get_rate_AB()
220:     print "k({}) {}".format(direction, kAB / knorm)210:     print "k(B<-A)", kAB / knorm
221:     211:     
222:     if True:212:     if True:
223:         fname = "{}.rates".format(out_prefix)213:         fname = "{}.mfpt".format(out_prefix)
224:         print "saving rates and mean first passage times for all minima to reach {} to file {}".format(destination, fname)214:         print "saving mean first passage times for all minima to reach B to file", fname
225:         mfpt = sorted([(m, t) for m, t in 215:         mfpt = sorted([(m, t) for m, t in 
226:                        calculator.mfpt_computer.mfpt_dict.iteritems()])216:                        calculator.mfpt_computer.mfpt_dict.iteritems()])
227:         with open(fname, "w") as fout:217:         with open(fname, "w") as fout:
228:             fout.write("#Rates and mean first passage times for each node to reach {}\n".format(destination)) 
229:             fout.write("#The rate is just the inverse of the mean first passage time\n") 
230:             fout.write("#minimum_index rate mean_first_passage_time\n") 
231:             for i, t in mfpt:218:             for i, t in mfpt:
232:                 mt = t * knorm219:                 fout.write("%d %.12g\n" % (i, t * knorm))
233:                 fout.write("{index} {rate} {mfpt}\n".format(index=i, rate=1./mt, 
234:                                                             mfpt=mt)) 
235: 220: 
236:     print "computing committor probabilities"221:     print "computing committor probabilities"
237:     sys.stdout.flush()222:     sys.stdout.flush()
238:     calculator.compute_committors()223:     calculator.compute_committors()
239:     print "computing steady state rate"224:     print "computing steady state rates"
240:     kSS = calculator.get_rate_AB_SS() / knorm225:     kSS = calculator.get_rate_AB_SS() / knorm
241:     print "kSS({}) {}".format(direction, kSS)226:     print "kSS(B<-A)", kSS
242:     227:     
243:     # print the committors228:     if True:
244:     # get the alternate definition of committors for the nodes in A and B229:         # get the alternate definition of committors for the nodes in A and B
245:     Acomm = calculator.get_alternate_committors(A, B)230:         Acomm = calculator.get_alternate_committors(A, B)
246:     Bcomm = calculator.get_alternate_committors(B, A)231:         Bcomm = calculator.get_alternate_committors(B, A)
247:     all_committors = calculator.committor_computer.committor_dict.copy()232:         all_committors = calculator.committor_computer.committor_dict.copy()
248:     all_committors.update(Acomm)233:         all_committors.update(Acomm)
249:     all_committors.update(Bcomm)234:         all_committors.update(Bcomm)
250:     235:         
251:     fname = "{}.committors".format(out_prefix)236:         fname = "{}.committors".format(out_prefix)
252:     print "saving committor probabilities for all minima to file", fname237:         print "saving committor probabilities for all minima to file", fname
253:     coms = sorted([(m, c) for m, c in all_committors.iteritems()])238:         coms = sorted([(m, c) for m, c in all_committors.iteritems()])
254:     with open(fname, "w") as fout:239:         with open(fname, "w") as fout:
255:         fout.write("#probability a trajectory starting from each node ends up "240:             for i, c in coms:
256:                    "in {B} before returning to {A}\n".format(B=destination, A=source))241:                 fout.write("%d %.12g\n" % (i, c))
257:         fout.write("#minimum_index committor_probability\n") 
258:         for i, c in coms: 
259:             fout.write("{} {}\n".format(i, c)) 
260:     242:     
261:     time_solve = calculator.mfpt_computer.time_solve + calculator.committor_computer.time_solve243:     time_solve = calculator.mfpt_computer.time_solve + calculator.committor_computer.time_solve
262:     print "time spent solving linear equations", time_solve, "seconds"244:     print "time spent solving linear equations", time_solve, "seconds"
263:     print "total time", time.clock() - tstart245:     print "total time", time.clock() - tstart
264: 246: 
265: 247: 
266: def main():248: def main():
267:     """the main loop"""249:     """the main loop"""
268:     tstart =  time.clock()250:     tstart =  time.clock()
269:     parser = argparse.ArgumentParser(description=description)251:     parser = argparse.ArgumentParser(description=description)
270: 252: 
271:     parser.add_argument("-d", type=str, default=".",253:     parser.add_argument("-d", type=str, default=".",
272:                         help="directory with the min.data, ts.data files")254:                         help="directory with the min.data, ts.data files")
273:     parser.add_argument("-T", type=float, default=1., help="Temperature")255:     parser.add_argument("-T", type=float, default=1., help="Temperature")
274:     parser.add_argument("--prefix", type=str, default="out", help="output prefix")256:     parser.add_argument("--prefix", type=str, default="rates.out", help="output prefix")
275:     args = parser.parse_args()257:     args = parser.parse_args()
276:     258:     
277:     print "temperature", args.T259:     print "temperature", args.T
278:     directory = args.d260:     directory = args.d
279:     print "reading from directory:", os.path.abspath(directory)261:     print "reading from directory:", os.path.abspath(directory)
280:     out_prefix = args.prefix262:     out_prefix = args.prefix
281:     263:     
282:     A = read_minA(directory+"/min.A")264:     A = read_minA(directory+"/min.A")
283:     B = read_minA(directory+"/min.B")265:     B = read_minA(directory+"/min.B")
284:     266:     
285:     print "\ncomputing rates from A to B"267:     print "\ncomputing rates from A to B"
286:     run(directory, args.T, A, B, out_prefix + ".AtoB", tstart, reverse=False)268:     run(directory, args.T, A, B, out_prefix + "AtoB", tstart)
287:     print "\ncomputing rates from B to A"269:     print "\ncomputing rates from B to A"
288:     run(directory, args.T, B, A, out_prefix + ".BtoA", tstart, reverse=True)270:     run(directory, args.T, B, A, out_prefix + "BtoA", tstart)
289:     271:     
290:     272:     
291:      273:      
292: 274: 
293: 275: 
294: if __name__ == "__main__":276: if __name__ == "__main__":
295:     main()277:     main()


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0