Multisurface Multimode Molecular Dynamical Simulation of

Article pubs.acs.org/JPCA

Multisurface Multimode Molecular Dynamical Simulation of Naphthalene and Anthracene Radical Cations by Using Nearly Linear Scalable Time-Dependent Discrete Variable Representation Method Basir Ahamed Khan,†,‡ Subhankar Sardar,*,§ Pranab Sarkar,‡ and Satrajit Adhikari§ †

Department of Physics, Krishnath College, Berhampore, West Bengal 742101, India Department of Chemistry, Visva-Bharati University, Santiniketan, West Bengal 731235, India § Department of Physical Chemistry, Indian Association for the Cultivation of Science, Jadavpur, Kolkata 700032, India ‡

S Supporting Information *

ABSTRACT: The major portion of the algorithm of the time-dependent discrete variable representation (TDDVR) method is recently parallelized using the shared-memory parallelization scheme with the aim of performing dynamics on relatively large molecular systems. Because of the astronomical importance of naphthalene and anthracene, we have investigated their radical cations as models for theoretical simulation of complex photoelectron spectra and nonradiative decay process using the newly implemented parallel TDDVR code. The strong vibronic coupling among the six lowest doublet electronic states makes these polynuclear hydrocarbons dynamically important. The aim of the present investigation is to show the efficiency of our current TDDVR algorithm to perform dynamics on large dimensional quantum systems in vibronically coupled electronic manifold. Both the sequential and the parallelized TDDVR algorithms are almost linear scalable for an increase in number of processors. Because a significant speed-up is achieved by cycling in the correct way over arrays, all of the simulations are performed within a reasonable wall clock time. Our theoretical spectra well reproduce the features of the corresponding experimental analog. The dynamical outcomes, for example, population, photoelectron spectra, and diffused interstellar bands, etc., of our quantumclassical approach show good agreement with the findings of the well-established quantum dynamical method, that is, multi configuration time-dependent Hartree (MCTDH) approach.

1. INTRODUCTION Femtosecond chemistry1−3 has received much attention in recent years due to its widespread applicability in molecular physics, spectroscopy, chemistry, and biology. The availability of powerful lasers has stimulated the design of new experimental techniques that make it possible to analyze femtosecond internal conversion processes and radiationless transitions. Theoreticians simultaneously developed a powerful molecular dynamics algorithm4−9 to treat this fundamental process in an exact quantum or mixed quantum mechanical way. The molecular quantum dynamics become interesting in the presence of conical intersection (CI),10−12 which is essentially responsible for the ultrafast decay processes. Such nonadiabatic or non Born−Oppenheimer (BO) phenomena13−22 are very frequent in polyatomics and need special attention for accurate quantum treatment. The large molecular systems possessing many degrees of freedom (DOFs) can be handled in a mixed quantum-classical way because all of them are not significant in dynamics. This feature is utilized in our first principle-based quantum-classical method, the timedependent discrete variable representation (TDDVR) method. The method has the following important dynamical features: (i) The TDDVR method has the capability to distinguish the quantum and classical regime of a molecular system, and thereby it treats the dynamically predominating modes © XXXX American Chemical Society

quantum mechanically and the less important modes classically with minor quantum correction as per requirement. (ii) Because the grid-points involved in dynamics are moving, the required number of basis functions is less, and they achieve convergence significantly fast. (iii) The different vibrational modes independently contribute to the TDDVR equation of motion (EOM) that facilitates the parallelization of the algorithm. The conical intersections, frequently present in polyatomics, among the potential energy surfaces (PESs) are the origin of nonadiabatic or vibronic coupling, and the role of such coupling in molecular dynamics is an important research topic10−12,14 in modern chemistry. Naphthalene and anthracene appeared to be the challenging molecules where complicated electronic nonadiabatic features prevail. Over the last few decades, both the experimentalists23−31 and the theoreticians32−38 have become very much interested in the photophysics of these radical cations, because the diffused interstellar bands (DIBs) of star Cernis 5230 and HD2811591731 are coincident within the strongest bands of these polynuclear hydrocarbons (PNHs). The origin and assignment of these DIBs are two issues that Received: July 25, 2014 Revised: November 15, 2014

A

dx.doi.org/10.1021/jp507459m | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

orthogonal but not normalized, which should be normalized before the initiation of dynamics. For rth DOF, the mrth basis function is represented in TDDVR scheme as a product of DVR basis function (Ξη(zr)) and plane wave (χ(qr, t)) as follows:

have fundamentally existed for a long time in astrophysical investigation.29−31 The removal of a valence electron from neutral naphthalene produces its doublet radical cation (Nap+) with the following six low-lying electronic states: D0 (X2Au), D1 (A2B3u), D2 (B2B2g), D3 (C2B1g), D4 (D2Ag), and D5 (E2B3g), which are vibronically coupled through 48 vibrational modes. In a similar fashion, the six lowest doublet PESs, X2B2g, A2B1g, B2Au, C2B2g, D2B3u, and E2Ag, of anthracene radical cation (Ant+) are entangled with 66 vibrational DOFs. As the electronic states of both PNHs are vibronically coupled, the BO theory fails to describe the underlying dynamics. The coupling terms in adiabatic condition are singular at the interaction region, and thus usage of diabatic representation of Hamiltonian is essential. Because the prime interest of our effort is to implement and establish an accurate and efficient parallel molecular dynamics method, the well-known Hamiltonians of both PNHs developed by Ghanta et al.33 are used as models to reproduce the molecular properties. The experimental and theoretical studies of the aromatic polycyclic hydrocarbon cations are most prevalent in current research because of their existence in the interstellar medium (ISM) showing the DIBs. Such DIBs originated due to complex entanglement of the electronic surfaces, which are incorporated in these realistic model Hamiltonians.33 As the quantum phenomena are predominating among the coupled PESs, a rigorous quantum mechanical study is always challenging for any theoretician, and thus one of the purposes of our present numerical effort is the theoretical simulation of these DIBs. The different sections of this Article are organized as follows: The theoretical background of the TDDVR method is briefly presented in section 2 with a few important equations. The next section elaborately describes the recent development of the parallel TDDVR algorithm. The computational advantages of the parallel TDDVR methodology in terms of speed-up and scalibility are illustrated in section 4. In section 5, we have discussed in detail the outcomes of our theoretical simulation in comparison with the results obtained by MCTDH method35,36 and experimental measurements.23−25,27 The few important features of the parallel TDDVR algorithm and its application to two real molecules are summarized in the final section.

ηr

Φmr (qr , t ) =

m1m2...mp

qm (t ) = r

ℏ zmr + qrc(t ) 2 Im A r

(3)

The width parameter (Ar) of the wavepacket is considered here as time independent for any mode (rth mode) to avoid stiffness that arises due to classical propagation. If the system Hamiltonians33 adopted for dynamics and the wave function (WF) expanded in TDDVR scheme (Ψk({qr}, t)) are placed in the time-dependent Schrödinger equation (TDSE), we obtain the quantum equation of motion (EOM) as presented below. The amplitude vector of the WF (am1m2···mp) for any PES is related to the TDDVR expansion coefficients by the relation, β m 1 m 2 ···m p = a m 1 m 2 ···m p ∏ r =p 1 (S m r m r ) 1/2 , where S m r ,m r ′ = ∑ η =ηr 0 Ξ η* (z m r )Ξ η (z m r ′). The time-dependent amplitude (β̇m1m2···mp,k) of the WF, which we have to evaluate at each time step, in the kth PES is ⎡

βṁ m .... m , k = 1 2

p

∑ ⎢⎢ ∑ r

⎧ 2 ⎪ iq ̇ −⎨ r + ⎪ 2ℏ ⎩ +

Im A r Tm̅ r mr′

2i ⎣ m1′m2′..m′p ⎫ ipq̇ c Fm̅ r mr ⎪ r ⎬βm m ..m , k 1 2 p 4ℏIm A r ⎪ ⎭

1 {Vkk(m1m2.... mp)βm m .... m , k + 1 2 p iℏ ⎤ βm m .... m , k ′}⎥ 1 2 p ⎥⎦

p

βm ′m ′..m′ , k ∏ δmr′mr′′ 1

2

p

r ′≠ r

∑ Vkk′(m1m2.... mp) k ′≠ k

(4)

The TDDVR formulation has the scope to represent all of the electronic states (k) properly, and the corresponding EOM can transfer the information of one surface to others when the wavepacket moves from upper to lower PESs. The explicit form of the kinetic matrices used in the TDDVR quantum EOM (eq 4) is presented in the Supporting Information, which have several important dynamical features. The component matrices (see eqs A1a,A1b of the Supporting Information) of eq 4 are required to evaluate once at the start of the dynamics as they are independent of time. Along with the above-mentioned

am1m2.... mp , k (t ) ∏ Φmr (qr , t ) r=1

(2)

The DVR functions are constructed with the harmonic oscillator eigen functions as primitive basis. The Gauss− Hermite (GH) basis (ψη(qr, t)) for rth normal mode is orthonormal, and its ground state is called a Gaussian wave packet (GWP). The grid-point, qmr, of the TDDVR basis function is prepared with the roots (zm) of the ηrth Hermite polynomial. It is worthwhile to mention that although the root points (zmr) are fixed, the grid-points of the TDDVR method move on the configuration space with time propagation due to the time-dependent parameter, qcr (t):

p

∑

r

η=0

2. THEORETICAL FRAMEWORK OF THE TDDVR METHODOLOGY The different versions of the TDDVR formulation were presented in detail in our last few publications,39−41 and therefore we present here only the important equations of the latest version of the method. In case of Nap+ and Ant+, the TDDVR scheme propagates the wavepacket (WP) on a six coupled electronic manifold, where the time-dependent DVR grid-points are used to propagate the WP with fixed width parameter. The diabatic wave function (Ψk({qr}, t)) for kth PES with total p number of vibrational modes can be expanded in product form with basis functions [{Φmr (qr, t)}] in TDDVR scheme as Ψk({qr }, t ) =

∑ Ξ*η (zm )ψη(qr , t )

(1)

The individual state wave function (Φr({qr}, t)) for different surfaces produces the system wave function (Ψk({qr}, t)). The set of TDDVR basis functions, Φm, for the rth mode is B



Article

for different (kth number) PESs. The delta function, δmr′mr′′ , in eq 4 describes that the multiplication between the WF vector and kinetic matrix for a particular mode (r) does not depend on the same operation of the other mode (r′). Such a special feature of the multiplication motivates us to parallelize the TDDVR algorithm, where this operation is performed by a group of processors either declared by programmer or fixed by the default parallel rules. At each time step, the total number of iterations are equally divided among threads so that each of them gets a specified number of iterations. In the recent parallelism, a major part of the TDDVR code is parallelized except for a few portions, that is, time-dependent grid point generation and initialization, etc. This theoretical computation code is parallelized on a multiple core cluster with shared memory using OpenMP-based architecture on the single function with multiple data (spmd) parallelism scheme. In the present version, many unscoped variables (in the potential energy matrix and the quantum propagation part) are declared as private for converting them in a parallel OpenMP-based application. Further, that large number of scratch variables are utilized in the previous parallelism scheme for controlling intermediate results, which are now scoped by private. When we are using this clause, some variables, which need to be shared, are explicitly controlled by the shared clause. The clause dynamic is included in the code as scheduled so that the iterations are divided into chunks with a particular size declared by the programmer. In the former versions, it has been observed that (i) the evaluation of the potential matrix takes a substantial computation time; (ii) the transfer of the total potential matrix, which is actually a very large dimensional array, to each of the computing processors for kinetic matrix− wave vector multiplication and then reverse back to the same matrix takes a considerable elapsed time. Therefore, to reduce such a waste of time for message passing, we have constructed the large dimensional potential matrix in each processor for the potential-vector multiplications. Currently, all of the data exchange among the processors is performed implicitly by accessing shared memory. Such modification in the algorithm further helps to reduce the overall memory of the total dynamical computation. To explore the workability and speedup of this newly implemented parallelized TDDVR algorithm, the two realistic multisurface multimode Hamiltonians of the PNHs are chosen as the testing ground.

quantum picture, the WP in the TDDVR scheme is dictated by classical EOMs:

qṙ c(t ) = pq c (t )

