Effective Fragment Potential Method for H-Bonding ... - ACS Publications

Jun 6, 2017 - Department of Condensed Matter Physics, National Research Nuclear University “MEPhI”, 31 Kashirskoe Highway, Moscow,. 115409, Russia...
0 downloads 0 Views 2MB Size
Subscriber access provided by UNIV OF ARIZONA

Article

Effective Fragment Potential Method for H-Bonding: How to Obtain Parameters for Non-Rigid Fragments Nikita O. Dubinets, and Lyudmila V. Slipchenko J. Phys. Chem. A, Just Accepted Manuscript • Publication Date (Web): 06 Jun 2017 Downloaded from http://pubs.acs.org on June 8, 2017

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry A 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 22

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

The Journal of Physical Chemistry

Effective Fragment Potential Method for H-bonding: How to Obtain Parameters for Non-rigid Fragments Nikita Dubinetsa,b and Lyudmila V. Slipchenkoc* a

National Research Nuclear University “MEPhI”, Department of Condensed Matter Physics, 31 Kashirskoe highway, Moscow, 115409, Russia b Federal Research Center “Crystallography and Photonics”, Photochemistry Center RAS, 7a Novatorov str., Moscow, 119421, Russia c Department of Chemistry, Purdue University, 560 Oval Drive, West Lafayette, IN 47907, United States * Corresponding author’s email: [email protected]

Abstract Accuracy of the Effective Fragment Potential (EFP) method was explored for describing intermolecular interaction energies in three dimers with strong H-bonded interactions, formic acid, formamide and formamidine dimers, which are a part of HBC6 database of non-covalent interactions. Monomer geometries in these dimers change significantly as a function of intermonomer separation. Several EFP schemes were considered, in which fragment parameters were prepared either for a fragment in its gas-phase geometry or recomputed for each unique fragment geometry. Additionally, a scheme in which gas-phase fragment parameters are shifted according to relaxed fragment geometries is introduced and tested. EFP data are compared against the coupled cluster with single, double and perturbative triple excitations (CCSD(T)) method in a complete basis set (CBS) and the symmetry adapted perturbation theory (SAPT). All considered EFP schemes provide a good agreement with CCSD(T)/CBS for binding energies at equilibrium separations, with discrepancies not exceeding 2 kcal/mol. However, only the schemes that utilize relaxed fragment geometries remain qualitatively correct at shorter than equilibrium intermolecular distances. The EFP scheme with shifted parameters behaves quantitatively similar to the scheme in which parameters are recomputed for each monomer geometry and thus is recommended as a computationally efficient approach for large-scale EFP simulations of flexible systems.

Introduction Intermolecular interactions play a leading role in determining structures and properties of biological macromolecules such as DNA and proteins, organic crystals, liquids and colloids. The challenge in theoretical modeling of non-covalent interactions arises due the large number of small-magnitude interactions, approximately scaling as N2 with the total number of atoms in the system. Additionally, non-covalent forces have a complicated make-up ranging from Coulombbased electrostatic components to van der Waals terms originating in long-range electron correlation and Pauli exclusion principle. Many-body effects also contribute to complexity associated with modeling of non-covalent interactions.

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry

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

Page 2 of 22

The Effective Fragment Potential (EFP) method is one of promising and computationally efficient ways to describe non-covalent interactions.1-7 EFP is a model potential based on first principles that represents the total intermolecular energy of the system as a sum of electrostatic (Coulomb), polarization (induction), dispersion, exchange-repulsion, and charge-transfer interactions. All EFP terms are based on a distributed approach, in which parameters are located at multiple points within a fragment. The EFP method has been previously applied to modeling structures and binding patterns in clusters and liquids8-12, including water-alcohol13-15 and waterbenzene systems16, DNA base pairs5, 17, to name a few. Additionally, EFP potentials can be successfully used for describing solute-solvent interactions in polarizable QM/MM-type schemes, as was shown in a number of recent work.18-25 Recently, the EFP methodology has been extended to biological macromolecules such as peptide sequences, by developing a scheme of splitting the peptide into smaller fragments and preparing parameters for individual fragments.26 However, one of the present limitations of the EFP method is that the formulism is based on rigid fragments. In the original EFP method, parameters for a rigid-geometry fragment are precomputed and used in all subsequent calculations of the system of interest. This approach works well for weakly interacting fragments, which is a common case in typical H-bonded and van der Waals complexes. However, rigidity of fragments becomes a drawback when stronger interactions are considered, such as in complexes with ionic or multiple hydrogen bonds or in molecules with flexible degrees of freedom, such as biological or organic polymers. Even working with small systems such as dimers, we have found several occasions when the gasphase structures differ significantly from the monomer structures adapted in the interacting complex, with a few examples being strongly H-bonded dimers such as formic acid dimer or adenine-thymine dimer.9 The issue of flexible fragments is elegantly addressed in a hybrid of the EFP and FMO (Fragment Molecular Orbital)27-30 methods, called the Effective Fragment Molecular Orbital (EFMO) method.6, 31-33 In EFMO, each fragment is described quantum-mechanically within the FMO framework, but non-covalent interactions between the fragments are described at the EFP level. The resulting approach is computationally more efficient and often more accurate than the original FMO method, in which both inter- and intra-fragment interactions are described quantum-mechanically. On the other hand, EFMO improves on the main shortcoming of the EFP by allowing intra-monomer flexibility, but for a price of significantly increasing the computational cost of EFP calculations. Specifically, recomputing EFP parameters for each unique fragment geometry is a significant computational expense of EFMO calculations. In this work we explore an alternative way of introducing intra-monomer flexibility in the EFP method, without a need of recomputing parameters of fragments. Can the EFP method be used “as is” for fragments with non-rigid geometries? How to consistently obtain EFP parameters for fragments that change their geometry along a considered chemical process? As a first step toward this goal, we measure accuracy of various EFP schemes on example of one of the databases developed to benchmark non-covalent interactions, namely the HBC6 database.34 The HBC6 database contains six doubly-H-bonded dimers based on formic acid (FaOO), formamide (FaON), and formamidine (FaNN) monomers computed at varying inter-monomer separations,

ACS Paragon Plus Environment

2

Page 3 of 22

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

The Journal of Physical Chemistry

and presents a challenge for a conventional EFP method in which geometries of each monomer are kept rigid. In the present work we validate the approach in which the EFP parameters are shifted based on geometric changes of each monomer, using homo-dimers from the HBC6 database as a test set.

