Intermolecular potential energy

Continuing the discussion from Intermolecular potential energy:

I have been working on this exact problem and I think neither the compute group/group nor the equivalent rerun plus neigh_modify exclude molecule/intra are the correct way of computing the intermolecular energy. My understanding is that bonded pairs do contribute to the cohesive energy via the E_long term and the correct way to treat the Coulomb interactions with 3D boundary conditions plus a long-range solver is to correct the real-space pairwise Coulomb contribution (i.e. in any coul/long-type pair style) with the special_bond command. This is evident for small molecules like water: there, E_intra = E_vdwl + E_coul + E_long. For bigger molecules, excluding the intramolecular contribution from the neighbour list will only work for pair styles without long-range electrostatics terms.

The solution I came up so far is to rerun the system with a new “custom” topology, where all the atoms belonging to the same molecule are connected by dummy bonds. This hack was mentioned here: Modeling a Partially Rigid Molecule - #7 by Efrem_Braun1

I attach a simple toy system (an array of linear molecules with 4 atoms and opposite dipoles) which I have used to test the behavior of compute group/group and rerun. I believe that the correct intermolecular energy is the PotEng in toy03.log: 762.06174 kcal/mol.

toy01.log (5.0 KB)
toy01.lt (1.7 KB)
toy02.in (1004 Bytes)
toy02.log (3.9 KB)
toy03.in (970 Bytes)
toy03.log (3.8 KB)

I wonder if there is a better way to compute the intermolecular energy.
Any comment is greatly appreciated!

Cheers,
Otello

Hello there, I have an important update on this topic. Along with Stephen Farr, we have developed a new compute inter for the intermolecular energy of systems with long-range electrostatics. It also works for systems without long-range potentials, but the existing compute group\group can also be used for those. You can find the source code on Stefan’s fork on GitHub.

(the repo also contains a new coarse-grained force model called MOLC, but this will be the topic of another thread.)

I have prepared an input deck with a sample of 1000 n-decane molecules.
intermolecular.tgz (5.3 MB)
The force field is the GROMOS-ATB, composed of the usual bonds, angles, dihedral terms, Lennard-Jones spheres and point charges.

These are the relevant commands for every simulation:

# decane_01
compute inter1 all group/group all pair yes kspace yes molecule inter
compute inter2 all inter kspace yes # the new compute.

# decane_02
neigh_modify exclude molecule/intra all
compute inter1 all group/group all pair yes kspace yes # include inter- and intra-molecular contributions
compute inter2 all inter kspace yes # the new compute.
rerun decane_01.dump dump x y z vx vy vz box yes

# decane_03
special_bonds lj/coul 0.0 0.0 0.0
compute inter1 all group/group all pair yes kspace yes # include inter- and intra-molecular contributions
compute inter2 all inter kspace yes # the new compute.
rerun decane_01.dump dump x y z vx vy vz box yes

For simulation decane_02, I have used the method suggested by Steve Plimton here.
For simulation decane_03, I have used a dummy molecular topology where all intramolecular interactions are scaled to zero by adding extra dummy bonds between all pairs of carbon atoms (suggested here).

This table reports the values (kcal/mol) computed in the last step:

Filename E_\text{vdwl}+E_\text{coul}+E_\text{long} c_inter1 c_inter2 Note
decane_01 -8111.2936 -21622.214 -12451.946 production run
decane_02 -21622.21664 -21622.217 -21622.217 mailing list reference 1
decane_03 -12451.9485 -12451.948 -12451.948 mailing list reference 2

In simulation 01, the intermolecular energy from compute group/group is almost twice the one computed by the new compute and corresponds exactly to the sum of E_vdwl+E_coul+E_long in simulation 02, e.g. by excluding intramolecular pairwise energy through neigh_modify. Conversely, by using a custom topology where the special_bonds turn off all intra-molecular interactions, the sum of E_vdwl+E_coul+E_long in simulation 03 nicely corresponds to compute inter in simulation 01.

In fact, compute 1 all inter is nothing but a modified compute 1 all group\group all where all the intramolecular interactions have been scaled to zero. The important bit is that the weighting function for the real part of the Coulomb summation is computed with factor_coul = 0 applied to every atom pairs belonging to the same molecule. This correction is necessary to account for the long-range solver.

