The quadrupole correction: from molecular electrostatic potential to

Subscriber access provided by University of Winnipeg Library

Molecular Mechanics

The quadrupole correction: from molecular electrostatic potential to free energies of halogen bonding Oleg Igorevich Titov, Dmitry A. Shulga, and Vladimir A. Palyulin J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00739 • Publication Date (Web): 21 Dec 2018 Downloaded from http://pubs.acs.org on December 22, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

The quadrupole correction: from molecular electrostatic potential to free energies of halogen bonding Oleg I. Titov, Dmitry A. Shulga, Vladimir A. Palyulin* Department of Chemistry, Lomonosov Moscow State University, Moscow, 119991 Russia Keywords: halogen bonding, electrostatic interactions, drug design, force fields, free energy

ABSTRACT

Recently several molecular mechanics models of halogen bonding have been published. They describe the electrostatic potential anisotropy near the heavy halogen atom (known as a σ-hole) in different ways, ranging from an all-atom multipole expansion to a single positive extra-point charge. However, the question of a reasonable balance between the accuracy and the simplicity of the model remains open. In this paper, we introduce the simplistic RESPQ electrostatics model built on the RESP charges complemented with fixed atomic quadrupoles. We show that it: 1) correctly describes the MEP anisotropy of aromatic halogen atoms, 2) improves the description of the halogen-water interaction energies in both halogen and hydrogen bonding cases, 3) provides an excellent estimations of solvation free energy differences of aromatic halogens, 4) is compatible (with the help of multipole charge cluster approximation) with contemporary molecular modeling packages.

1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Page 2 of 37

Introduction Over the last several years halogen bonding (XB), an attraction between the halogen atom and the

electron density donor, has evolved from a novel counterintuitive subtype of multipolar interactions into a recognizable and well-known interaction pattern.1-9 It has found its application in different areas of chemistry ranging from new materials to drug design.1, 4, 10-14 The origin of the halogen bonding, the σ-hole effect,2, 5 has been studied and generalized for the other heavy elements such as chalcogens and pnictogens.2, 15-17 This effect cannot be described by the simplistic atomic charges model,16 thus several molecular mechanics models have been proposed to describe the characteristic anisotropy of the halogen atoms most of the which18-24 are based on the extra-point charge approximation, where an additional positive charge is placed on the continuation of the C-Hal bond to represent the electron deficient area formed by the σ-hole, while the others25-29 take advantage of multipole expansion as the natural expansion of atomic charges. Some of them require a complete charge refit for every molecule in question,18, 20, 26, 28 some models use tabulated values or other fast schemes.22-24, 27 Earlier, we have published an approach describing halogen bonding by means of multipole expansion.30, 31 We have shown that the main features of anisotropy of halogen atoms can be described by complementing the atomic charges with the fixed symmetric atomic quadrupole placed on a halogen atom.31 Being fixed for the atomic types, these quadrupoles bring no additional degrees of freedom to the electrostatic model and thus seem to be very promising for the incorporation into the modeling of the halogen bonding, which is harder to achieve with extra-point charges. The known halogen bonding models differ in complexity significantly from a simple tabulated extra-point charge incorporation into a classical force field to a full reparametrization of the halogen atom parameters in a multipole-based polarizable force field. On one hand, the simple extra-pointbased models tend to focus solely on the halogen bonding description not paying specific attention to 2 ACS Paragon Plus Environment

Page 3 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


X-H hydrogen bonding features. They also tend to have the extra-point charge placed as far as 2.2 Å away from the halogen atom which may lead to the ESP fitting problems (such as being too close to the reference grid or “burying” the halogen atom)31 and numerical issues during the simulation as this charge approaches the van der Waals boundaries of the halogen. On the other hand, the complicated polarizable models give extremely accurate free energy results, though they are not used as widely as classical force fields. Additionally, the complexity of such models and the numerous parameters introduced make them prone to overfitting and poor transferability. Thus, in this work, we are aiming at creating a very simple model with minimal force field modifications which would give the most significant improvement to the free energy estimations. Previously we have found31 that fixed symmetric atomic quadrupoles describe MEP anisotropy of the halogen atom well enough (MEP RMSE is less than 0.7 kcal/mol). This simplistic and extremely transferable model opens the way for the accurate modeling of halogen interactions in drug-like systems. In this work, we are testing this model in applied molecular modeling by calculating solvation free energy difference with a classical GAFF force field.32 Halobenzene – water system was chosen for this test. On one hand, halobenzenes are the simplest molecules containing aromatic halogen atoms which are of particular interest for the medicinal chemistry. On the other hand, a water molecule, like halobenzenes, can participate in both halogen and hydrogen bonding (HB), meaning that correct description of halogen-water interaction requires a correct and consistent description of both of these interaction patterns, especially in the view that hydrogen bonding description within GAFF has been well tested by numerous molecular dynamics modeling experiments. Needless to say, water is present in abundance in any biological system and correct description of solvation is one of the crucial factors in the modeling of such systems. With this system, we would like to ensure that both HB and XB are described accurately and consistently. In what follows we, first, introduce two anisotropic electrostatic models and compare their performance with the atomic charge only model. Second, we show how the introduction of halogen 3 ACS Paragon Plus Environment


Page 4 of 37

atom anisotropy changes molecular mechanics interaction of the studied halobenzenes with a water molecule. And finally, we analyze the performance of the anisotropic model in the reproduction of the free energy difference of hydration.

2

Methods

