compute/fep dlambda

Hello LAMMPS developers/users and Professor Padua,

I would like to use the compute/fep module written for LAMMPS. First of all, please let me thank you for spending your time to write such nice tools.

I have one question about the fdti examples provided in LAMMPS.

In the “CH4hyd” or “CC-CO” example, the change frequency in the fix adapt/fep line is 125000 steps and the total number of steps is 2500000, so each change is 5% or 0.05, but you set variable dlambda to equal 0.002 in these fdti examples.

What am I missing here?

Actually, the only major difference between the fdti and the fep examples is dlambda, the change frequency in the fix adapt/fep line was set to be the same for both cases.



Thank you for using the fep code. LAMMPS is “materials” oriented and it’s such a great MD code I though could benefit of some more “chemical” tools.

When calculating a free energy difference using the FEP method, you perform a perturbation at the value \lambda_{i+1}, on a trajectory generated at \lambda_i (the activation parameter that defines your free energy route from the initial to final systems) Therefore the size of the \delta \lambda (the perturbation) and the step in the fix adapt/fep must match in the FEP method. This is so that the next trajectory is acquired at \lambda_{i+1} and then the perturbation will be performed for \lambda_{i+2}, and so on.

If you are using FDTI, then the idea is to use the \delta \lambda to obtain a numerical derivative of the free energy at a certain value of \lambda_{i}. In this case the size of the “perturbation" is not connected with the size of the step in the fix adapt/fep.
This \delta \lambda has to be sufficiently small to enable a good estimate of the derivative, but not too small otherwise we lose precision. I just used 0.002 as an example because this is what worked for me in certain applications.

In both types of calculation, the fix adapt defines for how long your trajectory is generated at a certain value of \lambda. You should allow for some initial equilibration period after each “step”. The length of the panels I used are just examples. You will certainly need longer ones for serious work (depending on your problem).

In my opinion FEP is a superior method if you’ll be using this code. Because the derivative of the Hamiltonian is not obtained analytically, and because it is so easy to go by small steps, I think FEP will be better. But the BAR method is a very efficient one, although more laborious to apply in multiple steps using this code.

Feel free to ask more questions.


Agilio Padua

Thank you for detailed reply. The explanation is clear.

I agree that BAR is indeed a very efficient method, if it is sufficient for
including only the \lambda=0 and \lambda=1 states. However, BAR does not
seem to be reliable for my particular systems as shown in the variance of
the results. I could do multiple steps with BAR using your module, but as
you said, it is tedious to do so with the tools that are available at the
moment. The multi-step FEP does seem to work well for me with a small
variance and little hysteresis.