Transcript Polyspherical Description of a N-atom system Christophe Iung LSDSMS, UMR 5636

Polyspherical Description of
a N-atom system
Christophe Iung
LSDSMS, UMR 5636
Université Montpellier II
e-mail : [email protected]
Collaboration avec
Dr. Fabien Gatti, Dr. Fabienne Ribeiro et G. Pasin
Pr. Claude Leforestier (Montpellier)
Pr. Xavier Chapuisat et Pr. André Nauts
(Orsay et Louvain La Neuve)
EXAMPLES OF SYSTEMS
H
nnCH
F
F
F
(H2O)n
FIT of a Potential Energy Surface
(PES) to describe the water solvent
CF3H
Intramolecular Energy transfer
in an excited System :
Dynamical Behaviour of an
Excited system :
Is it ergodic or selective?
Schrödinger ro-vibrational Equation
1- Born-Oppenheimer Approximation :
===> The Potential energy surface V can be expressed
in terms of (3N-6) internal coordinates that
describe the deformation of the molecular system
2- A Body-Fixed Frame (BF) has to be defined :
Tc = Tc(G : 3 coordinates) + Tc(rotation-vibration:3N-3 coordinates)
3- Ro-Vibrationnal Schrödinger Equation : an eigenvalue equation
H |Y> = (Tc+V) |Y(3N-3) internal coordinates> = Ero-vibrationnal|Y>
Problem to be Solved
1- Choice of the set of coordinates adopted to describe the system :
A crucial Choice
2- Expression of the Kinetic Energy Operator (KEO) Tc
3- Calculation and Fit of the Potential energy Surface (PES), V, a function
of the 3N-6 internal nuclear coordinates.
4- Definition of a working basis set in which the Hamiltonian is
diagonalized, this basis should contain 150000 states, for instance.
5- Schrödinger Equation to be solved
-Pertubative Methods (CVPT...)
- Variational method (VSCF, MCSCF, Lanczos, Davidson,...)
6- Comparison between the calculated and experimental spectrum
1-Choice of the set of coordinates
Curvilinear
Rectilinear
Large amplitude motions
Low energy spectrum
DY1
DY2
DZ1
Dq
DZ2
Very Simple Expression of
the Kinetic Energy Operator
Basis of the traditional
Spectroscopy
More Intricate expression of
the KEO
2- Expression of the KEO (Tc)
with
Px = - i h
i

,
xi
* We need an exact expression of the KEO adapted to
the numerical methods used to solve the Schrödinger equation.
* We have to know how to act this operator on vectors
of the working basis set.
3- Analytical Expression of the PES
calculated on a grid (of few thousands
points). (Fit of this function)
Potential Energy Surface
Coordinate 1
Coordinates 2
Outlines of the talk
1- KEO Expression
2.1 : Historical Expressions of the KEO
2.2 : More Recent (1990-2005) Strategies that provide KEO
operator
2- Polyspherical Parametrization of a N-atom System
(IJQC review paper on the web)
2.1 : Principle
2.2 : Application to the study of large amplitude motion
2.3 : Application to highly excited semi-rigid systems :
Jacobi Wilson Method
3- Direct Methods that solve the Schrödinger Equation
3.1 : Lanczos Method
3.2 : Block Davidson Method
4- Application to HFCO
1- Some Famous References
B. Podolsky, Phys. Rev. 32,812 (1928)
E.C. Kemble “The fundamental Principles of Quantum Mechanics”
Mc GrawHill, 1937
E.B. Wilson, J.C. Decius, P.C. Cross “Molecula Vibrations”
McGrawHill, 1955
H.M. Pickett, J. Chem. Phys, 56, 1715 (1971)
A. Nauts et X. Chapuisat, Mol. Phys., 55, 1287 1985
N.C. Handy, Mol. Phys., 61, 207 (1987)
X.G. Wang, E.L. Sibert et M.S Child, Mol. Phys., 98, 317 (2000)
Quantum Expression of KEO for J=0 in the
Euclidean Normalization
2Tc = (tpx)+ px
where pxi is the conjugate momentum associated with
the mass-ponderated coordinates
If a new set of curvilinear coordinates qi (i=1,…,3n-6) is introduced
q = J-1 x