2.1 Atomic quadrupoles and their representation in molecular dynamics Three electrostatic models were compared. The first one is the RESP charges, which are recommended by GAFF force field developers. These charges, along with AM1-BCC33 ones are the default choice for small molecules in AMBER family force fields. They serve as a baseline isotropic model for our comparison. The second approach, mdzq, takes advantage of the multicenter multipolar expansion by optimizing atomic dipole and quadrupole placed on the halogen atom along with atomic charges to fit the MEP provided. It is the most flexible one and demonstrates the lowest MEP error on this particular grid, however, our previous research indicates that this kind of approximation should not be needed. This approach provides the most detailed MEP and serves as the example of the unreasonably extended solution to the XB description problem thus effectively bracketing the optimal model. The last series of models, RESPQ, combines the previous approaches and serves as the proposed way of describing halogen anisotropy: RESP charges are used with fixed symmetric atomic quadrupoles (1.15, 1.59, and 2.76 a.u. for Cl, Br, and I respectively), placed on the halogen atom. These values were obtained earlier31 as the medians of quadrupole distribution of a large set of various disubstituted halobenzenes and thus are good ‘average’ quadrupoles. This model has no additional degrees of freedom compared to the RESP charges since these values are fixed and never changed, so the correction can be easily transferred to any molecule. All the electrostatic models discussed in this paper are based on the reference RHF/6-31G* MEP. The molecular geometry of halobenzenes was generated with the help of OpenBabel software34 (with --gen3D option). The reference MEP was estimated on 5 layers (1.6, 1.8, 2.0, 2.2, and 2.4 times 4 ACS Paragon Plus Environment

Page 5 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


larger than the van der Waals radius) of Connolly surfaces with 1.12 Bohr-2 point density at RHF/631G* level of theory in Firefly v 7.1G software.35 This quantum chemistry approximation is a de facto standard for electrostatic models of AMBER-based classical force fields without explicit polarization treatment.32, 33 The point density was increased from the standard one for better sampling. RESP charges were obtained with ‘resp’ program from the AmberTools version 1.536 package. Atomic multipoles were fit directly to the reference MEP by the same procedure as in our previous work31 with the help of Electrostatic Tools package v 0.4.0.37 For each halogen atom a local coordinate frame was introduced based on the local nuclei symmetry: the Z-axis was directed along C-Hal bond, the X-axis was located in the plane of the aromatic ring, while Y-axis was perpendicular to it. The atomic multipoles were fit to the reference MEP with their principal axes constrained to this chosen local frame to reduce numerical noise arising from irregularities of the reference grid and molecular geometries.

2.2 From atomic multipoles to multipole charge cluster Most of the modern molecular dynamics packages cannot work with atomic multipoles. In order for the model to be accessible to a wide range of researchers, one needs to convert the multipoles to extrapoint charges. There are several procedures described in the literature,38, 39 however they either make little use of molecular symmetry or use rather complicated algorithms. We have developed our own simple algorithm to convert atomic multipoles to multipole charge clusters (MCC).30 In our approach, the monopole, dipole, and quadrupole can be represented as a set of 1 to 7 atomic charges (Figure 1). To achieve this, a regular octahedron of six extra-point charges was added on every halogen atom.



Page 6 of 37

Figure 1. Extra-point charges added to the bromobenzene molecule to represent atomic dipole and quadrupole (MCC radius is increased to 1.0 Å) with local multipole frame axes directions labeled. The RESP model uses no extra-points, mdzq uses ±Z and ±X, while RESPQ needs only ±Z. The electrostatic potential of the atom associated with an atomic monopole (q), an atomic dipole (with components dx, dy, and dz), and a symmetric atomic quadrupole (with diagonal components Qxx, Qyy, and Qzz) for the purposes of molecular modeling can be represented as the potential of a cluster of charges with a central charge (q0) located at the atomic nucleus and an octahedron of six additional charges located at the distance r along the coordinate frame axes (q+x, q-x, q+y, q-y, q+z, q-z), as shown in Figure 1. Using the definition of dipole and traceless quadrupole moments one can write down 6 equations on these charge values (7 variables): 1 0 0 0

1 𝑟 0 0

1 −𝑟 0 0

0

𝑟2

𝑟2

1 2 0 − 𝑟 ( 2

1 − 𝑟2 2

1 0 𝑟 0 1 − 𝑟2 2

1 0 −𝑟 0 1 2 − 𝑟 2

𝑟2

𝑟2

1 0 0 𝑟 1 2 − 𝑟 2 1 − 𝑟2 2

1 𝑞0 𝑞 0 𝑞+𝑥 𝑑𝑥 0 𝑞−𝑥 𝑑 𝑦 −𝑟 𝑞+𝑦 = 1 2 𝑑𝑧 𝑞−𝑦 − 𝑟 𝑄𝑥𝑥 2 𝑞+𝑧 1 2 (𝑄𝑦𝑦 ) − 𝑟 ) ( 𝑞−𝑧 ) 2

(1)

This uncertainty can be resolved either “mathematically” by using singular value decomposition to pick a minimum norm solution out of an infinite pool of solutions or “physically” by the addition of the seventh equation to zero out the excessive charges based on molecular (and multipole) symmetry. We have chosen the second approach and created the following set of rules: argmin|𝑄𝑖𝑖 | ,

|𝒅| = 0

𝑖=𝑥,𝑦,𝑧

𝑞𝑘 = 0,

𝑘=

argmin|𝑄𝑖𝑖 | , 𝑑𝛼 = 0, 𝑑𝛽 = 0, 𝑑𝛾 ≠ 0, 𝛼, 𝛽, 𝛾 ∈ {𝑥, 𝑦, 𝑧} 𝑖=𝛼,𝛽

