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

1 Like

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.