# E_pair does not equal E_coul+E_vdwl+E_long?

I just completed a simulation on a long chain alkane and from the thermo output I noticed that the outputted E_pair value does not equal the sum of E_coul, E_vdwl and E_long.

From the thermo_style documenation:

Epair = pairwise energy (evdwl + ecoul + elong)

Example output:

Step Temp TotEng E_vdwl E_coul E_pair E_long Enthalpy Density
0.0 296.86229 20395.576 -3674.6346 6355.8932 -2812.4255 -5493.6841 22274.883 0.78288505

My styles are:

units real
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style harmonic
pair_style lj/cut/coul/long 10.0 10.0
pair_modify mix geometric
special_bonds lj/coul 0.0 0.0 0.5
kspace_style pppm 0.0001

There was also a warning regarding the system having a net charge:

WARNING: System is not charge neutral, net charge = -7.74 (src/kspace.cpp:313)

I am wondering if this has something to do with it, does this mean that the simulation results are unreliable?

E_pair = -2812.4255
E_vdwl + E_coul + E_long = -3674.634 + 6355.8932 - 5493.6841 = -2812.4255

This looks correct to me.

No, it has (obviously) nothing to do with it, but it is a warning that should be taken seriously.
Any lattice sum will diverge for non-neutral cells and such an object would explode instantly in reality. The Ewald/PPPM formulation thus ignores divergent terms from the net charge which is physically equivalent to embedding the system into a distributed counter charge. While this is mathematically “clean” it can be problematic. The issue is to understand where the net charge originates from. Usually this will not happen unless you have either incorrect computed partial charges on atoms or have incorrectly rounded/truncated them. For an alkane chain you should have a net zero charge. In fact, each repeat unit should be neutral. So it is very likely that there is an error or an inaccuracy in your system preparation.

Thank you, my post processing script was causing the issue with the pairwise interaction sum.

Regarding the warning, would a system of, for example, phenol have the same problem? Would the same styles and settings be applicable?
(I built the system with moltemplate using the OPLSAA ff)

If you have neutral molecules then the system should be neutral. If moltemplate does not produce a neutral system (and I can only guess it is because of rounding/truncation of assigned partial charges) then this is a bug in the software or the parameters/inputs used.

MD follows the GI-GO principle and if your input is already tainted by rounding errors, then this doesn’t bode well for the reliability of any results, especially if those may be looking at subtle differences.

Hi, I scrutinized the input file and found the mistake (using the wrong type of atom from the oplsaa file). Thank you for your help.