halogen atoms are commonly 
 found in inorganic , organic , and pharmacological 
 molecules .
it has been reported that 50% of compounds 
 in high - throughput drug screening contains halogens .
halogen bonds , referring to the noncovalent interactions between 
 halogens and electron - rich atoms , are both strong and tunable . 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  . halogen atoms and some elements of groups
iv  vi 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 , ligand design 
 is also increasingly taking advantage of the halogen - bonding phenomenon .
for example , a number of recent 
 studies reported halogen - based hiv reverse transcriptase inhibitors . 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 . 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 .
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 .
similar treatments of halogen electrostatics and the -hole 
 effect have been reported previously . 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 .
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 . 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 .
the total energy is given by1the 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 . to derive 
 the amoeba polarizable 
 force field parameters , poltype was used to estimate initial atomic 
 multipoles at the mp2/6 - 311 g * * level of theory using stone s 
 original
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  ) was determined by matching to the ab initio qm molecular polarizability 
 tensor for ch3cl ( 5.3 , 3.8 , 3.8  ) .
the 
 same cl atomic polarizability was applied to mono- to tetrachloromethanes 
 ( table 1 ) .
the van der waals ( vdw ) parameters 
 of c and h have been transferred from the ch- in methanol and methylamine .
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 in conjunction with the tinker 6 package .
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 cc 
 or co 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 .
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 .
ccl4 , chcl3 , ch2cl2 , and ch3cl  as 
 pure liquids and as single molecules solvated in water . for the pure 
 liquid
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:2where 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 .
forcebalance allows 
 the fitting of selected force field parameters to qm cluster energies 
 and/or liquid thermodynamic properties using newton 
the hfe for each chloromethane 
 was computed using the conventional thermodynamic cycle consisting 
 of discharging , van der waals decoupling , and gas - phase recharging 
 steps .
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:3where 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 147 potential , as described in our previous 
 publications .
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 . for each of 
 the chloromethanes that exist as a liquid at room temperature 
 ( all but ch3cl ) , we simulated a solution consisting of 
 395882 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 .
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:4 in eq 4 , obs is the 
 experimental c = o vibrational frequency of 
 acetophenone dissolved in one of the three chloromethane solvents , fsolv 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 ) and are also calibrated against a set of reference solvents .
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 . , 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 .
these observations are recapitulated by the amoeba atomic multipole 
 moments derived from qm calculations . 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 . 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 ) .
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 
 cc distances for homodimers and different co 
 distances for heterodimers .
the forcebalance program was then used 
 to fit the vdw parameters to the gas - phase dimer energy .
( a ) ccl4ccl4 , ( b ) ccl4water , 
 ( c ) chcl3chcl3 , ( d ) chcl3water , 
 ( e ) ch2cl2ch2cl2 , ( f ) ch2cl2water , ( g ) ch3cl  ch3cl , ( h ) ch3cl  water .
carbon , 
 orange ; cl , green ; h , white ; o , red . amoeba and qm calculated pairwise interaction energy for dimers 
 in gas phase .
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 . 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 chcl3h2o 
 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 chcl3h2o 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 .
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 / cm 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 ) .
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 
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 .
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 ) .
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 .
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 .
density and heat of vaporization 
 data were used in the parameterization ; the hydration results were 
 not .
hvap ( kcal / mol ) , heat of vaporization ;  ( g / cm ) , 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 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 ) .
the 
 amoeba force field accurately predicts solvent electric fields , and 
 so , by extension , can describe solvent - induced frequency shifts and 
 band - broadening . 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 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 .
peak frequency and linewidth ( cm ) , electric field mean , and standard deviation ( mv / cm ) . electric field calculated from 
  
 = ( 0.484  0.029 ) fsolv + 1703.6 .
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 , 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 
it is noticed from figure 4f that the closest clo distance in ccl4 solution is in the range of 34  , which is much longer 
 than the c  cl covalent bond ( 1.8  )
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 .
( 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 ) clcl 
 distance rdfs for ccl4 and ch3cl in pure liquid 
 simulations .
( 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 / gas indicates the correlation coefficient between qm and mm simulations : 
 ( a ) ccl4ccl4 , rmses for liquid 
 and gas are 0.27 and 0.23 kcal / mol , respectively .
( b ) cclh2o , rmses for liquid and gas are 0.16 and 0.20 kcal / mol , respectively . 

( c ) ch3clch3cl , rmses for 
 liquid and gas are 0.07 and 0.13 kcal / mol , respectively .
( d ) ch3clh2o , 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 / gas values
, it can be seen that the final liquid - optimized parameter 
 set in general yields very good correlation with qm results , with r = 0.850.99 and rmse = 0.070.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 ch3clch3cl 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 .
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 ccl4ccl4 dimer 
 configuration a ( in figures 5 and 6 ) , the intermolecular clcl 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 ccl4ccl4 in figure 4e shows that there are very few pairs with clcl 
 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(co 
 distance  3.2  , similar to the gas - phase structure ) 
 and d2(co 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 . among all the liquid phase 
 results computed by amoeba ,
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 ccl4water heterodimer interaction energies 
 were not in strong accordance with qm ( figure 5b ; rliquid = 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 .
conformations found in 
 liquid - phase simulations for which the amoeba 
 model has large errors .
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 
, 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 / cm 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 .
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 .
however , 
 the observation that ccl4 s electric fields were 
 modeled qualitatively incorrectly ( table 4 ) 
 and the problems with ccl4h2o 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 .
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 . 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 .
 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
Computational Methods
Results
and Discussion
Conclusion