[lammps-users] Unable to understand calculation of ecoul in thermo_style


I am trying to extract pairwise Columbic interaction energy of water (tip3p) using ecoul in thermo_style. A single water molecule was initialized, energy minimized, equilibrated in nve ensemble and finally in nvt ensemble. The value of ecoul was around +39.

I need to understand whether all 3 combinations of interaction between 3 atoms of water (H-H, O-H and O-H) are included in it, or is it only the non-bonded interaction between H-H?

To understand it, I changed charges of oxygen and hydrogen to zero, in different combinations, and performed same steps explained before.
When I assigned the charge of O as -0.830 and both Hydrogen as 0.0, the ecoul was 0.0 as expected.
But when the charge of O as -0.830, one H as 0.415 and other H as 0 was assigned, ecoul value was around 25. I was expecting it to be negative, as charge of atoms are opposite in this combination.
Then I assigned charge of O as 0.0 and both H as 0.415, then ecoul was -11, where I was expecting a positive sign.

I have included the different combinations of charges in atom section of data file and corresponding truncated output in details.docx and attached with this mail.

Can you please help me to understand how ecoul is calculated in LAMMPS and why the sign of ecoul here is not as expected.

The input LAMMPS script and data file are also attached.

Thanks in advance.

water1.txt (533 Bytes)

in.water_single (786 Bytes)

details.docx (15.8 KB)

You didn’t use it, but the command special_bonds (https://docs.lammps.org/special_bonds.html) is important to understand. By default, all 1-2 (bond), 1-3 (angle), and 1-4 (dihedral) interactions ignore contributions from LJ and Coulombic potentials. This means that your system’s Coulombic potential can only be coming from the kspace style PPPM, which is unnecessary for a system of one molecule anyway. You can see this by removing the kspace_style command and switching to pair_style lj/cut/coul/cut, for example.

It would be better to make a physically accurate system (a fully solvated box) in order to understand Coulombic interactions.

Michael Jacobs
Dobrynin Group
Department of Chemistry
UNC at Chapel Hill

You are making your life needlessly complicated by using periodic boundaries and long-range coulomb interactions.

Because of that, “ecoul” is only part of the coulomb interaction, the real-space part. you have to add to this the value of elong to get the total coulomb interaction corresponding to an infinitely long coulomb cutoff.

If you look at the sum of ecoul and elong you should be getting a value close to (but not identical to) zero.

That is because the (full) coulomb between the 1-2, 1-3, and 1-4 neighbors are excluded, i.e. the direct interactions within the same molecule.
The pair style only computes a (damped!) short range interaction that must be augmented with (inversely damped) reciprocal space contributions.
Due to periodic images and the long-range nature of coulomb interactions the oxygen atom is still interacting with the periodic images of itself and the hydrogen atoms and one hydrogen atom with the periodic images of the other and so on.
that and the fact that you are not converging the long-range sum in kspace to the full will result in a small remaining total coulomb energy when computing the sum of ecoul and elong.

I suggest you have a closer look into a proper text book on classical MD force fields and how long-range coulomb is computed.

BTW: you are NOT simulating TIP3P, because that would be a rigid water model, but your input does not contain any such constraints.

1 Like


Then I assigned charge of O as 0.0 and both H as 0.415, then ecoul was -11, where I was expecting a positive sign.

because you are using a kspace style the actual total colomb energy is plus infinity. you have an infinite number of charged cells.
as soon as you system is not neutral, then your total coulomb energy with kspace diverges. but due to the way how the sum is done, one can discard the diverging term, which is equivalent to adding a compensating charge (uniformly distributed) into the boundary surrounding the (infinite) crystal (yes, it makes physically no sense, but is mathematically possible).
with that setup, you can (still) compute forces and have essentially a bogus energy, but since you care about the forces and not the (absolute) energy you can still do MD.