Classical Force Fields Tailored for QM Applications: Is It Really a

Subscriber access provided by UNIV OF CAMBRIDGE

Article

CLASSICAL FORCE FIELDS TAILORED FOR QM APPLICATIONS: IS IT REALLY A FEASIBLE STRATEGY? Oliviero Andreussi, Ingrid Guarnetti Prandi, Marco Campetella, Giacomo Prampolini, and Benedetta Mennucci J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00777 • Publication Date (Web): 14 Sep 2017 Downloaded from http://pubs.acs.org on September 15, 2017

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

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

Page 1 of 35

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

Journal of Chemical Theory and Computation

CLASSICAL FORCE FIELDS TAILORED FOR QM APPLICATIONS: IS IT REALLY A FEASIBLE STRATEGY? Oliviero Andreussi1,2*, Ingrid G. Prandi3*, Marco Campetella3, Giacomo Prampolini4, Benedetta Mennucci3 1 Institute of Computational Science, Università della Svizzera Italiana, Via Giuseppe Buffi 13, CH-6904 Lugano, Switzerland 2 Theory and Simulations of Materials (THEOS) and National Centre for Computational Design and Discovery of Novel Materials (MARVEL), École Polytechnique Fédérale de Lausanne, Station 12, CH-1015 Lausanne, Switzerland 3 Department of Chemistry, University of Pisa, Via Giuseppe Moruzzi 3, I-56124 Pisa, Italy 4 CNR, UOS Pisa, Istituto di Chimica dei Composti Organometallici ICCOM CNR, Area Ric, Via Giuseppe Moruzzi 1, I-56124 Pisa, Italy * The first and second authors equally contributed to this work.

ABSTRACT Classical molecular dynamics is more and more often coupled to quantum mechanical based techniques as a statistical tool to sample configurations of molecular systems embedded in complex environments. Nonetheless, the classical potentials describing the molecular systems are seldom parameterized to reproduce electronic processes, such as electronic excitations, which are instead very sensitive to the underlining description of the molecular structure. Here we analyze the challenging case of the peridinin molecule, a natural apocarotenoid responsible for the light-harvesting process in the PCP antenna protein of dinoflagellates. Ground state structural and vibrational properties, as well as electronic transitions of the pigment are studied by means of quantum-mechanical static and dynamic calculations. Thereafter, classical molecular dynamics simulations are performed with a number of different force-fields, ranging from a popular, general purpose one to refined potentials of increasing level of complexity. From the comparison of classical results with their quantum mechanical counterparts, it appears that, while very poor results are obtained from standard transferrable force-fields, specifically tuned potentials are able to correctly characterize most of the structural and vibrational features of the pigment. Nonetheless, only an advanced parameterization technique is able to give a semiquantitative description of the coupling between vibrations and electronic excitations, thus suggesting that the use of classical MD in combination of QM calculations for the study of photoinduced processes, albeit possible, should be considered with care.

1. INTRODUCTION In the past years, hybrid computational methods, which combine quantum mechanical (QM) descriptions with classical models, have become standard techniques to study properties and processes of molecular systems embedded in different environments. While for homogeneous solvents, continuum models are successfully used,1,2 for more complex environments where anisotropies and heterogeneities become important, atomistic models such as those based on Molecular Mechanics (MM) are to be preferred.3–6 The resulting

ACS Paragon Plus Environment

1


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

Page 2 of 35