(5)

r

pq̇ c (t ) = − r

dV ({qr }, t ) dqr

+ fr (t )

(6)

The quantity f r(t) denotes the “quantum force” for rth DOF, which is rigorously derived (see eq B1 in the Supporting Information) by minimizing the integral, ∫ [−(iℏ∂Ψk*)/(∂t) − H({pqr},{qr})Ψk({qr},t)*] [(iℏ∂Ψk*({qr},t))/(∂t) − H({pqr}, {qr})Ψk({qr}, t)]∏r p= 1 dqr with respect to ṗqrc, according to the Dirac−Frankel variational principle.42 Although the term “quantum force” used here has no direct relation to that of the Bohmian dynamics,43 it has some resemblance to its actual definition. The Bohmian mechanics is conceptually appealing in the limit where the classical mechanics is well behaved. Such a method can be formulated by introduction of the quantum potential with a scaling factor starting from zero to unity, where the zero indicates the classical limit and the full quantum behavior is switched on for 1. The TDDVR method has the same feature, and it is exact in principle, although it propagates trajectories without using the quantum force as used in Bohmian mechanics. The TDDVR EOMs (eqs 4, 5, and 6) cannot be solved by propagating such trajectories; rather they execute a large dimensional matrix equation, where the size of the matrix equals the number of basis functions involved in the TDDVR expansion (see eq 1). The “quantum force”, involved in the classical EOMs, includes the information on the quantum equations through time-dependent expansion coefficients (am1m2···mr···mp,k(t)), and it is obtained variationally as the procedure of obtaining actual quantum force.

3. PARALLELIZATION SCHEME OF THE TDDVR ALGORITHM The parallelization is a specific tool to handle a multidimensional numerical problem. In parallel computation, such a problem can be divided into smaller ones in its algorithm, and then the different parts can be evaluated concurrently. A parallel architecture composite has several computing nodes, which consist of many processors linked with a single network. All of these threads are connected to a common memory, which can be utilized directly, and the shared data can be accessed by a network from one processor to other. The threads are identical and have equal memory access. For the implementation of our serial TDDVR code, we analyze the performance of the code by profiling, which describes all of the function as well as subroutine calls, the time required for each call, percentage of total time, number of subroutine calls, and many other informations. The nested loop-level parallelism scheme is one of the useful ways of parallelism, where one loop is nested over the other without interfering. We have illustrated in our previous publication41 that the TDDVR approach has enough flexibility to parallelize the substantial areas of its algorithm, and in the subsequent publications,44,45 a part of the code has been parallelized. In the TDDVR algorithm, 94.2% of total computation time is required for the evaluation of the multiplication of the kinetic matrix (T̅ m,m′) and the wave vector (βm1m2···mp) for the rth mode in TDDVR EOM (eq 4), which is the most time-consuming portion of the code. On the other hand, the rest of the time is spent on the βm1m2···mp calculation

4. NEARLY LINEAR SCALABILITY OF TDDVR ALGORITHM The radical cations of the naphthalene and anthracene are included in the present investigation due to their astronomical importance. The proceeding sections demonstrate the accuracy of our parallelized TDDVR methodology to simulate their overlapping PE spectra and DIBs with respect to the results of the well-established molecular dynamical methods, MCTDH33,35 and ML-MCTDH.36,46 Although the computation time exponentially increases with the vibrational modes involved in quantum dynamics, it is now possible to solve the numerical problems with many dimensions due to the considerable increment of the shared as well as distributed memory parallel computers, where many processors can be involved in a single work. A limited number of vibrational modes relevant to dynamics are involved in our treatment considering linear vibronically coupled multisurface multimode Hamiltonians.33 The efficiency and scalibility of both serial and C



Article

tasks. At the time of the parallel execution, the variables used in the code are stored in different storage locations of each worker processor. The utilization of that variable by a thread will refer to that private copy of the variable of the processor. Such memory location of the thread is inaccessible to others. Additionally, more than one iteration of a loop can be performed concurrently in a thread. Thus, multiple threads can invoke the subroutines separately and safely compute without disturbing other processors. Each thread waits awhile for others to complete their iterations to proceed to the next iteration. The recent modifications of the code make it possible to speedup the computations because many costly synchronizations have been eliminated. After completion of execution by all threads, the workers stop computation. Finally, the master thread serially executes to finish the first job and follows the remaining parallel loops. The computation times required to perform dynamical simulation for both PNHs with increasing number of total TDDVR basis sets (ℳ) are listed in Table 2. The speed-up (ℒi) is a well-known quantity used to measure the gain in parallel computation using i number of processors in comparison to a single thread. The speed-up is obtained when the computation time for serial calculation is divided by the time required to calculate the same calculation parallel with i number of threads. Table 2 illustrates the increasing speed-up as a function of processors (i) for inclusion of a greater number of basis sets. The left and right halves of Table 2 present the computation time in real scale for dynamical calculation of Nap+ and Ant+, respectively. In both cases, we have increased the total number of TDDVR basis functions involved in dynamics to solve TDSE with increasing number of processors. It is important to note the speed-up of the parallel TDDVR computation is independent of the number of basis functions involved in dynamics (see the ℒ from top to bottom in Table 2 for a fixed number of processors with increasing TDDVR basis sets). Such a feature is displayed in pictorial form in a separate figure (see Figure 1) for detailed understanding. In this figure, we observe that while the basis set is increased from low (2 × 107) to high (14 × 107), the speed-up of the calculation to solve the coupled differential equation solution remains the same and even slightly increases at high basis set. This is only possible when the major portion of the algorithm is parallelizable and sharing-synchronization overhead is least. Theoretical simulations of these large dimensional systems are computationally expensive as millions of coupled differential equations are solved simultaneously. Figure 2a illustrates the performance of the TDDVR method when the code is implemented using OpenMP-based parallelism scheme. The dashed (black) line in Figure 2a presents the exact linear speedup, whereas the speed-up obtained in the parallel TDDVR algorithm for Nap+ and Ant+ is presented by square (blue) and circle (red) points, respectively. The speed-up points show that an extensive parallelism (∼98% of algorithm) of the TDDVR method reduces the sharing and synchronization overhead such that the amount of ℒi due to parallelization closely equals the number of processors (see ℒi with i = 2, 4, and 6 in Figure 2a) involved in computation. Because distribution of the variables from master to slaves, and thereafter assembly of those dispersed data, requires more time with increasing number of processors, the speed-up gradually decreases with a higher number of processors (see ℒi with i = 8, 10, 12, and 16 in Figure 2a). Scalability analysis is an effective tool for predicting the performance of any algorithm-architecture-based methods,

parallelized TDDVR algorithm are illustrated in the current section. In the multisurface multimode molecular dynamics, a maximum of 13 DOFs are included quantum mechanically, whereas the majority of the vibrational modes are treated classically in both PNHs. Note that the classically treated modes are explicitly involved in the dynamics. The number of basis functions used to represent both quantum and classical EOMs of TDDVR method for Nap+ and Ant+ is presented in the upper and lower parts of Table 1, respectively. A Table 1. Optimum Number of Basis Functions Used for Important Vibrational Modes in TDDVR Calculationa vibrational modes Nap+ η1, η2, η3, η39, η41, η1, η2, η3, η32, η34, η1, η2, η3, η37, η39, η1, η2, η6, η37, η39, η1, η2, η3, η32, η37, Ant+ η1, η2, η3, η58 η1, η2, η3, η47, η50, η1, η2, η3, η39, η40, η1, η2, η7, η51, η55, η1, η2, η3, η47, η55,

TDDVR basis functions

η4, η6, η7, η34, η37, η42 η4, η6, η7, η18, η22, η42 η6, η7, η14, η27, η32, η41, η42 η7, η14, η18, η27, η32, η41, η42 η4, η6, η7, η14, η27, η39, η41, η42

9, 3, 3, 3, 5, 5, 3, 3, 5, 5, 5 7, 3, 3, 3, 3, 5, 5, 7, 5, 5, 5 7, 3, 3, 3, 5, 3, 3, 5, 5, 5, 5, 5 7, 3, 5, 5, 5, 5, 3, 5, 3, 3, 3, 5 7, 5, 3, 3, 3, 3, 4, 3, 3, 3, 5, 5, 3

η7, η9, η47, η50, η51,

9, 5, 3, 5, 5, 3, 5, 5, 5 7, 5, 3, 5, 5, 3, 3, 3, 3, 5, 5, 2 6, 4, 3, 3, 5, 3, 4, 3, 3, 3, 5, 5 7, 3, 5, 5, 3, 5, 3, 5, 3, 3, 3, 5 7, 3, 3, 3, 5, 3, 4, 3, 5, 3, 3, 3, 5

η7, η9, η23, η27, η29, η51, η58 η7, η9, η27, η29, η35, η55, η57, η58 η9, η29, η35, η40, η47, η57, η58 η7, η9, η29, η35, η40, η57, η58, η62

ℳ (×107) 4.10 7.441

initial state X2Au or A2B3u B2B2g

15.94

C2B1g

15.94

D2Ag

13.77

E2B3g

0.75 6.37

X2B2g or A2B1g 2 B Au

15.74

C2B2g

9.56

D2B3u

13.77

E2Ag

In case of Nap+, the 21 normal modes (ν1, ν2, ν3, ν4, ν6, ν7, ν10, ν11, ν14, ν18, ν22, ν27, ν28, ν29, ν32, ν34, ν37, ν39, ν41, ν42, ν46) with six electronic states are included to carry out dynamics. We have performed dynamics on six electronic surfaces entangled with 22 vibrational modes (ν1, ν2, ν3, ν4, ν6, ν7, ν10, ν11, ν14, ν18, ν22, ν27, ν28, ν29, ν32, ν34, ν37, ν39, ν41, ν42, ν46) for Ant+. The normal modes treated quantum mechanically are presented in the table, whereas the DOFs with one basis function (treated classically) are not included in the table. The integer numbers in the second column present the number of basis functions (ηi) used for ith DOF. The quantity ℳ indicates the total number of basis functions used (i.e., the number of coupled differential equations solved at each time step) in the TDDVR dynamics. The last column indicates the state from which the dynamics is initiated. a

SuperServer 8025C-3RB cluster (4 × Intel Xeon E7330 (Quad-core/6 M Cache/2.40 GHz/1066 MHz FSB)) is used for all numerical computations. The parallelized TDDVR algorithm currently solves more than 20 × 107 number of coupled differential equations at each time step for WP propagation. The present implementations depend upon the shared-memory parallelism scheme using OpenMP-based architecture. The work-sharing architecture developed in the parallelism of the TDDVR method is as follows: A particular thread (denoted as master) is utilized to execute the serial part before the initiation of the parallel portion of the TDDVR algorithm. As the program comes to the parallel region, the master thread is instructed to construct a group of processors (called slaves) as workers for simultaneously computing parallel D



Article

Table 2. Computational Advantage of the Parallelized TDDVR Algorithm over Serial One for Large Dimensional Calculationsa Nap+ ℳ i τe ℒi ℳ i τe ℒi ℳ i τe ℒi ℳ i τe ℒi ℳ i τe ℒi

1 120.5

2 61.6 1.95

1 376.1

2 192.2 1.95

1 679.4

2 348.3 1.95

1 1155

2 590 1.95

1 1973.8

2 998.3 1.97

19 136 250 4 8 31.6 16.6 3.80 7.23 57 408 750 4 8 99.4 52.2 3.78 7.20 103 335 750 4 8 178.3 93.8 3.80 7.23 172 226 250 4 8 304.4 158.8 3.79 7.26 287 043 750 4 8 516.1 269.4 3.82 7.32

Ant+ 12 11.6 10.33

16 9.4 12.76

1 91.1

2 46.6 1.95

12 37.2 10.10

16 28.8 13.01

1 487.7

2 248.8 1.95

12 65.5 10.36

16 57.6 13.15

1 700.5

2 357.7 1.95

12 111.6 10.34

16 87.7 13.16

1 1157.7

2 592.7 1.95

12 189.4 10.41

16 150 13.16

1 1518.1

2 776 1.95

