Daniele
Loco
*a,
Louis
Lagardère
bc,
Gérardo A.
Cisneros
d,
Giovanni
Scalmani
e,
Michael
Frisch
e,
Filippo
Lipparini
f,
Benedetta
Mennucci
f and
Jean-Philip
Piquemal
*agh
aSorbonne Université, CNRS, Laboratoire de Chimie Théorique, LCT, Paris, France. E-mail: jean-philip.piquemal@sorbonne-universite.fr; daniele.loco@lct.jussieu.fr
bSorbonne Université, CNRS, Institut Parisien de Chimie Physique et Théorique, IP2CT, Paris, France
cSorbonne Université, Institut des Sciences du Calcul et des Données, ISCD, Paris, France
dUniversity of North Texas, Department of Chemistry, TX, USA
eGaussian, Inc., Wallingford, CT, USA
fUniverisita di Pisa, Dipartimento di Chimica e ChimicaIndustriale, Pisa, Italy
gInstitut Universitaire de France, IUF, Paris, France
hThe University of Texas at Austin, Department of Biomedical Engineering, TX, USA
First published on 11th June 2019
In this work, we present a general route to hybrid Quantum Mechanics/Molecular Mechanics (QM/MM) Molecular Dynamics for complex systems using a polarizable embedding. We extend the capabilities of our hybrid framework, combining the Gaussian and Tinker/Tinker-HP packages in the context of the AMOEBA polarizable force field to treat large (bio)systems where the QM and the MM subsystems are covalently bound, adopting pseudopotentials at the boundaries between the two regions. We discuss in detail the implementation and demonstrate the global energy conservation of our QM/MM Born–Oppenheimer molecular dynamics approach using Density Functional Theory. Finally, the approach is assessed on the electronic absorption properties of a 16500 atom complex encompassing an organic dye embedded in a DNA matrix in solution, extending the hybrid method to a time-dependent Density Functional Theory approach. The results obtained comparing different partitions between the quantum and the classical subsystems also suggest that large QM portions are not necessary if accurate polarizable force fields are used in a variational formulation of the embedding, properly including the QM/MM mutual polarization.
The definition of the target of such systems is not always a straightforward task, especially when the target and the environment are linked together by a covalent bond. As a matter of fact, in biological matrices (e.g.: proteins, DNA) several covalent bonds between the QM and the classical subsystem can be involved. It is evident that in all these cases, the way the system is partitioned becomes crucial to achieve a correct description of the process of interest. For instance, peculiar strong specific interactions between a small portion of the environment and the target may require to extend the QM treatment to a larger part of the system.10–13
The modeling of mixed QM/MM bonds has been a long-lasting argument of debate in the community and several different approaches have been developed in the years to address this issue.3,4,14,15 All these models have to solve, in some way, three main problems. First, cutting a covalent bond results in a QM atom with an unrealistic, highly reactive, unsaturated valence configuration, which can introduce strong artifacts in the simulation. Such an atom has to be capped in order to complete its valence. Second, the QM density at the boundary can be overpolarized when an electrostatic embedding is applied, which can be aggravated when a polarizable embedding is used as also the classical part can overpolarize. Third, the bonding terms in the classical MM force fields involving atoms from both subsystems have to be carefully selected in order to avoid double-counting of interactions.
Three different classes of boundary schemes can be found in the literature. A first class is made by link-atom schemes, where an additional (usually hydrogen) atom is introduced in the QM subsystem to cap the missing valence. Such an additional atom is not part of the real system and the charge next to the QM subsystem may be blurred or redistributed.16–18 The second one is made by localized orbital schemes,1,15,19 which employ hybrid orbitals placed at the boundary to cap the QM region, keeping some of them frozen. Finally, the third class uses a capping potential to terminate the QM region.20–23 In this work we will focus on a specific version of this third class, the pseudobond scheme.20–22
This implementation is integrated in our polarizable QM/MM method, based on the AMOEBA polarizable (and multipolar) force field,24–29 which interfaces the Gaussian and Tinker packages to achieve molecular dynamics (MD) on the adiabatic ground state energy surface of the QM/MM system. This new implementation is extensively tested on a set of oligopeptides, for which we discuss the energy conservation along hybrid ground state trajectories. Three test systems of increasing complexity are considered to show the stability of the method: an alanine dipeptide, a 6 aminoacid oligopeptide (SAPPAS) and a 20 aminoacid oligopeptide (PEP). We then detail the applicability of our method to the computation of excitation properties on a large solvated system encompassing a chromophore intercalated in DNA for which we will compute electronic excitations within a sequential approach based on the hybrid QM/MM molecular dynamics.
In this work we apply the recently proposed PBs by Yang and co-workers.22,43 The use of PBs for QM/AMOEBA hybrid simulations has been previously reported and implemented in the LICHEM code.44,45 Here, we employ the same approach as developed by Kratz et al.44
The ECPs are of the form
(1) |
A remarkable advantage of the PB approach is that it does not require any additional term to be implemented in the electronic structure code, as long as such a code is able to handle ECPs, which is usually the case. On the other hand, a significant effort has been put into the design of a flexible QM/MM interface, capable of handling the double quantum and classical nature of the PB atom. Such a duality requires in fact the handling of its associated electronic density within the electronic structure code (Gaussian 09), and the inclusion of bonding interactions involving quantum atoms in the MD code (Tinker), which normally ignores such class of atoms in the energy computation.
(2) |
The QM/MM coupling term is further divided into two different contributions
(3) |
Finally, the last term in eqn (2), , includes all the electrostatic interactions between the classical sites and all the bonding interactions of the classical FF. The bonding interactions which are accounted for in include also those bonding terms arising from the QM and MM atoms connected across QM/MM boundary regions. To compute these bonding terms between classical and quantum atoms connected through a PB atom, the PB atom itself and the QM atoms involved are treated as classical atoms, using the atomic parameters defined in the classical FF for the atom type required. In Fig. 1 the bonding interactions are depicted in a schematic manner, representing respectively: bond stretching, angle bending decomposed into in-plane and out-of-plane components, bond-angle cross terms, and torsional rotation, including torsion–torsion coupling terms.
The term does not depend on either the electronic density or the induced dipoles, and for this reason it does not enter the SCF procedure. It is worth noting that the PB atoms do not need any classical multipoles or polarizabilities since the description of the PB atoms is included in the nuclear and electronic part of the QM subsystem and thus any electrostatic and polarization interactions involving the PB atoms are included in the hybrid terms of eqn (3).
Taking advantage of the variational formulation, the coupled QM/AMOEBA equations needed to plug the QM/MM scheme in the quantum chemistry method of choice have been derived, and an effective QM/AMOEBA Fock (or Kohn–Sham, depending on whether Hartree–Fock or DFT theory is concerned) matrix can be obtained as the gradient of the energy functional in eqn (2) with respect to the density matrix.33
Analytical derivatives of the QM/AMOEBA energy are obtained by differentiating the energy functional of eqn (2) with respect to both classical and quantum nuclear coordinates. The hybrid QM/MM forces can then be used to perform QM/MM Born–Oppenheimer molecular dynamics (BOMD).34 An extended Lagrangian approach (XL-BO), in the formulation proposed by Niklasson and co-workers,47–50 is applied to improve the guess for the electronic density at each integration time step of the dynamics to speed-up the convergence of the SCF procedure.34
The detailed derivation of the TDDFT equations in a polarizable environment can be found elsewhere.51 Here we only report the final modified TDDFT equations:52
(4) |
(5) |
To account for the further environment effects due to the relaxation of its polarization in response to change in the QM density in the excitation process, a so-called “state-specific” correction,53,54 completely equivalent to the so-called “corrected Linear Response” (c-LR) scheme originally developed for polarizable continuum models,55 was implemented in the QM/AMOEBA method.33 Within this framework, a relaxed density matrix is calculated for the excited state of interest and the corresponding excitation energy is corrected for the interaction of the corresponding electric field and the induced dipole moments within the environment.
The oligopeptides are arbitrarily partitioned in a QM portion, modeled with DFT, and an MM one, modeled with the AMOEBA polarizable potential. The QM portion consists of one of the alanine residues for ALA, two central prolines for SAPPAS and a four aminoacid sequence (tyrosine–histidine–tryptophan–aspartate) for PEP. The boundaries between the QM and the classical portions in the oligopeptide are treated with the adequate PB, depending on whether the PB formed on the QM side is a CPB–N or CPB–CO bond, where CPB indicates the carbon atom bearing the pseudopotential, and the parameterized basis set, and included in the QM portion as the boundary atom. Simulations where a QM/MM partition is performed within the PB scheme will be referred to in the rest of the manuscript with the tag QM/MM(PB).
When the PB approach is applied, the fixed multipoles and the atomic polarizability of the classical atoms directly bonded to the PB atom (MM1 sites according to Fig. 1) are set to zero and the point charges are distributed on the MM2 atoms to balance the total charge of the system, as suggested by Zhang and co-workers.20,21
ALA is the simpler and smaller system, and for this reason it is employed to fully explore the capabilities and the limits of the PB approach. For such a purpose two sets of simulations are performed on ALA: (i) gas-phase dynamics and (ii) dynamics in a small droplet of water, treated with the AMOEBA potential. For each set of simulations the QM/MM(PB) dynamics is compared with a full MM trajectory, where the ALA is described completely with the AMOEBA potential, and a full QM one, where the ALA is full DFT. The total energy conservation is compared between the different models used to describe the ALA system and also their effects on some structural parameters, whose adequate description could be crucial in the modeling of a more realistic chemical problem.
Energy conservation is discussed for the QM/MM(PB) dynamics of SAPPAS and PEP only in a water droplet (AMOEBA water), and compared with results obtained from a trajectory performed with a full MM oligopeptide. Those two systems, even though small compared to real-life MD applications, are indeed more complex than the simple ALA dipeptide. For this reason, and because the AMOEBA FF is parameterized only for simulations in solution, gas-phase tests on SAPPAS and PEP are of little interest to draw significant conclusions on the physics of the PB scheme.
All the trajectories are run for 2 ps at the B3LYP/6-31G(/AMOEBA) level of theory, in the NVE ensemble, with an integration time-step δt = 0.5 fs and initial velocities sampled from a Maxwell–Boltzmann distribution at 80 K. The integrator applied is the velocity Verlet. The low temperature is chosen to minimize the noise in the energy fluctuations.34
In Fig. 2 the partitioned system is depicted and the total energy Etot time-series of full QM, QM/MM(PB) and full MM trajectories is reported. Either in vacuo (@ QM(VAC), @ PB(VAC) and @ MM(VAC)) and in the small droplet of classical AMOEBA water (@ QM(WAT), @ PB(WAT) and @ MM(WAT)) the drift in the total energy is barely noticeable, even in the very small scale of the energy fluctuations due to the low temperature (∼±10−5–10−4 hartree (EH)). Both the drifts and the oscillations in the total energy appear slightly larger in the simulation in solution. This is due to the increased complexity in the dynamics, which stems from the presence of the solvent. Overall, the drift is still remarkably small (∼10−4 EH). Furthermore, no significant differences can be observed in the behaviors of the full QM, QM/MM(PB) and full MM simulations, meaning that the introduction of the PB scheme does not affect the total energy conservation significantly.
The “in vacuo” trajectories are analyzed more in detail from a structural point of view. To have a larger sampling the dynamics have been extended to 4 ps. The energy conservation of the elongated dynamics is still consistent with what was observed for the shorter ones and the equivalent analysis is reported in the ESI (Section S1†).
The three compared trajectories all start from the same geometry but exhibit different initial velocities, and thus they can, in principle, evolve in different regions of the PES. Since they evolve for a short time, and because of the relatively “simple shape” of the PES for the ALA system and the low temperature of the simulation, this effect is not expected to play a major role in the structural differences that can be detected.
First, two dihedral angles are analyzed as they are involved in bonding terms between quantum and classical atoms connected through the PB atom within the QM/MM(PB) dynamics. These dihedrals are reported in Fig. 3 and named CNCC (green) and NCCN (orange). We also report their distributions along the full QM (red), QM/MM(PB) (blue) and full MM (yellow) dynamics in the gas-phase. For the CNCC dihedral, where two of the involved atoms are in the MM region, the full MM and QM/MM(PB) distributions are almost overlapping.
Fig. 3 Comparison between selected structural features of ALA along a full MM (yellow), full QM (red), and the hybrid QM/MM(PB) (blue) dynamics. The distributions are drawn from 8000 snapshots extracted from 4 ps trajectories of ALA in the gas-phase. The structural parameters analyzed are: two hybrid QM/MM dihedrals occurring across the PB defined in the QM/MM partition (see Fig. 2), named CNCC (green) and NCCN (orange), the OCNC (blue) dihedral, which is fully QM in the QM/MM(PB) dynamics and the CCNC (cyan) one, which involves the PB atom and it is fully MM in the QM/MM(PB) dynamics. The length of the PB formed in the QM/MM partition scheme is also analyzed. The continuous lines are Gaussian fitting to the distribution. |
The scenario is quite different for the NCCN dihedral, where the two involved atoms are in the QM region (see Fig. 3). In this case the full QM distribution shows two main domains, centered around 50 and 90°, respectively, while both the QM/MM(PB) and the full MM span a smaller domain of dihedral angles, between 70 and 110°. The QM/MM(PB) and full MM distributions are mostly overlapping, even though the intensities of the single bins can vary a lot between them, and they also exhibit a similar shape, with two separate peaks in both cases centered around the intervals 70–80° and 90–100°.
The large difference between the full QM results and the other distributions can be certainly attributed to the mismatch between the DFT and the AMOEBA potentials, AMOEBA being fitted on a higher level of QM accuracy.27 The comparison also shows that the overall quality of the QM/MM dynamics is not significantly deteriorated with respect to the one based on the AMOEBA FF by the PB approach. To conclude the analysis of dihedrals in the QM/MM(PB) trajectory, two further ones, namely, one involving only QM atoms (OCNC, blue) and one only involving MM atoms and the pseudoatom (CCNC, cyan) are analyzed. In both cases the comparison between the distributions of such angles in the three different trajectories shows a good agreement, especially for the OCNC dihedral, where the QM/MM(PB) distribution strongly resembles a hybrid between the full MM and the full QM one. The OCNC case, however, points out that the effect of the PB, even though reasonably small, is not completely negligible in the close proximity of the QM/MM boundary.
A similar analysis has been carried out on the length of the bond which corresponds to the PB length defined in the QM/MM(PB) simulation, compared to the length of the corresponding bond in the full QM and full MM dynamics. The resulting distributions (Fig. 3) show a Gaussian behavior, with the full QM case (red) having a C–C bond which is on average longer than that of the corresponding PB (in blue), with an average difference of ∼0.05 Å. The full MM value (yellow) lies in between the QM/MM(PB) and the full QM. Considering the full QM simulation as a reference, the PB approach is reasonably in accordance. The differences in the bond lengths are well motivated by the effect of the smaller basis set for the valence electrons and the pseudopotential used to treat the core electrons.
It is worth reminding here that these analyses on the QM/MM boundaries represent a limiting case for the method since, as usual in QM/MM, the QM region of interest, where structural parameters can play a fundamental role, is always kept as far as possible from the QM/MM boundary region, in order to minimize the effects of the cut.
For the SAPPAS hexapeptide two PB atoms are defined, one on each of the two alanine units (see Fig. 4, top panel). In one case the PB replaces the α-C of a C–N bond in the quantum part, while in the other a C–CO bond is replaced, as for ALA. In Fig. 4 (Bottom panel) the SAPPAS total energy conservation along the full MM trajectory (yellow) is compared with that of the QM/MM(PB) dynamics (blue). The QM/MM(PB) energy time-series (SAPPAS @ PB(WAT)) slightly deviates from the full MM starting a positive drift after 0.5 ps of simulation. It suddenly disappeared with an opposite sign drift. Around 2 ps the drift has been removed. A longer 6 ps trajectory is analyzed in the ESI (Section S2†), showing a negative drift reappearing after 3 ps. The drift appears linear in time and of very small entity, with a negligible largest energy deviation from the average of ∼0.3 mEH (less than 0.2 kcal mol−1).
Finally, as a last numerical test, the PEP oligopeptide (see Fig. 4, top panel) trajectory is analyzed. The QM/MM(PB) total energy conservation (blue) is reported in comparison with that of the full MM dynamics (yellow) shifted by −0.2 mEH to avoid the superposition of the two time series. In this case, as in the ALA dynamics, no significant effect on the total energy conservation can be seen when the PB approach is applied. For PEP the fluctuation in the total energy is generally larger than those of ALA and SAPPAS, due to the larger dimension of the system (WAT simulations involve water boxes of: 630 molecules for ALA, 473 molecules for SAPPAS and 1053 molecules for PEP). Also for PEP a comparison between longer simulations (7.5 ps) is reported in the ESI (Section S2†).
Fig. 5 (A) Pictorial representation of the TO dye buried in the DNA double helix, embedded in a sphere of water (16500 atoms). The different colors represent the differences peculiar to each of the QM/MM partition schemes compared, highlighting the portions of the system included in the QM subsystem of the QM/MM (orange) and the QM/MM(PB) (blue) partition schemes, and the environment selection applied in the movE simulations which differs from that of the other schemes. See the text and Table 1 for a detailed definition of the QM/MM, QM/MM(PB) and movE schemes. (B) The DNA structure is highlighted in a ball-and-stick representation, leaving the first water solvation shell around the TO dye. The QM subsystem defined in this work, tagged as QM/MM(PB) and made up of the TO dye and the four closest nucleobases (NBs) is zoomed out and represented in licorice style in the blue square. The atoms bearing the pseudopotential on the QM/MM boundaries are represented in black. (C) NTOs and excitation energies relative to the ππ* bright excitation of the TO dye embedded in different environments, from vacuum to the QM/MM(PB) partition scheme. Also the results with a different approach to treat the QM/MM boundaries, namely the link atom (QM/MM capping H), are reported. All the NTOs and energies are computed on the initial configuration used for the MD2 simulation (see Table 2 and in the text). H and P mean the hole and particle, respectively. |
In a previous study by some of the authors58 the vibronic structure of the absorption band of TO in the biological matrix was simulated applying a sequential scheme based on the QM/AMOEBA MD performed with the Gaussian 09/Tinker interface. The TDDFT including explicit effects of the polarizable environment on the excitation properties was employed. In that study, the system was partitioned in such a way that only the TO dye was treated at the QM level, while the rest of the environment was represented with the AMOEBA FF. In the trajectories used to compute excitation energies, only the TO was allowed to move in a fully classical and frozen environment.58
Neither the role of the QM/MM partition scheme employed nor the effects of the frozen environment on the excitation properties have been investigated.
Aiming at shedding more light on these aspects, a new set of simulations is compared to the previous results, applying the PB scheme to assess how the excitation properties are affected by the definition of a larger QM subsystem which goes beyond the TO dye, including also the two pairs of nucleobases (NBs) nearest to the dye. This modified QM/MM partition scheme entangles together the effects on the dynamics of the system and on the resulting excitation properties. The PB approach is applied to the QM/MM boundaries, which are placed between the nucleobase and the sugar connecting the base to the DNA backbone (see Fig. 5). We refer to the new definition of the QM subsystem including the NBs as QM/MM(PB), while the previous one only including TO as QM/MM.
In order to better investigate how the environment relaxation affects the overall results, we performed a third type of MD simulation. Now, all the systems, not only the QM subsystem, is allowed to move; we label the test performed under these conditions as movE.
In Table 1 the notation used to define the different QM/MM partition schemes is reported, together with the composition of the QM and MM subsystems and whether the MM one is frozen or it is moving during the QM/MM MD simulation. In Fig. 5 a pictorial representation of the whole system and of the different schemes applied in this work is reported.
Simulation tag | QM/MM partition | Dynamics | |
---|---|---|---|
QM/MM | QM | TO dye | Moving |
MM | DNA + water | Fixed | |
QM/MM(PB) | QM | TO dye + four NBs | Moving |
MM | Rest of the DNA + water | Fixed | |
movE | QM | TO dye + four NBs | Moving |
MM | Rest of the DNA + water | Moving |
The structures coming from the QM/MM, QM/MM(PB) and movE dynamics are finally used to compute the excitation energies using the same QM–MM partition adopted for the dynamics. Since the present work does not aim at computing the whole absorption band shape of the embedded TO dye, instead of performing several hybrid QM/MM trajectories, and sampling the initial system configuration from a classical dynamics, a reduced set of starting geometries is extracted from the same classical trajectory used in ref. 58, obtaining two starting configurations, far in time from each other as schematically depicted in Fig. 6. The trajectories produced from the two initial configurations selected are named MD1 and MD2.
Fig. 6 Schematic representation of the procedure followed to compute the MD1 and MD2 trajectories. Two configurations of the TO dye embedded in the DNA in water solution have been extracted from the classical molecular dynamics (CMD, blue arrow) carried out in ref. 58. From these two configurations the MD1 and MD2 2 ps long trajectories have been run within the three different QM/MM partition schemes reported in Table 1: QM/MM, QM/MM(PB) and movE. Each trajectory, for a total of six simulations, has been sampled every 80 fs, removing the first 480 fs as an equilibration time, to extract system-environment configurations to perform TDDFT calculations. |
The B3LYP/6-31G/AMOEBA level of theory has been used for all molecular dynamics simulations. 2 ps runs in the NVE ensemble with a time step δt = 0.5 fs have been performed to keep consistency with the simulation proposed in ref. 58. No periodic boundary conditions were applied (i.e. water droplet), and the initial atomic velocities are assigned through the sampling of a Maxwell–Boltzmann distribution at 300 K.
First the effect of a larger QM subsystem is discussed. In these simulations the classical environment is kept frozen to be consistent with the previous study. The environment, composed of the DNA helix, counterions and water molecules is made up of about 16500 atoms, most of which are due to the large solvation box (see Fig. 5). In ref. 58 the B3LYP functional has been employed also for the TDDFT calculations which follow the dynamics. When the new QM/MM(PB) partition is applied, the use of the B3LYP functional results in the appearance of charge-transfer (CT) intruder states, delocalizing the TO ππ* transition of interest on the QM NBs. This effect has been investigated computing TDB3LYP/6-31+G(d)/AMOEBA excitations on the snapshots extracted from MD1 every 80 fs (see Fig. 6) and the results are reported in the ESI (Section S3†). Excitations computed on the same MD1 snapshots but restricting the QM subsystem at the TO dye only, where the CTs cannot take place in the TDDFT calculation, have been used for comparison and the nature of the excitation has been analyzed in terms of the corresponding Natural Transition Orbitals (NTO).59
In Section S4 of the ESI† we report the comparison between the B3LYP, CAM-B3LYP, M062X, PBE0, ωB97X and ωB97XD functionals, showing that the emergence of the intruder CT-like states can be cured by choosing a range-separated hybrid functional.60 For this reason the CAM-B3LYP exchange-correlation functional has been selected to compute excitation properties.
The results are reported in Table 2 in terms of the average excitation energies and the corresponding transition dipoles from MD1 and MD2 run both in the previous QM/MM and the new QM/MM(PB) schemes. The relative differences between the averaged values are also reported. It can be noticed that the differences between the excitation properties computed on MD1 and MD2 are much smaller in the QM/MM(PB) scheme than those of the QM/MM one. In particular, the average excitation energies in the QM/MM(PB) trajectories differ by only 0.04 eV. The same simulations have been carried out on a third TO-environment configuration extracted from the long classical MD run, giving results very close to MD2 both for the QM/MM and QM/MM(PB) schemes. Excitation energies and transition dipole moments have been reported in the ESI (Section S5†).
QM/MM | QM/MM(PB) | movE | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
〈Eexc〉 | |μ|2 | 〈BLA〉 | C corr | 〈Eexc〉 | |μ|2 | 〈BLA〉 | C corr | 〈Eexc〉 | |μ|2 | 〈BLA〉 | C corr | |
MD1 | 3.02 | 13.1 | −2.5 | −0.61 | 2.91 | 9.9 | −2.7 | −0.54 | 2.82 | 10.8 | 0.3 | −0.32 |
MD2 | 2.83 | 14.3 | −0.5 | −0.59 | 2.87 | 10.7 | −0.8 | −0.23 | 2.85 | 11.1 | 0.2 | 0.11 |
MD1–MD2 | 0.19 | −1.2 | −2.0 | 0.04 | 0.8 | −2.1 | −0.03 | −0.03 | 0.1 |
This analysis, based on the two sets of trajectories MD1 and MD2 and fostered by the third set of simulations reported in Section S5 of the ESI,† indicates that, as expected, allowing the local environment due to the closest NBs to move together with TO reduces the differences introduced by the different initial configurations of the environment.
If we now want to estimate the contribution of the QM–MM partition, separating it from the structural effect related to the NB geometry relaxation, a different comparison has to be done. To do that, we selected a single TO-environment configuration (namely one of the two initial configurations extracted from the classical MD trajectory), and we use it in combination with different embeddings: (i) no environment, defined as vacuum, (ii) the QM–MM partition scheme, where only TO is QM, (iii) the QM/MM(PB) scheme, where also the four closest NBs are treated at the QM level and (iv) substituting the PB approach with link atoms, namely the QM/MM capping H. As usually done when the link atom is used, the electrostatic parameters of the classical MM atom that has been substituted by the hydrogen atom have been distributed between the classical atoms connected to it. The hydrogen link-atom has been compared here with the PB scheme only for single point calculations, to provide a wider range of methodologies to the reader, also considering the large diffusion of the link-atom approach.3,4
The excitation energy, together with the NTOs of the corresponding states, has been reported in Fig. 5. Adding an environment on the selected structure has generally a non-negligible effect on the excitation properties. The 3 eV excitation energy found for the vacuum case is red shifted by roughly the same amount (0.15 eV) in the three different embeddings. As a result, we can observe how the QM–MM partition itself does not significantly affect the excitation properties of TO. This observation is confirmed by the comparisons of the effects of different environments, i.e. for instance the (ii) and (iii) embeddings, on the many structures extracted from the MD1 (see Section S7 of ESI†). Additionally, the same analysis has been performed modeling the classical environment with the AMBER99sb61 force field for DNA and the TIP3P model for water,62 finding that within a classical, non-polarizable description of the environment, the computed properties are more sensitive to the definition of the QM and MM subsystems (see Section S8 of ESI†).
A third set of simulations has been performed, where also the environment has been allowed to move together with the QM subsystem (moving environment scheme, movE). In this test, a selection of the water solvent has also been performed: two selection radii are tested, first including all the solvent molecules (and the counterions) within a radius of 15 Å from each TO atom, and a second selection with a radius of 23 Å. All classical atoms (roughly 4200 for the small radius and 8400 for the larger one) are moving and treated as polarizable with the AMOEBA potential. As no significant differences in the computed excitation energies have been found between the two selection radii (<10−2 eV), only the results on the larger radius simulations are shown in Table 2, while the comparison between the two tested radii is reported in the ESI (Section S6†). As observed for the QM/MM(PB) simulations, also in this case the relaxation of the environment has a larger effect on MD1, with a final red shift in the averaged excitation energy of 0.09 eV if the QM/MM(PB) value is used for comparison, and 0.2 eV with respect to that in the QM/MM scheme.
As a key parameter to describe the relation between the internal structure of the TO dye and its excitation properties, the bond length alternation (BLA) defined in ref. 58, and highlighted in Fig. 8 with the TO structure, has also been computed on the structures obtained from all six hybrid trajectories. Their average values are reported in Table 2. Also the correlation coefficients computed between the BLAs and the excitations computed on the sampled snapshots are reported. Negative correlation coefficients highlight an anticorrelation between the two quantities, which can also be observed from the time series reported in Fig. 7. Particularly for the QM/MM simulations a non-negligible anticorrelation can be observed, with a coefficient larger in absolute value than 0.5. For MD2 the correlation coefficient becomes significantly smaller also in the QM/MM(PB) scheme. The MD1 and MD2 BLA values start to be very close to each other only in the movE. In this case, in fact, the constraints due to the fixed environment, together with the differences in the system configurations related to the initial structures are reduced.
Fig. 7 TO ππ* excitation energies (red lines) and BLA (blue lines) fluctuations along the QM/MM, QM/MM(PB) and movE simulation schemes adopted for the MD1 and MD2 trajectories. On the x-axis is reported the snapshot number while on the y-axis is reported the difference between the considered properties at the actual snapshot and the average value of the property (either the excitation energy or BLA) along the trajectory. All average values are reported in Table 2. |
Fig. 8 Molecular structure of the TO dye, represented here in licorice style, with the atoms involved in the BLA definition highlighted in orange. |
Structural effects have been studied on the QM/MM dynamics of oligopeptides, partitioning the aminoacid sequences into QM and MM portions, linked by covalent bonds. The effects of different QM–MM partitioning schemes have been investigated on the dynamics and the resulting electronic excitation of a dye embedded in a DNA fragment. The comparison shows that an extended QM subsystem going beyond the dye has small effects on the excitation properties if the structure of the whole system is kept the same. This indicates that very large QM subsystems are not required if accurate polarizable force fields such as AMOEBA are used. We also show that a non-polarizable description of the environment using a classical force field may lead to computed properties which are more sensitive to the definition of the QM and MM subsystems. Overall, the importance of a proper handling of the coupled dynamics of the QM subsystem and the environment is also highlighted. If a discussion on the computational cost of such treatment has not been explicitly included in this work, we can add, as a final remark, that (i) Molecular Dynamics with polarizable force fields using distributed multipoles such as AMOEBA can be achieved efficiently on very large systems using modern algorithmics;32,63,64 (ii) the proper treatment of the polarizable environment represents only a small overhead of a few percent on the already expensive SCF procedure of the QM subsystem that dominates the total cost of the simulation.
These results demonstrate that large scale hybrid QM/MM molecular dynamics simulations of complex systems using advanced point dipole polarizable embeddings are now possible offering enhanced accuracy and new perspectives of applications in various fields. Finally, this work paves the way for a production implementation of a polarizable QM/MM scheme to appear in the distribution versions of Gaussian and Tinker/Tinker-HP, for which a more detailed technical description will be given in order to detail the general performance of the model.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9sc01745c |
This journal is © The Royal Society of Chemistry 2019 |