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