13 781 250 4 8 23.8 12.7 3.81 7.13 68 906 250 4 8 128.3 66.6 3.80 7.31 96 468 750 4 8 183.8 96.1 3.80 7.28 151 593 750 4 8 303.8 159.4 3.80 7.26 206 718 750 4 8 398.3 207.7 3.81 7.30

12 8.8 10.25

16 7.2 12.61

12 47.2 10.32

16 37.2 13.10

12 68.3 10.25

16 53.8 13

12 111.6 10.36

16 88.8 13.01

12 145.5 10.43

16 115 13.20

The quantities are demonstrated as follows: ℳ represents the number of TDDVR basis sets used in dynamics, ℒi represents speed-up employing i number of processors in parallelized TDDVR computation, and τe represents elapsed time (in hours) required to propagate the WP.

a

Figure 1. Speed-up (ℒi) of the parallelized TDDVR algorithm with increasing number of basis sets. The hollow (red) and filled (blue) points present ℒi’s for Nap+ and Ant+, respectively. Note that the speed-up of the parallelized TDDVR algorithm is independent of the number of basis functions used in the dynamical calculation. The same speed-up (ℒi) will be obtained for any dimensional computation.

(from top to bottom of Table 2) for fixed i for increment in ηm. Yet the shared memory parallelization of the TDDVR algorithm does not reduce the total computational cost of the calculation, but rather the computational time. Although the computation time of different versions of the MCTDH method, serial MCTDH,47−52 G-MCTDH,53 MLMCTDH,36,46,54−56 etc., is not directly comparable with that of the parallel TDDVR method because the same levels of data are not available, we have presented a few TDDVR calculated outcomes compared to recent ML-MCTDH36 computation results for Nap+ and Ant+ in Table 3. A careful investigation is performed to detect important vibrational modes to simulate

because it measures the efficiency of a parallel algorithm with respect to the number of threads. Switching from serial to parallel version, the TDDVR algorithm shows a closely linear speed-up because an effective sparse matrix-vector multiplication is involved in the formulation (note eq 4). Also, the number of multiplication for η number of basis function for both serial and parallel version is η∑mηm at each step of propagation as shown in Figure 2b, where ηm TDDVR basis functions are used to represent each mth mode. The most important feature of the master-worker parallel TDDVR architecture is that the computation time falls linearly (as seen in Table 2 from left to right), and the ℒi tends to increase E



Article

involvement of all vibrational DOFs improves a little in the dynamical results. By further illustration of the computational aspects of both dynamical methods, we have depicted the speed-up data in Figure 2a for the parallel TDDVR method in case of Nap+ and Ant+ with 30 × 107 TDDVR basis functions in comparison with the available parallel MCTDH speed-up data. In Figure 2a, the ℒi’s of shared memory (SM-MCTDH) and distributed memory (DM-MCTDH) parallelism of MCTDH algorithm55 are presented with (11,12,12,11,7,7) SPFs (i.e., Avector length of 853 776 coefficients). The MCTDH results are reported for a full 15-dimensional quantum-dynamical simulation of the protonated water-dimer. In this figure, both the SM-MCTDH and the DM-MCTDH show lower speed-up than the TDDVR parallelization. It can be concluded from (i) Figure 2 that speed-up of parallel TDDVR computation is greater than that of parallel (both SM and DM) MCTDH method;55 and (ii) Table 3 that the computation time required by the parallel TDDVR method is much less than that by the MCTDH method, and slightly greater than the ML-MCTDH method to reproduce the same level of dynamical results. Additionally, the speed-up of the parallel MCTDH method gradually falls with the number of coefficients, whereas that of the parallel TDDVR method does not change with an increase in the number of basis functions.

5. DYNAMICAL RESULTS AND DISCUSSIONS As several complex dynamical phenomena are involved in the electronic manifold of both PNHs, a six-state 48 mode Hamiltonian33 for Nap+ and 66 mode Hamiltonian33 for Ant+ are adopted to perform quantum-classical dynamics, where the first- and second-order coupling terms are used to represent the tuning10 as well as coupling modes. The normal-mode frequency values, their linear and quadratic coupling constants, and ionization potential values are utilized from ref 33. Many experimental groups23,27−31 have recorded the photo absorption spectra of naphthalene as well as anthracene, and revealed that the PE spectra originated due to the six lowest conical intersecting electronic states (D0−D5). The PESs are vibronically coupled through several vibrational modes differently, and therefore we have carefully found out the most important modes for spectral simulation in both cases. As an example, the ν37, ν39, ν41, ν42 play a crucial dynamical role in Nap+, because they couple between X̃ -Ã , B̃ -C̃ (not ν42), and D̃ -Ẽ states, whereas the B̃ -D̃ , C̃ -Ẽ PESs are coupled through ν27, ν28, etc. Similar tasks are also performed in case of Ant+, where the ν1, ν2, ν7, ν9, ν47, ν57, ν62 modes are most relevant to the dynamics. Because we have adopted the Nap+ and Ant+ as two testing examples for application of the newly implemented parallel TDDVR method, it is interesting to observe the TDDVR calculated dynamical results, and their analogy with the outcomes34−36 of the well-established MCTDH method.47−51 The authors of ref 33 have revealed that 29 out of 48 vibrational modes (for Nap+) are relevant to the dynamics, which are distributed in the irreducible representation as 7ag + 3au + 2b1g + 4b1u + 4b2g + 3b2u + 4b3g + 2b3u. We have determined that 21 out of those 29 DOFs are enough to reproduce the photoelectron band of Nap+, in which many DOFs can be controlled classically. On the other hand, it was previously reported33 that 31 out of 66 modes have a major contribution to the nuclear dynamics for Ant+, which are classified as 9ag + 1au + 5b3g + 2b3u + 3b2g + 5b2u + 1b1g + 5b1u. Almost 13 vibrational modes in case of Ant+ are observed as quantummechanically important to perform dynamics. We have

Figure 2. (a) The parallel TDDVR algorithm shows nearly linear speed-up (points connected by dotted lines) with increasing number of processors, whereas the dashed line presents the exactly linear speed-up. Almost 30 × 107 number of TDDVR basis functions are used to check the speed-up of the parallel TDDVR algorithm. A comparative speed-up of the SM and DM parallelization of MCTDH algorithm55 is presented with (11,12,12,11,7,7) SPFs with A-vector length of 853 776 coefficients. MCTDH results are reported for a full 15-dimensional quantum-dynamical simulation of the protonated water-dimer. (b) The number of multiplications is performed by TDDVR algorithm at each time step for a given number of basis sets.

the PE spectra of both PNHs. Thus, we have included quantum mechanically most important 9−14 DOFs in our calculation depending on the electronic states of the molecular systems. The computation time displayed in the left part of Table 3 is obtained from parallel TDDVR calculation (using 16 processors) with optimized basis functions, where the WP is propagated up to 400 fs. The right part of the same table illustrates the full dimensional ML-MCTDH36 calculation for Nap+ and Ant+ (takes 60−140 and 40−150 h, respectively) requires less elapsed time than the reduced dimensional MCTDH computation35 (takes 160−350 and 160−230 h, respectively). It is interesting to note that the CPU time required for parallel TDDVR calculation is 2−3 times less than the reduced dimensional MCTDH calculation35 to reproduce the same dynamical results. Our theoretical results successfully reproduce (which will be discussed in the next section) the outcomes of reduced dimensional MCTDH calculation35 and the features of 48D ML-MCTDH36 results within this short time scale. It is important to mention that TDDVR calculated results are not comparable with those of the full dimensional MCTDH computation36 because a higher number of DOFs are involved (more important DOFs are involved) there. Yet this F



Article

Table 3. Comparison of Computation Time for the MCTDH,35 ML-MCTDH,36,46 and TDDVR Algorithmsa initial state

method

CPU time (h:min)

TDDVR basis set (×107)

method

CPU time (h:min)

MCTDH coefficients

Nap+ X2Au

TDDVR

42:13

4.100

A2B3u

TDDVR

41:06

4.100

B2B2g

TDDVR

76:40

7.441

C2B1g

TDDVR

170:00

15.94

D2Ag

TDDVR

171:06

15.94

E2B3g

TDDVR

132:13

13.77

48D 29D 48D 29D 48D 29D 48D 29D 48D 29D 48D 29D

ML-MCTDH MCTDH ML-MCTDH MCTDH ML-MCTDH MCTDH ML-MCTDH MCTDH ML-MCTDH MCTDH ML-MCTDH MCTDH

92:06 283:47 96:09 351:41 59:35 208:15 120:34 338:35 143:09 245:41 116:17 163:57

59 459 252 540 59 459 252 540 48 004 252 540 81 799 252 540 83 613 252 540 44 525 252 540

66D 31D 66D 31D 66D 31D 66D 31D 66D 31D 66D 31D

ML-MCTDH MCTDH ML-MCTDH MCTDH ML-MCTDH MCTDH ML-MCTDH MCTDH ML-MCTDH MCTDH ML-MCTDH MCTDH

39:00 154:25 39:49 216:25 65:52 208:25 139:33 160:04 136:17 187:07 153:49 228:57

64 020 2 395 224 64 020 2 395 224 115 112 2 395 224 98 859 2 395 224 92 031 2 395 224 104 462 2 395 224

Ant+ 2

X B2g

TDDVR

7:46

0.75

A2B1g

TDDVR

8:00

0.75

B2Au

TDDVR

68:53

6.37

C2B2g

TDDVR

169:19

15.74

D2B3u

TDDVR

97:46

9.56

E2Ag

TDDVR

131:26

13.77

a

The WP is propagated up to 400 fs in all cases. All of the data of MCTDH simulation were taken from ref 36, where the calculations are run on the machine and CPU type (AMD OpteronTM Processor 246, 2000 MHz). However, all of the TDDVR calculations are performed on a SuperServer 8025C-3RB cluster using 16 processors. The fourth column presents the total number of optimized TDDVR basis functions used in the present calculation, whereas the seventh column shows the total number of time-dependent MCTDH coefficients propagated in each case.36

Figure 3. Convergence study of population dynamics with increasing number of TDDVR basis functions. For Nap+ the WP is initiated from (a) Ẽ 2B3g, (b) C̃ 2B1g, and (c) Ã 2B3u PESs, and the corresponding population profiles are presented. The Ã 2B3u state population profiles are displayed in panel (d) for Ant+, when the dynamics is performed locating the WF from the B̃ 2Au state.

G



Article

uppermost one. Both the ν6 and the ν7 modes play a distributive role in the spectral envelop, and therefore are essential in nonadiabatic dynamics for all electronic states. Other tuning as well as coupling modes, using single basis function, are treated classically in our investigation due to their negligible contribution. The five population profiles of Figure 3a have been calculated using five different sets of TDDVR basis functions. The dotted (black) line is obtained with the lowest TDDVR basis set (2.01 × 107), where the distribution of the primitive basis functions is presented in set Ia of Table 4. The same calculation is repeated

