Dear all,
I am trying to calculate the free energy of hydration of ethanol in TIP4P water using the
USER FEP method. I have a problem with the TIP4P and ethanol LJ pair coefficients.
If I run the simulation letting LAMMPS pick the i-j coefficients with the mix geometric
command then the simulation runs perfectly fine. If I choose them explicitly then
I needed to play around with some labels to get a decent run otherwise the carbon atom
and the oxygen atom of the water got way too close.
pair_coeff 1 5 lj/cut/tip4p/long/soft 0.1244 3.3151 1.0 #C_3 OW
Does anyone know why this happens and what is the correct way to choose the i-j 0.0 1.0
coefficients as above.?
This is part of the input file:
pair_style hybrid &
lj/cut/coul/long 10.0 10.0 &
lj/cut/tip4p/long/soft 5 6 5 3 0.125 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm/tip4p 1.0e-4
You should not use pair style hybrid, but the tip4p style for everything. Otherwise the coulomb contributions are not computed correctly.
Axel.
How come then the USER-FEP benchmark test for the hydration energy of methane in TIP4P water uses these
settings?:
created by fftool
print "Methane in TIP4P water"
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
special_bonds lj/coul 0.0 0.0 0.5
read_data data.lmp
<b>pair_style hybrid &
lj/cut/coul/long 10.0 10.0 &
lj/cut/tip4p/long/soft 3 4 2 2 0.125 1 0.5 10.0 10.0 10.0</b>
pair_modify tail no
kspace_style pppm/tip4p 1.0e-4
pair_coeff 1 1 lj/cut/coul/long 0.0660 3.5000 # C4H C4H
pair_coeff 1 2 lj/cut/coul/long 0.0445 2.9580 # C4H H
pair_coeff 2 2 lj/cut/coul/long 0.0300 2.5000 # H H
pair_coeff 3 3 lj/cut/tip4p/long/soft 0.1628 3.1644 1.0 # Ow Ow
pair_coeff 3 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # Ow Hw
pair_coeff 4 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # Hw Hw
pair_coeff 1 3 lj/cut/tip4p/long/soft 0.1036 3.3279 1.0 # C4H Ow
pair_coeff 1 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # C4H Hw
pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 1.0 # H Ow
pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # H Hw
I did not write those examples (keep in mind that large parts of LAMMPS are contributed), but the implicit use of the point M charge is only done inside the tip4p style. All other coulomb styles would consider those charges at the original location. If you need a soft/non-soft mix you would have to model tip4p with 4 explicit positions and fix rigid/small instead of fix shake.
Also, in hybrid pair styles you can have mixing of lj parameters only between the same substyle, even if two substyles have the same lj component in their potential.
Axel
Hi Adriano,
I think you are missing the point of using lj/cut/tip4p/long/soft here. During a FEP simulation, the interactions between the molecule and the solvent (in a hydration calculation for example) are progressively turned on or off (depending on the reaction path you take).
As explain in the manual, the lambda parameter you are talking about is a factor which softens the forces if set between 0 and 1. If at 0, the interaction is deactivated, if at 1 it is just like a LJ interaction. In the example, you will see that a lambda variable is defined using the ramp function and that this variable is used to scale the lambda parameter of the interaction and the charges of the methane using the adapt/fep fix. It might not be such a big deal if some atoms overlap while lambda is set to 0 initially, since they might not see each other. The other Lj parameter can be computed through mixing rule with usual caution. So using lj/cut/tip4p/long/soft between your ethanol atoms, lambda should be set to 1 and left untouched. Only cross interaction with waters (and charges) should progressively change.
Beside that, I’m unsure why using lj/cut/coul/long for ethanol-ethanol interaction is incorrect. But if Axel says so, there might be a good reason. The examples’ in files were written using automated tools (not mine) so that might explain the use of pair hybrid. By changing the example for the hydration of methane to using only lj/cut/tip4p/long/soft as mentioned above, I got results from fep01 case that were as close to the reference experimental value in the Readme file than the original example.
Regards,
Germain