Simulating Electron Dynamics in Polarizable ... - ACS Publications

Subscriber access provided by Hong Kong University of Science and Technology Library

Article

Simulating Electron Dynamics in Polarizable Environments Xiaojing Wu, Jean-Marie Teuler, Fabien Cailliez, Carine Clavaguera, Dennis R. Salahub, and Aurélien de la Lande J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00251 • Publication Date (Web): 24 Jul 2017 Downloaded from http://pubs.acs.org on July 27, 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 48

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

SIMULATING ELECTRON DYNAMICS IN POLARIZABLE ENVIRONMENTS

Xiaojing Wu1, Jean-Marie Teuler1, Fabien Cailliez1, Carine Clavaguéra1, Dennis R. Salahub2,3, Aurélien de la Lande1* 1: Laboratoire de Chimie Physique, CNRS - Université Paris Sud, Université Paris-Saclay. 15 avenue Jean Perrin, 91405, Orsay CEDEX, France. 2: Department of Chemistry, Centre for Molecular Simulation, Institute for Quantum Science and Technology and Quantum Alberta, University of Calgary, 2500 University Drive N. W., Calgary, Alberta, Canada T2N 1N4. 3: College of Chemistry and Chemical Engineering, Henan University of Technology, No 100, Lian Hua Street, High-Tech Development zone, Zhengzhou 450001, P. R. China *To whom correspondence should be addressed: [email protected] ABSTRACT We propose a methodology for simulating attosecond electron dynamics in large molecular systems. Our approach is based on the combination of Real-Time Time-Dependent-Density-Functional-Theory (RTTDDFT) and polarizable Molecular Mechanics (MMpol) with the point-charge-dipole model of electrostatic induction. We implemented this methodology in the software deMon2k that relies heavily on auxiliary fitted densities. In the context of RT-TDDFT/MMpol simulations, fitted densities allow the cost of the calculations to be reduced drastically on three fronts i) the Kohn-Sham potential, ii) the electric field created by the (fluctuating) electron cloud which is needed in the QM/MM interaction, and iii) the analysis of the fluctuating electron density on-the-fly. We determine conditions under which fitted densities can be used without jeopardizing the reliability of the simulations. Very encouraging results are found both for stationary and time-dependent calculations. We report absorption spectra of a dye molecule in the gas phase, in non-polarizable water and in polarizable water. Finally, we use the method to analyze the distance-dependent response of the environment of a peptide perturbed by an electric field. Different response mechanisms are identified. It is shown that the induction on MM sites allows excess energy to dissipate from the QM region to the environment. In this regard, the first hydration shell plays an essential role in absorbing energy. The herein presented methodology opens the

1 ACS Paragon Plus Environment


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 48

possibility of simulating radiation-induced electronic phenomena in complex and extended molecular systems. INTRODUCTION Recent years have seen a growing interest in the electron dynamics taking place in molecules when they are subjected to an external perturbation. This interest has been stimulated by progress in attosecond spectroscopy that now gives access to details on electron dynamics. The realm of sub-femtosecond electron dynamics involves fascinating processes such as ultrafast charge migration1, Auger decays and Intra Coulomb Decays2-5. These are not driven by nuclear dynamics but instead by electron correlation and energy redistribution6. This nascent research field has led to new debated concepts like attosecond chemistry, a possible promise of which would be the possibility to control chemical reactions by the control of electronic motion7. Electron dynamics is also important in the description of ultrafast nonadiabatic molecular dynamics. The relaxation pathways within molecules electronically excited or ionized by a photon or a high-energy particle are particularly rich and complex. They involve coupled electronnuclear dynamics8. On the computational side, much effort has been spent to devise simulation algorithms of electron dynamics. In the family of wave function approaches the TD-HF (Time-dependent-Hartree–Fock)9, TD-CI (Time-dependent-configuration

interaction)10-11

or

the

TD-MCSCF

(Time-dependent

Multi

Configurational Self Consistent Field)12-13 methods have been developed. Another popular approach for simulating electron dynamics relies on time-dependent Density Functional Theory (TDDFT). This approach is frequently referred to as Real-Time TDDFT (RT-TDDFT) to distinguish it from the Linear Response (LR-TDDFT) formalism. The latter relies on perturbation theory to simulate UV-visible absorption spectra14. Although not exempt from intrinsic limitations like the self-interaction-error15-16, a noticeable advantage of TDDFT is its excellent computational cost/accuracy ratio. TDDFT can be applied to molecular systems comprised of hundreds of atoms. TDDFT finds its root in the seminal work of Runge and Gross17. Under the Kohn-Sham framework that refers to a fictitious reference system of noninteracting electrons, the coupled time-dependent KS equations describe the time evolution of the KS molecular orbitals, hence the dynamics of the electron density of the real system. RT-TDDFT has been used to calculate static and dynamic polarizabilities and hyperpolarizabilities of molecules18, to simulate UV-visible spectra of molecules19-20 or of nanoparticles21, to simulate core-level near-edge X-ray absorption spectra22-23, to simulate photoelectron emission spectra24-25, electron conductance in electronic junctions26-27, photoinduced electron transfer28-29, magnetization dynamics in inorganic 2 ACS Paragon Plus Environment

Page 3 of 48

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


complexes30, attosecond dynamics following X-ray photoionization of gas molecules3 or charge migration following radiolysis of water5, 31. RT-TDDFT has also been coupled to mean-field (Ehrenfest) nuclear dynamics to simulate non-adiabatic processes32-34 with many interesting applications, for example to simulate the ultrafast dynamics of photoexcited metal complexes or in optimal control of chemical reactions8. Curiously enough, most implementations have been designed for molecular systems in the gas phase and not in contact with environments. If the systems of interest are periodic, periodic boundary conditions can be used to simulate infinite systems35. But in many cases, systems are not periodic and alternatives must be found. A challenge is to account for the electronic response of the environment due to the changes in the electronic structure of the molecule, and vice versa. The environment may be homogeneous (solutions) in which case a polarizable dielectric continuum (PCM) can be used. The environment may also be heterogeneous, as for example for extended biosystems (DNA, proteins, lipid membranes), nanoclusters or interfaces. In such cases hybrid QM/MMpol (i.e. using polarizable force fields)36 constitutes a method of choice to retain the atomistic details of the environment at moderate computational cost. In the linear response formalism, coupling between TDDFT and either PCM or polarizable QM/MMpol have been devised37-39. For explicit propagation in time of the TDDFT equations, a further challenge is to account, by definition, for the time-dependence of the environment’s response. Remote atoms should take longer times to respond than closer ones, for instance. Li and co-workers developed a combined RT-TDDFT/PCM40-41 method, with applications to charge transfer dynamics in bulk heterojunction models42. In their PCM model the dielectric constant of the environment was made timedependent even though the PCM was made stationary with the evolving potential created by the QM region. Corni et al.43-44 as well as Ding et al.45 later described a more general approach of RT-TDDFT/PCM calculations where both the PCM and QM region were propagated in time. Regarding hybrid RTTDDFT/MMpol approaches Dinh et al. reported a few years ago a coupling between RT-TDDFT and a polarizable force field (FF)46. Induction was introduced by distinguishing core from valence electrons on MM atoms. The average position of the core electrons + nucleus and the average position of valence electrons had the possibility to be different depending on the electrical environment, thereby creating induced dipoles. The authors reported insightful applications to sodium clusters deposited on metal surfaces46 with detailed analyses of their optical properties. In the computational set-up of Dinh et al. the RT-TDDFT engine itself relies on a grid-based implementation of DFT. This is quite different from the algorithms employed in the community of quantum chemists that generally rely on local basis sets



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 48