included 21 modes for Nap+ and 22 modes in case of Ant+ in dynamics (see Table 1), where many DOFs are tracked classically. The TDDVR methodology has the specialty to distinguish the quantum and classical nature of the vibrational modes by calculating the corresponding time-dependent coefficients of the WF, and thereby it treats them accordingly. Such nature obviously depends on the electronic state from which the dynamics is initiated because the electronic surfaces are entangled differently by the normal modes. In case of Nap+, this group of modes (ν10, ν11, ν14, ν18, ν22, ν27, ν28, ν29, ν32, ν46) is treated classically when the WF is initiated from the lowest electronic state, whereas another group of modes (ν10, ν11, ν18, ν22, ν28, ν29, ν34, ν46) is controlled classically for the Ẽ 2B3g surface (see Table 1). Similar phenomena are also observed for the other electronic surfaces of Nap+ and for Ant+. The dynamical observables are reported in terms of different types of spectra and nonradiative population transfer from upper electronic states to lower ones. The subsequent section illustrates the convergence of TDDVR calculated dynamical results to establish the accuracy of TDDVR algorithm with respect to the MCTDH methodology. We have also monitored the nuclear DOFs for reduced densities along with population dynamics to elucidate the internal features of population dynamics. For better understanding of different nonadiabatic phenomena, we have presented the dynamical results of both of the systems in a comparative manner in the proceeding sections. 5.1. Analysis of Convergence. The convergence study of the basis functions is an important criteria for any quantum dynamical method. Qualitative calculations could be possible with a small number of basis functions, but to obtain numerically exact results, it is necessary to use an optimum number of basis functions, otherwise the result would be erroneous. Therefore, the convergence of all of the dynamical findings from both population dynamics as well as photoelectron spectra is studied and presented elaborately in the current section. Because the observables originated due to different types of vibronic coupling of various DOFs entangled in PESs, we have treated different vibrational modes involved in the model Hamiltonians with varying TDDVR basis functions. The number of basis functions increased from lower to higher until the property of our interest converges. Although a large endeavor with several basis sets is carried out for both of the systems, we have presented only a few significant results because the demonstration of all of those data is beyond the scope of this Article. 5.1.1. Population Dynamics. a. Naphthalene. The present nuclear dynamics is carried out initiating the WP from the FC point of the Ẽ 2B3g, C̃ 2B1g, and Ã 2B3u electronic surfaces of sixsurface 21 mode Hamiltonian of Nap+, and the corresponding dynamical results are displayed in panels a, b, and c of Figure 3, respectively. The convergence of basis functions used in the TDDVR method is tested separately for each state in the population dynamics. For clear visualization, a particular region of all population profiles is enlarged in the inset of each panel. We have first determined the priority of all vibrational modes, then included 21 DOFs in molecular dynamics. After that, we carefully distributed the modes for quantum and classical treatment according to their contribution to the dynamics. The ν1 is always taken as a strongly active vibrational mode, whereas the ν2 is considered as a moderately significant one. The ν3 is reasonably important for C̃ and Ẽ electronic states, whereas ν4 is moderately active for the two lowest states (X̃ and Ã ) and the

Table 4. Number of Basis Functions Used for Different Vibrational Modes in TDDVR Dynamics for Nap+, Where Other DOFs Are Treated Classically with One TDDVR Basis Functiona vibrational modes ν1, ν2, ν6, ν7, ν11, ν14, ν18, ν25, ν27, ν32, ν37, ν39, ν41, ν42 ν1, ν2, ν6, ν7, ν14, ν18, ν25, ν27, ν32, ν37, ν39, ν41, ν42 ν1, ν2, ν3, ν4, ν6, ν7, ν14, ν27, ν32, ν37, ν39, ν41, ν42 ν1, ν2, ν3, ν4, ν6, ν7, ν14, ν27, ν32, ν37, ν39, ν41, ν42 ν1, ν2, ν3, ν6, ν7, ν14, ν27, ν32, ν37, ν39, ν41, ν42 ν1, ν2, ν3, ν4, ν6, ν7, ν14, ν25, ν27, ν32, ν39, ν41, ν42 ν1, ν2, ν6, ν7, ν14, ν18, ν27, ν32, ν37, ν39, ν41, ν42 ν1, ν2, ν6, ν7, ν14, ν27, ν28, ν32, ν37, ν39, ν41, ν42 ν1, ν2, ν6, ν7, ν14, ν18, ν27, ν32, ν37, ν39, ν41, ν42 ν1, ν2, ν6, ν7, ν14, ν18, ν22, ν25, ν27, ν37, ν39, ν41, ν46

distribution of basis functions

total basis functions

name of basis set

4, 3, 2, 2, 3, 3, 4, 2, 3, 3, 3, 3, 3, 4

2.01 × 107

7, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 3, 3 8, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 3, 3 7, 5, 3, 3, 3, 3, 4, 3, 3, 3, 5, 5, 3 8, 6, 4, 4, 3, 3, 3, 4, 4, 5, 5, 4 5, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3 7, 5, 3, 5, 3, 3, 3, 3, 3, 5, 5, 5 7, 3, 3, 5, 3, 3, 3, 5, 5, 5, 5, 5 8, 4, 4, 5, 5, 5, 3, 5, 4, 3, 3, 5 6, 3, 3, 5, 3, 3, 4, 4, 5, 3, 3, 3, 3

∼4.00 × 107

set IIa

∼11.0 × 107

set IIIa

∼14.00 × 107

set IVa

∼20.00 × 107

set Va

∼1.60 × 107

set Ib

∼9.50 × 107

set IIb

∼16.0 × 107

set IIIb

∼25.0 × 107

set IVb

∼9.44 × 107

set Vb

set Ia

a

The total number of basis functions used in dynamics are shown in the third column.

with a total of ∼4.00 × 107 TDDVR basis functions (see set IIa of same table), and the corresponding population profile [dashed (red) line of Figure 3a] makes a wide difference with the former one. The dot dashed (green) line of the same panel is obtained by solving more than 10.70 × 107 coupled differential equations with the TDDVR basis functions of set IIIa in Table 4. Indeed, we have increased the basis functions to more than 2 times that of 4 × 107, but it does not affect substantially the dynamics, because this increment is on the tuning modes, which are almost converged. The solid (blue) line is achieved using the total TDDVR basis set of ∼14.00 × 107 (for distribution of basis functions, see set IVa of Table 4). Although the Ẽ state population profile with lowest basis set (set Ia) decays more rapidly than any other profile, it is not the converged one because further increment of basis makes slower the ultrafast decay and ultimately converges. The solid (blue) line does not change with the increment of basis functions either on tuning vibrational mode or coupling one. Therefore, it is concluded that we have reached the optimized basis set for the present calculation. The increment of TDDVR basis functions to any vibrational DOF beyond optimized basis set does not influence the decaying process. As a case study, a large TDDVR basis set ∼20.00 × 107 (note the set Va) is attempted, H



Article

The dot dot dashed (magenta) line is obtained by solving almost 18.00 × 107 coupled differential equations [set Vc of Table 5 leads to 172 800 000 TDDVR basis functions], but the calculated profile closely follows the former one. This indicates that the solid (blue) line is well converged by the utilized TDDVR basis set. b. Anthracene. Because the detailed convergence study of the individual state population profiles of two PNHs is beyond the scope of this Article (as the study of the individual profile is space consuming), we have displayed only one panel (Figure 3d) to show convergence of population dynamics of Ant+. The WP (WP) dynamics is performed initiating the function from B̃ 2Au state with different TDDVR basis sets, and analogous population profiles of A2B1g state are presented in Figure 3d. The three population profiles of Figure 3d have been calculated with the following basis sets: (1) The set Id of Table 5 leads to ∼1.24 × 107 TDDVR basis functions; (2) the optimized TDDVR basis set (set IId) amounts to 6.37 × 107 coupled differential equations; abd (3) by increasing the basis set more abruptly [set IIId of Table 5 leads to ∼65 × 107 TDDVR basis functions] than the optimized set, the same calculation is repeated, but the profile remains almost unaltered. This indicates that the TDDVR basis functions used in the dynamics are truly converged. 5.1.2. Photoelectron Spectra. a. Naphthalene. The photoelectron spectra, containing six overlapping energy bands, of Nap+ are simulated with the help of the parallelized TDDVR method. Three low-lying electronic states (X̃ , Ã , and B̃ ) contribute to three lower energetic bands, and the next extremely overlapping higher energetic band is produced due to C̃ , D̃ , and Ẽ states. The ν6 and ν7 vibrational modes have a major contribution to the three peaks of X̃ state located at ∼8.177, ∼ 8.37, and ∼8.55 eV. The improvement of the broad structure of the Ã state (positioned at ∼8.93−9.27 eV) is (see Figure 4a) acquired with the major input of ν2 and ν7. The B̃ band of Nap+ has been the subject of major interest because, in addition to the presence of the diffuse interstellar bands, it consists of three new interstellar or circumstellar bands30 as measured by gas-phase laboratory spectroscopy at low temperatures. Such identification of new interstellar features supports that these simplest polycyclic aromatic hydrocarbons are the carriers of both DIBs and anomalous microwave emission. To capture the DIBs of the B̃ state, ν1 along with ν6 and ν7 play a crucial role. These characteristic features of each vibrational mode are reflected in the choice of optimized basis functions, where a considerable number of TDDVR basis functions are utilized for ν6, ν7, and ν1 when the dynamics is initiated from X̃ , Ã , and B̃ states (see Table 1). On the other hand, ν2 and ν6 have moderate contribution to the PE band of C̃ state. Unlike any other state band, ν3 is essential for Ẽ state spectra with ν1. The Ẽ state band, spreading over a broad region, can be acquired with the help of ν1, ν3, and ν7. It is a difficult task to separate out the individual contributions of C̃ , D̃ , and Ẽ states, because they construct a highly overlapping C̃ D̃ -Ẽ band structure. The nuclear dynamics is carried out on the model Hamiltonian33 of Nap+ using a different number of TDDVR basis sets locating the initial WF on each of the six electronic states. The resultant autocorrelation functions (C(t) = ⟨Ψ(t)|Ψ(0)⟩) are then utilized to calculate the individual state PE spectra using the low damping parameter (τ = 50 fs). These separate envelopes are convoluted with equal weight factor to capture the total PE spectra of Nap+. Those PE spectra are displayed in Figure 4a with three distinct line types.

and the computed population profile is displayed with a dot dot dashed (magenta) line. The last two profiles closely follow each other, indicating the solid (blue) line is the fully converged profile. In Figure 3b, the five population profiles have been calculated with the following basis sets. The dotted (black) line is obtained using a minimum basis set 1.5 × 107 functions (see set Ib of Table 4). The dashed (red) line is computed with a increment of the TDDVR basis set to almost ∼9.50 × 107 (see set IIb). The converged population profile (see dot dashed (green) line) is achieved with ∼16 × 107 TDDVR basis functions as presented in set IIIb. We have enormously increased the TDDVR basis functions to more than 25.00 × 107 according to set IVb of Table 4 to calculate the population profile as indicated by the solid (blue) line. It is gratifying to see the population profiles gradually converge with a systematic increase of the total TDDVR basis set. Interestingly, if we increase the basis functions to double the optimized one, the converged profile remains the same, indicating accurate quantum dynamics. On the contrary, the dot dot dash (magenta) line is obtained using a substantial basis set Vb (in Table 4), but because of the improper distribution of the TDDVR basis functions, the profile remains away from the converged one. Similarly, Figure 3c presents the characteristic population profiles of Nap+ when the dynamics is initiated from Ã 2B3u electronic state employing different TDDVR basis sets. The lowest number of TDDVR basis set (see set Ic of Table 5) is Table 5. Same as Table 4 with set Ic−set Vc for Nap+ and set Id−set IIId for Ant+ vibrational modes

distribution of basis functions

total basis functions

name of basis set

ν1, ν2, ν6, ν7, ν34, ν37, ν39, ν41, ν42 ν1, ν2, ν3, ν4, ν6, ν7, ν34, ν37, ν39, ν41, ν42 ν1, ν2, ν3, ν4, ν6, ν7, ν34, ν37, ν39, ν41, ν42 ν1, ν2, ν3, ν4, ν6, ν7, ν34, ν37, ν39, ν41, ν42 ν1, ν2, ν3, ν4, ν6, ν7, ν34, ν37, ν39, ν41, ν42 ν1, ν2, ν3, ν7, ν9, ν23, ν27, ν29, ν47, ν50, ν51, ν58 ν1, ν2, ν3, ν7, ν9, ν23, ν27, ν29, ν47, ν50, ν51, ν58 ν1, ν2, ν3, ν7, ν9, ν23, ν27, ν29, ν47, ν50, ν51, ν58

6, 4, 3, 3, 3, 3, 3, 3, 3 5, 3, 3, 3, 3, 3, 3, 3, 5, 5, 5 6, 4, 4, 4, 4, 4, 3, 3, 4, 4, 4 9, 3, 3, 3, 5, 5, 3, 3, 5, 5, 5 9, 4, 4, 4, 5, 5, 4, 4, 5, 5, 5 7, 5, 2, 5, 5, 3, 4, 4, 3, 2, 2, 2 7, 5, 3, 5, 5, 3, 3, 3, 3, 5, 5, 2 8, 6, 4, 6, 6, 4, 4, 3, 3, 6, 6, 3

∼0.31 × 107

set Ic

∼0.82 × 107

set IIc

∼2.14 × 107

set IIIc

∼4.10 × 107

set IVc

∼18.0 × 107

set Vc

∼1.24 × 107

set Id

∼6.37 × 107

set IId

