Car-Parrinello Molecular Dynamics on excited state surfaces

21 Pages·0.25 MB·English

by  Eric R. Bittner

Car-Parrinello molecular dynamics on excited state surfaces
Eric R. Bittner and D.S. Kosov∗
Department of Chemistry
University of Houston
Houston, TX 77204, USA
(February 1, 2008)
This paper describes a method to do ab initio molecular dynamics in electronically excited
systems within the random phase approximation (RPA). Using a dynamical variational treatment
oftheRPAfrequency,whichcorrespondstotheelectronicexcitationenergyofthesystem,wederive
coupled equations of motion for the RPA amplitudes, the single particle orbitals, and the nuclear
coordinates. These equations scale linearly with basis size and can be implemented with only a
9 single holonomic constraint. Test calculations on a model two level system give exact agreement
9 with analytical results. Furthermore, we examined the computational efficiency of the method by
9 modelingtheexcitedstatedynamicsofaone-dimensionalpolyenelattice. Ourresultsindicatethat
1 the present method offers a considerable decrease in computational effort over a straight-forward
n configurationinteraction(singles)plusgradientcalculationperformedateachnuclearconfiguration.
a Pacs: 71.15-m, 71.15.Pd,31.50.+w
J

2
2 I. INTRODUCTION

3
In recent years,there has been considerable advances towards fully ab initio molecular dynamics simulations using
v
1 inter- and intramolecular potentials derived from quantum many-body interactions. This activity has been largely
7 spurred by the introduction of hybrid quantum/classical dynamics treatments, such as the density functional the-
2 ory/moleculardynamics methods introduced by Car andParrinello1–3 (CPMD). In this approach,the quantum elec-
5 tronic orbitals, φ (1) , are coupled to the classical nuclear degrees of freedom, R , through a fictitious Lagrangian,
i n
0 which along wit{h an a}ppropriate set of Lagrange multipliers, (Λ ), leads to a set of equations of motion
ij
8
9
δE[ φ ,R ]
/ µφ¨(1)= { i} n + Λ φ (1),
t i − δφ∗(1) ik k
a i k
m X
MnR¨n = RE[ φi ,Rn], (1)
- −∇ { }
d
where E[ φ ,R ] is the Hohenberg and Kohn energy functional4. This approach has proven to be highly successful
n { i} n
in simulating chemical dynamics in the condensed phase.
o
c However, there are a number of distinct disadvantages to the method. First, the equations of motion involve both
: fast (the orbitals)and slow (the nuclei) degreesof freedom. Thus, the computationaltime-step must be rathersmall.
v
Secondly, imposing the holonomic constraintoforbital orthogonalityforces a single iterationof the approachto scale
i
X as M N2 where N is the number ofparticles andM the size ofthe electronicbasis. Finally,the methodcanonly
∼ ×
r treat systems in the electronic ground state.
a
Several groups have attempted to overcome each of these difficulties. The time-step problem has been attacked
via multiple time-scale methods5,6 and by using the electronic density instead of the Kohn-Sham7 orbitals as the
dynamical variable8. Several algorithms have been introduced to tackle the scaling problem3. In these cases, scaling
more favorable than N3 is achieved only at the cost of several assumptions and approximations.
∼
Recently, Alavi, Kohanoff, Parrinello, and Frenkel have developed a density functional based method for treating
systems with finite temperature electrons9,10. This approach has a number of attractive features, foremost being
the inclusion of fractional occupation of the electronic orbitals. Dynamics in this scheme are isothermal rather than
adiabatic,allowingforincoherenttransitionsbetweenelectronicallyexcitedstates. While the approachworkswellfor
small bandgap systems (such as metals or dense hydrogen), this approach is not useful for systems with larger band
gaps (∆E kT) or photo-excited systems.
≫
Inthispaper,wepresentanalternativeformulationofthe Car-Parrinellomethodfortreatingelectronicallyexcited
systems. We begin by introducing the general computational algorithm for the dynamical optimization of excited
electronicstatesinthepresenceof“classicallypropagating”ions. Thisgeneraltheoreticalschemeisdevelopedwithout
∗
Present address: Department of Chemistry, UMIST,Manchester, M60 1QD,UK
1
specific reference to the description of the electronic excited state. The principal results of this section are a set of
dynamical equations of motion for the electronic amplitudes and nuclear coordinates in a generalized phase space.
We then derive specific equations of motion for the electronic amplitudes and nuclear coordinates under the Random
Phase Approximation (RPA) for the electrons. Finally, we demonstrate the salient features of the approach through
a simplified SU(2) model of interacting electrons and via a realistic model of conjugated one-dimensional polymer
lattices.
II. CLASSICAL EQUATIONS OF MOTION FOR THE DYNAMICAL OPTIMIZATION OF EXCITED
QUANTUM STATES
We begin by developing a general dynamical theory for molecular dynamics in excited electronic states. We first
define the N-electron Fock space as being parameterized by a set of dynamical variables λ(τ) which include the
nuclear positions and any other dynamical variables and ensemble constraints placed upon the system (e.g. pressure,
temperature, volume etc...):
= n ,n ,...(λ) , n =N (2)
1 2 i
F (| i )
i
X
We describe the ground state wave function as an expansion on the basis in :
F
Ψ (λ) = C (λ)n ...n ...(λ) , Ψ (λ)Ψ (λ) =1 (3)
o i 1 i o o
| i | i h | i
1...i...
X
Ψ (λ) is not a Slater determinant in general case.
o
| We niext define the mapping operators,P†, which act on the ground state to produce a new state, ν(λ) , which is
ν | i
orthogonalto Ψ (λ) in the subspace of the Fock space:
o +
| i F
ν(λ) =P† Ψ (λ) , (4)
| i ν| o i
P Ψ (λ) =0, (5)
ν o
| i
ν(λ)Ψ (λ) = Ψ (λ)P Ψ (λ) =0, ν(λ)ν(λ) =1 (6)
o o ν o
h | i h | | i h | i
Ifthe Ψ (λ) is a true groundstate oftheN-particlesystem,the excitedstatesshouldlie inthe orthogonalsubspace
o
11.| i
+
F
Variationaldeterminationsoftheexcitedstatesarethussearchesintheorthogonalsubspacefortheoptimalmapping
operator P† and for the optimal basis of the Fock space n n ...n .. such that the energy functional
ν | 1 2 i i
ε = ν(λ)H ν(λ) ,
ν
h | | i
= Ψ (λ)H Ψ (λ)
o o
h | | i
+( ν(λ)H ν(λ) Ψ (λ)H Ψ (λ) ),
o o
h | | i−h | | i
= Ψ (λ)H Ψ (λ) + Ψ (λ) P , H,P† Ψ (λ) ,
h o | | o i h o | ν ν | o i
=E +ω , (7)
o ν (cid:2) (cid:2) (cid:3)(cid:3)
is at a variationalminimum. We write the excitationenergy, ε , as a functional of both the groundand excited state
ν
and introduce explicitly via the double commutator term the energy gap, ω , between the ground and the excited
ν
state.
According to the calculus of variations, we can derive equations of motion which dynamically optimize the excited
state system. We first define formally a set of “velocities” of the states in the Fock space
d d
ν˙ = ν , Ψ˙ = Ψ , (8)
o o
| i dτ| i | i dτ| i
where τ is not the true physical time but rather a parameter which characterizes the evolution of the amplitudes
through the Hilbert space subject to the orthogonalization constraints given above. The “classical” Lagrangian for
our system is given by
2
1 1 1
= µΨ˙ Ψ˙ + µ′ ν˙ ν˙ + mλ˙2(τ)
o o
L 2 | ih | 2 | ih | 2
ε [0 , ν ,λ]
ν
− | i | i
Λ( Ψ Ψ 1) Γ( ν ν 1), (9)
o o
− h | i− − h | i−
whereΛandΓareLagrangemultipliersarisingfromtheholonomicortho-normalizationconstraintforthegroundand
excitedstatesandwhichallow ν and Ψ to betreatedasindependent variables. The massesµandµ′ arefictitious
o
| i | i
masses andhaveno true physicalmeaning. This Lagrangianleads to the classicalEuler-Lagrangeequations12 for the
wave functions:
d ∂ ∂
L L =0, (10)
dτ ∂|Ψ˙oi − ∂|Ψoi
d ∂ ∂
L L =0. (11)
dτ ∂ ν˙ − ∂ ν
| i | i
The Euler-Lagrange equation dictates the following equations of motion for quantum dynamical optimization in the
Hilbert space:
δε
µΨ¨ = ν +ΛΨ ,
o o
| i −δ Ψ | i
o
h |
δε
µ′ ν¨ = ν +Γν ,
| i −δ ν | i
h |
δε
mλ¨ = ν. (12)
− δλ
This dynamical view is the starting point for the remainder of our calculations.
III. RANDOM PHASE APPROXIMATION FOR ELECTRONIC MOTIONS
For a molecular system with N electrons and n nuclei, the full N-body Hamiltonian for a fixed set of nuclear
positions R , is given as (atomic units are used throughout)
ν
{ }
1 Z 1
H = 2 ν + , (13)
2 ∇i − R r r r
ν i i i
i ν i | − | ij | − |
X XX X
N
1
= h(i)+ , (14)
◦ r r
i i
i ij | − |
X X
where h(i) is the one-body electron operator which includes the electron kinetic energy and the interaction between
◦
the electron and the nuclear cores. We define the single-particle Hamiltonian, H , through either the Hartree-Fock
1
(HF) approximation13 or via the density functional theory of Kohn and Sham7. In either case, H is a functional of
1
the single particle density
occ
ρ(12)= φ∗(1)φ (2). (15)
i i
i
X
The single particle states, φ (1) , are solutions of a nonlinear Schr¨odinger equation,
i
{ }
H [ρ,R]φ =E [R]φ , (16)
1 i i i
| i | i
whichmust be determinedat the startof the calculation. Here E [R]is the orbitalenergy. The HF spin-orbitalwave
i
functions and energies differ from the KS ones in the way the terms in the single-particle Hamiltonian are treated.
We will denote the Slater determinant wave-function constructed from occupied single-particle HF (of KS) orbitals
as HF and note that this state represents either the Hartree-Fock ground state or the Kohn-Sham ground state.
| i
Havingobtaineda basisofsingle-particlestates, wecanwrite the two-bodyHamiltonianinsecondquantizedform:
3
H =E + E :a† a :
HF i iσ iσ
iσ
X
1
+ 2 (ij|v|kl):a†iσa†jσ′alσakσ′ :, (17)
ijkl σσ′
XX
where E is the HF single-particle eigenenergy,
i
(ij v kl)= d1d2φ∗(1)φ∗(2)v(12)φ (1)φ (2), (18)
| | i j k l
Z
is the matrix element of Coulomb electron-electroninteractions,: : denotes a normal orderedproduct and E is
HF
the HF ground state energy. The operators a† (a ) create (anni·h·i·late) single electrons in the iσ HF spin-orbital
state, where σ = 1 is the spin index of the eilσectrioσns. { }
±2
Itisconvenientatthispointtomovetoaparticle/holerepresentationandintroducetheparticle/holequasi-particle
operator,
α HF =0, (19)
iσ
| i
through its action on the HF vacuum, where i = p,h represents particles (occupied states above the Fermi energy)
and holes (vacancies below the Fermi energy). These quasi-particles operators are related to the fermion operators
written above via
a† =θ(E E )α† +θ(E E )α , (20)
iσ i− f iσ f − i iσ
a =θ(E E )α +θ(E E )α† , (21)
iσ i− f iσ f − i iσ
where θ(x) is the Heaviside step function and E is the Fermi energy.
f
We canexpandH interms ofthe particle/holeoperatorsandneglectterms with anoddnumber ofparticle or hole
terms (i.e. anharmonic interaction terms of the form α†ααα and α†α†α†α),
∼ ∼
H =E + E α† α + E α α†
HF p pσ pσ h hσ hσ
pσ hσ
X X
1
+ (p p v h h )
1 2 3 4
2 | |
p1pX2h3h4
α† α† α† α† +h.c.
× p1σ p2σ′ h3σ h4σ′
Xσσ′ (cid:16) (cid:17)
+ (p1h2|v|h3p4) α†p1σα†h3σαh2σ′αp4σ′
p1hX2h3p4 Xσσ′
− (p1h2|v|p3h4) α†p1σα†h4σ′αh2σ′αp3σ
p1hX2p3h4 Xσσ′
+anharmonic terms, (22)
where p and h refer to particle and hole single particle states.
Note that each “harmonic” term in our expansion of H involves an even number of particle and hole operators.
We can form a new set of exciton operators by taking bilinear combinations of the particle/hole operators with total
angular momentum S and projection m ,
S
1 1
A† (Sm )= σ σ Sm α† α† ,
ph S h2 12 2| si pσ1 hσ2
σX1σ2
1 1
A (Sm )= σ σ Sm α α , (23)
ph S h2 12 2| Si hσ2 hσ1
σX1σ2
where s σ s σ Sm isaClebsch-Gordoncoefficient. Ournotationissuchthatσ denotestheactionoftimereversal
1 1 2 2 s
h | i
operator on the hole state,
4
α† =( )12+σα† . (24)
hσ − h−σ
The A† (Sm ),A (Sm ) serve to create (or destroy) particle/hole exciton with spin S and projection m .
{ ph S ph S } s
Writing just the “harmonic” part of H in terms of exciton operators,
H = E α† α + E α α†
2 p pσ pσ h hσ hσ
pσ hσ
X X
1
(p p v h h ) ( )S A† (Sm )A† (Sm )+h.c.
− 2 1 2| | 3 4 − p1h4 s p2h3 s
p1pX2h3h4 SXms (cid:16) (cid:17)
+ 1 (p h v h p ) ( )S +1 ( )S′ +1 A† (Sm )A (S′m′)
2 1 2| | 3 4 − − p1h3 s p4h2 s
p1hX2h3p4 SXmsSX′m′s(cid:0) (cid:1)(cid:16) (cid:17)
(p h v p h ) A† (Sm )A (Sm ), (25)
− 1 2| | 3 4 p1h4 s p3h2 s
p1hX2p3h4 SXms
where A (Sm ) = ( )S+msA (S m) is the time reversal exciton annihilation operator. This defines our “quasi-
ph s ph
− −
harmonic” exciton Hamiltonian.
The exact commutation relation between the exciton operators is
[ A (Sµ),A† (S′µ′)]
ph p′h′
1 1 1 1
= σ σ λµ σ′ σ′ λ′µ′
h2 12 2| ih2 12 2| i
σX1σ2σX1′σ2′
×(αhσ2α†h′σ2′δpp′δσ1σ1′ −α†p′σ1′αpσ2δhh′δσ2σ2′). (26)
We apply the Wick theorem to the right hand side of Eq. 26 and then neglect the normal ordering contribution to
this relation. This is the RPA14 and in this approximation the commutation relation between the exciton operators
is bosonic
[Aph(Sµ),A†p′h′(S′µ′)]≈δSS′δµµ′δhh′δpp′. (27)
We then perform a canonicalBogolubovtransformation15 for both singlet(S =0)and triplet (S =1) excitonbosons
Q†(00)= Xi A† (00) Yi A (00),
i ph ph − ph ph
ph
X
Q†(1m )= Xi A† (1m ) Yi A (1m ), (28)
i s ph ph s − ph ph s
ph
X
where Xi and Yi are the RPA amplitudes which are yet to be determined. These operators create and destroy
ph ph
harmonic excitations in the electronic degrees of freedom. Physically, the singlet RPA excitations correspond to
“breathing”modesoftheelectroncloudandtripleRPAexcitationscorrespondtodipolaroscillationsoftheelectronic
systemaboutthe nuclearframe. Inwhatfollows,wewillfocus ourattentionuponthe tripletexcitations,astheseare
the states created upon photo-excitation of a singlet ground state. Equations of motion for singlet RPA excitations
can be obtained as well.
The commutation relations for the RPA excitation creation operators goes as follows:
[Q ,Q† ]= Xi ∗Xi′ Yi ∗Yi′
µi µ′i′ ph ph− ph ph
ph
X
=δii′δµµ′. (29)
By defining the RPA vacuum as
Q RPA =0, (30)
µi
| i
and the RPA excited state as
Q† RPA = ω , (31)
µi| i | µii
5
the orthogonality condition of Eq. 29 is fulfilled automatically for single RPA excitations
RPAQ Q† RPA = RPA[Q ,Q† ]RPA
h | µi µ′i′| i h | µi µ′i′ | i
=δµµ′δii′. (32)
TheRPAexcitationsaretrueoscillationsofelectronicsystemsinthesensethattheHeisenbergequationsofmotion
for the excitation creation operators is precisely as we expect for a harmonic system,
[H ,Q† ]=Q† ω , (33)
2 µi µi i
where ω is the RPAfrequency. From this we can derive the RPA amplitudes, Xi and Yi , and RPAfrequencies ω .
i ph ph i
After a bit of algebra we arrive at the equations:
EphXpih+ (pp′|v|h′h)Ypi′h′ − (ph′|v|p′h)Xpi′h′
p′h′ h′p′
X X
=ω Xi , (34)
i ph
EphYi(ph)+ (pp′|v|h′h)Xpi′h′ − (ph′|v|p′h)Ypi′h′
p′h′ h′p′
X X
= ω Yi . (35)
− i ph
In matrix form:
A B X X
B A |Yi =ω | Yi . (36)
(cid:18) (cid:19)(cid:18) | i (cid:19) (cid:18)−| i(cid:19)
The diagonal blocks, A, are determined by Coulombic interactions between single particle/hole excitations.
Aphp′h′ =Ephδpp′δhh′ (ph′ v p′h). (37)
− | |
However, the off diagonal terms,
Bphp′h′ =(pp′ v h′h), (38)
| |
are related fluctuations about the HF vacuum and thus give rise to a net polarization of the HF vacuum. Thus, the
RPAvacuumdiffersfromtheHFvacuumbytheadditionofthesepolarizationfluctuations. Lastly,wenote,thatifwe
takeY 0 orassume thatB 0 andneglectthe fluctuationterms,the RPAequations reduceto the Tamm-Dankoff
→ ≈
or configuration interaction equations with single excitations.
The matrix RPA method canbe alsoderived as the small amplitude limit of the time dependent HF equation16,14.
We emphasize that the time-dependent density functional theory17,18 is equivalent to matrix RPA if we also assume
harmonic oscillation of the time-dependent density ρ(t) around equilibrium. The standard practical realization of
the time-dependent density functional theory is the time-dependent local density approximation19 where nonlocal
exchange interaction is approximated by local density dependent interaction in adiabatic limit. In this respect the
matrixRPAmethodwe useis advantageousbecause this methodcantreatnonlocalFockinteractionsinnaturalway.
We can transform this into an eigenvalue problem by multiplying the right hand side of this equation by the Pauli
spin matrix, σ ,
z
MΨ =ωσ Ψ , (39)
z
| i | i
where
X
Ψ = | i , (40)
| i Y
(cid:18) | i (cid:19)
is a two componentspinor state. The eigenvectorsolutionto this equationcanbe obtained by transformingthe RPA
matrix equation
σTMΨ =ω Ψ . (41)
z | i | i
Thus, Ψ is the right-handed eigenvector of the antisymmetric matrix
| i
6
A B
σTM= . (42)
z B A
(cid:18)− − (cid:19)
The RPA eigenstates are subject to the orthogonality condition
Xi ∗Xj Yi ∗Yj =δ . (43)
ph ph− ph ph ij
ph
X
Thus, we define ΨR as the right-handed eigenvector
Xi
|ΨRi i= |Yii , (44)
(cid:18) | i (cid:19)
and ΨL as the left-handed eigenvector
Xi
|ΨLi i= | Yii =σz|ΨRi i. (45)
(cid:18)−| i (cid:19)
Theinnerproductrelationbetweenthe leftandright-handedeigenvectorsproducesthedesiredorthogonalityrelation
for the RPA states.
ΨL ΨR =ΨL ΨR =δ . (46)
h i | j i i · i ij
These states are determined at the start of the calculation in a basis of single particle/hole states which spans the
Fock space.
The RPA frequency, ω , is determined by taking the expectation value of M between the left and right-handed
i
eigenvectors.
ω [X,Y,φ,R]= ΨL MΨR
i h i | | i i
A B X
=(X∗, Y∗) k (47)
k − k B A Yk
(cid:18)− − (cid:19)(cid:18) (cid:19)
A. Equations of motion for the excited states molecular dynamics in the RPA
The energy functional for the RPA equations is
E (i)= RPAQ HQ† RPA ,
ex h | µi µi| i
= RPAH RPA
h | | i
+ RPA Q , H,Q† RPA ,
h | µi µi | i
=E +ωh, h ii (48)
o i
where E is the RPA ground state energy and
o
ω = RPA Q , H,Q† RPA , (49)
i h | µi µi | i
h h ii
is the RPA frequency. In the previous section the RPA equation is derived via the Heisenberg equation (Eq.33) of
motion for the excitation creation operators. The Heisenberg equation is physically the same as the requirement of
a variational minimum for the RPA excited state. To demonstrate this, we multiply from the left with an arbitrary
variation of the RPA excited state:
RPAδQ H,Q† RPA =ω RPAδQ Q† RPA .
h | µi µi | i ih | µi µi| i
h i (50)
We can rewrite this equation as a variational of the double commutator in Eq. 49, because RPAQ† =0
h | µi
7
δ
RPA Q , H,Q† RPA ω ( RPAQ Q† RPA 1) =0 (51)
δQ h | µi µi | i− i h | µi µi| i−
µi
n h h ii o
where ω plays a role of the Lagrange multiplier and insures the normalization of the RPA excited states.
i
Neglecting vacuum polarization effects, i.e. approximating RPA by the HF Slater determinant 0 in Eq. 51
HF
we obtainthe so-calledequationof motion method which was|originiallyproposedby D.J. Rowe14 an|d widiely used in
the studies ofmolecularexcitedstatesinquantumchemistry20–22. Thesamematrix RPAequationcanbe derivedby
several methods, e.g. by Green function methods13,23, linear response function method24,25 Each of these methods
hasthestrongpointsinspecificaspectsbutisdeficientinotherrespects. AswedemonstratedbyEq.51ourmethodis
equivalenttotheRayleigh-RitzvariationalprinciplefortheRPAexcitationamplitudesandtheuseofthevariational
approachis pivotal to formulation classical equations of motion for the dynamical optimization of the excited states.
Finally, we assume that polarization of the ground state vacuum due to particle-hole interaction is weak, then the
RPA vacuum is approximately the same as the HF or KS vacuum and we can write the RPA ground state energy as
E E .
o HF
≈
Following these considerations, we can write the Lagrangian
1 1
= ν d1φ˙(1)2+ µ X˙k 2
L 2 | | 2 | ph|
i Z k ph
X X X
1 1
+ µ Y˙k 2+ m R˙2
2 | ph| 2 n n
k ph n
X X X
E [φ,R] ω [φ,X,Y,R]
HF i
− −
+ Λ ( d1φ∗(1)φ (1) δ )
ij i j − ij
ij Z
X
+ Γ ( Ψ Ψ δ ) (52)
kl k l ij
h | i−
kl
X
where the constraint
Γ ( Ψ Ψ δ )
kl k l kl
h | i−
kl
X
= Γ Xk Xl ∗ YkYl ∗ δ
kl ph ph − ph ph − kl
kl ph
X X
  (53)