(Gaussian or Slater atomic orbitals). There is thus a clear need to develop hybrid RT-TDDFT/MMpol schemes for the modelling of electron dynamics with local basis sets in extended molecular systems. The structure of the present article is as follows. First, we report an implementation of polarizable QM/MMpol based on the charge-induced dipole model36 of electronic induction in the software deMon2k47. Second, we describe our implementation of RT-TDDFT and its coupling with polarizable MM. In both modules density fitting techniques are used to reduce the computational cost drastically48-49. In section II we carefully test the reliability of substituting the Kohn-Sham density by the auxiliary density for propagating the electron dynamics or for calculating the QM/MMpol coupling interactions. Very encouraging results are obtained. We investigate in section III the time-dependent electronic response of molecules in vacuum and in solution.

I METHODOLOGY I.1 Auxiliary Kohn-Sham Density Functional Theory We start the methodology section by recalling the general DFT framework implemented in deMon2k. This program solves the Kohn-Sham DFT equations with KS molecular orbitals (MO) represented as Linear Combinations of Gaussian-Type Atomic Orbitals (LCGTAO). For simplicity, we will consider only

closed-shell molecules, but we mention that the methodologies presented in this work have been adapted to the open-shell case too. =

, = ,

⁄

∗ = 2

(1)

(2)

(3)

Greek letters are used both as indexes and as AO function names. The MO coefficients ( , hence the

density matrix ( ) and the electron density (), depend on time. When solving the stationary KS equations, this dependence would by definition vanish. But we keep here the more general RT-TDDFT 4 ACS Paragon Plus Environment

Page 5 of 48

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


formulation. Note that the MO coefficients are complex numbers in real time propagation. deMon2k relies heavily on the variational density fitting method originally introduced by Dunlap48 to avoid the

calculation of four-center electron repulsion integrals (ERIS). The fitted densities () are expressed as linear combinations of auxiliary basis functions : = ∑ . For computational efficiency the

auxiliary basis functions are Hermite Gaussian polynomials that are grouped by functions sharing the same exponents50. With this auxiliary density, the electronic energy expression for an isolated molecule

reads:

1 " = # + %‖' − * %‖+' + ",- ./ 2 ,

,

,*

(4)

The symbol ∥ stands for the coulomb operator (1/|3 − 3 |). As evident in Eq. 4, no four-centers ERIS are

needed but only two- and three- centers ERIS. # are the matrix elements of the core Hamiltonian,

encompassing the kinetic energy and the electron-nuclei attraction. The Kohn-Sham potential matrix

elements for an isolated molecule are obtained by differentiating the total energy with respect to the density matrix elements: 45* # ≡

7" 7",- ./ = # + %‖' + 7 7

(5)

deMon2k offers also the possibility to use the fitted density in the calculation of the exchange-

correlation (XC) energy, in which case ",- ./ is replaced by ",- ./ in Eq. 4 and 551. Now the KS potential does not depend explicitly on the KS density but only on the auxiliary density. We refer to this framework as Auxiliary DFT52.

I.2 Polarizable QM/MM in deMon2k Model of electronic induction. There are various ways to carry out DFT/MM calculations with deMon2k and we refer the interested reader to a recent review describing these alternatives53. In the present work we focus on the so-called in-deMon2k QM/MM by which both QM, here (TD)DFT, and MM calculations are carried out by deMon2k, without using program interfaces53. We expect advantages in terms of data passing management between the DFT and MM modules. The objective of the present work is to upgrade the preexisting QM/MM method to QM/MMpol. Electrostatic induction can be introduced in classical force fields in different ways36, 54-55, 56. We note, for example, that the Drude polarizable force5 ACS Paragon Plus Environment


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

Page 6 of 48

field has recently shown great promise for ions interacting with protein models, for which additive fixedcharge force fields come up short57-58. We have chosen here to consider the point-charge dipole model by which induction is simulated by the introduction of induced dipoles (8 , note that vectors are written

in bold) on every polarizable MM site i36-37. Each induced dipole is determined from the electric field 9

at the MM atom position. 9 stems from the electric field created by other MM permanent charges :

(9 ) and by other MM induced dipoles (9;< ). In QM/MMpol calculations one further adds the electric

?=> field created by the QM region (9=> ) and by the electron ), that is by the atomic nuclear charges (9 @

density (9 ). We also introduced the possibility to add an external electric field (9A,B ) to mimic for

instance the interaction with the electric part of an electromagnetic wave. The mathematical expressions for the various contributions are given by: 8 = C 9 = C D9 9

:

=

G∈>> GJ

:

FG 3G H G

9;< = − KG L G∈>> GJ

1 3 KG = H M − O PQ 3G 3G R

+ 9;< + 9=> + 9A,B E

(7)

(8) Q Q RQ

9=> = 9?=> + 9 = @

(6)

∈=>

R QRS R

(9)

T − V3 − U H | − |H 3

(10)

C is the polarizability of MM atom W. It is assumed to be isotropic. FG is the charge of MM atom X; G is

the vector between atoms W and X; KG is the dipole-dipole interaction tensor and M is the identity matrix. T is the nuclear charge of QM nucleus . The total induction energy is comprised of three terms that

;< reflect the interaction between the MM induced dipoles with i) the MM permanent charges ("YZ[>> ), ii) ;< ;< the atomic nuclei of the QM atoms ("YZ?=> ), iii) the electron cloud ("YZ@ ). ;< ;< ;< ;< "B5B = "YZ[>> + "YZ?=> + "YZ@

(11)

: Contributions (i) and (ii) are calculated as −1⁄2 ∑∈>> 8 . 9 and −1⁄2 ∑∈>> 8 . 9?=> respectively.

The last term depends on the electron density and is given by


Page 7 of 48

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


;< "YZ@ =

1 1 8 − U 8 . − V3 = \] ]^ H | − | | − |H 2 2 ∈>>

(12)

∈>> ,

QM/MMpol calculations must capture the interdependence between the electron cloud and the polarizable environment. For stationary DFT calculations this is done by regularly updating the induced

and the electrostatic potential created by the induced dipoles to be dipoles that depend on 9=> included in the KS Hamiltonian. We use a similar algorithm for RT-TDDFT propagation. We will come back

;< to this methodological point at the end of Section I . Note that besides "YZ@ , the interaction between

the QM and MM atoms also includes the electrostatic energy between the electron cloud and the MM permanent charges (F ). "[Z@ = U _A`a

∈>>

1 F . V = F bc cd | − | | − |

(13)

∈>> ,

Within this electrostatic embedding scheme, the potential created by permanent charges and induced dipoles on MM sites are obtained by differentiation of the respective interaction energies with respect to

