Importance of Accurate Structures for Quantum Chemistry Embedding

Subscriber access provided by NEW MEXICO STATE UNIV

Spectroscopy and Excited States

On the importance of accurate structures for quantumchemistry embedding methods: which strategy is better? Erik Rosendahl Kjellgren, Jógvan Magnus Haugaard Olsen, and Jacob Kongsted J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00202 • Publication Date (Web): 03 Jul 2018 Downloaded from http://pubs.acs.org on July 6, 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 32 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

On the importance of accurate structures for quantum-chemistry embedding methods: which strategy is better? Erik Rosendahl Kjellgren,† Jógvan Magnus Haugaard Olsen,†,‡ and Jacob Kongsted∗,† †Department of Physics, Chemistry and Pharmacy, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark ‡Hylleraas Centre for Quantum Molecular Sciences, Department of Chemistry, UiT The Arctic University of Norway, N-9037 Tromsø, Norway E-mail: [email protected] Abstract Quantum chemistry embedding methods have become a popular approach to calculate molecular properties of larger systems. In order to account for finite temperature effects — including both configurational and conformational averaging — embedding methods are often combined with molecular dynamics (MD) simulations either in a direct or sequential manner. One of the decisive factors for a successful application of embedding methods is that that the underlying structures provided by the MD simulation are accurate – if not this will result in low-quality prediction of the molecular properties in question. Here we investigate different approaches for generating a set of molecular structures to be used in subsequent embedding calculations ranging from classical MD using a standard molecular mechanics (MM) force field to combined quantum mechanics/molecular mechanics (QM/MM) MD. Overall, we find an

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

intermediate approach relying on classical MD followed by a constrained QM/MM geometry optimization to be a fairly accurate and very cost-effective approach although this procedure naturally leads to underestimation of e.g. spectral bandwidths.

1

Introduction

In recent years methods based on quantum–classical approaches have become an important part of the computational chemists toolbox. In broad terms these quantum–classical methods may be divided into two groups based on the physical approach taken for the description of the classically treated part of the system. In the first approach this part is described based on a dielectric continuum model, 1,2 while in the second approach a discrete representation — based on molecular mechanics — is used. 3,4 In both cases the classical description is coupled to a quantum mechanical description of the smaller part of the total system, and either ways of describing the classical region possesses pros and cons. However, the great advantage of such quantum–classical approaches is their reduced scaling with respect to system size as compared to a regular full quantum-mechanical treatment of the whole system, which quickly becomes impossible or in the best case very expensive from a computational point of view. Thus, quantum–classical approaches enable efficient calculations of all kinds of chemical systems where a division into a core region (treated using quantum mechanics) and its environment (treated classically) is possible. One example of such a system is the case of a solute in a solvent. For this, one is usually not interested in evaluating all the details of the electronic structure of the solvent, but rather the solute is of interest. From a computational point of view the dielectric continuum description of such a system is very appealing since in this case it is possible to introduce the effect of the solvent in one single calculation. This is, on the other hand, very different from the explicit solvent models. The reason for this is that in the explicit solvent models averaging of the solvent configurations have to be included explicitly, which often requires coupling to a sampling technique such as molecular dynamics (MD) simulations. Thus, compared to the dielectric continuum method, the use 2


Page 2 of 32

Page 3 of 32 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 an explicit solvent model is much more involved since in this case there is a need to couple to MD simulations generating the structures to be used for the property calculations of the solute-solvent system. Often the generation of structures is decoupled from the quantum– classical calculations in the way that one first performs the MD simulation and then extracts a sufficient number of structures to be used in the following quantum–classical calculations. This approach is often denoted as the sequential procedure. 5 If the aim of a given study is to calculate spectroscopic properties of the solute, e.g. properties of either electric or magnetic character or to address excited states and their properties, the underlying quality of the structural parameters of especially the solute molecule becomes of utmost importance, and can significantly influence the final spectroscopic properties, see for example refs 6 and 7 for recent discussions. At the same time the MD simulations need to be run for a sufficient time to ensure adequate sampling. While it is in principle possible to satisfy both of these requirements, it will in practice lead to very expensive calculations. A possible solution to this is to rely on specially optimized force fields. 6,8 However, the procedure of deriving such optimized parameters is usually rather involved and needs to be repeated for each individual system. Thus, this way of proceeding will add a further layer of complexity to the overall computational approach. Another alternative is to use quantum-mechanical force fields that are based on efficient fragmentbased quantum-mechanical methods which thus circumvent the scaling issue of conventional quantum-mechanical methods. 9–29 While such methods are very promising, they are still being actively developed and have not yet reached a level of maturity that allows straightforward application of the methods. 30 In this paper we will investigate the possibility of using structures for quantum–classical embedding calculations that are based on automated procedures without the need for performing any force field adjustment or re-parametrization. The three methods explored are based on either classical MD using a general molecular mechanics (MM) force field, combined quantum mechanics / molecular mechanics (QM/MM) MD, 3,4,31,32 or an intermediate

3



Page 4 of 32

approach where structures from the classical MD are geometry optimized based on a constrained QM/MM approach. Using the structures extracted from the QM/MM MD as a reference, we study both differences in the produced structures as well as how these differences in the structures affect calculation of molecular properties of the solute based on embedding calculations using the polarizable embedding (PE) model. 33 In this context we will consider two different kinds of molecular properties as exemplifying cases namely excited electronic states (excitation energies and corresponding oscillator strengths) and NMR shielding constants. We have chosen these two different molecular properties since these first of all are rather sensitive to the quality of the underlying molecular structure and since these properties represent common spectroscopic quantities to be obtained from quantum chemical embedding approaches. We note that the intermediate approach relying on partial geometry optimization – in this case of only the solute – naturally removes some of the thermal effects introduced by the MD and one central question of this paper is to quantify to what degree such a removal of thermal effects will impact the final spectroscopic properties. We note that the approach relying on (partial) geometry optimization of structures generated by classical MD has been found useful in other types of applications, e.g. for post-simulation optimizations in relation to pKa predictions, 34 for calculating ligand binding affinities 35 and for inclusion of the vibrational structure in spectral shapes of molecules in solution. 36 The specific molecules studied are shown in Figure 1 (PNP1-3) and consist of a set of three molecules of increasing complexity regarding structural parameters. H2 O2 N