Computational details Three H-bonded dimers from the HBC6 database were examined, formic acid (FaOO) dimer, formamide (FaON) dimer, and formamidine (FaNN) dimer. For each dimer, a set of fully optimized structures at the CCSD(T)/aug-cc-pVDZ level exists for inter-dimer separations between 3.4 Å and 10 Å. Interaction energies with CCSD(T) in the complete basis set (CBS) and with the symmetry adapted perturbation theory (SAPT) are available for the dataset.34 SAPT2+3(CCD)35/aug-cc-pVTZ results are used as reference values in this work. We computed a set of EFP potential energy slices with 3.4 Å – 10.0 Å intermonomer separations for each dimer. At distances shorter than 4.4 Å, the step was 0.1 Å; for distances between 4.6 Å and 7.0 Å, the step was 0.2 Å – 0.6 Å. Finally, for separations longer than 7.0 Å, the step was 1.0 Å. EFP parameters are computed in a single electronic structure calculation (so called MAKEFP job) in the GAMESS quantum chemistry package36-37. The parameters include (1) electrostatic multipoles (charges, dipoles, quadrupoles, octopoles) distributed at atoms and bond mid-points, (2) static and dynamic polarizability tensors computed at centroids of localized molecular orbitals (LMO), (3) localized wave function, Fock matrix, atomic basis set. Effectively, two sets of coordinate points are used in EFP calculations, coordinates of atoms and related coordinates of bond mid-points, and coordinates of LMO centroids. Electrostatic multipoles and atomic basis set are stored at the former points, while polarizabilities are kept at the latter ones. Additionally, coordinates of the LMO centroids are used in expressions for exchange-repulsion and chargetransfer terms and overlap-based damping functions for electrostatics and dispersion.38 Even though most of the EFP parameter components are tensor properties that are directionally dependent (with exception of charges and dynamic polarizabilities that are scalars), we do not exploit directional dependence of the parameters in this work. In other words, tensors (dipoles, quadrupoles, octopoles, polarizabilities) are translated to new positions but not rotated. This approach can be considered as a simplified version of the procedure used in Amoeba and other force fields employing dipole and quadrupole moments in which tensors are also rotated following a rotation of a local coordinate system associated with neighboring atoms.39 EFP energy terms have different dependence on the basis set used for computing corresponding parameters. For example, working expressions for the EFP exchange-repulsion energy are obtained in the assumption of a complete basis set.40-41 Static and dynamic polarizabilities are known to converge slowly with the basis set, such that a large diffuse basis is needed for convergence of polarization and dispersion terms. Electrostatic multipoles (obtained from a distributed multipole analysis, DMA) become unstable in a diffuse basis set and are better computed in a compact basis set.42 With these considerations, we typically recommend to generate EFP parameters in a mixed basis, using 6-31+G*43-44 or 6-31G* for electrostatic multipoles and 6-311++G(3df,2p)45-46 for other terms (exchange-repulsion, polarization,

ACS Paragon Plus Environment

3

The Journal of Physical Chemistry

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

Page 4 of 22

dispersion, charge-transfer).9 In the present work, we use the mixed 6-31+G*/6-311++G(3df,2p) basis in one of the employed EFP schemes and the small 6-31+G* basis in three other schemes in order to simplify both computations and analysis. In order to understand what is the most optimal way to compute EFP interaction energies and energy components for fragments with non-rigid geometries, we employed four schemes for generating EFP parameters. These are: - Scheme 1, “optimized parameters”: A complete set of EFP parameters is generated for each monomer’s structure in the 6-31+G* basis set. This scheme is expected to be the most rigorous, but a computational cost of recomputing parameters at each unique geometry of a fragment might be impractical for large-scale computations. - Scheme 2, “gas phase parameters”: The monomer structures extracted from the dimers separated by 10.0 Å are considered as the “gas phase” structures of the monomers. EFP parameters computed for these “gas phase” monomers in 6-31+G* basis are used in all calculations, i.e., EFP fragments are kept unchanged. This is a standard way of performing EFP calculations. - Scheme 3, “mixed-basis gas phase parameters”: The scheme is similar to Scheme 2, but the mixed basis set is used: 6-31+G* for electrostatic multipoles (for electrostatic term) and the 6-311++G(3df,2p) basis for all other terms. This scheme was used extensively in the past.9, 13 - Scheme 4, “shifted parameters”: This scheme is a hybrid of Schemes 1 and 2. First, “gas phase” parameters of the fragments are computed as in Scheme 2. However, these parameters are shifted to geometric points corresponding to fully optimized fragments, as in Scheme 1. Practically, it means that in the “gas phase” parameter file, coordinates of atoms, bond midpoints and LMO centroids are adjusted to match the coordinates of the corresponding fragment from Scheme 1. Scheme 4 is tried for the first time in this work. If successful, Scheme 4 might bring a computationally-efficient alternative to Scheme 1, in case non-rigid fragments are to be modeled. In all considered schemes, short-range damping parameters for electrostatic, polarization, and dispersion terms were employed. Overlap-based dampings for electrostatics and dispersion, and Gaussian-like damping for polarization with default damping parameters were used, as described in Ref. 38. When the schemes with non-gas-phase fragment geometries are used (i.e., Schemes 1 and 4), one also needs to account for a deformation energy contribution to the total interaction energy, i.e.: ௜௡௧ ஺ ஻ ஺ ஻ = ‫ ܧ‬௜௡௧ + ൫‫ܧ‬௥௟௫ − ‫ܧ‬௚௔௦ ൯ + (‫ܧ‬௥௟௫ − ‫ܧ‬௚௔௦ ), (1) ‫ܧ‬௥௟௫ ஺/஻ ஺/஻ where ‫ܧ‬௥௟௫ − ‫ܧ‬௚௔௦ is the deformation energy required to bring monomer A or B from its equilibrium geometry to that in the complex. ‫ ܧ‬௜௡௧ is the interaction energy for the dimer ௜௡௧ computed at its relaxed (optimized) geometry, and ‫ܧ‬௥௟௫ is the interaction energy of the dimer with respect to the infinitely separated monomers in their equilibrium geometries. Note that if the gas phase geometries of the fragments are used, the relaxed interaction energy is equal to the interaction energy as directly obtained from EFP calculations. In this work, deformation energies were computed at the HF/6-31+G* level of the theory, and the data were extracted from MAKEFP parameter calculations. For reference, these deformation energies are provided in Table S1 of Supporting Materials.

ACS Paragon Plus Environment

4

Page 5 of 22

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

The Journal of Physical Chemistry