the electronic density. The calculation of electrostatic integrals (be ` ed) required in Eq. 13 was

optimized by Alvarez et al.53, 59. To take advantage of these algorithmic developments we represent each

induced dipole 8 by two charges of opposite sign (±g ) separated by 0.5 bohr and centered around the Y@

MM atom positions60-62. This way the potential created by the induced dipoles (# ) are included in the

KS potential via a set of point charges, the calculation of which is performed efficiently in deMon2k53, 59. Y@ #

;< 7"YZ@ 8 ≡ = bc H cd | 7 − | ∈>>

= hg ij ∈>>

1

c − + 0.5

8 c ‖L ‖

jm − g ij

1

c − − 0.5

8 c ‖L ‖

jmn

(14)

As for any QM/MM scheme, a critical point of DFT/MMpol and RT-TDDFT/MMpol calculations is to set the boundaries between MM and QM regions. Setting boundaries across polar groups may deteriorate the efficiency of the hybrid energies and the derived properties. The choice of QM/MM portioning is the responsibility of the user. Our QM/MMpol implementation doesn't rely on interfaces between QM and 7 ACS Paragon Plus Environment


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

Page 8 of 48

MM software and both DFT and MM calculations are done within deMon2k. This a critical advantage to reach efficient DFT/MMpol calculations. Indeed, passing information between the DFT and MM branches of the same program is rapid compared to I/O operations. The electric field and the induced dipoles are vectors that can be stored in RAM (Random Access Memory), and one can easily restart convergence of induced dipoles at every new SCF cycle or RT-TDDFT time step from previous steps. Finally, our indeMon2k QM/MMpol uses a direct algorithm similar to the calculation of electron repulsion integrals (ERIS)50, 63 . At every SCF cycle or RT-TDDFT step o, instead of calculating the QM electric field 9 from the current density matrix ; , we increment it from 9

∆ = ; − ;Z , namely 9

=>,;

= 9

=>,;Z

=>,;Z

− ∑Y,r ΔYr bLe

etd. `s

=>,;

and the difference density This procedure has the

advantage that one can screen many terms of the sum if ΔYr is below a predefined threshold. This helps to decrease the computational time in SCF calculations when the density is close to convergence or in

RT-TDDFT simulations when the density evolves slowly. Direct SCF procedures have been used by other groups in the context of QM/MMpol calculations64-65. Polarization catastrophe. A well-known pitfall of polarizable force fields is the risk of "polarization catastrophe" the origin of which has been exposed by Thole56, 66-67. This term defines a divergence of the polarization energy that happens when adjacent dipoles align on the same line in head-to-tail configurations. Most polarizable MM implementations avoid the polarization catastrophe by damping the electric fields at short distance. For MM atoms bonded in 1-2, 1-3 or 1-4 positions to a given polarizable MM site, the electronic field can be simply ignored. In deMon2k, choice is given to the user to set-up these parameters. For non-bonded atoms electric damping is achieved by the modification of the dipole interaction tensor with two distance-dependent screening functions (uA and uB ):

uA 3uB KG = H M − O PQ 3G 3G R

Q Q RQ

R QRS R

(15)

Three alternatives of the screening function have been implemented in deMon2k following previous proposals reported in the literature. One is the linear scheme67: v = 3G ⁄w with w = xyC CG z uA = |

⁄{

1.0 Wu 3G > w 4v − 3v Wu 3G < w H

(16) (17)


Page 9 of 48

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


uB = |

1.0 Wu 3G > w v Wu 3G < w

(18)

another is the exponential scheme: v = x3G ⁄yC CG z uA = 1 − uB = 1 −

⁄{

(19)

v + v + 1 exp −v 2

vH v + + v + 1 exp −v 6 2

(20) (21)

while the third one is the Tinker-exponential form: = 3G ⁄yC CG z

⁄{

uA = 1 − exp −xH

uB = 1 − 1 + xH exp −xH

(22) (23) (24)

In these expressions x is a unitless parameter that depends on the force field. The higher the x the faster the field damping with distance. In the water model of the AMOEBA force field x is for example set to 0.39.68

On the use of fitted densities. We have tested various alternatives for estimating 9 . The most correct

way is to calculate it from the Kohn-Sham density (Eq. 2 and 10). Since 9

=>

=>

needs to be evaluated on

every MM site at every SCF cycle and at every RT-TDDFT step, this task can become computationally

expensive. A tempting alternative is to replace by the auxiliary density . Because the number of

auxiliary basis functions is typically four to five times lower than the number of products of atomic basis

functions, substituting by is expected to drastically reduce the cost of calculation of 9=> . Such a

substitution is however not guaranteed to yield reliable results. This is because auxiliary fitted densities are not designed to reproduce but to provide auxiliary densities from which approximate electronic

repulsion interactions can be computed with reduced computational cost. That said, we recently showed that electrostatic multipoles obtained either from the KS density or from the fitted density are very similar, provided sufficiently flexible auxiliary basis sets are used69. This is an encouraging result. Indeed,

we may expect that if the intrinsic multipoles on QM atoms extracted from are similar to those 9

ACS Paragon Plus Environment


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

Page 10 of 48

extracted from , so will be the electric fields generated by and . We will thus test the accuracy of stationary and time-dependent DFT/MMpol calculations when replacing by .

I.3 Electron dynamics equations-of-motion We now move to the description of our RT-TDDFT implementation. It is largely based on algorithmic developments reported previously by other groups in the last decade9, 20, 26, 70. We thus refer the reader to the original publications. We insist here on the specificities of the implementation in deMon2k and on the novel features that we have introduced, notably the coupling between RT-TDDFT and the QM/MMpol just described. Runge and Gross developed the many body wave function TD Schrödinger equation into the single-particle TD density Kohn-Sham (TDKS) equation with an effective Hamiltonian #t uniquely described by the TD electron density .17 W

7 = #. / 7

(25)

where # is the time-dependent Kohn-Sham operator which is a functional of the charge density. It includes the KS potential of the isolated molecule the matrix elements of which are given by Eq. 5 and

the interaction potentials of the electron cloud with external electric fields # A,B . In QM/MMpol the KS [@

operator also includes the perturbation from the MM charges (# ) and MM induced dipoles (# # = # 45* + #

[@

+#

Y@

+ # A,B

Y@

).

(26)

# 45* includes the contribution from the XC potential. We make the adiabatic approximation and

consider only the spatial dependence of the XC potential, neglecting its temporal non-locality. Eq. 25 can

be recast using the density matrix into a Liouville-von Neumann type of equation, which reads in the case of a non-orthogonal basis set: W

7 = # − # 7

(27)

is the overlap matrix in the atomic orbital basis set. As in Ref. 9 and 20 we transform # and to the

orthogonal MO basis leading to # and . Since the MO are orthogonal Eq. 27 simplifies to:


Page 11 of 48

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

W


7 = .# , ′ / 7

(28)

We will use primes to denote matrices in the molecular orbital (MO) basis and no primes to denote matrices in the atomic orbital (AO) basis. The formal solution of Eq. 28 can be expressed as:

= , : : , :

(29)

Where U is the evolution operator, which can be discretized into small time steps ∆ , ;Z

, : = + Δ ,

BB

+ ∆ , = |−W U B

(30) # ′ V

(31)

is the time-ordering operator, ensuring that operators associated with later times always appear to

the left of those associated with earlier times. Many schemes have been proposed to evaluate the propagator in RT-TDDFT and we refer the reader to recent reviews describing the physical conditions that propagators should fulfill70-71. We have implemented in deMon2k the Euler and second-order Magnus propagators. Euler propagation. To solve Eq. 28 by applying Lagrange's Mean Value Theorem we obtain: + ∆ = − W.# , / ∗ ∆

(32)

The propagation of the density matrix requires only the value of the density matrix and Kohn-Sham matrix at the current time. These are easy to obtain; however, this propagation scheme doesn't guarantee the preservation of the norm of the KS wave-function which can lead to divergence of the electronic propagation. We found the Euler propagation to be unstable in most of our applications and it will not be considered any further in this article. Magnus propagation. A convenient solution to Eq. 31 is given by a Magnus expansion72: B∆B

|−W U B

# V = = ⋯

(33)

Where Ω is a series of nested commutator integrals: 11 ACS Paragon Plus Environment


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

BB

Ω + Δ , = −W U B

BB

Ω + Δ , = −W U B

# V

(34)

V U V #′ , # ′ B

Page 12 of 48

(35)

Stopping at second order: W = Ω , this integral can be evaluated using a quadrature formula26: Ω + Δ , ≃ −W# ¡ +

Δ

¢ ∗ Δ

2

(36)

This is equivalent to the well-known split-operator method70.

Three algorithms have been implemented to calculate the matrix exponential entering eq. 33. The first one is based on the diagonalization of the £ matrix:

= ¤ , with £ = ¥

(37)

another is based on a Taylor expansion of the exponential:

1 = £ ; o!

(38)

; :

while the third one is based on the Baker-Campbell-Hausdorff (BCH) scheme20, 73. + Δ = +

1 1 1 .£, / + W, .W, P’t/ + ©£, £, .£, /ª + ⋯ 1! 2! 3!

(39)

Note that the latter scheme assumes the Kohn-Sham matrix is Hermitian. Other methods based on polynomial Chebychev expansion or Krylov subspace projections have been considered by other groups to evaluate the matrix exponents70, 74. Propagation with the Magnus scheme requires the knowledge of

the KS matrix at later time #′ + iterative algorithm70 #′ +

B

B ,

which is unknown. Two methods have been implemented. In the

is first extrapolated from the knowledge of #′ at earlier times. is then

propagated from to + Δ by the Magnus propagator and the resulting density matrix is used to build

the KS potential at + Δ . A new KS potential at time +

B

is interpolated from the potential at and

+ Δ . The propagation of is repeated with this new KS potential. The iterations are continued until

convergence. This is a robust but time-consuming procedure. An alternative is the two-step predictorcorrector scheme proposed by Van Voorhis and co-workers26.


Page 13 of 48

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


Application of an external electric field. We continue this section with the mathematical definitions of the external electric fields that can be applied in deMon2k. One option is to apply a Gaussian shaped pulse:

³ 9 = «a¬, . exp.− − : /2 / cos ± ²

(40)

³ = ´, Q´, R̂ is the Where : is the center of the pulse, is the pulse width, ± the field pulsation, ²

polarization vector, and «a¬, is the maximum field strength.

³ < : «a¬, cos± ² 9 = ¶ : ³ «a¬, cos± ² ≥ :

(41)

A drawback of the Gaussian shaped pulse or of the linear ramp is the possible introduction of spurious ¹º» static field effects which arise if the zero pulse area condition (ZPAC,¸B : 9 = 0) is not fulfilled 75-76.

B

Care must be taken to avoid such effects, for example by setting the center of the Gaussian pulse

sufficiently far from the initial time. In fact Gaussian pulses are not convenient in practical applications because of the shallow decay of Gaussian functions which require long simulation times to ensure the ZPAC. To alleviate this inconvenience some authors proposed to use squared sinusoidal functions75, 77. We implemented the following one in deMon2k:

¿

³ 9 = «a¬, . sin ¾ Â . Ã . ÃyÀ_Á*4A − z. cos ± ² À_Á*4A

(42)

where À_Á*4A is the duration of the pulse and Ã is the Heavyside function. Finally we also implemented

an infinitely narrow kick in the first step of the RT-TDDFT simulation. As illustrated below kick perturbations are useful to simulate absorption spectra. The applied field excites the molecular system through the coupling with the electrostatic dipole. The corresponding potential term is added to the KS matrix. " ¬__ = −8Ä ∙ 9Ä

8Ä = TÆ ÇÆ − %||' Æ

,

(43) (44)



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

A,B # ≡

7" A,B = %||'9 7

Page 14 of 48

(45)

Several analysis tools of the electron dynamics have been implemented. They will be introduced in the following applications sections as needed. On the coupling between the electron and the induced dipole intrinsic dynamics. The full hybrid QM/MMpol Hamiltonian is given by Eq. 46. #

=>/>>

= # => + # AaÈA< + # >>

(46)

where # => is given by eq. 5, # AaÈA< collects the coupling between the QM and MM region (embedding energy) which including the interaction between the electrons and QM nuclei with the permanent and

induced dipoles on MM sites. Finally # >> is the energy of the MM part computed with molecular

mechanics force field. The latter term also hold a dependence in time because of the interaction

between the time-dependent MM induced dipoles and the MM permanent charges. If one is interested in solving the time-independent KS equations to determine the stationary states of the system of interest, a common procedure is to relax the MM induced dipoles at every SCF cycle. The MM dipoles are then injected in the next SCF cycle to calculate a new embedding potential. The convergence threshold for converging the MM dipole moments is tightened along with the SCF convergence to reach, at global convergence, a user-defined value, typically 10-8 to 10-10 D. On the other hand if one is interested in the time-dependent solutions of the KS equations more subtle algorithms are needed because of the time dependence of each terms of eq. 46. In principle one needs to set up the coupled equations of motion for the overall system. This is not a trivial task because of the composite quantumclassical nature of the system. One may think of coupling RT-TDDFT for the electron cloud to a fictitious dynamics of the MM induced dipoles, in the spirit of what is done for molecular dynamics simulations with MMpol78-80. Here we consider a simpler scheme in which we make the assumption that the MM dipoles completely relax at each RT-TDFT step. In other words, we look for the stationary polarization state of the environment along with the non-stationary propagation of the electron cloud. For sufficiently small time steps this mixed stationary-non-stationary scheme is certainly a valid


Page 15 of 48

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


approximation to the real dynamics. As shown in the validation section we found this approximation to be acceptable for standard electron dynamics simulations with 1as time steps or below.

II PERFORMANCE AND VALIDATION II.1 RT-TDDFT propagation with density fitting As explained in the methodology section our implementation of RT-TDDFT relies on the use of fitted densities. In particular, one has the choice in deMon2k of using either the KS or the fitted density to calculate the (time-dependent) XC energies and potentials. These two approaches are referred to as BASIS and AUXIS respectively. Using the fitted densities usually induces a reduction of the computational

cost by a factor of 10, which is clearly advantageous. Yet it remains to be tested if can be used safely in RT-TDDFT propagation. As test cases, we consider two molecules: carbon monoxide and cysteine, a sulphur containing amino acid taken in the non-zwitterionic form. The propagation has been run using the Magnus propagator with an integration time step of 1as and diagonalization (Eq. 37). The simulations have been carried out with the TZVP-FIP281 basis set and the PBE functional82. Auxiliary basis sets are generated by an automatic procedure implemented in deMon2k that depends on the atomic orbital basis set. The GEN-An auxiliary function sets contain groups of auxiliary functions with s and spd angular momenta. The index n determines the number of auxiliary function sets, i.e. the number of these sets increases with increasing n. We have considered the GEN-A2 and GEN-A3 auxiliary function sets, as well as the GEN-A2* and GEN-A3* that are supplemented by f and g auxiliary functions. As a general rule of thumb the larger the auxiliary basis set the more accurate the DFT based energies and properties. An adaptive grid of accuracy 10-7 Ha has been used to integrate the XC potential and energies51. In Figure 1 we report the fluctuations of the x-component of the dipole moment of the two molecules when subjected to a constant electric field of intensity 0.01 a.u along the z-axis. Similar conclusions looking at the y- or z- components of the induced dipoles can be drawn. On the left-hand-side of Figure 1 we analyze the sensitivity of the propagation to the electron density used to calculate the XC potential, using the GEN-A2* auxiliary basis set. During the 10 fs of the propagation we find no important differences between the BASIS (orange) and AUXIS (blue) approaches. The simulation corresponding to the graphs on the right hand side have been obtained with the AUXIS approach, but with different auxiliary basis sets. For CO2 all the simulations give similar electronic 15 ACS Paragon Plus Environment


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

Page 16 of 48

evolution. The electronic response of CO2 is well captured by each auxiliary basis set. For cysteine the results are clearly more contrasted. Taking GEN-A3* as the reference auxiliary basis set, we find that the simulations with GEN-A2 and GEN-A3 are clearly different. On the other hand with GEN-A2* the simulation nicely reproduces the evolution of the induced dipole moment. This is an encouraging result. Provided sufficiently flexible auxiliary basis sets are chosen, one can rely on density fitting techniques to reduce the computational cost of the RT-TDDFT propagations similar to what is done in stationary auxiliary DFT or auxiliary perturbation theory calculations. Our RT-TDDFT implementation thus takes advantage fully of the optimized density fitting algorithms already implemented in deMon2k. We will come back to the code performance at the end of section II.

Figure 1

II.2 RT-TDDFT to calculate static polarizabilities We continue the validation section with the calculation of polarizabilities of molecules. Static and dynamic polarizabilities, as well as hyper-polarizabilities, can be calculated with standard DFT either by finite-field methods83-84, by the coupled-perturbed KS approach85-86 or by the auxiliary density perturbation theory87-89. These are actually recommended approaches for computing these properties at modest computational cost. Here they are used to test the validity of our RT-TDDFT module. After converging the stationary ground state of the molecule in the absence of an external field, an electric

field 9A,B is applied and the response of the electron density is simulated by RT-TDDFT. The resulting

induced dipole is related to the applied electric field vector via the polarizability tensor. The external electric field may either be constant or time-dependent, giving access to static or dynamic polarizabilities respectively. Focusing here on the static case, the polarizability tensor elements are given by83:

2 1 CG «G = L̅ y«G z − L̅ y−«G z − L̅ y2«G z − L̅ y−2«G z 3 12

(47)

where L̅ y«G z denotes the average contribution along W= , Q, R of the induced dipole when an external

field «G has been applied along X. The uncertainty in the knowledge of L̅ is given by L̅ = ⁄ÊËAÌÌ

where is the standard deviation of the sample of dipole moments calculated along the RT-TDDFT propagation. Ë is the number of terms in the sample and AÌÌ is the statistical chain efficiency. The

latter has been evaluated with the Coda package of the R project for statistical computing90-91. The


Page 17 of 48

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


calculations have been done with the PBE functional82, the TZVP-FIP1 basis set81 which has been optimized for electric properties calculations and the GEN-A2* auxiliary basis set. Fields «G of 0.01 a.u.

have been applied. This value is also used by default in the finite field method in deMon2k84. The RT-

TDDFT propagations have been run for 10 fs with an integration time step of 1 as. The propagatorcorrector Magnus scheme with diagonalization has been used. As can be seen in Table 1 the static polarizability tensors calculated from simulation of the electron density subjected to a perturbation by RT-TDDFT match nicely those obtained by the finite field method. The agreement between both approaches validates our implementation. 3

Table 1: Static polarizability tensors (bohr ) computed with RT-TDDFT and with a finite field method.

Molecule BENZENE

PHENOL

CYSTEINE

Finite field difference

RT-TDDFT

CG

X

Y

Z

X

Y

Z

X

83.5

0.0

0.0

83.6 ± 0.6

0.0 ± 0.0

0.0 ± 0.0

Y

0.0

83.5

0.0

0.0 ± 0.6

83.6 ± 0.6

0.0 ± 0.0

Z

0.0

0.000

44.5

0.0 ± 0.0

0.0 ± 0.6

44.3 ± 0.1

X

97.1

0.0

1.9

98.1 ± 0.4

0.0 ± 0.6

2.3 ± 0.0

Y

0.0

46.8

0.0

0.0 ± 0.0

46.8 ± 0.0

0.0 ± 0.0

Z

2.0

0.0

87.1

2.0 ± 0.1

0.0 ± 0.0

87.2 ± 0.2

X

78.7

13.3

-8.0

78.6 ± 0.0

13.4 ± 0.0

-8.1 ± 0.0

Y

13.3

72.6

-2.4

13.4 ± 0.0

72.4 ± 0.0

-2.4 ± 0.0

Z

-8.0

-2.4

85.2

-8.0 ± 0.0

-2.3 ± 0.0

85.0 ± 0.0

II.3 Absorption spectra in the gas phase A further validation of our RT-TDDFT implementation is now sought by comparing molecular absorption spectra calculated by LR- and by RT-TDDFT. Indeed, the electronic spectrum of a molecule is encoded in the evolution of the molecular dipole simulated by RT-TDDFT after a molecule is perturbed by an 17 ACS Paragon Plus Environment


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

Page 18 of 48

infinitely narrow electric field. The Fourier transform of the dipole signal gives access to the polarizability tensor in the frequency domain which in turn yields the dipole strength function (Absorption spectrum). More details about this procedure can be found in many recent publications19-21. Coumarin has been chosen here as a test case because it is a solvatochromic dye, the dipole moment of which is strongly modified in the excited state corresponding to the first electronic absorption band. Calculations have been done with the DZVP/GEN-A2* combination of atomic and auxiliary basis sets and with the PBE XC functional. The LR-TDDFT absorption spectrum has been simulated by assigning a Lorentzian function of width 0.25 eV centered at each excited state energy and with an amplitude proportional to the oscillator strength of the transition. The first 493 (singlet) excited states have been included in the construction of the spectrum, spanning an energy window 18 eV wide. The LR-TDDFT calculation has been carried out with deMon2k within the framework of Auxiliary Density Perturbation Theory (ADPT)

92

. For the RT-

TDDFT spectrum three 15 fs length simulations have been carried out, each with a different orientation of the initial kicking electric field. The strength of the field was set to 0.005 a.u.. A time step of 1as was used in the propagation (using the PC Magnus propagator and the diagonalization technique). To

construct the spectrum the molecular dipole was damped by an exponential function ( ZB⁄ with

= 180 x. .) to broaden the absorption peaks in the RT-TDDFT spectrum. In both the LR- and RT-TDDFT

approaches fitted densities are employed to evaluate the Coulomb and the XC integrals. The spectra are depicted on Figure 2 The lowest energy transition with significant oscillator strength is found at 4.12 eV (300 nm). This excitation energy is close to the first band of the absorption spectra of coumarin in isopentane (310 nm)93. As illustrated on Figure 3, this electronic excitation mainly corresponds to ππ* transitions from the HOMO-2 and HOMO toward the LUMO with coefficients of 0.53 and 0.29 respectively. The agreement between the LR- and RT-TDDFT spectra is excellent over the entire range of energies. Both the positions of the maximum and the relative amplitude obtained match with the two types of TDDFT implementations. Our results validate the use of the Auxiliary DFT framework52 for constructing absorption spectra from RT-TDDFT simulations. As usual though, the choice of the auxiliary basis is critical and must include polarization functions (GEN-An*) to reach good accuracy.

Figure 2 Figure 394,95


Page 19 of 48

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


II.4 QM/MMpol framework in stationary and non-stationary cases (DFT+density fitting)/MMpol: stationary calculations. We now present the first results of QM/MMpol calculations with the in-deMon2k QM/MM module. To this end we consider a peptide (Tyr-Gly-Gly-PheMet) treated by DFT immersed in a box of 4,030 polarizable POL396 water molecules. The full system was previously equilibrated by classical MD simulations (data not shown). Five geometries extracted from the classical MD simulations have been calculated at the QM/MMpol level. For the latter we have used the PBE functional and the DZVP-GGA atomic basis set in combination with the GEN-A2, GEN-A2* or GENA3* auxiliary sets. The fitted density has been used to calculate both the Coulomb and XC potentials52. A grid of high accuracy is used to integrate the XC contributions (10-7 Ha). The induced dipoles have been updated at every SCF cycle by an iterative procedure until the Root-Mean-Square between two successive cycles is below 10-9 D. To avoid polarization catastrophe we have used the Tinker-exponential forms of field attenuation (Eqs 22-24). For the calculation of the electric field created by the electron density (Eq. 10) we consider two options. The first one is to use the KS density, the other is to use the auxiliary fitted density, respectively referred to as FBASIS and FAUXIS. Table 2: Induction energies and timings of stationary QM/MMpol calculations. QM/MMpol

QM/MMpol

(FBASIS)

(FAUXIS)

GEN-A2

GEN-A2*

GEN-A2

GEN-A2*

GEN-A3*

0.03

0.00

0.04

0.02

0.02

2.74

0.00

26.94

1.49

1.56

6.35

0.00

55.61

2.97

3.06

3.64

0.00

28.70

1.46

1.50

QM electric field

1060 (61%)

657 (43%)

3 (≈0%)

3 (≈0%)

5 (≈0%)

MM dipole iterations

236 (14%)

196 (13%)

246 (36%)

163 (18%)

154 (13%)

dipole embedding

167 (10%)

133 (8.8%)

169 (25%)

106 (12%)

121 (10%)

Total induction

1463 (85%)

986 (65%)

418 (61%)

272 (30%)

280 (23%)

20.9

29.6

8.4

16.1

20.0

a

Induction energies (kcal/mol) ;< ÎÏÐ"YZ[>> ) ;< ÎÏÐ"YZ?=> ;< ÎÏÐ"YZ@ ;< ÎÏÐ"B5B

Timings (s)

b

SCF (s/cycle)

a: We report the RMSDs over five conformers of the peptide taking the FBASIS/GEN-A2* as reference method. b: total time spent in SCF divided by the number of SCF cycles to reach convergence. For the timing, the numbers within brackets represent the percentage of the total time spent to compute the given contribution.



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

Page 20 of 48

Let us first consider the QM/MMpol-FBASIS scheme. The timings reported in Table 2 indicate that including induction significantly increases the cost of the calculation. It represents 65% of the total time spent in the SCF module with GEN-A2*. This is for the most part the calculation of the electric field created by the electron density on the MM sites that leads to this timing. 30 s are needed on average for each SCF cycle. For comparison, a non-polarizable DFT/MM calculation on the same system requires 12 s per SCF cycle. There is indeed a clear cost to the inclusion of induction. When switching to the less accurate GEN-A2 auxiliary set the computational cost is decreased to 21 s/SCF cycle, but more SCF cycles are needed to reach convergence so that 85% of the time is spent in the induction modules. In an attempt to reduce this supplementary cost, we consider the QM/MMpol-FAUXIS scheme. With GEN-A2* the cost drops to 16.1 s per SCF cycle, which is a just a little more than a non-polarizable DFT/MM (12 s). This reduction comes from the calculation of the electric field created by the electron density that has dropped to almost zero. The remaining time for induction is spent in the induced dipole iterations and in the calculation of the induced-dipole electrostatic potential. To assess whether the FAUXIS is reliable we report in Table 2 the different contributions to the polarization energy (Eq. 11). For each contribution, the Root-Mean-Square-Deviations are calculated taking the FBASIS/GEN-A2* as reference. At SCF ;< ;< + "YZ@ + "WoV convergence we find that the RMSD of induction energies ("YZ?=> L−FÏÏ ) differs by less than

1.5 kcal/mol between the FBASIS and FAUXIS approaches with GEN-A2*. A similar value is found for FAUXIS/GEN-A3*. We can thus consider the FAUXIS/GEN-A2* combination to be an excellent approximation of FBASIS/GEN-A2*, without further need to go to GEN-A3*. On the other hand, the comparison is less encouraging with GEN-A2. In this case, the RMSD of the induction energy between the

FAUXIS and FBASIS approaches is 29 kcal/mol. This is a rather large value: the fitted density can be used safely in lieu of the KS density in QM/MMpol calculations as long as sufficiently flexible auxiliary basis sets are used (GEN-A2* or larger). In such cases the computational advantage of the FAUXIS approach is significant. We make two final remarks. First, rather expectedly, the induction terms that are the most sensitive to the choice of the auxiliary basis are those between the MM induced dipoles and the QM

;< region. The induction energy between the dipoles and the MM charges ("YZ[>> ) is always the same.

Second, we note that the SCF process converges more rapidly with larger auxiliary basis sets, somewhat lowering the increase in computational cost due to the greater number of integrals to compute.

(DFT+density fitting)/MMpol: non-stationary calculations. We now examine the sensitivity of the RTTDDFT/MMpol simulation to the method chosen to calculate the electric field generated by the QM region. After SCF convergence the electronic density of the peptide is perturbed by a Gaussian shaped


Page 21 of 48

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


electric field centered at 10 as and with standard deviation 1 as. Three field strengths have been tested: 0.001, 0.01 and 0.1 a.u. The simulations have been run for 1 fs with the propagator-corrector Magnus scheme and a time-step of 1 as. The BCH formula (Eq. 39) has been used to calculate the exponential of the complex matrices. 30 terms have been used in the expansion of Eq. 39. We report in Figure 4 the evolution of the peptide dipole moment (Top), and the difference of polarization energy with respect to their values at the initial time. Both the FAUXIS and FBASIS schemes are tested. Each simulation has been carried out with either the GEN-A2 or GEN-A2* auxiliary basis set. As expected the stronger the intensity of the electric field perturbing the electron density at the beginning of the propagation, the larger the response of the system. This can be seen on the evolution of the peptide dipole moment, that exhibits larger amplitudes with the stronger field (0.1 a.u.), and also on the fluctuations of the polarization energy. With the weakest field the subsequent variations of the polarization energy are small (of the order of 10-2 kcal/mol). The agreement between the FAUXIS and FBASIS approaches is always very satisfactory. With the GEN-A2* auxiliary basis set the results between the FAUXIS and FBASIS simulations are even indistinguishable whatever the strength of the electric pulse that perturbs the system.

Figure 4

Coupling between electron and MM dipoles in RT-TDDFT/MMpol simulations. As explained in the methodology section a central aspect of the present RT-TDDFT/MMpol implementation is the assumption that the MM induced dipoles respond instantaneously to the electronic motion taking place in the QM region. Technically, this means that the MM dipoles are fully converged at every RT-TDDFT step. We tested the suitability of this strategy by repeating the previous simulations of the solvated peptide with shorter time steps of 0.75, 0.5, 0.25 and 0.1 as. For short enough time steps the decoupling approximation is certainly valid. With a time step of 0.1 as we found, indeed, that the MM induced dipoles evolve very smoothly. Figure 5 depicts the differences of QM/MMpol total energy, of polarization energy, and of embedding energy as a function of time taking the 0.1 as time-step simulation as reference. The initial perturbing field was set to 0.001 a.u. Clearly, the larger the time step the larger the difference. For the total energy and for the embedding energy the maximum error is of the order of a few thousandths of a kcal/mol with a 1 as time step. It is an order of magnitude smaller for the polarization energy. These values are rather small compared to the variations of the total energy in these simulations caused by the initial perturbation with the external electric field (around 0.06 kcal/mol). Interestingly the energy errors fluctuate around zero. This suggests that the simulations with time steps 21 ACS Paragon Plus Environment


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

Page 22 of 48

larger than 1as eventually depart from the reference trajectory but do not diverge from it. We also found that the peptide dipole as well as the water dipoles of the first solvation layer (i.e. those mainly impacted by the electron dynamics taking place on the peptide) were within 1.0E-5D from those of the reference trajectory. This is a very small value. In simulations in which the initial perturbing electric field strength was increased to 0.01 a.u. the same trends are obtained albeit with a factor of ten in the amplitudes of the errors (Figure S1). This again seems acceptable in view of the overall total energy change (6 kcal/mol). For an even stronger perturbing field (0.1 a.u.) the errors in total, polarization energy and embedding energies are of the order of a kcal/mol, a tenth of a kcal/mol and a hundredth of a kcal/mol, respectively (Figure S2). These values are quite high, but again much smaller than the fluctuations of the total energy of the molecule (around 50 kcal/mol). We note that such electric fields are extremely strong and would trigger non-linear effects like ionization. Altogether these tests justify the nonstationary/stationary coupling scheme between RT-TDDFT for the QM part and stationary MMpol for the environment although one should be careful to adapt the propagation time step to the amplitude of the electronic fluctuation that take place in the QM region. The most suitable time step might depend on the particular systems of interest. In principle though there should be a time step beyond which the decoupling between electrons and MM dipoles ceases to be valid. When we increased the time step (2 or 5 as) the electronic propagation was not stable anymore and diverged in a few steps. RT-TDDFT propagations are usually very sensitive to discontinuities that may arise in the time-dependent KS potential. Hence, a plausible explanation for the numerical instabilities observed in RT-TDDFT/MMpol simulations for the largest time steps may stem to potential discontinuities caused by significant variations of MM induced dipoles between two propagation steps. Interestingly, sudden instability of electron dynamics propagation may thus well be a sign of the breakdown of the decoupling hypothesis between the electron cloud dynamics and induced MM dipoles. Further work will be needed to examine this point in more detail.

Figure 5

II.5 Absorption spectra from RT-TDDFT/MMpol simulations We report on Figure 6 the absorption spectrum for the coumarin molecule solvated by a 30 Å radius sphere of water molecules. The system was equilibrated in a previous step by a classical MD simulation (data not shown) followed by a few-steps geometry optimization of the coumarin molecule before the 22 ACS Paragon Plus Environment

Page 23 of 48

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


RT-TDDFT/MMpol simulations. This partial optimization was intended to avoid too large distortions of the molecule, which would cause unreasonable displacements of the electronic absorption bands. We followed a similar protocol as for the gas phase case (see II.3) to build the spectra. In most simulations the QM region encompasses the coumarin molecule while the water environment is described by the force field. However in order to define a reference spectrum we have also carried out RT-TDDFT/MMpol simulations including the first hydration layer of coumarin in the QM region (12 water molecules). As seen on Figure 6 the environment has a significant effect on the spectrum. This is especially noticeable for the band around 4 eV that corresponds to a transition having charge transfer character. The center of this band is red-sifted by 0.22 eV when comparing the gas phase and the reference spectra (see inset). The non-polarizable TIP3P model leads to a red-shift of 0.15 eV, while the polarizable RTTDDFT/MMpol (POL3) spectrum exhibits a more pronounced red-shift of 0.18 eV, that is closer to the reference. Electrostatic induction thus permits a slight improvement on the position of the most displaced absorption bands. We also remark that the absorption spectra calculated with polarizable water combined with either the FBASIS or the FAUXIS option (see above) are indistinguishable. This is a further element showing that as long as sufficiently flexible auxiliary basis sets are chosen (GEN-A2* here), density fitting techniques can be safely used in TDDFT/MMpol calculations.

Figure 6

II.6. Numerical stability and performances We conclude section II with some notes on the numerical stabilities of RT-TDDFT simulations and on the computational performance of the implementation in deMon2k. It is known that simulating electron dynamics by RT-TDDFT can be difficult in terms of numerical stability. This numerical stability is highly system-dependent and the best propagation scheme has to be sought for each new molecular system of interest. The list of parameters impacting the numerical stability encompasses not only those directly related to the propagation schemes (propagator, matrix exponentiation method, size of time-step…) but also other more general DFT parameters: the quality of atomic basis set, of auxiliary basis set, the grid quality for XC potential numerical integration... In particular, we found the initial conditions of the electronic propagation to be extremely important. For instance, in most of our calculations the preliminary stationary SCF needs to be converged with tolerance criteria below 10-10 Ha for the total electronic energy and 10-7 for the charge density error (respectively defined by the TOL and CDF options 23 ACS Paragon Plus Environment


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 24 of 48

of the SCFTYPE keyword in deMon2k). These convergence thresholds are tight compared to those customarily used in stationary DFT calculations. We also found the method for determining the fitted density coefficients to be important. To solve the sets of inhomogeneous systems of linear equations associated with density fitting, one can either use an analytical97 or numerical approach98. The former generally led to more stable RT-TDDFT simulations and has been used throughout. In Figure 7 we report the computational timing to carry out a 1.5 fs RT-TDDFT/MMpol electron dynamics simulation on the solvated peptide just described. Following our previous conclusions, this simulation has been run with the AUXIS and FAUXIS approaches with the DZVP-GGA/GEN-A2* combination of basis sets. We have used the predictor-corrector Magnus propagator with the BCH expansion (30 terms) and a time step of 1as. A grid of accuracy 10-7 Ha has been used for the XC contributions. The time-dependent auxiliary density was integrated at every RT-TDDFT step to extract intrinsic atom multipoles (charges,

dipole, quadrupole) according to the Hirshfeld scheme. With GEN-A2* multipoles extracted from are very close to those extracted from

69

.The simulation took 12 h on 48 processors with the message

passing interface protocol.

Figure 7

The most time-consuming part of the simulation corresponds to the matrix multiplications (of which there are almost 400 000). These are needed i) to transform the density and KS matrices between the MO and AO representation, ii) in the BCH approximation that involves nested commutators. In second position is the cost of including induction stemming from the embedding of the QM region by induced dipoles (8%) and by the convergence of MM induced dipoles by the iterative procedure (18%). This situation could be improved in the future by adopting more advanced simulation algorithms that avoid the iterative procedure used here to converge the MM induced dipoles99. The XC potential represents 24% of the overall computational cost while the calculation of the Coulomb contribution is almost negligible. This remarkable performance is possible thanks to the mixed scheme implemented in deMon2k to store short range ERIS in RAM and to compute long range ERIS by double asymptotic expansions100. Finally we note that the repeated analysis of the time-dependent electron density induced only a small supplementary cost to the calculation69. In summary, we think our RT-TDDFT/MMpol is efficient enough to tackle the simulation of electron dynamics in large molecular systems.


Page 25 of 48

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


III. ELECTRONIC POLARIZATION DYNAMICS III.1 Electron dynamics within an isolated molecule. In this subsection we analyze the electron dynamics taking place when a coumarin molecule in the gas phase is submitted to a monochromatic laser field. We choose a wavelength that corresponds to an excitation energy of 4.1547 eV (see section II.3). We have propagated the electron density for 25 fs with a time step of 10 as. A Gaussian shaped electric field centered at 9.5 fs and full-width at half maximum of 3.53 fs ( = 1.5uw) was applied (Eq. 40). The electric field orientation was set in the molecular plane

with an intensity of 10-4 a.u. (51 mV/nm). The frequency ± of the field was set to correspond to an excitation energy of Δ" = 0.1528 #x. To follow the evolution of the electron density we have

charges F and dipole moments (8Ñ ). To define the latter three distinct population schemes have been computed the molecular time-dependent dipole moment and the time-dependent intrinsic atomic

tested, namely the Hirshfeld101, Becke102 and Voronoi Deformation Density103 schemes. Only the

Hirshfeld results are shown in the main text while the others are given in SI. Because repetitive population analyses along the RT-TDDFT simulations may become time-consuming, we have extracted

atomic charges and dipoles from the analysis of the auxiliary density () which is cheaper than from the

KS density69. We have previously shown this alternative to be reliable as soon as the GEN-A2* auxiliary basis (or larger) is chosen, as is the case in this application.

The set of time-dependent atomic charges (FÑ ) permits the definition of a total charge-derived dipole

moment 8[ = ∑Ñ FÑ Ñ . The variations in time of 8[ provide information on charge transfer between

atoms in the course of the simulation. On the other hand, the sum of the intrinsic atomic dipoles 8_5* = ∑Ñ 8Ñ provides information on the internal polarization of the atoms during the simulation. The

sum of 8[ and 8_5* gives rise to the full molecular dipole moment. As seen in Figure 8 the molecular

dipole starts to oscillate at around 7fs as a result of the application of the electric field. Regular

oscillations are seen in subsequent times, the period of which (0.99 fs) corresponds well to the applied electric field energy (4.1547 eV). The molecular system thus undergoes Rabi oscillations between the ground and the targeted excited state. We observe soft beatings on a few fs that correspond to the Gaussian shaped pulse. Looking at the decomposition of the dipole moment we find that the oscillations

are essentially due to 8[ , i.e. they are caused by charge transfers between atoms. The sum of the

internal polarization 8Ñ (8_5* ) fluctuates much less, as well as each intrinsic dipole moment (Figure 8, right). Of course, the separation between 8[ and 8_5* is arbitrary and depends on the chosen population 25 ACS Paragon Plus Environment


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 26 of 48

scheme. With Voronoi Deformation Density, we found very similar results as for Hirshfeld (Fig S3). On the other hand with the Becke partitioning scheme a very different separation between polarization and charge transfer is obtained (Fig S4). Actually the Becke scheme is not recommended for extracting electrostatic multipoles because it may produce “non-chemical” charges in many cases (e.g. hydrogen always have charges around -0.5e). It should probably be avoided for analyzing time-dependent electron density. That said we found fluctuations of the total dipole moment come from fluctuations of the charge-derived dipole moments whatever the chosen partitioning scheme.

Figure 8

The isosurfaces of the deformation electron density shown in Figure 9 render a pictorial view of ultrafast dynamics over three Rabi oscillations. For this simple excitation process involving two electronic states the overall shape of the deformation density can be rationalized looking at the shape of the MOs involved in the process. Especially the population of the LUMO can be identified looking at Figure 9.

Figure 9

III.2 Dynamics of the response of the environment In this section, we analyze the electronic response of the environment of a central molecule after perturbation of the latter by an external electric field. To this end we take the same system as in II.4, namely a methionine enkephalin solvated in a box of POL3 water molecules. After tight SCF convergence, the central peptide was perturbed by a Gaussian shaped electric pulse centered at 20 as with 3 as width. The field strength was set to either 0.001, 0.01 or 0.1 a. u. Note the latter corresponds to a very strong intensity. Following our conclusions from the validation section we use the AUXIS and FAUXIS approaches in combination with the GEN-A2* auxiliary basis. The simulation was conducted for 3 fs with a time step of 3 as using the predictor-corrector-Magnus/BCH propagator. We report on Figure 10 the

variation of the induced dipoles on MM atoms with respect to the initial time (∆L = L − L0) and

their normalized auto-correlation functions (Ò , ACF). Both quantities are averaged by hydration layers

as indicated by the angular brackets 〈… 〉.

We start by considering the upper graphs that correspond to perturbing field strength of 0.001 a.u.. As expected the longer the distance between the water molecules and the peptide, the smaller the impact on the induced dipoles. The first hydration layer is the one that experiences the highest variations of 26 ACS Paragon Plus Environment

Page 27 of 48

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


induced dipoles. As evident from the black curve in Figure 10, Top-Left, the average induced dipoles undergo damped oscillations. These are caused by energy dissipation in the MM environment, which is possible thanks to the use of a polarizable FF. Dissipation is very pronounced for the first hydration layer but it is also seen for the outer hydration layers. The induced dipoles for molecules pertaining to the inner hydration layer completely lose correlation within a few tens of as, while beyond 15 Å, the average ACF remains close to 0.8 at 200 as. The characteristic response time is distance dependent. This characteristic time increases for each successive hydration shells. Some of the averaged ACF become

negative which is to be related to the oscillatory nature of the variations of 〈∆L〉. We finally remark that the response of MM induced dipoles is not fully instantaneous but also exhibits relaxation components over hundreds of attoseconds. When the strength of the initial perturbing field is increased to 0.01 a.u. the amplitude of oscillation of the average induced dipoles is larger by a factor of ten. This is true for each hydration shell. When the field strength is further increased to 0.1 a.u. a further increase of response amplitudes is observed for 〈∆L〉. The ACFs exhibit more complex evolutions with the increase of perturbing field strength. For the

weakest perturbing field strength (0.001a.u.) we already mentioned that the response was distance dependent (Top-Right). For a perturbing field strength of 0.01 a.u. the response of induced dipoles is not distance dependent within the first 50as, apart for water molecules situated beyond 15 Å (Middle-Right). Only after this time a scattering of the average ACFs becomes apparent. Finally, for a perturbing field strength of 0.1 a.u. all the average ACFs but one (again for water molecules situated beyond 15 Å) are almost superimposed (Bottom-Right). The response mechanism of MM induced dipoles is therefore not distance dependent at all within 15 Å. All these results reflect subtle response mechanisms that deserve a detailed analysis. We recall that the induced dipoles are determined by the electric field created by the other MM atoms :

(9

=> ;< + 9;< ) and by the QM region (9 ) (Eqs. 6-10). In the present RT-TDDFT simulations only 9

and 9=> can account for the variations of the MM induced dipoles since the nuclei are fixed. We also

recall that we employ here a mixed non-stationary/stationary RT-TDDFT/MMpol scheme (see section II). Accordingly the response of the MM induced dipoles caused by variations of 9=> are expected to be

enhanced in our scheme compared to what they would be in fully dynamical simulation. Nonetheless the

average 〈∆L 〉 and the associated ACF extracted from a 0.1as time-step simulation were found to be

very similar to the graphs shown on Figure 10, thereby indicating that the artificial enhancement of the dipole relaxation due to the RT-TDDFT/MMpol coupling scheme is moderate. The response of MM 27 ACS Paragon Plus Environment


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 28 of 48

dipoles should be less pronounced as the distance 3 increases because of the decay of 9

=>

with

should be associated with a distance. On the other hand, the response caused by variation of 9;<

certain delay, because it requires the other induced dipoles to be affected. For example, the induced dipoles of water molecules situated between 6 and 9 Å will be affected by induced dipoles of innermost hydration waters only when their induced dipoles have varied. The prevalence of one mechanism over the other should depend on the relative strength of 9;< and 9=> .

For the inner hydration layer (< 3 Å) 〈Ò 〉 is almost identical whatever the initial perturbation field. It

decays to 0.3 in around 50as (although the variations of induced dipole amplitudes are different for each

perturbing electric field). For this hydration layer the source of variation of MM induced dipoles is =>

primarily9 . It is the time-dependent field created by the electron cloud of the peptide that determines the response of MM induced dipoles. The oscillations of the MM induced dipoles essentially

follow that of the peptide dipole moment (Fig. S5). For the outer hydration layers, the response

and 9=> mechanism depends on the relative importance of 9;< , the latter being itself dependent on

the perturbing field strength. For the strongest perturbing field (0.1 a.u.) the response mechanism of the MM induced dipoles is completely imposed by 9=> for all hydration shells (except for water molecules

beyond 15 Å). This explains why the average ACFs are almost superimposed. The amplitude of the response decays with distance but the speed at which the induced dipoles vary is the same. In this

regime 9=> ≫ 9;< so that 9=> imposes the response mechanism: the MM induced dipoles within 15 Å

follow the variations of the peptide dipole moment. For the intermediate perturbing field (0.01 a. u.)

9=> dominates the response mechanism for the shorter distances (

Simulating Electron Dynamics in Polarizable ... - ACS Publications

Recommend Documents