C2

H1 C1

C3

O3

C4

C6 C5

H3

PNP1

H4

H3

H2 C3

O2 O1

N

C5 C1

H1

H6 C8

C4

C2

O3

H5 C9

H4

H7

H8

C3

O3

C12 C11

C6

H3

H2 C10 O1

C7

H8

H7

N O2

PNP2

C8

C4 C5

C2 C1 H1

C10 C11

H6

C14

C7 H5

C6 H4

C12 O1

C9

H10

C13 H9

PNP3

Figure 1: Structures of the investigated molecules: 4-nitrophenolate (PNP1), 4’-nitro-4biphenylate (PNP2), and 4-[(E)-2-(4-nitrophenyl)vinyl]phenolate (PNP3). 4


Page 5 of 32 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


This paper is organized in the following way. In section 2 we describe the computational procedures in detail. Section 3 contains the analysis of the structural and embedding results for PNP1 followed by analysis of the corresponding results for PNP2 and PNP3. Finally, section 4 contains the conclusions and recommendations based on the main findings of this work.

2 2.1

Computational Details Classical MD simulations

The partial charges used in the force-field description of PNP1-3 were derived following a standard protocol, i.e. the molecular geometries were first optimized at the HF/6-31G* 37,38 level followed by an electrostatic potential (ESP) calculation at the same level of theory. For these calculations Gaussian09 39 was used. The atomic charges were assigned from the calculated ESP by a restrained electrostatic potential 40 (RESP) calculation using Antechamber. 41 Solvation of the PNP1-3 molecules in explicit water was performed using xleap. The solutes were described using GAFF 42 while the water solvent was modeled using TIP3P. 43 The charge of the total system was neutralized with sodium ions. In the solvation process, water molecules were placed around the solute molecule in a cube such that all atoms in the solute molecule were at least 13.0 ˚ A away from the box edge. The MD simulations were performed using Amber. 44 To prepare for the MD simulations, a minimization was performed with a non-bonded interaction cutoff of 10.0 ˚ A. The subsequent equilibration was performed in two steps. The first equilibration was run in the NPT ensemble, and the second equilibration was run in the NVT ensemble. The final production step was also run in the NVT ensemble. The Andersen thermostat 45 and Berendsen barostat 46 were used. Each equilibration step ran for 400 ps while the production step ran for 1.2 ns, all using a 2.0 fs time step. All hydrogen atoms were constrained using SHAKE. 47 Particle mesh Ewald (PME) 48 was used in all simulations for the electrostatics. For further analysis a snapshot was taken out every 5



Page 6 of 32

5000’th step from the production run using CPPTRAJ. 49

2.2

QM/MM geometry optimizations

From each of the classical MD simulations, 120 snapshots were extracted at regular time intervals. Based on these snapshots, QM/MM geometry optimizations of the solute molecules PNP1-3 within a fixed MM region were performed using QSite. 50,51 In these QM/MM optimizations the QM region consisted of the solute molecule while the MM region consisted of the water molecules and sodium ion. The QM region was treated using the BLYP 52,53 functional combined with the 6-31G* basis set. The OPLS2005 54 force field was used for the MM region. For PNP1 a few of the snapshots had the sodium ion and a few of the water molecules included in the quantum region according to a distance criterion set as default by the Schrödinger package. This included water and the sodium ion, if hydrogen bonds connected these to the PNP1 molecule. For PNP1 the 400 water molecules closest to the solute were included in the MM region and the rest was not considered in the calculation. The closest water molecules were defined by the distance between the atoms in the PNP molecule and the oxygen atoms of the water molecules. For PNP2 the 470 closest water molecules were included and for PNP3 the closest 610 water molecules were included.

2.3

QM/MM Born-Oppenheimer MD simulations

The QM/MM Born-Oppenheimer MD (BOMD) simulations were performed using the CPMD 55–57 code interfaced with the GROMOS 58 program. A total of 12 different starting structures were acquired from the classical MD simulation and a trajectory was produced from each of these starting structures. The starting structures were taken equally distributed from the classical MD trajectory (0.12 ns between each structure). The solutes were treated as the quantum region modeled by Kohn-Sham density functional theory (DFT) using a planewaves basis (energy cutoff at 80.0 Ry) with Martin-Troullier pseudopotentials 59 and the BLYP functional. The MM region consisted of all water molecules and the sodium which 6


Page 7 of 32 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


were described by the same force field as used in the classical MD simulation. The time step was set to 10.0 a.u. and 44000 steps were taken providing a trajectory of 10.6 ps. A Nosé-Hoover thermostat 60,61 with a frequency of 3400.0 cm−1 was used for the QM region, while for the MM region a frequency of 3000.0 cm−1 was used. Both thermostats were used to maintain a temperature of 298.0 K. The Hockney Poisson solver 62 was used to mimic an isolated system. The cell size of the quantum system was equal in all directions (the cell size was 30.4 a.u. for PNP1, 38.7 a.u. for PNP2 and 43.8 a.u. for PNP3). The electrostatic coupling between the QM and MM region is taken explicitly into account for MM atoms within 20.0 a.u. from the QM region, while atoms further use the long-range electrostatic coupling of Laio et al. 63

2.4

Excitation energies and oscillator strengths

To calculate excitation energies and associated oscillator strengths, the polarizable embedding 33 (PE) model was used to describe environment effects. In brief the PE model accounts for environmental effects by assigning to the solvent molecules a set of distributed multipoles and polarizabilities giving rise to an embedding potential to be included in the optimization of the electronic structure of the solute molecule. The calculations were performed with the solute molecule (PNP1-3) in the quantum region, and the classical region was chosen to be the water molecules and potentially a sodium ion. The embedding-potential parameters were calculated using the CAM-B3LYP 64 functional and 6-31+G* basis set. The calculations were performed at the DFT level also using the CAM-B3LYP functional and the 6-31+G* basis set. These PE-DFT calculations were performed using the Dalton program 65,66 employing PElib 67 and Gen1Int. 68,69 The generation of spectra was performed on the basis of the 5 lowest-lying excited states employing eq. (SI9) in the supporting information. We note that in this generation of absorption spectra we do not have an explicit account of the vibronic structure. Recently, different strategies have been discussed for inclusion of vibrational effects into band shapes based on combinations of explicit and implicit solvent models. Even if 7



