[lammps-users] cutoff in lj/gromacs/coul/gromacs


I can't figure out how is Coulomb energy computed in pair style

When I put two point charges +/-e 4A apart and use pair_style:
   pair_style lj/charmm/coul/charmm 9.0 10.0
I get E_coul = -33.2064 kcal/mol, which is correct.

However, when I use;
  pair_style lj/gromacs/coul/gromacs 9.0 10.0
I get E_coul = -19.2376 kcal/mol, which is incorrect.

Changing the cutoff radii to 90A and 100A, changes E_coul to -31.8095. I
was looking in pair_lj_gromacs_coul_gromacs.cpp, but I cannot find any
place where the inner cutoff radius would be used incorrectly. Bellow is
in.lammps and data.lammps that reproduce this behavior.

Did I miss any required settings or is there a bug?

Thank you,


units real
atom_style full
boundary f f f
dielectric 2.5
special_bonds lj/coul 0.0 0.0 1.0

#pair_style lj/gromacs/coul/gromacs 9.0 10.0
pair_style lj/charmm/coul/charmm 9.0 10.0

read_data data.lammps
thermo_style multi
run 0


NaCl 2 atoms
           2 atoms
           2 atom types
          -1.400493 5.759659 xlo xhi
          -1.759659 1.759659 ylo yhi
          -1.759659 1.759659 zlo zhi


    1 35.4530 # Cl
    2 22.9897 # Na

Pair Coeffs

    1 0.28330 3.51932
    2 0.50000 2.80099


       1 999 2 1.00000 0.00000 0.00000 0.00000
       2 999 1 -1.00000 4.00000 0.00000 0.00000

this seems odd - I'll take a look


The code in pair_lj_gromacs_coul_gromacs.cpp has
an extra term with coulsw5 in the Coulombic energy
expression that I don't understand. It is producing
the odd behavior in this case. If this is a bug,
then it would only be affecting the output energy, not
the forces, hence it wouldn't alter dynamics.

I need to ask the person who implemented this - he's
likely gone for the Christmas break ... so stay tuned.


I think Aidan just fixed this issue - see the 26Aug11 patch.
The problem was the energy shift term not being
added in correctly (across the entire range of the potential)
for lj/gromacs. It was OK for lj/gromacs/coul/gromacs, but that
confused me in my earlier reply.