 Software
 Open access
 Published:
An automated calculation pipeline for differential pair interaction energies with molecular force fields using the Tinker Molecular Modeling Package
Journal of Cheminformatics volume 16, Article number: 96 (2024)
Abstract
An automated pipeline for comprehensive calculation of intermolecular interaction energies based on molecular forcefields using the Tinker molecular modelling package is presented. Starting with nonoptimized chemically intuitive monomer structures, the pipeline allows the approximation of global minimum energy monomers and dimers, configuration sampling for various monomer–monomer distances, estimation of coordination numbers by molecular dynamics simulations, and the evaluation of differential pair interaction energies. The latter are used to derive Flory–Huggins parameters and isotropic particle–particle repulsions for Dissipative Particle Dynamics (DPD). The computational results for force fields MM3, MMFF94, OPLSAA and AMOEBA09 are analyzed with Density Functional Theory (DFT) calculations and DPD simulations for a mixture of the nonionic polyoxyethylene alkyl ether surfactant C_{10}E_{4} with water to demonstrate the usefulness of the approach.
Scientific Contribution
To our knowledge, there is currently no open computational pipeline for differential pair interaction energies at all. This work aims to contribute an (at least academically available, open) approach based on molecular force fields that provides a robust and efficient computational scheme for their automated calculation for small to mediumsized (organic) molecular dimers. The usefulness of the proposed new calculation scheme is demonstrated for the generation of mesoscopic particles with their mutual repulsive interactions.
Introduction
The quantitative description of nonbonding interactions between molecules is fundamental to understanding and designing chemical processes in materials and life sciences [1, 2]. In contrast to covalent bonding within molecules, nonbonding intermolecular interactions comprise dispersed variations of electromagnetic interactions like dipole/dipole, dipole/induced dipole, induced dipole/induced dipole (van der Waals) interactions, hydrogen bonding, (partial) charge interactions, π–π, cation/anion–π or polar πeffects. Each different spatial configuration of two molecules can be assigned a corresponding nonbonding intermolecular interaction energy. Its determination is challenging because the intermolecular interactions are generally small compared to covalent bonding and especially to the total energy of a molecular configuration.
For an accurate quantitative description, a quantum chemical treatment with a suitable model chemistry is commonly advised, e.g., application of Density Functional Theory (DFT) with an appropriate combination of functional and basis set: the complexation energy (i.e., the energy difference between a specific dimer configuration and the two monomer molecules that form it) can then quantify the intermolecular interaction. In particular, SymmetryAdapted Perturbation Theory (SAPT) allows direct computation of nonbonding intermolecular interactions (i.e., without the need to calculate the total energies of monomers and dimer) and provides a physically meaningful decomposition of its contributing (electrostatics, induction, dispersion, shortrange repulsion) terms [2]. Recent DFTSAPT approaches have demonstrated a comparatively fast calculation in combination with remarkable accuracy for small to mediumsized dimers up to the adenine–thymine base pair [3, 4].
When multiple spatial dimer configurations need to be sampled, quantum chemical approaches become increasingly expensive due to their considerable computational complexity, as a single calculation can easily take minutes or even longer. As a purely classical alternative, molecular force fields may be employed instead: they allow intermolecular energy calculations within fractions of a second for a specific spatial dimer configuration, i.e., an acceleration by several orders of magnitude. This comes at the expense of accuracy, as a given force field can easily lead to deceptively erroneous results.
This work aims at providing a robust automated molecular forcefield based calculation pipeline for comprehensive estimation of mutual intermolecular energies of a set of small to mediumsized monomer molecules. For a chosen force field, the monomer molecules must be provided with an initial, chemically intuitive spatial geometry with associated atom types (where for small molecules a planar 2D geometry provided by any 2D structure editor seems to be sufficient). Then, each monomer start geometry is globally optimized to its minimum energy conformer. Several monomer–monomer distances with multiple configurations at each distance are sampled to obtain a nearglobal minimum energy dimer. This dimer is then locally optimized towards its global minimum. The sampled configurations can be averaged by Boltzmann weights to get mean nonbonded dimer interactions at different distances. The calculations can be extended to obtain differential molecule pair interaction energies, which describe the excess intermolecular interaction of two different molecules in comparison to the average interaction of the molecules themselves (see Eqs. 1 and 2 in “Methods” section below), where molecular dynamics (MD) simulations are used to estimate the coordination numbers of neighboring molecules.
The resulting interaction energies may be useful for different purposes. Based on the comprehensive sampling, suitable spatial start configurations can be obtained for more elaborate (e.g., quantum chemical) refinement procedures. Pairwise nonbonded intermolecular interactions can be considered as moleculepair descriptors for Cheminformatics tasks like molecular similarity estimation. Differential molecule pair interaction energies play an important role in statistical thermodynamics, e.g., for the quantitative estimation of excess quantities that determine the properties of mixtures [5]. In particular, they can be related to Flory–Huggins interaction parameters for polymer models (Eq. 9), which in turn can be used to describe isotropic particle–particle repulsions for mesoscopic simulation approaches (Eq. 10) [6]. The latter application is particularly studied in this work.
The backbone of the calculation pipeline is implemented using the Java programming language. All forcefieldbased calculations are performed with the Tinker Molecular Modelling package [7]. The Tinker package is modularized into standalone, taskbased executables (marked in italics in “Methods” section below), which fit well into the Java backbonedriven parallelized computational scheme that fully exploits the computational capabilities of multicore workstations. All force fields provided by Tinker can be used for the calculation pipeline.
Methods
The forcefieldbased intermolecular energy \({E}_{ij}^{C}\left(r\right)\) between two molecules \(i\) and \(j\) is determined by different noncovalentbonding contributions (van der Waals and partial charge interactions, hydrogen bonding etc.) and depends on the intermolecular distance \(r\) between the centers of the molecules as well as their relative spatial configurations \(C\). For each specific distance \({r}_{fix}\) there is a minimum energy configuration \({C}_{min}\) with \({E}_{ij}^{{C}_{min}}\left({r}_{fix}\right)\le {E}_{ij}^{C}\left({r}_{fix}\right)\). The global minimum energy dimer is characterized by a distinct distance \({r}_{min}\) so that \({E}_{ij}^{min}={E}_{ij}^{{C}_{min}}\left({r}_{min}\right)\le {E}_{ij}^{C}\left(r\right)\), i.e. \({r}_{min}\) is the distance between the centers of two molecules \(i\) and \(j\) when the dimer geometry corresponds to the global energy minimum. If different spatial configurations with a fixed intermolecular distance \({r}_{fix}\) are averaged, an averaged intermolecular energy for this distance \(\langle {E}_{ij}\rangle \left({r}_{fix}\right)\) is obtained: The corresponding minimum energy distance \({r}_{\langle E\rangle ,min}\) with \({\langle {E}_{ij}\rangle }^{min}=\langle {E}_{ij}\rangle \left({r}_{\langle E\rangle ,min}\right)\le \langle {E}_{ij}\rangle \left(r\right)\) does not necessarily coincide with minimum distance \({r}_{min}\) of the global minimum energy dimer. Among the concrete averaged configurations with a fixed intermolecular distance \({r}_{fix}\) the configuration \({C}^{*}\) with the minimal intermolecular energy is denoted \({E}_{ij}^{{C}^{*}}\left({r}_{fix}\right)\ge {E}_{ij}^{{C}_{min}}\left({r}_{fix}\right)\).
A differential pair interaction energy describes the excess intermolecular interaction of molecules \(i\) and \(j\) in comparison to the average interaction of the molecules themselves. This may be defined with respect to the global minimum energy dimers
or corresponding averages (see Eqs. 5 and 6 below for calculation details)
A positive (negative) differential pair interaction energy indicates a less (more) favorable intermolecular interaction compared to the individual ones. When differential pair interaction energies are applied to lattice models, all vertices of the lattice have a fixed number of surrounding neighbours (e.g. \(Z=4\) in two dimensions or \(Z=6\) in three dimensions). In contrast, continuum models can (and usually do) have different coordination numbers \({Z}_{ij}\), where the number \({Z}_{ij}\) of molecules \(j\) surrounding molecule \(i\) can (and usually does) differ from the number \({Z}_{ji}\) of molecules \(i\) that surround molecule \(j\). Equations 1 and 2 can be extended accordingly to
and
respectively, where superscript Z denotes the coordination number extension in contrast to Eqs. 1 and 2, \({E}_{ij}^{min}={E}_{ji}^{min}\) and \({\langle {E}_{ij}\rangle }^{min}={\langle {E}_{ji}\rangle }^{min}\) but \({Z}_{ij}\) may be different from \({Z}_{ji}\), i.e. \({Z}_{ij}\ne {Z}_{ji}\).
Thus, the estimation of differential pair interaction energies requires two steps: (I) A calculation scheme to obtain the different (averaged) molecule pair interaction energies \({E}_{ij}^{min}\) (\({\langle {E}_{ij}\rangle }^{min}\)), and (II) a procedure to estimate the different coordination numbers \({Z}_{ij}\). In the following, the concrete implementation of a corresponding automated calculation pipeline for a selected molecular force field using the Tinker Molecular Modelling package is described. All dimerrelated calculations are started with the global minimum energy conformers of the two constituent monomer molecules (see following “Global minimum energy monomers” section), where multiple dimer configurations are analyzed with these global minimum energy monomers (see following “Global minimum energy dimers” section) to obtain an approximate configuration for the global minimum energy dimer. The latter is achieved by optimizing this final approximate dimer configuration without constraints using all atomic degrees of freedom, so that the monomers are no longer confined to their individual global minimum energy conformers.
Global minimum energy monomers
The global minimum energy conformers are derived in advance from conformer search procedures: STARTING from an initially defined chemically intuitive geometry of a monomer molecule, a first geometry improvement is achieved with Tinker optimize using the Optimally Conditioned Variable Metric (OCVM) optimization technique [by default with a root mean square (RMS) gradient of 0.01 kcal/mole/Å] [8]. The resulting locally optimized geometry is then transferred to a lowmode conformational search (LMOD) with Tinker scan to find the minimum energy conformer (where by default, the RMS gradient for the LMOD procedure is set to 0.0001 kcal/mole/Å, the energy threshold for local minima is set to 100 kcal/mole, torsion angles are automatically selected, and five eigenvectors are used for the local search) [9]. For small molecules with O(10) number of atoms, the detected minimum energy conformers usually coincide with the global minimum energy conformers.
Global minimum energy dimers
To approximate the global minimum energy dimer of a pair of molecules \(i\) and \(j\), the centers of the molecules are positioned at different distances ranging (by default) from 3 to 16 Å in steps of 0.5 Å where the initial relative configuration of the two molecules is arbitrary. For each distance, a configuration sampling procedure is performed, which is sketched in Fig. 1. A (unit) sphere around each center is constructed with a number of \({N}_{sphere}\) evenly spaced points being generated on each sphere using a Fibonacci lattice as an adequate approximation (in comparison to a latitude–longitude lattice the surface points generated by a Fibonacci lattice are more evenly spaced with a smaller axial anisotropy) [10]. By rotating the molecules around their centers so that two adjacent spherical surface points and the centers of both molecules lie on a straight axis, \({N}_{sphere}\times {N}_{sphere}\) configurations are generated for which the corresponding interaction energies are determined by Tinker analyze (with settings to only compute the nonbonding interactions for a significantly accelerated performance). In addition, for each single configuration the second molecule is rotated for a number \({N}_{rot}\) of angles around the axis with corresponding interaction energy calculations, so that in total \({N}_{sphere}\times {N}_{sphere}\times {N}_{rot}\) spatial configurations are sampled for a single fixed distance between the monomer molecules (with each monomer being constrained to its individual global minimum energy conformer). The distance resolution is then refined around the distance with the detected minimum interaction energy dimer by reducing the step size from 0.5 to 0.1 Å, and then again from 0.1 to 0.01 Å (compare Fig. 2), so that a final near minimum energy dimer configuration is evaluated for the latter resolution. This resulting configuration \({C}^{*}\) is then optimized with Tinker optimize without constraints using all atomic degrees of freedom (i.e. the monomers are no longer confined to their individual global minimum energy conformers) to arrive at the forcefield dependent approximation for \({E}_{ij}^{min}\) of the two molecules \(i\) and \(j\), the global minimum energy dimer using the specified forcefield.
Configuration sampling
If configuration sampling is desired for a specific distance \({r}_{fix}\) of the dimer molecules, the (already) obtained interaction energies of the \({N}_{sphere}\times {N}_{sphere}\times {N}_{rot}\) configurations for this distance may be averaged with Boltzmann weights \({w}_{C}\left({r}_{fix}\right)\) for a defined temperature (with \({\text{k}}_{\text{B}}\), the Boltzmann constant, and \(\text{T}\), thermodynamic temperature) using \({E}_{ij}^{min}\):
Then \(\langle {E}_{ij}\rangle\) can be evaluated from the different \(\langle {E}_{ij}\rangle \left(r\right)\). Figure 2 shows the quantitative results for the acetic acid dimer using the Merck molecular force field (MMFF94) [11].
Coordination numbers
The coordination numbers \({Z}_{ij}\) are estimated by MD simulations. The simulation box construction is based on a pure water box at 298 K to get consistent results. A water molecule has a van der Waals volume of \({V}_{vdW}^{{H}_{2}O}=17.35\) Å^{3}, in a pure water box it occupies a volume of \({V}_{box}^{{H}_{2}O}=30.00\) Å^{3} at 298 K due to its density [12] and molar mass [13]. This relation is mapped to other molecules \(X\) with
so that the edge length \(a\) of a cubic simulation box of \(N\) molecules \(X\) is given by
The van der Waals volumes are approximated with the VABCVolume [14] descriptor of the Chemistry Development Kit (CDK) [15, 16]. A simulation box with a defined number \(N\) of molecules \(j\) (default is 400) and defined edge length \(a\) is created using Tinker xyzedit. Then a single molecule \(i\) is added to the box, where Tinker xyzedit automatically removes molecules \(j\) to keep the defined density of the simulation box. The (possibly unfavorable) start configuration is optimized with Tinker minimize to avoid atomic contacts that could lead to instabilities. The following MD simulation uses Tinker dynamics with a step size of one femtosecond and an Andersen thermostat [17] for temperature equilibrium (default is 298 K). 10,000 (default) initial steps are used for box equilibration, followed by several hundred thousand steps with data recording (default is 400,000). For each recorded simulation step the number of molecules \(j\) surrounding the single molecule \(i\) is analyzed. This is done either by a “bruteforce” counting approach or, alternatively, by the cellindex method with periodic boundary conditions [18]. The “bruteforce” approach calculates the distances between each atom of the single molecule \(i\) and each atom of all molecules \(j\). Alternatively, the cellindex method only considers the (drastically reduced) distances between each atom of single molecule \(i\) and the atoms of solvent molecules \(j\) in neighbouring cells. The criterion for including a molecule \(j\) as a relevant neighbor for the coordination number \({Z}_{ij}\) is based on the distance between its atoms and those of molecule \(i\). If the distance of the respective atoms is less than or equal to the sum of their van der Waals radii plus an arbitrary “catch radius” (default is 1 Å), the molecules are considered neighbors. In Fig. 3, a snapshot of a simulation step is illustrated. For each simulation step, the number of neighboring molecules \(j\) is determined. The average over all recorded steps is used to estimate the coordination number \({Z}_{ij}\), see Fig. 3.
Flory–Huggins and mesoscopic repulsion parameters
Differential pair interaction energies may be directly utilized to estimate Flory–Huggins interaction parameters \({\chi }_{ij}\) by
to describe polymer solutions [5], with \(\Delta {\langle {E}_{ij}\rangle }^{Z}\) being defined in Eq. 4.
For “bridging the gap between atomistic and mesoscopic simulation” (Groot and Warren [6]), the interacting molecules can be identified with the particles of “bottomup” mesoscopic Dissipative Particle Dynamics (DPD), where the microscopic Flory–Huggins interaction parameters \({\chi }_{ij}\) can be directly related to mesoscopic isotropic particle–particle repulsions \({a}_{ij}\left(T\right)\) (expressed in units of \({\text{k}}_{\text{B}}\text{T}\), with \({\rho }_{DPD}\) being the dimensionless DPD density, refer to [6] for details)
which determine the conservative DPD forces \({\underline{F}}_{ij}^{C,DPD}\):
with \({\underline{F}}_{ij}^{C,DPD}\), \({\underline{F}}_{ij}^{C,Bond}\), soft repulsive DPD force and harmonic bond force on particle i exerted by particle j; \({a}_{ij}\), maximum isotropic repulsion between particles \(i\) and \(j\); \({\underline{r}}_{ij}={\underline{r}}_{i}{\underline{r}}_{j}={r}_{ij }{\underline{r}}_{ij}^{0}\); \({\underline{r}}_{ij}^{0}\), unit vector. The numerical factor (3.4965) in Eq. 10 is derived from Eq. 24 in reference [6] where the inverse value (0.286) is given.
In interplay with the dissipative (frictional) forces \({\underline{F}}_{ij}^{D}\) and random forces \({\underline{F}}_{ij}^{R}\) the conservative forces \({\underline{F}}_{ij}^{C}\) guide the trajectories \({\underline{r}}_{i}\left(t\right)\) of the DPD particles according to Newton’s equation of motion:
with \({m}_{i}\), \({\underline{r}}_{i}\), mass and position vector of particle \(i\); \(t\), time; \({\underline{F}}_{i}\), total force on particle \(i\) exerted by other particles \(j\); \(N\), total number of particles in simulation; \({\underline{F}}_{ij}^{C}\), \({\underline{F}}_{ij}^{D}\), \({\underline{F}}_{ij}^{R}\), conservative, dissipative, and random force on particle \(i\) exerted by particle \(j\).
Thus, the described calculation pipeline may be applied to construct a forcefieldbased particle set for DPD simulations.
Pipeline code availability
The pipeline code is written in Java and openly available at [19]. A dedicated installer executable for the Java pipeline code, which comprises a full Java runtime environment, is available at [20] for the Windows operating system. For Linux operating systems a zip file is available at [20]. According to its licensing terms the Tinker executables for optimize, scan, analyse etc. have to be downloaded from its website [21] into a specified pipeline directory: a detailed instruction how to perform this is provided at [22]. A set of standalone Mathematica notebooks [23] for result visualizations is provided at [24]. A test pipeline is available at [19] to ensure proper installation.
Pipeline calculation performance
Calculation of a full single differential pair interaction energy for the force fields MM3, MMFF94 and OPLSAA with default settings (\({N}_{sphere}\times {N}_{sphere}\times {N}_{rot}=144\times 144\times 16=\text{331,776}\) dimer configurations for each fixed distance to approximate the intermolecular interaction energies, 10,000 equilibration steps and 400,000 simulation steps for the MD simulations to estimate the coordination number with 400 molecules in the box) takes several hours, where the AMOEBA09 force field requires a multiple. Since the pipeline supports comprehensive calculation parallelization for a set of monomer molecules, a single differential pair interaction energy can be obtained on average in less than an hour on a 64core AMD Ryzen™ Threadripper™ PRO 5995 CPU workstation [25].
DFT calculations for result evaluation
DFT calculations are performed with Gaussian 16 [26] and analyzed with GaussView 6 [27]. All molecular geometries are optimized using the dispersioncorrected wB97XD functional [28] with the 6–311++G(d,p) basis set where counterpoise calculations are used to obtain basis set superposition error (BSSE) corrected interaction energies. All Gaussian jobs files used are openly available at [29].
DPD simulations for result evaluation
All DPD simulations of this study are performed with the MFsim simulation system [30] using the Jdpd simulation kernel [31]. All constructed particle sets and MFsim simulation jobs are openly available at [32].
Results and discussion
To demonstrate the applicability of the different steps of the calculation pipeline, several small molecules are selected: Water (abbreviated H2O), methane (Me), ethane (Et), methanol (MeOH), dimethyl ether (Me2O) and the acetic acid dimer (HAc). The calculation results for this molecule set are evaluated and compared with alternative approaches and experimental results. All averaged energies and MD simulations refer to a temperature of \(T=298 K\). All intermolecular energy calculations were performed with \({N}_{sphere}\times {N}_{sphere}\times {N}_{rot}=144\times 144\times 16=\text{331,776}\) dimer configurations for each fixed distance of the molecule centers.
Acetic acid dimer
The acetic acid dimer is stabilized by two hydrogen bonds and has a planar geometry. Calculation results with the MMFF94 force field are shown in Figs. 2 and 4. Figure 2 depicts the minimal energy configuration \({C}^{*}\) energies \({E}_{ij}^{{C}^{*}}\left({r}_{fix}\right)\) for each specific distance \({r}_{fix}\) in red, exhibiting a single minimum at r_{fix} = 4.91 Å with \(E_{ij}^{C*}\) (4.91 Å) = − 15.9 kcal/mole. The corresponding averaged intermolecular energies \(\langle {E}_{ij}\rangle \left({r}_{fix}\right)\) are shown in blue, where the single minimum distance coincides at r_{fix} = 4.91 Å with \(\left\langle {E_{ij} } \right\rangle\) (4.91 Å) = − 15.4 kcal/mole in this case. The two hydrogen bonds of the minimal energy configuration \({C}^{*}\) have a length of 1.69 Å and a distance of 2.63 Å between the corresponding donor and acceptor oxygen atoms (see Fig. 4), which is close to the experimental values of 2.68 Å [33].
With the final optimization of the near minimal energy configuration \({C}^{*}\) the global MMFF94 forcefieldbased minimum energy configuration \({C}_{min}\) with \({E}_{ij}^{min}={E}_{ij}^{{C}_{min}}\left({r}_{min}\right)=17.6\) kcal/mole is obtained, which is 1.7 kcal/mole below the sampled minimal energy \({C}^{*}\) configuration: the finally optimized dimer keeps the distance of r_{min} = 4.91 Å but shows a planar geometry with a slightly reduced hydrogen bond length of 1.63 Å and a distance of 2.62 Å between the donor and acceptor oxygen atoms (see Fig. 5).
Mutual dimer calculations
Table 1 presents the mutual dimer calculations: it comprises \(\langle {E}_{ij}\rangle\) and \({E}_{ij}^{{C}^{*}}\) for the detected minima (compare Fig. 2) as well as the global forcefieldbased energy minimum \({E}_{ij}^{min}\) for force fields MMFF94, AMOEBA09, MM3 [34] and OPLSAA [35, 36] with the water models TIP3P [37] and TIP5P [38]. As expected, the nonpolar pure alkyl (methane and ethane) dimers exhibit only small interaction energies, the acetic acid dimer with two hydrogen bonds shows the largest interaction, and the polar dimers are in between. There are clear differences between the force fields, with the OPLSAA (TIP5P) interactions for the alcohol–water dimers being the strongest. On average, MM3 differs most significantly from the other force fields.
For the HAc–HAc and the Me–Me dimer the results for force fields OPLSAA, AMOEBA09 and MMFF94 were compared with corresponding DFT single point calculations (denoted DFT sp). In addition, the spatial configurations with \({E}_{ij}^{min}\) were used as start geometries for DFT geometry optimizations (denoted DFT opt), see Table 2. The DFT calculations indicate that the automated pipeline leads to acceptable near minimum energies and corresponding spatial configurations—with individual exceptions: in contrast to the MMFF94 and OPLSAA force fields, AMOEBA09 results in an eclipsed minimum energy configuration for the Me–Me dimer instead of a staggered one (see Fig. 6, this eclipsed configuration is maintained by the DFT geometry optimization), but this finding has no significant influence on the subsequent investigations due to its low energetic effect (see Table 2). The difference in interaction energies between the DFT opt and the OPLSAA force field result for the HAc–HAc dimer is most significant.
Coordination numbers \({Z}_{ij}\)
Table 3 contains the mutual coordination numbers \({Z}_{ij}\) (at \(T=298 K\)), where a single molecule i is surrounded by 400 molecules j, using the MMFF94, OPLSAA (with TIP3P and TIP5P water models), and AMOEBA09 force fields. For comparison, static packing [40] results are included which were taken from previous research [41] and obtained with the commercial Blends software of Materials Studio [42] using the Condensedphase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) force field [43].
Figures 7 and 8 show the similarity of trends in coordination number assignment. While the MDderived coordination numbers are similar for the different force fields used, the results for the COMPASS static packing approach are significantly reduced by about 50%. The spread of the values derived from MD (Fig. 7) is significantly higher than that of the static packing approach. Interestingly, linear mapping of the coordination numbers to the interval [0,1] yields an approximate overlap of the results (Fig. 8).
Repulsion parameters \({a}_{ij}\)
For the different force fields, particle sets for mesoscopic DPD simulations with the isotropic mutual repulsions \({a}_{ij}\) for methane (Me), ethane (Et), methanol (MeOH), ethanol (EtOH), dimethyl ether (Me2O), and water (H2O) were constructed. The offdiagonal \({a}_{ij}\) values of a particle set were linearly scaled with MFsim [30] so that the maximum absolute deviation between the smallest \({a}_{ij}\) value and the diagonal values \({a}_{ii}=24.83\) (for a thermodynamic temperature of 298 K) is 20. For the OPLSAA force field the water models TIP3P and TIP5P were considered, denoted OPLSAA (TIP3P) and OPLSAA (TIP5P). An additional OPLSAA (TIP5P) particle set, denoted OPLSAA (TIP5P \({Z}_{ij}=1\)), was created which is solely based on the minimal averaged intermolecular energy \(\langle {E}_{ij}\rangle\) with a fixed coordination number \({Z}_{ij}=1\) for all dimers. With force field MM3, interaction energies can be calculated only, so a combined particle set was created using MM3 for interaction energy calculation and the MMFF94 force field for MDbased coordination number estimation, denoted MM3/MMFF94. The particle set from [41, 44], based on the COMPASS force field, is used for comparison.
Figure 9a–f display the different repulsion parameters \({a}_{ij}\). The red dashed line indicates the diagonal value of 24.83. A crucial difference between the different particle sets is the water–methanol (H2O–MeOH) repulsion: for the COMPASS, OPLSAA (TIP5P) and OPLSAA (TIP5P, \({Z}_{ij}=1\)) force fields, this repulsion is the smallest one (and thus the base value for offdiagonal repulsion parameter scaling), whereas for force fields OPLSAA (TIP3P), AMOEBA09, MMFF94, and MM3/MMFF94 the water–dimethyl ether (H2O–Me2O) repulsion is minimal. Interestingly, the different water models TIP3P and TIP5P of the OPLSAA force field exhibit this crucial difference, emphasizing their relevance. Note, that the differences of the nonwater repulsions for OPLSAA (TIP3P) and OPLSAA (TIP5P) are caused by the different scaling due to these different base values for offdiagonal repulsion parameter scaling. The COMPASS force field exhibits an extraordinary difference for the water–methane (H2O–Me) and water–ethane (H2O–Et) repulsions, which is otherwise not visible. A significant difference between the OPLSAA (TIP5P) and OPLSAA (TIP5P, \({Z}_{ij}=1\)) force field is the water–dimethyl ether (H2O–Me2O) repulsion. These obvious differences between the force fields lead to different levels of usefulness for DPD simulation approaches, where even a single difference in dimer interaction can become a decisive factor.
DPD simulations
The created particle sets were used for DPD simulations of mixtures of water with the nonionic polyoxyethylene alkyl ether surfactant C_{10}E_{4}, where “C_{10}” indicates the number of 10 carbon atoms in the alkyl chain of the lyophobic part, and “E_{4}” represents the number of 4 lyophilic ethylene oxide units [44]. A stable lamellar L_{α} phase is formed by a C_{10}E_{4}/water mixture around 298 K with a C_{10}E_{4} mass fraction of 0.75. The performance of the particle sets for the different force fields may be evaluated by monitoring the emergent formation of C_{10}E_{4} bilayers from initial random mixing in the simulation box [44]. The DPD simulations are carried out with the settings outlined in [44], using the SPICES 9Me–4Me2O–MeOH fragmentation for C_{10}E_{4} [45, 46], an integration step size of 0.04, and a deactivated periodic boundary in zdirection (to induce bilayer formation in the xyplane).
For particle set OPLSAA (TIP5P \({Z}_{ij}=1\)), a stacked bilayer superstructure emerges at simulation step 62,000, see Fig. 10 and Table 4, which was even below the COMPASS particle set from [44] with 116,000 steps, where the emerged bilayer structure corresponds well to the one reported in [44], see Fig. 11. The OPLSAA (TIP5P) particle set required more than twice as many simulation steps (132,000) and the OPLSAA (TIP3P) particle set more than tenfold as many (846,000 steps). The AMOEBA09 and MMFF94 particle sets show a bilayer superstructure formation, but the bilayers do not align parallel to the xyplane (as induced by the periodic boundary conditions) but parallel to the yzplane. Using the hybrid MM3/MMFF94 particle set, no bilayer was formed within 1,000,000 simulation steps. For the specified DPD simulation task, the OPLSAA (TIP5P \({Z}_{ij}=1\)) particle set can be regarded as the most suitable choice, which is in good agreement with experimental findings.
Conclusions
The outlined automated comprehensive calculation of intermolecular interaction energies based on molecular forcefields shows satisfactory results for small molecule interactions and can even be successfully used to estimate mesoscopic simulation parameters. Special care should always be taken, as individual force fields can lead to erroneous results. Therefore, despite all automation, manual checking of the results is still essential.
The new calculation pipeline can be easily extended to additional force fields (which may require their conversion into the Tinker format). Therefore, calculating differential pair interaction energies will advance with the progress in improving the underlying molecular force fields. Due to the modularized pipeline approach, alternative modeling packages can also be used for certain computational tasks if they can provide the required specific functions.
The robustness of the outlined computational pipeline can be seen as a crucial advance, as the estimation of differential pair interaction energies along alternative computational paths often led to ambiguous results, which in particular prevented the construction of consistent particle–particle repulsions for larger mesoscopic particle sets [47]. The construction of mesoscopic sets with dozens of particles (including hundreds or thousands of mutual particle repulsions) is now within practical reach.
Availability and requirements

