# [lammps-users] Reax and compute pe/atom

HI All,
I have a question on the compute pe/atom. I am using ReaxFF.

When I have a single O2 molecule, and I use the following commands:

group o2 type 2
compute o2pepa o2 pe/atom
compute o2pe o2 reduce sum c_o2pepa

I obtain that o2pe is equal to pe, i.e., to the total pe of the
system. I assume that is the energy of the O2 molecule using the
ReaxFF.

However, when I place another atom (specifically a Pt atom, for which
I have the parameters) very very far away from the O2 molecule, and I
re-run, I obtain that o2pe does not equal the energy computed before.
It seems that the presence of another atom, even though very far away,
influences the energy of the O2 molecule. Is that due to the Coulombic
energy term even at such large distances?

Moreover, when I turn the Pt atom to an oxygen atom, I do get that
o2pe becomes equal to the initially computed O2 pe.

Best regards,
Paolo Valentini

I would use compute pe/atom and dump these values to a dump
file for a few snapshots. The sum of the per-atom eng should
equal the total PE of the system. If you stare at the dump
file you should be able to see what portions have been
assigned to what atoms. If the Pt atom is close enough
to interact in some way with the O2, then it will change the
O2 per-atom eng. I don't recall what the cutoffs are in ReaxFF,
but they are finite, even for Coulombic.

Steve

When I have 3 atoms (i.e., O2 and a Pt atom), the sum of their pe/atom
energies does equal
the sum of the total potential energy. I do that using the following
commands (type 2 atoms are oxygen atoms, in
the parameter file the third element is oxygen and the seventh is platinum):

pair_style reax 10.0 1.0e-6
pair_coeff * * ffield.reax 7 3

neighbor 2.5 bin
neigh_modify delay 0 every 5 check no

group o2 type 2
compute o2pepa all pe/atom
compute o2pe all reduce sum c_o2pepa

The output is:

Step Temp E_vdwl E_coul PotEng KinEng TotEng o2pe
0 0 -130.02729 -6.6400871 -136.66738
0 -136.66738 -136.66738

However, even though the Pt atom is very far away from the O2
molecule (I am pretty sure it is outside any cutoff radius), the sum
of the pe/atom of the
two oxygen atoms does not equal (or it is not very close) to the
energy computed when the Pt
atom is not included in the system (i.e., when I only have the oxygen
molecule, where E_vdwl is roughly -130).
The commands I am using are the following:

pair_style reax 10.0 1.0e-6
pair_coeff * * ffield.reax 7 3

neighbor 2.5 bin
neigh_modify delay 0 every 5 check no

group o2 type 2
compute o2pepa o2 pe/atom
compute o2pe o2 reduce sum c_o2pepa

And the output is:

Step Temp E_vdwl E_coul PotEng KinEng TotEng o2pe
0 0 -130.02729 -6.6400871 -136.66738
0 -136.66738 -155.7718

where I would expect E_vdwl to be very close (if not equal) to o2pe.

Am I using the compute command correctly? Or maybe it has to do with
the parameters file?

Regards,
Paolo Valentini

I don't know about the interaction lengths, but this
command

compute o2pepa o2 pe/atom

will give you 0.0 for the energy of any atom not in the
group o2pepa - so I would use group all if you are trying
to account for all energy in the system. Again I would
print these energies individually in a dump file and see
what changes when you add/subtract the Pt atom.

Steve