these schemes usually rely on the approximation of treating the solute rigid in the MD simulation (see e.g. Refs. 70 and 71) or to rely on a partial geometry optimization of the solute in a frozen solvent cage 36 such approaches can surely be seen as more complete compared to the route taken in our study since in the following we completely neglect any account of vibrational structure in the predicted absorption spectra.

2.5

NMR shielding constants

For the calculation of NMR shielding constants the definition of the classical region as well as the construction of the embedding potentials follow the same strategy as used for the calculation of the excitation energies. The calculation of the NMR shielding constants was performed at the DFT level using the KT3 72 functional and the pcS-1 73,74 basis set. The calculations were performed using the Dalton program employing PElib and Gen1Int.

3

Results and Discussion

In this section, we present and discuss the results of the calculations performed for PNP1-3. In the first part we focus on PNP1, which is easier to analyze in detail because of its simpler molecular structure, and thereafter we move on to discuss PNP2 and PNP3.

3.1

PNP1

In Figure 2, we present the results of the five lowest excitation energies of PNP1 in aqueous solution obtained from calculations on structures sampled from the classical (MM) MD trajectory, the QM/MM geometry optimizations of the structures from the MM-MD trajectory, or the structures from the QM/MM-MD trajectories. The classical MD structures were sampled every 10 ps, and the QM/MM-MD structures sampled every 1 ps. The excitation energies are given both as a function of snapshot number (Fig. 2a) and as histograms (Fig. 2b). It can immediately be seen that the structures sampled by MM-MD give rise to 8


Page 8 of 32

Page 9 of 32

000' 4000RSWLPL]HG 40000'

&RXQW

