Strange behaviour of the EDIP potential

Hi, everyone!
I’ve noticed a strange behaviour of the current EDIP potential, calculating the pressure by two equivalent ways (calculating pressure by the Lammps command and calculating PV for each atom, summarizing and dividing on volume as it was shown in Lammps example for “compute stress/atom” ) and recieving different results. That’s strange because all the other potentials I’ve ever used (tersoff, sw, eam and meam) show exactly the same correct results,

LAMMPS (17 Oct 2012)
variable element index Si
variable el0 index 5.43
variable ep index 1.0075e6

#Initialization

units metal
dimension 3
boundary pp pp pp

atom_style atomic
neighbor 2 bin
neigh_modify delay 0

#Atom definition

lattice diamond 5.43
Lattice spacing in x,y,z = 5.43 5.43 5.43
region r_all block 0 10 0 10 0 10 units lattice
create_box 1 r_all
Created orthogonal box = (0 0 0) to (54.3 54.3 54.3)
1 by 1 by 2 MPI processor grid
create_atoms 1 box units lattice
Created 8000 atoms

mass 1 28.0855
pair_style edip
pair_coeff * * Si.edip Si

timestep 0.001
velocity all create 2500.0 87287

compute PV_ten_l all stress/atom
variable PV_l atom c_PV_ten_l[1]+c_PV_ten_l[2]+c_PV_ten_l[3]
compute p all reduce sum c_PV_ten_l[1] c_PV_ten_l[2] c_PV_ten_l[3]
variable press2 equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol)

thermo_style custom step temp pe etotal press v_press2 vol
thermo 100

fix 0nve all nve
fix 0T all temp/berendsen 2500 2500 1
run 3200
Memory usage per processor = 3.87498 Mbytes
Step Temp PotEng TotEng Press press2 Volume
0 2500 -37199.627 -34614.747 17521.179 17521.179 160103.01
100 1356.0998 -35880.988 -34478.846 4901.7935 4270.1586 160103.01
200 1094.2925 -35481.153 -34349.707 103.11081 -2042.9558 160103.01
300 1291.8962 -35563.75 -34227.991 2892.8914 1086.305 160103.01
400 1387.9292 -35546.915 -34111.863 3608.2855 1841.7892 160103.01
500 1475.1584 -35525.665 -34000.422 4387.668 2260.5719 160103.01
600 1535.0627 -35481.06 -33893.879 4467.0018 2260.2553 160103.01
700 1548.5237 -35392.636 -33791.537 3747.0632 1638.4163 160103.01
800 1618.2619 -35367.17 -33693.965 4389.1028 1183.5905 160103.01
900 1624.9449 -35279.089 -33598.974 3477.4914 -286.74412 160103.01
1000 1648.6682 -35212.565 -33507.921 3414.7402 -1142.8671 160103.01
1100 1659.2395 -35137.727 -33422.153 2763.4888 -1504.5809 160103.01
1200 1729.4129 -35127.341 -33339.211 3429.686 -1435.1799 160103.01
(on the fifth and sixth column the results must be the same results)

I’m looking forward to recieving your help.
Best regards,
Anton Rudenko

These 2 values should agree. It's curious they agree on the 1st timestep.
When they don't agree it is typically b/c the pair style is not
correctly tallying the per-atom virial when each component of the
force is computed. Since this is a USER-MISC pair style,
I'll refer you to the author. Luca - can you check your calls
to the Pair::tally() functions to see if they are all in place for
each force component?

Steve

I'll have a look.

luca

Some time ago a bug was detected in lines 469-472 about the sign of the over/under-coordination force forceModCoord_ij. If the a negative sign of this force is written on these three lines, then it should be also imposed on the call of the tally function a few lines below (line 484).

  Fernan

Yes, confirmed. The negative sign must be imposed also in the tally
call in line 484.

Here is the patch:

--- USER-MISC/pair_edip.cpp 2012-10-26 05:10:42.000000000 +0200
+++ USER-MISC/pair_edip.cpp 2012-11-22 12:06:35.000000000 +0100
@@ -482,7 +482,7 @@

         evdwl = 0.0;
         if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0,
- forceModCoord_ij, dr_ij[0], dr_ij[1], dr_ij[2]);
+ -forceModCoord_ij, dr_ij[0], dr_ij[1], dr_ij[2]);
     }
   }

After this patch, the test posted by Антон Руденко <[email protected]...>
gives correct results:

LAMMPS (22 Oct 2012)
  using 1 OpenMP thread(s) per MPI task
Lattice spacing in x,y,z = 5.43 5.43 5.43
Created orthogonal box = (0 0 0) to (54.3 54.3 54.3)
  1 by 1 by 1 MPI processor grid
Created 8000 atoms
Setting up run ...
Memory usage per processor = 5.07967 Mbytes
Step Temp PotEng TotEng Press press2 Volume
       0 2500 -37199.627 -34614.747 17521.179
17521.179 160103.01
     100 1356.0998 -35880.988 -34478.846 4901.7935
4901.7935 160103.01
     200 1094.2925 -35481.153 -34349.707 103.11081
103.11081 160103.01
     300 1291.8962 -35563.75 -34227.991 2892.8914
2892.8914 160103.01
     400 1387.9292 -35546.915 -34111.863 3608.2855
3608.2855 160103.01
     500 1475.1584 -35525.665 -34000.422 4387.668
4387.668 160103.01
     600 1535.0627 -35481.06 -33893.879 4467.0018
4467.0018 160103.01
     700 1548.5237 -35392.636 -33791.537 3747.0632
3747.0632 160103.01
     800 1618.2619 -35367.17 -33693.965 4389.1028
4389.1028 160103.01
     900 1624.9449 -35279.089 -33598.974 3477.4914
3477.4914 160103.01
    1000 1648.6682 -35212.565 -33507.921 3414.7402
3414.7402 160103.01

Loop time of 12.4296 on 1 procs (1 MPI x 1 OpenMP) for 1000 steps with
8000 atoms

luca

just posted a patch for this - thanks for the quick fix.

Steve