[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

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

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

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.
  2. 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?
  3. 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?
  4. If I understand correctly, ‘compute fep’ uses FEP, but is it possible to switch to MBAR for eg.?

The LAMMPS version I’m using is 22Aug2018.

Thank you very much.

Regards,
Kunal Mavani