Total energy calculation

I’ve been playing with a very simple problem, a Lennard-Jones 7 cluster in 2D, and I noticed something I do not understand about the energy that is computed. In particular, I do not understand how it’s counting bonds.

I set up the cluster with the commands:

boundary p p p
units lj
atom_style atomic
atom_modify map array sort 0 0.0
dimension 2

region box block 0 20 0 20 0 10 units box
create_box 1 box

create_atoms 1 single 10 10 0
create_atoms 1 single 11 10 0
create_atoms 1 single 10.5 10.866025403784 0
create_atoms 1 single 9.5 10.866025403784 0
create_atoms 1 single 9 10 0
create_atoms 1 single 9.5 9.1339745962156 0
create_atoms 1 single 10.5 9.1339745962156 0

mass 1 1.0
pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0
pair_modify shift yes
neighbor 0.3 bin
neigh_modify delay 0 every 1 check yes

fix 99 all enforce2d
run 0

thermo 1
min_style cg
minimize 1.e-14 1.e-14 1000 1000

And if I look at the output, I see:

Step Temp E_pair E_mol TotEng Press
0 0 -0.099697525 0 -0.099697525 0.35246486
1 0 -1.7219978 0 -1.7219978 0.021817974
2 0 -1.7276334 0 -1.7276334 -0.015209308
3 0 -1.7415426 0 -1.7415426 0.0020181014
4 0 -1.741734 0 -1.741734 0.00045648444
5 0 -1.7417403 0 -1.7417403 -0.00028734981
6 0 -1.7417442 0 -1.7417442 -7.997489e-05
7 0 -1.7417445 0 -1.7417445 -6.8336763e-06
8 0 -1.7417445 0 -1.7417445 -5.1660289e-08
9 0 -1.7417445 0 -1.7417445 -2.9622362e-12
10 0 -1.7417445 0 -1.7417445 -1.1339129e-12

The thing that I want to highlight here is that the E_pair and TotEng are the same value, of about -1.74. Now, if I think about my cluster of 7 atoms, and assume that near equilibrium for this, all the pair distances are close to the equilibrium value of 2^(1/6). In this case, the energy of nearest neighbor pair bonds is about -1. Each of the atoms on the exterior of the heptameter has 3 nearest neighbors, and the atom in the center has 6, so the total energy would be about -24, which is quite a bit different than the TotEng value.

To put it another way, I would have assumed that

TotEng = sum_{i<j} phi(r_{ij})

But that seems not to be the case.

-gideon

With units lj the thermo_output is normalized, so TotEng is the average total energy. Since there is twelve bonds, that amounts to -1 * 12 / 7.0 which is about -1.714285. The additional energy is probably from the second-nearest neighbor interactions.

Ah, I just realized where you got 24 from: You were counting each bond twice. The center atom has 6 bonds and each outer has 3, but if you correct for the factor 2 (which goes automatically if you sum over only i<j) you get 12 too.

I’ve been playing with a very simple problem, a Lennard-Jones 7 cluster in
2D, and I noticed something I do not understand about the energy that is
computed. In particular, I do not understand how it’s counting bonds.

​please be careful with your nomenclature.​
there are no *bonds* being computed here. ​you have only *non-bonded*
interactions in your model, i.e. you only use a "pair style". LAMMPS will
build a neighborlist and then loop over all *pairs* within the cutof​f and
considering periodic boundary conditions.

so, to avoid confusion, i recommend you talk about pair-wise interactions
instead.

axel.

Thanks. The issue was cleared up for me by Stefan, having pointed out that LJ units normalizes extensive quantities.