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