[lammps-users] Coulombic energy


I noticed something puzzling to me when playing with cutoff values in the calculation of the Coulombic energy. When I use pair style lj/cut/coul/long along with PPPM, I get consistent values for the Coulombic (E_coul+E_long) energy while changing the cutoff values. However, when I use lj/cut/coul/cut the Coulombic (E_coul) energy is all over the place. (The results are obtained from 1 step of simulation only.)

cutoff E (pppm) E (no pppm)
7 -76089 55231
10 -76143 -100749
12 -76149 -105903
15 -76112 -45346
17 -76122 -60270
20 -76113 -36366

I repeated the simulation for the same system using Amber. The results using particle mesh Ewald (PME) are almost identical to PPPM. The results without PME are about 10% lower but they do not show any appreciable change with cutoff value.

This could be a system specific as I have not seen this situation with homogenous systems. My system consists of a box of water (ice) sandwiched between a mixture of two types of small molecules. The input script I use is below.

I appreciate any input.

units real
neigh_modify delay 1 every 1 page 100000 one 20000
atom_style full
bond_style harmonic
special_bonds amber
angle_style harmonic
dihedral_style harmonic

pair_style lj/cut/coul/cut 12.0 22.0
#pair_style lj/cut/coul/long 12.0 22.0
#kspace_style pppm 1e-4

pair_modify mix arithmetic
special_bonds amber

boundary p p p
read_data data.inter

group coat type 1 2 3 4 5 6 7 8 9 12 13 14
group ice type 10 11

fix 2 coat shake 0.000001 20 0 m 1.0008
fix 3 coat nvt 375.0 375.0 100.0

thermo 1
thermo_style multi
dump 2 all atom 1 inter.dump
dump_modify 2 image yes scale yes
dump 3 all custom 1 inter_vdw.dump tag eng fx fy fz
restart 1 inter.restartB inter.restart
run 1


lj/cut/coul/cut is a very simple truncated potential, so it is not surprising that you’re seeing the energy jumping all over the place. I assume that the cutoff potential you’re using in AMBER is either a shifted or switched coulomb potential. You might try using lj/charmm/coul/charmm, which is a switched potential and probably more similar to what AMBER is doing. But for production runs, I’d definitely recommend sticking with a coul/long potential that properly deals with the long-range electrostatics. It is good to know that you’re getting good agreement with AMBER when doing long-range electrostatics.