Finally, by assuming that the intramolecular energy of a molecule in the liquid and gas phases is the same, the enthalpy of vaporisation becomes:
\Delta H_\text{vap}(T) = -E_\text{inter}(T) +k_B T
from which the following enthalpies of vaporization at 300 K (kJ/mol) are computed:

exp (NIST) c_inter1 c_inter2
51.2(1) 93.1(1) 54.7(1)

I encourage the community to test this compute on cohesion energy, Hansen solubility parameters, enthalpies, and likes.
I will publish a paper giving more technical details, but I hardly doubt that any referee will dig into the source code. For me, this is the preferred way to present this work to this community for scrutiny, hoping that you will find it useful. Once validated, this work could be merged into the official LAMMPS distribution.

Best wishes,
Otello

2 Likes

Otello,

I installed LAMMPS from the source your provided but “compute inter” was not supported in the compiled package. Can you please check and let me know if I missed any specific setting/package while compiling the LAMMPS?

Best,
Amir

Hi, this is likely happening because the molc package was not compiled. Try this:

make yes-asphere yes-extra-dump yes-kspace yes-misc yes-molecule yes-rigid yes-molc

@Lukun
This is an unrelated inquiry. Please see the forum guidelines and post this as a new topic with a suitable subject line and more detailed information about your system.

okay.

Hello, I’m glad I found this topic.

I am interested in computing Hildebrand solubility parameters.

I think there is an other approach where we can use rerun and neigh_modify exclude molecule/inter instead of neigh_modify exclude molecule/intra. Then we substract E_{inter+intra} (from 1st run) and E_{intra} (from rerun) to obtain E_{inter}. What do you think ? I got a mean value of -12454.986 kcal/mol.
I updated the .tgz file you provided with a decane_04.in.

The new compute inter sounds interesting as it computes E_inter on the fly. I had to modified it a little because I am using the last version of LAMMPS and compiling with cmake. It looks like it is working as I obtained the same value as you.
It is the first time I look at a c++ code so I am not an expert and could have made mistakes. I also put the slightly modified MOLC package in the archive molc.tgz.
Does this seem correct to you?

intermolecular.tgz (5.3 MB)
molc.tgz (54.7 KB)

Thank you for porting the MOLC package to the latest dev, it seems correct to me.

The method you suggest also works fine, and does not require any new compute. To summarise, the inter-molecular (aka cohesion) energy for a molecular system with long-range electrostatics can be obtained from the equation:
E_\text{inter}=E_\text{pair}-E_\text{intra}
where E_\text{pair}=\sum E_\text{vdwl}+E_\text{coul}+E_\text{long} and E_\text{intra} is computed as follows:

compute intra all group/group all pair yes kspace yes molecule intra
variable enb equal evdwl+ecoul+elong
variable inter equal v_enb-c_intra
compute inter2 all inter kspace yes # the new compute.
compute inter3 all group/group all pair yes kspace yes molecule inter

The intermolecular energy computed in this way is numerically identical to that computed with the new compute inter command:

  v_inter       c_inter2       c_inter3          v_enb        c_intra
-12353.99     -12353.937     -21524.214     -8025.8405      4328.1496

To me this result means that the intermolecular energy computed with computed group/group is not correct, as it is equivalent to a rerun with neigh_modify exclude molecule/intra all. Another way to see that something does not add up is by summing c_intra+c_inter3, which does not return v_enb.

@stamoor, can we take this offline? Or, better, in the Open Review? I still need to publish this stuff, and your revision is definitely the best I can hope for.

@hothello I’m trying to get up to speed on the discussion. Are you saying that exclusions are not being properly handled with compute group/group? From your post, it seems like there are extraneous exclusions when there shouldn’t be any? For long-range electrostatics (i.e. kspace), there is no concept of exclusion. Every point charge interacts with every other charge (including itself!). So any exclusions have to be handled inside the short-range pair style, including those from kspace. In other words, extra interactions from kspace are subtracted out in the pair style.

It is possible the logic for compute group/group with exclusions is incorrect in the coul/long pair styles.