QM/MM formulations however require to be used in combination with statistical tools allowing for a proper sampling of the configurational space of the full system. In most cases, a sequential classical-QM strategy7,8 is used, where sampling is performed by means of Molecular Dynamics (MD) at a classical level, due to the large dimensions of the system of interest and/or the long time window that has to be explored. As a result, the quality of the final QM description strongly relies on the accuracy through which the initial classical MD describes both intra- and inter- molecular interactions. The QM/MM calculations of the property of interest are in fact performed on configurations extracted from the classical MD trajectory, generally without any further refinement. If the internal structure of the QM subsystem and/or the relative position and orientations of the QM and the MM subsystems are not accurately predicted, the QM/MM calculations will lead to artifacts.9 Even in the absence of an embedding medium, the most delicate aspect in a sequential MD-QM strategy is to achieve a reliable description of the intra-molecular degrees of freedom of the investigated molecule. As a matter of fact, a not accurate description of the geometry of the molecule and its vibrational modes is a real risk when photo-induced processes are investigated, since the QM calculations will be largely affected by the way the classically described nuclear degrees of freedom couple with the electronic excitations. An effective way of quantifying these couplings is through the frequency dependent spectral density (SD) function10. If we assume a linear coupling of the electronic degrees of freedom and the nuclear modes and treat them as harmonic oscillators, the spectral density can be computed using a Fourier transform of the excitation energy autocorrelation function, which can be in turn computed straightforwardly through the aforementioned sequential classical-QM approach.11–15 In such application, the nuclear degrees of freedom are treated using the classical MD whereas the excitation energies are calculated quantum mechanically along the classical trajectory. It is evident that the quality of the results will be almost completely determined by the quality of the force field (FF) used in the MD.9 The simulation of the SD can thus become an extremely important validation property for assessing the quality of MM FF to be used in quantum chemical based approaches for photo-induced processes. In this work, a very challenging molecule, the apocarotenoid peridinin (PID), is used to investigate limits and potentials of this strategy. Peridinin represents the absorbing unit used in the peripheral water-soluble light-harvesting (LH) complex of most photosynthetic dinoflagellates, the peridinin-chlorophyll-complex (PCP), which constitute the main part of oceanic plankton. The uniqueness of PCP as LH antenna, is that carotenoid molecules (the PIDs) are preponderant over the chlorophylls, and they act as primary absorbing units. The key structural features of PID are a lactone and an allene group conjugated to the polyene chain, conferring special spectroscopic properties to the molecule16–25. These unique properties are exactly the reason for our choice. Indeed the presence of a polar group (the lactone) within the polyene chain makes the coupling of the vibrational modes with the electronic degrees of freedom extremely challenging for classical FF based simulations. Thereafter, classical MD simulations are performed with three different FFs, ranging from a popular, general purpose, one to more refined models, derived according to two different strategies from Density Functional Theory (DFT) data, purposely computed on the target chromophore. These two newly developed FFs, which represent extremely refined formulations specifically developed for PID, are here used as “optimal” FF descriptions. The


2

Page 3 of 35

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


results obtained from QM calculations over trajectories obtained with such refined FFs, when compared with fully QM descriptions, should give us an insight about the limits and potentials of the sequential classical-QM strategy. Here two fully QM descriptions are used as benchmarks, namely a static model based on a single QM optimized structure and a trajectory from an ab initio Born-Oppenheimer MD simulation. As the scope of this study is the analysis of the coupling between electronic and internal nuclear degrees of freedom arising form a classical description of the internal degrees of freedom, the whole investigation has been limited to peridinin in vacuum. From the comparison of classical results with their quantum mechanical counterparts, it appears that, while very poor results are obtained from standard transferrable FF, specifically tuned potentials are able to correctly characterize most of the structural and vibrational features of the pigment. Nonetheless, only an advanced parameterization technique is able to give a semiquantitative description of the coupling between vibrations and electronic excitations. The article is structured as follows: in Section 2 we summarize the different approaches exploited to generate and sample configurations of PID, followed by a detailed description of methods and parameterizations of the different classical Force Fields adopted; in Section 3 we report the results obtained for the different simulations in terms of the coupling between the lowest bright excitation of PID and its vibrational modes; this analysis is accompanied by an investigation of structural and vibrational properties.

2. METHODS AND PARAMETERIZATIONS 2.1 STRUCTURES AND TRAJECTORIES The following structures or ensembles of structures were considered in this work: • Structures obtained by geometry optimization of the eight different molecules present in the crystallographic structure of the PCP monomer (PID611, PID612, PID613, PID614, PID621, PID622, PID623, PID624, as extracted from the 1PPR structure16 found in the Protein Data Bank). Calculations were performed at the DFT level with the hybrid B3LYP functional and the Gaussian 6-311g(d,p) basis set, exploiting the Gaussian 09 program.26 In the following these results are reported with the label GAU-OPT. • Single structure obtained by geometry optimization of the first peridinin molecule (PID611) in the crystallographic structure of the PCP monomer. Also in this case, calculations were performed at the DFT B3LYP level of theory, but exploiting the CP2K package,27 using a mixed basis set composed by plane-waves and localized functions, consisting of plane-waves within a cutoff of 300 Ry and a double zeta plus polarization basis set (DZVP) with GTH pseudo-potentials to handle core electrons. A simulation cell of 30 Angstrom with the removal of periodic boundary artifacts was adopted. In the following these results are reported with the label CP2K-OPT. • Structures sampled from an ab initio Born-Oppenheimer MD (BOMD) simulation performed with the CP2K package27, starting from the crystallographic structure of PID611 and using the same parameters as reported above. Simulations in the canonical NVT ensemble were performed using the thermostat from Bussi and


3


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

