Modeling organochlorine compounds and the σ-hole effect using a polarizable multipole force field.

The charge distribution of halogen atoms on organochlorine compounds can be highly anisotropic and even display a so-called σ-hole, which leads to strong halogen bonds with electron donors. In this paper, we have systematically investigated a series of chloromethanes with one to four chloro substituents using a polarizable multipole-based molecular mechanics model. The atomic multipoles accurately reproduced the ab initio electrostatic potential around chloromethanes, including CCl4, which has a prominent σ-hole on the Cl atom. The van der Waals parameters for Cl were fitted to the experimental density and heat of vaporization. The calculated hydration free energy, solvent reaction fields, and interaction energies of several homo- and heterodimer of chloromethanes are in good agreement with experimental and ab initio data. This study suggests that sophisticated electrostatic models, such as polarizable atomic multipoles, are needed for accurate description of electrostatics in organochlorine compounds and halogen bonds, although further improvement is necessary for better transferability.

Introduction

Halogen atoms are commonly found in inorganic, organic, and pharmacological molecules.[1−3] It has been reported that ∼50% of compounds in high-throughput drug screening contains halogens.[4,5] Halogen bonds, referring to the noncovalent interactions between halogens and electron-rich atoms, are both strong and tunable.[6−8] When attached to a strong electron-withdrawing group, halogen atoms may display an electron-depleted region on the outermost portion of the molecular surface, which has been referred to as the “σ-hole”.[2,9−14] Halogen atoms and some elements of groups IV–VI[13,15,16] with a σ-hole can thus form strong interactions with negatively charged sites or electron-donor-rich groups (e.g., Lewis bases, π-electrons, and anions). Such strong interactions are often referred to as halogen bonds (X-bonds) and offer alternatives to other common classical interactions, such as hydrogen bonds (H-bonds), that play important roles in the supramolecular chemistry of biosystems and nanomaterials. In addition to areas such as crystal engineering and solid-state materials,[17−19] ligand design is also increasingly taking advantage of the halogen-bonding phenomenon.[3,20,21] For example, a number of recent studies reported halogen-based HIV reverse transcriptase inhibitors.[1,22−24] In another biochemical application, Ho and co-workers utilized bromine-substituted uracil to promote the assembly of four-stranded DNA junctions where the halogen atoms facilitated the noncovalent bonding by acting as electrophilic sites.[25,26]

Given the importance of halogens and halogen bonds, it is crucial to develop accurate molecular mechanics (MM) models to capture their electronic structure and molecular interactions. However, because of their complicated charge distribution and high polarizability, it is difficult to model the halogen atoms using atomic point charges with spherically symmetric potentials.[6,24,27,28] The partial negative charges typically assigned to halogen atoms make their electrostatic interactions with electron donors repulsive instead of attractive. To deal with such anisotropic charge distributions and the σ-hole effect, Jorgensen et al. introduced off-center charged sites to halogen atoms in the OPLS-AA force field to capture the halogen bonds, which led to improvement in predictions for the density, heat of vaporization, relative hydration, and binding free energy.[24,29] Similar treatments of halogen electrostatics and the σ-hole effect have been reported previously.[2,30−33]

Beyond the fixed-charge model, a more physically appealing approach to model the complicated anisotropic electronic structure of halogens is to incorporate higher-order multipole moments. Among the available methods, the atomic multipole optimized energetics for biomolecular applications (AMOEBA) force field, which features atomic-based multipole moments (up to quadrupole moments) and inducible dipole-based polarizability, ought to be a good platform for the development of a halogen model. The AMOEBA force field has been developed for water, common single-atom ions, common organic molecules, and proteins.[34−37] In addition, the importance of an explicit representation of electronic polarization in molecular mechanics models has been demonstrated by various groups over the past decade.[34,36,38−49]

In this paper, we investigated a polarizable multipole-based molecular mechanics model for a series of chlorine-substituted methanes, including CCl4 (carbon tetrachloride), CHCl3 (chloroform), CH2Cl2 (dichloromethane), and CH3Cl (chloromethane). These simple chlorine-containing compounds allow us to focus on chlorine itself before investigating more complex chlorine-containing drugs. The classical force field was parametrized by using a combination of ab initio quantum mechanics (QM) and experimental data. Transferability of the force field was tested by computing the chlorocompounds’ hydration free energies (HFE) and solvent reaction fields on a reference solute. Detailed analysis of the energetics of homo- and heterodimers using the resulting force field along with ab initio calculations gives further insight into the intermolecular interactions of halogen atoms. Overall, the calculated hydration free energies and solvent electric fields agreed with experiment more satisfactorily for the less substituted compounds, and tetrachloromethane’s properties were the most difficult to reproduce.

