One Reply to “compute fep for LAMMPS”

  1. 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

Leave a Reply

Your email address will not be published.