Compute pe/atom versus compute group/group

LAMMPS LAMMPS General Discussion
1

Hello,

I am using LAMMPS version dated 23 Jun 2022 to study a phosphate crystal (Na4Fe3P4O15) that comprises of sodium, iron, phosphorous, and oxygen atoms. I used the following commands to calculate the total pair potential energy of sodium atoms in two ways:

compute s sodium pe/atom pair
compute spe sodium reduce sum c_s
compute sod_all sodium group/group all pair yes kspace no

On outputting c_sod_all and c_spe, I was expecting c_sod_all to be twice c_spe, since in the former (compute group/group), the energy for each pair is accumulated while in the latter (compute pe/atom), for each pair, the energy for that pair is split in half and added to each atom participating (see: computre group/group vs pe/atom).

However, the output I am getting is as follows:
Step c_spe c_sod_oxy c_sod_e_oxy c_sod_all
360000 -12343.569 -214911.13 179884.75 -35026.38
where c_sod_all is quite far from 2*c_spe. I do not understand what I am doing wrong. Any help is appreciated. Please find my complete code below:

#Initialization
read_restart <restart.file>

Potential Parameters

pair_style hybrid/overlay coul/long 15 morse 15 lj/cut 15
pair_coeff * * coul/long
kspace_style pppm 1e-5
pair_coeff <list of pair coeffs>

Simulation Specifics

group sodium type 1:4
group iron type 5:7
group phos type 8:11
group oxy type 12:26
group excp_oxy type 1:11 # Except oxygen

Computes

compute s sodium pe/atom pair
compute spe sodium reduce sum c_s
compute sod_oxy sodium group/group oxy pair yes kspace no
compute sod_e_oxy sodium group/group excp_oxy pair yes kspace no
compute sod_all sodium group/group all pair yes kspace no

thermo 400
thermo_style custom step c_spe c_sod_oxy c_sod_e_oxy c_sod_all
rerun dump.data skip 1 dump x y z
~

2

You have to think about both sodium-sodium (Na-Na) interactions and sodium-other (Na-X) interactions.

With a compute NaPE sodium pe*, all Na-X interactions contribute only half their potential energy (since Na is in the group but X isn’t), but all Na-Na interactions contribute the full potential energy (since both Na are in the group).

With a compute NaGG sodium group/group all, both Na-X and Na-Na interactions contribute all their potential energy.

So you should expect the result of the first to be half the result of the second, plus half of all Na-Na interactions. For your results, this means a total Na-Na interaction energy of about +10400 kcal/mol. You would know whether this is a reasonable value for your system, but it makes sense that this number is positive since it would be dominated by medium-range Coulomb repulsion.

*which should give the same result as summing over a compute pe/atom.

3

Thank you for the explanation! It makes a lot more sense now. Thus, since all Na-X (X=Na,Fe,P,O) pair-wise interactions contribute all their PE when compute NaGG sodium group/group all is used, the appropriate way to calculate the PE of all the sodium atoms due to all pair-wise interactions (excluding kspace) would be to halve the output obtained from this command.