Computational Methods

AMOEBA Force Field

The AMOEBA model has been described in detail in previous publications.[34,42,44] The total energy is given byThe terms in 1 include the valence contributions corresponding, respectively, to the bonds, angles, bond-angle cross couplings, out-of-plane, and torsional energies. The long-range electrostatic interactions, including both permanent and polarizable components, are treated with the particle-mesh Ewald (PME) algorithm.

Gas Phase Calculations

To derive the AMOEBA polarizable force field parameters, POLTYPE was used to estimate initial atomic multipoles at the MP2/6-311G** level of theory using Stone’s original “distributed multipole analysis” (DMA) procedure.[35,50] The resulting multipole moments (with monopoles fixed) were then further refined by fitting to the electrostatic potential computed at the MP2/6-311++G (2d, 2p) level of theory to overcome the potential instability of directly applying DMA to bigger basis sets with diffuse functions. The atomic polarizability of Cl (2.5 Å3) was determined by matching to the ab initio QM molecular polarizability tensor for CH3Cl (5.3, 3.8, 3.8 Å3). The same Cl atomic polarizability was applied to mono- to tetrachloromethanes (Table 1).

Table 1

Calculated Molecular Polarizability (Å3)

	methods	αxx	αyy	αzz
CCl4	QMa	10.23	10.22	10.22
	MMb	10.47	10.47	10.47
CHCl3	QM	9.14	9.14	6.56
	MM	9.21	9.20	7.08
CH2Cl2	QM	7.99	5.75	5.15
	MM	7.83	6.13	5.54
CH3Cl	QM	5.29	3.82	3.82
	MM	5.32	4.07	4.07

QM results were obtained at the MP2/aug-cc-pVTZ level.

MM results are calculated from interactive atomic dipole polarizabilities.[34]

The van der Waals (vdW) parameters of C and H have been transferred from the CH- in methanol and methylamine.[37] The Cl vdW parameters were first determined by fitting to the ab initio potential energy profiles of molecular dimers evaluated at the MP2/aug-cc-pVQZ level of theory. The fitting procedure was carried out using the ForceBalance program[51] in conjunction with the TINKER 6 package.[52] All ab initio calculations in this study were performed using the Gaussian09 program.[53] The geometry of each dimer pair was optimized using MP2/aug-cc-pVTZ in gas phase. The potential energy profile for each pair was then generated by displacing the molecules along the C···C or C···O axis. Single-point energy calculations at the MP2/aug-cc-pVQZ level of theory with basis set superposition error (BSSE) correction were applied to obtain the interaction energy for each generated structure. Because the dispersive energy may contribute significantly to the interaction energy of chloride compounds, the quality of MP2/aug-cc-pVQZ interaction energy was also examined by comparing with the CCSD(T)/CBS result for a CH3Cl-water dimer. The CCSD(T)/CBS interaction energy was estimated following the extrapolation approach of Hobza and co-workers.[54] The MP2 correlation energy was extracted to CBS from aug-cc-pVTZ and aug-cc-pVQZ, and then the correlation energy was adjusted from the MP2 to the CCSD(T) level by adding the difference between the two methods using the aug-cc-pVDZ basis set. The calculated MP2/aug-cc-pVQZ and CCSD(T)/CBS dimer interaction energy values are −1.071 and −1.100 kcal/mol, respectively. In addition, for the purpose of calibrating the classical model, the HF and MP2 energies seemed to converge reasonably well at the aug-cc-pVQZ basis set by comparing to the CBS results (see Supporting Information (SI) Table S1). Thus, MP2/aug-cc-pVQZ has been used to obtain all the gas-phase dimer energies in this study.

Liquid Calculations