∼65.0 × 107

set IIId

used to calculate the profile as displayed by the dotted (black) line. Next, two higher basis sets are attempted: first ∼0.82 × 107 (see set IIc) and second ∼ 2.14 × 107 (note set IIIc), where the corresponding profiles are indicated by dashed (red) and dot dashed (green) lines, respectively. The solid (blue) line is obtained with more than 4.00 × 107 TDDVR basis functions (for the distribution of basis functions, see set IVc of Table 5). The noticeable feature of the Ã state population profiles is that the decaying rate of Ã state population is boosted with an increase in the number of basis functions and ultimately approaches convergence. It observed that the solid (blue) line remains unaltered with the increment of basis functions. We increase the functions more than 4 times than the 4.00 × 107, but such a hike does not significantly affect the ultrafast decay. I



Article

Figure 4. Convergence study of the PE spectra for (a) Nap+ and (b) Ant+ with increasing TDDVR basis sets. The theoretical PE spectra are calculated with τ = 50 fs and τ = 45 fs, respectively. The PE spectra obtained with more than optimized basis sets are almost identical to the optimized one, and therefore are not displayed in this figure to avoid crowding.

bands originated due to the ground (X̃ ), first excited (Ã ), and second excited (B̃ ) electronic states, and a higher energetic band is developed due to C̃ , D̃ , and Ẽ states. Although it was previously mentioned in ref 35 that nine tuning vibrational modes are relevant to dynamics of Ant+, we have analyzed that the five DOFs are sufficient to reproduce the PE spectra, and out of them all do not equally contribute to the band. The X̃ state band is a composite of three main peaks located at ∼7.17, ∼7.35, and ∼7.53 eV with diminishing intensity (see Figure 4b), which can be achieved with major input of ν6 or ν7 and ν9 modes. Further, the ν3 and ν9 are needed to catch the three closed positioned peaks at ∼8.44, ∼ 8.54, and ∼8.63 eV of the Ã state. The B̃ band has drawn the large attention of astronomers because of its analogous new broad interstellar band of the stars Cernis 52,29 a likely member of the very young star cluster IC 348 and HD281159. Details on this issue have been discussed in section 5.2.3. The weak coupling of B̃ state with the lower electronic state (Ã ) functions in a major role to construct such DIBs. The B̃ state band consists of only one peak with an additional collapsed structure, which can be obtained by the ν7, ν9 with a limited contribution of ν2. The significant progressions of these vibrational modes are reflected in the choice of basis set as presented in Table 1. The overlapping structure of the C̃ , D̃ , and Ẽ states of Ant+ is

Theoretical PE spectra, consisting of the contribution of the six electronic states, of Figure 4a are calculated with three different numbers of TDDVR basis sets. The dot dash (green) line is obtained with 107 number of TDDVR basis functions. We have increased the basis set more than six times (∼6.05 × 107), and thereby the calculated PE spectra [dashed (blue) line] substantially differ from the former one at high energy region but not at low. The solid (red) line is obtained with the optimized TDDVR basis set (>15.9 × 107). Because there is no significant difference between the spectra displayed by dashed (blue) line and solid (red) line, we have increased the TDDVR basis functions (∼34.00 × 107) to compute the same band, but it gives almost the same result (not displayed in the figure to avoid overlapping) as the optimized one. Therefore, similar to the converged population profiles, further increment in the TDDVR basis functions beyond the optimized basis set does not show any significant change in the converged PE spectra as shown in Figure 4a. Note that the TDDVR calculated spectral profiles readily converge at low energy region with a lower number of basis functions, but a greater number of TDDVR basis functions are required to converge at higher energy region. b. Anthracene. Similar to Nap+, we have quantum-classically simulated the PE spectra of Ant+, where three lower energetic J



Article

Figure 5. Diabatic population profiles resulting from six-surface 21 mode dynamics for Nap+ initiating the WP from (a) E2B3g, (b) D2Ag, (c) C2B1g, (d) B2B2g, (e) A2B3u, and (f) X2Au surfaces, respectively. In all cases, the optimized basis sets are used to propagate WP.

expanded over an energy range ∼10.18−10.90 eV, but this nonuniform hump is not so diffused as Nap+. The ν1, ν7, and ν9 make significant contributions to C̃ state band consisting of many peaks over the region ∼10.38−10.84 eV, when the individual spectra are calculated with τ = 45 fs. The D̃ state spectra consist mainly of a peak at ∼10.64 eV with additional low intensity peaks at lower and higher energy range, where the ν7 and ν8 modes have the dominant progression. The Ẽ state band, containing many structures in the region ∼10.31−10.89 eV, can be obtained with a significant contribution of modes ν1, ν7, and ν9. The convergence study of the TDDVR basis sets for the simulation of PE spectra of Ant+ is also performed similar to Nap+. In Figure 4b, (i) the lowest basis set (∼ 1.50 × 107) is used to calculate the dot dash (green) line; (ii) the dashed (blue) line is the result of dynamics with a medium basis set (∼6.00 × 107); and (iii) the solid (red) line is obtained with the optimized basis set (∼ 15.7 × 107). It is clearly detected in those spectral profiles that with increasing basis sets, the broadening and peaks gradually converge, but the number of peaks remains the same with identical positions. The optimized spectral profiles for Ã and B̃ states are of slightly greater intensity comparative to the other state profiles. The increment of the TDDVR basis functions beyond optimized basis set does not show any significant difference of the converged PE spectra as usual. Table 1 presents the optimized number of basis

functions required in TDDVR dynamics, which is substantially lower than the basis set used by Ghanta et al.,35 but the number of basis functions distributed on various modes in TDDVR dynamics is similar to the dynamics of the MCTDH method. 5.2. Workability of the Parallelized TDDVR Method. Wave packet dynamics on a multisurface multimode Hamiltonian is challenging for any quantum dynamical method because several CIs prevailed among the PESs that demand an accurate computational framework. With this aim, we have performed dynamics on two such molecular systems, Nap+ and Ant+, to explore their internal dynamical features. These polycyclic aromatic hydrocarbons are considered as the natural carriers of DIBs in the interstellar medium, and therefore have special interest in astronomical research. The femtosecond internal conversion dynamics of a WP initially prepared on different PESs for both systems is presented here. We are presenting our TDDVR calculated time-dependent observeables in the current section in a comparative manner with the results computed by the well-established quantum dynamical method, MCTDH.47−51 5.2.1. Ultrafast Population Decay Dynamics in Coupled Electronic Manifold. As the analysis of population is an effective tool to understand the femtosecond internal conversion dynamics of any molecular system, we have presented the population dynamics in detail for both PNHs. K



Article

Figure 6. Contour diagram of the WP on two important DOFs (q1 and q2) for Nap+. The snapshots are taken at different time intervals when the WP is initiated from (i) Ẽ 2B3g state [panels (a), (b), and (c) are contours on D̃ surface], [panels (d) and (f) are on C̃ PES], [panel (e) is on B̃ surface]; (ii) D̃ 2Ag state [panels (g), (h), and (i) are contour images of WP on D̃ surface], [(j) and (k) are on C̃ potential surface]; (iii) C̃ 2B1g state [(l) for state D̃ ], [panels (m), (n), and (o) display WP on B̃ surface].

strongly coupled and their CIs exist within ∼0.58 eV, a major portion (45%) of the Ẽ state population sharply decays to D̃ state within ∼25 fs, and that of ∼35% population quite slowly transfers to the C̃ state. Only 5−10% of Ẽ state population transfers to the B̃ state. The intrinsic features of this population dynamics can be understood with the help of the reduced density plots of the WP on two important coordinates. However, the actual dynamics is much more complicated; it would be understandable from the following discussion. Let us discuss the WP on two important vibrational DOFs (q1, q2) planes in each PES. At the start of the dynamics, the WP exists at (q1 = 0, q2 = 0) on Ẽ surface. The solid lines in Figure 6a−c indicate the D̃ -Ẽ seam. Because the (0,0) point exists in the vicinity of the D̃ -Ẽ seam, the packet immediately reaches the CI, after starting the dynamics, to decay to the lower D̃ PES. Therefore, the D̃ state becomes substantially populated within t

To explore the dynamical features originating from the ultrafast population transfer, we have separately carried out quantumclassical dynamics initiating the WP from all of the PESs. The investigation is further enriched with the study of the diabatic reduced densities, ζk(qm, qn, t) = ∫ Ψ*k ({qm}, t)Ψk({qn},t)∏l≠m,n dql of the time-dependent WF on different PESs l = D̃ 0, D̃ 1, ..., D̃ 5, where the kth component of the total WF Ψ({q},t) is presented by Ψk({q}, t). The Ψl({q},t) consists of the contributions of six low-lying electronic surfaces of both PNHs. a. Napthalene. We are discussing the diabatic population dynamics of Nap+ when the initial WP is propagated individually from the Ẽ 2B3g to X̃ 2Au potential surfaces to carry out dynamics. Figure 5a−f displays the time-dependent population obtained from that dynamics. Panel a presents the electronic population of six electronic states when the initial WP is prepared for Ẽ state. Because the C̃ -D̃ -Ẽ states are L



Article

Figure 7. Diabatic population dynamics of Ant+ using six-state 22 mode Hamiltonian initiating the WP from (a) Ẽ 2Ag, (b) D̃ 2B3u, (c) C̃ 2B2g, (d) B̃ 2Au, and (e) Ã 2B1g and X̃ 2B2g electronic surfaces, respectively. Because of very weak coupling between X̃ state with other states, the rate of the decaying process is very slow as seen in panel (f).

= 3 fs (see the D̃ state profile in Figure 5a). On the D̃ surface, the function shifts toward positive q2 direction (see Figure 6a), and then (t = 50 fs) translates back toward negative q2 (see Figure 6b). At longer time (t = 110 fs), the WP moves again toward positive q2 direction (see Figure 6c). Indeed, the WP oscillates about q2 coordinate, which is reflected on the D̃ state population profile that shows an upward shift at t ≈ 3 fs, downward movement at t ≈ 50 fs, and again upward shift at t ≈ 110 fs. We observe that the WP on the Ẽ electronic surface is gradually decreasing and remains almost stationary on its position. As the D̃ -Ẽ and C̃ -D̃ seams are energetically close, the WF reaches to the next lower C̃ state in a similar fashion. Reaching the C̃ surface, the WP shifts toward negative q2 axis and, interestingly, splits into two parts (see Figure 6d): one reaches to the lower B̃ surface (see Figure 6e) and the other gradually shifts toward positive q2 axis on the same PES (see Figure 6f). In summary, the gradual shift of the WP toward the positive q2 axis is reflected on the continuous increasing population profile (with a saturated value ∼35%) of C̃ state in Figure 5a. Diabatic population dynamics, initiating the WP from D̃ state, of six electronic manifold are shown in Figure 5b. Beneath the D̃ PES, the C̃ state and B̃ electronic states are rapidly populated after initiation of dynamics. Because the C̃ -D̃ CI is lying less than 0.06 eV from D̃ state and the B̃ -D̃ states are

strongly coupled33 through ν27 and ν28 modes, an ultrafast population transfer occurs from D̃ to lower B̃ and C̃ states. A considerable amount of population (more than 5%) is transferred to the Ẽ state, which indicates that the WP can suitably access the upper lying D̃ -Ẽ CI. The gradually increasing population of Ã state from 0% (at ∼20 fs) to almost 5% (at 150 fs) is due to the vibronic coupling between this state with B̃ and C̃ states because Ã state cannot access the D̃ -Ẽ CI. After a sharp fall, the D̃ state population profile is slowly decaying with a nonuniform hump at ∼100 fs. At the beginning, the WP is located at (q1 = 0, q2 = 0) on the D̃ surface. Up to t ≈ 60 fs, the WP shifts toward the positive q2 direction [see Figure 6g]. It then backs [see Figure 6h] to almost the starting position to further shift beyond [see Figure 6i] the D̃ -C̃ seam. After t = 100 fs, because it moves far from the D̃ -C̃ seam, the decay rate becomes slower. The characteristic features of the WP on D̃ surface are continuous shifting toward positive q2, then a back shift and further a permanent shift toward positive q2 is reflected in the population profile of the D̃ state, where the population continuously decays up to 60 fs (from vicinity of D̃ C̃ seam to go far), then it slightly increases (back toward D̃ -C̃ seam) and remains almost uniform with slight decaying (far from D̃ -C̃ seam). When the WP at D̃ state remains near the D̃ C̃ seam, it rapidly decays to lower C̃ state without movement. The WP on the C̃ surface shows a feature similar to that of the M