px=t(J-1) pq
where J is the matrix which relies
the cartesian coordinates to the new set of coordinates qi
The determinant of J is the Jacobian of the transformation denoted by J
dtEuclide = dx1 dx2… dx3N-6= J dq1 dq2… dq3N-6
Tc expression of the KEO for J=0
in Euclidian normalization
If 2Tc = (tpx)+ px
and
px=t(J-1)pq
2Tc = (tpq)+ J-1 t(J-1) pq
2Tc = (tpq)+ g pq
 det(g)=J-2
What is the adjoint of pqi ?
It depends on the normalisation chosen
In an Euclidean Normalization
(pqi)+ = J-1 pqi J

where J est the Jacobian
2Tc = J-1 tpq J g pq
Démonstration de (pq)+=J-1 pq J en normalisation euclidienne
Définition de l’adjoint de pqi ?
<(pqi)+ j | f >= < j | pqi f >
Or
... pq (J j f) dq1 dq2… dq3n-6 = 0
si (J j f) s ’annule sur les bornes d ’intégration
... pq (J j* f) dq1 … = ... pq (J j*) f dq1…dq3n-6 + ... J j* pq ( f) dq1…dq3n-6
d ’où
... (J-1pq J j)* f J dq1… dq3n-6 = ... j* pq ( f) J dq1... dq3n-6 =0
... (J-1pq J j)* f dtEuclide = ... j* pq ( f) dtEuclide =0
d ’où
(pq)+ = J-1 pq J
Other way to find 2Tc = J-1 tpq J g pq
Let use the expression of the Laplacian in spherical coordinates :


D = J
J g k ,l
qk
ql
k ,l =1
3 n -6
-1
where J is the Jacobian
2Tc = -h/2p D,
This expression can be re-expressed by
2Tc = (tpq)+ g pq
Quantum Expression of Tc for J=0
in Wilson Normalization dtWilson =dq1 dq2… dq3n-6
This normalization can be helpful to calculate some
integrals.
(jEuclide)* ÂEuclide fEuclide dtEuclide =
(jEuclide)* ÂEuclide fEuclide J dq1 … dq3n-6 =
(jWilson)* ÂWilson fWilson dtWilson
(jWilson)* ÂWilson fWilson dq1 … dq3n-6
(J0.5jEu)* (J0.5 ÂEuJ-0.5) (J0.5 fEu) dtWilson =
(jWilson)*
ÂWilson
(jWilson)* ÂWilson fWilson dtWilson
fWilson
2TcW = J0.5 TcEuJ-0.5 = J0.5 J-1 tpq J g pq J-0.5 = J-0.5tpqJ g pq J-0.5
2 TcWilson = J-0.5 tpq J g pq J-0.5
2 TcWilson = J-0.5 tpq J g pq J-0.5
OR
2TcEuclide = J-1 tpq J g pq
Rectilinear Description
J, g do not depend on q
2Tc =tpq g pq
No problem for Tc
but problem for the fit of V
and for the
physical meaning of q
Curvilinear Description
J, g depend on q
Problem of no-commutation
More Intricate expression
To find and to act on a basis
But easy fit of V
et better physical meaning
of q
Different strategies developed :
Application of the Chain Rule
Handy et coll. (Mol. Phys., 61, 207 (1987))
Starting with the expression with cartesean coordinates
: 2Tc = (tpx)+ px
The chain rule is acted (with the kelp of symbolic calculation)
and provides :
2 Tc = S gkl pk pl + S hk pk
in Euclidean Normalization
 qj
hk =  2
i =1  xi
3N
with
2
Other normailization can be used…
But it results more intricate expression of the KEO Tc
Other formulation :
Pickett expression: JCP, 56,1715 (1972)
Starting from 2 TcWilson = J-0.5 tpq J g pq J-0.5
One can find
2 TcWilson =
tp
q
g pq + V’
V’ « extrapotential term » that depends
on the masses. It can be treated with the potential
1  ln J  ln J 1  2 ln J 1 g kl  ln J
V ' =  - g kl
2
qk ql
4  qk
4 qk  ql
k ,l =1 8
3n -6
This formulation has be exploited by E.L. Sibert et coll. in his
CVPT perturbative formulation: J. Chem. Phys., 90, 2672 (1989)
Ideal features of a KEO expression
1- Compact Expression of the KEO :
larger is the number of terms, larger is the CPU time
2- Use of a set of coordinates adapted to describe the motion of atoms
in order
*to reduce the coupling between these coordinates
* to define a working basis set such that the Hamiltonian matrix is sparse
3- The numerical action of the KEO must be possible
and not too much CPU time consuming
4- The expression should be general and should allow to treat a large
variety of systems
2- Polyspherical Parametrization
The N-atom system is parametrized by (N-1) vectors described by their
Spherical Coordinates ((Ri,i, ji), i=1,...,N-1)
The General Expression of the KEO is given in terms of either
1- the kinetic momenta
associated to the vectors