Molecular dynamics (MD) simulations were performed on all four molecules—CCl4, CHCl3, CH2Cl2, and CH3Cl—as pure liquids and as single molecules solvated in water. For the pure liquid, a box of 216 chloromethane molecules was used. Liquid density and heat of vaporization were evaluated. The liquid simulations were first equilibrated through a 50 ps trajectory in the NPT ensemble at 298 K and 1 bar, followed by another 500 ps production run for data collection. To estimate the heat of vaporization, we used the equation below:[34]where Eliq is the averaged potential energy of single molecule in the liquid box and Egas is the energy of the single molecule in the gas phase and is calculated by running a 500 ps simulation for one molecule using a 0.1 fs time step. By computing the liquid density and heat of vaporization for each chloromethane from molecular dynamics simulations, ForceBalance was applied to iteratively optimize the vdW parameters for Cl.[51,55] ForceBalance allows the fitting of selected force field parameters to QM cluster energies and/or liquid thermodynamic properties using Newton–Raphson and other optimization algorithms.

The HFE for each chloromethane was computed using the conventional thermodynamic cycle consisting of discharging, van der Waals decoupling, and gas-phase recharging steps.[42,48] The charging/discharging step handles the polarization effect in the AMOEBA force field through the scaling of solute atomic polarizability. Thus, the final HFE can be expressed as:where ΔAdischarging(aq) and ΔAdecoupling(aq) are the free energy changes due to turning off the electrostatic and vdW between the solute and environment respectively, and ΔArecharging(vac) corresponds to the intramolecular electrostatic interactions in vacuum. The vdW annihilation used a soft-core buffered 14–7 potential, as described in our previous publications.[42,45] The discharging schedules are λ = 0, 0.1, 0.2, ... 1; decoupling schedules are λ = 0.0, 0.2, 0.4, 0.5, 0.55, 0.6, 0.625, 0.65, 0.7, 0.75, 0.8, 0.9, 0.95. The Bennett acceptance ratio (BAR) method was used to evaluate the free energy changes.[56]

For each of the chloromethanes that exist as a liquid at room temperature (all but CH3Cl), we simulated a solution consisting of 395–882 chloromethane molecules (to fill a 45 Å cubic box) and 1 acetophenone molecule and calculated the electric field the solvent exerts onto the C=O bond of acetophenone using methods previously described.[57] The solvent box was equilibrated for 100 ps in an NPT ensemble, and production dynamics were carried out for 500 further ps, during which the solvent field was calculated every 10 fs. The calculated electric fields can be compared with experimental values evaluated from vibrational frequencies through the linear vibrational Stark equation:In eq 4, ν̅obs is the experimental C=O vibrational frequency of acetophenone dissolved in one of the three chloromethane solvents, F⃗solv is the electric field the chloromethane solvent exerts onto the C=O bond, ν̅0 is a reference frequency associated with zero-electric field, and Δμ⃗probe is the C=O vibration’s Stark tuning rate. ν̅0 and Δμ⃗probe correspond to the vibration’s gas-phase frequency and difference dipole moment (measured in Stark spectroscopy)[58] and are also calibrated against a set of reference solvents.

Results and Discussion

Gas-Phase Study

To examine the electronic structure of Cl atoms and the σ-hole effect, we evaluated the electrostatic potential surface of each of the four chloromethanes (Figure 1). As suggested by Scholfield et al.,[2] the electron-withdrawing ability of the atoms and functional groups that are bonded to Cl affect the size of the σ-hole on Cl, which is consistent with our QM calculations. As shown in Figure 1d, in CH3Cl, the electrostatic potential around Cl in CH3Cl is mostly negative (red), but the outermost portion of Cl’s surface facing away from the central carbon is actually neutral (green). This is in accordance with the fact that Cl is more electronegative than C and, thus, bears a partial negative charge in a C–Cl bond. However, the covalent bond pulls the electron distribution toward the center of the bond and leaves a neutral patch on the outermost surface of Cl along the C–Cl axis. From CH2Cl2 to CCl4 (Figures 1), this “patch” becomes larger and increasingly positive (blue) as the electron-withdrawing forces generated by additional Cl atoms increase. The Cl on CCl4 has the largest σ-hole. These observations are recapitulated by the AMOEBA atomic multipole moments derived from QM calculations.

Figure 1

Ab initio molecular electrostatic potential surfaces calculated at the MP2/6-311G++ (2d, 2p) level. (a) CCl4, (b) CHCl3, (c) CH2Cl2, (d) CH3Cl. The electrostatic potential is mapped on the surface of molecular electron density at 0.001 au contours. Coloring scheme: red (<−12.55 kcal/mol), yellow (−5.02 kcal/mol), green (0 kcal/mol), light blue (5.02 kcal/mol), and blue (>12.55 kcal/mol).