Article

Figure 8. Contour plots of the diabatic WF for Ant+ initiating the WP from (i) D̃ PES; (ii) B̃ state [panels (f), (g), and (h) are pictures on B̃ surface], [(i) and (j) for state Ã ]. Panels (a) and (b) present the snapshot of the WF on D̃ surface, whereas panels (c), (d), and (e) present for state C̃ in former case (i).

upper state. After gaining a substantial population on the C̃ state, the WP slowly moves toward positive q2 direction (see Figure 6j), then backs (see Figure 6k) to remain static (where it decays to lower states) at almost equilibrium point. When we have initiated dynamics from C̃ state, the ultrafast relaxation occurs mainly (more than 70%) to the upper lying D̃ surface rather than lower energetic B̃ state (maximum 20%), which is unexpected because access of the WP to C̃ -D̃ CI is more difficult than low-lying C̃ -B̃ CI. Further, the slower decay rate of the WP on the C̃ state in the absence of D̃ state indicates the major role of D̃ state in case of C̃ state non adiabatic decay dynamics. Almost 5% population growth of the Ẽ state at the beginning of the dynamics (∼2 fs) indicates a high impact of the upper lying electronic states to the internal conversion process of lower PESs. Because the C̃ and D̃ states are very close in energy, at the start of dynamics from C̃ state

the population rapidly transfers to the D̃ state, where the WP remains in the vicinity of the C̃ -D̃ seam [see Figure 6l] throughout the dynamics. On the contrary, the WP on B̃ state moves to the positive q2 axis [see Figure 6m] showing a node, then comes to an equilibrium position [see Figure 6n] to shift again to the positive q2 axis [see Figure 6o]. Such oscillations of the WP hamper the decaying process on the B̃ state. Figure 5d presents the six state population profiles, while the dynamics is performed initiating the WP from B̃ state. The slow decaying process of the B̃ state to the lower Ã and X̃ states indicates the presence of a weak vibronic interaction among the states through B̃ -C̃ and B̃ -Ã CIs. Such emission helps to produce the diffuse interstellar bands of B̃ state, which makes naphthalene an astronomically important one. This phenomenon supports the fluorescence emission of this state. A little population transfer to the other excited electronic states, C̃ , D̃ , N



Article

and Ẽ states, inspired us to treat the corresponding coupling vibrational modes (ν14, ν27, ν28, ν37, ν39, and ν41) classically (see Table 1). The time-dependent six state population dynamics locating the WP from Ã state of Nap+ is displayed in Figure 5e. In the panel, only the Ã and X̃ state profiles are visible, but other state populations are too low to observe. Actually, the present system can be divided into two subspaces: the X̃ -Ã states and other B̃ C̃ -D̃ -Ẽ states. The upper state CIs are lying too high to access the WP of Ã or X̃ , and thus effectively no population transfer occurs from the lower to upper subspace. Because the X̃ -Ã CI exists almost ∼0.11 eV lower than the Ã state, the WP rapidly (within ∼20 fs) decays (almost 98.2%) to the lowest state, X̃ . Unlike any other states, when the WP dynamics is performed on the X̃ state, a very low population transfers to the Ã state (see Figure 5f) and ∼0% to other energy states, because the WP on X̃ cannot access higher energetic CIs, even not the X̃ -Ã CI. b. Anthracene. We have also performed WP dynamics for Ant+ locating the initial WP separately from upper (Ẽ 2Ag) to lowest state (Ã 2B1g), and the corresponding population profiles are presented in Figure 7a−f, respectively. Because the C̃ and D̃ electronic states are coupled to the Ẽ state with almost equal competence and there is no substantial interaction between Ẽ state with lower states, the Ẽ state population is almost equally divided to those states as shown in Figure 7a. The diabatic populations of six electronic states are shown in Figure 7b for an initial location of the WF on the D̃ state. The maximum of the D̃ state population (∼83%) rapidly transfers to the lower energetic C̃ state through C̃ -D̃ CI lying very close to the D̃ state minimum. It has been observed that the molecular dynamics of Ant+ is more complex than Nap+, and all of the calculated reduced densities cannot properly explain the corresponding population profiles. Although the population of the D̃ state is decaying and that of C̃ state is increasing with almost similar fashion, the WP in the D̃ surface shows mild oscillations, but that of C̃ state shows a large oscillation with time. Starting from (q1 = 0, q2 = 0), the WP on the D̃ PES remains intact with mild oscillation in its position (see Figure 8a) even at longer times (see Figure 8b). The C̃ state profile (see Figure 7b) shows a down−up−down at t = 170, 190, and 200 fs, where we have taken the snapshots of corresponding WP. At t = 170 fs, the WP shifts toward the positive q2 direction [see Figure 8c]. It then returns back to its original position at t = 190 fs [see Figure 8d] and finally escapes from the D̃ -C̃ CI [see Figure 8e]. When the dynamics is initiated from C̃ state, only a few population transfers to the lower B̃ state and upper D̃ state. The decay rate of the C̃ state in Ant+ is much slower than that of the Nap+ as shown in Figure 7c. The time evolution of the WP from the B̃ state of Ant+ is shown in panel d of Figure 7. The B̃ state population rapidly decays to two lower electronic states, Ã and X̃ , unlike the case of Nap+ (in Nap+ the decay process is rather slow). The low energy difference between Ã state minimum and Ã -B̃ CI (∼0.1 eV) helps to ultrafast decay the B̃ state population to the Ã state. Because the B̃ -X̃ CI is upper lying from the B̃ state minimum by about 6.68 eV, the Ã state population transfer through this CI is not possible. Thus, the X̃ state is only populated by Ã state population, similar to Nap+, via X̃ -Ã CI. The WP on the B̃ state remains almost in the same place, but on the Ã state it shows oscillations. The general feature of the WP is that it remains in almost the same place [see Figure 8f,g] through which the dynamics is initiated; second, the part of that

WP on the lower surface shows oscillations [see Figure 8h, i], where one or more nodes may appear [see Figure 8j]. The time-dependent population profiles are presented in Figure 7e when the initial WP is located on the FC point of Ã state. The Ã state slowly transfers population to the X̃ state unlike for the Nap+, where the same transfer is rapid (see Figure 5e). As the X̃ -Ã CI exists ∼0.72 eV above the Ã state, the transfer of population from the X̃ state to any upper state is quite difficult. Thus, a few of the X̃ state population transfer to the Ã state as shown in Figure 5f. 5.2.2. Simulation of Complex Photoelectron Spectra. The proper understanding of the contribution of the upper electronic manifold to the lower ones is a contemporary research topic in modern chemistry. The key element of analyzing nonadiabatic dynamics in polyatomic molecules is the study of the vibronic structure of the electronic spectra. Therefore, a quantum-classical dynamical study is carried out to examine the complicated spectral features underlying six electronic states of the two PNHs. The highly overlapping structures in the PE spectra of Nap+ and Ant+ are found to originate due to electronic transitions in the six lowest states (D̃ 0−D̃ 5). The photoelectron spectral data of Nap+ and Ant+ were experimentally recorded by da Silva et al.27 as displayed in Figures 9a and 10a, respectively, which are theoretically simulated by the parallelized TDDVR methodology. For this purpose, the TDDVR propagation scheme is utilized to propagate the nuclear WP on the configuration space to calculate the autocorrelation function, C(t) = ⟨Ψ*(t/2)|Ψ(t/ 2)⟩. The present form of the C(t) is used here for its various advantages.44 The WF has components in six different PESs, which are combined to generate C(t). The photoelectron spectra of the Nap+ and Ant+ are obtained by Fourier transforming the C(t) using the relation I(ω) ∝ ω∫ ∞ −∞C(t) exp(iωt) dt. The PE spectra are calculated by propagating WP for six different initial conditions. The resolution and broadening of the experimental spectra are incorporated into the calculated spectrum by damping C(t) with h(t) = exp[−((| t|)/(τ))2] with an appropriate parameter, τ. Further, the Gibbs phenomenon generated in the calculated spectra is reduced by multiplying a suitable function.45 To compute the PE spectra, we assumed that the neutral ground state ionizes with equal probability to all excited PESs. Because the PE spectra are generated from the contribution of six electronic surfaces, we have carried out dynamics initiating the WP on each of those specific states to calculate the corresponding PE bands separately for each PNH, and finally have convoluted them to obtain the theoretical spectra of the system. a. Napthalene. The main contribution in the X̃ state PE band comes from the vibrational modes ν6 and ν7, where the peaks are ∼0.185 and ∼0.202 eV spaced in energy that correspond very well with corresponding frequencies,33 ω6 (0.1847 eV) and ω7 (0.2014 eV). The ν1 and ν2 modes mainly support the generation of the broad structure of the Ã state. On the other hand, the most important DIB, the irregular B̃ band that makes naphthalene astronomically valuable, comes mainly due to the ring breathing mode (ν1). The tuning modes, ν6 and ν7, are essential to obtain these three lower energetic (X̅ 2Au, Ã 2B3u, and B̃ 2B2g) bands. To capture these bands, technically we have performed three calculations starting with an initial WP Ψ(0) on those electronic states of six-state 21 mode Hamiltonian. Numerical propagation of Ψ(t) then is carried out to acquire the corresponding C(t)’s, which are utilized to compute the individual state partial bands. Equal weighted O



Article

Figure 9. Photoelectron spectrum for the X̃ -Ã -B̃ -C̃ -D̃ -Ẽ electronic manifold of Nap+. The intensity (in arbitrary units) is plotted as a function of energy of the electronic states. (a) Experimental gas-phase photoelectron emission spectrum was recorded by da Silva Filho et al.27 (b) Results obtained by the WP propagation method within the TDDVR scheme considering 21 vibrational degrees of freedom. τ = 64 fs is used here. (c) The PE spectra calculated by MCTDH method involving six-state with 29 modes.35

Figure 10. (a) Experimental photoelectron bands28 pertinent to the coupled X̃ -Ã -B̃ -C̃ -D̃ -Ẽ electronic states of Ant+ within the binding energy range of 6.63−11.16 eV, whereas (b) and (c) panels display the theoretically calculated envelope using TDDVR and MCTDH35 methodology, respectively. τ = 30 fs is used to compute the spectra by our approach. The PE spectra calculated by MCTDH method are involved with six states with 31 DOFs.

It has been clearly seen that in the C̃ -D̃ -Ẽ band the effect of nonadiabatic coupling effect is severe. Although the effects of upper electronic states to the Ẽ band are excluded from our study, it was observed that the ν3+ν4+ν6 modes play a major role to generate this spectra. The vibronic coupling modes such as ν37 + ν39 + ν41 + ν42 also have a high impact to construct the irregular structure of Ẽ PE band. Although at high energy region the TDDVE calculated C̃ -D̃ -Ẽ band shows some finer dissimilarity with the MCTDH calculated one, the TDDVR calculated total PE spectra of Nap+ show peak by peak correspondence with the experimental and MCTDH calculated spectra. b. Anthracene. Wavepacket propagation within the TDDVR framework is initiated on the six-state 22 mode Hamiltonian separately for each electronic surface to simulate the corresponding experimental spectra27 (see Figure 10a). The ν7 and ν9 are found relevant to the X̃ state band. Corresponding peaks, appearing ∼0.19 eV energetically spaced (see Figure 10b for TDDVR simulated spectra), show well agreement with the frequency value (0.1986 eV) of mode ν9. In case of Ã state PE spectra, the tuning modes ν3 and ν9 have significant contributions. Electronic transition of B̃ band has a paramount importance because of its correlation with the observed DIB of