Two versions of SAPT calculations are reported. Both are performed at the SAPT2+3(CCD)/aug-cc-pVTZ level. One version is referred to as “SAPT: optimized geometries” and uses optimized geometries of the dimers at each separation. These data are based on the original values provided in Ref. 34 but are modified by adding the deformation energy term (Eq. 1) computed at the CCSD(T)/CBS level, which corresponds to the following ௜௡௧ ௜௡௧ energy difference: ∆‫ܧ‬௥௟௫ = ‫ܧ‬௥௟௫ − ‫ܧ‬௖௣ (from Eq. 4 in Ref. 34). Another set of SAPT values, called “SAPT: gas phase geometries” is based on the gas-phase geometries of the monomers, i.e., geometries that the monomers possess in a dimer separated by 10.0 Å. These data were computed in the present work. SAPT energies computed using the optimized geometries of the monomers can be used as a reference for the EFP values computed using the optimized geometries of the monomers, i.e., Schemes 1 and 4. Similarly, the “gas-phase” SAPT correlates with the EFP employing gas-phase parameters, i.e., Schemes 2 and 3. Note that the energy components in both SAPT and EFP methods are uncorrected for relaxation energies, such that the energy components for the “optimized geometry” sets do not sum up to the total energies.

Results and discussion Total interaction energies Total EFP interaction energies and energy components for three dimers are shown in Figs. 1-3. CCSD(T)/CBS data from Ref. 34 are provided for reference. SAPT2+3(CCD)/aug-cc-pVTZ energies, computed as described in Computational Section, are also shown. Comparison of equilibrium distances, total binding energies and interaction energies at a shorter 3.5 Å separation for all EFP schemes and CCSD(T) and SAPT data are reported in Table 1. Additionally, Table 1 contains mean absolute errors (MAE) of the total interaction energies of the EFP schemes with respect to CCSD(T)/CBS, computed over the whole range of considered separations in the dimers.

ACS Paragon Plus Environment

5

The Journal of Physical Chemistry

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

Page 6 of 22

3.4Å – 10.0Å

FaOO - FaOO

Figure 1. Total interaction energies and energy components (kcal/mol) in FaOO dimer.

As follows from Figs. 1-3 and Table 1, all EFP schemes consistently overestimate equilibrium distances by 0.0-0.2 Å with respect to CCSD(T)/CBS. We observed slight overestimation of equilibrium separations by EFP in our previous work as well. All EFP schemes produce equilibrium separations within 0.1 Å from each other.

ACS Paragon Plus Environment

6

Page 7 of 22

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

The Journal of Physical Chemistry

3.4Å – 10.0Å

FaON - FaON

Figure 2. Total interaction energies and energy components (kcal/mol) in FaON dimer.

We can derive the following observations from analysis of the total binding energies. Scheme 1, in which parameters are recomputed for each monomer geometry, is the most binding from the EFP curves. Scheme 2 with the gas-phase parameters is always less binding than Scheme 1, by 0.1 - 1.4 kcal/mol. Further, Scheme 3 with the mixed-basis gas-phase parameters is less binding than Scheme 2, by 1.6 - 1.8 kcal/mol. We can expect that if the parameters for Schemes 1 and 4

ACS Paragon Plus Environment

7

The Journal of Physical Chemistry

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

Page 8 of 22

were prepared using a mixed basis, these Schemes would produce similarly less attractive (i.e., with interaction energies smaller by a magnitude) curves. We have shown previously that the mixed-basis scheme with the gas-phase parameters, i.e., Scheme 3, is reliable in description of various systems at or near equilibrium separations.9 Behavior of Scheme 4, which is a hybrid scheme in which the gas phase parameters are shifted to relaxed geometries of the fragments, is a topic of a special interest. For all considered dimers Scheme 4 is less binding than Scheme 1, by 1.2 - 1.3 kcal/mol. With respect to gas-phase Scheme 2, Scheme 4 can be either slightly overbinding, as in FaNN dimer, or under-binding, as in FaON and FaOO dimers. With respect to CCSD(T)/CBS reference values, all considered EFP schemes show a good agreement for the values of binding energies. Specifically, Scheme 1 (optimized parameters) overbinds the reference by 0.0 – 2.0 kcal/mol. Scheme 2 (gas phase parameters) is within +0.4 (-0.9) kcal/mol, i.e., can be both underbinding and overbinding, with respect to CCSD(T)/CBS. Scheme 3 (mixed-basis gas-phase parameters) always underbinds by 0.8-2.1 kcal/mol. Finally, Scheme 4 with shifted parameters varies from 1.2 kcal/mol underbinding in FaOO dimer to 0.8 kcal/mol overbinding in FaNN dimer. Errors in EFP binding energies are comparable to the corresponding errors of SAPT with respect to CCSD(T), which constitute 0.7-1.5 kcal/mol for the considered dimers. These errors are typical for this level of SAPT (SAPT2+3(CCD)/aug-ccpVTZ) in describing H-bonded complexes.35 Slightly different variants of SAPT such as SAPT2+(3)δMP2 or SAPT2+(CCD)δMP2 in aug-cc-pVTZ basis provide a better description of H-bonded and other non-covalent interactions.35 To summarize, if only interaction energies near equilibrium distances are in question, gas-phase Scheme 2 produces the closest agreement with the CCSD(T)/CBS values. However, situation is quite different if interaction energies at shorter separations are taken into account. In Table 1 we show two metrics that shed light on short-range behavior of different EFP schemes, namely, interaction energies at 3.5 Å separation, which is 0.3 – 0.7 Å shorter than the equilibrium distances in the considered dimers, and MAE of the EFP interaction energies computed over the whole distance range, i.e., between 3.4 and 10.0 Å. Both metrics reveal deficiencies of gasphase-based Schemes 2 and 3. It should be noted also that prediction of interaction energies at short intermolecular separations is a challenging test for methods based on the perturbation theory, and significantly higher discrepancies are expected, as will be obvious from analysis of both EFP and SAPT data. As follows from Table 1, all EFP schemes underestimate interaction energies at short separations. One reason, common to all EFP schemes, is an overestimation of the equilibrium distances, which results in shifting the EFP curves to the right with respect to the reference curves and in increasing discrepancies in energies in the region with repulsive potentials. Accordingly, in the dimers where EFP overestimates the equilibrium distances the most (FaOO dimer, overestimation by 0.2 Å), discrepancies in the interaction energies at 3.5 Å are the biggest, 16-17 kcal/mol for Schemes 1 and 4. On the other hand, in FaNN dimer, where the EFP equilibrium distance is within 0.0-0.1 Å from the reference, the short-range interaction energies by Schemes 1 and 4 are within 3 kcal/mol from the CCSD(T) values. Another reason of significant discrepancies in the interaction energies, which explains problems in schemes with

ACS Paragon Plus Environment

8

Page 9 of 22

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

The Journal of Physical Chemistry