And the (N-1) radial conjugate momenta pRi = -i
Ri
===> adapted to the description of large amplitude motion
OR
2- the momenta
(p R , pq , pj )
i
i
i
conjugated with the polyspherical
coordinates ((Ri,i, ji),i=1,...,N-1)
===> adapted to the description of highly excited semi-rigid systems
Development of this parametrization
First description of its interest :
X. Chapuisat et C. Iung , Phys. Rev. A,45, 6217 (1992)
Review papers :
F. Gatti et C. Iung,J. Theo. Comp. Chem.,2 ,507 (2003) et
C. Iung et F. Gatti, IJQC (sous presse)
Orthogonal Vectors : F. Gatti, C. Iung,X. Chapuisat JCP, 108, 8804 (1998), and
108, 8821 (1998)
M. Mladenovic, JCP, 112, 112 (2000)
NH3 Spectroscopy : F. Gatti et al , JCP, 111, 7236, (1999) and 111, 7236, (1999)
Non Orthogonal Vectors : C. Iung, F. Gatti, C. Munoz, PCCP, 1, 3377 (1999)
M. Mladenovic, JCP, 112, 1082 (2000);113,10524(2000)
Semi-Rigid Molecules : C. Leforestier, F. Ribeiro, C. Iung 114,2099 (2001)
F. Gatti, C. Munoz and C.Iung : JCP, 114, 8821 (2001)
X. Wang, E.L.Sibert and M. Child : Mol. Phys, 98, 317(2000)
H.G Yu, JCP,117, 2020 (2002);117,8190(2002)
HF trimer : L.S. Costa et D.C. Clary, JCP, 117,7512 (2002)
Diatom-diatom collision : E.M. Goldfield,S.K. Gray, JCP, 117,1604(2002)
S.Y. Lin and H. Guo, JCP, 117, 5183(2002)
“ORTHOGONAL” SET OF VECTORS
BF Gz
H
C
H
O
F
C
H
Jacobi Vectors
Radau Vectors
2T = 
i
Pi + Pi
i
Polyspherical Coordinates :
R3, R2, R1,
q1, q2
et  (out-of-plane dihedral angle)
O
Non Orthogonal Set of Vectors
H
C
O
BF Gz
H
Valence Vectors
2T = 
i, j
Pi + Pj
Mij
Polyspherical Description : R3, R2, R1, q1, q2 and 
- M matrix determination M (Trivial)
- Dramatic Increase of term number… CPU can dramatically
increase
Determination de la Matrice M
Any set of vectors can be related to a set of Jacobi vectors :
R = A r Jacobi
if 2T = 
i
(Pi Jac ) + Pi Jac
i
 M = At  A
= 
i, j
Pi + Pj
Mij
La Matrice M est une matrice très facile à déterminer et dépendant
des masses
Elle permet de généraliser les résultats obtenus avec les vecteurs
orthogonaux
Developed expressions of the KEO
+
Starting from 2Tc = P M P;
*kinetic momentum Li associated with
Ri and the radial momenta

(Li , PR = - i
)
Ri
i
By using
PR Ri
L R
Pi =
+ i 2 i
Ri
Ri
i
Pi is substituted by
*Conjugate radial and angular
momenta



