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:
- 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.
- 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?
- 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?
- 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