Project name: Mesoscopic Interaction Parameter Estimation with Tinker for Java (MIPET4Java)

Project home page: https://github.com/zielesny/MIPET

Current version: v1.0.0.0

Operating system(s): Windows (× 64), Linux (× 64)

Programming language: Java

Other requirements: Java v17.0.4 or higher, Tinker Molecular Modeling Package v8.10.2

Licence: GPL3.0

Any restrictions to use by nonacademics: While the backbone code of this project is not restricted to academic use, the Tinker Molecular Modeling Package is subject to corresponding restrictions, see [21] for license details.
Availability of data and materials
The source code of the calculation pipeline and the description for downloading the Tinker Molecular Modeling Package are available on GitHub at https://github.com/zielesny/MIPET.
Abbreviations
 AMOEBA09:

Atomic Multipole Optimized Energetics for Biomolecular Simulation force field 09
 BSSE:

Basis set superposition error
 COMPASS:

Condensedphase Optimized Molecular Potentials for Atomistic Simulation Studies
 CDK:

Chemistry Development Kit
 DFT:

Density Functional Theory
 DPD:

Dissipative Particle Dynamics
 LMOD:

Lowmode conformational search
 MD:

Molecular dynamics
 MM3:

Molecular mechanics force field 3
 MMFF94:

