Dear all,
I am attemping to selectively scale I J long range (kspace) electrostatic calculations in a water-methane simulation using the fix adapt command.
The general approach is to selectively scale the I and J interactions that correspond only to methanol -water non-bonded terms, while leaving intra-methanol and water-water terms unscaled. In particular, I would like to have enough control to do only electrostatic terms first, followed by the LJ terms.
My ultimate goal is to use LAMMPS to calculate the free energy of solvation for a methane molecule in water, following an approach similar to that described here (using NAMD):
http://www.ks.uiuc.edu/Training/Tutorials/namd/FEP/tutorial-FEP.pdf
Keeping long range electrostatics is important, so I have added “epsilon” and “scale” to the extract method in pair_lj_charmm_coul_long.cpp and .h (and also in ./src/KSPACE), recompiled, and run some tests.
I have found that e_vdwl and e_coul scale linearly with my linearly changing lambda_LJ and lambda_coul value, which is great. However, e_long does not change at all. I suspect that this may be due to the fact that it is not possible to decompose the kspace calculation into individual I and J interactions.
#Code snippet for lambda values and scaling#
variable lambda equal 1.0
variable lj equal -1*v_lambda+1
variable dlj equal -1
variable elec equal -2*v_lambda+1
variable delec equal -2
if “${elec} <= 0.0” then &
“variable elec equal 0.0” &
“variable delec equal 0.0”
#Methane atom types are 1 and 2, water types are 3 and 4.
#Therefore need to exclude 1-3, 1-4, 2-3, 2-4 interactions.
fix 2 all adapt 0 pair lj/charmm/coul/long epsilon 1 3 v_lj pair lj/charmm/coul/long scale 1 3 v_elec scale yes
fix 3 all adapt 0 pair lj/charmm/coul/long epsilon 1 4 v_lj pair lj/charmm/coul/long scale 1 4 v_elec scale yes
fix 4 all adapt 0 pair lj/charmm/coul/long epsilon 2 3 v_lj pair lj/charmm/coul/long scale 2 3 v_elec scale yes
fix 5 all adapt 0 pair lj/charmm/coul/long epsilon 2 4 v_lj pair lj/charmm/coul/long scale 2 4 v_elec scale yes
My first question: Is this correct, or have I missed a way to scale the e_long component selectively for I and J interactions?
As a work-around, I have thought of scaling the charges directly by using “pair_style hybrid/overlay” to create a clean separation between each group:
#Code snippet for hybrid/overlay#
pair_style hybrid/overlay lj/charmm/coul/charmm 8 10 0.0 1e-20 coul/long 10 coul/long 10
kspace_style pppm 1e-4
read_data solvate1.data
#Methane atom types are 1 and 2, water types are 3 and 4.
#Therefore need to create separation between 1-3, 1-4, 2-3, 2-4 interactions and all the rest.
pair_coeff 12 12 coul/long 1
pair_coeff 34 34 coul/long 1
pair_coeff 1 3 coul/long 2
pair_coeff 1 4 coul/long 2
pair_coeff 2 3 coul/long 2
pair_coeff 2 4 coul/long 2
However, this doesn’t seem to allow specification of different charges for coul/long 1 and coul/long 2 or to allow me to selectively scale coul/long 2.
To summarize, the selective scaling of the I J kspace calculation is the real sticking point in this. I would be greatful for any solution that allows me to do this selectively while keeping long range electrostatics in the system.
Kind regards,
Patrick