Package to perform free energy calculations with LAMMPS. Includes core-softed versions of Lennard-Jones and Coulomb potentials, a technique to avoid singularities in many free-energy calculation routes. This has been integrated in the LAMMPS code.
Laboratoire de Chimie
Dear Professor,
I am using your fdti.py tool to post process the result from the example CH4hyd, fdti01.
Below is your fdti.py file:
#!/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[1])
hderiv = float(sys.argv[2])
line = sys.stdin.readline()
while line.startswith("#"):
line = sys.stdin.readline()
v = 1.0
tok = line.split()
if len(tok) == 4:
v = float(tok[3])
lo = -rt * math.log(float(tok[2]) / v)
i = 0
sum = 0.0
for line in sys.stdin:
tok = line.split()
if len(tok) == 4:
v = float(tok[3])
hi = – rt * math.log(float(tok[2]) / 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