Surface hopping simulations on photoexcitation dynamics of conjugated polymer

Highlights

• Derive necessary equations for the surface hopping method based on PPP Hamiltonian and CIS formalism.

• Find that an excited state may be falsely populated using the Min-Cost method, and give a solution to this problem.

• Study the contribution of Z matrix to the nonadiabatic couplings and forces.

Abstract

Surface hopping is a very active research area in the recent decades of years. In this work, we derive all necessary equations for a surface hopping method based on Pariser-Parr-Pople (PPP) Hamiltonian and configuration interaction singles (CIS). The surface hopping method is then used to study the dynamics of a cis-polyacetylene chain after it is photoexcited to a high-lying excited state. We find that an excited state may be falsely populated because it is impossible to detect the trivial crossings between noninteracting excited states at time zero using the traditional Min-Cost method.This problem can be solved by screening the initial lattice geometries that do not match a certain condition. By comparison, we discuss the contribution of molecular orbital response with respect to atom displacements (Z matrix) to the nonadiabatic coupling terms and forces, and find that the contribution of Z matrix can be neglected, to some degree. This can greatly save the computational cost.

Introduction

In recent decades, surface hopping method aroused much attention in scientific community and was in rapid developing in the field of nonadiabatic molecular dynamics [1]. Surface hopping method is the so-called semiclassical (mixed quantum-classical) approach, where the nuclear motions are described by classical mechanics and the electrons quantum mechanically. The nuclei evolve on a certain potential energy surface (PES), but the transitions between different adiabatic states are allowed. Ehrenfest method is another nonadiabatic molecular dynamics methods, where the nuclei evolve on an effective PES which is an average of all adiabatic states involved weighted by their populations (also called mean-field method).

The relaxation process of a molecular system from a high-lying excited state to lower excited state is controlled by the strength of nonadiabatic coupling terms (NACTs) between different excited states. Some excited states have strong nonadiabatic couplings, some are noninteracting, depending on the characters of their wave functions and the energy level difference between them. When two noninteracting excited states are approaching degeneracy, the nonadiabatic coupling between them presents sharp peaks that are strongly localized in time. They actually cross each other, called “trivial" crossings. These trivial crossings may lead to “unphysical" hoppings and long range energy transfer in numerical simulations, and they must be avoided.

Different methods have been proposed to deal with the trivial crossing problem [2], [3], [4], [5], [6], [7], [8], [9]. Our group proposed a method in which the excited states are classified into two groups according to their transition density matrices and different group of excited states are noninteracting [9]. Because of this, the system can separately evolve within one group of excited states, and the trivial crossings between noninteracting states are completely eliminated. However, this method only applies to those systems whose excited states can be classified into noninteracting groups, such as a single polymer chain. For those systems composed by interacting polymer chains, it is not applicable. Another method is the so-called Min-Cost algorithm, which is used in the present work [4]. The key point of this algorithm is that the identities of excited states must be followed during dynamics based on the characters of their wave functions, instead of the energy-ordering criterion. In other words, at each time step the excited states are tracked according to a maximum overlap principle so that the old states can be mapped to the new states. However, there is a problem here that can very easily be ignored by people. At time zero, we are not able to track the excited states like at any other time and have to assign them based on energy-ordering criterion.This may lead to fake populations of an excited state during dynamical simulations. We will discuss this problem and give a solution to it in the present work.

People often have to consider how to reduce the computational load because the computational cost is a bottleneck for most surface hopping simulations. In surface hopping simulations, the most time-consuming parts are calculating the NACTs and the forces exerting on atoms. The higher the level of excited-state calculations, the more expensive the calculating NACTs and forces. In general, configuration interaction singles (CIS) is thought to be the zeroth-order approximation for excited-state calculations. It is generally used in surface hopping approaches [10], [11]. The time-dependent density functional theory (TDDFT) is also an efficient approach for calculating excited states of complex molecular systems, which has similar computational cost with CIS. TDDFT is also widely used in surface hopping approaches [12], [13]. However, even for CIS and TDDFT, surface hopping simulations on conjugated polymers are still expensive.

In order to speed up the surface hopping simulations, people have developed many approximation methods, especially on calculation the NACTs since they are the most time-consuming step. Pittner et al. proposed a scheme to neglect all NACTs between other states than the current state.[14] For a surface hopping simulation involving N_s states, taking into account that NACT σ_JJ = 0 and σ_IJ = σ_JI, N_s(N_s − 1)∕2 coupling terms should in principle be evaluated. Using this scheme, the number of coupling terms that should be computed is reduced to N_s − 1. Fabiano et al. neglected the second term of the nonadiabatic coupling vector expression, since they thought the first term dominate the expression in the region of a conical intersection [2]. Nelson et al. thought that the contribution of MO response to the NACTs, reflected by the Z matrix, is small. They truncated the summation up to the number of involved excited states when calculating the Z matrix [11]. Lobaugh and Rossky used approximate CIS wave functions to calculate NACTs and forces. These wave functions are generated using only 50 single excitation configurations [10].

In this work, we fully introduce the surface hopping method which is based on one-dimensional Pariser-Parr-Pople (PPP) Hamiltonian and CIS formalism. This method is easy to implement and has low computational cost [9], [15], [16]. In this method, we deduced analytic expressions for the forces and the NACTs. Using them, the computation speed is greatly accelerated. In order to further reduce the computational cost, we discussed the contribution of Z matrix to see if it is possible to neglect it.

In Section II, we fully describe the implementation of the surface hopping method. In Section III, the results of numerical simulations of photoexcitation dynamics in cis-polyacetylene are presented. Finally, in Section IV we conclude our investigation.