•

•

•

Page 4 of 35

coworkers,28 with a target temperature of 300 K and a thermostat time constant of 1 ps. A simulation timestep of 0.5 fs was adopted, for a total simulation time of 32.5 ps. For the collection of ensemble averages, only configurations from the second half of the trajectory were considered, whereas the first half of the simulation was used for equilibration. For the calculation of the IR spectrum and of the spectral density, configurations were sampled every 2 fs, for a total of about 8000 structures. In the following these results are reported with the label CP2K-BOMD. Structures sampled from classical MD simulations performed with the Amber suite of programs and the general Amber force-field (GAFF)29. Atomic charges were computed using the RESP procedure with a HF/6-31G(d) level of theory, averaged over the first four peridinin structures (PID611-PID614), optimized at the B3LYP/631G(d) level of theory. After a preliminary minimization, simulations in the NVT ensemble using a Langevin thermostat with a target temperature of 300 K and a 1 ps relaxation time were performed in open boundary conditions. This choice of thermostat was dictated by the lack of better alternatives to enforce a true microcanonical ensemble. A timestep of 1 fs was used throughout, with an initial equilibration run of 1 ns followed by a production run of 5 ns, from which 1000 uncorrelated configurations were obtained by sampling every 5 ps. For the calculation of the IR spectrum and of the spectral density, a further simulation with a reduced timestep of 0.5 fs and a total of 20 ps of simulation time was performed, from which 10000 structures were extracted (one every 2 fs). In the following these results are reported with the label GAFF-CLMD Structures sampled from classical MD simulations performed with the Amber suite of programs30 and a specifically tailored force-field, developed starting from firstprinciples simulations performed on the first four peridinin molecules (PID611PID614) at the B3LYP/6-311G(d,p) level of theory. The Amber functional form of the force-field was adopted, and an atom-specific parameterization was used, similar to the strategy adopted by Prandi et al.31. Details of the simulations are the same as in the previous GAFF-CLMD case. In the following these results are reported with the label STAFF-CLMD. Structures sampled from classical MD simulations performed with the Gromacs software (version 4.6.1)32 and a tailored force-field generated with the help of the JOYCE program,33,34 applied to first-principles simulations performed at the B3LYP/6311G(d,p) level of theory on the minimal energy configuration obtained among the different peridinin molecules (corresponding to PID613). Simulations were performed in the canonical NVT ensemble, with a target temperature of 300K enforced by the velocity rescale thermostat of Bussi et al.35 with 0.1 ps time constant. A simulation timestep of 1 fs and no periodic boundary conditions were adopted. The hydrogen bonds were constrained with LINCS algorithm36. Following 5 ns of equilibration, a correlated 32 ps MD was performed and 8000 snapshots were extracted every 4 fs for the IR spectrum and spectral density analyses. For the uncorrelated MD we have continued the correlated simulation for other 5 ns and extracted 1000 frames from the MD trajectory (1 frame every 5 ps). In the following these results are reported with the label JOYFF-CLMD.


4

Page 5 of 35

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


For all the above structures, electronic excitation energies were computed through the Gaussian 09 program using the Tamm-Dancoff (TDA) formulation of the TD-DFT approach, coupled with the hybrid B3LYP density functional and a 6-311G(d,p) basis set. IR vibrational spectra were computed either using linear response calculations on optimized structures (GAU-OPT and CP2K-OPT) or from the autocorrelation of QM electric dipoles, as computed along the trajectories by means of DFT calculations at the B3LYP/6-311G(d,p) level of theory, still exploiting the Gaussian 09 program.

2.2 FORCE-FIELD PARAMETERIZATION In the standard Amber formulation, the generalized force-field (GAFF) is partitioned into bonded and non-bonded terms. The former are described through sums of bond stretching, angle bending and two types of dihedral torsions (proper and improper), with the following expression = ∑ − + ∑ − + ∑ ∑ , 1 + (1)

cos (%& − ', )) + ∑*++ & − & , with , , , and being the force constants; and are the bond length and bond angle respectively; and are the corresponding bond and angle equilibrium values; & corresponds to the dihedral angle, % is its periodicity number, ', is its phase, and & is the equilibrium value of the improper dihedral angle.45 Non-bonded terms are instead modeled as a sum of pair-wise interactions, namely longrange electrostatic Coulomb and short-range Lennard-Jones (LJ) interactions expressed as ,-

:

2 / /0 1 30 =. 2 + . 460 7 8 9 1 30 4 30 40

40

Classical Force Fields Tailored for QM Applications: Is It Really a

Recommend Documents