Kspace ewald (coul/long potential) generates strange forces

LAMMPS LAMMPS General Discussion
1

Dear LAMMPS users,

I found that long-range Coulomb potential (lj/cut/coul/long) with kspace_style ewald generates forces between atoms, even if the atoms do not have interactions.
Is this an intended behavior of LAMMPS?

- Dimer calculation with charge but with zero interaction. - Note that there are three types written in the input script but the actual system only includes type 1 and 2. - The Coulomb potential is only between type 1-3 and 2-3. Thus, there should be zero forces. - lj/cut with eps=0 is overlay, and this is a dummy potential that gives zero energy and zero force.

The tested code is very simple with only two atoms, as follows,

atom_style charge
atom_modify map array sort 0 0
units           metal
boundary        p p p

read_data lmp_manual_in_mm1_2.data

pair_style hybrid/overlay lj/cut 12.0 lj/cut/coul/long 12.0 12.0
pair_coeff * * lj/cut 0 0
pair_coeff 2 3 lj/cut/coul/long 0 0
pair_coeff 1 3 lj/cut/coul/long 0 0
kspace_style ewald 1.0e-4

mass 1 22.989769274086011563440479221754
mass 2 35.449999990880691314032446825877
mass 3 12.010999996910237896940998325590

set atom 1 charge 0.8
set atom 2 charge -0.8
set atom 3 charge 0


thermo_style    custom step etotal temp lx vol press pxx pyy pxy
thermo          1

dump 4a all custom 1 dump.ewald_mm1 id type q x y z fx fy fz

#fix             f1 all nve

run 0

data file is shown below. (The name is “lmp_manual_in_mm1_2.data”)

2 	 atoms 
3  atom types
0.0                      20  xlo xhi
0.0                      20  ylo yhi
0.0                      20  zlo zhi


Atoms 

     1   1 0                      10                      10                      10
     2   2 0                      10                      10                      14

The dump file of the result, dump.ewald_mm1 shows

ITEM: TIMESTEP
0
ITEM: NUMBER OF ATOMS
2
ITEM: BOX BOUNDS pp pp pp
0.0000000000000000e+00 2.0000000000000000e+01
0.0000000000000000e+00 2.0000000000000000e+01
0.0000000000000000e+00 2.0000000000000000e+01
ITEM: ATOMS id type q x y z fx fy fz
1 1 0.8 10 10 10 1.89175e-18 1.57951e-18 0.0795596
2 2 -0.8 10 10 14 -1.89297e-18 -1.8637e-18 -0.0795596

This gives some forces in z direction that are not negligible.
If I calculate with “lj/cut/coul/cut”, short-range interaction without Ewald sum, the forces are just zero.
Would you please kindly tell me if there is a reason why force appears even if there is no atoms that have non-zero interaction in the system?

I checked this behavior in lammps-stable_29Sep2021 and 2022-09-18 version of windows executable (LAMMPS-64bit-latest.exe).

I would be happy to hear your opinion.

Best regards,

K. HISAMA

2

There is nothing strange about this. What you see is exactly the difference between using a cutoff and a long-range Coulomb solver which has effectively no cutoff.

3

Thank you for your response.
However, there should be no interaction between atom type 1 and 2.
There should be zero forces, definitely.

pair_coeff 2 3 lj/cut/coul/long 0 0
pair_coeff 1 3 lj/cut/coul/long 0 0

Sorry for my bad explanation.

4

You are only turning off the Lennard-Jones interactions. You still have Coulomb interactions and atoms with charges and thus there should be forces due to the charges from the Kspace style. That one always computes interactions between all atoms in the system. It is not a pair-wise computation, so how you assign pair styles in the hybrid pair style is irrelevant.

This can be easily confirmed by changing your thermo_style command to:

thermo_style    custom step etotal epair evdwl ecoul elong

Which produces the following output:

Per MPI rank memory allocation (min/avg/max) = 9.339 | 9.339 | 9.339 Mbytes
   Step         TotEng         E_pair         E_vdwl         E_coul         E_long    
         0  -1.5575404     -1.5575404      0              0             -1.5575404    
Loop time of 2.21e-06 on 1 procs for 0 steps with 2 atoms

As you can see, the contributions are from the long-range solver.

5

Thank you for your quick answer.
Now I see that it is not a pair-wise computation.