the gas-phase parameters, is caused by usage of unrelaxed monomer geometries. For the considered dimers, the latter discrepancies become significant at separations shorter than the equilibrium, which will be rationalized by analysis of interaction energy components. Specifically, EFP Schemes 2 and 3 are more repulsive (less binding) than the CCSD(T) reference by more than 25 kcal/mol in all dimers. The only exception is the case of FaNN dimer in which Scheme 3 matches CCSD(T) values within 2.5 kcal/mol. However, this happens by a wrong reason: this EFP curve shows an unphysical kink and becomes attractive at very shorter distances. This unphysical behavior is probably produced by a failure of a multipole approximation at short separations. Clearly, Scheme 3 should not be used for any realistic modeling of FaNN dimer or related systems due to a potential collapse of molecular structures. Overall, at short separations, EFP schemes 1 and 4 are significantly closer to the CCSD(T) values than rigid-geometry Schemes 2 and 3, improving the latter ones by as much as 30 kcal/mol. SAPT computed using the gas-phase geometries of monomers underbinds considered dimers by 0.4-0.8 kcal/mol (at equilibrium), with discrepancies increasing at shorter separations to as much as 17 kcal/mol at 3.5 Å in FaNN dimer (see Table 1). SAPT computed at the relaxed monomer geometries (and corrected by the deformation energies at the CCSD(T)/CBS level) overestimates equilibrium binding energies by 0.7-1.5 kcal/mol, while discrepancies at 3.5 Å are within 3-6 kcal/mol. Thus, it is clear that the geometry relaxation is essential for accurate description of complexes at shorter than equilibrium distances. It is very encouraging that a similar behavior is observed between the EFP and SAPT methods, suggesting that a necessity to relax monomer geometries at close inter-monomer separations is a universal feature rather than an artifact of the EFP formalism. Another observation, identical for SAPT and EFP, is that the equilibrium binding energies computed using optimized geometries are always lower (more binding) than those computed using gas-phase monomer geometries. The second metrics that helps to understand accuracy of EFP at short separations is MAE computed for a whole separation range. Obviously, this metrics is dominated by short-range contributions where the absolute values of interaction energies and discrepancies in energies become significant. Comparing MAEs for different EFP schemes, we observe a very similar situation to what follows from comparing the interaction energies at short range. Namely, Schemes 2 and 3 show very large disagreement with the CCSD(T) reference, between 5 and 8 kcal/mol. On the other hand, Schemes 1 and 4 produce MAEs of around 1 kcal/mol for FaNN, 2 kcal/mol for FaON, and 3 kcal/mol for FaOO dimers. Again, similar to the case of the interaction energies, MAEs are the best for the FaNN dimer in which EFP predicts equilibrium geometry most accurately. Comparing Schemes 1 and 4, Scheme 1 produces better MAEs in two out of three cases, but the discrepancies between both schemes do not exceed 0.8 kcal/mol. Thus, we conclude that for describing strongly H-bonded systems at short intermolecular distances, it is essential to incorporate methods that account for internal changes in monomer geometries. Moreover, Scheme 4 with the gas-phase parameters shifted according to optimized geometry of the monomer produces very good agreement with Scheme 1 and with the reference CCSD(T) values. For comparison, optimized-geometry SAPT produces MAE in the range of 0.8-1.6 kcal/mol, while SAPT based on the gas-phase geometries is in error of 0.7-3.6 kcal/mol.

ACS Paragon Plus Environment

9

The Journal of Physical Chemistry

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

Page 10 of 22

3.4Å – 10.0Å

FaNN - FaNN

Figure 3. Total interaction energies and energy components (kcal/mol) in FaNN dimer.

ACS Paragon Plus Environment

10

Page 11 of 22

FaON

FaOO

Table 1. Equilibrium distances (Å), binding energies at equilibrium separation (kcal/mol), total interaction energies at 3.5 Å separation (kcal/mol) computed by different EFP schemes, SAPT2+3(CCD)/aug-cc-pVTZ, and the reference CCSD(T)/CBS method, as well as mean absolute deviations (MAE) of EFP and SAPT interaction energies with respect to CCSD(T)/CBS over the total range of considered inter-monomer separations (kcal/mol).

FaNN

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

The Journal of Physical Chemistry

EFP Scheme 1: optimized parameters EFP Scheme 2: gas phase parameters EFP Scheme 3: mixed basis parameters EFP Scheme 4: shifted parameters SAPT: optimized geometries SAPT: gas phase geometries CCSD(T)/CBS EFP Scheme 1: optimized parameters EFP Scheme 2: gas phase parameters EFP Scheme 3: mixed basis parameters EFP Scheme 4: shifted parameters SAPT: optimized geometries SAPT: gas phase geometries CCSD(T)/CBS

Equilibrium distance, Å

Binding energy, kcal/mol

Mean Interaction absolute error energy at 3.5 (MAE), Å, kcal/mol kcal/mol

3.98

-16.15

5.04

3.43

3.99

-15.78

17.58

5.73

4.02

-14.01

13.93

5.17

4.01

-14.97

6.04

3.63

3.72

-17.64

-15.66

1.00

3.84

-15.36

-8.49

0.72

3.81

-16.15

-11.42

4.12

-15.58

7.17

2.01

4.14

-15.46

40.51

8.21

4.21

-13.81

35.51

7.70

4.16

-14.29

9.66

2.54

4.00

-15.28

-4.26

0.76

4.06

-14.18

5.36

1.45

4.05

-14.60

-1.43

-16.26

4.22

1.34

-14.89

39.91

5.94

EFP Scheme 1: 4.18 optimized parameters EFP Scheme 2: gas 4.29 phase parameters

ACS Paragon Plus Environment

11

The Journal of Physical Chemistry

EFP Scheme 3: mixed basis parameters EFP Scheme 4: shifted parameters SAPT: optimized geometries SAPT: gas phase geometries CCSD(T)/CBS

4.30

-13.22

4.45

4.61

4.22

-15.04

8.71

0.54

4.14

-15.05

1.03

1.56

4.23

-13.54

24.19

3.63

4.21

-14.29

6.95

Energy components A more detailed understanding of different behavior in various EFP and SAPT schemes can be reached by analysis of interaction energy components. Energy component plots are shown in Figs. 1-3. Additionally, Table 2 provides comparison of numerical values of the interaction energy terms at equilibrium separations in all dimers. As is well known, hydrogen-bonding complexes are dominated by the electrostatic and polarization (induction) interactions counterpoised by the exchange-repulsion term, while dispersion plays a somewhat secondary role. This is true for three considered dimers, even though relative magnitudes of each term differ from one dimer to another.