convolution of those partial bands produces the total spectra as displayed in Figure 9b within the energy range from 7.95 to 10.75 eV. TDDVR simulated spectra for Nap+ are presented in Figure 9b in comparison with the experimental (Figure 9a) and MCTDH calculated (Figure 9c) results. It is worthwhile mentioning that for better simulation of the experimental PE spectra, the TDDVR simulated PE spectra are gripped with τ = 64 fs (low value) because the parameter τ controls the broadening and not the positions of the peaks. These X̃ , Ã , and B̃ bands are well separated from each other, indicating that there is substantial energy separation among the states. On the other hand, the excited C̃ , D̃ , and Ẽ states show overlapping structure in the experimental PE spectra27 because they are energetically close and also corresponding CIs are not well separated. The WP dynamics has been carried out for each state separately, and we have calculated each band individually using the damping parameter τ = 64 fs. The main contribution of the C̃ state spectra comes from the tuning modes ν2+ν3+ν6, whereas the vibronic coupling has a significant effect on this band through ν27. The primary structure of the D̃ band comes from ν1+ν7 with two additional peaks from the mode ν6. In addition, the coupling mode, ν42, substantially affects this band. P



Article

astronomical experiments in many observatories. A large number of theoretical attempts32−35 were made to simulate such DIBs to explore the underlying mechanism of subpicosecond decay dynamics. An attempt of quantum-classical dynamical simulation on two PNHs is made to explore such long-standing issues of interstellar physics. The detail dynamical study in quantum framework is restricted by the huge number of strongly coupled nuclear DOFs. Many PESs entangled by vibrational modes in such molecules are vibronically coupled, and the well-known adiabatic or BO approximation is not sufficient to interpret the resulting femtosecond population decay from upper to lower electronic surfaces. Such a nonadiabatic phenomenon gives rise to compound structure in the photoabsorption spectra. Table 1 displays the numerical details of the present calculations. The current dynamical treatment includes the six (most relevant ν1, ν2, ν3, ν6, ν7, and ν9,) totally symmetric vibrational modes and other strongly coupled DOFs. The strongest coupling between the X̃ -B̃ states by ν47 is considered along with ν40, ν44, ν45. The coupling in this case is considered by ν23, ν27, ν29 between Ã -B̃ states. Vibronic coupling is expected to play a major role in the appearance of DIBs. Therefore, such coupling is explicitly included in the dynamics for precise description of the location and shape of the B̃ band of Ant+ in 707.0−711.0 nm range of the recorded DIB. TDDVR calculated time-dependent autocorrelation function is Fourier transformed to calculate the spectrum as presented in Figure 11a along with the corresponding MCTDH analog (note Figure 11b) and the experimental resonance enhanced multiphoton dissociation (REMPD)25 and Ar-matrix spectroscopy24 results (see the inset of Figure 11b). The peak positions of 579, 636.22, and 709 nm obtained by the TDDVR method show good accordance with the 582.4, 639, and 708.5 nm positions calculated by the MCTDH method (see Figure 11b).35 We have also calculated the B̃ state spectra (see Figure 12a) of Nap+, which shows close proximity with the 29D MCTDH spectra (see panel b of Figure 12). The experimental gas-phase photoelectron spectra (see Figure 12d) and absorption spectra based on the MIS technique (see Figure 12c) are also presented in Figure 12 for comparison. Although we have treated the Hamiltonian in reduced dimension, the TDDVR calculated spectra show peak by peak corresponding with the 29D MCTDH spectra.35 Yet the two tiny peaks at ∼10.09 and ∼10.15 eV of the result calculated by MLMCTDH36 are absent in our result because we have not included all of the DOFs in our calculation.

Cernis 52, which is described quite in detail in the next section. The vibrational modes ν2, ν7, and ν9 have dominant progressions to accure this band. Peak spacing of ∼0.77, ∼0.184, and ∼0.19 eV for ν2, ν7, and ν9 modes in the TDDVR simulated spectra shows well agreement with the corresponding ab initio frequency values33 of 0.0772, 0.1823, and 0.1986 eV, respectively. To obtain six X̃ 2B2g-Ã 2B1g-B̃ 2Au-C̃ 2B2g-D̃ 2B3u-Ẽ 2A1g state spectra of Ant+, dynamics is performed using the Hamiltonian of Ant+ initiating the WP separately on those electronic surfaces to calculate partial state spectra. After equal weighted convolution of all of those bands, we obtain the TDDVR calculated theoretical PE spectra as presented in Figure 10b in comparison with the experimental (Figure 10a) and MCTDH calculated (Figure 10c) spectra. It is worth mentioning that TDDVR simulated envelope is convoluted with τ = 30 fs. For C̃ state spectra, the major contribution comes from ν1, ν7 and one of ν8 or ν9. As the ν8 and ν9 are dynamically equivalent, we have included one of them in our study. The ν7 and ν8 show dominant progressions for D̃ state band with peak spacing 0.185 and 0.188 eV, which show good agreement with the corresponding ab initio33 frequencies 0.1823 and 0.1890 eV. For Ẽ state spectra, the major contribution comes from ν1, ν7, and ν9. Note that the ν1, ν2, and ν9 vibrational modes are essential to obtain total PE spectra of Ant+. Although our TDDVR simulated X̃ -Ã -B̃ state envelopes are almost similar to the MCTDH ones within the drawing accuracy, the TDDVR calculated C̃ -D̃ -Ẽ band is quite broad as compared to its exact quantum mechanical analog. In summary, our X̃ -Ã -B̃ -C̃ -D̃ -Ẽ state spectra of Ant+ properly reproduce all of the peaks of the experimental and MCTDH band. 5.2.3. Diffuse Interstellar Bands of Nap+ and Ant+. The DIBs are a special absorption feature present in the IR or visible range of various stars in the interstellar clouds. In astrophysical investigation, the origin and detection of such bands is so far not established. This investigation is fundamentally important to properly understand the evolution of the DIBs in ISM. The electronic transitions of Nap+ and Ant+ situated in the spectral range, where these DIBs appear, are considered as origin of such DIBs. Spectroscopic measurements of the star Cernis 5230 and also HD2811591731 persuade a discovery of three new DIBs in 580−680 and 708.68−709.08 nm appeared in the strongest bands of the PE spectra of Nap+ and Ant+, respectively, in a low temperature experiment. Iglesias-Groth et al.29 have reported the discovery of a broad interstellar or circumstellar band at 708.8 nm that exists in the range of experimental uncertainty of the sharp spectral region of Ant+. A broad interstellar feature is explored at ∼708.8 ± 0.2 and ∼708.494 nm of the spectra of Cernis 5230 and HD281159,31 respectively. The B̃ band of the PE spectra of Ant+ (in the ISM condition made in laboratory) carries the same feature. Theoretical simulation of such DIBs is a great challenge because the underlying mechanism of subpicosecond decay is so far not explored. The construction of Hamiltonian33 using the rigorous ab initio calculation and a detail theoretical spectral study35 was performed previously by Ghanta et al. In our present theoretical investigation, we have simulated such DIBs of Ant+ utilizing the same Hamiltonian used to reproduce the spectral characteristics of those stars and laboratory measurements26 with the help of TDDVR methodology. In the experimental photoelectron spectra30,31 of both Nap+ and Ant+, the B̃ band is the most intense one, and the same band is found in the detectable spectral region of the

6. SUMMARY A state-of-the-art algorithmic development to parallelize the TDDVR method and its application to the two realistic model Hamiltonians of naphthalene and anthracene is presented in this Article. The few important characteristic features make the proposed method dynamically more efficient than the other contemporary methodologies. Because the TDDVR method uses time-dependent grid points, the required number of basis functions is very few to converge a particular dynamical property. As the method can distinguish the quantum and classical nature of the vibrational modes of a multidimensional molecular system, it treats the significant contributing modes quantum mechanically and the dynamically less important modes classically. The present implementation of the TDDVR algorithm shows almost linear speed-up because the matrix multiplication between kinetic matrix and the wave vector involved in the algorithm is highly sparse, and therefore, Q



Article

Figure 12. Comparison between the TDDVR calculated (panel a) B̃ band of Nap+ with 29D MCTDH35 (dashed line in panel b) and 48D ML-MCTDH36 (solid line in panel b) simulated bands, where the corresponding experimental band23 (panel c) is obtained by the MIS techanique. The experimental PE spectra are reported in ref 27. The intensity in arbitrary units is plotted as a function of energy in eV. The theoretical spectrum calculated by TDDVR is captured with τ = 74 fs. Figure 11. (a) The B̃ state REMPD band of Ant+ is simulated by parallel TDDVR method. The experimental REMPD25 and Armatrix24 spectroscopy results are presented in the inset of panel (b). (b) The complete B̃ band calculated by MCTDH method: The spectra captured by two states with 10 vibrational modes (solid line) and six states with 31 modes (dashed line) are shown. The intensity in arbitrary units is plotted as a function of the wavelength in nanometers.

results of the well-established MCTDH47−51 and MLMCTDH36,46 methods.

■

ASSOCIATED CONTENT

S Supporting Information *

Explicit form of the matrices used in the TDDVR quantum and classical equation of motions; such time-dependent matrices have special dynamical features. This material is available free of charge via the Internet at http://pubs.acs.org.

although the dimensions of the matrices are huge, the number of effective multiplication is reasonable. The focus of our formulation to develop an accurate and efficient numerical algorithm with tracing has three main points: (1) the accommodation of many electronic surfaces with vibrational modes; (2) to reproduce different types of experimental spectra and nuclear dynamics; and (3) to carry out multidimensional calculation within a reasonable real time. The numerical performance of quantum-classical dynamical simulation by the TDDVR approach is tested on two realistic model Hamiltonians of Nap+ and Ant+, where 6 electronic states with 21 vibrational DOFs for Nap+ and 22 for Ant+ fabricate a rich structure in the photoabsorption spectra. Such systems provide us a scope to explore the workability of the newly implemented parallel TDDVR algorithm with respect to the observables obtained from MCTDH methodology and experimental measurement. In section 5, we have demonstrated satisfactorily the accuracy and workability of our parallelized TDDVR methodology to simulate the overlapping band structures of PE and REMPD spectra of these PAHs33,35 with respect to the

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS S.S. acknowledges the Council of Scientific and Industrial Research (CSIR), New Delhi, for Doctoral Fellowship and Indian Association for the Cultivation of Science for Super computaion facility in CRAY super computer (AMD Opteron(tm) Processor 23). S.A. acknowledges the Department of Science and Technology (DST, Government of India) for financial support through project no. SR/S1/PC-13/20082011. R



■

Article