In Table 2, the calculated electrostatic charges (monopoles) indicate that net charge on Cl become less negative going from mono- to tetrachloromethane. Furthermore, the z component of the dipole moment on Cl (pointing along the Cl–C bond vector) decreases from CH3Cl to CCl4, suggesting diminishing charge separation along the Cl–C bond, while the positive σ-hole gets larger. In addition, there are large quadrupole components on the Cl atoms, relative to those of C and H (see the parameters in the SI). The small root-mean-square difference between the electrostatic potentials calculated from QM and atomic multipoles (Table 2) suggests the distributed atomic multipole moments are able to represent these complex charge distributions. In Table 1, QM and force field calculated molecular polarizability values show very good agreement across the four compounds. The atomic polarizability values for C and H were transferred from the existing AMOEBA parameter set; the same Cl atomic polarizability was used for all four compounds here (all parameters are given in the SI). The Cl atomic polarizability is almost twice that of the C atom, and the molecular polarizability steadily increases from CH3Cl to CCl4.

Table 2

Electrostatic Parameters for Cl in Each of the Chlorocompounds

		Cl multipole (au)a	RMSE of electrostatic potential (kcal/mol)b
CCl4	monopole	–0.07475	0.1270
	dipole	0.00000 0.00000 0.01306
	quadrupole	–0.68662
		0.00000 −0.68672
		0.00000 0.00000 1.37334
CHCl3	monopole	–0.13452	0.2189
	dipole	0.00000 0.00000 0.06939
	quadrupole	–0.80136
		0.00000 −0.79348
		0.00000 0.00000 1.59484
CH2Cl2	monopole	–0.18091	0.0741
	dipole	0.00000 0.00000 0.10467
	quadrupole	–0.79232
		0.00000 −0.90414
		0.00000 0.00000 1.69646
CH3Cl	monopole	–0.22443	0.0665
	dipole	0.00000 0.00000 0.17165
	quadrupole	–0.88958
		0.00000 −0.88997
		0.00000 0.00000 1.77955

The atomic multipoles were derived from ab initio calculations at the MP2/6-311++G(2d,2p) level.

The RMSE is based on a comparison to the ab initio electrostatic potential at the same level.

We examined a series of homodimers of CCl4, CHCl3, CH2Cl2, and CHCl3, as well as heterodimers with a water molecule (Figure 2) in which the geometries were randomly chosen near local minima on the dimer energy surface. Because POLTYPE already produced the electrostatic and valence parameters based on the QM calculations on the monomers, the QM calculations on dimers allow us to quickly estimate the vdW parameters (each chloromethane has its own Cl parameters). The use of heterodimers with water ensures that the parameters are transferable to another environment. Geometry optimizations of these dimers were carried out at the MP2/aug-cc-pVTZ level. The BSSE-corrected association energy at the MP2/aug-cc-pVQZ level was obtained at different C···C distances for homodimers and different C···O distances for heterodimers. The ForceBalance program was then used to fit the vdW parameters to the gas-phase dimer energy.[51]

Figure 2

Homo- and heterodimer geometries in gas phase. (a) CCl4–CCl4, (b) CCl4–water, (c) CHCl3–CHCl3, (d) CHCl3–water, (e) CH2Cl2–CH2Cl2, (f) CH2Cl2–water, (g) CH3Cl–CH3Cl, (h) CH3Cl–water. Carbon, orange; Cl, green; H, white; O, red.

By optimizing the vdW parameters, AMEOBA was able to reproduce the ab initio interaction energies of homodimers reasonably well (Figure 3). However, for heterodimers, the agreement is generally worse, particularly for the CHCl3–H2O dimer. The difference between the QM and AMOEBA interaction energy is ∼1.0 kcal/mol at the minimum energy distance. In Figure 3, it was also noticed that neither vdW parameter d nor ε appears to display a consistent “chemical” trend going from CCl4 to CH3Cl, as observed for the electrostatic parameters (Table 1 and Table 2). As noted earlier, these configurations were randomly chosen and might not represent the most important configurations in the liquid phase. For example, in the CHCl3–H2O dimer (Figure 2d), the H atom instead of Cl atom of CHCl3 is facing the water O atom. Although the dimer data is likely insufficient for determining the final vdW parameters, these simple gas-phase calculations provide a set of starting parameters for subsequent examination of liquid state properties using molecular dynamics simulations.