FaOO

Table 2. Total interaction energies with and without deformation energy (DE), as well as energy components (kcal/mol) computed at the CCSD(T)/CBS equilibrium distances by different EFP schemes and SAPT.a

EFP Scheme 1: optimized parameters EFP Scheme 2: gas phase parameters EFP Scheme 3: mixed basis parameters EFP Scheme 4: shifted parameters SAPT: optimized geometries SAPT: gas phase geometries EFP Scheme 1: optimized parameters

Fa O N

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

Page 12 of 22

Electrostat. Energy, kcal/mol

ExchangeRepulsion Energy, kcal/mol

Dispersion Energy, kcal/mol

Induction Energy, kcal/mol

Total Interaction Energy with DE, kcal/mol

Total Interaction Energy without DE, kcal/mol

-31.81

29.36

-4.99

-9.92

-14.11

-17.36

-27.19b

27.00

-4.98

-8.18

-

-13.35

-27.35 b

30.94

-6.38

-8.47

-

-11.27

-30.01

28.09

-4.88

-8.93

-12.48

-15.73

-32.34

41.52

-10.36

-19.28

-17.47

-20.45

-28.21

38.52

-9.98

-15.69

-

-15.37

-26.26

21.40

-4.98

-8.12

-15.02

-17.96

ACS Paragon Plus Environment

12

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

The Journal of Physical Chemistry

FaNN

Page 13 of 22

a

EFP Scheme 2: gas phase parameters EFP Scheme 3: mixed basis parameters EFP Scheme 4: shifted parameters SAPT: optimized geometries SAPT: gas phase geometries EFP Scheme 1: optimized parameters EFP Scheme 2: gas phase parameters EFP Scheme 3: mixed basis parameters EFP Scheme 4: shifted parameters SAPT: optimized geometries SAPT: gas phase geometries

-22.00 b

18.69

-4.80

-6.47

-

-14.58

-22.22 b

23.06

-6.13

-6.69

-

-11.97

-24.50

20.43

-4.91

-7.27

-13.31

-16.25

-27.31

31.50

-8.76

-12.81

-15.28

-17.39

-23.26

27.41

-8.12

-10.16

-

-14.13

-27.80

23.77

-5.78

-9.30

-16.26

-19.12

-23.97 b

23.04

-5.89

-7.76

-

-14.58

-24.16 b

26.59

-7.28

-8.11

-

-12.96

-26.34

22.64

-5.67

-8.52

-15.04

-17.90

-27.60

32.81

-9.23

-13.15

-15.00

-17.16

-24.49

31.19

-9.04

-11.20

-

-13.54

computed at intermolecular separations of 3.8 Å for FaOO, 4.0 Å for FaON, 4.2 Å for FaNN;

b

differences in electrostatic energies between Schemes 2 and 3 are due to differences in the values of the overlap-based electrostatic damping, which is computed in the same basis as the exchange-repulsion term, i.e., 6-31+G* for Scheme 2 and 6-311++G(3df,2p) for Scheme 3.

Generally, one can notice that the EFP energy components computed with schemes 1 and 4 (both based on the relaxed fragment geometries) behave quite similarly and agree well with the relaxed SAPT energy components. In majority of cases, shifted-parameters Scheme 4 produces energy components that are slightly smaller by magnitude (i.e., less attractive for electrostatics, polarization, and dispersion and less repulsive for exchange-repulsion) than “optimized parameters” Scheme 1. However, for all energy terms except electrostatics, where Scheme 1 matches SAPT values very closely, the differences between Scheme 1 and Scheme 4 are smaller than the deviations of both schemes from the relaxed SAPT curves. Thus, energy component analysis justifies our conclusion based on comparing total interaction energies, that Scheme 4 with shifted parameters is a valid substitute for Scheme 1 and can be safely used when geometric flexibility of fragments is of importance. On the other hand, Schemes 2 and 3 with the gas phase parameters qualitatively match SAPT energy component curves with the unrelaxed geometries. Again, the EFP energy terms underestimate attractiveness of the electrostatics, polarization and dispersion interactions, and underestimate repulsiveness of the exchange-repulsion interactions, with respect to SAPT. This results in a partial error cancellation and an overall good accuracy of the EFP method. As expected, large basis set (in Scheme 3) contributes to increase in magnitude of polarization, dispersion and exchange-repulsion energies in all dimers. This increase is the biggest (in absolute values) for the exchange-repulsion term, but in relative values it is the most significant for the

ACS Paragon Plus Environment

13

The Journal of Physical Chemistry

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

Page 14 of 22

dispersion term. As a result, and in agreement with our previous work, the mixed-basis scheme produces less attractive total interaction energies than the small-basis schemes. It is noteworthy to compare energy component curves between the schemes with the optimized and gas-phase parameters. In all dimers, the schemes with relaxed geometries produce less repulsive (“softer”) exchange-repulsion curves at short separations. Similarly, the relaxedgeometry schemes result in less attractive (again “softer”) energy curves for electrostatics, polarization, and dispersion. These trends are the same for SAPT and EFP. Moreover, as was already discussed above, Scheme 4 with shifted parameters follows these trends as well. It is interesting though that in most cases, the gas-phase and relaxed-geometry energy component curves cross, and the crossing point occurs near or at slightly shorter distances than equilibrium separations (3.7 – 4.1 Å). This results in a peculiar effect that at the equilibrium, optimizedgeometry schemes produce more attractive values of electrostatics, polarization and dispersion, and more repulsive values of exchange-repulsion than the corresponding gas-phase-geometry schemes. Again, the trend is the same in SAPT and EFP. Indeed, even positions of the crossing points are very similar in both methods and do not differ by more than 0.1 Å in most cases, suggesting that the underlying physics is captured similarly by both methods. Turning now to Table 2, we can inspect numerical values of all EFP and SAPT schemes at equilibrium separations. The best agreement between the EFP and SAPT energy components in these dimers is observed for electrostatic terms, where the observed discrepancies are of the order of 1 kcal/mol, which is about 3-4% of the absolute values of electrostatic energies. The agreement is significantly worse for the other terms. For example, SAPT and EFP values for exchange-repulsion differ by 8-12 kcal/mol, which is about 30% of the exchange-repulsion energy, if the small basis set is employed for computing the exchange-repulsion parameters. If exchange-repulsion parameters are prepared in a larger basis, as in Scheme 3, the agreement with SAPT improves, such that deviations do not exceed 4-8 kcal/mol. For dispersion, the disagreement between EFP and SAPT schemes is 3-5 kcal/mol that constitutes as much as 50% of the total dispersion energy. Increase of the basis set in the EFP calculations diminishes the discrepancy in dispersion energies to 2-4 kcal/mol. In case of polarization, differences between EFP and SAPT at equilibrium are as big as 4-9 kcal/mol, or 30-50% of the total induction energy. Larger basis set in EFP (Scheme 3) produces only a slight improvement of the results, of the order of 0.5 kcal/mol compared to Scheme 2 with a small basis and gas phase parameters. Shifting of the EFP parameters: justification and analysis While the above discussion suggests that the EFP Scheme 4 with the gas-phase-based but shifted parameters provides a close agreement with the parameter-optimized Scheme 1, it is worth taking a closer look at explicit values of some of the parameters, in order to better understand a magnitude of the applied approximation. First, we analyze geometrical changes in the monomers that happen due to dimerization. Fig. 4a shows changes in the most important geometrical parameters of the monomers as a function of separation between the monomers in the dimer. It is clear the molecules undergo significant and non-linear geometrical rearrangements upon complexation. For example, changes in O-H and N-H bond lengths of the bonds involved in H-bonding are 0.05 Å and more in the considered range of intermonomer separations. Variations in O(N)-C-O(N) valence angles

