Modeling and simulating the excited-state dynamics of a system with condensed phases: A machine learning approach
MModeling and simulating excited-statedynamics of a system in condensed phases:Machine Learning approach
Seiji Ueno ∗ , † and Yoshitaka Tanimura ∗ , ‡ † HPC Systems Inc., Japan ‡ Department of Chemistry, Kyoto University, Kyoto, Japan
E-mail: [email protected]; [email protected]
Abstract
Simulating quantum irreversible dynamics of exciton and electron transfer prob-lems possesses a non-trivial challenge. Because irreversibility of the system dynamicsresults from quantum thermal activation and dissipation caused by a surrounding en-vironment, it is necessary to include infinite environmental degrees of freedom in thesimulation. Because the capability of full quantum dynamics simulation that includesthe surrounding molecular degrees of freedom is limited, a practical approach is to em-ploy a system-bath model, in which dynamics of excitons or electrons are described by asystem Hamiltonian, while the other degrees of freedom, arising from the environmentalmolecules, are described by a harmonic oscillator bath (HOB) plus a system-bath in-teraction. By extending a previous study of a machine-learning approach for molecularliquids [J. Chem. Theory Comput. 2020, 16, 2099], here we construct a system-bathmodel for exciton and electron transfer problems. We determine both the system andsystem-bath interaction parameters that include the bath spectral distribution, using a r X i v : . [ phy s i c s . c h e m - ph ] F e b he electronic excitation energies obtained from Quantum Mechanics/Molecular Me-chanics (QM/MM) simulations as a function of time. Using analytical expressions ofoptical response functions, we calculate linear and two-dimensional electronic (2DES)spectra for indocarbocyanine dimer in methanol. From these results, we demonstratethe capability of our approach to elucidate the non-equilibrium exciton dynamics of aquantum system in a non-intuitive manner. Quantum dynamics play a significant role in many chemical physics and biochemicalphysics problems. Frequently studied problems of this kind are exciton and electron transferprocesses those are involved in photosynthetic, electron transfer,
DNA, andphotovoltaic systems.
In these problems, the environments (baths), for example, proteinand solvation play a central role: These baths are complex and strongly coupled to a molec-ular system of interest at finite temperature. Recent theoretical works have demonstratedthat the system and baths are quantum mechanically entangled (bathentanglement) thatare essential to properly understand the quantum dynamics displayed by the system.
For example, it has been shown that the optimal condition for excitation energy transfer inlight-harvesting complexes is realized under the nonMarkovian system-bath interactions ina strong coupling regime, in which the noise correlation time of the bath is comparable tothe time scale of the system dynamics. To perform high accuracy simulations with reducedcomputational cost, some approaches have utilized machine learning methods to make mod-els to reproduce open quantum dynamics, or estimate chemical properties for classicalmolecular dynamics.
Although irreversibility of the system dynamics results from quan-tum thermal activation and dissipation caused by a surrounding environment, it is difficultto conduct the quantum molecular dynamics simulation that exhibits such a characteristicfeature arising from the macroscopic degrees of freedom.Thus we introduce a system-bath model in which dynamics of excitons or electrons are2escribed by a system Hamiltonian, while the other degrees of freedom, arising from theenvironmental molecules, are described by a harmonic oscillator bath (HOB). The HOB,whose distribution takes a Gaussian form, exhibits wide applicability to simulate bath effects,despite its simplicity. This is because the influence of the environment can in many casesbe approximated by a Gaussian process, due to the cumulative effect of the large numberof environmental interactions. In such a situation, the ordinary central limit theorem isapplicable, and hence the Gaussian distribution function is appropriate.
The distinctivefeatures of the HOB model are determined by the spectral distribution function (SDF) ofthe coupling strength between the system and the bath oscillators for various values offrequency. By choosing the form of the SDF, the properties of the bath can be adjustedto represent vareous environments consisting of, for example, solid state materials andprotein molecules.
Because the SDF can be different for different form of a system-bathcoupling, it is not easy to find an optimized Hamiltonian associated with optimized SDF, inparticular for a bath describing a fluctuation of site-site interaction energy.In the previous study, we employed a machine-learning approach to construct a system-bath model for intermolecular and intramolecular modes of molecular liquid using atomictrajectories obtained from molecular dynamics (MD) simulations. In this study, we ex-tend the previous approach in order to investigate an exciton or electron transfer problemthat is characterized by electronic states embedded in molecular environment using Quan-tum Mechanics/Molecular Mechanics (QM/MM) calculations for the atomic coordinates ofmolecules. In particular, we are aiming at the exciton transfer process of the photosynthesisantenna system to investigate how natural systems can realize such highly efficient yields,presumably by manipulating quantum mechanical processes. As a demonstration, we con-struct a model Hamiltonian of an indocarbocyanine dimer compound. The accuracy of amodel is examined by calculating linear and two-dimensional electronic spectra.This paper is organized as follows. In Sec. 2, we introduce an exiciton or electron transfermodel that is coupled to a harmonic heat-bath. We then describe the machine learning3pproach that we use to determine the system parameters, the system-bath interaction,and the SDFs on the basis of QM/MM simulations. In Sec. 3, we present results for anindocarbocyanine dimer model constructed from the analysis of the QM/MD trajectories.Linear absorption and two-dimensional spectra are calculated from the analytical expressionsof linear and nonlinear response functions. Section 6 is devoted to concluding remarks.
We consider a situation in which an exciton or electron transfer system interacts withmolecular environments that give rise to dissipation and fluctuation in the system. TheHamiltonian of the system is expressed as ˆ H S = (cid:88) j (cid:126) ω j | j (cid:105)(cid:104) j | + (cid:88) j (cid:54) = k (cid:126) ∆ jk | j (cid:105)(cid:104) k | , (1)where the j th exiciton or electron states with energy (cid:126) ω j are represented by the bra and ketvectors as | j (cid:105) and (cid:104) j | . The interaction energy between the j th and k th states is describedby (cid:126) ∆ jk . In our model, each state is coupled to a different molecular environment (labeledby a ) that are treated as N a harmonic oscillators. The total Hamiltonian is then given by H tot = H S − (cid:88) a N a (cid:88) j =1 α aj ˆ V a ˆ x aj + (cid:88) a N a (cid:88) j =1 (cid:20) (ˆ p aj ) m aj + 12 m aj ( ω aj ) (ˆ x aj ) (cid:21) , (2)where the momentum, position, mass, and frequency of the j th oscillator in the a th bath aregiven by ˆ p aj , ˆ x aj , m aj , and ω aj , respectively. The system part of the system-bath interacting is4xpressed as ˆ V a = (cid:88) j,k V ajk | j (cid:105)(cid:104) k | , (3)where V ajk is the coupling constant for the a th bath between the j and k states. The a thheat bath can be characterized by the spectral distribution function (SDF), defined by J a ( ω ) ≡ N a (cid:88) j =1 (cid:126) ( α aj ) m aj ω aj δ ( ω − ω aj ) , (4)and the inverse temperature β ≡ /k B T , where k B is the Boltzmann constant. Various en-vironments, for example, those consisting of nanostructured materials, solvents, and proteinmolecules, can be modeled by adjusting the form of the SDF. For the heat bath to be an un-limited heat source possessing an infinite heat capacity, the number of heat bath oscillators N a is effectively made infinitely large by replacing J a ( ω ) with a continuous distribution. Theabove model has been frequently used in the analysis of photosynthetic systems, electrontransfer, DNA, and solar battery systems.
We consider pigments in a molecular system, whose electric excitation or exciton statesare described by Eq.(1). The electric states of the pigments depend on a configuration ofthe surrounding atoms at time t . The time evolution of the excited states of the system andenvironmental molecules are described by QM/MM simulations. Because our aim of con-structing the system-bath model is to perform full quantum simulation of the entire system,we should use quantum MD simulation to provide date on the basis of all atomic coordi-nates. In practice, however, it is impossible to treat large environmental degrees of freedom,quantum mechanically accurately. Fortunately, we are expected to have reasonable SDFsfor quantum simulation even we evaluated them using the classical MD simulation, becausewe utilize the ensemble of the molecular trajectories, which exhibit Gaussian distributionwhere the difference between the quantum and classical trajectories is expected to be minor,5nd because the dynamics of harmonic oscillators are identical in both classical and har-monic cases, as both classical and quantal Liouvillian for the j th oscillator is expressed as ˆ L j = − ( p j /m j )( ∂/∂x j ) − ( mω )( ∂/∂p j ) for the j th oscillator. We thus use the classical MDsimulation to acquire the atomic coordinates of the pigments and the molecular environment.We then conduct quantum chemistry calculations to obtain the desired electronic states, mosttypically the HOMO and LUMO states of the pigments as a function of time. The excitedenergy of the j th pigment is denoted by (cid:15) jj ( t ) and the interaction energy between the j thand k th pigment that includes the bath induced fluctuation is denoted by (cid:15) jk ( t ) . They areobtained using any kind of numerical programs for quantum chemistry calculations. In acase if the main system is too large to evaluate the entire electronic states, we evaluate thesite energy (cid:15) jj ( t ) and the interaction energy (cid:15) jk ( t ) separately. From the calculated (cid:15) jj ( t ) and (cid:15) jk ( t ) , we evaluate the system-bath coupling strength in ˆ V ajk and its SDF, in addition to theexcitation energy (cid:126) ω j and the interaction energy ∆ jk , on the basis of the machine learningapproach.While the SDFs evaluated from the MD simulation are temperature dependent, those forthe HOB are temperature independent: We eliminate the temperature dependence of opti-mized parameters, assuming that the sampled MD trajectories exhibit canonical ensembleat finite temperature. For n exciton or electronic excitation sites, we express the simulated data in terms of (cid:15) jk ( t ) as describing the excited and site-site interaction energies of interest obtained fromthe QM/MM simulation. The learning Hamiltonian is then expressed as H ( t ) = n (cid:88) j,k =1 (cid:15) jk ( t ) | j (cid:105)(cid:104) k | . (5)6e then attempt to reproduce the trajectories of (cid:15) jk ( t ) for the total Hamiltonian Eq.(2) withEqs. (1) and (3). Although the system-bath model considers an infinitely many degrees offreedom, here we employ a finite number of bath oscillators to estimate the SDFs. Then,sampling for the machine leering training is considered as the average of the classical bathoscillators for a choice of the system and system-bath parameters. The site energy andinteraction energy can be expressed as (cid:15) jj ( t ) = (cid:126) ω j − δ jj ( t ) (6)and (cid:15) jk ( t ) = (cid:126) ∆ jk − δ jk ( t ) , (7)respectively, where δ jk ( t ) is expressed in terms of the linear function of the bath coordinatesas δ jk ( t ) = (cid:88) a α ajk x ajk ( t ) . (8)Here the a th bath coordinate for the jk site is described as a function of time as x ajk ( t ) = A ajk sin (cid:0) φ ajk + ω ajk t (cid:1) , (9)where A ajk and φ ajk are the amplitude and phase of the jk th bath oscillator. The phase φ ajk is randomly chosen to avoid recursive motion of the oscillator. Here, we assume thatthe influences of individual bath modes are all independent and that the correlations of thefluctuations among different modes can be ignored, although we can treat such correlatedmodes separately by introducing additional baths.From Eqs. (8) and (9), δ jk ( t ) can be expressed as δ jk ( t ) = (cid:88) a c ajk sin (cid:0) φ ajk + ω ajk t (cid:1) , (10)7here c ajk = α ajk A ajk (11)and we treat the system-bath coupling parameters as the product of these two variables. Thebath parameters and the system-bath interactions are expressed as a set of latent variablesin the machine learning context, defined as θ = ( { ω j } , { ∆ jk } , (cid:8) c ajk (cid:9) ) , (12)where { ... } is the set of system and bath parameters. The trajectories of (cid:15) jj ( t ) and (cid:15) jk ( t ) ,obtained from the QM/MM calculations, are described as the vibrational motions of thepigment molecular and the surrounding molecules. We then assume that the probabilitydistribution of the pure state energy λ i is determined from a Gaussian process and aredescribed by the set of the bath parameters α ajk and φ ajk , by optimizing the probabilitydistribution given by P ( λ i | θ ) = (cid:90) (cid:89) k,j,a dφ ajk P ( λ i | θ ; φ ajk ) P ( φ ajk ) , (13)which is marginalization of the phase of the oscillators φ ajk , where P ( φ ajk ) is the uniformdistribution of [0 , π ) . The marginalization is introduced in order to avoid a trapping in alocal minimal sate during the gradient method. The probability distribution is expressed inthe normal form as P ( λ i | θ ; φ ajk ) ∝ exp (cid:2) − σ ( λ i − E i ) (cid:3) , (14)where E i ≡ E i ( θ ; φ ajk ) is the predicted energy as a function of the parameter set θ and initialphase φ ajk for the model Hamiltonian Eq. (5) and σ is the error width. Our purpose of em-ploying a machine learning method is to choose the optimal parameter set in E i ( θ ; φ ajk ) thatmakes the probability distribution for given data λ i maximize. Among several optimizationmethods, here we use the maximum likelihood method (MLE), where the loss function is8xpressed in terms of the negative log of the probability as L = (cid:88) i ( λ i − E i ) . (15)To find the maximum value of L , we employ the Adam gradient method for optimization ofthe parameter set as θ ← θ + γ ∂L∂θ , (16)where γ is the learning rate. In this way, we obtain the J jk element of the SDF for the state jk . Because the energy distribution of each bath oscillator E ajk = m ajk (cid:0) ω ajk (cid:1) is assumed toobey an canonical ensemble, the oscillator amplitude can be expressed as (cid:10) A ajk (cid:11) = 1 (cid:113) πβm ajk (cid:0) ω ajk (cid:1) . (17)Using Eqs. (11) and (17) for Eq.(4), we obtain J jk ( ω ) = N a (cid:88) a =1 πβ (cid:126) ω ajk ( c ajk ) δ ( ω − ω ajk ) . (18)Because J jk ( ω ) rapidly changes in time in accordance with the structure change of thepigment molecules and environment, we evaluate J jk ( ω ) by averaging over the differentsample trajectories. In the mathematical aspect, c jk is the frequency domain expression ofthe time domain data and J jk ( ω ) can be obtained by averaging the power spectra c jk fromthe Wiener-Khinchin theorem.It should be noted that the absolute intensity of J jk ( ω ) cannot be determined in theframework of this study, because we do not evaluate the dipole moment of this complex ma-terial for simplicity: We evaluate the intensity of J jk ( ω ) from the width of the experimentallyobtained linear absorption spectrum. 9 Numerical demonstration
Figure 1: The molecule structure of indocarbocyanine dimer. Two pigments are connectedby methylene chains. The gray/blue/white atoms represent the carbon/nitrogen/hydrogen,respectively. The red square represents the pigment 1, whereas the blue square representsthe pigment 2.
We now demonstrate our approach for a dimer of identical indocarbocyanine molecules. Figure 1 displays the structure of the pigment molecule. The ground and the excited statesof each pigment are expressed as | (cid:105) j and | (cid:105) j for j = 1 and 2, respectively. The groundstate energies are chosen to be zero. The system Hamiltonian is then expressed as ˆ H = ω ( | (cid:105) (cid:104) | + | (cid:105) (cid:104) | )+ ∆ ( | (cid:105) (cid:104) | + | (cid:105) (cid:104) | ) , (19)where ω is the excitation energy of a pigment and ∆ is the interaction energy between thedimer. By diagonalizing H , we obtained the eigenvalues ω k for the k = + and − eigenstates10 (cid:105) = ( | (cid:105) | (cid:105) + | (cid:105) | (cid:105) ) / √ and | −(cid:105) = ( | (cid:105) | (cid:105) − | (cid:105) | (cid:105) ) / √ , respectively, as ω ± = ω ± ∆ . (20)The fluctuation of excitation energies and interaction energy as functions of time that arisefrom the intramolecular motions of the pigment and the intermolecular motions of surround-ing molecules are expressed as δω ± ( t ) and δ ∆( t ) , respectively. They are evaluated from thequantum chemistry calculations for given atomic trajectories of the entire molecular systemevaluated from a MD simulation.In our model, because each exciton state is localized and the effects of the environmentalmodes are exciton site specific, we employ individual heat-bath for a description of energyfluctuation in each exciton state. Although it is possible to introduce a global heat-bath forlow-frequency environmental modes coupled to multiple exciton states, here we ignore sucheffects. Then the fluctuation of the excitation energy and interaction energies are expressedas δω ± ( t ) = , (cid:88) m w ± ,m ( t ) (cid:88) a c aω m sin( ω a t + φ aω m ) , (21)and δ ∆( t ) = w ∆ ( t ) (cid:88) a c a ∆ sin( ω a t + φ a ∆ ) , (22)where c abm is the amplitude (scaled by α as described in Eq. (11)) and φ abm is the initialphase of the a th oscillator for the state index b = 11 (or 22) and . We introduced thelocalization weight functions w ± ,m ( t ) and w ∆ ( t ) to describe the pigment specific effects ofthe environments in the delocalized exciton states representation that is obtained from thediagonalization of the pigment based Hamiltonian expressed as Eq. 19. They are evaluatedfrom the electronic states of the pigment m = 1 and on the basis of the atomic orbitals11AO) obtained from quantum chemistry calculations.Thus the targeting eigenenergies to be described by the system-bath model, λ ± ( t ; θ ) , areexpressed as λ ± ( t ; θ ) = ω + δω ± ( t ) ± (∆ + δ ∆( t )) , (23)where θ is the set of parameters θ = (cid:0) ω , (cid:8) c a ± ,m (cid:9) , ∆ , { c a ∆ } (cid:1) . As learning data, we computethe exciton energy E ± ( t ) , the molecular orbital (MO) coefficients for each exciton state, andwavefunctions (atomic orbital (AO) coefficients for each MO) from the quantum chemistrycalculations for given atomic coordinates as a function of time, while the movements ofall atoms in the system are evaluated from the classical MD simulation. Using these, weoptimize the set of parameter θ . In order to evaluate the weight function w k,m ( t ) , we calculatethe exciton and hole populations p exk,m ( t ) and p hk,m ( t ) that are obtained as the summation ofthe absolute square of the AO coefficients, which are evaluated from the AO coefficientsinvolved in the MO in the pigment m for the excited state k . The weight function is thenevaluated as w k,m ( t ) = p exk,m ( t ) p hk,m ( t ) and w ∆ ( t ) = (cid:80) k = ± (cid:0) p hk, ( t ) p exk, ( t ) + p hk, ( t ) p exk, ( t ) (cid:1) . Asthe definitions indicated the exciton is localized when w ± ,m is close to 1, whereas the excitonstates are mixed among the pigments when w ∆ is close to 1.To optimize the system and bath parameter set, we minimize the loss function L = (cid:88) n (cid:88) t L n ( t )= (cid:88) n (cid:88) t (cid:2) ( λ − ( t ; θ n ) − E n − ( t )) + ( λ + ( t ; θ n ) − E n + ( t )) (cid:3) , (24)where E n − ( t ) and E n + ( t ) are the lowest ( | −(cid:105) ) and 2nd lowest ( | (cid:105) ) excitation energies, andthe index n indicates the n th sample. Using the MLE method, we optimize c aω m and c a ∆ for each time series as a sample set. In order to apply the machine learning algorithm,time series of the tuple ( E n − ( t ) , E n + ( t ) , w nk,m ( t )) are regarded as the input feature variables.12n the indocarbocyanine case, the two pigments are symmetric, and the SDF of the bathfor each pigment is considered to be identical. Therefore we use the averaged value c aω =( c aω + c aω ) / . We then evaluated J jj ( ω )( j = 1 , and J ( ω ) , which are J ( ω ) for ω and ∆ , from c aω and c a ∆ using Eq. (18), respectively. Commonly used approach to evaluate the SDFs for (cid:15) ij ( t ) utilizing the Fourier trans-formation of the auto-correlation function expressed as F [ (cid:104) δ(cid:15) ij (0) δ(cid:15) ij ( t ) (cid:105) ] , where δ(cid:15) ij ( t ) ≡ (cid:15) ij ( t ) − (cid:104) (cid:15) ij (cid:105) . In the actual calculation, the time series (cid:15) nij ( t ) , where n is the index of a sample,are evaluated as the average of the auto-correlation function expressed as C ij ( t ) = 1 N (cid:88) n (cid:10) δ(cid:15) nij (0) δ(cid:15) nij ( t ) (cid:11) , (25)where N is the number of the total sample. We then obtain the SDF as J ij ( ω ) = F [ C ij ( t )] . (26)Alternatively, using the Wiener-Khinchin’s theorem for stationary random processes, we canobtain the SDF as an average of power spectrum P nij ( ω ) = F [ (cid:15) ij ( t )] as J ij ( ω ) = 1 N (cid:88) n P nij ( ω ) . (27)Although this Fourier based approach is simple and straightforward, the obtained SDFs arenot necessary the optimal choice for the system-bath Hamiltonian to describe the QM/MMdata, because the exciton and interaction energies are mutually depend on each other andthus J ij ( ω ) and J ik ( ω ) cannot be evaluated separately. In the machine learning approach,however, it is possible to optimize not only J ij ( ω ) and J ik ( ω ) , but also ω and ∆ withoutassuming the explicit relationships between the SDFs and system parameters. Moreover, if13ecessary, we can introduce additional conditions for optimization of the SDFs and systemparameters, as we employed w k,m ( t ) to take into account the effects of the exciton localizationfor the indocarbocyanine dimer. Step 1: Classical MD
We prepared a system consists of an indocarbocyanine dimer molecule with 1024 methanolmolecules for solvent. The classical MD simulations were carried out with the GROMACSsoftware package.
The condition for preparation MD was set to be 1atm at 300K withthe NPT ensemble. The equilibrium MD was carried out for 20 ps in the NVT ensemblefollowed by sampling MD with 5 ps in the NVE ensemble. These equilibrium MD and sam-pling MD were repeated 100 times. The entire MD simulation was done with the time stepof 0.1 fs.
Step 2: Quantum Chemistry Calculation
To obtain the sample trajectories of the excitation energies, we conducted ZINDO calcu-lation and natural transition orbital analysis for 1 femto-second period of one sampleusing the ORCA software package. We then obtained 100 samples of ( E − ( t ) , E + ( t ) , w ( t )) in 5 ps length. Step 3: Machine Learning
We extracted 1000 points of data in the interval 2 fs. They were taken from the 5 pslength of data obtained from step 2 with shifting the starting time in each 1 step. Weexcluded a sample if E + − E − is smaller than × − eV. We then extracted the 721 seriesof ample in the 2 ps duration from the 100 samples in the 5 ps length of data. These samplingdate was used as the input feature values in the machine learning calculations.14o perform learning calculations, we developed Python codes using the Tensorflow li-brary. The training was performed until the value of the loss functions converges withthe learning rate γ = 1 × − . The initial values of the targeting optimize variables forthe amplitudes of SDFs were set as c amω = 1 × − and c a ∆ = 1 × − , and the excitonand interaction energies were set as ω = ( (cid:104) E + (cid:105) + (cid:104) E − (cid:105) ) / and ∆ = ( (cid:104) E + (cid:105) − (cid:104) E − (cid:105) ) / . Theinitial phases φ ab were randomized 10 times for each series of sample. The loss functions wereaveraged for each 10 samples as a minibatch, while the parameters were optimized everyminibatch. For 721 samples, each 1 epoch contained 73 iterations. To minimize the lossfunction, we employed the Adam algorithm. Step 4: Spectrum Calculation
We assumed the dipole operator for the indocarbocyanine dimer was given by ˆ µ + ˆ µ = µ ( | (cid:105) (cid:104) | + | (cid:105) (cid:104) | + | (cid:105) (cid:104) | + | (cid:105) (cid:104) | ) creating the transition between the triplet state, | (cid:105) , | (cid:105) , and | (cid:105) , while the optical transitions from these states to the singlet state | −(cid:105) wereforbidden. Thus the optical transitions in the present system were described by the three-level system with the eigenenergies 0, Ω + , and ω . This allowed us to employ the analyticalexpressions of the linear and nonlinear response functions, as presented in appendix A. Wethen calculated the linear absorption and 2D electronic spectroscopy signals using the line-shape functions. Representative examples of the prepared dataset are plotted in Fig. 2. The abrupt changeof the exciton energies in Fig. 2(b) occurs due to the exciton transfer between the pigment1 and 2 that takes place in the time duration of 10-100 fs. As illustrated in Fig. 2(c), thedifference of the exciton energies E + − E − exhibits minima in accordance with the excitontransfer processes. As illustrated by the red-circles in Fig. 2(c), although such minimal pointare significantly sharper than the minimal point arise from the fluctuation of energies, it is15igure 2: Samples of data used for a learning calculation for (a) the excitation energy E k for k = ± , (b) the weight functions w k,m ( t ) for the pigment m = 1 (blue) and m = 2 (orange)for k = + (solid line) and k = − (dotted line), respectively, and w ∆ (green dashed curve),and (c) the difference of the energy levels E ± . The panel (c) is plotted to illustrate therelationship between the energies and the weight functions.Figure 3: The learning curve of the model.16igure 4: The SDFs for the exciton energy J ( ω ) ( = J ( ω ) ) (blue) and the interactionenergy J ( ω ) (orange) obtained from the machine learning approach.17 δω ( t ) δω ( t ) δ∆ ( t ) + - Ω + Ω - Figure 5: Energy level diagram of a dimer system that undergoes the random fluctuationsof the excited energy and coupling strength described by δω ( t ) and δ ∆( t ) , respectively.For the description of pure dephasing, only the difference of the energies involved in theoptical excitation is important: The frequency fluctuation between | (cid:105) and | (cid:105) is given by δω ( t ) + δ ∆( t ) , whereas | (cid:105) and | (cid:105) is given by δω ( t ) − δ ∆( t ) . For perfectly uncorrelatedfluctuations, we treat δω ( t ) and δ ∆( t ) independently.18ot easy to separate the effects of exciton transfer from the fluctuation of energies arisesfrom the environmental motions. By introducing the localization weight functions w ± ,m ( t ) and w ∆ ( t ) in Eqs. (21) and (22) to eliminate the effects of non-environmental origin in thelearning trajectories, we can stabilize and enhance the efficiency of the machine learningprocess.In Fig. 3, we depict the learning curve of the loss function defined as Eq. (24). Afterrandom samplings of φ , the loss function converged monotonically to a certain positive value,which demonstrated the efficiency of the present algorithm. The initial parameter values ofthe excitation energy and the interaction energy were set by ω = 17736 cm − and ∆ =1004 cm − , whereas the optimized values of the excitation energy and the interaction energywere given by ω = 18060 cm − and ∆ = 845.0 cm − , which are closer to the experimentalobtained values.In Fig. 4, we display the results of the SDFs for the excitation energy J ( ω ) ( = J ( ω ) )and the interaction energy J ( ω ) . Various intermolecular modes below 2000 cm − are ob-served as the prominent sharp peaks near 446, 562, 1530, 1790, 1840, and 1920 cm − . In theregion over 2000 cm − , only two tiny peaks near 3000 cm − and 3840 cm − are observed. Thenormal mode analysis (B3LYP/def-SV(P)) indicates these peaks under 3300 − arise fromthe intramolecular modes of the indocarbocyanine dimer, whereas the peak at 3850 cm − arises from a molecular vibration of solvent methanol molecules. We found that each sharppeak can be fitted by the Brownian spectral distribution, whereas the broadened back-ground peak in the range from 0 to 2000cm − corresponding to the intramolecular modes isfitted by the Drude-Lorentz distribution. The intensities of the peaks in J ( ω ) are consid-erably weaker than those in J ( ω ) : Only the peaks near 456, 562, 1840 and 1920 cm − areidentified. As we expected, the positions of intermolecular peaks are determined from theclassical MD, whereas the heights of these are predominately determined from the quantumchemistry calculation.To verify the description of obtained SDFs and system parameters, we computed lin-19ar absorption and two-dimensional electronic spectra (2DES), where the experimentallyobtained spectra are available. In general, these spectra should be calculated in the frame-work of open quantum dynamics taking into account the complex interaction between theexciton sites. For a demonstration purpose, however, here we employ the analytical ex-pressions of response functions ignoring the transition to the singlet state that is usuallyforbidden. The details of calculations are presented in Appendix A.Linear absorption spectrum calculated from Eqs. (A.1 ) and (A.2 ) is presented in Fig.6. Here, the Gaussian spectral peak, λ exp (cid:2) − (( ω − ω c ) /γ ) (cid:3) , was fitted on the vicinity ofthe peak, where the the amplitude, central frequency, and width were chosen as λ = 348 . , ω c = 18750cm − , and γ = 466cm − . Note that because we did not try to determined theamplitude of dipole operator for simplicity, we could not determine the intensities of SDFs:We scale the intensities of J ( ω ) to fit the experimentally obtained signal. As presented inFig. 6, we observe the single broadened absorption peak at ω + ∆ corresponding to thetransition between | (cid:105) and | (cid:105) , while the transition between | (cid:105) and | −(cid:105) is forbidden.(see Fig. 5). Although the experimentally observed linear absorption spectrum exhibit − phonon sideband peak near ω = 19500cm − , here we only observe it as the asymmetry ofthe Gaussian peak in the high frequency region.The 2D correlation electronic spectra calculated using the analytical expressions of theresponse function (Eqs. (A.5 ) -(A.7 )) are presented in Fig 7. At t = 0 fs only onepeak stretched near the ω = ω line arising from the | (cid:105) → | (cid:105) transition is observed.At t = 20 fs, then the peak is elongated to the low-frequency ω direction due to shift ofthe eigenenergy caused by the heat-bath induced exciton-exciton interaction described by J ( ω ) . Because the system-bath interaction here we considered is non-Markovian whoseeffects appear only after the period longer than the inverse correlation time of noise, wedo not observe such heat-bath effects for small t . Then, around t = 40 fs, the off-diagonal peak near ( ω , ω ) = (17800 , in units of cm − corresponding to the tran-sition between | (cid:105) → | (cid:105) is observed, whereas the peak along the ω = ω line shifts to20igure 6: Linear absorption spectrum calculated from Eqs. (A.1 ) and (A.2 ) with the lineshape function Eq. (A.3 ) for the system parameters and SDFs obtained from the machinelearning approach. The dotted line is the fitted Gaussian peak centered at − ,indicating that the calculated peak is asymmetric due to the − phonon transition near − . ( ω , ω ) = (21000 , due to the transition | (cid:105) → | (cid:105) arises from the exciton-excitoninteraction described by ∆ and J ( ω ) . As t goes, the intensities of these two peaks changeoscillatory, as the results of the population transition among | (cid:105) , | (cid:105) , and | (cid:105) caused by ∆ and J ( ω ) . This phenomena was also observed experimentally. The appearance of thisoscillatory feature at finite period in t indicates the importance of the off-diagonal heat-bath, whose modeling is not easy in the framework of the existing approach. At t = 300 fs the peak elongated in the ω = ω direction due to the inhomogeneous broadening arisesfrom the baths described as J ( ω ) and J ( ω ) . We introduced a machine leaning approach to constructing a model to analyze the dy-namics of exciton or electron transfer processes in a complex environment on the basis ofthe energy eigenstates evaluated from QM/MM simulations as a function of time. The keyfeature of the present study is a system-bath model, in which the primary exciton/electron21igure 7: 2DES calculated from Eqs. (A.5 ) -(A.7 ) with the line shape function Eq. (A.3 )for the system parameters and SDFs obtained from the machine learning approach. Thewaiting time t for each signal is displayed on the left top of each pane. The peak intensityof the signal was normalized for each t . 22ynamics are described by a system Hamiltonian expressed in terms of discretized energystates, while the other degrees of freedom are described by a harmonic heat-bath that ischaracterized by SDFs. An optimized system-bath Hamiltonian obtained from the machinelearning approach allows us to conduct time-irreversible quantum simulations not feasiblefrom a full quantum MD simulation approach.Here, we demonstrated the above features by calculating linear and nonlinear opticalspectra for the indocarbocyanine dimer system in the methanol environment. in which thequantum entanglement between the system and bath plays a central role. The calculatedresults explain the experimental results reasonably well: We found that the heat-bath forthe exciton-exciton interaction plays a key role for descripting exciton transfer process inthis system. Here, we ignore the transition between the singlet and triple states due to anapplicability of the analytical expression, if necessary, we can treat it explicitly using theHEOM formalism.
Finally, we briefly discuss some extensions of this study. As shown in the previous pa-per, the machine learning approach can be applied not only for a system described byelectronic states but also a system described by reaction coordinates, which is useful forinvestigating chemical reaction processes characterized by potential energy surfaces. Bycombining the previous and present approaches, we can further investigate a system de-scribed as not only electronic states but also molecular configuration space in a frameworkof the system-bath model, for example, photoisomerization, molecular motor, and non-adiabatic transition problems. In this way, we may construct a system-bath model for en-tire photosynthesis reaction processes consisting of photoexcitation, exciton transfer, electron transfer, and proton transfer processes, including the conversion processes,such as exciton-coupled electron transfer and electron coupled proton transfer processes. Further theoretical and computational efforts that include providing learning dates fromaccurate and large quantum simulations, improving learning algorithms, and developingaccurate and efficient open quantum dynamics theory to treat a complex system-bath model23ust be made. We leave such extensions to future studies, in accordance with progress intheoretical techniques.
Acknowledgement
The authors are thankful to professor Yuki Kurashige for helpful discussions concerningthe QM/MM simulations for an indocarbocyanine dimer system. Financial support fromHPC Systems Inc. is acknowledged.
A One-dimensional and two-dimensional spectra
Linear and nonlinear optical spectra can be expressed in the Fourier transformationform of the response functions. In the present dimer case, we can analytically expressthe response functions in terms of a line shape function including the contribution from anexiciton-exiciton interaction. The linear absorption spectrum is given by S ( ω ) = (cid:90) ∞ dte iωt R (1) ( t ) − c.c., (A.1)where R (1) ( t ) = (cid:104) [ µ ( t ) , µ (0)] (cid:105) is the 1D response function expressed in terms of the transitiondipole moment µ ( t ) . For a coupled dimmer system, the analytical expression for the responsefunction for Eq.(A.1 ) is expressed as R (1) ( t ) = iµ (cid:126) exp[ i Ω + t − g − ( t ) − g − ( t )] − c.c. (A.2)24here the line shape function, g a ± ( t ) , for the SDF, J a ( ω ) , with a = 11 and 12 is given by g a ± ( t ) ≡ (cid:90) t d (cid:48) t (cid:90) t (cid:48) dt (cid:90) dω π × J a ( ω ) (cid:20) coth (cid:18) β (cid:126) ω (cid:19) cos( ωt ) ± i sin( ωt ) (cid:21) , (A.3) ttt (v) (vi) (vii) (viii) ttt (i) (ii) (iii) (iv) Figure A.0: Double-sided Feynman diagrams for the third-order response functions R (3) ( t , t , t ) . In each diagram, time runs from bottom to top, and t i represents the timeintervals for the i th sequence between successive laser–system interactions. The left linerepresents the time evolution of the ket, whereas the right line represents that of the bra.The black, red and blue lines represent that the system is in the | (cid:105) , | (cid:105) , and | (cid:105) states,respectively. The complex conjugate paths of these, which can be obtained by interchangingthe ket and bra diagrams, are not shown here.For the coupled dimmer system, the third-order response function R (3) ( t , t , t ) = (cid:104) [ µ ( t ) , [ µ ( t ) [ µ ( t ) , µ (0)]]] (cid:105) (A.4)25s also evaluated in the analytical form as R (3) ( t , t , t ) ≡ i µ (cid:126) (cid:88) α =1 exp[ Q α (t)] − c.c. (A.5)where Q ( t ) = − i Ω + t − i Ω + t − f ( t , t , t ) − f ( t , t , t ) ,Q ( t ) = i Ω + t − i Ω + t − f ( t , t , t ) − f ( t , t , t ) ,Q ( t ) = i Ω + t − i Ω + t − f ( t , t , t ) − f ( t , t , t ) ,Q ( t ) = − i Ω + t − i Ω + t − f ( t , t , t ) − f ( t , t , t ) ,Q ( t ) = − i Ω + t + i Ω − t − f ( t , t , t ) − f ( t , t , t ) ,Q ( t ) = i Ω + t + i Ω − t − f ( t , t , t ) − f ( t , t , t ) ,Q ( t ) = i Ω + t + 2 iω t + i Ω − t − f ( t , t , t ) − f ( t , t , t ) ,Q ( t ) = − i Ω + t − iω t − i Ω + t − f ( t , t , t ) − f ( t , t , t ) (A.6)26ith f a ( t , t , t ) = − g a − ( t ) − g a + ( t ) − [ g a + ( t ) − g a + ( t ) − g a − ( t ) + g a − ( t )] ,f a ( t , t , t ) = − g a + ( t ) − g a + ( t )+ [ g a − ( t ) − g a − ( t ) − g a + ( t ) + g a + ( t )] ,f a ( t , t , t ) = − g a + ( t ) − g a − ( t )+ [ g a + ( t ) − g a + ( t ) − g a + ( t ) + g a + ( t )] ,f a ( t , t , t ) = − g a − ( t ) − g a − ( t ) − [ g a − ( t ) − g a − ( t ) − g a − ( t ) + g a − ( t )] ,f a ( t , t , t ) = − g a − ( t ) − g a + ( t ) − [ g a + ( t ) − g a + ( t ) − g a − ( t ) + g a − ( t )] ,f a ( t , t , t ) = − g a + ( t ) − g a + ( t )+ [ g a − ( t ) − g a − ( t ) − g a + ( t ) + g a + ( t )] ,f a ( t , t , t ) = − g a + ( t ) − g a − ( t )+ [ g a + ( t ) − g a + ( t ) − g a + ( t ) + g a + ( t )] ,f a ( t , t , t ) = − g a − ( t ) − g a − ( t ) − [ g a − ( t ) − g a − ( t ) − g a − ( t ) + g a − ( t )] . (A.7)Here, t ≡ t + t , t ≡ t + t , and t ≡ t + t + t . Because the fluctuation in the | (cid:105) and | (cid:105) states are described by J ( ω ) + J ( ω ) and J ( ω ) + J ( ω ) = 2 J ( ω ) , respectively,the line shape function g a ± ( t ) in Eq.A.7 is now expressed as g ± ( t ) = g ± ( t ) + g ± ( t ) , and g ± ( t ) = 2 g ± ( t ) . By using the third-order diagrams, the pump–probe spectrum and photonecho spectrum are, for example, calculated from the Q ( t ) , Q ( t ) , and Q ( t ) elements, andthe Q ( t ) , Q ( t ) , and Q ( t ) elements, respectively.Although the change of exciton population can be explored by pump-probe spectroscopy,27f we wish to investigate not only population dynamics but also the system-bath coherence,two-dimensional electronic correlation spectroscopy is a better choice. This spectrum can becalculated from I (corr) ( ω , t , ω )= I (NR) ( ω , t , ω ) + I (R) ( ω , t , ω ) , (A.8)where the non-rephasing and rephasing parts of the signal are defined by I (NR) ( ω , t , ω ) =Im (cid:90) ∞ d t (cid:90) ∞ d t e iω t e iω t R (3) ( t , t , t ) , (A.9)and I (R) ( ω , t , ω ) =Im (cid:90) ∞ d t (cid:90) ∞ d t e iω t e − iω t R (3) ( t , t , t ) . (A.10)28 eferences (1) May, V.; Kühn, O. Charge and Energy Transfer Dynamics in Molecular Systems ; Wiley-VCH, 2003.(2) Schröter, M.; Ivanov, S.; Schulze, J.; Polyutov, S.; Yan, Y.; Pullerits, T.; Kühn, O.Exciton–vibrational coupling in the dynamics and spectroscopy of Frenkel excitons inmolecular aggregates.
Physics Reports , , 1–78.(3) Kreisbeck, C.; Kramer, T.; Aspuru-Guzik, A. Scalable High-Performance Algorithm forthe Simulation of Exciton Dynamics. Application to the Light-Harvesting Complex IIin the Presence of Resonant Vibrational Modes. J. Chem. Theory Comput. , ,4045–4054.(4) Lee, M. K.; Coker, D. F. Modeling Electronic-Nuclear Interactions for Excitation En-ergy Transfer Processes in Light-Harvesting Complexes. J. Phys. Chem. Lett. , ,3171–3178.(5) Kreisbeck, C.; Kramer, T. Long-Lived Electronic Coherence in Dissipative ExcitonDynamics of Light-Harvesting Complexes. J. Phys. Chem. Lett. , , 2828–2833.(6) Strümpfer, J.; Schulten, K. Light harvesting complex II B850 excitation dynamics. J.Chem. Phys. , , 225101.(7) Strümpfer, J.; Schulten, K. The effect of correlated bath fluctuations on exciton transfer. J. Chem. Phys. , , 095102.(8) Strümpfer, J.; Schulten, K. Excited state dynamics in photosynthetic reaction centerand light harvesting complex 1. J. Chem. Phys. , , 065101.(9) Olbrich, C.; Strümpfer, J.; Schulten, K.; Kleinekathöfer, U. Quest for Spatially Cor-related Fluctuations in the FMO Light-Harvesting Complex. J. Phys. Chem. B , , 758–764. 2910) Ishizaki, A.; Fleming, G. R. On the adequacy of the Redfield equation and relatedapproaches to the study of quantum dynamics in electronic energy transfer. J. Chem.Phys. , , 234110.(11) Raszewski, G.; Saenger, W.; Renger, T. Theory of Optical Spectra of PhotosystemII Reaction Centers: Location of the Triplet State and the Identity of the PrimaryElectron Donor. Biophys. J. , , 986–998.(12) Renger, T.; Klinger, A.; Steinecker, F.; Schmidt am Busch, M.; Numata, J.; Müh, F.Normal Mode Analysis of the Spectral Density of the Fenna–Matthews–Olson Light-Harvesting Protein: How the Protein Dissipates the Excess Energy of Excitons. J. Phys.Chem. B , , 14565–14580.(13) Müh, F.; Plöckinger, M.; Ortmayer, H.; Schmidt am Busch, M.; Lindorfer, D.;Adolphs, J.; Renger, T. The quest for energy traps in the CP43 antenna of photo-system II. J. Photochem. Photobiol., B , , 286–300.(14) Müh, F.; Plöckinger, M.; Renger, T. Electrostatic Asymmetry in the Reaction Centerof Photosystem II. J. Phys. Chem. Lett. , , 850–858.(15) Gelzinis, A.; Abramavicius, D.; Ogilvie, J. P.; Valkunas, L. Spectroscopic properties ofphotosystem II reaction center revisited. J. Chem. Phys. , , 115102.(16) Fujihashi, Y.; Fleming, G. R.; Ishizaki, A. Impact of environmentally induced fluctu-ations on quantum mechanically mixed electronic and vibrational pigment states inphotosynthetic energy transfer and 2D electronic spectra. J. Chem. Phys. , ,212403.(17) Novoderezhkin, V. I.; Romero, E.; Dekker, J. P.; van Grondelle, R. Multiple Charge-Separation Pathways in Photosystem II: Modeling of Transient Absorption Kinetics. Phys. Chem. Chem. Phys. , , 681–688.3018) Novoderezhkin, V. I.; Romero, E.; van Grondelle, R. How exciton-vibrational coherencescontrol charge separation in the photosystem II reaction center. Phys. Chem. Chem.Phys. , , 30828–30841.(19) Gelzinis, A.; Valkunas, L.; Fuller, F. D.; Ogilvie, J. P.; Mukamel, S.; Abramavi-cius, D. Tight-binding model of the photosystem II reaction center: application totwo-dimensional electronic spectroscopy. New J. Phys. , , 075013.(20) Adolphs, J.; Renger, T. How Proteins Trigger Excitation Energy Transfer in the FMOComplex of Green Sulfur Bacteria. Biophys. J. , , 2778–2797.(21) Garg, A.; Onuchic, J. N.; Ambegaokar, V. Effect of friction on electron transfer inbiomolecules. J. Chem. Phys. , , 4491.(22) Wolynes, P. G. Dissipation, tunneling, and adiabaticity criteria for curve crossing prob-lems in the condensed phase. J. Chem. Phys. , , 1957.(23) Yan, Y. J.; Sparpaglione, M.; Mukamel, S. Solvation dynamics in electron-transfer,isomerization, and nonlinear optical processes: a unified Liouville-space theory. J. Phys.Chem. , , 4842–4853.(24) Tanaka, M.; Tanimura, Y. Quantum Dissipative Dynamics of Electron Transfer Reac-tion System: Nonperturbative Hierarchy Equations Approach. J. Phys. Soc. Jpn. , , 073802.(25) Tanaka, M.; Tanimura, Y. Multistate electron transfer dynamics in the condensedphase: Exact calculations from the reduced hierarchy equations of motion approach. J.Chem. Phys. , , 214502.(26) Sim, E.; Makri, N. Path Integral Simulation of Charge Transfer Dynamics in Photo-synthetic Reaction Centers. J. Phys. Chem. B , , 5446–5458.3127) Sim, E. Determination of the Electron Transfer Mechanism through Decomposition ofthe Density Matrix. J. Phys. Chem. B , , 19093–19095.(28) Dijkstra, A. G.; Tanimura, Y. Correlated fluctuations in the exciton dynamics andspectroscopy of DNA. New J. Phys. , , 055005.(29) Gélinas, S.; Rao, A.; Kumar, A.; Smith, S. L.; Chin, A. W.; Clark, J. Ultrafast Long-Range Charge Separation in Organic Semiconductor Photovoltaic Diodes. , ,6.(30) Zirzlmeier, J.; Lehnherr, D.; Coto, P. B.; Chernick, E. T.; Casillas, R.; Basel, B. S.;Thoss, M.; Tykwinski, R. R.; Guldi, D. M. Singlet fission in pentacene dimers. Proc.Natl. Acad. Sci. U.S.A. , , 5325–5330.(31) Tamura, H.; Burghardt, I. Ultrafast Charge Separation in Organic Photovoltaics En-hanced by Charge Delocalization and Vibronically Hot Exciton Dissociation. J. Am.Chem. Soc. , , 16364–16367.(32) Huix-Rotllant, M.; Tamura, H.; Burghardt, I. Concurrent Effects of Delocalization andInternal Conversion Tune Charge Separation at Regioregular Polythiophene–FullereneHeterojunctions. J. Phys. Chem. Lett. , , 1702–1708.(33) Oviedo-Casado, S.; Urbina, A.; Prior, J. Magnetic field enhancement of organic photo-voltaic cells performance. Sci. Rep. , , 4297.(34) Cainelli, M.; Tanimura, Y. Exciton Transfer in Organic Photovoltaic Cells: A Role ofLocal and Nonlocal Electron-Phonon Interactions in a Donor Domain. J. Chem. Phys. , , 034107.(35) Tanimura, Y. Stochastic Liouville, Langevin, Fokker–Planck, and Master Equation Ap-proaches to Quantum Dissipative Systems. J. Phys. Soc. Jpn. , , 082001.3236) Tanimura, Y. Numerically “exact” approach to open quantum dynamics: The hierar-chical equations of motion (HEOM). J. Chem. Phys. , , 020901.(37) Dral, P. O.; Barbatti, M.; Thiel, W. Nonadiabatic Excited-State Dynamics with Ma-chine Learning. J. Phys. Chem. Lett. , , 5660–5663.(38) Hartmann, M. J.; Carleo, G. Neural-Network Approach to Dissipative Quantum Many-Body Dynamics. Phys. Rev. Lett. , , 250502.(39) Flurin, E.; Martin, L.; Hacohen-Gourgy, S.; Siddiqi, I. Using a Recurrent Neural Net-work to Reconstruct Quantum Dynamics of a Superconducting Qubit from PhysicalObservations. Phys. Rev. X , , 011006.(40) Smith, J. S.; Isayev, O.; Roitberg, A. E. ANI-1: an extensible neural network potentialwith DFT accuracy at force field computational cost. Chem. Sci. , , 3192–3203.(41) Sauceda, H. E.; Chmiela, S.; Poltavsky, I.; Müller, K.-R.; Tkatchenko, A. Molecularforce fields with gradient-domain machine learning: Construction and application todynamics of small molecules with coupled cluster forces. J. Chem. Phys. , ,114102.(42) Sifain, A. E.; Lubbers, N.; Nebgen, B. T.; Smith, J. S.; Lokhov, A. Y.; Isayev, O.;Roitberg, A. E.; Barros, K.; Tretiak, S. Discovering a Transferable Charge AssignmentModel Using Machine Learning. J. Phys. Chem. Lett. , , 4495–4501.(43) Chmiela, S.; Tkatchenko, A.; Sauceda, H. E.; Poltavsky, I.; Schütt, K. T.; Müller, K.-R.Machine learning of accurate energy-conserving molecular force fields. Sci. Adv. , , e1603015.(44) Chmiela, S.; Sauceda, H. E.; Müller, K.-R.; Tkatchenko, A. Towards exact moleculardynamics simulations with machine-learned force fields. Nature Comm. , , 1–10.(45) Kampen, N. V. Stochastic Processes in Physics and Chemistry ; Elsevier, 1981.3346) Chen, L.; Zhao, Y.; Tanimura, Y. Dynamics of a One-Dimensional Holstein Polaronwith the Hierarchical Equations of Motion Approach.
J. Phys. Chem. Lett. , ,3110–3115.(47) Dunn, I. S.; Tempelaar, R.; Reichman, D. R. Removing instabilities in the hierarchicalequations of motion: Exact and approximate projection approaches. J. Chem. Phys. , , 184109.(48) Ueno, S.; Tanimura, Y. Modeling and simulating excited-state dynamics of a system incondensed phases: Machine Learning approach. J. Chem. Theory Comput. , ,2099–2108.(49) Halpin, A.; Johnson, P. J. M.; Tempelaar, R.; Murphy, R. S.; Knoester, J.; Jansen, T.L. C.; Miller, R. J. D. Two-dimensional spectroscopy of a molecular dimer unveils theeffects of vibronic coupling on exciton coherences. Nature Chem. , , 196–201.(50) Berendsen, H.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passingparallel molecular dynamics implementation. Comput. Phys. Comm. , , 43–56.(51) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E.GROMACS: High performance molecular simulations through multi-level parallelismfrom laptops to supercomputers. SoftwareX , , 19–25.(52) Bekker, H.; Berendsen, H.; Dijkstra, E.; Achterop, S.; Van Drunen, R.; Van derSpoel, D.; Sijbers, A.; Keegstra, H.; Reitsma, B.; Renardus, M. Gromacs: A parallelcomputer for molecular dynamics simulations. Physics computing. 1993; pp 252–256.(53) Ridley, J.; Zerner, M. An intermediate neglect of differential overlap technique forspectroscopy: Pyrrole and the azines. Theoret. Chim. Acta , , 111–134.(54) Thompson, M. A.; Zerner, M. C. A theoretical examination of the electronic structure34nd spectroscopy of the photosynthetic reaction center from Rhodopseudomonas viridis. J. Am. Chem. Soc. , , 8210–8215.(55) Martin, R. L. Natural transition orbitals. J. Chem. Phys. , , 4775.(56) Neese, F. The ORCA program system. WIREs Comput. Mol. Sci. , , 73–78.(57) Abadi, M.; Agarwal, A.; Barham, P.; Brevdo, E.; Chen, Z.; Citro, C.; Corrado, G. S.;Davis, A.; Dean, J.; Devin, M.; Ghemawat, S.; Goodfellow, I.; Harp, A.; Irv-ing, G.; Isard, M.; Jia, Y.; Jozefowicz, R.; Kaiser, L.; Kudlur, M.; Levenberg, J.;Mane, D.; Monga, R.; Moore, S.; Murray, D.; Olah, C.; Schuster, M.; Shlens, J.;Steiner, B.; Sutskever, I.; Talwar, K.; Tucker, P.; Vanhoucke, V.; Vasudevan, V.;Viegas, F.; Vinyals, O.; Warden, P.; Wattenberg, M.; Wicke, M.; Yu, Y.; Zheng, X.TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems. arXiv:1603.04467 [cs] ,(58) Tanimura, Y. Reduced hierarchy equations of motion approach with Drude plus Brow-nian spectral distribution: Probing electron transfer processes by means of two-dimensional correlation spectroscopy. J. Chem. Phys. , , 22A550.(59) Tanimura, Y.; Mukamel, S. Femtosecond pump–probe spectroscopy of intermolecularvibrations in molecular dimers. J. Chem. Phys. , , 1981–1984.(60) Ikeda, T.; Tanimura, Y. Probing photoisomerization processes by means of multi-dimensional electronic spectroscopy: The multi-state quantum hierarchical Fokker-Planck equation approach. The Journal of Chemical Physics , , 014102.(61) Ikeda, T.; Dijkstra, A.; Tanimura, Y. Modeling and analyzing a photo-driven molecularmotor system: Ratchet dynamics and non-linear optical spectra. J. Chem. Phys. , , 114103. 3562) Ikeda, T.; Tanimura, Y. Phase-space wavepacket dynamics of internal conversion viaconical intersection: Multi-state quantum Fokker-Planck equation approach. Chem.Phys. , , 203.(63) Chen, L.; Shi, Q. Quantum rate dynamics for proton transfer reactions in condensedphase: The exact hierarchical equations of motion approach. J. Chem. Phys. , , 134505.(64) Shi, Q.; Zhu, L.; Chen, L. Quantum rate dynamics for proton transfer reaction in amodel system: Effect of the rate promoting vibrational mode. J. Chem. Phys. , , 044505.(65) Zhang, J.; Borrelli, R.; Tanimura, Y. Proton tunneling in a two-dimensional poten-tial energy surface with a non-linear system–bath interaction: Thermal suppression ofreaction rate. J. Chem. Phys. , , 214114.(66) Sakamoto, S.; Tanimura, Y. Exciton-Coupled Electron Transfer Process Controlled byNon-Markovian Environments. J. Phys. Chem. Lett. , , 5390–5394.(67) Zhang, J.; Borrelli, R.; Tanimura, Y. Probing photoinduced proton coupled electrontransfer process by means of two-dimensional electronic-vibrational spectroscopy. J.Chem. Phys. , , xxx.(68) Mukamel, S. Principles of Nonlinear Optical Spectroscopy ; Oxford University Press,1999.(69) Tanimura, Y.; Mukamel, S. Real-time path-integral approach to quantum coherenceand dephasing in nonadiabatic transitions and nonlinear optical response.
Phys. Rev.E , , 118–136. 36 raphical TOC Entryraphical TOC Entry