LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Question about Ewald summation in 1, 4 pair interactions
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Question about Ewald summation in 1, 4 pair interactions

From: Axel Kohlmeyer <akohlmey@...24...>
Date: Fri, 6 Oct 2017 17:03:27 -0400

On Fri, Oct 6, 2017 at 6:34 AM, J. Rene-Espinosa <jr752@...690...> wrote:

Dear all

I am finding the following problem comparing results from GROMACS and
Lammps. When I am simulating TIP4P in both packages I am obtaining exactly

​please always report the version of LAMMPS you are using.​
the same potential energy and same density for several thermodynamic states
(i. e. 298K 1 bar, 298K 1000 bar 298K 2000 bar ...). Everything OK.

But when I am simulating octanol and I have intramolecular interactions,
like the 1,4 pair interactions, (I am using OPLS/AA), then I am not obtaining exactly the same densities nor the same potential energies.
If I compare term by term the contributions in both packages, I get the
same energies for bonds, angles and dihedrals. I also obtain the same LJ,
interaction (including 1,4), but not the Coulombic one, and that is provoking that I am not obtaining the same densities and energies.

I set special bonds as 0.0 0.0 0.5 as is defined in OPLS/aa both in GROMACS and in Lammps, for the LJ interaction is working nice but not for the Coulombic part although I set it for both Lj/coul.

If I remove the kspace style in Lammps and I set Coul/cut but using a long cut-off, in my case I have used 20 Amstrong, then I recovered a very similar densities both in Lammps and Gromacs for several thermodynamic states.

So this indicates that the Ewald summation is not working properly when is dealing with the 1,4 pair interactions. When I switch off Ewald and with a long cut-off (20 A) then is when I recovered very similar values between Gromacs and Lammps.

​the question ​here is now: is the ewald summation in LAMMPS or the one in Gromacs deviating from "correct" result?

​i suggest the gromacs developers on this subject. i know that some of them are involved in a project validating different MD codes​ for consistent results, and they should be interested and checking out your case with their tools. i am quite confident, that for any recent version of LAMMPS the handling of excluded coulomb interactions is correct.
several years ago, we had a review of that code, where the exclusions for coulomb were revised and thoroughly tested.


There must be something that probably I am missing.
Any help would be very useful for me.
Thank you very much.

I attached the input file where I have all the parameter defined for the case of octanol.
By the way, I am not using constraints neither in Gromacs nor in Lammps,

units           real
dimension       3
boundary        p p p
atom_style      full
#log            log.${simname}.txt

# Defining the 1-4 interaction
# OPLS weights the lj and coulomb interactions by 0.5

bond_style      harmonic
angle_style     harmonic
dihedral_style  opls
improper_style  none
pair_style      lj/cut/coul/long 14.000000 14.000000
pair_modify     mix geometric

special_bonds   lj/coul 0.0 0.0 0.5 angle no dihedral no
#pair_modify pair lj/cut/coul/long   special lj/coul 0.0 0.0 0.5

kspace_style    pppm 1.0e-5


#Pair Coeffs

pair_coeff           1           1         0.00000000         0.00000000
pair_coeff           2           2         0.16988631         3.12000000
pair_coeff           3           3         0.02997994         2.50000000
pair_coeff           4           4         0.06595586         3.50000000
pair_coeff           5           5         0.06595586         3.50000000

#Bonds coeffs

bond_coeff 1 5.52630171e+02       0.9450
bond_coeff 2 3.19785994e+02       1.4100
bond_coeff 3 3.39772619e+02       1.0900
bond_coeff 4 2.67820770e+02       1.5290

# Angle Coeffs

angle_coeff 1  5.49632177e+01        108.5
angle_coeff 2  3.49765931e+01        109.5
angle_coeff 3  4.99665616e+01        109.5
angle_coeff 4  3.29779306e+01        107.8
angle_coeff 5  3.74749212e+01        110.7
angle_coeff 6  5.83110012e+01        112.7

#Dihedral Coeffs

dihedral_coeff 1  2.12137771e-16    -0.00000000e+00     4.49699054e-01     -0.00000000e+00
dihedral_coeff 2  -3.55766934e-01    -1.73882201e-01     4.91670249e-01     -0.00000000e+00
dihedral_coeff 3  2.38845897e-06    -0.00000000e+00     4.67686539e-01     -0.00000000e+00
dihedral_coeff 4  1.70985598e+00    -4.99665616e-01     6.62557323e-01     -0.00000000e+00
dihedral_coeff 5  -0.00000000e+00    -0.00000000e+00     2.99799369e-01     -0.00000000e+00
dihedral_coeff 6  1.29913060e+00    -4.99665616e-02     1.99866246e-01     -0.00000000e+00

neighbor        2.000000  bin
neigh_modify    delay 0 every 1 check yes one 4000

velocity all create 298.0 128872 mom yes rot yes dist gaussian

fix             2 all npt temp 298.0 298.0 500 iso 1.0 1.0 500
restart         1000 old_config_1.dat old_config_2.dat

thermo 10
thermo_style custom step temp density etotal pe ke evdwl ecoul elong epair ebond eangle edihed
thermo_modify  flush yes

dump 1 all xyz 10000
timestep 0.3
run 5000000
write_data      data_restart_init

Check out the vibrant tech community on one of the world's most
engaging tech sites,!
lammps-users mailing list

Dr. Axel Kohlmeyer  akohlmey@...92......
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.