ACS Paragon Plus Environment

14

Page 15 of 22

constitute 10-15 degrees. Fig. 4b illustrates relative changes in O(N)-H separation (this is a distance between a hydrogen on one monomer and an oxygen or nitrogen on another one, which are involved in a H-bond). At distances shorter than 3.8 Å (in FaOO dimer), 4.0 Å (in FaON dimer) or 4.2 Å (in FaNN dimer), we notice significant deviations between O(N)-H distances in fully optimized dimers vs dimers built from gas-phase monomers (positive values correspond to lengthening of O(N)-H distance in the optimized dimer with respect to the dimer constructed from gas-phase monomers). Comparing this observation with the analysis of energy components presented in Figs. 1-3, we notice that these distances match the distances at which optimized and gas-phase calculations start to differ. Thus, we conclude that differences between schemes with gas-phase and dimer-optimized monomers are governed by geometrical changes in monomer geometries.

Bond%Length,%Å%%%%%%%%%

aa

Angle,%°%%%%%%

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

The Journal of Physical Chemistry

b

Figure 4. (a) Variation of O-H (in FaOO) and N-H (in FaON and FaNN) bond lengths (in Å) of the bonds involved in H-bonding (black lines) and variation of C-O(N)-H and O(N)-C-O(N) angles (red and green curves, respectively, in degrees) as a function of separation between monomers in the dimer and in a gas-phase optimized monomer. (b) Changes in O(N)-H separation (in Å), i.e., the distance between the hydrogen on one monomer and the oxygen or nitrogen on another one that are involved in a H-bond, as a function of separation between monomers in the dimer. Positive values correspond to lengthening of O(N)-H distance in the dimer-optimized geometry with respect to the dimer geometry constructed from gas-phase monomers.

ACS Paragon Plus Environment

15

The Journal of Physical Chemistry

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

Page 16 of 22

In Fig. 5, we show the values of charges (Fig. 5a), dipoles (Fig. 5b) and polarizabilitities (Fig. 5c) at or near atoms directly involved in H-bonds, computed at various geometries of the monomers. That is, even though charges, dipoles, and polarizabilities are always computed for a monomer, the monomer’s geometry corresponds to that in the dimer characterized by a particular intermonomer distance. Specifically, the values of partial charges and magnitudes of dipole moments, obtained from the DMA, are considered at O, N, H atoms involved in H-bonding and at bond-mid-points between either O-H or N-H bonds. For this analysis, charges and dipoles at atoms also include neighboring charges/dipoles at bond-mid-points. Since distributed polarizability tensors are located at LMO centroids, two polarizabilities that contribute the most to H-bonding are considered. One polarizability corresponds to the N-H or O-H bond and is located approximately mid-way between N and H or O and H atoms. The other one corresponds to a lone-pair on N or O atoms that are acceptors of H-bond. For this analysis, isotropic values of polarizabilities were computed and plotted: 1 ߙത = ൫ߙ௫௫ + ߙ௬௬ + ߙ௭௭ ൯ 3 Note that full anisotropic polarizability tensors are used for computing polarization term in EFP. As follows from Fig. 5a, the values of partial charges do not change by more than 0.05 a.u. within the range of considered monomer geometries, which is a relatively small deviation to absolute values of partial charges of 0.3-0.7 a.u. The changes in dipole moments (Fig. 5b) are even smaller, with the only exception of the dipole on N atom in FaNN dimer, where the change is 0.1 a.u. A kink observed in that curve as well as in some other places originates in peculiarities of DMA that finds different solutions upon a small change in the geometry of the molecule. Values of the total dipole moments of the monomers are provided in Supporting Materials. Overall, analysis of partial charges and dipole moments is in agreement with the observation that the main reason in discrepancies in electrostatic energies between the gas-phase and relaxed fragments originates from changes in fragment geometries rather than fragment parameters. Turning to the analysis of polarizabilities (Fig. 5c), one can notice that polarizabilities on the lone pairs are not sensitive to geometrical changes. On the other hand, polarizabilities on OH and NH bonds show noticeable deviations with respect to the corresponding gas-phase values, of 0.20.3 a.u.3 or 15% in relative values. These deviations result in less than 10% discrepancy between the schemes with the optimized and shifted parameters, while the scheme with the gas-phase parameters produces dramatic failures in FaON and FaNN dimers. Dispersion interactions are computed using dynamic polarizabilities, with the main contribution originating from the zero-frequency (static) polarizability, such that we should expect qualitatively similar behavior due to applying approximate polarizabilities in polarization and dispersion terms. Indeed, we notice that discrepancies between optimized and shifted parameters schemes for dispersion are very small and do not exceed ~5% of the absolute values of the dispersion energies.

ACS Paragon Plus Environment

16

Page 17 of 22

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

The Journal of Physical Chemistry

a

b

c

Figure 5. (a) Partial charges, (b) dipole moments and (c) polarizabilities in atomic units at positions of atoms, bond-mid-points and LO centroids involved in H-bonding.

ACS Paragon Plus Environment

17

The Journal of Physical Chemistry

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

Page 18 of 22