insures the orthogonality of the RPA excited states. Solving the Euler-Lagrange equations for the single particle
amplitudes, φ (1),
k
δE [φ,R] δω [φ,X,Y,R]
νφ¨(1)= HF i
i − δφ∗(1) − δφ∗(1)
i i
+ Λ φ (1) (54)
ij j
j
X
TheequationsforthesingleparticleorbitalsaresimilartotheCPequationswrittenabove(Eq.1)andincludeaterm
duetothe electronicexcitations. Inthe caseofacollectiveexcitation,inwhichthe RPAexcitationisdelocalizedover
a number of single particle states (or in other worlds, the oscillator strength of the RPA excitations is distributed
more or less uniformly over a number of particle-hole excitations) the variation of the energy of collective excitation
ω with an infinitesimal change in one of the single particle orbitals will be very small. Because of this, the RPA
excitationvariableswillevolveonaslowertime-scalethanthesingleparticlevariablesandwecaninvokeanadiabatic
separation between the RPA variables and the single particle states. That is to say
δω δE
k HF
, (55)
δφ ≪ δφ
j j
so that if one deals with collective electron motions, Eq. 54 can be approximated with very high degree of accuracy
by the CP equations.
8
δE [φ,R]
νφ¨(1)= HF + Λ φ (1). (56)
i − δφ∗(1) ij j
i j
X
Thisshouldbeareasonableapproximationformanyrealmolecularsystemswhereonehasdelocalizedvalenceelectrons
and relatively strong residual Coulomb interaction between them. We emphasize that our computational scheme is
not restricted by the case of collective electronic excitation, i.e. by the approximation (Eq. 55). Taking into account
the explicit expression26 for the functional derivatives of ω in respect to the molecular orbitals φ our dynamical
k i
equations can be straightforwardlyused to do molecular dynamics on “non-collective” excited state surfaces.
The dynamical equations for the RPA amplitudes are determined to be (assuming only a single RPA excitation in
the system),
δω [φ,X,Y,R]
µX¨k = k +Γ Xk
ph − δXk∗ kk ph
ph
=EphXpkh+ ((pp′|v|h′h)Ypk′h′ −(ph′|v|p′h)Xpk′h′)
p′h′
X
+Γ Xk , (57)
kk ph
δω [φ,X,Y,R]
µY¨k = k Γ Yk
ph − δYk∗ − kk ph
ph
=EphYpkh+ (pp′|v|h′h)Xpk′h′ −(ph′|v|p′h)Ypk′h′
p′h′
X(cid:0) (cid:1)
Γ Yk. (58)
− kk ph
or in matrix form
X¨ A+Γ B X
µ ||Y¨ii!=(cid:18) B A−Γ(cid:19)(cid:18)||Yii(cid:19) (59)
One important point to consider is that the functional derivatives of E with respect to the particle states (i.e.
HF
single particlestates abovethe Fermienergy)vanishsince the HF density matrixdoes not include these states. As we
seen below, we need these states and their energies to compute the Coulomb matrix elements in the RPA dynamical
equations. However, all is not lost since we have at hand the states below the Fermi energy and hence the ground
state density. We can thus reconstruct and simply diagonalize H [ρ] to obtain the single particle excited states.
1
Finallyforthe nuclearpositionswehaveclassicalequationsofmotiononthe excitedstatepotentialenergysurface,
δE [φ,X,Y,R]
µR¨ = ex
n − δR
n
δE [φ,R ] δω [φ,X,Y,R ]
HF n k n
=
− δR − δR
n n
ThedetailsofthederivationofeachofthefunctionalderivativesaboveareprovidedintheAppendixA.Thederivation
of the molecular dynamics equations of motion in the RPA compose the central result of this paper.
IV. NUMERICAL AND ANALYTICAL EXAMPLES
A. The dynamical optimization of the excited state on a two level model
Inordertogetanindicationofhowourapproachworks,weapplyourschemetoamodeltwo-levelsystemconsisting
of N distinguishable electrons labeled by p = 1,2,...,N. Each of these can occupy one of two orbitals, a lower and
upper, having energies 1ε(R) and 1ε(R) and distinguished by g =1 and g =2 respectively. The coupling constant
2 −2
V(R) does not depend on any quantum number. The gapbetween orbitalε(R) and the strength of interaction V(R)
depends on a classical variable R. The HF approximation is obtained by setting V = 0, such that the lower level is
fully occupied and the upper is empty, i.e. N =Ω where Ω is the degeneracy of the levels. This model encapsulates
9
the salient effects one expects to see in physically realistic systems. The model Hamiltonian possesses the SU(2)
symmetry27 and can be written in terms of quantum quasi-spin operators and the classical variable R:
1
H =ε(R)Jˆ V(R) Jˆ Jˆ +Jˆ Jˆ
z + + − −
− 2
MR˙2 1 (cid:16) (cid:17)
+ + K(R R )2 (60)
0
2 2 −
where the operators Jˆ ,Jˆ ,Jˆ are defined as follows:
+ − z
Ω
1
Jˆ = a† a a† a ,
z 2 2p 2p− 1p 1p
Xp=1(cid:16) (cid:17)
Ω
Jˆ = a† a
+ 2p 1p
p=1
X
†
Jˆ = Jˆ . (61)
− +
(cid:16) (cid:17)
Above a† ,a are particle creation and annihilation operators on the lower (g =1) or the upper (g =2) level. K is
gp gp
the spring constant which does not allow to classical system to move far from the equilibrium.
Place the classical variable at the position
R=R +q, (62)
0
whichisthesumofthe equilibriumpositionR andthedisplacementq. Assumingthe displacementfromequilibrium
0
is small, we can expand the Hamiltonian in q
dε
ε(R)=ε(R )+ q+O(q2) (63)
0
dR
(cid:12)R0
(cid:12)
V(R)=V(R0)+ dV (cid:12)(cid:12) q+O(q2). (64)
dR
(cid:12)R0
(cid:12)
(cid:12)
Truncating this at first order gives a linear interactionbetw(cid:12)een classical and quantum systems. Neglecting the terms
of O(q2) we get the following quantum many-body Hamiltonian which depends parametricallyon the displacement q
of the classical variable:
dε
H(q)= ε(R )+ q Jˆ
0 z
dR
(cid:12)R0 !
(cid:12)
1 d(cid:12)V
V(R )+ (cid:12) q Jˆ Jˆ +Jˆ Jˆ
0 + + − −
− 2 dR
(cid:12)(cid:12)R0 !(cid:16) (cid:17)
Mq˙2 1 (cid:12)
+ + Kq2 (cid:12) (65)
2 2
In the RPA, the operators Jˆ ,Jˆ behave like pure bosonic operators
+ −
Jˆ Jˆ
− +
, =1 (66)
"√N √N#
Using these bosons, we construct the RPA excitation creation operator (in the two-level model there is only one
possible RPA excitation)
1
Q† = (XJˆ YJˆ ) (67)
+ −
√N −
The orthogonalityofthe RPAexcited state imposes the followingnormalizationconditionon the RPA amplitudes:
10