𝛼, {0,

(2)

𝑑𝛼 = 0, 𝑑𝛽 ≠ 0, 𝑑𝛾 ≠ 0, 𝛼, 𝛽, 𝛾 ∈ {𝑥, 𝑦, 𝑧} 𝑑𝛼 ≠ 0 , ∀𝛼 ∈ {𝑥, 𝑦, 𝑧}


Page 7 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Since the orientation of the quadrupole is always fixed to the coordinate frame the dipole vector decides which charge is excessive: in case the atom has no dipole the charge on the axis with the smallest quadrupole component is set to zero; in case the dipole has two zero components, the charge on the axis with the smallest quadrupole component and a zero dipole component is set to zero; in case the dipole has only one zero component, the charge on the corresponding axis is set to zero; and finally, in case the multipoles have a dipole with nonzero components along all the axes, the central charge is zeroed out. As the last step, the extra-points with zero charges are dropped out of the model. Thus, the conversion from atomic multipoles to the MCC is performed by solving a linear system of equations. The procedure described allowed us to automatically pick the minimal charge model needed for any kind of atom. It was implemented with the help of Electrostatic Tools package.37 We have arbitrarily chosen a value of 0.1 Å for the MCC radius. It is small enough to provide high precision of multipole reproduction and large enough to avoid numerical problems during simulation, however, we believe that the choice of this radius should not matter much. The resulting charge values are available in the Supporting Information.

2.3 PES calculations All potential energy surface (PES) calculations were performed with frozen geometries of isolated molecules (as generated by OpenBabel). Ab initio PES scans were conducted with ORCA version 4.0.1 software package.40 Two levels of theory have been used. The first one is the RI-MP2/aug-ccpVTZ level of theory (aug-cc-pVTZ-PP for Br and I atoms) with a basis set superposition error (BSSE) corrected by a counterpoise procedure.41 A similar approach (with SCS-MP2) has been shown to give good estimations of halogen bonding energies “in a vacuum” at a relatively low computational cost.8 The second one, RI-MP2/6-31G* with no BSSE corrections, has been used in the GAFF validation procedure.32 It implicitly captures a certain portion of the polarization of molecules in a water medium. Van der Waals intermolecular energies were calculated with ‘sander’ program from 7 ACS Paragon Plus Environment


Page 8 of 37

Amber 11 suite. Electrostatic interaction energies were calculated with in-house software based on the Electrostatic Tools version 0.4.037 package which is available on our website.

2.4 Free energy calculation and error estimation procedures Molecular dynamics simulations were carried out in GROMACS software package, although since GROMACS does not support GAFF force field internally the topology of the bromobenzene (RESP electrostatics) solvated by 10 Å layer of TIP3P water was prepared in AmberTools 1.5 and then converted with amb2gmx.pl script. The resulting files were modified to include 4 additional extrapoints (maximum number needed for these simulations). The system was heated with Langevin thermostat up to 300K over 50 ps, then the Parrinello-Rahman barostat was engaged for another 50 ps and, after that, the system was pre-equilibrated for another 50 ps. The resulted “equilibrated” geometry and topology files were used as a starting point for all free energy runs. The thermodynamic integration was performed using Gaussian integration formula with 7 points as recommended in AMBER manual.42 A 3.5 ns simulation was carried out for each point, but only the last 3 ns were used for the analysis leaving the initial 500 ps for system equilibration. The molecular dynamics simulation was performed with 1 fs time steps with aromatic hydrogens of the benzene molecule unfrozen to prevent any issues with one of the six hydrogens being different from the others when PhH→PhF system was considered. PME was used for non-bonding van der Waals and electrostatic interactions with 10 Å cutoffs for condensed phase simulations. To minimize the integration errors the integration was performed between the closest systems, which assumes the least perturbation of the system, for example

to

obtain

ΔΔG(PhF(RESP)→PhCl(RESPQ))

two

calculations

were

performed

(PhF(RESP)→PhCl(RESP) and PhCl(RESP)→PhCl(RESPQ)). Since no atoms were dissolved completely both the separation of van der Waals and electrostatic mutations and softcore potentials were not necessary and thus were not used. Average ∂U/∂λ values and error estimates (through block


Page 9 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


averaging technique) were extracted with the gmx analyze program. The reference solvation free energy values were calculated from the data published in Abraham’s excellent paper.43 The thermodynamic integration protocol was designed to minimize the perturbation of the system for each step. In every experiment, we change the molecule a little gradually going from H to F, then from F to Cl, etc. Since the changes are small, we believe 7 middle points are enough for sufficient sampling. Moreover, such approach seems to be more consistent than the one, where a molecule is completely dissolved, or mutated from PhH to PhI directly with a fixed number of points.

3

Results and Discussion

3.1 The force field As noted above we have shown previously that fixed symmetric atomic quadrupole is the only addition to classical isotropic atomic charges that is needed to describe the anisotropic features of aromatic halogens. The purpose of this paper is to demonstrate its applicability to free energy estimations and thus any kind of MD simulations needed for the tasks of medicinal chemistry. To achieve this goal, we compare the proposed approach with two similar electrostatic models. GAFF force field was chosen as the base for the comparison not only because it is a general-purpose force field for small molecules, which is widely used in conjunction with AMBER family force fields to model ligand-protein interactions. GAFF was also built based on the same ab initio electrostatic approximations as our enhanced electrostatic models and thus quadrupolar correction proposed has both a sound theoretical basis and practical compatibility. It is crucial that the only difference in the force field models compared is the electrostatic part. All other force field parameters were taken from GAFF as is without any adjustments in order to estimate the influence of the anisotropic electrostatics on the halogen properties.



Page 10 of 37

3.2 Electrostatic models Three electrostatic models were compared: the first one is the isotropic RESP charges, the second approach, mdzq, is a fully optimized anisotropic multipole model (mdzq stands for monopole, zdirected dipole and quadrupole), and the third one is simplistic and anisotropic RESPQ model which adds fixed symmetric atomic quadrupoles for halogen atoms to the RESP atomic charges. The performance (as the ability to reproduce ab initio electrostatic potential) of these models on the halobenzene set of molecules is summarized in Table 1. While the RESP charges demonstrate good performance for light molecules, the quality of MEP reproduction deteriorates slowly with the increase of the atomic number of a halogen atom due to the increasing magnitude of the σ-hole effect.3 The mdzq model being the most flexible one demonstrates the best and uniform MEP reproduction for all the molecules of the series with benzene and iodobenzene having the largest errors. The RESPQ model demonstrates the performance with the quality lying in between the two extremes described above with errors gradually increasing from Cl to I but still being significantly smaller than those of a classical approach. Keeping in mind that this is a compromise that does not require any additional fitting compared to RESP charges these results look promising. Table 1. Root mean square errors (RMSE, kcal·mol-1·e-1) of the molecular electrostatic potential of halobenzenes for the models studied. RESP

mdzq

RESPQ

PhH

0.829

0.829

0.829

PhF

0.680

0.662

0.680

PhCl

1.097

0.602

0.965

PhBr

1.276

0.618

1.058

PhI

2.084

0.793

1.574

Additional dof*

0

2

0

*dof stands for “degrees of freedom” 10 ACS Paragon Plus Environment

Page 11 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2 shows the reference MEP and its reproduction by the models discussed in more details with bromobenzene molecule used as an example. The reference potential shows a prominent σ-hole of the bromine atom located on the continuation of the C-Br bond axis. The negative ring of the bromine atom is located roughly at the 90° angle to the C-Br bond and merges with the negative area of the aromatic π-system of the benzene ring. RESP charges fail to reproduce this anisotropy of the bromine atom predicting an excessively negative potential for the σ-hole area and an excessively positive potential for the negative ring. The most complex mdzq model describes the potential almost perfectly. It slightly overestimates the potential in the σ-hole area but the negative ring looks good. The RESP model in terms of RMSE shows the performance in between the RESP and md zq approaches however in terms of a spatial arrangement of the errors the same conclusion cannot apply. It has even smaller errors in the σ-hole area than the mdzq model, however, this comes at the cost of two positive spots on the “side” view graph located roughly at 135° to the C-Br bond. It also predicts a too negative potential around the C-Hal bond and near the α-H atoms which can be attributed to the charges of the α-hydrogens being too low in RESP (which were transferred to RESPQ model) to help to describe the negative ring of the bromine atom.



Page 12 of 37

Figure 2. The reference MEP and MEP reproduction errors for the models discussed for the bromobenzene molecule. 12 ACS Paragon Plus Environment

Page 13 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


It should be noted that further investigation of the iodobenzene molecule shows that iodine atom cannot be described with just atomic dipoles and quadrupoles (thus the RMSE rises from 0.6 to 0.8 kcal·mol-1·e-1 in mdzq model; Table 1): the MEP errors around it show a complex octahedral pattern unachievable by quadrupoles. The MEP maps for iodobenzene and other molecules are available in the Supporting Information.

3.3 Potential Energy Surface Given the differences in the MEP reproduction observed, one should expect a significant change in interaction energies. Indeed, the addition of quadrupole modifies the model potential energy surface (PES) of halobenzene – water system significantly (Figure 3). We show two ab initio references: the original one used for GAFF validation (MP2/6-31G*) and the more elaborate one which was proven to give relatively good estimates of XB energies (MP2/aug-cc-pVTZCP). The high-level ab initio calculations show attraction in two cases: when water approaches halogen with its oxygen atom at the C-X···O angle of 180° (Figure 3, XB) and when water approaches halogen with its hydrogen atom at the angle of 90° (Figure 3, HB). The low-level ab initio reference actually shows the similar interaction pattern but with expected overestimation of binding due to a small basis set and thus a large BSSE. Also, one can notice the polar flattening effect44 of the halogen atom: the repulsion for 90° and 180° approaches to bromine atom happens at significantly different separations (which is especially evident on the Br···O contacts). The default model of RESP charges fails to describe this pattern showing attraction every time the halogen is approached with hydrogen and repulsion when water touches halogen with its oxygen atom. The advanced electrostatic models give a qualitatively correct description: both XB and HB approaches have minima and other approaches show repulsion, however quantitatively they show significant variations.

Atomic quadrupole adds the expected heavy halogen character to the

interaction: on one hand, hydrogen contact becomes more favorable in the “negative ring” region in 13 ACS Paragon Plus Environment


Page 14 of 37

MEP and almost unfavorable for the σ-hole area, on the other hand, oxygen approach becomes almost attractive for the σ-hole area and more repulsive otherwise.

Figure 3. Interaction energies in bromobenzene – water system: a σ-hole approach (left column) vs a “negative ring” approach (right column) for TIP3P water molecule facing hydrogen atom (top row) vs oxygen atom (bottom row). The difference between the reference RESP model and the anisotropic (multipole) ones is great, however, the multipole approaches behave very similarly. In the areas of largest deviations, the curves for the anisotropic models lie within 0.2 kcal/mol of each other, while the error of MEP reproduction differs significantly. This can be explained by the overfitting of mdzq model: all the important features


Page 15 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of the MEP (Figure 2) can be described by the simple RESPQ model and no additional precision is needed for the PES description at the molecular mechanics level. It is worth mentioning that the anisotropic models proposed are aimed at consistently enhancing the halogen atom model, meaning it should explain both HB and XB well. Perhaps surprisingly, most of the difference between isotropic (RESP) and anisotropic approaches (RESPQ and mdzq) are attributed to HB, which is in good accordance with the fact that water is a good HB donor and a relatively poor XB acceptor. However, for a different system (with more electron-withdrawing substituents in halobenzene molecule, or with a better XB-acceptor as a solvent), XB may play a more important role in intermolecular interactions, and quadrupole correction should still be able to describe the PES correctly. As the result, we see that the RESPQ provides an adequate balance between the accuracy of PES reproduction and the complexity of the model.

3.4 Free energy reproduction Given good improvements to PES reproduction, one can expect significant changes in free energy estimates. One of the advantages of halobenzene – water system is that these compounds are studied extensively, and the experimental data is easily available. Thermodynamic integration approach was employed for solvation free energy difference calculation in the series of halobenzenes and the results of the molecular dynamics simulation were compared with the experimental data. For this experiment, advanced electrostatics was used only for Cl, Br, and I atoms since H and F are generally not believed to possess any significant anisotropic properties, which can also be seen from MEP reproduction errors of RESP model (Table 1, Figure 2). Thermodynamic integration method is mathematically well defined and physically precise, thus it should give an ideal result for an ideal force field and ideal calculation procedures. However, the real force fields contain systematic errors (thus the point of this study) and computational methods contain 15 ACS Paragon Plus Environment


Page 16 of 37

two major sources of errors: the numerical integration of ∂U/∂λ curve and ensemble averaging for each λ point (for the integration procedure). Integration errors The integration error is heavily dependent on the magnitude of the derivatives of the ∂U/∂λ function. Large derivatives cause the function to fluctuate severely, thus, requiring more sampling points for accurate integration. The ∂U/∂λ curves are smooth and even the edge points (λ=0.02544 and λ=0.97455) show no sign of singular peculiarities (Figures 4 and 5). Thus, we can state that the errors of integration should be negligible and the 7-point gaussian integration formula provides the necessary accuracy.

Figure 4. The ∂U/∂λ plot for PhBr (RESP charges) to PhI (RESP charges) transformation. Another important observation is the linearity of the ∂U/∂λ derivative for quadrupole addition transformation (Figure 5). This means that quadrupole does not bring any features to the integrand function, so it has no effect on the precision of the calculation and thus models employing atomic quadrupoles can be used in any of well-established simulation protocols without modification.


Page 17 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. The ∂U/∂λ plot for quadrupole introduction transformation for PhBr molecule. Averaging errors The second source of errors, averaging of ∂U/∂λ values over the trajectory, was addressed in a standard way with the block-averaging method. However, this procedure sometimes can overestimate or, even worse, underestimate the errors due to the peculiarities of the dataset, so to doublecheck the results additional calculations were performed to obtain the excessive ΔΔG values and close thermodynamic cycles. As can be seen in Figure 6 all thermodynamic cycles are closed with free energy change lying within the error margins. For example, the PhBr(RESP)→PhI(RESP) energy difference

is

0.2±0.03

kcal/mol

which

is

in

accordance

with

PhBr(RESP)→PhBr(RESPQ)→PhI(RESPQ)→PhI(RESP) (-0.38 - 0.34 + 0.94 = 0.22±0.03 kcal/mol). Thus, we can state that our calculations and error estimates are accurate enough to reliably compare the electrostatic models’ performance within the protocol described.



Page 18 of 37

Figure 6. Thermodynamic integration paths and solvation free energy difference results (in kcal/mol; doubled standard error is shown as error estimates). The rows correspond to different molecules while the columns correspond to different models tested. The arrows showing the transformation paths and the corresponding values are colored for the ease of reading. Experimental values are calculated from Abraham’s data.43 Free energy results Finally, once the reliability of the simulation procedure has been established, the actual results can be analyzed (Figure 7). The performance of the new models looks very promising (note, that H and F atom models were not changed). While the classical atomic charges approach struggles to describe the F < Cl < Br < I trend, the anisotropic models give a very accurate prediction of solvation differences in this series.


Page 19 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7. The results of solvation free energy differences calculations for the models studied compared to the experimental43 estimates of solvation free energies of benzene and halobenzenes. All the calculated energies were obtained as the sum of the experimental free energy of solvation of phenylfluoride and the corresponding simulated ΔΔG values. The error bars represent thermodynamic integration error estimates (doubled standard error). It is worth mentioning, that in our study we modified only molecules containing Cl, Br, and I atoms as donors of halogen bonds. The electrostatics of hydrogen and fluorine is believed to be described well and we see no evidence of the opposite. However, GAFF – as our baseline for the classical force field – fails to correctly assess the difference between hydrogen and fluorine, so even with the almost perfect description in the F, Cl, Br, I series RESPQ fails to correctly estimate the solvation free energy difference between unsubstituted and halogen-substituted aromatics. We believe that the issue might be in GAFF nonbonding parameters. The investigation of van der Waals potential of X···O of the force field (Figure 8) shows that a fluorine atom, on one hand, is considered to be larger than an aromatic hydrogen and has a significantly deeper energy well. On the other hand, in classical QSPR models giving excellent estimations of solvation energies45, 46 the fluorine atom is often treated as 19 ACS Paragon Plus Environment


Page 20 of 37

hydrogen. Further refinement of GAFF parameters for aromatic H and F atoms is needed and, for now, additional caution should be taken for investigation of systems with H - Hal divergencies.

Figure 8. GAFF van der Waals potentials for X···O pairs (X=H,F,Cl,Br,I). Despite the mdzq model being optimized to precisely reflect MEP of the corresponding molecules, it shows the same performance as RESPQ model which differs from the original RESP model only by the quadrupoles of the fixed values. Moreover, RESPQ gives a smoother transition in the heavy halogen series and demonstrates slightly smaller errors compared to the experiment. Combining its good performance with the fact that RESPQ does not contain any additional degrees of freedom compared to RESP charges makes RESPQ the model of choice for modeling of heavy halogens.

4

Conclusion A new quadrupole-based approach for halogen bonding description was tested in halobenzene-water

interaction simulations. It was shown, that the suggested RESPQ electrostatics model built by the addition of fixed atomic quadrupoles to heavy halogen atoms, provides a consistently good description of intermolecular interactions between halogen atom and both hydrogen bond donors and halogen


Page 21 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


bond acceptors. The solvation free energy differences for F-, Cl-, Br-, and I-substituted benzenes are in excellent agreement with the experimental values. The quadrupoles used in RESPQ values were obtained by fitting electrostatic models to an ab initio reference MEP and thus are independent of the force field. So, a similar approach should work for any force field, although only GAFF was tested by us. This approach of quadrupole addition should work for any kind of numerical experiment in the field of medicinal chemistry ranging from a direct molecular dynamics investigation of a protein-ligand system to an incorporation into scoring functions and QSAR/QSPR models.

ASSOCIATED CONTENT Supporting Information. The final atomic charges used, the description of the multipole optimization procedure and the conversion from atomic multipoles to the multipole charge cluster procedure, the MEP maps for the reference potential, its prediction, and prediction errors for all molecules (PDF). This information is available free of charge via the internet at http://pubs.acs.org. AUTHOR INFORMATION Corresponding Author *Corresponding

author

phone:

+7

(495)

939-3969,

fax:

+7

(495)

939-0290

e-mail: [email protected] Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. 21 ACS Paragon Plus Environment


Page 22 of 37

Funding Sources This work has been funded by the Russian Foundation for Basic Research (RFBR grant 18-0301065a). ABBREVIATIONS MCC, multipole charge cluster; MEP, molecular electrostatic potential; PES, potential energy surface; HB, hydrogen bonding; XB, halogen bonding; RMSE, root mean square error REFERENCES 1.

Paulini, R.; Müller, K.; Diederich, F., Orthogonal Multipolar Interactions in Structural

Chemistry and Biology. Angew. Chem., Int. Ed. 2005, 44 (12), 1788-1805. 2.

Clark, T.; Hennemann, M.; Murray, J. S.; Politzer, P., Halogen Bonding: the σ-Hole. J. Mol.

Model. 2007, 13 (2), 291-296. 3.

Politzer, P.; Lane, P.; Concha, M. C.; Ma, Y.; Murray, J. S., An Overview of Halogen

Bonding. J. Mol. Model. 2007, 13 (2), 305-311. 4.

Bissantz, C.; Kuhn, B.; Stahl, M., A Medicinal Chemist’s Guide to Molecular Interactions. J.

Med. Chem. 2010, 53 (14), 5061-5084. 5.

Clark, T., σ-Holes. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 3 (1), 13-20.

6.

Desiraju, G. R.; Ho, P. S.; Kloo, L.; Legon, A. C.; Marquardt, R.; Metrangolo, P.; Politzer,

P.; Resnati, G.; Rissanen, K., Definition of the Halogen Bond (IUPAC Recommendations 2013). Pure and Appl. Chem. 2013, 85 (8), 1711-1713. 7.

Politzer, P.; Murray, J. S.; Clark, T., Halogen Bonding and Other σ-Hole Interactions: A

Perspective. Phys. Chem. Chem. Phys. 2013, 15 (27), 11178-89. 22 ACS Paragon Plus Environment

Page 23 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


8.

Wolters, L. P.; Schyman, P.; Pavan, M. J.; Jorgensen, W. L.; Bickelhaupt, F. M.; Kozuch,

S., The Many Faces of Halogen Bonding: a Review of Theoretical Models and Methods. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4 (6), 523-540. 9.

Kolář, M. H.; Hobza, P., Computer Modeling of Halogen Bonds and Other σ-Hole

Interactions. Chem. Rev. 2016, 116 (9), 5155-5187. 10. Lu, Y.; Liu, Y.; Xu, Z.; Li, H.; Liu, H.; Zhu, W., Halogen Bonding for Rational Drug Design and New Drug Discovery. Expert Opin. Drug Discovery 2012, 7 (5), 375-383. 11. Scholfield, M. R.; Zanden ,C. M. V.; Carter, M.; Ho, P. S., Halogen Bonding (X‐Bonding): A Biological Perspective. Protein Sci. 2012, 22 (2), 139-152. 12. Gilday, L. C.; Robinson, S. W.; Barendt, T. A.; Langton, M. J.; Mullaney, B. R.; Beer, P. D., Halogen Bonding in Supramolecular Chemistry. Chem. Rev. 2015, 115 (15), 7118-7195. 13. Cavallo, G.; Metrangolo, P.; Milani, R.; Pilati, T.; Priimagi, A.; Resnati, G.; Terraneo, G., The Halogen Bond. Chem. Rev. 2016, 116 (4), 2478-2601. 14. Li, B.; Zang, S.-Q.; Wang, L.-Y.; Mak, T. C. W., Halogen Bonding: A Powerful, Emerging Tool for Constructing High-Dimensional Metal-Containing Supramolecular Networks. Coord. Chem. Rev. 2016, 308, 1-21. 15. Murray, J. S.; Lane, P.; Clark, T.; Politzer, P., σ-Hole Bonding: Molecules Containing Group VI Atoms. J. Mol. Model. 2007, 13 (10), 1033-1038. 16. Politzer, P.; Murray, J. S.; Concha, M. C., σ-Hole Bonding Between Like Atoms; A Fallacy of Atomic Charges. J. Mol. Model. 2008, 14 (8), 659-665.



Page 24 of 37

17. Murray, J. S.; Lane, P.; Politzer, P., Expansion of the σ-Hole Concept. J. Mol. Model. 2009, 15 (6), 723-729. 18. Ibrahim, M. A. A., Molecular Mechanical Study of Halogen Bonding in Drug Discovery. J. Comput. Chem. 2011, 32 (12), 2564-2574. 19. Ibrahim, M. A. A., AMBER Empirical Potential Describes the Geometry and Energy of Noncovalent Halogen Interactions Better than Advanced Semiempirical Quantum Mechanical Method PM6-DH2X. J. Phys. Chem. B 2012, 116 (11), 3659-3669. 20. Rendine, S.; Pieraccini, S.; Forni, A.; Sironi, M., Halogen Bonding in Ligand-Receptor Systems in the Framework of Classical Force Fields. Phys. Chem. Chem. Phys. 2011, 13 (43), 1950819516. 21. Kolář, M.; Hobza, P., On Extension of the Current Biomolecular Empirical Force Field for the Description of Halogen Bonds. J. Chem. Theory Comput. 2012, 8 (4), 1325-1333. 22. Jorgensen, W. L.; Schyman, P., Treatment of Halogen Bonding in the OPLS-AA Force Field: Application to Potent Anti-HIV Agents. J. Chem. Theory Comput. 2012, 8 (10), 3895-3901. 23. Soteras Gutiérrez, I.; Lin, F.-Y.; Vanommeslaeghe, K.; Lemkul, J. A.; Armacost, K. A.; Brooks, C. L.; MacKerell, A. D., Parametrization of Halogen Bonds in the CHARMM General Force Field: Improved Treatment of Ligand–Protein Interactions. Bioorg. Med. Chem. 2016, 24 (20), 48124825. 24. Lin, F.-Y.; MacKerell, A. D., Polarizable Empirical Force Field for Halogen-Containing Compounds Based on the Classical Drude Oscillator. J. Chem. Theory Comput. 2018, 14 (2), 10831098.


Page 25 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


25. Mu, X.; Wang, Q.; Wang, L.-P.; Fried, S. D.; Piquemal, J.-P.; Dalby, K. N.; Ren, P., Modeling Organochlorine Compounds and the σ-Hole Effect Using a Polarizable Multipole Force Field. J. Phys. Chem. B 2014, 118 (24), 6456-6465. 26. Bereau, T.; Kramer, C.; Meuwly, M., Leveraging Symmetries of Static Atomic Multipole Electrostatics in Molecular Dynamics Simulations. J. Chem. Theory Comput. 2013, 9 (12), 54505459. 27. Bereau, T.; Andrienko, D.; von Lilienfeld, O. A., Transferable Atomic Multipole Machine Learning Models for Small Organic Molecules. J. Chem. Theory Comput. 2015, 11 (7), 3225-3233. 28. El Hage, K.; Bereau, T.; Jakobsen, S.; Meuwly, M., Impact of Quadrupolar Electrostatics on Atoms Adjacent to the Sigma-Hole in Condensed-Phase Simulations. J. Chem. Theory Comput. 2016, 12 (7), 3008-3019. 29. Hédin, F.; El Hage, K.; Meuwly, M., A Toolkit to Fit Nonbonded Parameters from and for Condensed Phase Simulations. J. Chem. Inf. Model. 2016, 56 (8), 1479-1489. 30. Titov, O. I.; Shulga, D. A.; Palyulin, V. A.; Zefirov, N. S., Description of Halogen Bonding on the Basis of Multicenter Multipole Expansion. Dokl. Chem. 2013, 450 (1), 139-143. 31. Titov, O. I.; Shulga, D. A.; Palyulin, V. A.; Zefirov, N. S., Perspectives of Halogen Bonding Description in Scoring Functions and QSAR/QSPR: Substituent Effects in Aromatic Core. Mol. Inf. 2015, 34 (6-7), 404-416. 32. Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A., Development and Testing of a General AMBER Force Field. J. Comput. Chem. 2004, 25 (9), 1157-1174.



Page 26 of 37

33. Jakalian, A.; Jack, D. B.; Bayly, C. I., Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: II. Parameterization and Validation. J. Comput. Chem. 2002, 23 (16), 1623-41. 34. O'Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R., Open Babel: An Open Chemical Toolbox. J. Cheminf. 2011, 3 (1), 33. 35. Granovsky, A. A. Firefly version 8.2. www http://classic.chem.msu.su/gran/firefly/index.html. 36. Case, D. A.; Brozell, S. R.; Cerutti, D. S.; Cheatham, D. E. I.; Cruzeiro, V. W. D.; Darden, T. A.; Duke, R. E.; Ghoreishi, D.; Gohlke, H.; Goetz, A. W.; Greene, D.; Harris, R.; Homeyer, N.; Izadi, S.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Lin, C.; Liu, J.; Luchko, T.; Luo, R.; Mermelstein, D. J.; Merz, K. M.; Miao, Y.; Monard, G.; Nguyen, H.; Omelyan, I.; Onufriev, A.; Pan, F.; Qi, R.; Roe, D. R.; Roitberg, A.; Sagui, C.; Schott-Verdugo, S.; Shen, J.; Simmerling, C. L.; Smith, J.; Swails, J.; Walker, R. C.; Wang, J.; Wei, H.; Wolf, R. M.; Wu, X.; Xiao, L.; York, D. M.; Kollman, P. A., AmberTools 1.5. University of California, San Francisco, 2011. 37. Titov, O. I.;

Shulga, D. A.; Palyulin, V. A. Electrostatic Tools v 0.4.0.

http://qsar.chem.msu.ru/elec_tools/. 38. Devereux, M.; Raghunathan, S.; Fedorov, D. G.; Meuwly, M., A Novel, Computationally Efficient Multipolar Model Employing Distributed Charges for Molecular Dynamics Simulations. J. Chem. Theory Comput. 2014, 10 (10), 4229-4241. 39. Unke, O. T.; Devereux, M.; Meuwly, M., Minimal Distributed Charges: Multipolar Quality at the Cost of Point Charge Electrostatics. J. Chem. Phys. 2017, 147 (16), 161712. 40. Neese, F., The ORCA Program System. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2 (1), 73-78. 26 ACS Paragon Plus Environment

Page 27 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


41. Boys, S. F.; Bernardi, F., The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19 (4), 553-566. 42. Case, D. A.; Darden, T. A.; Cheatham, T. E. I.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M.-J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER12 User Manual. http://ambermd.org/doc12/Amber12.pdf. 43. Abraham, M. H.; Andonian-Haftvan, J.; Whiting, G. S.; Leo, A.; Taft, R. S., Hydrogen bonding. Part 34. The Factors That Influence the Solubility of Gases and Vapours in Water at 298 K, and a New Method for its Determination. J. Chem. Soc., Perkin Trans. 2 1994, (8), 1777-1791. 44. Sedlak, R.; Kolář, M. H.; Hobza, P., Polar Flattening and the Strength of Halogen Bonding. J. Chem. Theory Comput. 2015, 11 (10), 4727-4732. 45. Antipin, I. S.; Arslanov, N. A.; Palyulin, V. A.; Konovalov, A. I.; Zefirov, N. S., Solvation Topological Index - A Topological Model for the Description of Dispersional Interactions. Dokl. Akad. Nauk SSSR 1991, 316 (4), 925-927. 46. Antipin, I. S.; Arslanov, N. A.; Palyulin, V. A.; Konovalov, A. I.; Zefirov, N. S., Prediction of Nonspecific Solvation Enthalpy for Organic Non-Electrolytes. Dokl. Akad. Nauk SSSR 1993, 331 (2), 173-176.



Page 28 of 37

TABLE OF CONTENTS GRAPHIC


Page 29 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Extra-point charges added to the bromobenzene molecule to represent atomic dipole and quadrupole (MCC radius is increased to 1.0 Å) with local multipole frame axes directions labeled. The RESP model uses no extra-points, mdzq uses ±Z and ±X, while RESPQ needs only ±Z. 84x33mm (300 x 300 DPI)

ACS Paragon Plus Environment


Figure 2. The reference MEP and MEP reproduction errors for the models discussed for the bromobenzene molecule. 177x224mm (300 x 300 DPI)


Page 30 of 37

Page 31 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Interaction energies in bromobenzene – water system: a σ-hole approach (left column) vs a “negative ring” approach (right column) for TIP3P water molecule facing hydrogen atom (top row) vs oxygen atom (bottom row). 169x135mm (300 x 300 DPI)



Figure 4. The ∂U/∂λ plot for PhBr (RESP charges) to PhI (RESP charges) transformation. 84x84mm (300 x 300 DPI)


Page 32 of 37

Page 33 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. The ∂U/∂λ plot for quadrupole introduction transformation for PhBr molecule. 84x84mm (300 x 300 DPI)



Figure 6. Thermodynamic integration paths and solvation free energy difference results (in kcal/mol; doubled standard error is shown as error estimates). The rows correspond to different molecules while the columns correspond to different models tested. The arrows showing the transformation paths and the corresponding values are colored for the ease of reading. Experimental values are calculated from Abraham’s data. 84x50mm (300 x 300 DPI)


Page 34 of 37

Page 35 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7. The results of solvation free energy differences calculations for the models studied compared to the experimental43 estimates of solvation free energies of benzene and halobenzenes. All the calculated energies were obtained as the sum of the experimental free energy of solvation of phenylfluoride and the corresponding simulated ΔΔG values. The error bars represent thermodynamic integration error estimates (doubled standard error). 84x84mm (300 x 300 DPI)



Figure 8. GAFF van der Waals potentials for X···O pairs (X=H,F,Cl,Br,I). 84x84mm (300 x 300 DPI)


Page 36 of 37

Page 37 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table of Contents and Abstract Graphic 82x44mm (300 x 300 DPI)


The quadrupole correction: from molecular electrostatic potential to

Recommend Documents