Sequentially turning off coulombic and LJ interactions

LAMMPS LAMMPS General Discussion
1

Hello,

I am recently working on an alchemical transformation that involves turning off coulombic and LJ interactions sequentially between a cation and a framework. I have been setting up the input according to the fep/CH4hyd/fep10 example in LAMMPS. I wonder why it is necessary to scale down the charges while the coul lambda value is being scaled down as shown below in the first example scripts.

In a second example, fep/SCPEhyd, the decoupling was done in two steps: q and LJ sequentially. In this case, however, only the charges are being scaled down while the coul/long pair-wise interactions are not changed. Could anyone justify the differences between the two cases?

Thanks a lot !!

Best,
Jianming

fep10/CH4hyd

fix ADAPT all adapt/fep 100000 &
  pair lj/cut/tip4p/long/soft lambda 1*2 3*4 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.24*v_dlambda
variable dq2 equal  0.06*v_dlambda

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

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

fep/SPCEhyd example

fix ADAPT all adapt/fep 100000 &
  atom charge 1 v_qH &
  atom charge 2 v_qO &
  after yes

thermo_style custom step etotal ke pe evdwl ecoul elong temp press density v_lambda v_qO v_qH

variable dlambda equal -0.05
variable dqH equal  0.4238*v_dlambda
variable dqO equal -0.8476*v_dlambda

compute FEP all fep ${TK} &
  atom charge 1 v_dqH &
  atom charge 2 v_dqO &
  volume yes
2

Sorry, but what you are asking about is not about LAMMPS but about the method. That makes it off-topic for this forum category. Even though you are using LAMMPS, it seems to be working fine as it is. There are many review articles and text books about this kind of free energy calculations. I know for example a good review book edited by Chris Chipot (you should find it when you search for his name and “free energy”) that is worth reading.

3

I assume that since Coulomb repulsion is “softer” and (the energy) only goes as r^{-1}, while LJ repulsion is “harder” and goes as r^{-12}, the /soft styles (which use lambda) are necessary for LJ repulsion but not for Coulomb repulsion.

Hence, I guess, the first example script uses lambda since the LJ component is changing, while the second example script does not use lambda since the LJ component is not changing (in the second step when the LJ component is changed, lj/cut/soft and lambda are used).

Note this is only an explanation; I am not saying these scripts are correct (or wrong). Test on your own system and understand for yourself what is happening.