Figure 3

AMOEBA and QM calculated pairwise interaction energy for dimers in gas phase. Energy is in kcal/mol, and distance, in Å. The Cl’s vdW for CCl4: d (diameter) = 3.62 Å, ε = 0.4 kcal/mol, and RMSE (QM vs AMEOBA) = 0.51 kcal/mol. The Cl’s vdW for CHCl3: d = 3.06 Å, ε = 1.62 kcal/mol, and RMSE = 1.43 kcal/mol. The Cl’s vdW for CH2Cl2: d = 3.49 Å, ε = 0.5 kcal/mol, and RMSE = 1.02 kcal/mol. The Cl’s vdW for CH3Cl: d = 3.67 Å, ε = 0.20 kcal/mol, and RMSE = 1.24 kcal/mol.

Liquid-Phase Simulations

The gas-phase QM study is important in the sense that it provides detailed information on electrostatic and intermolecular interactions in an isolated environment, although it is challenging to directly compare these results to experimental measurements. We next optimized the gas-phase QM derived Cl’s vdW parameters for each of the chloromethanes using the experimental liquid density and heat of vaporization data and the ForceBalance method. Given that there are two free parameters, it is not surprising that the final AMOEBA results are in excellent agreement with experimental values for both density and heat of vaporization. From initial optimization, it was found that the resulting Cl diameters, d, across all four chloromethane are very similar. Thus, we decided to use the same d parameter and reoptimized the vdW energy depth parameter ε for each Cl type in the four chloromethanes. The RMSE is 0.01 g/cm3 for the density and 0.12 kcal/mol for the heat of vaporization. Interestingly, the resulting optimal ε parameter systematically decreases from mono- to tetrachloromethane (Table 3). Recall that the Cl charge and dipole moment also follow a similar trend (Table 2).

Table 3

Final vdW Parameters Fitted from Liquid-Phase Simulation and Comparison of Calculated Liquid Properties to Experimenta,b

	CCl4	CHCl3	CH2Cl2	CH3Cl
Cl, d (Å)	3.898	3.898	3.898	3.898
Cl, ε (kcal/mol)	0.319	0.340	0.362	0.413
ΔHvap_cal	7.71 (±0.02)	7.35 (±0.01)	6.84 (±0.01)	4.34 (±0.01)
ΔHvap_expt	7.74[65]	7.50c	6.82[66]	4.52[65]
ρcalc	1.593 (±0.017)	1.490 (±0.017)	1.330 (±0.017)	0.920 (±0.019)
ρexpt	1.584[65]	1.480[65]	1.327[67],c	0.911[65]
ΔAhyd_cal	1.99	–0.49	–0.73	–0.26
ΔAhyd_expt[68,69]	0.08	–1.08	–1.31	–0.55

Density and heat of vaporization data were used in the parameterization; the hydration results were not.

ΔHvap (kcal/mol), heat of vaporization; ρ (g/cm3), liquid density at room temperature; ΔAhyd (kcal/mol), hydration free energy; subscript with “cal”, calculated values; subscript with “expt”, experimental reference.

This value is measured at 293 K[70]

In addition to the neat liquid simulation, we evaluated the hydration free energy of each of the chloromethanes to examine the transferability of the resulting parameters and the potential energy model. The HFE of small molecules is an important physical property in many chemical and biological processes, such as protein–ligand binding.[59] In biological force field development, the calculation of HFE is considered as a critical validation step. Extensive efforts have been made to improve the fixed-charge force fields to reduce the HFE error of small organic molecules to ∼1.0 kcal/mol.[42,60−64] Previously, we have shown that for a group of common organic molecules, the HFE calculated by AMOEBA is in excellent agreement with experimental data (RMSE ∼ 0.4 kcal/mol).[37] The computed HFE for each chloromethane using the alchemical free energy methods as described in the Computational Methods section is shown in Table 3. The errors are well within 1 kcal/mol for mono-, di-, and trichloromethane, indicating reasonable transferability of the model from neat liquid to a water environment. However, the error for tetrachloromethanes is 1.91 kcal/mol. Overall, there is a clear trend of increasing error from mono- to tetrachloromethanes. The large error suggests that the model for tetrachloromethanes needs to be further examined for the cause of the transferability issues.