The Naphthalene Cation (C10H8+). J. Chem. Phys. 1991, 94, 6964− 6977. (24) Szczepanski, J.; Martin, V. Electronic And Vibrational Spectra Of Matrix Isolated Anthracene Radical Cations: Experimental And Theoretical Aspects. J. Chem. Phys. 1993, 98, 4494−4511. (25) Roland, D.; Specht, A.; Blades, M.; Hepburn, J. Resonance Enhanced Multiphoton Dissociation Of Polycyclic Aromatic Hydrocarbons Cations In An RF Ion Trap. Chem. Phys. Lett. 2003, 373, 292− 298. (26) Sukhorukov, O.; Staicu, A.; Diegel, E.; Rouillé, G.; Henning, T.; Huisken, F. D2 ← D0 Transition Of The Anthracene Cation Observed By Cavity Ring-down Absorption Spectroscopy In A Supersonic Jet. Chem. Phys. Lett. 2004, 386, 259−264. (27) da Silva Filho, D. A.; Friedlein, R.; Coropceanu, V.; Ö hrwall, G.; Osikowicz, W.; Suess, C.; Sorensen, S. L.; Svensson, S.; Salaneckb, W. R.; Brédas, J. L. Vibronic Coupling In The Ground And Excited States Of The Naphthalene Cation. Chem. Commun. 2004, 1702−1703. (28) Sánchez-Carrera, R. S.; Coropceanu, V.; da Silva Filho, D. A.; Friedlein, R.; Osikowicz, W.; Murdey, R.; Suess, C.; Salaneck, W. R.; Brédas, J. L. Vibronic Coupling In The Ground And Excited States Of Oligoacene Cations. J. Phys. Chem. B 2006, 110, 18904−18911. (29) Iglesias-Groth, S.; Manchado, A.; Rebolo, R.; Hernández, J. I. G.; Garcia-Hernández, D. A.; Lambert, D. L. A Search For Interstellar Anthracene Towards The Perseus Anomalous Microwave Emission Region. Mon. Not. R. Astron. Soc. 2010, 407, 2157−2165. (30) Iglesias-Groth, S.; Manchado, A.; Garcia-Hernández, D. A.; Hernández, J. I. G.; Lambert, D. L. Evidence For The Naphthalene Cation In A Region Of The Interstellar Medium With Anomalous Microwave Emission. Astrophys. J. 2008, 685, L55−L58. (31) Galazutdinov, G.; Lee, B.; Song, I.; Kazmierczak, M.; Krelowski, J. A Search For Interstellar Naphthalene And Anthracene Cations. Mon. Not. R. Astron. Soc. 2011, 412, 1259−1264. (32) Reddy, V. S.; Ghanta, S.; Mahapatra, S. First Principles Quantum Dynamical Investigation Provides Evidence for The Role Of Polycyclic Aromatic Hydrocarbon Radical Cations In Interstellar Physics. Phys. Rev. Lett. 2010, 104, 111102(1)−(4). (33) Ghanta, S.; Reddy, V. S.; Mahapatra, S. Theoretical Study Of Electronically Excited Radical Cations Of Naphthalene And Anthracene As Archetypal Models For Astrophysical Observations. Part I. Static Aspects. Phys. Chem. Chem. Phys. 2011, 13, 14523− 14530. (34) Reddy, V. S.; Mahapatra, S. Photostability Of Electronically Excited Polyacenes: A Case Study Of Vibronic Coupling In The Naphthalene Radical Cation. J. Chem. Phys. 2008, 128, 091104(1)− (4). (35) Ghanta, S.; Reddy, V. S.; Mahapatra, S. Theoretical Study Of The Electronically Excited Radical Cations Of Naphthalene And Anthracene As Archetypal Models For Astrophysical Observations. Part II. Dynamics Consequences. Phys. Chem. Chem. Phys. 2011, 13, 14531−14541. (36) Meng, Q.; Meyer, H.-D. A Multilayer MCTDH Study On The Full Dimensional Vibronic Dynamics Of Naphthalene And Anthracene Cations. J. Chem. Phys. 2013, 138, 014313(1)−(12). (37) Hall, K. F.; Boggio-Pasqua, M.; Bearpark, M. J.; Robb, M. A. Photostability Via Sloped Conical Intersections: A Computational Study Of The Excited States Of The Naphthalene Radical Cation. J. Phys. Chem. A 2006, 110, 13591−13599. (38) Sato, T.; Tokunaga, K.; Tanaka, K. Vibronic Coupling in Naphthalene Anion: Vibronic Coupling Density Analysis for Totally Symmetric Vibrational Modes. J. Phys. Chem. A 2008, 112, 758−767. (39) Adhikari, S.; Billing, G. D. A Time-Dependent Discrete Variable Representation Method. J. Chem. Phys. 2000, 113, 1409−1414. (40) Barkakaty, B.; Adhikari, S. Time-Dependent Discrete Variable Representation Method In A Tunneling Problem. J. Chem. Phys. 2003, 118, 5302−5318. (41) Puzari, P.; Sarkar, B.; Adhikari, S. Quantum-Classical Dynamics Of Scattering Processes In Adiabatic And Diabatic Representations. J. Chem. Phys. 2004, 121, 707−721.

REFERENCES

(1) Manz, J.; Wö ste, L. Femtosecond Chemistry; Wiley-VCH: Weinheim, 1994; Part 1, Chapter 2, pp 14−128. (2) Chergui, M. Femtochemistry, Ultrafast Chemical and Physical Processes in Molecular Systems; World Scientific: Singapore, 1996. (3) Chergui, M.; Zewail, A. H. Electron and X-Ray Methods of Ultrafast Structural Dynamics: Advances and Applications. ChemPhysChem 2009, 10, 28−43. (4) Kosloff, R. Propagation Methods For Quantum Molecular Dynamics. Annu. Rev. Phys. Chem. 1994, 45, 145−178. (5) Cerjan, C., Ed. Numerical Grid Methods And Their Application To Schrödinger’s Equation; Kluwer Academic Publishers Group: Dordrecht, The Netherlands, 1993. (6) Leforestier, C.; Bisseling, R. H.; Cerjan, C.; Feit, M. D.; Friesner, R.; Guldenberg, A.; Hammerich, A.; Jolicard, G.; Karrlein, W.; Meyer, H.-D.; Lipkin, N.; Roncero, O.; Kosloff, R. A Comparison Of Different Propagation Schemes For The Time Dependent Schrö dinger Equation. J. Comput. Phys. 1991, 94, 59−80. (7) Sim, E.; Makri, N. Time-Dependent Discrete Variable Representations For Quantum Wave Packet Propagation. J. Chem. Phys. 1995, 102, 5616−5625. (8) Gelman, D.; Schwartz, S. D. Tunneling Dynamics With A Mixed Quantum-Classical Method: Quantum Corrected Propagator Combined With Frozen Gaussian Wave Packets. J. Chem. Phys. 2008, 129, 024504(1)−(7). (9) Thompson, A. L.; Punwong, C.; Martinez, T. J. Optimization Of Width Parameters For Quantum Dynamics With Frozen Gaussian Basis Sets. Chem. Phys. 2010, 370, 70−77. (10) Köppel, H.; Domcke, W.; Cederbaum, L. S. Multimode Molecular Dynamics Beyond The Born-Oppenheimer Approximation. Adv. Chem. Phys. 1984, 57, 59−246. (11) Yarkony, D. R. Diabolical Conical Intersections. Rev. Mod. Phys. 1996, 68, 985−1013. (12) Domcke, W.; Yarkony, D. R.; Köppel, H. Conical Intersections: Electronic Structure, Dynamics And Spectroscopy. Advanced Series in Physical Chemistry; World Scientific: Singapore, 2004; Vol. 15, Chapter 2, pp 41−127. (13) Herzberg, G.; Longuet-Higgins, H. C. Intersection Of Potential Energy Surfaces In Polyatomic Molecules. Discuss. Faraday Soc. 1963, 35, 77−82. (14) Baer, M. Beyond Born-Oppenheimer: Electronic Nonadiabatic Coupling Terms And Conical Intersections; Wiley-Interscience: Hoboken, 2006; Chapter 4, pp 84−104. (15) Baer, M.; Englman, R. A Modified Born-Oppenheimer Equation: Application To Conical Intersections And Other Types Of Singularities. Chem. Phys. Lett. 1996, 265, 105−108. (16) Baer, M. Nonadiabatic Effects In Molecular Adiabatic Systems: Application To Linear Plus Quadratic E⊗e System. J. Chem. Phys. 1997, 107, 10662−10666. (17) Adhikari, S.; Billing, G. D. The Geometric Phase Effect In Chemical Reactions: A Quasiclassical Trajectory Study. J. Chem. Phys. 1997, 107, 6213−6218. (18) Baer, M.; Lin, S. H.; Alijah, A.; Adhikari, S.; Billing, G. D. Extended Approximated Born-Oppenheimer Equation. I. Theory. Phys. Rev. A 2000, 62, 032506(1)−(8). (19) Adhikari, S.; Billing, G. D. Non-Adiabatic Effects in Chemical Reactions: Extended Born-Oppenheimer Equations And Its Applications. Adv. Chem. Phys. 2002, 124, 143−196. (20) Sarkar, B.; Adhikari, S. Extended Born-Oppenheimer Equation For A Three-State System. J. Chem. Phys. 2006, 124, 074101(1)−(18). (21) Yonehara, T.; Takahashi, S.; Takatsuka, K. Non-BornOppenheimer Electronic And Nuclear Wavepacket Dynamics. J. Chem. Phys. 2009, 130, 214113(1)−(14). (22) Yonehara, T.; Hanasaki, K.; Takatsuka, K. Fundamental Approaches To Nonadiabaticity: Toward a Chemical Theory Beyond The Born-Oppenheimer Paradigm. Chem. Rev. 2012, 112, 499−542. (23) Salama, F.; Allamandola, L. J. Electronic Absorption Spectroscopy Of Matrix-Isolated Polycyclic Aromatic Hydrocarbon Cations. I. S



Article

(42) Dirac, P. A. M. Note on Exchange Phenomena In The Thomas Atom. Proc. Cambridge Philos. Soc. 1930, 26, 376−385. (43) Bohm, D.; Hiley, B. J.; Kaloyerou, P. N. An Ontological Basis For The Quantum Theory. Phys. Rep. 1987, 144, 321−375. (44) Sardar, S.; Paul, A. K.; Sharma, R.; Adhikari, S. The Multistate Multimode Vibronic Dynamics Of Benzene Radical Cation With A Realistic Model Hamiltonian Using A Parallelized Algorithm Of The Quantumclassical Approach. J. Chem. Phys. 2009, 130, 144302(1)− (12). (45) Sardar, S.; Puzari, P.; Adhikari, S. Multi-State Multi-Mode Nuclear Dynamics On Three Isomers Of C6H4F2+ Using Parallelized TDDVR Approach. Phys. Chem. Chem. Phys. 2011, 13, 15960−15972. (46) Manthe, U. A Multilayer Multiconfigurational Time-Dependent Hartree Approach For Quantum Dynamics On General Potential Energy Surfaces. J. Chem. Phys. 2008, 128, 164116. (47) Meyer, H.-D.; Manthe, U.; Cederbaum, L. S. The MultiConfigurational Time-Dependent Hartree Approach. Chem. Phys. Lett. 1990, 165, 73−78. (48) Manthe, U.; Meyer, H.-D.; Cederbaum, L. S. Wave-Packet Dynamics Within The Multiconfiguration Hartree Framework: General Aspects And Application To NOCl. J. Chem. Phys. 1992, 97, 3199−3213. (49) Beck, M. H.; Jäckle, A.; Worth, G. A.; Meyer, H.-D. The Multiconfiguration Time-dependent Hartree (MCTDH) Method: A Highly Efficient Algorithm For Propagating Wavepackets. Phys. Rep. 2000, 324, 1−105. (50) Meyer, H.-D.; Worth, G. A. Quantum Molecular Dynamics: Propagating Wavepackets And Density Operators Using The Multiconfiguration Time-Dependent Hartree Method. Theor. Chem. Acc. 2003, 109, 251−267. (51) Meyer, H.-D.; Gatti, F.; Worth, G. A. Multidimensional Quantum Dynamics: MCTDH Theory And Applications; Wiley-VCH Verlag GmbH & Co. KGaA: Weinheim, 2009. (52) Raab, A.; Worth, G.; Meyer, H.-D.; Cederbaum, L. S. Molecular Dynamics Of Pyrazine After Excitation To The S2 Electronic State Using A Realistic 24-Mode Model Hamiltonian. J. Chem. Phys. 1996, 110, 936−946. (53) Burghardt, I.; Giri, K.; Worth, G. A. Multimode Quantum Dynamics Using Gaussian Wavepackets: The Gaussian-Based Multiconfiguration Time-Dependent Hartree (G-MCTDH) Method Applied To The Absorption Spectrum Of Pyrazine. J. Chem. Phys. 2008, 129, 174104(1)−(14). (54) Vendrell, O.; Meyer, H.-D. Multilayer Multiconfiguration TimeDependent Hartree Method: Implementation And Applications To A Henon-Heiles Hamiltonian And To Pyrazine. J. Chem. Phys. 2011, 134, 044135(1)−(16). (55) Vendrell, O.; Brill, M.; Gatti, F.; Lauvergnat, D.; Meyer, H.-D. Full Dimensional (15-Dimensional) Quantum-Dynamical Simulation Of The Protonated Water-Dimer III: Mixed Jacobi-Valence Parametrization And Benchmark Results For The Zero Point Energy, Vibrationally Excited States, And Infrared Spectrum. J. Chem. Phys. 2009, 130, 234305(1)−(13). (56) Meng, Q.; Faraji, S.; Vendrell, O.; Meyer, H.-D. Full Dimensional Quantum-Mechanical Simulations For The Vibronic Dynamics Of Difluorobenzene Radical Cation Isomers Using The Multilayer Multiconfiguration Time-Dependent Hartree Method. J. Chem. Phys. 2012, 137, 134302(1)−(20).

T


Multisurface Multimode Molecular Dynamical Simulation of

Recommend Documents