Conclusions A question of flexibility of fragment geometries in the Effective Fragment Potential (EFP) method is considered on the example of strongly H-bonded dimers FaOO-FaOO, FaON-FaON, and FaNN-FaNN. EFP schemes in which parameters are pre-computed at the gas-phase fragment geometry or recomputed for each new geometry of the monomer were probed. Additionally, a scheme, in which gas-phase fragment parameters are shifted according to the optimized geometry of the monomer, is introduced. The EFP interaction energies along varying inter-monomer distances were compared with the reference CCSD(T)/CBS calculations and SAPT results. Similarly to EFP, two versions of SAPT were considered, one using relaxed fragment geometries and another one based on rigid gasphase geometries. Generally, all EFP and SAPT schemes provide good agreement with the reference CCSD(T) calculations near equilibrium geometries. However, the accuracy of the schemes that employ gas-phase monomer geometries dramatically deteriorates at shorter than equilibrium inter-monomer separations, while the schemes utilizing relaxed monomer geometries remain reliable in that region. The EFP scheme with shifted parameters closely follows the scheme in which fragment parameters are fully reoptimized for each unique fragment geometry. Analysis of EFP and SAPT energy components provides an interesting insight on the reasons of failure of rigid-geometry schemes at short intermonomer distances. Energy components in relaxed-geometry schemes show “softer” slopes toward decreasing intermonomer separations, as compared to rigid-geometry schemes in which energy component curves fall off steeply, which results in overly repulsive total potentials. As this behavior is consistent between the EFP and SAPT methods, we believe this is a general phenomenon rather than an artifact of a particular computational scheme. This study shows that the new EFP scheme in which fragment parameters are shifted according to the changes in fragment geometries is a reasonable substitute for the optimized-geometry scheme. This scheme with shifted parameters will be exploited and tested in more complex systems in future work. Supporting Information EFP deformation energies, graphs of the EFP interaction energies in 6-31+G*/aug-cc-pVTZ mixed basis, and EFP total dipole moments of the monomers are provided in Supporting Information section.

Acknowledgements LVS acknowledges support of the National Science Foundation under grants no. CHE-1465154 and CHE-1450088. ND and LVS acknowledge support by the Russian Science Foundation, project no. 14-43-00052. This research was supported in part through computational resources provided by Information Technology at Purdue University.

ACS Paragon Plus Environment

18

Page 19 of 22

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

The Journal of Physical Chemistry

References 1. Day, P. N.; Jensen, J. H.; Gordon, M. S.; Webb, S. P.; Stevens, W. J.; Krauss, M.; Garmer, D.; Basch, H.; Cohen, D. An Effective Fragment Method for Modeling Solvent Effects in Quantum Mechanical Calculations. J. Chem. Phys. 1996, 105, 1968-1986. 2. Gordon, M. S.; Freitag, M. A.; Bandyopadhyay, P.; Jensen, J. H.; Kairys, V.; Stevens, W. J. The Effective Fragment Potential Method: A QM-Based MM Approach to Modeling Environmental Effects in Chemistry. J. Phys. Chem. A 2001, 105, 293-307. 3. Gordon, M. S.; Slipchenko, L. V.; Li, H.; Jensen, J. H. The Effective Fragment Potential: A General Method for Predicting Intermolecular Forces. Annu. Rep. Comput. Chem. 2007, 3, 177-193. 4. Gordon, M. S.; Smith, Q. A.; Xu, P.; Slipchenko, L. V. Accurate First Principles Model Potentials for Intermolecular Interactions. Annu. Rev. Phys. Chem. 2013, 64, 553-578. 5. Ghosh, D.; Kosenkov, D.; Vanovschi, V.; Williams, C. F.; Herbert, J. M.; Gordon, M. S.; Schmidt, M. W.; Slipchenko, L. V.; Krylov, A. I. Noncovalent Interactions in Extended Systems Described by the Effective Fragment Potential Method: Theory and Application to Nucleobase Oligomers. J. Phys. Chem. A 2010, 114, 12739-12754. 6. Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation Methods: A Route to Accurate Calculations on Large Systems. Chem. Rev. 2011, 112, 632-672. 7. Ghosh, D.; Kosenkov, D.; Vanovschi, V.; Flick, J.; Kaliman, I.; Shao, Y.; Gilbert, A. T. B.; Krylov, A. I.; Slipchenko, L. V. Effective Fragment Potential Method in Q-Chem: A Guide for Users and Developers. J. Comput. Chem. 2013, 34, 1060-1070. 8. Smith, Q. A.; Gordon, M. S.; Slipchenko, L. V. Benzene−Pyridine Interactions Predicted by the Effective Fragment Potential Method. J. Phys. Chem. A 2011, 115, 4598-4609. 9. Flick, J. C.; Kosenkov, D.; Hohenstein, E. G.; Sherrill, C. D.; Slipchenko, L. V. Accurate Prediction of Non-Covalent Interaction Energies with the Effective Fragment Potential Method: Comparison of Energy Components to Symmetry-Adapted Perturbation Theory for the S22 Test Set. J Chem. Theory Comp. 2012, 8, 2835-2843. 10. Kemp, D. D.; Gordon, M. S. Theoretical Study of the Solvation of Fluorine and Chlorine Anions by Water. J. Phys. Chem. A 2005, 109, 7688-7699. 11. Kemp, D. D.; Gordon, M. S. An Interpretation of the Enhancement of the Water Dipole Moment Due to the Presence of Other Water Molecules. J. Phys. Chem. A 2008, 112, 4885-4894. 12. Brorsen, K. R.; Pruitt, S. R.; Gordon, M. S. Surface Affinity of the Hydronium Ion: The Effective Fragment Potential and Umbrella Sampling. J. Chem. Phys. B 2014, 118, 1438214387. 13. Hands, M. D.; Slipchenko, L. V. Intermolecular Interactions in Complex Liquids: Effective Fragment Potential Investigation of Water–Tert-Butanol Mixtures. J. Phys. Chem. B 2012, 116, 2775-2786. 14. Rankin, B. M.; Hands, M. D.; Wilcox, D. S.; Fega, K. R.; Slipchenko, L. V.; Ben-Amotz, D. Interactions between Halide Anions and a Molecular Hydrophobic Interface. Faraday Discuss. 2013, 160, 255-270. 15. Adamovic, I.; Gordon, M. S. Methanol−Water Mixtures:  A Microsolvation Study Using the Effective Fragment Potential Method. J. Phys. Chem. A 2006, 110, 10267-10273. 16. Slipchenko, L. V.; Gordon, M. S. Water-Benzene Interactions: An Effective Fragment Potential and Correlated Quantum Chemistry Study. J. Phys. Chem. A 2009, 113, 2092-2102.

ACS Paragon Plus Environment

19

The Journal of Physical Chemistry

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

Page 20 of 22