Solvent Fields

Previous work has demonstrated that a carbonyl vibration’s frequency reports on the local electric field created by the surrounding solvent according to a simple linear model (eq 4).[71] The AMOEBA force field accurately predicts solvent electric fields, and so, by extension, can describe solvent-induced frequency shifts and band-broadening.[57] Drawing on this concept, we evaluated the electric fields that the chloromethane solvents exert on the test solute acetophenone by recording the infrared spectrum of acetophenone dissolved in CH2Cl2, CHCl3, and CCl4 (10 mM) and mapping the peak frequency and linewidth to the mean electric field and electric field standard deviation using linear models calibrated against seven nonhalogen-containing reference solvents. The results are shown in Table 4.

Table 4

Comparison of Simulated Electric Fields to Experiment

	CCl4	CHCl3	CH2Cl2
peak frequencya	1691.2	1683.3d	1684.6d
expected mean fieldb	–25.6 ± 1.4	–41.9 ± 2.3	–39.2 ± 2.2
simulated mean field	–11.4 ± 0.2	–36 ± 2	–33 ± 1
linewidtha	8.9	16.1d	11.4d
expected field std devc	10.9 ± 1.1	20.9 ± 2.0	14.4 ± 1.5
simulated field std dev	5.8	24	15

Peak frequency and linewidth (cm–1), electric field mean, and standard deviation (MV/cm).

Electric field calculated from ν̅ = (0.484 ± 0.029) Fsolv + 1703.6.

Field standard deviation calculated from LW = (0.714 ± 0.08)σF + 1.14.[57]

Data from Fried et al.[72]

The average electric fields for the chloromethane solvents are significantly greater than expected from continuum models (such as the Onsager reaction field or Poisson–Boltzmann equation), implying the presence of H-bonds and X-bonds between chloromethane molecules and the C=O vibrational probe. The AMOEBA model reproduces the average solvent field exerted by CH2Cl2 and CHCl3 reasonably well, including the small difference between them. Agreement among simulations and experimental values for the dispersion of these two solvents’ electric fields is similarly strong. The good agreement between AMOEBA and experimentally determined electric fields suggests the electrostatic moments assigned to these molecules, derived from gas phase calculation, give a correct description of their electrostatic interactions in the condensed phase, as well.

CCl4 exerts smaller but nonetheless significant electric fields, although MD simulations substantially underestimate them. Because CCl4 possesses no H-bonding capacity but is expected to donate the strongest X-bonds among the solvents studied, the increased discrepancy is consistent with the possibility that the AMOEBA model renders electrostatics of X-bonds insufficiently attractive.

Transferability of Classical Model

The general philosophy of the AMOEBA model is a classical potential energy function applicable in different chemical and physical environments. The introduction of the inducible atomic dipole, which allows the electrostatics of a molecule to respond to its instantaneous environment, is a step forward in the development of a more transferable force field. In the past, for many ionic and organic molecular systems,[35,37,73−76] we were able to achieve such transferability by iteratively optimizing the vdW parameters between the gas-phase cluster energy and liquid experimental properties. In this study, the electrostatic parameters were derived from QM; vdW parameters were initially derived from QM dimer energetics and optimized in liquid simulations.

To examine the transferability of the final parameters back to the gas-phase environment, about 20 new dimer configurations were extracted from the neat liquid and hydration free energy simulation trajectories of CCl4 and CH3Cl, respectively. By computing the radial distribution function (RDF) for relevant atom pairs (Figure 4a to d), we were able to identify the first solvation shell, from which the homo- and heterodimers were chosen randomly for the subsequent QM–force field comparisons. It is noticed from Figure 4f that the closest Cl···O distance in CCl4 solution is in the range of 3–4 Å, which is much longer than the C–Cl covalent bond (1.8 Å). Similar contact distances between the Cl···O have been reported for chlorobenzene.[24] This observation, together with the evidence that our clasical model can well reproduce the dimer interaction energy as shown next (Figure 5), indicates that the so-called halogen bonds, formed by the positive σ-hole on Cl and electron rich O, are of noncovalent nature and can be modeled by a combination of electrostatic and vdW inteaction as can most hydrogen bonds.

