Interfacing CRYSTAL/AMBER to Optimize QM/MM Lennard-Jones Parameters for Water and to Study Solvation of TiO2 Nanoparticles
Asmus Ougaard Doh, Daniele Selli, Gianluca Fazio, Lorenzo Ferraro, Jens Jørgen Mortensen, Bartolomeo Civalleri, Cristiana Di Valentin
AArticle
Interfacing CRYSTAL / AMBER to Optimize QM / MMLennard–Jones Parameters for Water and to StudySolvation of TiO Nanoparticles
Asmus Ougaard Dohn * , Daniele Selli , Gianluca Fazio , Lorenzo Ferraro ,Jens Jørgen Mortensen , Bartolomeo Civalleri and Cristiana Di Valentin Faculty of Physical Sciences and Science Institute, University of Iceland, 107 Reykjavík, Iceland Dipartimento di Scienza dei Materiali, Università di Milano-Bicocca, via Cozzi 55, 20125 Milano, Italy;[email protected] (D.S.); [email protected] (G.F.); [email protected] (L.F.);[email protected] (C.D.V.) CAMD, Department of Physics, Technical University of Denmark, 2800 Kongens Lyngby, Denmark;[email protected] Dipartimento di Chimica, Università di Torino and NIS Centre of Excellence, Via P. Giuria 7, I-10129 Torino,Italy; [email protected] * Correspondence: [email protected]; Tel.: + // / journal / molecules / special_issues / QM------------------------------------------------------------------------------------------------------------------------------------------------Received: 1 October 2018; Accepted: 6 November 2018; Published: 13 November 2018 !" ! !" Abstract:
Metal oxide nanoparticles (NPs) are regarded as good candidates for many technologicalapplications, where their functional environment is often an aqueous solution. The correct descriptionof metal oxide electronic structure is still a challenge for local and semilocal density functionals,whereas hybrid functional methods provide an improved description, and local atomic function-basedcodes such as CRYSTAL17 outperform plane wave codes when it comes to hybrid functionalcalculations. However, the computational cost of hybrids are still prohibitive for systems of real sizes,in a real environment. Therefore, we here present and critically assess the accuracy of our electrostaticembedding quantum mechanical / molecular mechanical (QM / MM) coupling between CRYSTAL17and AMBER16, and demonstrate some of its capabilities via the case study of TiO NPs in water.First, we produced new Lennard–Jones (LJ) parameters that improve the accuracy of water–waterinteractions in the B3LYP / TIP3P coupling. We found that optimizing LJ parameters based on watertri- to deca-mer clusters provides a less overstructured QM / MM liquid water description than whenfitting LJ parameters only based on the water dimer. Then, we applied our QM / MM couplingmethodology to describe the interaction of a 1 nm wide multilayer of water surrounding a sphericalTiO nanoparticle (NP). Optimizing the QM / MM water–water parameters was found to have little tono e ↵ ect on the local NP properties, which provide insights into the range of influence that can beattributed to the LJ term in the QM / MM coupling. The e ↵ ect of adding additional water in an MMfashion on the geometry optimized nanoparticle structure is small, but more evident e ↵ ects are seenin its electronic properties. We also show that there is good transferability of existing QM / MM LJparameters for organic molecules–water interactions to our QM / MM implementation, even thoughthese parameters were obtained with a di ↵ erent QM code and QM / MM implementation, but with thesame functional.
Keywords: QM / MM; multiscale; nanoparticles; force field parameters; water; titanium dioxide;geometry optimization; molecular dynamics
Molecules , , 2958; doi:10.3390 / / journal / molecules olecules , , 2958 2 of 24
1. Introduction
After the pioneering work of M. Levitt and A. Warshel [1], the fundamental QM / MM methodologyof joining explicit electronic structure to classical potential functions is constantly being both appliedand expanded upon [2–7]. While much work is being put into methods that go beyond the electrostaticembedding-type of additive QM / MM, which explicitly couples the two subsystems via a Coloumbterm but leaves the MM charges fixed, this form of coupling still seems prevalent within additiveQM / MM applications [8], most likely due to its generally accepted trade-o ↵ between accuracy ande ciency. Some of the authors involved in this work have previously demonstrated how electrostaticembedding QM / MM can be used to predict solvation dynamics of bond formation in a model catalystsystem, which was later confirmed by experiment [9,10]. However, the QM / MM implementationemployed previously [11] relies on a DFT code that at the time of writing does not yet support thecalculation of forces from hybrid functionals. In contrast, the MPP-CRYSTAL17 [12–14] code has ane cient implementation of the exact exchange term [12]. This capability is crucial for successfullymodeling semiconducting oxide systems, such as TiO nanoparticles (NPs) in aqueous solution [15].Metal oxide nanoparticles NPs are widely used in many technological applications, and especially TiO nanomaterials are very popular because of their abundance, nontoxicity, high stability under a varietyof conditions, and biocompatibility. Their uses ranges from traditional ones, such as in the pigmentindustry or in cosmetics, to more advanced fields, such as photocatalysis, photoelectrochemistry,fuel production, and nanomedicine [16–19]. Many, if not all, of the applications take place in aqueousenvironments, where water is often a main actor in the physics and chemistry of such processes [20,21].However, the modeling of TiO nanoparticles is computationally very demanding [22–25]: apartfrom requiring hybrid functionals, realistic models consist of 500–1000 atoms for NPs of at least 2 nmof diameter, which cannot be treated periodically. Therefore, such simulations require the use ofparallelized, e ciently performing codes, such as MPP-CRYSTAL17. However, when also includingthe aqueous environment of the NP (easily around 3000 atoms), current computational power still doesnot allow for a full, exact-exchange-including DFT treatment. Thus, we here report on the developmentand benchmarking of an electrostatic embedding QM / MM framework that allows for simulations ofsystems of these sizes and geometries. The highly arched curvature of spherical NP surface is knownto partially dissociate the first water monolayer around it [26,27]. In addition, the curved surface haslong-range physical e ↵ ects on surrounding water, as recently reported [28]. Thus, when preparingsimulations of NPs, it is important to assess how much water should be included in the simulation toachieve the desired realism. To describe the aforementioned dissociation of the first monolayer, it isnecessary to include it in the QM subsystem in the multiscale divisioning, when working with simple,rigid, point charge-based force fields. This means that there will almost exclusively be water close tothe QM / MM intersection, which again motivates a thorough analysis / optimization of the water–waterQM / MM coupling accuracy.The paper is structured as follows: In Section 2, we present the results from the water–waterQM / MM LJ re-parameterization and implementation (Section 2.1), before turning to solvation of aTiO NP (Section 2.2). In Section 3.1, we analyse the consequences of choosing various LJ fittingstrategies, assess the transferability of QM / MM LJ parameters between di ↵ erent electronic structurecodes, and the extent of the influence of water–water LJ QM / MM interactions. Section 3.2 discusses theconsequences of including more water in the NP simulations. Finally, the computational details of thevarious sections are presented in Section 4.2, before Section 5 wraps up the conclusions from our work.The appendices contain further benchmarks on our QM / MM implementation.The updated QM / MM implementation and CRYSTAL17 / AMBER16 interface are available in theo cial version of the Atomic Simulation Environment (ASE) [29,30] (see https: // wiki.fysik.dtu.dk / ase / ). olecules , , 2958 3 of 24
2. Results / MM Water–Water Interactions
By comparing QM / MM interaction energies with their pure QM counterparts on a system thatonly contains one molecular species, we can assess the accuracy of our novel electrostatic embeddingframework [11,31]. Since the method is focused on describing molecules and larger systems suchas nanoparticles in aqueous solution, we start by focusing on neat water. Some of the main criteriaof success for achieving an accurate QM / MM method are that: (1) the QM / MM coupled interactionenergies do not over- or under-bind when compared to neither the chosen QM nor the MM models:and (2) the various possible combinations of QM / MM geometries and regions are as similar as possible,so that there will be no orientation-induced artifacts in the total energies and forces [2]. The QM / MMcoupling model used in electrostatic embedding (see Section 4.1) leaves open two main routes for thetuning of the accuracy of the coupling: (1) through how the additional electrostatic term in the externalpotential in the DFT code is constructed [11,32]; and (2) via a re-parameterization of the Lennard–Jonespotential that describes non-electrostatic interactions between the atoms in the QM and MM subsystem(Equation (4) in Section 4.1). Here, we focus on the second option. We have employed two strategiesfor optimizing QM / MM LJ parameters: The
Dimer Fit (DF) is based on optimizing the LJ parametersto minimize the di ↵ erence between the multiscale interaction energy curves of the water dimer andthe average of the pure QM and pure MM interaction energy curves, a methodology akin to previouswork on other systems [33]. The other strategy attempts to take into account that the structure of liquidwater cannot be completely understood only analyzing the water dimer, and thus, the Cluster Fit (CF)is a global fit on interaction energies on water clusters of size 3–10 molecules per cluster. The details ofthe fits can be found in Section 4.2.2, which the very interested reader might benefit from reading first,before continuing to the results section. Apart from minimizing the errors in structure and energyintroduced by partitioning the system, the fitting e ↵ ort also leads to insights into the transferability ofthe molecular characteristics of the water dimer into the modeling of liquid water.Hybrid functionals provide a better description of the water liquid structure, compared to theirGGA relatives [34–36], although often semi-empirical corrections to GGA functionals are employed [37].In contrast, for the MM subsystem, there are many known improvements to the TIP3P descriptionof water, either within the same general model of static point charges [38], or even extending thedescription to account for the polarizability of water [39–41], and / or flexible molecules [41,42]. However,employing more advanced MM water potentials is outside the scope of this work because: (1) TIP3Pis still ubiquitous in many QM / MM simulations, especially of very large total systems, where theadded computational expense of going beyond pairwise potentials can become prohibitive. (2) Whenincluding the first solvation shell (or more) in the QM region and focusing on explicit solvation e ↵ ectson the solute in particular, and not so much the solute–solvent interactions themselves, the outer layersof water become less important, thus making TIP3P an adequate choice.2.1.1. Water DimerThroughout this section, we use the following terminology: “QM / MM” means that thehydrogen-donating molecule is in the QM region, while the hydrogen-accepting is described withTIP3P, and vice versa for “MM / QM”. We address both configurations collectively as “the multiscaleconfigurations”, to avoid confusion.Figure 1 reports the structure of the water dimer considered and relevant quantities assessed forour QM / MM implementation.The hydrogen bonding energy of the water dimer is calculated with the three di ↵ erent sets of LJparameters, and compared to the pure B3LYP and pure TIP3P results in Figure 2. Using the TIP3PLJ parameters with our choice of functional and basis set within CRYSTAL17 produces overboundmultiscale dimer binding curves for both configurations, but most prominently for the QM / MMconfiguration, where the electronic density of the hydrogen-donating molecule is modeled explicitly. olecules , , 2958 4 of 24 At very short distances, dominated by Pauli repulsion, the di ↵ erences in binding energies betweenthe multiscale and the B3LYP curves are large, indicating that a re-evaluation of the LJ parametersare warranted. For the multiscale binding curves using the DF LJ parameters (see Section 4.2.2 fordetails on the fitting strategies), the QM / MM curve becomes almost identical to the TIP3P result,even though the fit target is the average of the QM and the MM binding curves. The resulting MM / QMcurve is also brought closer to its QM counterpart, and all overbinding has been removed from themodel. The CF strategy seemingly produces binding curves that are still overbound around theirminima, but it interestingly decreases the di ↵ erence between the O-O distances with the two multiscaleminimum values, e ↵ ectively making the multiscale model more consistent with respect to geometricdi ↵ erences. See Appendix A for an analysis of the various relaxed multiscale geometries for all possibleconfigurations of the water dimer. (cid:436)(cid:437) (cid:443) r OO H D O A O D H A H-DonatingMolecule
QM/MM: H-Donating is QMMM/QM: H-Accepting is QM
H-AcceptingMolecule
Figure 1.
Illustration of the geometric quantities used in the water benchmarks. For the dimer testsand fits, we use the shorthand nomenclatures “QM / MM” and “MM / QM” to discern between the twopossible subsystem configurations. Three angles provide insight into the hydrogen bonding networkstructure of liquid water, and are sampled in the molecular dynamics simulations: The donor angles ↵ = \ ( O D , H D , O A ) and = \ ( H D , O D , O A ) , and the acceptor angle ✓ = \ ( H D , O A , H A ) . E n e r g y ( k c a l / m o l ) B3LYPTIP3PB3LYP/TIP3P LJ: TIP3PTIP3P/B3LYP LJ: TIP3P E n e r g y ( k c a l / m o l ) B3LYPTIP3PB3LYP/TIP3P LJ: DFTIP3P/B3LYP LJ: DF
O-O Distance (˚A) E n e r g y ( k c a l / m o l ) B3LYPTIP3PB3LYP/TIP3P LJ: CFTIP3P/B3LYP LJ: CF
Figure 2.
Water dimer binding curves (counterpoise corrected), calculated using LJ parameters from:TIP3P ( top ); the dimer fit ( middle ); and the cluster fit ( bottom ). The starting geometry is obtainedfrom the S22 database, where the water dimer was optimized using CCSD(T) [43], and all otherintramolecular geometries are kept constant while scanning the O-O distance. olecules , , 2958 5 of 24 ↵ erence in interaction QM / MM energies and their pure B3LYP counterparts: DD E = D E QMMMint D E B3LYPint n (1)evaluated over a large dataset of water clusters ranging from n = / MM LJ parameters for water–water interactions.Figure 3 shows box plots representing the distributions in interaction energy di ↵ erences with respectto the pure B3LYP results. The blue patches represent the maximum di ↵ erence in interaction energiesfrom the pure B3LYP and pure TIP3P results. Using the TIP3P LJ parameters results in parts of theQM / MM DD E -values falling below the blue patches, meaning that there is a significant overbindingwith respect to the single-description results. We also observe a trend of the binding systematicallybeing tightest for the multiscale configurations that have the most QM / MM interactions, resulting in“u-shaped” curves. When using the DF LJ parameters, this u-shape has disappeared, but, especially forthe 8- and 9-mer clusters, the QM / MM interaction energies give significantly underbound clusters,as the DD E -values end up over the upper limit of the blue patches. The CF parameters retain theu-shape, indicative of the same systematic increase in binding with number of QM / MM LJ terms,but all QM / MM distributions are now within the total di ↵ erence spanned by the limits of the two puredescriptions. In all cases, spread of energy di ↵ erences is largest when the multiscale description isfurthest removed from the reference (i.e., 1 QM water vs. pure QM water), as one would expect for awell-behaved interface. The fact that the TIP3P and CF results are more similar is most likely becausethe TIP3P parameters are already optimized to reproduce bulk phase properties [2], which is discussedfurther in Section 3.1. Number Of QM Water Molecules E ( k c a l / m o l ) Underboundw.r.t B3LYPOverbound w.r.t B3LYP max( E TIP3P E B3LYP ) / n LJ: TIP3P LJ: DF LJ: CF
Figure 3.
Boxplot of the di ↵ erence in interaction energies with respect to the pure B3LYP results.The horizontal axis is firstly divided divided into eight subsets based on the total number of moleculesin each size of water cluster, and then further divided into subsets defined by how many of the totalamount of molecules are in the QM subsystem. Each box represents the distribution of interaction energydi ↵ erences from all possible QM / MM combinations with that number of QM molecules. The boxesspans 50% (the IQR) of the distribution of energy di ↵ erences, while the whiskers represent ± ↵ erent sets of LJ parameters: The original TIP3Pparameters (blue), the dimer fit (DF) parameters (red), and the cluster fit (CF) parameters (green).The blue patches behind the box plots show the maximum di ↵ erence between the two pure descriptions. olecules , , 2958 6 of 24 As mentioned in Section 1, it is important to explore the e ↵ ects of including increasing amounts ofwater around NPs, since the nanoparticle / water interfacial e ↵ ects are quite long range. Here, we havestudied three systems: (A) the bare NP (TiO ) · O (Figure 4a); (B) the NP with a first watermonolayer (ML) adsorbed on the surface containing 134 molecules, ⇠
20% of which are dissociated,as discussed and detailed in Section 4.2.5 and in [28] (Figure 4b); and (C), comprised of systemB, with a molecular mechanic (MM) region composed by 824 surrounding waters added around it(Figure 4c). The figure shows the structures after geometry optimizations. Systems A and B wereoriginally prepared for a previous study by some of us [28]. The geometry optimization for system Cwas carried out twice, once with the water in the MM region described with the CF LJ parameters,and once with the TIP3P LJ parameters, to assess any possible e ↵ ects on the NP. Figure 4.
Optimized structure with B3LYP: ( a ) for the bare TiO nanoparticle; ( b ) for the nanoparticlewith an adsorbed water monolayer composed by 134 molecules ( ⇠
20% of which are dissociated);and ( c ) with QM / MM for the nanoparticle with a water multilayer of 958 molecules, using CF LJparameters. Titanium, oxygen and hydrogen atoms of the NP are in cyan, red and white, respectively.Oxygen atoms of molecular water are in blue, while oxygen and hydrogen atoms from the dissociatedwater molecules are in green and yellow, respectively. olecules , , 2958 7 of 24 First, the e ↵ ect of water around the NP has been investigated with respect to the induced structuralchanges in the optimized geometries. In Figure 5, the simulated EXAFS spectra of the bare NP, of theNP with an adsorbed water monolayer and of the one with a water multilayer are compared withthat of bulk anatase. The distribution is quite broad for the bare NP, while it is reduced upon waterML adsorption. The change is most drastic going from no water at all to some water, but, when themultilayer is added, there is only a further tendency of the distances, in particular of Ti-O ax , to be centredat the bulk value, as the grow-in of a new peak at this distance indicates. This is true for calculationsperformed with both the CF and TIP3P LJ parameters, which produce almost identical spectra. Figure 5.
Distances-distribution (simulated EXAFS) computed with B3LYP for the bare NP (black),the NP with the adsorbed water monolayer (ML, blue), computed with QM / MM for the NP with thewater multilayer using the CF LJ (multiL, red) and using the TIP3P LJ (multiL, dashed green). Dashedblack lines correspond to the anatase bulk Ti—O, Ti ---
Ti and T ---
O distances.
The e ↵ ect on the NP electronic properties of the QM / MM relaxation has also been investigated.Figure 6 compares the total (DOS) and projected density of states (PDOS) on oxygen species forthe bare NP, the NP with a water monolayer and with a water multilayer around it. The electronicstructure of the NP in vacuum has been evaluated both with PBE and B3LYP functionals (Figure 6a,b,respectively). First, when switching from a GGA functional such as PBE to a functional including exactexchange, we observe an expected widening of the bandgap from 2.44 to 4.14 eV [15]. When a ML ofwater is adsorbed on the NP surface (Figure 6c), we observe a 0.1 eV change in the band gap value,and a decrease of the work function by + ↵ ect is quite reduced, with a shift of + + olecules , , 2958 8 of 24 Figure 6.
Total (DOS) and projected density (PDOS) of states on di ↵ erent oxygen species of: ( a ) theNP in vacuum calculated with the PBE functional; ( b ) the NP in vacuum calculated with the B3LYPfunctional; ( c ) the NP with a QM water monolayer; ( d ) with a QM / MM water multilayer and the CFLJ parameters; and ( e ) with a QM / MM water multilayer and the TIP3P LJ parameters. In each panel,the band gap and the energy shift of the valence and conduction band edges with respect to the vacuumNP model in a vacuum are given in eV. The zero energy is set to the vacuum level. Panels ( b , c ) are fromprevious work by some of the authors [28].
3. Discussion / MM LJ Re-Parameterization and Parameter Transferability Tests
Table 1 compares the values of OO and ✏ OO , used in the non-electrostatic term of the QM / MMcoupling (Equation (4) in Section 4.1) obtained from our two di ↵ erent methods of fitting water–waterQM / MM LJ parameters, and shows that optimizing LJ parameters only on the water dimer givesrise to a pairwise potential that is very di ↵ erent from both the TIP3P and CF potentials: while thevan der Waals (vdW) radius is increased significantly, the potential well depth is over 10 times moreshallow. On the contrary, the CF methodology reaches its optimum values by reducing the vdWradius and deepening the potential. The cluster- (Figure 3) and liquid-water benchmarks (Figures A3 olecules , , 2958 9 of 24 and A4 in Appendix C) present a good case for the necessity of reaching beyond dimer systems whengenerating QM / MM parameters for solvents. This might not be particularly surprising, since the fullthree-dimensional hydrogen bonding network of water cannot be realized by molecular systems smallerthan the water hexamer. Furthermore, the MM parameters are usually optimized to best describecondensed-phase systems that behave systematically more polarized, which again means that MM gasphase potentials might have overemphasized interactions [2]. The description of polarization in theQM system is more transferable between phases, therefore it should not have as overemphasized gasphase interactions as the MM description. Thus, when fitting the QM / MM interaction (partially) to theMM interaction (and partially to the QM interaction which is not correspondingly under-interacting),the over-emphasis carries over into the obtained parameters.
Table 1.
The two sets of QM / MM LJ parameters generated in this work, compared to the TIP3P LJparameters. See Section 4.1 for their definition.
Type OO ( Å ) ✏ OO (kcal / mol) TIP3P 3.15061 0.1521Dimer Fit (DF) 3.89048 0.0122Cluster Fit (CF) 3.10031 0.2629
Remembering that the DF method with its order of magnitude lower ✏ -value is just short ofremoving the attractive part of the LJ potential altogether, it is interesting to briefly revisit Figure 3.The two datasets with normal ✏ -values are seen to systematically bind hardest when the number ofQM / MM LJ interactions is largest, giving these u-shaped collections of box plots. Since this featureis absent from the DF dataset, we can attribute this type of QM / MM overbinding (within the limitsspanned by the pure QM pure MM di ↵ erence) to the LJ potential. Since performing a fit on the entiredataset of clusters does not remove this behaviour, it must be inherent in the actual formulation of theinteraction, and thus, it is necessary to improve on that, or go beyond electrostatic embedding QM / MMaltogether, if one wishes to increase the accuracy even further.It should be noted that when it comes to QM / MM methodologies that involve fitting of one sortor another, many much more “modern” strategies involving machine learning routines have beenformulated [6], but how to apply such strategies to explicitly interacting QM and MM subsystems(e.g., electrostatic embedding QM / MM) seems not yet quite realized. The fitting methodology presentedhere is very basic in comparison, but has the advantage of only requiring re-evaluations of the pairwiseLJ potential, which means that, if one were to replace the MM force field, no new QM simulationswould have to be made to redo the fit, and, when the new parameters are obtained, there is no extracost to the following simulations.In the same vein as our water–water LJ re-parameterization, we have also tested a previousQM / MM re-parameterization of LJ parameters of organic molecule / water dimers to assess the moregeneral accuracy of our novel CRYSTAL17 / AMBER16 coupling, the details of which can be found inAppendix B. The QM / MM LJ parameters for this section were obtained from the work of Freindorfet al. When transferring the QM / MM LJ parameters optimized by Freindorf et al. to CRYSTAL17,the RMSD values increase from the results in the paper by 0.11 kcal / mol, 0.01 Å, and 0.98 degrees,for the energy, hydrogen bond distances and angles, respectively. Taking all possible calculationaldi ↵ erences mentioned in Appendix B into account, we conclude that simply transferring the previouslyproduced QM / MM parameters to a new QM / MM implementation can be done without a dramatic lossof accuracy. Lastly, we also note that we find no systematicity in any of the errors based on whetherthe role of the water molecule in the dimers is to act as a hydrogen-donor or -acceptor, which againindicates a robust QM / MM description, which is not significantly and / or systematically sensitive to thenature of the hydrogen bonding geometry over the QM / MM division. olecules , , 2958 10 of 24 For the simulated EXAFS spectra of the three NP systems (Figure 5), we saw that, while the mostdrastic changes happen when going from naked to ML hydrated NP, there were only fine di ↵ erencesbetween the QM monolayer and the QM / MM multilayer simulations, in particular associated tothe emergence of a peak at the bulk anatase Ti–O ax distance in the QM / MM multilayer simulation.The e ↵ ects of the QM / MM LJ re-parameterization on the NP structure were negligible, indicating thatthe inclusion of the first layer of water in the QM subsystem makes the NP structural properties robustto changes in the water–water QM / MM LJ coupling terms. This conclusion is thus also an input to theongoing debate about the impact of the somewhat ad hoc treatment of non-electrostatic interactionsover the QM / MM border. This study seems to show that, as long as this border is far enough removedfrom the main subject of the study, inaccuracies within the QM / MM non-electrostatic potential are oflesser concern.The monolayer / multilayer di ↵ erences could indicate that, when fully immersed in water,the nanoparticle undergoes a slow process of recrystallization, but this would have to be confirmedwith molecular dynamics and thermal sampling. An experimental EXAFS of Rajh et al. [44] reportedthe partial restoration of the octahedral Ti coordination when enediol ligands adsorb to the surface ofspherical TiO nanoparticles.With regards to the electronic properties of the NP in di ↵ erent surrounding conditions, it appearsevident that inclusion of exact exchange is crucial for providing the correct electronic structure.The exact exchange term partially corrects for the electronic self-interaction error that is present instandard GGA functionals, such as PBE, and causes a large underestimation of the band gap, as isevident from Figure 6. We wish to comment on the fact that the DFTB +/ AMBER approach used inprevious work [28] is clearly faster than the current approach that does not employ the Tight-Bindingapproximation. The previous approach allows introducing a much larger number of water molecules todescribe the bulk of water around the nanoparticle. However, it cannot correctly capture the electronicproperties of the system being based on standard density functional theory methods.When adding more water, up to a 1 nm thick layer, we do not observe a relevant e ↵ ect on theband gap, but we observe an e ↵ ect on the NP’s work function with a shift of about 0.2 eV with bothCF and TIP3P LJ parameters, due to the dipole orientation of the water molecules in the multilayer.In particular, one should consider that the surface dipole moment for a naked NP is negative andpointing outwards. In the first adsorbed water monolayer (Layer I in Figure 7), all molecules arebound through the water oxygen with the OH bonds directed towards the vacuum. This creates anopposite (positive) dipole moment (i.e., pointing towards the center of the NP), which destabilizes theband states by + ↵ ect because the average radial dipole moment, created by the multilayer,is negative, causing a shift back of 0.2 eV (compare the position of the top of the valence band ofFigure 6d with that in Figure 6c). It must be re-iterated that these observations pertain to a specificconfiguration of water molecules in the multilayer, as obtained through the relaxation from a startingrandom distribution (see Section 4.2.5). To assess whether these changes are carried over to a thermalaverage, one should run molecular dynamics trajectory / ies with extended sampling of configurations,which is however beyond the scope of this work and could be object of future studies. olecules , , 2958 11 of 24 Figure 7.
Top view of the cross-section of the water multilayer on the NP model after QM / MMgeometry optimization. Oxygen atoms of the adsorbed water molecules and OH groups are color-codedaccording to their distance from the closest titanium atom. Titanium, oxygen, and hydrogen atomsof the nanoparticle are shown as cyan, red, and white spheres, respectively. Hydrogen bonds arerepresented by dashed black lines.
4. Materials and Methods / MM In electrostatic embedding QM / MM methodologies, the total energy of the entire system is a sumof the two subsystems, plus an interaction term: E TOT = E QM + E QM / MM + E MM (2)The MM region is here described using a point charge-based force field, with q i denoting the chargesand R i their spatial coordinates and n ( r ) the electronic density of the QM subsystem. The interactionenergy term E QM / MM (in atomic units) is: E QM / MM = N MM X i = q i Z n ( r ) | r R i | d r + N MM X i = N QM X ↵ = q i Z ↵ | R ↵ R i | + E NES (3)with Z ↵ being the atomic number of atom ↵ , running over all atoms in the QM subsystem. i runsover all N MM charges of the MM subsystem. In CRYSTAL17, the point charges are then formallyconsidered “atoms” with zero mass, nuclear charge as given in the input, and the default assigned basisset is a single s gaussian with an exponent of 1,000,000. Thus, the point charges are included in theexternal potential, and the electronic density is converged under their influence. Likewise, CRYSTAL17calculates the electrostatic forces from the density of the QM subsystem on the point charge centres. olecules , , 2958 12 of 24 The non-electrostatic electronic interactions, E NES , are in this work included througha Lennard–Jones pair potential between the atomic centers of each subsystem: E NES = N MM X i N QM X ↵ ✏ i ↵ "✓ i ↵ | R ↵ R i | ◆ ✓ i ↵ | R ↵ R i | ◆ (4)with ✏ i ↵ defining the well-depth, and i ↵ the equillibrium distance of the potential. Thus, nonon-electrostatic interactions are directly taken into consideration in the evaluation of the electronicdensity of the QM subsystem. The re-parameterizations via the DF and the CF strategies are carriedout for ✏ i ↵ and i ↵ of Equation (4), for the atom types (elements) Oxygen and Hydrogen, combined viathe Lorenz–Berthelot combining rules. For the QM subsystems in all of the QM / MM simulations carried out for this work, the choiceof functional has been (the CRYSTAL17 implementation of) B3LYP, comprised of Becke’s exchangefunctional [45], 20% HF exchange, and the Lee-Yang-Parr correlation functional [46]. For the benchmarks,we opted for the CRYSTAL17 basis set 6-31G(d,p) [47], which resembles the basis set used in the workwhere we obtained the non-water QM / MM LJ parameters from [33]. The Anderson ConvergenceAccelerator was used for all CRYSTAL simulations [13], and the CRYSTAL standard convergence criteriafor single point calculations were used. For water simulations including more than a single QM watermolecule, the counterpoise correction scheme was employed to avoid basis set superposition errors:when calculating monomer energies of single molecules in clusters, ghost atoms were added at thepositions of the other molecules, so that the same quality of basis set was kept constant, when evaluatinginteraction energies. As in previous work [11], the ASE-native TIP3P implementation [30,48] was usedfor the pure water benchmarks, thus avoiding the process of generating AMBER input files for these6000 + calculations.4.2.1. Water Dimers and ClustersThe interaction energy, or binding energy, is defined as D E int = E cluster P nm E m , where E m is theenergy of monomer m in the cluster, and n is the total number of molecules in the cluster. Thus, e.g., fora QM / MM water dimer using TIP3P as the MM model, only the QM monomer has a monomer energy.To transferably and generally minimize over- or underbinding of QM / MM systems in relation tothe total di ↵ erence between pure B3LYP and TIP3P models, we go beyond the water dimer, and extendthe water dataset to encompass binding energy calculations of all N possible QM / MM combinations of n -molecule water clusters: N = X n = n X k = n ! k ! ( n k ) ! c n (5)where k is the number of QM molecules, and c n is the number of clusters with n total molecules in thedataset. The water cluster geometries were obtained from the work by the Bates and Tschumper [49],and the set of (H O) n , n = / MM single point energy calculations. First, a dataset of interaction energy di ↵ erences DD E = D E QMMMint D E B3LYPint was produced using the TIP3P LJ parameters. The resulting datasetwill then consist of a distribution of interaction energies per k , for each n . This means quite a fewdistributions will have to be held up against each other, which is most easily visualized through abox-plot, where the distribution is simplified by showing the interquartile range (IQR), or the middle50% of the distribution as a box, as done in Figure 3. olecules , , 2958 13 of 24 Since only electrostatic interactions are included in the external potential in CRYSTAL17, one cansimply subtract the non-electrostatic contributions to the interaction energies, and re-calculate the LJterm with any given LJ parameters, to achieve new total interaction energies, without having to redothe more expensive DFT calculations. This method was employed to obtain fitted LJ parameters thatminimized the total the over- and underbinding of the entire water cluster dataset, as described in thefollowing section.4.2.2. Fitting StrategiesThe dimer fit minimizes the di ↵ erence between the QM / MM dimer binding curve and theaverage of the pure QM and the pure MM also in the 2.1 < r < / MM non-electrostatic coupling term (Equation (4)). The MM-MMLJ interactions are not modified, i.e., these interactions are still evaluated using the TIP3P parameters.The fit was done using the curve_fit tool from the Scipy.optimize submodule for python. The intervalis chosen to avoid giving weight to the very repulsive region at r < ↵ erence, , summed over all i clusterconfigurations, and normalized with the number of molecules in each cluster N i : = X i ⇣ D E QMMMint, i D E avgint, i ⌘ N i D E avgint, i (6)where D E avgint, i is the average of the pure QM and pure MM interaction energy of cluster configuration i .The fit was carried out by using the basin hopping algorithm implemented in the Scipy python module,starting from TIP3P LJ parameters, but also allowing for non-zero LJ parameters on the hydrogenatoms, for a total of 4 degrees of freedom in the fit. The initial step size was 1 ⇥ eV or Å. Keepingin mind the physical meaning of the - and ✏ -parameters in the LJ potential, the following constraintswere employed: 1 Å < OO < HH < , with theparameters CFOO = ✏ CFOO = / MM liquid water simulations, a cubic box with 27.96 Å sides containing729 water molecules was equilibrated at 300 K with the TIP3P force field, using a Langevin thermostatwith a friction coe cient of 0.02 a.u. After the equilibration was completed, the simulation wascontinued to produce initial configurations for spawning QM / MM MD simulations at 1 ps intervalsfrom the equilibrated MM MD trajectory. Three sets of QM / MM MD simulations were prepared:(1) using the TIP3P LJ parameters; (2) using the DF LJ parameters; and (3) using the CF parameters.The simulation setups were otherwise identical: The QM subsystem was defined to include a singlewater molecule, and re-equilibration was carried out for 4 ps for each of the spawned QM / MM MDtrajectories before sampling of, e.g., RDFs, was carried out. The RDFs were calculated using the VMDprogram [51]. The total sampled time for the CF LJ parameter-runs were 182 ps, 123 ps for the DF LJparameter-runs, and 56 ps for the TIP3P runs. The results can be found in Appendix C. olecules , , 2958 14 of 24 / Water Dimer BenchmarksTo benchmark the accuracy of our QM / MM methodology against pure DFT results on organicmolecules, we used a subset of the database used by Freindorf et al. [33]. All dimers consisting of asingle organic molecule and a water molecule from the original database were recreated, as well astwo new systems, namely glycine and aspartic acid, were included in our dataset. Their geometrieswere then optimized with B3LYP / / MM geometryoptimizations were performed, with the MM subsystem consisting of the water molecule, using theoriginal TIP3P LJ parameters, since the new parameters are only optimized for QM / MM water–waterinteractions. In these calculations, the MM atoms were fixed at their positions obtained from theprevious calculations of the dimers, performed with the single-description QM level of theory.The geometry relaxations were carried out until the maximum force magnitude within the systemswas below 0.01 eV / Å. The results can be found in Appendix B.4.2.5. Nanoparticle SimulationsThe anatase TiO spherical nanoparticle (NP) model used in this work was obtained from asimulated annealing process described in a previous work of some of us [24]. The model has astoichiometry of (TiO ) · O and it is characterized by a diameter of about 2.2 nm (see Figure 4a).A water monolayer around the NP was modeled by adding to the undercoordinated Ti atoms of thesurface water molecules. The procedure is described in detail in ref. [28]. The most stable watermonolayer around the NP is constituted by 134 water molecules and has an extent of dissociation ↵ = ↵ is defined as: ↵ = n OH,H n tot , n tot = n H O + n OH,H (7)where n H O is the number of molecular and n OH,H of dissociated water molecules in the watermonolayer. Finally, a spherical water shell of about 1 nm thickness was added around the NP with thewater monolayer using the PACKMOL code [52]. The water density inside the shell is approximatelyof 1 g / cm corresponding to 824 water molecules (see Figure 4c). For all the DFT calculations in thissection, the all–electron basis sets used are O 8–411(d1), Ti 86–411(d41), and H 511(p1), as definedin [23]. In the case of the bare NP and the NP with the water monolayer, forces were relaxed to lessthan 0.02 eV / Å using the CRYSTAL17 code. For the bare NP, we also used the generalized gradientapproximated (GGA) functional PBE [53], for comparison to B3LYP results. In the case of the NPsurrounded by the multilayer, forces were relaxed to 0.05 eV / Å to decrease the overall computationalcost, using CRYSTAL17 / AMBER16 QM / MM approach. For the QM part, we used the B3LYP functional,while to describe the QM / MM non-electrostatic interaction we made two simulations, one using CFand one using TIP3P LJ parameters.To simulate the direct-space extended X-ray adsorption fine structure (EXAFS) spectra,Gaussian convolution of peaks ( = eq = ax = --- Ti = --- Ti = --- O = cients ( = cients in thelinear combination of atomic orbitals (LCAO) of each molecular orbital, since summing the squares ofthe coe cients of all the atomic orbitals centred on a certain atom type results (after normalization) inthe relative contribution of each atom type to a specific eigenstate. olecules , , 2958 15 of 24
5. Conclusions
With this work, we have made a QM / MM electrostatic embedding implementation that combines(MPP-)CRYSTAL17 and AMBER16, and tested it thoroughly. We have produced new water–waterQM / MM LJ parameters for B3LYP / TIP3P, which improve the water–water coupling over the QM / MMborder, and could help reduce problems in, e.g., future adaptive QM / MM methodologies, where thecontents of the QM subsystem is dynamically updated. Through this part of the study, we have alsomade it clear that for solvents such as water, where the water–water interactions over the borderare important, it is not always su cient to optimize new LJ parameters from dimer interactions.This is of extra importance when the fit includes matching an average of QM and MM potentials,where the latter is already optimized to reproduce bulk phase properties, and thus have overemphasisedgas-phase interactions.We have also analysed the transferability of LJ parameters between di ↵ erent DFT codes andQM / MM implementations, and found that the QM / MM vs. pure-QM errors did not increase significantly,even though the parameters were produced with a di ↵ erent DFT code.Finally, we have demonstrated that the implementation can handle geometry optimizations ofTiO nanoparticles immersed in water. In the same vein of assessing how large a role the QM / MMLJ parameters play in the total picture, we have shown that changing them has a negligible e ↵ ect onthe NP. This is perhaps not surprising, since the di ↵ erences do not directly involve the NP itself: theMM-to-QM-water forces would have to a ↵ ect the QM water shell structure rather drastically for this toa ↵ ect the NP.In addition, as the studies on neat QM / MM water have shown, if attempts to reduce computationalcosts further were to be made by describing all water layers with force fields, it is likely thatre-parameterization of Ti and NP oxygen–water LJ parameters would be necessary, apart from alsohaving to use a force field for water that can model dissociation.Our analysis showed that the major changes to the nanoparticle take place already after addingthe first water monolayer, and that adding further water layers can have a small e ↵ ect on the structuralproperties, here seen in terms of a finer recrystallization. For the specific water configurationsstudied here, we observed a more evident e ↵ ect on the electronic properties of the NP. Future workcould be focused on studying if and how these e ↵ ects could change with thermal sampling of thesolvated nanoparticles. Author Contributions:
D.S., G.F. and J.J.M. wrote the ASE interface to CRYSTAL17, with assistance from A.O.D.and B.C. A.O.D. carried out the water benchmarks and the water-LJ QM / MM fitting work. L.F. and D.S. carriedout the benchmark on organic molecule / water dimers. D.S. and C.D.V. carried out the NP simulations, and wrotethe sections of the paper pertaining to those. B.C. wrote the details of the charge embedding in CRYSTAL17.A.O.D wrote the rest of the paper. C.D.V. edited the manuscript and coordinated the work. Funding:
This research was funded by the Icelandic Research Fund (grant 174244-051) and VILLUM FONDEN,the European Research Council (ERC) under the European Union’s HORIZON2020 research and innovationprogramme (ERC Grant Agreement No [647020]).
Acknowledgments:
A.O.D. Would like to thank Jónsson, H. for discussions about fitting strategies. C.D.V. isgrateful to Lara Ferrighi, Massimo Olivucci, and Stefano Motta for fruitful discussions. A.O.D. Acknowledgesfunding from the Icelandic Research Fund (grant 174244-051) and VILLUM FONDEN. The project has receivedfunding from the European Research Council (ERC) under the European Union’s HORIZON2020 research andinnovation programme (ERC Grant Agreement No [647020]).
Conflicts of Interest:
The authors declare no conflict of interest.
Abbreviations
The following abbreviations are used in this manuscript:ADF Angular Distribution FunctionAIMD Ab Initio Molecular DynamicsBOMD Born–Oppenheimer Molecular Dynamics olecules , , 2958 16 of 24 CF Cluster FitDF Dimer FitDOS Density of StatesGGA Generalized Gradient ApproximationHF Hartree-FockIQR Interquartile RangeLCAO Linear Combination of Atomic OrbitalsLJ Lennard–JonesML MonolayermultiL MultilayerNP NanoparticlePDOS Projected Density of StatesQM / MM Quantum Mechanical / Molecular MechanicalRDF Radial Distribution FunctionRMSD Root Mean Square DeviationvdW van der Waals
Appendix A. Further Water QM / MM Benchmarks
Table A1.
Fully optimized geometries of the water dimer, using the various multiscale- andsingle-description configurations, and their corresponding interaction energies. The QM atomswere completely unconstrained while the MM atoms where constrained to the rigid TIP3P geometry.The optimizations were carried out until the maximum force component on any unconstrained atomwas below 0.01 eV / Å. Configuration r OO (Å) \ ( ↵ ) (deg.) \ (OOH) (deg.) D E int (kcal / mol) B3LYP 2.8803 51.7 6.3 / TIP3P (LJ: TIP3P) 2.6414 15.9 0.4 / TIP3P (LJ: DF) 2.7039 21.6 2.0 / TIP3P (LJ: CF) 2.7414 17.8 0.4 / B3LYP (LJ: TIP3P) 2.8578 61.3 0.1 / B3LYP (LJ: DF) 2.8502 72.2 3.2 / B3LYP (LJ: CF) 2.7859 61.0 0.1 / cc-pVQZ a a Multireference calculations data from the S22 dataset [43].
Table A1 shows results from performing geometry optimizations using the three di ↵ erent LJparameter sets. Note how the multiscale geometries most closely resemble each of the two single-modelgeometries depending what model is used to describe the lone pair-donating (hydrogen accepting)molecule in the water dimer. The CF LJ parameters provide the smallest di ↵ erence in optimal O-Odistance between the QM / MM and MM / QM configuration.
Appendix B. Organic Molecule / Water Dimer Benchmark
Appendix B.1. Structures and Nomenclature
Figure A1 provides a visual overview of the benchmarked dimers, and the hydrogen-bondinggeometries that the QM / MM interface needs to describe. olecules , , 2958 17 of 24 METHANE
CH HH H
Figure A1.
Structure and nomenclature of the di ↵ erent H-bond evaluated for di ↵ erent organicmolecules / water pairs: acetic acid / water (1,2,3,4), methylamine / water (1,2), acetone / water (1),dimethyl ether / water (1), methanol / water (1,2), methane / water (1), imidazole / water (1,2), glycine / water(1,2,3,4,5), and aspartic acid / water (1,2,3,4,5). It is debatable whether the methane–water interactioncan be called a hydrogen bond, but, since it is included in the original dataset, we also include it here. Appendix B.2. Benchmark Results
Table A2.
Comparison of the hydrogen bond binding energy ( E HB in kcal / mol), distance ( r HB in Å) andangle ( HB in degrees) as calculated with B3LYP and B3LYP / TIP3P. For the structure and nomenclatureof the hydrogen bond evaluated for di ↵ erent organic molecules / water and amino acids / water pairsrefer to Figure A1 in Appendix B.1. Molecule HB B3LYP B3LYP / TIP3P E HB r HB HB E HB r HB HB Acetic Acid 1 10.50 1.782 156.9 9.34 1.972 160.82 10.50 1.979 136.0 9.34 2.198 135.03 5.94 1.930 160.0 6.03 1.987 160.24 3.58 2.067 153.6 3.97 2.081 153.7 olecules , , 2958 18 of 24 Table A2.
Cont . Molecule HB B3LYP B3LYP / TIP3P E HB r HB HB E HB r HB HB Methyl Amine 1 2.68 2.188 158.7 3.20 2.175 166.32 7.90 1.908 171.2 6.59 2.055 173.4Acetone 1 6.31 1.907 164.6 6.47 1.948 164.9Dimethyl Ether 1 5.24 1.907 175.6 5.05 1.952 166.1Methanol 1 5.58 1.936 173.9 5.93 1.946 172.42 5.96 1.898 177.3 6.04 1.932 176.8Imidazole 1 6.50 1.953 179.2 6.77 1.964 178.52 7.08 1.949 179.4 6.05 2.092 177.0Methane 1 0.38 2.566 156.7 0.74 2.660 163.6Glycine 1 10.50 1.776 156.6 9.29 1.959 160.32 10.50 1.995 135.0 9.29 2.204 134.43 6.34 1.943 156.8 7.05 2.019 156.74 6.34 2.186 146.1 7.05 2.167 156.05 3.26 2.045 152.2 3.80 2.101 152.9Aspartic acid 1 6.56 2.155 135.1 7.03 2.321 134.12 6.56 2.563 111.2 7.03 2.820 108.13 5.66 2.178 149.1 6.70 2.167 148.84 5.66 2.477 126.0 6.70 2.479 126.25 7.29 1.893 159.9 7.26 1.968 159.8
Table A3.
Root-mean-square deviations (RMSD) for hydrogen bond binding energies ( E HB inkcal / mol), distances ( r HB in Å), and angles ( HB in degrees) between the pure B3LYP referenceand the B3LYP / TIP3P calculations.
RMSD ( E HB ) RMSD ( r HB ) RMSD ( HB ) Appendix B.3. Analysis
Figure A2 shows the correlation of hydrogen bond energies, distances, and angles betweenthe QM / MM simulations and the pure QM reference. The QM / MM angles and energies are evenlydistributed around their QM reference, and no systematic di ↵ erences can be observed for this dataset.For the hydrogen bond lengths, the correlation still seems linear, but with a tendency for the QM / MMmodel to systematically slightly overestimate the bond lengths. The total RMSD values for the datasetare 0.77 kcal / mol, 0.121 Å, and 3.7 degrees, for the binding energy, H-bond length, and H-bondangle, respectively. olecules , , 2958 19 of 24 E QMHB (kcal/mol) E Q M / MM H B ( k c a l / m o l ) H Donating WaterH Accepting Water r QMHB (˚A) r Q M / MM H B ( ˚A )
100 125 150 175 QMHB Q M / MM H B Figure A2. QM / MM vs. QM Correlations between hydrogen bond: binding energy ( top ); distance( middle ); and angle ( bottom ). The dashed black lines represent the pure QM reference values.
The authors report that their hydrogen bond energies are systematically too big, while the hydrogenbond distances are systematically too short [33]. Curiously, when we use their LJ parameters in ourCRYSTAL17 / AMBER16 QM / MM implementation, we only see, as shown in Figure A2, a systematicdeviation of the hydrogen bond distances, and, in our case, the QM / MM bonds are too long . There arethree main di ↵ erences between the work of Freindorf et al. and the work carried out here: firstly,the DFT code used by Freindorf et al. is Q-Chem, which is di ↵ erent from CRYSTAL17 in many aspects,and could have an implementation of B3LYP that is not exactly the same as in CRYSTAL17. The authorsdid not mention if any reduced-scaling algorithms were employed, which we assume must meanthat no such methods were used. However, some programs apply, e.g., the resolution of identityapproximation by default, which could a ↵ ect geometry optimizations, but, since no auxiliary basissets are reported, this is most likely not the case. Secondly, there is not a complete 1:1 correspondencebetween the basis sets available to the two codes, e.g., CRYSTAL17 uses 5d functions in the 6-31G(d,p)basis set, where many other codes would use 6d. Lastly, the geometry optimizations could have beencarried out with another type of optimization, and to di ↵ erent convergence criteria. olecules , , 2958 20 of 24 Appendix C. QM / MM Liquid Water
Figure A3 compares QM / MM water RDFs to the RDF obtained from a pure TIP3P run. Sincehybrid functionals such as B3LYP are still relatively costly, we have not made a pure B3LYP simulationfor comparison. Instead, we have shown a result from the literature [36], which was restricted to32 molecules and short trajectories due to cost, and run at 350 K, which impedes direct comparison.It is a well known strategy [54–56] to increase the temperature in AIMD simulations to counteractboth the overstructuring of liquid water by many DFT functionals, and to nullify the lack of protonquantum e ↵ ects in the dynamics. The overbinding in the TIP3P LJ parameter dimer binding curve isobserved to carry over to liquid water, and results in an overstructuring: the amplitude of both the firstpeak, and the first well is the largest of all the datasets. The first peak of the DF and CF datasets arealmost of the same height, with a peak magnitude of 2.92 and 2.86 for the CF and DF set, respectively.Both the first well and the second peak of the DF dataset have higher amplitudes than the CF dataset.The CF radius of the first solvation shell seems to lie in between the pure B3LYP and pure TIP3P result,whereas the shell is thinner for the two other QM / MM simulations.
Figure A3.
Pairwise O-O, O-H, and H-H radial distribution functions of liquid water at 300 K.The QM / MM RDFs are sampled from a single B3LYP water molecule in a box of 728 MM molecules.Three types of QM / MM simulations were made, one for each of the QM / MM LJ parameter sets: TIP3P,DF, and CF. The pure B3LYP data, simulated at 350 K, and for 32 water molecules was obtainedfrom [36]. The higher temperature is known to lower the amplitude of the O-O peaks for both classical,pairwise force fields [54], and for DFT methodologies [55,56].
Angular distributions of the hydrogen bond accepting and donating water molecules are centralquantities for characterizing the hydrogen bond network in liquid water. Figure A4 shows the sampleddistributions for the three QM / MM water datasets produced in this work. Again, as was the case forthe radial distributions, the dimer fitting methodology has reduced the water over-structuring relativeto the TIP3P-structure, but not to the same degree as the cluster fit methodology. Especially for the ↵ (and the connected ) angle, not only the amplitude has been decreased towards the TIP3P result,but also the peak position has moved. olecules , , 2958 21 of 24
100 120 140 160 180 (degrees) ( ) B3LYP/TIP3P (LJ: CF)B3LYP/TIP3P (LJ: DF)B3LYP/TIP3P (LJ: TIP3P)TIP3P
25 50 75 100 125 150 175 (degrees) ( ) (degrees) ( ) Figure A4. QM / MM and pure MM angular distribution functions. See Figure 1 for definitions. The ↵ and angles are most a ↵ ected by the change of QM / MM LJ parameters to the CF set, whereas the DFparameters provide angular orientations that are similar to the structure from the TIP3P parameters.
References
1. Warshel, A.; Levitt, M. Theoretical Studies of Enzymatic Reactions: Dielectric, Electrostatic and StericStabilization of the Carbonium Ion in the Reaction of Lysozyme.
J. Mol. Biol. , , 227–249. [CrossRef]2. Lin, H.; Truhlar, D.G. QM / MM: What have we learned, where are we, and where do we go from here?
Theor. Chem. Acc. , , 185–199. [CrossRef]3. Senn, H.M.; Thiel, W. QM / MM methods for biomolecular systems.
Angew. Chem. Int. Ed. , , 1198–1229. [CrossRef] [PubMed]4. Bulo, R.E.; Michel, C.; Fleurat-Lessard, P.; Sautet, P.; Heyden, A.; Lin, H.; Truhlar, D.G.; Pezeshki, S.; Lin, H.;Riahi, S.; et al. Multiscale Modeling of Chemistry in Water: Are We There Yet? J. Chem. Theory Comput. , , 2231–2241. [CrossRef] [PubMed]5. Duster, A.W.; Wang, C.H.; Garza, C.M.; Miller, D.E.; Lin, H. Adaptive quantum / molecular mechanics:What have we learned, where are we, and where do we go from here? WIRES Comput. Mol. Sci. , , e1310. [CrossRef]6. Zhang, Y.J.; Khorshidi, A.; Kastlunger, G.; Peterson, A.A. The potential for machine learning in hybridQM / MM calculations.
J. Chem. Phys. , , 241740. [CrossRef] [PubMed]7. Steinmann, C.; Reinholdt, P.; Nørby, M.S.; Kongsted, J.; Olsen, J.M.H. Response properties of embeddedmolecules through the polarizable embedding model. Int. J. Quantum Chem. , , e25717. [CrossRef] olecules , , 2958 22 of 24
8. Morzan, U.N.; Alonso de Armiño, D.J.; Foglia, N.O.; Ramírez, F.; González Lebrero, M.C.; Scherlis, D.A.;Estrin, D.A. Spectroscopy in Complex Environments from QM–MM Simulations.
Chem. Rev. , , 4071–4113. [CrossRef] [PubMed]9. Dohn, A.O.; Jónsson, E.O.; Kjær, K.S.; van Driel, T.B.; Nielsen, M.M.; Jacobsen, K.W.; Henriksen, N.E.;Møller, K.B. Direct Dynamics Studies of a Binuclear Metal Complex in Solution: The Interplay BetweenVibrational Relaxation, Coherence, and Solvent E ↵ ects. J. Phys. Chem. Lett. , , 2414–2418. [CrossRef][PubMed]10. van Driel, T.B.; Kjær, K.S.; Hartsock, R.W.; Dohn, A.O.; Harlang, T.; Chollet, M.; Christensen, M.;Gawelda, W.; Henriksen, N.E.; Kim, J.G.; et al. Atomistic characterization of the active-site solvationdynamics of a model photocatalyst. Nat. Commun. , , 13678. [CrossRef] [PubMed]11. Dohn, A.O.; Jónsson, E.Ö.; Levi, G.; Mortensen, J.J.; Lopez-Acevedo, O.; Thygesen, K.S.; Jacobsen, K.W.;Ulstrup, J.; Henriksen, N.E.; Møller, K.B.; et al. Grid-Based Projector Augmented Wave (GPAW)Implementation of Quantum Mechanics / Molecular Mechanics (QM / MM) Electrostatic Embedding andApplication to a Solvated Diplatinum Complex.
J. Chem. Theory Comput. , , 6010–6022. [CrossRef][PubMed]12. Erba, A.; Baima, J.; Bush, I.; Orlando, R.; Dovesi, R. Large-Scale Condensed Matter DFT Simulations:Performance and Capabilities of the CRYSTAL Code. J. Chem. Theory Comput. , , 5019–5027.[CrossRef] [PubMed]13. Dovesi, R.; Saunders, V.R.; Roetti, C.; Olando, R.; Zicovich-Wilson, C.M.; Pascale, F.; Civalleri, B.; Doll, K.;Harrison, N.M.; Bush, I.J.; et al. CRYSTAL17 User’s Manual ; University of Torino: Torino, Italy, 2017.14. Dovesi, R.; Erba, A.; Orlando, R.; Zicovich-Wilson, C.M.; Civalleri, B.; Maschio, L.; Rérat, M.; Casassa, S.;Baima, J.; Salustro, S.; et al. Quantum-mechanical condensed matter simulations with CRYSTAL.
Wiley Interdiscip. Rev. Comput. Mol. Sci. , , e1360. [CrossRef]15. Labat, F.; Baranek, P.; Domain, C.; Minot, C.; Adamo, C. Density functional theory analysis of thestructural and electronic properties of TiO rutile and anatase polytypes: Performances of di ↵ erentexchange-correlation functionals. J. Chem. Phys. , , 154703. [CrossRef] [PubMed]16. Bai, Y.; Mora-Seró, I.; De Angelis, F.; Bisquert, J.; Wang, P. Titanium Dioxide Nanomaterials for PhotovoltaicApplications. Chem. Rev. , , 10095–10130. [CrossRef] [PubMed]17. Ma, Y.; Wang, X.; Jia, Y.; Chen, X.; Han, H.; Li, C. Titanium Dioxide-based Nanomaterials for PhotocatalyticFuel Generations. Chem. Rev. , , 9987–10043. [CrossRef] [PubMed]18. Schneider, J.; Matsuoka, M.; Takeuchi, M.; Zhang, J.; Horiuchi, Y.; Anpo, M.; Bahnemann, D.W.Understanding TiO Photocatalysis: Mechanisms and Materials.
Chem. Rev. , , 9919–9986.[CrossRef] [PubMed]19. Rajh, T.; Dimitrijevic, N.M.; Bissonnette, M.; Koritarov, T.; Konda, V. Understanding TiO Photocatalysis:Mechanisms and Materials.
Chem. Rev. , , 10177–10216. [CrossRef] [PubMed]20. Diebold, U. Perspective: A Controversial Benchmark System for Water-oxide Interfaces: H O / TiO (110). J. Chem. Phys. , , 040901. [CrossRef] [PubMed]21. Mu, R.; Zhao, Z.j.; Dohnálek, Z.; Gong, J. Structural Motifs of Water on Metal Oxide Surfaces. Chem. Soc. Rev. , , 1785–1806. [CrossRef] [PubMed]22. De Angelis, F.; Di Valentin, C.; Fantacci, S.; Vittadini, A.; Selloni, A. Theoretical Studies on Anatase and LessCommon TiO Phases: Bulk, Surfaces, and Nanomaterials.
Chem. Rev. , , 9708–9753. [CrossRef][PubMed]23. Fazio, G.; Ferrighi, L.; Di Valentin, C. Spherical versus Faceted Anatase TiO Nanoparticles: A ModelStudy of Structural and Electronic Properties.
J. Phys. Chem. C , , 20735–20746. [CrossRef]24. Selli, D.; Fazio, G.; Di Valentin, C. Modelling Realistic TiO Nanospheres: A Benchmark Study of SCC-DFTBagainst DFT.
J. Chem. Phys. , , 164701. [CrossRef] [PubMed]25. Selli, D.; Fazio, G.; Di Valentin, C. Using Density Functional Theory to Model Realistic TiO Nanoparticles,Their Photoactivation and Interaction with Water.
Catalysts , , 357. [CrossRef]26. Shirai, K.; Fazio, G.; Sugimoto, T.; Selli, D.; Ferraro, L.; Watanabe, K.; Haruta, M.; Ohtani, B.; Kurata, H.;Di Valentin, C.; et al. Water-Assisted Hole Trapping at Highly Curved Surface of Nano-TiO Photocatalyst.
J. Am. Chem. Soc. , , 1415–1422. [CrossRef] [PubMed] olecules , , 2958 23 of 24
27. Li, G.; Li, L.; Boerio-Goates, J.; Woodfield, B.F. High Purity Anatase TiO Nanocrystals:Near Room-Temperature Synthesis, Grain Growth Kinetics, and Surface Hydration Chemistry.
J. Am.Chem. Soc. , , 8659–8666. [CrossRef] [PubMed]28. Fazio, G.; Selli, D.; Seifert, G.; Di Valentin, C. Curved TiO Nanoparticles in Water: Short (Chemical) andLong (Physical) Range Interfacial E ↵ ects. ACS Appl. Mater. Interfaces , , 29943–29953. [CrossRef][PubMed]29. Bahn, S.R.; Jacobsen, K.W. An object-oriented scripting interface to a legacy electronic structure code. Comput. Sci. Eng. , , 55–66. [CrossRef]30. Larsen, A.; Mortensen, J.; Blomqvist, J.; Castelli, I.; Christensen, R.; Dulak, M.; Friis, J.; Groves, M.;Hammer, B.; Hargus, C.; et al. The Atomic Simulation Environment—A Python library for working withatoms. J. Phys. Condens. Matter , , 273002. [CrossRef] [PubMed]31. Hunt, D.; Sanchez, V.M.; Scherlis, D.A. A quantum-mechanics molecular-mechanics scheme for extendedsystems. J. Phys. Condens. Matter , , 335201. [CrossRef] [PubMed]32. Laio, A.; VandeVondele, J.; Rothlisberger, U. A Hamiltonian Electrostatic Coupling Scheme for HybridCar–Parrinello Molecular Dynamics Simulations. J. Chem. Phys. , , 6941–6947. [CrossRef]33. Freindorf, M.; Shao, Y.; Furlani, T.R.; Kong, J. Lennard–Jones parameters for the combined QM / MM methodusing the B3LYP / / AMBER potential.
J. Comput. Chem. , , 1270–1278. [CrossRef] [PubMed]34. Gillan, M.J.; Alfè, D.; Michaelides, A. Perspective: How good is DFT for water? J. Chem. Phys. , , 130901. [CrossRef] [PubMed]35. Head-Gordon, T.; Hura, G. Water Structure from Scattering Experiments and Simulation. Chem. Rev. , , 2651–2670. [CrossRef] [PubMed]36. Todorova, T.; Seitsonen, A.P.; Hutter, J.; Kuo, I.W.; Mundy, C.J. Molecular Dynamics Simulation of LiquidWater: Hybrid Density Functionals. J. Phys. Chem. B , , 3685–3691. [CrossRef] [PubMed]37. Seitsonen, A.P.; Bryk, T. Melting temperature of water: DFT-based molecular dynamics simulations withD3 dispersion correction. Phys. Rev. B , , 184111. [CrossRef]38. Horn, H.W.; Swope, W.C.; Pitera, J.W.; Madura, J.D.; Dick, T.J.; Hura, G.L.; Head-Gordon, T. Development ofan Improved Four-Site Water Model for Biomolecular Simulations: TIP4P-Ew. J. Chem. Phys. , , 9665–9678. [CrossRef] [PubMed]39. Wikfeldt, K.T.; Batista, E.R.; Vila, F.D.; Jónsson, H. A transferable H O interaction potential based on asingle center multipole expansion: SCME.
Phys. Chem. Chem. Phys. , , 16542–16556. [CrossRef][PubMed]40. Medders, G.R.; Babin, V.; Paesani, F. Development of a “First-Principles” Water Potential with FlexibleMonomers. III. Liquid Phase Properties. J. Chem. Theory Comput. , , 2906–2910. [CrossRef] [PubMed]41. Cisneros, G.A.; Wikfeldt, K.T.; Ojamäe, L.; Lu, J.; Xu, Y.; Torabifard, H.; Bartók, A.P.; Csányi, G.; Molinero, V.;Paesani, F. Modeling Molecular Interactions in Water: From Pairwise to Many-Body Potential EnergyFunctions. Chem. Rev. , , 7501–7528. [CrossRef] [PubMed]42. Babin, V.; Leforestier, C.; Paesani, F. Development of a “First Principles” Water Potential with FlexibleMonomers: Dimer Potential Energy Surface, VRT Spectrum, and Second Virial Coe cient. J. Chem.Theory Comput. , , 5395–5403. [CrossRef] [PubMed]43. Jure˘cka, P.; ˘Sponer, J.; ˘Cern ý , J.; Hobza, P. Benchmark database of accurate (MP2 and CCSD(T) completebasis set limit) interaction energies of small model complexes, DNA base pairs, and amino acid pairs. Phys. Chem. Chem. Phys. , , 1985–1993. [CrossRef] [PubMed]44. Rajh, T.; Chen, L.X.; Lukas, K.; Liu, T.; Thurnauer, M.C.; Tiede, D.M. Surface Restructuring of Nanoparticles:An E cient Route for Ligand-Metal Oxide Crosstalk. J. Phys. Chem. B , , 10543–10552. [CrossRef]45. Becke, A.D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. , , 5648–5652. [CrossRef]46. Lee, C.; Yang, W.; Parr, R.G. Development of the Colle-Salvetti correlation-energy formula into a functionalof the electron density . Phys. Rev. B , , 785–789. [CrossRef]47. Gatti, C.; Saunders, V.; Roetti, C. Crystal-field e ↵ ects on the topological properties of the electron-densityin molecular-crystals. The case of urea. J. Chem. Phys. B , , 10686–10696. [CrossRef]48. Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of simple potentialfunctions for simulating liquid water. J. Chem. Phys. ,79