Merck molecular force field 94
 OCVM:

Optimally Conditioned Variable Metric
 OPLSAA:

Allatom optimized potentials for liquid simulations force field
 RMS:

Root mean square
 SAPT:

SymmetryAdapted Perturbation Theory
References
Engel T, Gasteiger J (2018) Chemoinformatics: basic concepts and methods. WileyVCH, Weinheim
Stone AJ (2013) The theory of intermolecular forces, 2nd edn. Oxford University Press, Oxford
Heßelmann A, Korona T (2014) Intermolecular symmetryadapted perturbation theory study of large organic complexes. J Chem Phys 141:094107. https://doi.org/10.1063/1.4893990
Heßelmann A (2018) DFTSAPT intermolecular interaction energies employing exactexchange KohnSham response methods. J Chem Theory Comput 14:1943–1959. https://doi.org/10.1021/acs.jctc.7b01233
Hill TL (1986) An introduction to statistical thermodynamics. Dover Publications, New York
Groot RD, Warren PB (1997) Dissipative particle dynamics: bridging the gap between atomistic and mesoscopic simulation. J Chem Phys 107:4423–4435. https://doi.org/10.1063/1.474784
Rackers JA, Wang Z, Lu C et al (2018) Tinker 8: software tools for molecular design. J Chem Theory Comput 14:5273–5289. https://doi.org/10.1021/acs.jctc.8b00529
Davidon WC (1975) Optimally conditioned optimization algorithms without line searches. Math Program 9:1–30. https://doi.org/10.1007/BF01681328
Kolossváry I, Guida WC (1999) Lowmode conformational search elucidated: application to C39H80 and flexible docking of 9deazaguanine inhibitors into PNP. J Comput Chem 20:1671–1684. https://doi.org/10.1002/(SICI)1096987X(19991130)20:15%3c1671::AIDJCC7%3e3.0.CO;2Y
González Á (2010) Measurement of areas on a sphere using fibonacci and latitudelongitude lattices. Math Geosci 42:49–64. https://doi.org/10.1007/s110040099257x
Halgren TA (1996) Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94. J Comput Chem 17:490–519. https://doi.org/10.1002/(SICI)1096987X(199604)17:5/6%3c490::AIDJCC1%3e3.0.CO;2P
Riddick JA, Bunger WB, Sakano T, Weissberger A (1986) Organic solvents: physical properties and methods of purification, 4th edn. Wiley, New York
PubChem Water. https://pubchem.ncbi.nlm.nih.gov/compound/962. Accessed 5 Feb 2024
Zhao YH, Abraham MH, Zissimos AM (2003) Fast calculation of van der Waals volume as a sum of atomic and bond contributions and its application to drug compounds. J Org Chem 68:7368–7373. https://doi.org/10.1021/jo034808o
Chemistry Development Kit. https://cdk.github.io/. Accessed 5 Feb 2024
Willighagen EL, Mayfield JW, Alvarsson J et al (2017) The Chemistry Development Kit (CDK) v2.0: atom typing, depiction, molecular formulas, and substructure searching. J Cheminform 9:33. https://doi.org/10.1186/s1332101702204
Andersen HC (1980) Molecular dynamics simulations at constant pressure and/or temperature. J Chem Phys 72:2384–2393. https://doi.org/10.1063/1.439486
Allen MP, Tildesley DJ (2017) Computer simulation of liquids, 2nd edn. Oxford University Press, Oxford
MIPET. https://github.com/zielesny/MIPET. Accessed 5 Feb 2024
MIPETv1.0.0.0. https://github.com/zielesny/MIPET/releases/tag/MIPET. Accessed 5 Feb 2024
Tinker Molecular Modeling Package. https://dasher.wustl.edu/tinker/. Accessed 5 Feb 2024
MIPET/README.md at main zielesny/MIPET. https://github.com/zielesny/MIPET/blob/main/README.md. Accessed 5 Feb 2024
Wolfram Mathematica: Modern Technical Computing. https://www.wolfram.com/mathematica/. Accessed 5 Feb 2024
MIPET/Visualizaton Mathematica notebooks at main zielesny/MIPET. https://github.com/zielesny/MIPET/tree/main/Visualizaton%20Mathematica%20notebooks. Accessed 5 Feb 2024
(2022) AMD Ryzen^{TM} Threadripper^{TM} PRO 5995WX Prozessor. https://www.amd.com/de/products/cpu/amdryzenthreadripperpro5995wx. Accessed 5 Feb 2024
Frisch MJ, Trucks GW, Schlegel HB, et al (2016) Gaussian 16 Rev. C.01. https://gaussian.com/gaussian16/. Accessed 5 Feb 2024
Dennington R, Keith TA, Millam JM (2019) GaussView Version 6. https://gaussian.com/gaussview6/. Accessed 5 Feb 2024
Chai JD, HeadGordon M (2008) Longrange corrected hybrid density functionals with damped atom–atom dispersion corrections. Phys Chem Chem Phys 10:6615–6620. https://doi.org/10.1039/B810189B
MIPET/Gaussian job files at main zielesny/MIPET. https://github.com/zielesny/MIPET/tree/main/Gaussian%20job%20files. Accessed 5 Feb 2024
van den Broek K, Daniel M, Epple M et al (2020) MFsim—an open Java allinone richclient simulation environment for mesoscopic simulation. J Cheminformatics 12:29. https://doi.org/10.1186/s13321020004329
van den Broek K, Kuhn H, Zielesny A (2018) Jdpd: an open java simulation kernel for molecular fragment dissipative particle dynamics. J Cheminformatics 10:25. https://doi.org/10.1186/s1332101802787
MIPET/C10E4water bilayer formation study at main zielesny/MIPET. https://github.com/zielesny/MIPET/tree/main/C10E4water%20bilayer%20formation%20study. Accessed 5 Feb 2024
Derissen JL (1971) A reinvestigation of the molecular structure of acetic acid monomer and dimer by gas electron diffraction. J Mol Struct 7:67–80. https://doi.org/10.1016/00222860(71)900081
Allinger NL, Yuh YH, Lii JH (1989) Molecular mechanics. The MM3 force field for hydrocarbons. 1. J Am Chem Soc 111:8551–8566. https://doi.org/10.1021/ja00205a001
Jorgensen WL, TiradoRives J (1988) The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin. J Am Chem Soc 110:1657–1666. https://doi.org/10.1021/ja00214a001
Jorgensen WL, Maxwell DS, TiradoRives J (1996) Development and testing of the OPLS allatom force field on conformational energetics and properties of organic liquids. J Am Chem Soc 118:11225–11236. https://doi.org/10.1021/ja9621760
Jorgensen WL, Chandrasekhar J, Madura JD et al (1983) Comparison of simple potential functions for simulating liquid water. J Chem Phys 79:926–935. https://doi.org/10.1063/1.445869
Mahoney MW, Jorgensen WL (2000) A fivesite model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions. J Chem Phys 112:8910–8922. https://doi.org/10.1063/1.481505
Matthews GP, Smith EB (1976) An intermolecular pair potential energy function for methane. Mol Phys 32:1719–1729. https://doi.org/10.1080/00268977600103031
Fan CF, Olafson BD, Blanco M, Hsu SL (1992) Application of molecular simulation to derive phase diagrams of binary mixtures. Macromolecules 25:3667–3676. https://doi.org/10.1021/ma00040a010
Truszkowski A, Epple M, Fiethen A et al (2013) Molecular fragment dynamics study on the water–air interface behavior of nonionic polyoxyethylene alkyl ether surfactants. J Colloid Interface Sci 410:140–145. https://doi.org/10.1016/j.jcis.2013.07.069
BIOVIA Materials Studio  BIOVIA  Dassault systèmes®. https://www.3ds.com/productsservices/biovia/products/molecularmodelingsimulation/bioviamaterialsstudio/. Accessed 5 Feb 2024
Sun H (1998) COMPASS: an ab initio forcefield optimized for condensedphase applications overview with details on alkane and benzene compounds. J Phys Chem B 102:7338–7364. https://doi.org/10.1021/jp980939v
Bänsch F, Steinbeck C, Zielesny A (2023) Notes on molecular fragmentation and parameter settings for a dissipative particle dynamics study of a C10E4/water mixture with lamellar bilayer formation. J Cheminformatics 15:23. https://doi.org/10.1186/s1332102300697w
van den Broek K, Daniel M, Epple M et al (2018) SPICES: a particlebased molecular structure line notation and support library for mesoscopic simulation. J Cheminformatics 10:35. https://doi.org/10.1186/s1332101802947
Zielesny A (2021) SPICES—A particlebased Molecular Structure Line Notation and Support Library for Mesoscopic Simulation. https://github.com/zielesny/SPICES. Accessed 5 Feb 2024
Truszkowski A (2015) Simulation von Peptiden, Proteinen und Biomembranen mit Molekularer Fragmentdynamik (MFD). Dissertation, University of DuisburgEssen, Germany
Acknowledgements
The authors would like to thank Tim Clark for valuable initial discussions as well as his “dooropener” support, and especially Andreas Heßelmann for important methodological advice. Jay Ponder provided extensive help and advice on the Tinker package, for which we are very grateful. Karina van den Broek and Matthias Epple kindly supported the initial concept, Veit Hucke thankfully contributed the rotation matrix code. JanMathis Hein and Martin Urban generously helped with the molecular coding and the construction of start geometries. Jonas Schaub thankfully revised the GitHub repository and tested the installation procedure. Support from the Carl Zeiss Foundation and the Open Access Publication Fund of the Westphalian University of Applied Sciences is gratefully acknowledged. Lastly, the authors would like to thank the communities that created the open software libraries used in this work.
Funding
Open Access funding enabled and organized by Projekt DEAL. This work was supported by the CarlZeissFoundation.
Author information
Authors and Affiliations
Contributions
FB and MD designed, developed, and tested the calculation pipeline. HL helped with all force field and molecular dynamics related issues. CS and AZ conceived the study. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Competing interests
AZ is cofounder of GNWI—Gesellschaft für naturwissenschaftliche Informatik mbH, Dortmund, Germany.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Bänsch, F., Daniel, M., Lanig, H. et al. An automated calculation pipeline for differential pair interaction energies with molecular force fields using the Tinker Molecular Modeling Package. J Cheminform 16, 96 (2024). https://doi.org/10.1186/s13321024008905
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s13321024008905