Scripts:mkwt.py

From UC_Chemistry

This program takes a 'biasing energy file' and computes weights -- c*exp(-W/kT). It can optionally exclude frames (assigning a weight of zero), specify a column other than -1 (=last) for the weight and/or compute W = en-d by subtracting a second file from the first.

Note: Make sure kT is correct for your system!

Source

#!/usr/bin/python

import sys
from ucgrad import *

#beta = 1./(Boltzmann*298.15)
beta = 4.184/(Boltzmann*298.15)

def main(argv):
        if len(argv) < 3:
                print "Usage: %s [-i <indicator file>] [-n wt energy col.] [-d diff. file (W = en - d)] <energy file> <weight>"%(argv[0])
                return 1

        i = 1
        n = -1
        ind = None
        diff = 0
        while i < len(argv):
                if argv[i] == '-n':
                        n = int(argv[i+1])
                        i += 2
                elif argv[i] == '-i': # flush violators to zero
                        ind = 1-(read_matrix(argv[i+1])[:,-1] > 1.0e-10)
                        i += 2
                elif argv[i] == '-d':
                        diff = read_matrix(argv[i+1])[:,-1]
                        i += 2
                else:
                        break

        M = read_matrix(argv[i]) # Future output file.
        if ind == None:
                ind = ones(len(M))

        M[:,1] = ind * exp((M[:,n]-diff)*beta)
        # re-normalize to make ea. sample count at most once.
        M[:,1] /= max(M[:,1])
        write_matrix(argv[i+1], M[:,:2])
        print "Effective/qualifying/total number of samples = %f/%d/%d"%(\
                        sum(M[:,1]), sum(ind), len(ind))

        return 0

if __name__=="__main__":
        main(sys.argv)