1import math 2import sys 3 4# read FILE with CVs and weights 5FILENAME_ = sys.argv[1] 6# number of CVs for FES 7NCV_ = int(sys.argv[2]) 8# read minimum, maximum and number of bins for FES grid 9gmin = []; gmax = []; nbin = [] 10for i in range(0, NCV_): 11 i0 = 3*i + 3 12 gmin.append(float(sys.argv[i0])) 13 gmax.append(float(sys.argv[i0+1])) 14 nbin.append(int(sys.argv[i0+2])) 15# read KBT_ 16KBT_ = float(sys.argv[3*NCV_+3]) 17OUTF_ = sys.argv[3*NCV_+4] 18 19def get_indexes_from_index(index, nbin): 20 indexes = [] 21 # get first index 22 indexes.append(index%nbin[0]) 23 # loop 24 kk = index 25 for i in range(1, len(nbin)-1): 26 kk = ( kk - indexes[i-1] ) / nbin[i-1] 27 indexes.append(kk%nbin[i]) 28 if(len(nbin)>=2): 29 indexes.append( ( kk - indexes[len(nbin)-2] ) / nbin[len(nbin) -2] ) 30 return indexes 31 32def get_indexes_from_cvs(cvs, gmin, dx): 33 keys = [] 34 for i in range(0, len(cvs)): 35 keys.append(int( round( ( cvs[i] - gmin[i] ) / dx[i] ) )) 36 return tuple(keys) 37 38def get_points(key, gmin, dx): 39 xs = [] 40 for i in range(0, len(key)): 41 xs.append(gmin[i] + float(key[i]) * dx[i]) 42 return xs 43 44# define bin size 45dx = [] 46for i in range(0, NCV_): 47 dx.append( (gmax[i]-gmin[i])/float(nbin[i]-1) ) 48 49# create histogram 50histo = {} 51 52# read file and fill in histogram 53for lines in open(FILENAME_, "r").readlines(): 54 riga = lines.strip().split() 55 # check format 56 if(len(riga)!=NCV_+1 and len(riga)!=NCV_+2): 57 print(FILENAME_,"is in the wrong format!") 58 exit() 59 # read CVs 60 cvs = [] 61 for i in range(1, NCV_+1): 62 cvs.append(float(riga[i])) 63 # get indexes 64 key = get_indexes_from_cvs(cvs, gmin, dx) 65 if(len(riga)==NCV_+2): 66 # read weight 67 w = float(riga[NCV_+1]) 68 else: w = 1.0 69 # update histogram 70 if key in histo: histo[key] += w 71 else: histo[key] = w 72 73# calculate free-energy and minimum value 74min_fes = 1.0e+15 75for key in histo: 76 histo[key] = -KBT_ * math.log(histo[key]) 77 if(histo[key] < min_fes): min_fes = histo[key] 78 79# total numbers of bins 80nbins = 1 81for i in range(0, len(nbin)): nbins *= nbin[i] 82 83# print out FES 84log = open(OUTF_, "w") 85# this is needed to add a blank line 86xs_old = [] 87for i in range(0, nbins): 88 # get the indexes in the multi-dimensional grid 89 key = tuple(get_indexes_from_index(i, nbin)) 90 # get CV values for that grid point 91 xs = get_points(key, gmin, dx) 92 # add a blank line for gnuplot 93 if(i == 0): 94 xs_old = xs[:] 95 else: 96 flag = 0 97 for j in range(1,len(xs)): 98 if(xs[j] != xs_old[j]): 99 flag = 1 100 xs_old = xs[:] 101 if (flag == 1): log.write("\n") 102 # print value of CVs 103 for x in xs: 104 log.write("%12.6lf " % x) 105 # print FES 106 if key in histo: 107 fes = histo[key]-min_fes 108 log.write(" %12.6lf\n" % fes) 109 else: 110 log.write(" Infinity\n") 111log.close() 112