#!/usr/bin/env python # fep.py - calculate free energy from compute fep results import sys import math if len(sys.argv) < 2: print("Free Energy Perturbation") print("usage: fep.py temperature < out.fep") sys.exit() rt = 0.008314 / 4.184 * float(sys.argv[1]) v = 1.0 sum = 0.0 for line in sys.stdin: if line.startswith("#"): continue tok = line.split() if len(tok) == 4: v = float(tok[3]) sum += math.log(float(tok[2]) / v) print(-rt * sum)