([FLWDWLRQHQHUJ\>H9@

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


6QDSVKRW

([F ([F 000' ([F ([F 4000RSWLPL]HG 40000' ([FLWDWLRQHQHUJ\>H9@

(a)

([F

(b)

Figure 2: Excitation energies of PNP1 calculated using either MM-MD structures, QM/MMoptimized structures, or QM/MM-MD structures. (a) shows the excitation energies of PNP1 as function of snapshot number and (b) shows the distribution of the excitation energies on the basis of frequency.

9



a much higher spread of the excitation energies compared to the QM/MM-optimized and QM/MM-MD structures. The lowest excitation energy seems to be especially sensitive to the quality of the structures. Here the excitation energies in the MM-MD structures range from around 1.5 to 3.5 eV, whereas the excitation energies in the other structures are clustered around 3.0 to 3.5 eV. The lowest spreads are observed for the QM/MM-optimized structures. This is not unexpected considering that the geometry optimization will remove finite-temperature effects within the PNP1 molecule. The solvent effect on the geometry is, however, still retained to some degree because the water solvent is kept frozen during the geometry optimization. In the following, we will focus on the lowest intense π → π ∗ transition which is dominated by a transition from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO). This transition is characterized as an intramolecular charge-transfer (ICT) transition where electron charge is transferred from the phenolate to the nitro group. The standard deviation of this excitation energy is calculated to be 0.37 eV based on the MM-MD structures compared to only 0.12 and 0.16 eV for the QM/MMoptimized and QM/MM-MD structures, respectively. This is a rather substantial difference and to investigate this discrepancy further we turn our attention to the oscillator strengths. The oscillator strengths are presented in Figure 3 for each snapshot (Fig. 3a) and as a function of the dihedral angle between planes defined by the nitro- and phenolate groups (Fig. 3b).1 First we examine the variation of the oscillator strength with different snapshots shown in Figure 3a. Just as observed for the excitation energies, we here also find larger fluctuations of the oscillator strengths based on the MM-MD structures compared to those obtained from either QM/MM-optimized and QM/MM-MD structures although the difference is less pronounced here. More interesting perhaps is the fact that the oscillator strengths from the MM-MD structures are noticeably lower and even close to zero for some snapshots. In terms of numbers, we calculated the average MM-MD oscillator strength to 1

Note that in Fig. 3b we always use the lowest intense transition which is the second state for some of the QM/MM-MD structures as seen in Fig. 3a.

10


Page 10 of 32

VW([F

2VFLOODWRUVWUHQJWK

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


QG([F

UG([F

000' 4000RSWLPL]HG 40000' 6QDSVKRW

WK([F

WK([F

2VFLOODWRUVWUHQJWK

Page 11 of 32

000'

4000RSWLPL]HG

40000'

(a)

&&12GLKHGUDODQJOH

(b)

Figure 3: Panel (a) shows the oscillator strength of PNP1 for the 120 snapshots. In the MM-MD and the QM/MM-MD there is change in order of the excited states in some of the snapshots, i.e. the first and the second excitations are changed. Panel (b) shows the oscillator strength of the lowest excitation against the dihedral angle C2-C3-N-O2, see Fig. 1. In panel (b) the transitions have been ordered such that the lowest excitation is defined as the excitation with the largest oscillator strength.

11



be 0.32 compared to 0.50 and 0.48 for QM/MM-optimized and QM/MM-MD, respectively. The observed differences in excitation energy and oscillator strengths can be explained by (at least) two key structural features. GAFF assigns the same atom-type to all carbons in the phenyl group thus leading to almost equal average carbon-carbon bond lengths in the MM-MD structures as can be seen in Figure 4. This is in contrast to QM/MM-MD which predicts different average bond lengths for the different carbon-carbon bonds. It is also apparent from Figure 4 that the QM/MM optimization leads to bond-length alternations that are consistent with QM/MM-MD, although nearly all fluctuations are removed as expected.

&RXQW

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 32

000' 4000RSWLPL]HG 40000'

&&

&&

&&

%RQGOHQJWK>c@

Figure 4: Histograms for the three different C-C bond lengths in PNP1. The results are based on statistics from 120 snapshots. Another important structural feature is the torsion of the nitro group. As the torsion increases, the conjugation will decrease, thus, affecting both the excitation energy and oscillator strength. Specifically, we expect that the oscillator strength will decrease for increasing torsion. The degree of the torsion can be measured by the dihedral angle between the planes defined by the phenyl and nitro group, respectively. This torsion has a rather big effect 12


Page 13 of 32

on the oscillator strength as can be seen in Figure 3b where the dihedral angle is plotted against the oscillator strength. As expected, we find a negative correlation between the torsion angle of the nitro group and the oscillator strength. It is evident that the dihedral angle term is too flexible in the MM-MD where it ranges from 0 to 90 degrees, i.e. full rotation seems to occur, whereas the QM/MM-MD structures only vary between 0 and 30 degrees. Furthermore it can be seen that some of the very low oscillator strengths observed for the MM-MD structures mainly occur at large torsions of the nitro group. Again we find that the QM/MM optimization is largely able to correct this flaw, although there are still some dihedral angles present that are not explored in the QM/MM MD.

$EVRUSWLRQ

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


000' 4000RSWLPL]HG 40000'

313


Figure 5: Predicted absorption spectra for PNP1 in water solution. The blue lines represent the different snapshots based on the 5 lowest excitations. The red curve is a calculated spectrum from the blue sticks. A normalized Gaussian has been imposed on each blue stick, with σ = 0.1 eV, resulting in the final convolution. The MM-MD has two low points around 1.5 eV that are not shown. Finally, we compare the absorption spectra produced from the calculations on MM-MD, QM/MM-optimized, or QM/MM-MD structures, respectively (see Fig. 5). The spectra are

13



produced by Gaussian convolution using a broadening factor of 0.1 eV, see the Supporting Information eq. (SI9) for more details. As expected, based on the analysis above, the MMMD structures result in a spectrum with much lower maximum absorption compared to the two other types of structures. This is partly because of the generally lower oscillator strengths and partly because of the larger span of the excitation energies. The QM/MMoptimized structures yield a spectrum with a peak that is somewhat more intense than the QM/MM-MD spectrum owing to the fact that the excitation energies are clustered more tightly. The energies at maximum absorption are found to be 3.31, 3.37, and 3.28 eV, for the MM-MD, QM/MM-optimized, and QM/MM-MD structures, respectively. If the excitation energies are instead obtained as averages based on the calculated excitation energies we obtain 3.15, 3.35, and 3.29 eV, respectively. Thus, the excitation energies are almost the same for the QM/MM-optimized and QM/MM-MD structures regardless of how they are obtained. The average excitation energy for the MM-MD structures, on the other hand, is somewhat lower than the maximum absorption energy found from the spectrum. This indicates that there is some correlation between low excitation energy and low oscillator strength for these structures. The excitation energy based on the QM/MM-optimized structures is higher than the excitation energy obtained from the QM/MM-MD structures.

3.2

PNP2 and PNP3

Histograms of the five lowest excitation energies in PNP2 and PNP3 based on MM-MD, QM/MM-optimized, or QM/MM-MD structures are presented in Fig. 6. Just as for PNP1 we observe a large spread of the excitation energies based on the MM-MD structures and a significantly lower spread of the excitation energies based on QM/MM-optimized structures. However, contrary to the PNP1 case, we here see that the spread of the excitation energies based on QM/MM-MD structures is closer to the spread of the MM-MD excitation energies. Indeed, the standard deviations of the lowest intense excitation energy in PNP2 are 0.40, 0.21, and 0.37 eV, for the MM-MD, QM/MM-optimized, and QM/MM-MD structures, 14


Page 14 of 32

313 ([F ([F 000' ([F ([F 4000RSWLPL]HG 40000' ([FLWDWLRQHQHUJ\>H9@

([F

&RXQW

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


&RXQW

Page 15 of 32

313 ([F ([F 000' ([F ([F 4000RSWLPL]HG 40000' ([FLWDWLRQHQHUJ\>H9@

(a)

([F

(b)

Figure 6: (a) Excitation energy distribution for the 120 snapshots of PNP2. (b) Excitation energy distribution for the 120 snapshots of PNP3.

15



respectively, and in PNP3 they are 0.40, 0.16, and 0.38 eV, for the MM-MD, QM/MMoptimized, and QM/MM-MD structures, respectively. Thus it appears that the QM/MM

000'

4000RSWLPL]HG

(a)

&&&&GLKHGUDODQJOH (b)

40000'

2VFLOODWRUVWUHQJKWV

optimization produces structures that result in a too narrow spread of the excitation energies.

&&12GLKHGUDODQJOH

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 32

(c)

Figure 7: Shows the two dihedral angles C3-C2-N-O2 and C8-C7-C5-C4 for the 120 snapshots of PNP2. The color scale shows the oscillator strength for the lowest and intense excitation (excitations ordered by magnitude of oscillator strengths). Continuing to the oscillator strengths we expect that internal torsions also have a large effect on the oscillator strengths for PNP2 and PNP3 just as we observed for PNP1. For this analysis, we will only consider PNP2 which has two important torsions, namely the torsion of the nitro group as in PNP1, and the additional torsion of the nitrophenyl group relative to the phenolate group. The heatmap plots in Figure 7 show the oscillator strength of the lowest excitation energy as a function of the two dihedral angles. It is clear from these plots that the MM-MD samples a wide range of dihedral angles that are not visited in the QM/MM-MD, and, furthermore, that the QM/MM optimization corrects this flaw to some degree. It can also be seen that the oscillator strength falls rapidly towards zero when any of the dihedral angles exceeds about 40◦ . From the discussion of PNP1 this reduction in oscillator strength as a function of these dihedral angels could be expected and would also be expected to apply for PNP3. In the region below 30◦ for both dihedrals, we find that the oscillator strengths based on the MM-MD structures are lower than those obtained from the QM/MM-optimized and QM/MM-MD structures similar to what we observed for PNP1. 16


Page 17 of 32

This suggests that there are also other structural differences that have a significant impact on the oscillator strength, e.g. the aromatic carbon-carbon bond lengths.

313

000'

4000RSWLPL]HG

$EVRUSWLRQ

$EVRUSWLRQ

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


40000'


313 000'

4000RSWLPL]HG

40000'

(a)


(b)

Figure 8: Predicted absorption spectra for PNP2 and PNP3 in water solution. The blue lines represent the different snapshots based on the 5 lowest excitations. The red curve is a calculated spectrum from the blue sticks. A normalized Gaussian has been imposed on each blue line, with σ = 0.1 eV. Figures 8a and 8b display the calculated spectra of PNP2 and PNP3, respectively, for the three different types of structures. The main peak in the PNP2 spectrum based on QM/MM-MD structures is located at 2.35 eV which has a shoulder on the blue side. The shoulder stems from the same transition as the main peak (which can be seen from Fig. SI2 in Supporting Information). It arises because of the larger spread of the excitation energies and the chosen broadening factor of the Gaussian functions. Besides the main peak (and the shoulder) there is also a low intensity peak centered at approximately 4.3 eV. It is clear, by comparison, that the MM-MD spectrum does not reproduce the QM/MM-MD reference very well. Besides the much lower intensity of the main transition, the energy at maximum 17



absorption is shifted to around 3 eV, and, furthermore, it has shoulders on both sides. The low intensity peak is, however, reasonably reproduced albeit slightly too intense. The QM/MM optimization, on the other hand, results in a spectrum with a clear peak of similar intensity as the main peak in the QM/MM-MD spectrum although it is slightly shifted to higher energy (the energy at maximum absorption is 2.56 eV) and it has no shoulder. The low intensity peak is also shifted a similar amount, which thus indicates that the main part of the shift is most likely due to the use of different functionals. It is clear that the QM/MM optimization removes some of the finite-temperature effects resulting in a smaller spread of the excitation energies, but the main feature is still reasonably reproduced. For PNP3, however, neither the MM-MD nor the QM/MM-optimized structures result in spectra that reproduce the QM/MM-MD reference very well. For PNP3, the QM/MM-MD structures lead to an even wider spread of the main transition than what we saw for PNP2, thus, leading to two peaks with energies 2.11 and 2.80 eV at maximum absorption. Again these peaks arise purely because of the somewhat arbitrary choice of broadening. There is also a low intensity peak at around 4 eV similar to the one in the PNP2 spectrum. The MM-MD structures result in a spectrum where the main transition leads to a single peak located at 2.65 eV, and also a relatively intense peak at 3.90 eV. It is clear that the QM/MM optimization leads to structures that are too ”optimal” thus leading to a single, very intense peak at 2.38 eV. Where the QM/MM geometry optimizations in the case of PNP1 and PNP2 lead to spectra that are closer to the QM/MM-MD reference, this is clearly not the case here. Let us finally comment on the use of PCM for incorporating the effects of the solvent. As shown in SI (Table SI15) for all of the PNPs the mean excitation energy computed using a single geometry optimized structure in combination with PCM is further away from the QM/MM-MD references compared to results based on QM/MM-optimized structures. This is especially the case for PNP3 as expected due to its increased number of torsional degrees of freedom.

18


Page 18 of 32

Page 19 of 32 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


3.3

NMR shielding constants

In addition to the consideration of excitation energies and transition strengths we also analyze the impact of the structural differences between MM-MD, QM/MM-optimized and QM/MMMD in relation to calculation of NMR parameters, e.g. chemical shielding constants. In Table 1 the results from the carbon NMR shielding constant calculations are presented. The RMSD is defined with the results based on the QM/MM-MD structures as the reference, and have been calculated for the averaged shielding constants for all the carbon atoms in each of the three PNPs. The presented standard deviations, σ, represent the average standard deviation for the carbon shielding constants. Table 1: RMSD of the averaged carbon shielding constants calculated using structures based on QM/MM-MD as the reference. In addition the average standard deviation of the calculation carbon NMR shielding constants have been tabulated. Results are in ppm. PNP1 PNP2 PNP3

RMSDMD 9.39 6.59 7.69

RMSDQM/MM 3.53 2.88 4.58

σMD 5.33 5.01 5.28

σQM/MM 2.40 2.78 2.64

σQM/MM−MD 6.45 5.67 6.00

The RMSD detailed in Table 1 supports the trend observed for the excitation energies, i.e. the QM/MM-optimized results are closer to the corresponding results based on the QM/MMMD reference. This is the case for all three PNP molecules. The improvement reflected in the RMSD does, however, not reflect an improvement for every single carbon shielding constants. Some of the shielding constants were found to be slightly worsened by performing the QM/MM geometry optimization. The average carbon NMR shielding constants and their standard deviation can be seen in the supporting information Tables SI11, SI12 and SI13. The average standard deviations of the shielding constants also suggest that the QM/MM geometry optimization removes some the dynamic features of the structures. Unlike for the excitation energies it can be noted that the standard deviations from the QM/MM-MD structure based shielding constants are larger than that of the classical MD based structures which most likely is a result of using constrained hydrogen atoms in the classical MD. 19



4

Conclusion

Classical MD simulations were performed for three different solute solvent systems on the basis of a general MM force field. Based on structural analysis of snapshots extracted from these MM-MD simulations, we investigated to which degree this protocol could reproduce structures comparable with QM/MM-MD simulations. Further it was investigated if QM/MM geometry optimization would improve the quality of the MM-MD structures, using again the QM/MM-MD structures as the reference. When comparing the differences in the absorption spectra based on the MM-MD and QM/MM-MD structures, it was shown that the differences were mainly due the dihedral angle at the nitro groups for PNP1, and due to two different dihedrals for PNP2. We find that QM/MM geometry optimization in general improve the dihedral angles, but at the cost of other structural features suffering a great loss of fluctuation as seen, e.g., for the bond lengths in Figure 4. The shape of the predicted absorption spectra also support the analysis. The spectra based on the MM-MD simulations were flattened out which indicates that a lot of different conformations were present, whereas QM/MM geometry optimization possesses the narrowest spectra. This is also directly reflected from the standard deviation of the excitation energies. For PNP1 the computational cost of QM/MM-MD was about 700 times higher compared to the MD with the QM/MM geometry optimization. For PNP2 and PNP3 the increase in cost was about 1400 and 2300 times, respectively. We also performed calculations of the carbon chemical shielding constants – a property that is well known to be very sensitive to the quality of the underlying molecular structures. The analysis based on the calculation of these shielding constants supports the finding based on the electronic excitation energies and absorption strengths, i.e. the use of QM/MM-optimization generally leads to results that are in better agreement with the ones based on the QM/MM MD but, as for the case of the transition properties, QM/MM-optimization leads to distributions that are too narrow as compared to the corresponding distributions based on QM/MM MD structures.

20


Page 20 of 32

Page 21 of 32 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


In conclusion we find that even though neither the protocol based on MM-MD alone nor the one based on a QM/MM geometry optimization lead to a simultaneous reproduction of the structures and properties based on QM/MM-MD simulations, use of QM/MMoptimized structures for use in embedding calculations represent a cost-efficient alternative to full QM/MM-MD, preserving many of the key structural parameters from QM/MM-MD; however it should be clear that this approach is expected to underestimate e.g. spectral bandwidths. As an alternative, specialized force fields could be constructed but at the cost of repeating such parametrization for each solute molecule.

Acknowledgement Computation/simulation for the work described in this paper was supported by the DeIC National HPC Centre, SDU. J. M. H. O. thanks the Carlsberg Foundation for financial support (grant id. CF15-0823).

Supporting Information Available MD trajectories of all PNPs. Oscillator strengths for PNP2 and PNP3. Oscillator strengths vs. dihedral angles in PNP3. Spectra of all PNPs. Mean excitation and standard deviation for all PNPs. Gaussian convolution equation. Mean shield constants and standard deviations of carbon for all PNPs.

This material is available free of charge via the Internet at

http://pubs.acs.org/.

References (1) Cramer, C. J.; Truhlar, D. G. Implicit Solvation Models: Equilibria, Structure, Spectra, and Dynamics. Chem. Rev. 1999, 99, 2161–2200.

21



(2) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999–3094. (3) Lin, H.; Truhlar, D. G. QM/MM: what have we learned, where are we, and where do we go from here? Theor. Chem. Acc. 2006, 117, 185–199. (4) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem. Int. Ed. 2009, 48, 1198–1229. (5) Coutinho, K.; Canuto, S. Advances in Quantum Chemistry; Elsevier, 1997; pp 89–105. (6) Andreussi, O.; Prandi, I. G.; Campetella, M.; Prampolini, G.; Mennucci, B. Classical Force Fields Tailored for QM Applications: Is It Really a Feasible Strategy? J. Chem. Theory Comput. 2017, 13, 4636–4648. (7) Schwabe, T. Assessing Molecular Dynamics Simulations with Solvatochromism Modeling. J. Phys. Chem. B 2015, 119, 10693–10700. (8) Cerezo, J.; Santoro, F.; Prampolini, G. Comparing classical approaches with empirical or quantum-mechanically derived force fields for the simulation electronic lineshapes: application to coumarin dyes. Theor. Chem. Acc. 2016, 135, 143. (9) van Severen M. C.; Gresh N.; Parisel O.; Cisneros G. A.; Darden T. A.; Piquemal J., C. R. U. S. Progress Towards Accurate Molecular Modeling of Metal Complexes Using Polarizable Force Fields. In Theory and Applications of Computational Chemistry 2008, 185192. (10) Gresh, N.; Claverie, P.; Pullman, A. Theoretical studies of molecular conformation. Derivation of an additive procedure for the computation of intramolecular interaction energies. Comparison withab initio SCF computations. Theor. Chim. Acta 1984, 66, 1–20.

22


Page 22 of 32

Page 23 of 32 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


(11) Piquemal, J.-P.; Chevreau, H.; Gresh, N. Toward a Separate Reproduction of the Contributions to the Hartree-Fock and DFT Intermolecular Interaction Energies by Polarizable Molecular Mechanics with the SIBFA Potential†. J. Chem. Theory Comput. 2007, 3, 824–837. (12) Gresh, N.; Cisneros, G. A.; Darden, T. A.; Piquemal, J.-P. Anisotropic, Polarizable Molecular Mechanics Studies of Inter- and Intramolecular Interactions and LigandMacromolecule Complexes. A Bottom-Up Strategy. J. Chem. Theory Comput. 2007, 3, 1960–1986. (13) Piquemal, J.-P.; Gresh, N.; Giessner-Prettre, C. Improved Formulas for the Calculation of the Electrostatic Contribution to the Intermolecular Interaction Energy from Multipolar Expansion of the Electronic Distribution. J. Phys. Chem. A 2003, 107, 10353–10359. (14) Piquemal, J.-P.; Cisneros, G. Many-Body Effects and Electrostatics in Biomolecules; Pan Stanford, 2016; pp 269–299. (15) Piquemal, J.-P.; Cisneros, G. A.; Reinhardt, P.; Gresh, N.; Darden, T. A. Towards a force field based on density fitting. J. Chem. Phys. 2006, 124, 104101. (16) Duke, R. E.; Starovoytov, O. N.; Piquemal, J.-P.; Cisneros, G. A. GEM: A Molecular Electronic Density-Based Force Field for Molecular Dynamics Simulations. J. Chem. Theory Comput. 2014, 10, 1361–1365. (17) Cisneros, G. A.; Piquemal, J.-P.; Darden, T. A. Quantum Mechanics/Molecular Mechanics Electrostatic Embedding with Continuous and Discrete Functions. J. Phys. Chem. B 2006, 110, 13682–13684. (18) Gao, J. Toward a Molecular Orbital Derived Empirical Potential for Liquid Simulations. J. Phys. Chem. B 1997, 101, 657–663.

23



(19) Xie, W.; Orozco, M.; Truhlar, D. G.; Gao, J. X-Pol Potential: An Electronic StructureBased Force Field for Molecular Dynamics Simulation of a Solvated Protein in Water. J. Chem. Theory Comput. 2009, 5, 459–467. (20) Xie, W.; Gao, J. Design of a Next Generation Force Field: The X-POL Potential. J. Chem. Theory Comput. 2007, 3, 1890–1900. (21) Song, L.; Han, J.; lin Lin, Y.; Xie, W.; Gao, J. Explicit Polarization (X-Pol) Potential Using ab Initio Molecular Orbital Theory and Density Functional Theory†. J. Phys. Chem. A 2009, 113, 11656–11664. (22) Han, J.; Truhlar, D. G.; Gao, J. Optimization of the explicit polarization (X-Pol) potential using a hybrid density functional. Theor. Chem. Acc. 2012, 131 . (23) Wang, Y.; Sosa, C. P.; Cembran, A.; Truhlar, D. G.; Gao, J. Multilevel X-Pol: A Fragment-Based Method with Mixed Quantum Mechanical Representations of Different Fragments. J. Phys. Chem. B 2012, 116, 6781–6788. (24) Han, J.; Mazack, M. J. M.; Zhang, P.; Truhlar, D. G.; Gao, J. Quantum mechanical force field for water with explicit electronic polarization. J. Chem. Phys. 2013, 139, 054503. (25) Giese, T. J.; Huang, M.; Chen, H.; York, D. M. Recent Advances toward a General Purpose Linear-Scaling Quantum Force Field. Acc. Chem. Res. 2014, 47, 2812–2820. (26) Giese, T. J.; Chen, H.; Huang, M.; York, D. M. Parametrization of an Orbital-Based Linear-Scaling Quantum Force Field for Noncovalent Interactions. J. Chem. Theory Comput. 2014, 10, 1086–1098. (27) Giese, T. J.; Panteva, M. T.; Chen, H.; York, D. M. Multipolar Ewald Methods, 2: Applications Using a Quantum Mechanical Force Field. J. Chem. Theory Comput. 2015, 11, 451–461. 24


Page 24 of 32

Page 25 of 32 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


(28) Giese, T.; York, D. Many-Body Effects and Electrostatics in Biomolecules; Pan Stanford, 2016; pp 3–31. (29) Giese, T. J.; Chen, H.; Dissanayake, T.; Giamba¸su, G. M.; Heldenbrand, H.; Huang, M.; Kuechler, E. R.; Lee, T.-S.; Panteva, M. T.; Radak, B. K.; York, D. M. A Variational Linear-Scaling Framework to Build Practical, Efficient Next-Generation Orbital-Based Quantum Force Fields. J. Chem. Theory Comput. 2013, 9, 1417–1427. (30) Giese, T. J.; York, D. M. Quantum mechanical force fields for condensed phase molecular simulations. J. Phys. Condens. Matter 2017, 29, 383002. (31) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227–249. (32) Brunk, E.; Rothlisberger, U. Mixed Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulations of Biological Systems in Ground and Electronically Excited States. Chem. Rev. 2015, 115, 6217–6263. (33) Olsen, J. M.; Aidas, K.; Kongsted, J. Excited states in solution through polarizable embedding. J. Chem. Theory Comput. 2010, 6, 3721–3734. (34) Giese, T. J.; York, D. M. A GPU-Accelerated Parameter Interpolation Thermodynamic Integration Free Energy Method. J. Chem. Theory Comput. 2018, 14, 1564–1582. (35) Kongsted, J.; Ryde, U. An improved method to predict the entropy term with the MM/PBSA approach. J. Comput.-Aided Mol. Des. 2008, 23, 63–71. (36) Zuehlsdorff, T. J.; Isborn, C. M. Combining the ensemble and Franck-Condon approaches for calculating spectral shapes of molecules in solution. J. Chem. Phys. 2018, 148, 024110.

25



(37) Petersson, G. A.; Bennett, A.; Tensfeldt, T. G.; Al-Laham, M. A.; Shirley, W. A.; Mantzaris, J. A complete basis set model chemistry. I. The total energies of closed-shell atoms and hydrides of the first-row elements. J. Chem. Phys. 1988, 89, 2193–2218. (38) Petersson, G. A.; Al-Laham, M. A. A complete basis set model chemistry. II. Openshell systems and the total energies of the first-row atoms. J. Chem. Phys. 1991, 94, 6081–6090. (39) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian09 Revision E.01. Gaussian Inc. Wallingford CT 2009. (40) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 1993, 97, 10269–10280. (41) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graphics Modell. 2006, 25, 247–260. 26


Page 26 of 32

Page 27 of 32 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


(42) 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, 1157–1174. (43) Jorgensen, W. L.; Madura, J. D. Quantum and statistical mechanical studies of liquids. 25. Solvation and conformation of methanol in water. J. Am. Chem. Soc. 1983, 105, 1407–1413. (44) Case, D. A.; Babin, V.; Berryman, J. T.; Betz, R. M.; Cai, Q.; Cerutti, D. S.; Cheatham, T. E.; Darden, T. A.; Duke, R. E.; Gohlke, H.; Goetz, A. W.; Gusarov, S.; Homeyer, N.; Janowski, P.; Kaus, J.; Kolossváry, I.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Luchko, T.; Luo, R.; Madej, B.; Merz, K. M.; Paesani, F.; Roe, D. R.; Roitberg, A.; Sagui, C.; Salomon-Ferrer, R.; Seabra, G.; Simmerling, C. L.; Smith, W.; Swails, J.; Walker,; Wang, J.; Wolf, R. M.; Wu, X.; Kollman, P. A. Amber 14. University of California, San Francisco, 2014. (45) Andersen, H. C. Molecular dynamics simulations at constant pressure and/or temperature. J. Chem. Phys. 1980, 72, 2384–2393. (46) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684– 3690. (47) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327–341. (48) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089–10092. (49) Roe, D. R.; Cheatham, T. E. PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theory Comput. 2013, 9, 3084–3095. 27



(50) Murphy, R. B.; Philipp, D. M.; Friesner, R. A. A mixed quantum mechanics/molecular mechanics (QM/MM) method for large-scale modeling of chemistry in protein environments. J. Comput. Chem. 2000, 21, 1442–1457. (51) Philipp, D. M.; Friesner, R. A. Mixed ab initio QM/MM modeling using frozen orbitals and tests with alanine dipeptide and tetrapeptide. J. Comput. Chem. 1999, 20, 1468– 1494. (52) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098–3100. (53) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785–789. (54) Banks, J. L.; Beard, H. S.; Cao, Y.; Cho, A. E.; Damm, W.; Farid, R.; Felts, A. K.; Halgren, T. A.; Mainz, D. T.; Maple, J. R.; Murphy, R.; Philipp, D. M.; Repasky, M. P.; Zhang, L. Y.; Berne, B. J.; Friesner, R. A.; Gallicchio, E.; Levy, R. M. Integrated Modeling Program, Applied Chemical Theory (IMPACT). J. Comput. Chem. 2005, 26, 1752–1780. (55) Marx, D.; Hutter, J. Ab initio molecular dynamics: Theory and implementation. Modern methods and algorithms of quantum chemistry 2000, 1, 141. (56) Andreoni, W.; Curioni, A. New Advances in Chemistry and Material Science with CPMD and Parallel Computing. Parallel Comput. 2000, 26, 819–842. (57) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics; Cambridge University Press, 2009; http://www.cambridge.org/gb/knowledge/isbn/item2327682/. (58) Christen, M.; H¨ unenberger, P. H.; Bakowies, D.; Baron, R.; B¨ urgi, R.; Geerke, D. P.; Heinz, T. N.; Kastenholz, M. A.; Kräutler, V.; Oostenbrink, C.; Peter, C.; Trzes-

28


Page 28 of 32

Page 29 of 32 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


niak, D.; van Gunsteren, W. F. The GROMOS software for biomolecular simulation: GROMOS05. J. Comput. Chem. 2005, 26, 1719–1751. (59) Troullier, N.; Martins, J. L. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B 1991, 43, 1993–2006. (60) Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 52, 255–268. (61) Hoover, W. G. Canonical dynamics: equilibrium phase-space distributions. Phys. Rev. A 1985, 31, 1695. (62) Hockney, R. Potential calculations and some applications. Methods Comput. Phys. 9: 135-211(1970). 1970, (63) Laio, A.; VandeVondele, J.; Rothlisberger, U. A Hamiltonian electrostatic coupling scheme for hybrid Car–Parrinello molecular dynamics simulations. J. Chem. Phys. 2002, 116, 6941–6947. (64) Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange–correlation functional using the Coulomb-attenuating method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51–57. (65) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; Dalskov, E. K.; Ekström, U.; Enevoldsen, T.; Eriksen, J. J.; Ettenhuber, P.; Fernández, B.; Ferrighi, L.; Fliegl, H.; Frediani, L.; Hald, K.; Halkier, A.; Hättig, C.; Heiberg, H.; Helgaker, T.; Hennum, A. C.; Hettema, H.; Hjertenæs, E.; Høst, S.; Høyvik, I.-M.; Iozzi, M. F.; Jans´ık, B.; Jensen, H. J. Aa.; Jonsson, D.; Jørgensen, P.; Kauczor, J.; Kirpekar, S.; Kjærgaard, T.; Klopper, W.; Knecht, S.; Kobayashi, R.; Koch, H.; Kongsted, J.; Krapp, A.; Kristensen, K.; Ligabue, A.; Lutnæs, O. B.; Melo, J. I.; Mikkelsen, K. V.; Myhre, R. H.; Neiss, C.;

29



Page 30 of 32

Nielsen, C. B.; Norman, P.; Olsen, J.; Olsen, J. M. H.; Osted, A.; Packer, M. J.; Pawlowski, F.; Pedersen, T. B.; Provasi, P. F.; Reine, S.; Rinkevicius, Z.; Ruden, T. A.; Ruud, K.; Rybkin, V. V.; Salek, P.; Samson, C. C. M.; de Merás, A. S.; Saue, T.; Sauer, S. P. A.; Schimmelpfennig, B.; Sneskov, K.; Steindal, A. H.; SylvesterHvid, K. O.; Taylor, P. R.; Teale, A. M.; Tellgren, E. I.; Tew, D. P.; Thorvaldsen, A. J.; Thøgersen, L.; Vahtras, O.; Watson, M. A.; Wilson, D. J. D.; Ziolkowski, M.; ˚ Agren, H. The Dalton quantum chemistry program system. WIREs Comput. Mol. Sci. 2014, 4, 269–284. (66) Dalton, a molecular electronic structure program, Release Dalton2017.alpha (2016). ”see http://daltonprogram.org/”, commit 8f11cac. (67) Olsen, J. M. H.; List, N. H.; Steinmann, C.; Steindal, A. H.; Nørby, M. S. PElib:

The

Polarizable

Embedding

library

(version

1.3.5).

2018;

https://doi.org/10.5281/zenodo.1209644. (68) Gao, B.; Thorvaldsen, A. J.; Ruud, K. GEN1INT: A unified procedure for the evaluation of one-electron integrals over Gaussian basis functions and their geometric derivatives. Int. J. Quantum Chem. 2011, 111, 858–872. (69) Gao, B. Gen1Int (version 0.2.1). 2012; http://gitlab.com/bingao/gen1int. (70) Zale´sny, R.; Murugan, N. A.; Gel’mukhanov, F.; Rinkevicius, Z.; O´smialowski, B.; Bartkowiak, W.; ˚ Agren, H. Toward Fully Nonempirical Simulations of Optical Band Shapes of Molecules in Solution: A Case Study of Heterocyclic Ketoimine Difluoroborates. J. Phys. Chem. A 2014, 119, 5145–5152. (71) Cerezo, J.; Ferrer, F. J. A.; Prampolini, G.; Santoro, F. Modeling Solvent Broadening on the Vibronic Spectra of a Series of Coumarin Dyes. From Implicit to Explicit Solvent Models. J. Chem. Theory Comput. 2015, 11, 5810–5825.

30


Page 31 of 32 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


(72) Keal, T. W.; Tozer, D. J. A semiempirical generalized gradient approximation exchangecorrelation functional. J. Chem. Phys. 2004, 121, 5654–5660. (73) Jensen, F. Polarization consistent basis sets: Principles. J. Chem. Phys. 2001, 115, 9113–9125. (74) Jensen, F.; Helgaker, T. Polarization consistent basis sets. V. The elements Si–Cl. J. Chem. Phys. 2004, 121, 3463–3470.

31



Graphical TOC Entry

32


Page 32 of 32

Importance of Accurate Structures for Quantum Chemistry Embedding

Recommend Documents