# Potential energy

Hi Lammps User.

I am trying to compute, during my simulation, energy differences betwen steps where the only difference is that I’ve added a particle to the simulation box, all the other particles remain in their original coordinates. In this situation the difference in energy between the system with N particles and the one with N+1 particle is only the interaction energy between the particle that I’m adding and the rest of the particles. As Lammps has a compute pe / atom command, I compared the two alternatives, computing energy differences of the systems and computing just the energy of the new particle. Unfortunately I’ve found differences.
In order to find what was going on I simulated a simple system of two water separated by their hydrogen bond distance. The model of water used was a coarse grain model based on a modify Stillinger-Weber potential
these are the results

PotEng (default thermo output) -5.4790852
compute all pe -5.4790852
compute all pe/atom (water 1) -2.7395426 !!!
compute all pe/atom (water 2) -2.7395426 !!!
sum of water 1 and water 2 -5.4790852

if I look the PotEng, I got a quite nice hydrogen bond energy around 5 Kcal/mol, but here is what I don’t understand, As I get from the manual the pe/atom is the energy of each particular atom interacting with the rest of the atoms in the simulation so as there are only two atoms here the energy of the water1-water2 should be the hydrogen bond energy -5.47 ,the same for the water2-water1 and that should be the total PotEng. So why does the pe/atom energy appears divided by two ?

I recall, I’m using a Stillinger-Weber with parameters for water molecules, these results where obtained by doing a minimization.

I don't understand exactly what you are doing,
but you appear to have a 2 atom system.
The PotEng printed out by thermo may or may
not be divided by the number of atoms in the system.
That depends on the setting of thermo_modify norm,
and the default is different for different units. See the
doc pages for details.

The sum of compute pe/atom over all atoms should
equal the total PotEng. It's not clear from your email
whether you are listing one atom's pe and how you
printed that out.

Since your values seem to be a difference of 2x for
a 2 atom system, I'm assuming some comboination
of the above explains it.

Steve