Re: [lammps-users] compute fep + fix adapt/fep + charge neutrality

Dear All,

I am trying to find out the free energy of CO2 (1 molecule) in a solvent in an NVE ensemble. My script looks las below.

variable lelj equal 1-${Nfac}*step #change lambda based on time

You could use ramp() to have a parameter changing, say, from 0 to 1 along the run.
Study the examples in lammps/examples/USER/fep (download the latest lammps version, where these were updated)

fix f1 all adapt/fep ${Nsol} pair lj/cut/coul/long/soft lambda 14 5 v_lelj pair lj/cut/coul/long/soft lambda 5 67 v_lelj pair lj/cut/coul/long/soft lambda 1*7 8 v_lelj atom charge 5 v_lelj atom charge 8 v_lelj

compute c1 all fep ${Temp} pair lj/cut/coul/long/soft lambda 14 5 v_lelj pair lj/cut/coul/long/soft lambda 5 67 v_lelj pair lj/cut/coul/long/soft lambda 1*7 8 v_lelj atom charge 5 v_lelj atom charge 8 v_lelj

fix f2 all ave/time 1 100 1000 c_c1[1] file energy.dat

I guess this is just a test run because statistics seem poor. Also, you should write at least the c_ID[1] and c_ID[2] elements of the compute because the c_ID[2] contains the averaged Boltzmann factor. Please read the documentation page for compute fep.

fix 2 all ave/time 1 1 1000 v_lelj file lambda2.dat

this will just print the variable. It is easier just to add it to thermo.

I am using ‘fix adapt/fep’ as I want to change the value of ‘lambda’ for pair style ‘lj/cut/coul/long/soft’ and atom charge every ‘Nsol’ time steps (here 5000).

I am facing some issues:

  1. When I add ‘atom charge’ for carbon (type 5) and oxygen (type 8) of CO2 to ‘fix adapt/fep’ and ‘compute fep’ I get a warning (not otherwise) ‘WARNING: System is not charge neutral, net charge =3’. I am not able to resolve this.

C and O atoms do not have the same charge (qO = - qC / 2), so you should ensure that charges vary during the simulation keeping the overall charge neutral (otherwise you’ll have errors from the long-range Coulomb. Again, look at the examples.

  1. I store value of lambda every 1000 timesteps in the file lambda2.dat. I see that the value is gradually decreasing. But for ‘compute fep’ I want to change the value at every ‘Nsol’ timesteps (here 5000). How can I be sure that ‘fix adapt/fep’ is providing a new value to ‘compute fep’ only every 5000 timesteps?

That’s what fix adapt/fep does. If you want to be sure, you can test with a very simple system (e.g. 2 LJ atoms at fixed positions in a box) and just print the energy terms with thermo. You can learn a lot by playing with this.

  1. If I want to change lambda, not as a function of time, how could I do that? I was thinking of ‘variable index’, but in that case, I would want ‘next’ to be a function of time, is that possible somehow?

Up to you, the possibilities are endless. As long as you understand what fix adapt and compute fep do.

  1. If I understand correctly, ‘compute fep’ uses FEP, but is it possible to switch to MBAR for eg.?

Compute fep calculates the delta Epot and the Boltzmann factor. You can use this to calculate free energies with FEP, finite-difference thermodynamic integration, or one-step BAR (which can work quite well for simple problems). Again, study the examples. MBAR is more sophisticated.

The LAMMPS version I’m using is 22Aug2018.

Should be fine for compte fep and fix adapt/fep. However, the examples are outdated.

Thank you very much.

You’re welcome.
Agilio Padua

Hi,

Dear Agilio Padua,

Thank you for your response. It was a good suggestion to go through the examples. I had two questions though. I would be glad if you could help me out.

  1. If I perturb the ‘lambda’ value for ‘pair lj/cut/coul/long/soft’, the ‘lambda’ value for Coulombic interactions is also perturbed right? In that case, is it really necessary to perturb the charge of atoms separately (I suppose the care is being taken by ‘lambda’ and may cause redundancy!)?

Yes, lambda also acts on the real-space part of Coulomb in the combined pair lj/cut/coul/long/soft.
And yes, you need to adapt and perturb the charges too in order to affect the k-space part of Coulomb.

  1. I went through the fep10 and fdti10 script for the hydration of methane, as my case is similar (energy of solvation of CO2 in a solvent). I’m just adding the example script (fep10) here for reference as I’d some doubts regarding the calculation of free energy.
    #########################################################################
    variable lambda equal ramp(1.0,0.0)
    variable q1 equal -0.24v_lambda
    variable q2 equal 0.06
    v_lambda

fix ADAPT all adapt/fep 100000 &
pair lj/cut/tip4p/long/soft lambda 12 34 v_lambda &
atom charge 1 v_q1 &
atom charge 2 v_q2 &
after yes

fix PRINT all print 100000 “adapt lambda = {lambda} q1 = {q1} q2 = ${q2}”

variable dlambda equal -0.05
variable dq1 equal -0.24v_dlambda
variable dq2 equal 0.06
v_dlambda

compute FEP all fep ${TK} &
pair lj/cut/tip4p/long/soft lambda 12 34 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2

fix PRINT1 all print 100000 “adapt dlambda1 = {dlambda} dq1 = {dq1} dq2 = ${dq2}”

fix FEP all ave/time 20 4000 100000 c_FEP[1] c_FEP[2] file fep10.lmp

#########################################################################

a) I see that the variables for ‘fix ADAPT’ (lambda, q1, q2) and ‘compute FEP’ (dlambda, dq1, dq2) are different. How in this case is ‘fix ADAPT’ influencing ‘compute FEP’? What would be the purpose of ‘dlambda’, ‘dq1’, and ‘dq2’ (all are constants although I see in fep10.lmp file c_FEP[1] and c_FEP[2] has different values for different timesteps) ?