Figure 4

Radial distribution functions from liquid-phase simulation. (a) CCl4 neat liquid simulation, 20 neighboring CCl4 molecules were found around the reference CCl4 within 11 Å (second peak). (b) CCl4 in solvent, with 21 water neighboring molecues around within 10 Å. (c) CH3Cl liquid simulation, 21 neighboring CH3Cl molecules were found within 7.6 Å. (d) CH3Cl in solvent, with 21 neighboring water molecules found within 5.4 Å. (e) Cl···Cl distance RDFs for CCl4 and CH3Cl in pure liquid simulations. (f) Cl···O distance RDFs for CCl4 and CH3Cl in solvent.

Figure 5

Correlations between QM and AMOEBA dimer interaction energies. The label “liquid” (open circles) refers to results from the vdW parameter set optimized using the liquid properties, and “gas” (stars) refers to the results from the parameter set determined from QM dimer energy profile in gas phase (Figure 3). Rliquid/gas2 indicates the correlation coefficient between QM and MM simulations: (a) CCl4···CCl4, RMSEs for liquid and gas are 0.27 and 0.23 kcal/mol, respectively. (b) CCl···H2O, RMSEs for liquid and gas are 0.16 and 0.20 kcal/mol, respectively. (c) CH3Cl···CH3Cl, RMSEs for liquid and gas are 0.07 and 0.13 kcal/mol, respectively. (d) CH3Cl···H2O, RMSEs for liquid and gas are 0.13 and 0.22 kcal/mol, respectively.

For a total of 83 dimers, we evaluated the association energy by using both QM (BSSE corrected MP2/aug-cc-pVQZ) and the AMOEBA force field. The correlation between the QM and force field results using the initial gas-phase determined vdW parameters and final, liquid-optimized vdW parameters is shown in Figure 5.

By inspecting the Rliquid/gas2 values, it can be seen that the final liquid-optimized parameter set in general yields very good correlation with QM results, with R2 = 0.85–0.99 and RMSE = 0.07–0.27 kcal/mol. In addition, the liquid-optimized parameters did better than the gas-phase QM energy-derived parameter set, especially for the CH3Cl-CH3Cl dimer. During the initial vdW parametrization based on the gas-phase QM dimer energy (Figure 3), one single conformation was chosen randomly for each dimer so that it is possible that we used a local minimum-energy configuration for CH3Cl···CH3Cl that is uncommon in the liquid state.

We further examined a few structures (Figure 6A, B1, B2, B3, D1, and D2) that correspond to the outliers labeled in Figure 5. Most of these “outlier” structures observed in liquid-phase simulations (Figure 6 B1, B2, B3, D1, and D2) are very different from those used in the gas-phase dimer vdW parametrizations (Figure 2). In CCl4···CCl4 dimer configuration A (in Figures 5 and 6), the intermolecular Cl···Cl distance (3.4 Å) is similar to that of the gas-phase QM optimized structure (3.3 Å, Figure 2a). On the other hand, the RDF of CCl4···CCl4 in Figure 4e shows that there are very few pairs with Cl···Cl distances less than 3.5 Å in the liquid state. In other words, configuration A is a rare case in the liquid state but close to the optimal structure in the gas phase. It is therefore understandable that the liquid-optimized vdW parameter set did worse on this structure than the parameter set directly fitted to the gas-phase dimer energy profile. A similar conclusion can be drawn by comparing D1(C···O distance ∼ 3.2 Å, similar to the gas-phase structure) and D2(C···O distance ∼ 4.5 Å, prevalent in liquid simulation). For the B1, B2, and B3 configurations extracted from liquids, the H atom of water is facing the Cl of CCl4, and in the gas-phase dimer we used above (Figure 2b), the O atom of water is facing the Cl. The overall worse performance of the gas-phase parameters on these structures is a reults of the inferior transferability when only one configuration is used in the preliminary parametrization.

Figure 6

Conformations found in liquid-phase simulations for which the AMOEBA model has large errors. The labels of the structures correspond to those in Figure 5.

