[lammps-users] Erroneous non-bonded energies

LAMMPS LAMMPS Mailing List Mirror
1

Hi!

The error I describe here concerns the basic things, so it seems very strange.

First, a "formal" part. My system is GNU/Linux Fedora 6. LAMMPS of
February 08, 2008 was compiled with g++ with stub MPI, without FFT
library and also without any optional packages (just default set). No
compiler optimization keys were used.

Next, the problem.
Let us calculate with LAMMPS the non-bonded energies in two simplistic systems.
Each system comprises just two equal atoms in periodic box. The atoms
have charges +1 and -1.
The systems differ in that in the first one there are no chemical
bonds whereas in the second one the two atoms are connected with a
bond.

Here are the input files for the above systems:

1)"non-bonded system"

===== in.AA =====

#...

units real
atom_style full

special_bonds 0.0 0.0 0.0
kspace_style ewald 0.00001
pair_style lj/cut/coul/long 12.0

read_data data.AA

thermo_style custom pe evdwl ecoul elong
thermo 1000

run 0

===== end of in.AA ====

===== data.AA =====

#...

      2 atoms

      1 atom types

   0.0000 4.0000 xlo xhi
   0.0000 4.0000 ylo yhi
   0.0000 8.0000 zlo zhi

Masses

      1 1.0000

Pair Coeffs

      1 1.0000 4.0000

Atoms

      1 1 1 -1.0000 2.0000 2.0000 2.0000
      2 1 1 1.0000 2.0000 2.0000 6.0000

===== end of data.AA ====

2) "bonded system"

===== in.A-A =====

#...

units real
atom_style full

special_bonds 0.0 0.0 0.0
kspace_style ewald 0.00001
pair_style lj/cut/coul/long 12.0
bond_style harmonic

read_data data.A-A

thermo_style custom pe evdwl ecoul elong ebond
thermo 1000

run 0

===== end of in.A-A ====

===== data.A-A =====

#...

      2 atoms
      1 bonds

      1 atom types
      1 bond types

   0.0000 4.0000 xlo xhi
   0.0000 4.0000 ylo yhi
   0.0000 8.0000 zlo zhi

Masses

      1 1.0000

Pair Coeffs

      1 1.0000 4.0000

Bond Coeffs

      1 100.0000 4.0000

Atoms

      1 1 1 -1.0000 2.0000 2.0000 2.0000
      2 1 1 1.0000 2.0000 2.0000 6.0000

Bonds

      1 1 1 2

===== end of in.A-A ====

Now let us see the results (only energies here; full log.lammps files
see in attachment).

For the "non-bonded system" we have:

     PotEng E_vdwl E_coul E_long
-72.347992 -8.0582398 11.535252 -75.825004

These values are correct; I've reproduced them in MatLab with good accuracy.

Now let's predict the Coulomb and van der Waals energies for "bonded system".
They should differ from the above values by energies of corresponding
interactions between bonded atoms.
For van der Waals this energy is zero. So E_vdwl should remain -8.0582398.
For Coulomb this energy is 332.063709 * 1 * (-1) / 4.0 = -83.015927.
So E_coul should become 94.551179 kcal/mol. (E_long should not differ)

Now let us see the LAMMPS output:

     PotEng E_vdwl E_coul E_long E_bond
  1791.1817 -2.8680738 1869.8748 -75.825004 0

Now the question: WHAT'S WRONG ?

Am I missing something or it's a bug?

Best regards,

Vadim

PS: The same behavior takes place with ewald/n solver with pair_style
lj/coul long long 12.0 (when LAMMPS is built with optional package
USER-EWALDN)

log.lammps.AA (1.18 KB)

log.lammps.A-A (1.26 KB)

2

Read the doc page for the special_bonds command.
When 2 atoms are bonded the weighting coeff
is 0.0, you are turning off the vdwl and Coulomb
interaction between those 2 atoms, hence the
energy will change.

Steve

3

This is an interesting problem. The issue is that
your cutoff is bigger than the box length and so there
are many images of the single molecule interacting
with itself. LAMMPS is eliminating the pairwise interaction
for atom 1 with all images of atom 2. This probably
isn't the right thing to do, but it's not clear to me how
to figure out which ones to eliminate. For example, if
you have a bond length of 4 and a box size of 4, then
there can be many periodic images of the atom J that I is
bonded to that are close to I. Which one is it bonded to?

So I need to think about this one a bit. The safe thing
is not to have a cutoff distance greater than the box size
minus a bond length.

Steve