Hi Patrick,
Hi Agilio,
I have been experimenting with your USER-FEP package, and I now have a few more questions. I hope you can answer them if it is not too much trouble!
First, I am interested in using 2 lambda values - one for the LJ terms and one for the Coulombic terms. The goal of this is to scale out the charge interactions at about twice the rate of the LJ terms. However, I notice that in your methane hydration example, you only use one lambda value, and never explicitly scale the coulombic interaction. Why is this? Does scaling the 'lambda variable' in the soft potential do this automatically?
It is possible to use several lambda values, evolving in different ways. You will need several variables to steer those using fix adapt/fep.
You can decouple the van der Waals and Coulomb parts by using the hybrid pair styles, e.g. lj/cut/soft and coul/long/soft
And yes, one consequence of using the soft-core potential is that the “activation” of the sites is carried out through the lambda parameter and not through epsilon or sigma or the charge. Check the formulas of the soft core pair functions in the doc page.
Sometimes you will want to change a charge (or an epsilon or sigma) when a certain site is being transformed from one type to another but not vanishing. In these cases there’s no need for soft core potentials and you can use fix adapt/fep on the charge (or the epsilon or the sigma) directly.
Second, I notice that you do scale the charges in the other examples, but you scale them directly rather than using the 'scale' variable as in the standard fix adapt command. If I did this with the methane hydtration example, I think it would interfere with normal intra-solvent dynamics - perhaps this is why you didn't do it?
I think the standard fix adapt will scale all the charges in the system (this is useful to compute the electrostatic energy in an ionic material, which is what that fix was written for).
In solvation calculations, we may want charges of the solute to change, but not those of the solvent. So, fix adapt/fep allows you to only modify charges on selected atom types.
Third, for the fdti method, which is the one I am interested in using, you still use the compute fep command, but with a very small dlambda value - approximately 0.002. I read about this in you manual entry, but I wonder if there is any rhyme or reason to the value you chose? Could it be smaller?
In FDTI one uses a FEP calculation to evaluate numerically the small \delta G corresponding to a small \delta\lambda variation, and from there estimate the derivative of the free energy at that point. I just chose 0.002 as a compromise between linearity and resolution. I haven’t tested values much different from these in the examples. If you use larger than this, might as well do a FEP. If smaller, you may have precision problems.
I’d go for the traditional FEP method instead of FDTI. In the past I used FDTI, but at the time I didn’t have a nice tool to scan the lambda values along a simulation. So, it would have been cumbersome to run 50 or more trajectories at different lambdas to do FEP properly. FDTI may be better when you can only do a small number of evaluations at different lambda points. But it may be less precise than FEP.
So, because with one run using fix adapt/fep you can evaluate as many lambda points as you want, I found FEP to be better than FDTI, both in hysteresis (comparing 0->1 and 1->0 routes) and in accuracy.
You may want to try BAR method, which is surprisingly good! Maybe for your application you cannot get away with 1-step BAR, and in these cases you can do a multistep BAR (you’ll have to run several simulations independently, at different values of lambda, but you may need only one, maybe very few intermediate values).
I recommend FEP over FDTI. BAR is excellent if 1-step works, a bit more laborious if MBAR required.
Agilio