Pressure Calculation , Reax/c

Hi All,
I have been trying to run a simulation with Reax/c Lammps and stand alone (original ) reax code.
When I compare energies for these two simulations at step 1 they are almost identical. However when I compare pressures for these two cases
they are different. i.e -1.5 GPa for Lammps and 3.9 GPa for reax code.
Is it possible to have a bug in Lammps/Reax/c pressure calculation part?
Regards.

Hi All,
I have been trying to run a simulation with Reax/c Lammps and stand alone
(original ) reax code.
When I compare energies for these two simulations at step 1 they are
almost identical. However when I compare pressures for these two cases
they are different. i.e -1.5 GPa for Lammps and 3.9 GPa for reax code.
Is it possible to have a bug in Lammps/Reax/c pressure calculation part?

since you have an alternate code, ​you are in a better position to track
this down than us who don't.
i suggest you set up a set of tiny test systems, and compare the results
for the virial of the two codes with what you get from F \dot r.

axel.

Hi Alex,
I tested couple of systems to diagnose the problem. Would you mind leading me where to look?

Test System -1

5

1.0

4.010000 0.000000 0.000000

0.000000 4.010000 0.000000

0.000000 0.000000 4.010000

Direct

1 0.333333 0.333333 0.333333 0 Ba 138

2 0.833333 0.833333 0.833333 0 Ti 48

3 0.333333 0.833333 0.833333 0 O 16

3 0.833333 0.333333 0.833333 0 O 16

3 0.833333 0.833333 0.333333 0 O 16

Box: 24.06 24.06 24.06 (6x6x6 : 320 atoms)

Pressure

x y z

-5.556 -5.556 -5.556 LAMMPS

2.5179 2.5179 2.5179 ADF

Test System-2

Same with System 5 , box : 44.06 44.06 44.06 20A vacuum in each direction

Pressure

x y z

0.6717 0.6717 0.6717 LAMMPS

1.800 1.800 1.8000 ADF

The only difference between Test System 1 and Test System 2 is 20 A vacuum in each direction. ADF-Reax pressures make sense in this regard. However LAMMPS calculates -5 Gpa for the periodic system and 0.6 GPa for the system with vacuum. The problem maybe with the periodic boundary conditions or ghost atoms?
Thanks in advance.

in.debug (4.85 KB)

Hi Alex,
I tested couple of systems to diagnose the problem. Would you mind leading
me where to look?

​a) i would first compare the forces. if the forces don't agree
sufficiently well, there is little purpose to try and debug ​the
virial/pressure.
b) the virial contributions are added to the global (or per atom)
accumulators with a variety of functions that have the word "tally" in
common. see here:

USER-REAXC]$ grep -n _tally *.cpp
reaxc_bond_orders.cpp:195: system->pair_ptr->v_tally(i,fi_tmp,delij);
reaxc_bond_orders.cpp:214: system->pair_ptr->v_tally(j,fj_tmp,delji);
reaxc_bond_orders.cpp:231: system->pair_ptr->v_tally(k,fk_tmp,delki);
reaxc_bond_orders.cpp:233: system->pair_ptr->v_tally(k,fk_tmp,delkj);
reaxc_bond_orders.cpp:251: system->pair_ptr->v_tally(k,fk_tmp,delki);
reaxc_bond_orders.cpp:253: system->pair_ptr->v_tally(k,fk_tmp,delkj);
reaxc_bonds.cpp:99:
system->pair_ptr->ev_tally(i,j,natoms,1,ebond,0.0,0.0,0.0,0.0,0.0);
reaxc_bonds.cpp:129:
system->pair_ptr->ev_tally(i,j,natoms,1,estriph,0.0,0.0,0.0,0.0,0.0);
reaxc_hydrogen_bonds.cpp:178:
system->pair_ptr->ev_tally3(i,j,k,e_hb,0.0,fi_tmp,fk_tmp,delij,delkj);
reaxc_multi_body.cpp:95:
system->pair_ptr->ev_tally(i,i,system->n,1,e_lp,0.0,0.0,0.0,0.0,0.0);
reaxc_multi_body.cpp:121:
system->pair_ptr->ev_tally(i,j,system->n,1,e_lph,0.0,0.0,0.0,0.0,0.0);
reaxc_multi_body.cpp:208:
system->pair_ptr->ev_tally(i,i,system->n,1,eng_tmp,0.0,0.0,0.0,0.0,0.0);
reaxc_nonbonded.cpp:180:
system->pair_ptr->ev_tally(i,j,natoms,1,pe_vdw,e_ele,
reaxc_nonbonded.cpp:293:
system->pair_ptr->ev_tally(i,j,natoms,1,e_vdW,e_ele,
reaxc_nonbonded.cpp:336:
system->pair_ptr->ev_tally(i,i,system->n,1,0.0,en_tmp,0.0,0.0,0.0,0.0);
reaxc_torsion_angles.cpp:467:
system->pair_ptr->ev_tally(j,k,natoms,1,eng_tmp,0.0,0.0,0.0,0.0,0.0);
reaxc_torsion_angles.cpp:469:
system->pair_ptr->v_tally4(i,j,k,l,fi_tmp,fj_tmp,fk_tmp,delil,deljl,delkl);
reaxc_valence_angles.cpp:393:
system->pair_ptr->ev_tally(j,j,system->N,1,eng_tmp,0.0,0.0,0.0,0.0,0.0);
reaxc_valence_angles.cpp:395:
system->pair_ptr->v_tally3(i,j,k,fi_tmp,fk_tmp,delij,delkj);

so these are the individual contributions to the stress tensor, which are
then converted to pressure.

that is about all the help i can give on this subject. beyond that, you
should try and discuss with the USER-REAXC author (cc'd).

axel.

Also pay special attention to the settings in the control file as one small different setting can affect the reulsts.

Ray

Have you checked if the charges are identical in both simulations? Which energy terms have different values?

Also, the example you’ve provided is incomplete - it lacks the force field you use.

Michal