17. Smith, Q. A.; Gordon, M. S.; Slipchenko, L. V. Effective Fragment Potential Study of the Interaction of DNA Bases. J. Phys. Chem. A 2011, 115, 11269-11276. 18. Arora, P.; Slipchenko, L. V.; Webb, S. P.; Defusco, A.; Gordon, M. S. Solvent-Induced Frequency Shifts: Configuration Interaction Singles Combined with the Effective Fragment Potential Method. J. Phys. Chem. A 2010, 114, 6742-6750. 19. DeFusco, A.; Ivanic, J.; Schmidt, M. W.; Gordon, M. S. Solvent-Induced Shifts in Electronic Spectra of Uracil. J. Phys. Chem. A 2011, 115, 4574-4582. 20. DeFusco, A.; Minezawa, N.; Slipchenko, L. V.; Zahariev, F.; Gordon, M. S. Modeling Solvent Effects on Electronic Excited States. J. Phys. Chem. Lett. 2011, 2, 2184-2192. 21. Minezawa, N.; Gordon, M. S. Optimizing Conical Intersections of Solvated Molecules: The Combined Spin-Flip Density Functional Theory/Effective Fragment Potential Method. J. Chem. Phys. 2012, 137, 034116. 22. De Silva, N.; Minezawa, N.; Gordon, M. S. Excited-State Hydrogen Atom Transfer Reaction in Solvated 7-Hydroxy-4-Methylcoumarin. J. Chem. Phys. B 2013, 117, 15386-15394. 23. Kosenkov, D.; Slipchenko, L. V. Solvent Effects on the Electronic Transitions of PNitroaniline: A Qm/Efp Study. J. Phys. Chem. A 2010, 115, 392-401. 24. Slipchenko, L. V. Solvation of the Excited States of Chromophores in Polarizable Environment: Orbital Relaxation Versus Polarization. J. Phys. Chem. A 2010, 114, 8824-8830. 25. Ghosh, D.; Isayev, O.; Slipchenko, L. V.; Krylov, A. I. The Effect of Solvation on Vertical Ionization Energy of Thymine: From Microhydration to Bulk. J. Phys. Chem. A 2011, 115, 6028–6038. 26. Gurunathan, P. K.; Acharya, A.; Ghosh, D.; Kosenkov, D.; Kaliman, I.; Shao, Y.; Krylov, A. I.; Slipchenko, L. V. Extension of the Effective Fragment Potential Method to Macromolecules. J. Chem. Phys. B 2016, 120, 6562-6574. 27. Fedorov, D. G.; Kitaura, K. Extending the Power of Quantum Chemistry to Large Systems with the Fragment Molecular Orbital Method. J. Phys. Chem. A 2007, 111, 6904-6914. 28. Fedorov, D. G.; Slipchenko, L. V.; Kitaura, K. Systematic Study of the Embedding Potential Description in the Fragment Molecular Orbital Method. J. Phys. Chem. A 2010, 114, 8742-8753. 29. Komeiji, Y.; Ishikawa, T.; Mochizuki, Y.; Yamataka, H.; Nakano, T. Fragment Molecular Orbital Method-Based Molecular Dynamics (Fmo-Md) as a Simulator for Chemical Reactions in Explicit Solvation. J. Comput. Chem. 2009, 30, 40-50. 30. The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems. Fedorov, D. G.; Kitaura, K., Eds.; CRC Press: 2009. 31. Steinmann, C.; Fedorov, D. G.; Jensen, J. H. Effective Fragment Molecular Orbital Method: A Merger of the Effective Fragment Potential and Fragment Molecular Orbital Methods. J. Phys. Chem. A 2010, 114, 8705-8712. 32. Pruitt, S. R.; Steinmann, C.; Jensen, J. H.; Gordon, M. S. Fully Integrated Effective Fragment Molecular Orbital Method. J. Chem. Theory Comput. 2013, 9, 2235-2249. 33. Bertoni, C.; Gordon, M. S. Analytic Gradients for the Effective Fragment Molecular Orbital Method. J. Chem. Theory Comput. 2016, 12, 4743–4767. 34. Thanthiriwatte, K. S.; Hohenstein, E. G.; Burns, L. A.; Sherrill, C. D. Assessment of the Performance of Dft and Dft-D Methods for Describing Distance Dependence of HydrogenBonded Interactions. J. Chem. Theory Comput. 2011, 7, 88-96.

ACS Paragon Plus Environment

20

Page 21 of 22

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

The Journal of Physical Chemistry

35. Parker, T. M.; Burns, L. A.; Parrish, R. M.; Ryno, A. G.; Sherrill, C. D. Levels of Symmetry Adapted Perturbation Theory (Sapt). I. Efficiency and Performance for Interaction Energies. J. Chem. Phys. 2014, 140, 094106. 36. Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J., et al. General Atomic and Molecular Electronic-Structure System. J. Comp. Chem. 1993, 14, 1347-1363. 37. Gordon, M. S.; Schmidt, M. W. Advances in Electronic Structure Theory: Gamess a Decade Later. In Theory and Applications of Computational Chemistry Dykstra, C. E.; Frenking, G.; Kim, K. S.; Scuseria, G. E., Eds. Elsevier: 2005; pp 1167-1190. 38. Slipchenko, L. V.; Gordon, M. S. Damping Functions in the Effective Fragment Potential Method. Mol. Phys. 2009, 107, 999-1016. 39. Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P. The Polarizable Atomic Multipole-Based Amoeba Force Field for Proteins. J. Chem. Theory Comput. 2013, 9, 4046-4063. 40. Jensen, J. H.; Gordon, M. S. An Approximate Formula for the Intermolecular Pauli Repulsion between Closed Shell Molecules. Mol. Phys. 1996, 89, 1313-1325. 41. Jensen, J. H.; Gordon, M. S. An Approximate Formula for the Intermolecular Pauli Repulsion between Closed Shell Molecules. Ii. Application to the Effective Fragment Potential Method. J. Chem. Phys. 1998, 108, 4772-4782. 42. Stone, A. J. The Theory of Intermolecular Forces. Oxford University Press: Oxford: 1996. 43. Hariharan, P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theor. Chim. Acta 1973, 28, 213-222. 44. Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. V. R. Efficient Diffuse Function-Augmented Basis Sets for Anion Calculations. III. The 3-21+G Basis Set for First-Row Elements, Li–F. J. Comput. Chem. 1983, 4, 294-301. 45. Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Self Consistent Molecular Orbital Methods. XX. A Basis Set for Correlated Wave Functions. J. Chem. Phys. 1980, 72, 650-654. 46. Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self Consistent Molecular Orbital Methods 25. Supplementary Functions for Gaussian Basis Sets. J. Chem. Phys. 1984, 80, 3265-3269.

ACS Paragon Plus Environment

21

The Journal of Physical Chemistry

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

Page 22 of 22

TOC Graphic

ACS Paragon Plus Environment

22