## One Reply to “compute fep for LAMMPS”

1. Dezhao Huang says:

Dear Professor,
I am using your fdti.py tool to post process the result from the example CH4hyd, fdti01.
#!/usr/bin/env python
# fdti.py – integrate compute fep results using the trapezoidal rule

import sys, math

if len(sys.argv) < 3:
print ("Finite Difference Thermodynamic Integration (Mezei 1987)")
print ("Trapezoidal integration of compute_fep results at equally-spaced points")
print ("usage: fdti.py temperature hderiv < fep.lmp")
sys.exit()

rt = 0.008314 / 4.184 * float(sys.argv)
hderiv = float(sys.argv)

while line.startswith("#"):

v = 1.0
tok = line.split()
if len(tok) == 4:
v = float(tok)
lo = -rt * math.log(float(tok) / v)

i = 0
sum = 0.0
for line in sys.stdin:
tok = line.split()
if len(tok) == 4:
v = float(tok)
hi = – rt * math.log(float(tok) / v)
sum += (hi + lo) / (2 * hderiv)
lo = hi
i += 1

print (sum / i) # int_0^1: divide by i == multiply by delta

By comparing your script with the lammps compute FEP webpage, I could not figure out why you have hderiv here sum += (hi + lo) / (2 * hderiv) to calculate the area. Although I know hderiv is the charge perturbation value in the simulation.
Can you further elaborate ?
I would be really appreciated!
Best,
******************************************************************
Dezhao Huang
Ph.D. candidate at University of Notre Dame
MOlecular/Nano-Sacle Transport & Energy Research (MONSTER) Laboratory
Department of Aerospace and Mechanical Engineering
B14 Fitzpatrick Hall
Notre Dame, IN 46556-5710
Tel: (573)-578-0869
dhuang2@nd.edu