The fix adapt modifies the parameters every 100 000 steps. That’s all it does. This sets the stages (equilibrium states that are ensemble-averaged) on which the perturbations are to be computed.

Then compute fep (through fix ave/time) performs the perturbations. In this example we have 20 stages: lambda goes from 1 to 0 over a total run of 2 000 000 steps, so the adapt changes parameters every 100 000 steps, for a delta of -0.05. In this manner, the delta G obtained at each stage corresponds to the transformation from one stage to the next (see the formula for the multi-stage FEP in the doc page of compute fep).
So, in this instance of the FEP method the duration of each stage in fix adapt is related to the delta used in the compute fep and to the duration of the run.

Actually, you don’t have to use successive stages as in these examples. Sometimes it’s better to use a separate run for just for the adapt, from which you save configurations for each stage. (Maybe these need to be equilibrated more at each stage, or maybe the run with the adapt was long enough.) Then take those as starting configurations and run independent trajectories at each stage to compute the perturbations, in which you can do +delta and -delta in the same run (2 compute feps). The choice depends on the size of the problem and on the machines available.

b) The code for ‘fdti10’ and ‘fep10’ is almost the same except for the value of ‘dlambda’ (which I suppose is a major determining factor I don’t know how, although maybe answer to a) helps). The output files have the ‘Ediff’ and ‘Boltzman factor’ values. So as I understand I could directly use the formula on the doc page (compute fep) and find out the free energy using FEP or FDTI right? For FDTI I could not understand the way of determination of the weight factor though ‘w’?

The idea of FDTI (see paper by Mezei, J Chem Phys 1987) is to use FEP on a very small delta around each stage in order to approximate the derivative (of the energy with respect to lambda) by a numerical derivative (finite difference). I picked +/-0.002 for the delta lambda because it is small enough to evaluate the derivative and not too small to cause too much numerical error (precision). But this may depend on the system.

The weights are those of some numerical integration (quadrature) scheme: you can use the trapezium rule (weights 1/2 1 1 1 … 1/2 ) or Simpson’s rule for example. More sophisticated quadratures may require non-evenly spaced points.

PS: In my previous mail I wanted to ask about BAR but mistakenly I wrote MBAR. Sorry about that. Also, right now it’s only a test run for me. I will run longer in other cases.

BAR can be surprisingly efficient. But in my group we use mostly FEP for bigger solutes. It is important to check if perturbations going up 0->1 or down 1->0 give similar profiles and close values for the free energy difference.

Agilio

Dear Agilio Padua,

Thank you for your response. It was a good suggestion to go through the examples. I had two questions though. I would be glad if you could help me out.

  1. If I perturb the ‘lambda’ value for ‘pair lj/cut/coul/long/soft’, the ‘lambda’ value for Coulombic interactions is also perturbed right? In that case, is it really necessary to perturb the charge of atoms separately (I suppose the care is being taken by ‘lambda’ and may cause redundancy!)?

  2. I went through the fep10 and fdti10 script for the hydration of methane, as my case is similar (energy of solvation of CO2 in a solvent). I’m just adding the example script (fep10) here for reference as I’d some doubts regarding the calculation of free energy.