PR = - i h
, P = -i h
, Pj = - i h
Ri
i
ji
Obtained by the
substitution of the angular momentum
i
{
i
i
* Lix = - sin j i p i + cos j i pj i
* Liy = cos j i p i + sin j i pj i
* Liz = pj i
A BF (Body Fixed) frame has to be defined to introduce
the total angular momentum (full rotation) vector J
Choix du Body Fixed
The (Gz)BF is chosen parallel to RN-1 ; LN-1 is substituted by
This requires 2 Euler rotations (,)
The last Euler rotation () can be chosen by the user
In general, RN-2 is taken parallel to the plane (Gxz)BF
But other choice can be done :
N atoms = 3N-3 degrees of freedom
•Kinetic Momenta Li (i-1,...,N-2)
(2N-5) angles (jN-1=j N-2=q N-1=0)
•the (N-1) radial conjugate
momenta
• the full rotation J (3 angles)
*(3N-6) conjugate momenta
(jN-1=j N-2=q N-1=0)
{PRi (i = 1,..., N), P i (i = 1, ..., (N -1)), Pj i (i = 1,..., (N - 2))
*the full rotation J (3 angles)
By taking into account the fact that
RN-1 and RN-2 are linked to the BF frame
(problem of no-commutation of the operator
that depends on vectors RN-1 and RN-2 )
It results in general expression of the KEO
with
One finds that :
The problem of no-commutation are such that
KEO developed expression for a system described
by a set of (N-1) orthogonal vectors
KEO developed expression for a system described
by a set of (N-1) orthogonal vectors
General Expression of Tc in terms of the conjugate momenta
Associated with the polyspherical coordinates
Expression used to study semi-rigid systems
F. Gatti, C. Munoz, C. Iung, JCP, 114, 8821 (2001)
The expression of the KEO are known…
How can we use them
for instance for semi-rigids systems ?
1- Orthogonal Coordinates provides rather simple expression of KEO…
However, these coordinates does not necessary describe a real
deformation of the system
2- Interesting coordinates, such valence coordinates, are not
‘orthogonal’ The KEO expression is intricate
Two sets of coordinates can be used…
This is the idea of the Jacobi-Wilson Method
Definition of “Curvilinear Normal Coordinates”,Qi,
In terms of polyspherical coordinates qj :
0
Hvib
1 3N- 6
0
=
 (qn Fn,m qm + pn Gn,m pm )
2 n,m
where Gn,m = G(qeq ) and
0
  2V 
Fn,m = 

q
q
 n mq
eq
Normal Modes
Defined in terms of
Polyspherical coordinates
-1
Q= L q
Polyspherical
Coordinates
Pq is substituted by (tL) PQ in Tc
Advantages :
Simplicity of Tc in terms of polyspherical coordinates
Physical Interest of the Normal Modes
Jacobi-Wilson Method
(C. Leforestier, A. Viel, C. Munoz, F. Gatti and C. Iung, JCP, 114, 2099 (2001))
JACOBI-WILSON STRATEGY
H
C
O
JACOBI
F
Jacobi Vector
• Description Polyspherique
• Simple Expression of the KEO T
• Normal Mode Coordinates :
WILSON
Q = L-1 q
• Definition of a working basis set :
1 ...3N-6
H = H 0 + DV + DT
DIAGONALIZATION
Application to HFCO et H2CO
Up to10000 cm-1
Improvement of the zero-order basis set
On can take into account to the diagonal anharmonicity:
2
2

0,
Hvib = 2 + V (Q )Q   = Qeq
2  Q
0
H = H + DV + DT
DV = V - Vanharm
h2
DT = 2
3 N -6

n,m


0
(Gn,m - Gn,m )
qn
qm
H Matrix calculation

q
DV
semi-analytical estimation of its action
pseudo spectral scheme used
DG
Spectral Representation :
 (Q1 ...Q6 ) =  Yn ...n n (Q1 )  ...  n (Q6 )
Grid Representation:
n1 ...n6
1
6
1
6
 (Q1a ...Q6 f ) =  Yn ...n  n (Q1a )  ...  n (Q6 f )
n1 ...n6
1
6
1
6
Ideal features of a method that provides
eigenstates and eigenvalues which can be located
in a dense part of the spectrum
•Application to a large variety of systems ;
•Use of huge basis set ;
•Obtention of eigenvalues and eigenstates;
•Control of the accuracy of the results ;
•Small CPU time, Small memory requirement ;;
•Easy to use and to adopt ;
•Specific Calculation of energies in a given part of
the spectrum ;
LANCZOS METHOD
Iterative Construction of the Krylov subspace
generated by {un, n=O,N} :
1- Initialization : A first guess vector u0 is chosen
2- Propagation : The following vector un+1 is calculated
n+1 un+1 = (H – n) un – n un-1
with
n = <un|H|un>
n+1= <un+1|H|un>
LANCZOS FEATURES
• Lanczos Method:







H
dim B


  0 1 0

   

1
2
 1

 0

 


N
N
 
dim N < dim B
• Avoid the determination of the full H matrix.
• The convergence is slower when the state density
increases
DIAGONALISATION DE H
Ouverture de Fenêtres en
énergie
300
r ( pour 100 cm-1)
250
Méthode de
Lanczos
200
150
100 Diagonalisation
directe
50
0
0
5000
10000
Energies en cm-1
E0
15000
One has to open some window energy
• Spectral Transform .
Lanczos applied to G=(Eref°-H)-1. or exp(-(H-Eref)2)
Modified Block-Davidson Algorithm to calculate
a set of b coupled eigenstates
Method based on one parameter : e which sets the accuracy
F. Ribeiro, C. Iung, C. Leforestier JCP in press
C. Iung and F. Ribeiro JCP in press
Prediagonalization step in order to reduce
the off-diagonal terms
The working basis set Banh is divided into :
Banh = P°  Q .
Where P° contains the zero-order states which play a significant rôle
in the calculation performed :
H is diagonalized in P°, et this new basis set {u°i ,E°i } is used during
the Davidson scheme
H°=
E°1
0
0
E°i
P°
0
Eanh1
0
0
Q
0
Eanhq
Determination of the Block of states
We can defined the block of states using the
second-order perturbation
 ,  = 
2
qQ
u H uq
u  H uq
E - E
E - E
States such that
un a given block :
0
0
q
0
 ,   0.1
2
0
are retained
APPLICATION TO HFC0
• Faible barrière de
dissociation (14000 cm-1)
HFCO
HF + CO
• Mode de déplacement
hors du plan très découplé à
haute énergie.
• Forte densité d ’états
Selectivity of the energy transfer in HFCO
whose out-of-plane mode is excited
6 modes
2981 cm-1 : CH stretch
1837cm-1 : C=0 stretch
1347 cm-1 : HCO bend
1065 cm-1 : CF stretch
662 cm-1 : FCO bend
1011 cm-1
Excitation of the
Out-of-plane mode
In-plane
modes
O
C
H
F
Out of plane mode
Moore et coll. have studied the highly excited
out-of-plane overtones (nnout-of-plane, n=14,…,20) :
they predict the localization of energy in these states
How can we understand that
a highly excited state can be localized in one mode while the
state density is large for Eexc=14000-20000cm-1?
CONVERGENCE OF THE DAVIDSON SCHEME
-État 60
-État 120
-État 180
Lanczos
Error
On the
Energy (cm-1)
Nombre
d ’itérations
Davidson
Davidson
Iteration Number
• Between 1 and 60 Davidson iterations are required to calculate
each state.
. DE is not a correct indicator of the convergence
Convergence criterion
• IIq(M)II constitue an excellent
Indicator of the convergence
And the accuracy of the
eigenenergy and eigenstate.
qM = ( Er(MM ) - Hˆ )Yr(MM )
Erreur (cm-1)
•I
: IIq(M)IImax < 10 cm-1
Pour DE < 0.01 cm-1
IIq(M)IImax < 50 cm-1
Pour DE < 0.5 cm-1
Davidson Iteration Number
Numerical Cost
8000
7000
BLOC-DAVIDSON
6000
LANCZOS
Nombre
d ’actions de H
5000
4000
3000
2000
1000
0
2000
2500
3000
3500
4000
4500
5000
Energie d ’excitation (cm-1)
5500
6000
6500
Determination of an Active Space specificaly built
to study a given state
Example : state |10n6> in HFCO
(State n° 1774, in a 100,000 state primitive basis set
Caracterized by vmax(in plane)=8 and Emax=32,000cm-1)
We begin with a Davidson calculation performed
on |10n6>(0) in a 7,000 state basis set defined
by vmax(in plane)=4 and Emax=24,000cm-1
The Davidson scheme in this small basis set converges and
provides an estimation of the eigenstate studied: |10n6>(1)
The largest |v1,…,v6>°contributions of this
estimated eigenstate |10n6>(1) are retained
in a small (368 states) « active space » Po used
in the calculation in the large (100,000 states) basis set
Application to the calculation of highly
Excited overtones in HFCO
State 10n6: State n°1700
Zero Order Energy
P0 dimension
Prediagonalized
Energy of the guess
Converged Energy
Overlap with the
guess
Main Contributions
of the eigenstate
v1averaged
v2averaged
v3averaged
v4averaged
v5averaged
v6averaged
6n
6439.6
1
6439.6
8n
8641.9
1
8641.9
10n
10849.9
368
10127.1
6018.4
0.77
7984.9
0.64
9948.8
0.92
0.1 4n
+0.77 6n
+0.28(n1+4n)
-0.22(n3+n)
+0.1(n2+4n)
+0.1(n4+n)
-0.1(n1+n)
0.27
0.08
0.14
0.06
0.02
5.47
- 0.2 n
+ 0.64 8n
-0.30(n1+n)
-0.24(n3+n)
-0.20(n2+n)
+0.1(n4+n)
-0.1(n1+n)
0.47
0.11
0.31
0.20
0.23
6.75
Similar coefficients
- 0.31 n
obtained for different
+ 0.51 10n
overtones :
-0.2 (n1+n)
The nature of the
-0.23 (n3+10n)
coupling
-0.20 (n2+n)
is identical
+0.1 (n4+10n)
for these overtones
-0.1 (n1+10n)
0.69 CH stretch
The CH stretch
0.17
is the more
0.37 HCO bend
coupled mode
0.15
0.04
8.57
Main features of this new
Prediagonalized Block-Davidson Scheme
•It can be coupled to any method which can provide the action of the
Hamiltonian on a vector
Huge basis set can be used (more than 100 000 states)
•Calculation of the eigenstates and eigenenergies
•The accuracy of the results can be controled (with ||qM||)
•Low memory cost
•Faster and more efficient than Lanczos
•Very easy to use because it depends only on one parameter e
•It is adapted to calculate a series of coupled states
Conclusions
The development of a general method to calculate high
excited ro-vibrational state is crucial
Different approachs have to be exploited :
• The Jacobi-Wilson method coupled with the Davidson
algorithm presents interesting advantages.
•It allows the specific calculation of eigenstates
associated with highly excited states
It can be improved by using a MC-SCF or SCF treatment
• However a lot of work has to be done… improve the fit of
V, use a fit of the KEO in order to reduce the CPU time…
3 – MCTDH Method (Time Dependant method)
A fit of the global PES has to be performed :
A factorized form of H is required
Spectrum calculation
•By Fourrier transform of the survival probability
•By filtered diagonalization
References :
H-D Meyer, U. Manthe and L. Cederbaum, Chem. Phys. Lett. 165, 73 (1990)
M. Beck, A. Jaeckle, G. Worth and H-D Meyer, Physics Reports, 324,1 (2000)
C. Iung, F. Gatti and H-D Meyer, J. Chem. Phys., 120, 6992 (2004)
The MCTDH Approach
•Primitive Basis set : { |v1,…,v9>°,vi=0,1,…,vimax }
•The MCTDH « Active Space » is generated
by the configurations
: { ji1(1)(Q1,t).. … ji9(2)(Q9 ,t); ij=0,1,…,ijmax ; j=1,…,9} :
(time-dependant functions which are adapted to the
location of the wave-packet describing the system)
It is efficient if ijmax< vimax
Time Dependent Coefficient and Functions to estimate
max
max
3
i3
Y (Q ,..., Q , t ) =  ... A , , (t ) j ( ) (q
, t)
1
9
j1 j 2 j3
j
j

 =1

j1
j3
i1
Time dependent coefficient to optimize Time Dependent 3D functions to optimize
.
i AJ =
 <  J I H I L  AL
L
Projection on the
Space generated by
Functions jj() (Qk)
. ( )
ij
= (1 - P
( )
) (r
( )
Density Matrix
-1
) <H
( )
j
( )
Mean Field
Hamiltonian
Application to the calculation of highly
excited overtones in HFCO (|nn6>)
Spectrum obtained with filtered Diagonalization
Fraction of Energy in the different Normal Modes
10n6
10n
- 0.31 n
+ 0.51 10n
+0.1 (n4+10n)
-0.23 (n3+10n)
-0.20 (n2+n)
-0.2 (n1+n)
-0.1 (n1+10n)
20n6
The CH stretch
is an
energy reservoir