Among all the liquid phase results computed by AMOEBA, CCl4’s HFE and solvent field have the greatest errors. The root of this poor agreement likely lies in the force field’s inaccurate description of CCl4’s X-bonding to the O atoms of water and acetophenone. The dielectric constant of CCl4 (2.2) is close to that of n-hexane (1.9), yet it exerts more than twice the electric field (−25 vs −11 MV/cm) on acetophenone. X-bonding interactions are responsible for the additional electrostatic attraction, yet this effect is clearly not captured in our electric field simulations (Table 4); more subtly, it might also explain why the calculated HFE of CCl4 was overly endergonic (Table 3) and why CCl4–water heterodimer interaction energies were not in strong accordance with QM (Figure 5b; Rliquid2 = 0.85). Strong agreement with CCl4’s experimental heat of vaporization, density (Table 3), and homodimer interaction energies (Figure 5a) may have been possible because CCl4 does not have a good X-bond acceptor itself, so the effect is muted in neat systems.

Conclusion

Given the abundance of organochlorine compounds—or halogen compounds in general—in areas of chemistry, biology, and drug discovery, it is important to have an accurate classical force field for modeling such compounds and their interaction with other molecules. However, the traditional fixed atomic charge force field is inadequate for treating the complex electrostatics in halogen compounds, especially when the “σ-hole” effect is strong.[24]

In this study, we applied the AMOEBA model to investigate a series of chloromethanes, including CCl4, CHCl3, CH2Cl2, and CH3Cl. The atomic multipole framework is a natural choice for complicated charge distribution seen on the halogen atoms. The electrostatic parameters were derived from QM calculations of the monomers, and the vdW parameters were optimized against liquid densities and heats of vaporization. The electrostatic potentials of monomers were well described by atomic multipoles, including the large charge separation observed on the Cl atom. The QM gas-phase dimer energies were reasonably reproduced by AMOEBA for homodimers, but less satisfactory for heterodimers formed with a water molecule. The neat liquid density and heat of vaporization, which were used as the vdW parameter optimization targets, were well reproduced by AMOEBA simulations. The RMSE is 0.01 g/cm3 for the simulated density and 0.12 kcal/mol in heat of vaporization. The hydration free energy and solvent field, which were not used in the parametrization process, agree quite well for the less substituted chloromethanes; however, the predicted tetrachloromethane HFE showed a larger deviation from the experimental data (1.91 kcal/mol). By examining a number of dimer configurations from the liquid state simulations, we show that, in general, the AMOEBA interaction energies are very well correlated with the QM results. The chloromethane–H2O heterodimer interaction energies suggest that X-bonds, formed between the positive σ-hole on Cl and the negative charge on O, can be reliably treated by a combination of electrostatic and vdW interactions, just like H-bonds.[77] However, the observation that CCl4’s electric fields were modeled qualitatively incorrectly (Table 4) and the problems with CCl4–H2O interaction energies probably reflect deeper issues with the AMOEBA potential function rather than the parametrization. In particular, the problems mentioned hint at a cancellation of error in which electrostatic interactions are made insufficiently attractive while vdW interactions are under-repulsive.

In diagnosing this energy decomposition problem, we point out that Cl atom has much larger vdW parameters (diameter and energy well depth) than the other elements (C, O, and H) in this study. The vdW interaction energys for unlike pairs (e.g., Cl–H or Cl–O) are computed by using simple combining rules, which are known to be problematic.[78,79] This problem would be more severe when multiple Cl atoms are in close contact with a very dissimilar atom such as O or H in water. On the other hand, it has been shown that the electrostatic energy by point multipoles (or point charges) tends to underestimate the electrostatic attraction for large diffuse molecular species.[80−300] Moreover, an inadequate vdW combination rule would prevent molecules from adopting the correct spatial orientations, and electric fields (especially from quadrupole moments, such as on Cl) depend very sensitively on distances. In developing force fields that can correctly decompose interaction energies into electrostatic and vdW components (and not rely on error cancelation), training data that encode information about specific intermolecular interactions is particularly incisive. To this end, vibrational Stark shifts and heterodimer interaction energies appeared to be more discriminating than the other properties analyzed.

To the best of our knowledge, this study represents the first example of using vibrational Stark measurements to assess the parametrization of a force field, and we expect they will be of further use moving forward, since they isolate the effect of electrostatic interactions from other intermolecular forces. Nonetheless, the current parametrization procedure reported herein appears satisfactory for modeling most organochlorine compounds and can be extended to chlorine-containing molecules of interest.