- Email: [email protected]
The determination of magnetisabilities using density functional theory
Volume 2 17. number 3
CHEMICAL PHYSICS LETTERS
14 January 1994
The determination of magnetisabilities using density functional theory Susan M. Colwell and Nicholas C. Handy University Chemical Laboratory, LensJeld Road, Cambridge CB2 IE W, UK Received 24 June 1993; in final form 23 November 1993
The theory for the calculation of magnestisabilities using current density functional theory is presented, following on from the original theory of Vignale, Rasolt and Geldart. It is shown that a Hellmann-Feynman identity holds for the first derivative of the energy with respect to a magnetic field. The magnetic Hessian which arises in the coupled perturbed Kohn-Sham equations contains an extra term involving the current density operator as well as the density itself. It is now in a form for computation.
1. Introduction The application of density functional methods to molecular problems has made significant advances recently with the development of new functionals, and it is a field that is growing rapidly. Calculations of energies, structures and vibrational frequencies are well developed, and in general are more accurate than Hartree-Fock self consistent field (SCF) calculations at similar cost. SCF calculations of electric and magnetic properties are known to be inaccurate, and correlated methods, particularly for frequency-dependent properties, are complicated and expensive. The result is that only a limited number of such calculations have been performed, and those only on small molecules. It is therefore important to consider density functional theory (DFT) methods for obtaining these properties. The authors have recently presented DFT calculations of polarisabilities, cr, and hyperpolarisabilities, /3, of molecules [ 1 ] using coupled Kohn-Sham theory and Senatore and Subbaswamy [ 21 have reported calculations of cr and second hyperpolarisabilities, y, for atoms. Many other workers have used finite field DFT calculations to obtain electric properties of molecules. Although there have been calculations of magnetic properties using DFT (e.g. Salahub et al. [ 3 ] ) these have not been based on a full coupled Kohn-Sham treatment. The difficulty, which is widely acknowledged, is that
with standard DFT using an exchange-correlation functional which only depends on the density, the corresponding exchange-correlation part of the magnetic Hessian vanishes. Therefore it is necessary to use a density functional theory which introduces the current density and specifically introduces an exchange-correlation functional which depends on both densities. The basic theory for this current density functional theory (CDFT) has been presented by Vignale, Rasolt and Geldart [ 41, and they have given the specific form of the Kohn-Sham equations. In this Letter, we discuss the evaluation of magnetisabilities. It is necessary to derive the appropriate coupled perturbed Kohn-Sham (CPKS) theory for current density functional theory. We present an explicit expression for the magnetic Hessian within this theory using the local density approximation. We find that a Hellmann-Feynman identity holds, and hence the expressions for second-order properties are straightforward.
2. The Schriidinger equation in an external electromagnetic field In this section we present dard electromagnetic theory so that subsequent results can context. The non-relativistic
0009-26141941% 07.00 0 1994 Elsevier Science B.V. All rights reserved. SSDZ OOOOS-2614(93)El399-2
an exposition of stanfor completeness, and more easily be put into Hamiltonian for a sys271
14 January 1994
CHEMICAL PHYSICS LETTERS
Volume 2 17, number 3
tern of N electrons and Nz nuclei moving in an external electromagnetic field (E, B) is given (within the Born-Oppenheimer approximation) by
where 3 is the Fock operator. This leads to the definition of the Fock matrix, which is diagonal at convergence if canonical orbitals are used. (6)
~,,=(PIhIq)+2(pqIjj)-(~jI~q)=~P~,,
(1)
where veXt= -e@, and @is the scalar electrostatic potential such that E= -V@-dA/dt, and A is the magnetic vector potential such that B=VxA. For a uniform field VeXt=eE.r and, in the Coulomb gauge
where we use the standard repeated suffix notation, with i, j, ... being occupied orbitals, a, b... being virtual orbitals and p, q, ... being any molecular orbitals. The expression for the RHF energy is given by fT=2(ilhli)+2(iiljj)-(ijlji).
(7)
In the presence of an electromagnetic field, the energy expression is still as above, but now
(V.A=O), A=bBxr. We define the density and current operators acting on the ground state wavefunction P(rr, r2, .... rN) by j?= U*(r) Y(r)
=N
J
!P(r, r,, .... rN) r2, .... rN) dr, drs...drN ,
x Y(r, f 3p=-
2
[
!P(r)V!P(r)-V!P(r)Y(r)]
(2)
(8)
=ho@+h@, where the V-A as an operator, the derivatives external field series,
is not to be interpreted as div A, but i.e. the term is V. (A@). Now consider of the energy with respect to a real 9’. Express the energy as a Taylor
(3)
(note that our definitions do not include any factors of e), i.e. the dependence of !Pon all but one spatial variable has been integrated out. The ground state energy is then given by
J
p(r)
J=$+ +e
J
Write the derivative
(;,4*+V,,>dr
jp,(r)-Adr ,
where 8. is the ground of an external field.
of a molecular
orbital as
(10) (4) or alternatively
state energic in the absence
@p=@; + @u;“,@: + $%YJ$@
+...
(11)
so that uL;,=@pl;9=o.
3. Coupled-perturbedHartree-Fock theory for electric and magnetic properties Consider closed-shell restricted Hartree-Fock (RHF ) and a system of iV= 2M electrons. The wavefunction is a single Slater determinant of M doubly occupied orbitals &, which are obtained as the M lowest eigenfunctions of the eigenvalue equation 9@i=ci@i, 272
(5)
(12)
Orthonormality of orbitals implies that 0 is skew Hermitian. Also invariance with rotations within the occupied and virtual spaces mean that without loss of generality we can chose &&=o.
(13)
Differentiation of the energy with respect to the field P and using eq. (6) yields the usual HellmannFeynman identity,
f$
=2(il$
Ii) .
d.+@ /=o
=2(il
dh df”
LO I‘) 1 ,
(15)
where I i) now refers to @p,i.e. the unperturbed orbitals. Second-order properties are given by d26 dp’d9
9 ‘=O
)
P=O
(22) (23)
Im ( U$)H2ai,bj = 0
and so U is purely real. The energy is expressed as - $EAE,~A, +... ,
#u;= -2Pi
Ia)
(24)
(25)
is the permanent dipole moment, and from ( 16) (YAP = - 4 U$ P$
+2iIm(U{i) dh
U&)H,+j=O
Pti+Re(
where
dh IO +(‘I -I d.9” ;9=o >
dFA 9=. -I
where P’;, is the dipole integral. The CHF equations give
6= 8. -E&
d2h Ii) 9=0 =2(i’ dY”dF6’ 99=.
+2Re(U$)
-(‘I
(21)
(14)
First-order properties are therefore obtained from db
14 January 1994
CHEMICAL. PHYSICS LETTERS
Volume 2 17. number 3
la)
IO . >
(16)
To evaluate these properties we need the U matrices, which are obtained from the coupled Hartree-Fock equations. These equations are obtained by differentiating the condition Fai=O 3
with respect to Ta and evaluating at 9=0
(17) to give
is the polarisability. For a uniform magnetic field, the first-order perturbation to the Hamiltonian is pure imaginary (- (i!ze/2m) (A-V+V*A), to be interpreted as before), and YA=BA, (rl g
IB=ols)
(@,I~~B~o~v+v~~~B~oldJs,
=-g i#ie
=- G
I
(18)
-where Hlai,bjand Hzai,bjare the electric and magnetic Hessians, respectively, as first presented by Stevens, Pitzer and Lipscomb [ 5 1, and are given by H,,i,aj=4(aiIbj)-(ajlbi)-(abIji) + (e=-+5ai,bj, Hzai,bj=(aj(bi)-(~blii)+(~“-t~)~=ni,bj
V@), dr=ihr;l,,
2@(rx
(27)
where iN = eL, and N is real and skew Hermitian. The CHF equations give Re(U$)H1,,,=O
+ i Im ( U&)Hzni,bj= 0 ,
(26)
(28)
3
Nb+Im(U&)Hm,bj=O
3
(29)
and so U is purely imaginary. The energy is expressed as ~=$-E~ml-fEnE,~~r+...,
(30)
where (19) (20)
for perturbations out of real orbitals. For a uniform electric field, the perturbation to the Hamiltonian ( =eE*r) is real, 9@=EA,
m,=-2iNi=0
(31)
is the magnetic moment (which vanishes within RHF), and &- --2Qf/‘++Im(U$)N$
(32) 273
Volume 2 17, number
3
is the magnetisability,
CHEMICAL
PHYSICS
14 January
LETTERS
where
1 e2 p(r,Mr2)
J[p] = -4ne, 2 s is the Coulomb energy is
These expressions are equivalent Stevens et al. [ 5 1.
Colwell et al. [ 1 ] and others have implemented coupled Kohn-Sham theory (which is exactly analogous to the above but using Kohn-Sham theory rather than RHF) and successfully calculated static polarisabilities and hyperpolarisabilities. However, a straightforward application of ordinary density functional theory to magnetic properties does not work because the exchange correlation part of the magnetic Hessian vanishes. Partly in order to overcome this kind of problem, Vignale, Rasolt and Geldart [ 41 have introduced current density functional theory where the exchange-correlation energy is a functional, not only of the density p(r) but also of the currentj, (r). They have also formulated the theory in terms of the gauge invariant vector y(r) , where
Ij,(r)/Ar)l
,
(34)
and we follow their theory in this section, before applying it in the next section to get new formulae for electric and magnetic properties. In the presence of an electromagnetic field the current density functional theory expression for the energy of the system previously.discussed (4) becomes
(35)
(37) 2
energy. The non-interacting
kinetic
where x’[p, j,] is a single Slater determinant of the M doubly occupied non-interacting orbitals which are the lowest eigenfunctions of the appropriate orbital equation. The generalised Hohenberg-Kohn theorem states that gvA is minimised when p and j, take the correct values for the potentials A and VeXt, and this minimum value is the true ground state energy for the system. The appropriate orbital equation for the non-interacting orbitals must have the form
(
&
where the current P=
V&r)
( -ifiV+eA,e)2+
>
@i=Ei@i,
(39)
and density are now given by
ye; 3
(41 and the effective potentials found. The energy is
+
Aeff and
Vcrr are to be
I w(r) dr+ I pV,,(r) dr+&W,l
which can be written
+ &
s
+&,[pJ,l 3
where u(r) is the nuclear 274
dr
’
,
(42)
as
(A’-A&)pdr+e2
j, * (A -4d
dr
+J[PI + s pv(r) dr+ s pK,,(r> dr
where
+Ecc[P,jpl
dr
to those given by
4. Current density functional theory
v(r)=Vx
Ic-r2l
1994
attraction
potential,
.
(43)
(36)
Taking the functional derivative of (43) with respect to p, (6&,/6p) lip, gives rise to the definition
and
of v,ff,
Volume 2 17, number 3
CHEMICAL PHYSICS LETTERS
I4 January 1994
be substituted back into the orbital equation (39) and the energy expression (43). + T/H(r) + ~x’x,(r)3
(44)
and the functional derivative with respect to j,,
(Wc/6jp)[P,j,l Ipgives 4ff(r) =A(r) +4,(r) ,
(45)
where
5. Energy derivatives in current Kohn-Sham
We now rewrite this in standard Kohn-Sham [ 6,7 ] language. The standard ground state energy expression is &=2(iIhIi)+2(iiIjj)+Ex,[p],
(46) (47)
theory
(53)
where the usual exchange term of HF theory is replaced by the exchange-correlation energy functional. Within current density functional theory this becomes, from eq. (43 ) , 8=2(iIh,,Ii)+2(iiIjj)+Ex,,[p,
v]
(48) (54)
The orbital equation becomes
where we define
- ~V2+v(r)+VH(r)+Vs(r) +
e*
~A*+Kx,-
ifie 2m
(A-V+ V-A)
>
+i=~igi,
- E
(49)
and, given E,, [p, j,] , this can in principle be solved to give the orbitals, and hence the density, current and energy. Vignale, Rasolt and Geldart [ 41 now recast the problem in terms of the gauge invariant vector y(r)=Vx v,(r)/p(r)]. It can be shown that Exc[~Jpl =E&,
~1 ,
(50)
i.e. the exchange-correlation energy is actually a functional of p and v and so is gauge invariant. Taking the functional derivatives of (50) with respect to p and v gives new expressions for V,,(r) and A,,(r), SE
V,,(r) = $
[p,v]l,-e’*,
(A,,~V+V*A,,)
=h+h,,
= i Vx $$
[P,Vllp,
(55)
The Kohn-Sham matrix is now
(56)
In the rest of this section we consider the derivatives of the energy with respect to a real external field .@, and show that a Hellmann-Feynman identity holds. We follow a similar procedure to the RHF case (section 2) with the same definition of 0 etc. Working within the LDA approximation, i.e. assuming there are no Vp or V-j, dependent terms in the potential, dE m =2(il $Fj li)+2(il
s
Ii)
(51) +2&[FY-(ilV,,Ir)]+
eA,,(r)
.
(52) -e
giving a gauge invariant expression for V, which can
s
dlip(r)~&l df”
s$$dr dr ’
(57)
275
Volume 2 17, number 3
CHEMICAL PHYSICS LETTERS
where
(58)
substituting Hermitian, CL!?
m
Note
(59) Then
in (57) and remembering one obtains
8 is skew
Ii).
(65)
Note that this is dh/d@ not dh,,/dP, i.e. we have verified that the Hellmann-Feynman identity holds for current density functional theory, as we would expect from the Kohn-Sham condition Fy ~0. Therefore first- and second-order electric and magnetic properties are given by equations which are formally the same as in the RHF case ( 15) and ( 16), but I i) now refers to @p, i.e. the unperturbed KohnSham orbitals. To evaluate these we need the U matrices, which are now obtained from the coupled-perturbed Kohn-Sham equations, described in the next section.
6. Coupled-perturbed current Kohn-Sham
but d& dp
=2(il$
14 January 1994
(61)
The coupled-perturbed Kohn-Sham equations obtained by differentiating the condition Fy=O
and
(62)
theory are
(66)
with respect to FA and evaluating entiating (56) gives
at 9~0.
Differ-
so
(63)
Then
dV
+@I xc dp”
14) )
and it can be seen that
276
(67)
14 January 1994
CHEMICAL PHYSICS LETTERS
Volume 2 17, number 3
choosing oij= oa:ob=0 as before
dl/,,
(69)
(” d3@ 6=. IO
where (77)
For perturbations out of real orbitals j, is first order in 4 and so is A,, and so the last term in (68) can be ignored when evaluating at T= 0. Also, from (62) (71)
dhxc _ li)=4iRe(&j)e
(al s
9-O x
k
d(Axca) =lp
X
dp
@bbQj
dr-4iIm(
obj>e
9F=o
J
k~ia~l,=,E~&~$).
(78)
So the coupled perturbed Kohn-Sham equations have the same form as the CPHF equations ( 18 ), but now Hlai,bj and Hzoi,bj, the electric and magnetic Hessians respectively, are given by
(72)
hence +4e
JjaiylF=?ejdr
(79)
(73)
(74)
If the orbitals are initially real, jj,: ikj, is pure imaginary, i.e. kjr= -kv is real. Hence di, dP
p=o
= -4Im(
Od)kj,
(75)
,??tF*[p, v] =EkFA[p]o+
and therefore dv, =-41m(Uti)VX d@ y=o
where the j matrices have been reintroduced in place of ik. Vignale et al. [ 41 give an expression for the LDA exchange-correlation energy in the case of a weak field (which may be the applicable case). It appears that it can be written
kjr 0
p
,
(76)
J
F(p”‘)
1v(r) I’dr
and sufficient data is given to enable F(p’/‘)
(81) to be 277
Volume 2 17, number 3
CHEMICAL PHYSICS LETTERS
fitted. Specifically for this form of the functional, setting F(Y)
=g(P) 3
(82)
(83)
E:,,(r) = ~xo(r> +&T’(P) I y(r) 12,
-L(r)
= j Vx
KWpPl
(84)
14 January 1994
where the j matrices have been reintroduced in place of ik. HI and Hz are now real, and so the same deductions can be made as in the RHF case. So for purely electric fields U is real and one recovers the usual CKS equations. For purely magnetic fields U is pure imaginary, and one obtains
hence, as v=O when Pc=O
&
=o.
dv
(85)
.%=O
As differentiations w.r.t. r and w.r.t. p or v do not commute, but differentiations w.r.t. r and w.r.t. 9 do, it is easier to obtain dA/dp and dA/dv by differentiating (84) w.r.t. 9 directly:
dJ4 --dP
dp
--+x(2gv)+~vxv
- dPt”
+
;v(2‘&)xv+
fvx(2g-&J,
1 (86)
hence
%I,_,= $+d$J =-41m(Ubj)
fVx[2gVx(4>].
(87)
Therefore, from ( 69 ) Ii)=-4iIm(Ubj)f?
-4ej
2g(Vx$)-(Vx$)dr],
(91)
enabling Ubjto be found, and hence the magnetisability to be evaluated using eq. (32) as all the other quantities are known. Although more work needs to be done before we can present any results from calculations on specific molecules, all the quantities in eq. (9 1) except g are calculated as a matter of routine in a standard DFT code. g is given in numerical form by Vignale et al. [ 41, and needs to be fitted to a simple function, but we do not envisage any difficulties with this. We believe that the theory as presented here is complete, and hence that it is appropriate to present it at this stage. Clearly there are questions about the numerical accuracy of this method in practice, and we can make no predictions about this yet. Our key interest is to extend this work to non-linear optics to derive working expressions for the determination of frequency-dependent polarisabilities. Now that we have obtained the form of the magnetic Hessian we anticipate that this will be more straightforward.
(a;;;TVx,2gVx(&$‘),dr References =-4iIm(LTb,)ej2g(VxF)-(VxF)dr.
(88) Hence now Hlai,bj=(~a-Ei)~=i,,+4(ailbj)
(89)
-4e/ 278
2g(Vx$)-(Vx%)dr,
(90)
[ 1 ] SM. Colwell, C.M. Murray, N.C. Handy and R.D. Amos, Chem. Phys. Letters 210 (1993) 261. [ 2 ] G. Senatore and K.R. Subbaswamy, Phys. Rev. A 34 ( 1986) 3619. [ 31 V.G. Malkin, O.L. Malkina and D.R. Salahub, Chem. Phys. Letters 204 (1993) 87. [4] G. Vignale, M. Rasolt and D.J.W. Geldart, Advan. Quantum Chem. 21 (1990) 235; G. Vignale and M. Rasolt, Phys. Rev. Letters 59 ( 1987) 2360; Phys. Rev. B 37 (1988) 10685. [ 51R.M. Stevens, R.M. Pitzer and W.N. Lipscomb, J. Chem. Phys. 38 (1963) 550. [ 61 P. Hohenberg and W. Kohn, Phys. Rev. 136 ( 1964) 864. [7] W. Kohn and L.J. Sham, Phys. Rev. A 140 (1965) 1133.
CHEMICAL PHYSICS LETTERS
14 January 1994
The determination of magnetisabilities using density functional theory Susan M. Colwell and Nicholas C. Handy University Chemical Laboratory, LensJeld Road, Cambridge CB2 IE W, UK Received 24 June 1993; in final form 23 November 1993
The theory for the calculation of magnestisabilities using current density functional theory is presented, following on from the original theory of Vignale, Rasolt and Geldart. It is shown that a Hellmann-Feynman identity holds for the first derivative of the energy with respect to a magnetic field. The magnetic Hessian which arises in the coupled perturbed Kohn-Sham equations contains an extra term involving the current density operator as well as the density itself. It is now in a form for computation.
1. Introduction The application of density functional methods to molecular problems has made significant advances recently with the development of new functionals, and it is a field that is growing rapidly. Calculations of energies, structures and vibrational frequencies are well developed, and in general are more accurate than Hartree-Fock self consistent field (SCF) calculations at similar cost. SCF calculations of electric and magnetic properties are known to be inaccurate, and correlated methods, particularly for frequency-dependent properties, are complicated and expensive. The result is that only a limited number of such calculations have been performed, and those only on small molecules. It is therefore important to consider density functional theory (DFT) methods for obtaining these properties. The authors have recently presented DFT calculations of polarisabilities, cr, and hyperpolarisabilities, /3, of molecules [ 1 ] using coupled Kohn-Sham theory and Senatore and Subbaswamy [ 21 have reported calculations of cr and second hyperpolarisabilities, y, for atoms. Many other workers have used finite field DFT calculations to obtain electric properties of molecules. Although there have been calculations of magnetic properties using DFT (e.g. Salahub et al. [ 3 ] ) these have not been based on a full coupled Kohn-Sham treatment. The difficulty, which is widely acknowledged, is that
with standard DFT using an exchange-correlation functional which only depends on the density, the corresponding exchange-correlation part of the magnetic Hessian vanishes. Therefore it is necessary to use a density functional theory which introduces the current density and specifically introduces an exchange-correlation functional which depends on both densities. The basic theory for this current density functional theory (CDFT) has been presented by Vignale, Rasolt and Geldart [ 41, and they have given the specific form of the Kohn-Sham equations. In this Letter, we discuss the evaluation of magnetisabilities. It is necessary to derive the appropriate coupled perturbed Kohn-Sham (CPKS) theory for current density functional theory. We present an explicit expression for the magnetic Hessian within this theory using the local density approximation. We find that a Hellmann-Feynman identity holds, and hence the expressions for second-order properties are straightforward.
2. The Schriidinger equation in an external electromagnetic field In this section we present dard electromagnetic theory so that subsequent results can context. The non-relativistic
0009-26141941% 07.00 0 1994 Elsevier Science B.V. All rights reserved. SSDZ OOOOS-2614(93)El399-2
an exposition of stanfor completeness, and more easily be put into Hamiltonian for a sys271
14 January 1994
CHEMICAL PHYSICS LETTERS
Volume 2 17, number 3
tern of N electrons and Nz nuclei moving in an external electromagnetic field (E, B) is given (within the Born-Oppenheimer approximation) by
where 3 is the Fock operator. This leads to the definition of the Fock matrix, which is diagonal at convergence if canonical orbitals are used. (6)
~,,=(PIhIq)+2(pqIjj)-(~jI~q)=~P~,,
(1)
where veXt= -e@, and @is the scalar electrostatic potential such that E= -V@-dA/dt, and A is the magnetic vector potential such that B=VxA. For a uniform field VeXt=eE.r and, in the Coulomb gauge
where we use the standard repeated suffix notation, with i, j, ... being occupied orbitals, a, b... being virtual orbitals and p, q, ... being any molecular orbitals. The expression for the RHF energy is given by fT=2(ilhli)+2(iiljj)-(ijlji).
(7)
In the presence of an electromagnetic field, the energy expression is still as above, but now
(V.A=O), A=bBxr. We define the density and current operators acting on the ground state wavefunction P(rr, r2, .... rN) by j?= U*(r) Y(r)
=N
J
!P(r, r,, .... rN) r2, .... rN) dr, drs...drN ,
x Y(r, f 3p=-
2
[
!P(r)V!P(r)-V!P(r)Y(r)]
(2)
(8)
=ho@+h@, where the V-A as an operator, the derivatives external field series,
is not to be interpreted as div A, but i.e. the term is V. (A@). Now consider of the energy with respect to a real 9’. Express the energy as a Taylor
(3)
(note that our definitions do not include any factors of e), i.e. the dependence of !Pon all but one spatial variable has been integrated out. The ground state energy is then given by
J
p(r)
J=$+ +e
J
Write the derivative
(;,4*+V,,>dr
jp,(r)-Adr ,
where 8. is the ground of an external field.
of a molecular
orbital as
(10) (4) or alternatively
state energic in the absence
@p=@; + @u;“,@: + $%YJ$@
+...
(11)
so that uL;,=@pl;9=o.
3. Coupled-perturbedHartree-Fock theory for electric and magnetic properties Consider closed-shell restricted Hartree-Fock (RHF ) and a system of iV= 2M electrons. The wavefunction is a single Slater determinant of M doubly occupied orbitals &, which are obtained as the M lowest eigenfunctions of the eigenvalue equation 9@i=ci@i, 272
(5)
(12)
Orthonormality of orbitals implies that 0 is skew Hermitian. Also invariance with rotations within the occupied and virtual spaces mean that without loss of generality we can chose &&=o.
(13)
Differentiation of the energy with respect to the field P and using eq. (6) yields the usual HellmannFeynman identity,
f$
=2(il$
Ii) .
d.+@ /=o
=2(il
dh df”
LO I‘) 1 ,
(15)
where I i) now refers to @p,i.e. the unperturbed orbitals. Second-order properties are given by d26 dp’d9
9 ‘=O
)
P=O
(22) (23)
Im ( U$)H2ai,bj = 0
and so U is purely real. The energy is expressed as - $EAE,~A, +... ,
#u;= -2Pi
Ia)
(24)
(25)
is the permanent dipole moment, and from ( 16) (YAP = - 4 U$ P$
+2iIm(U{i) dh
U&)H,+j=O
Pti+Re(
where
dh IO +(‘I -I d.9” ;9=o >
dFA 9=. -I
where P’;, is the dipole integral. The CHF equations give
6= 8. -E&
d2h Ii) 9=0 =2(i’ dY”dF6’ 99=.
+2Re(U$)
-(‘I
(21)
(14)
First-order properties are therefore obtained from db
14 January 1994
CHEMICAL. PHYSICS LETTERS
Volume 2 17. number 3
la)
IO . >
(16)
To evaluate these properties we need the U matrices, which are obtained from the coupled Hartree-Fock equations. These equations are obtained by differentiating the condition Fai=O 3
with respect to Ta and evaluating at 9=0
(17) to give
is the polarisability. For a uniform magnetic field, the first-order perturbation to the Hamiltonian is pure imaginary (- (i!ze/2m) (A-V+V*A), to be interpreted as before), and YA=BA, (rl g
IB=ols)
(@,I~~B~o~v+v~~~B~oldJs,
=-g i#ie
=- G
I
(18)
-where Hlai,bjand Hzai,bjare the electric and magnetic Hessians, respectively, as first presented by Stevens, Pitzer and Lipscomb [ 5 1, and are given by H,,i,aj=4(aiIbj)-(ajlbi)-(abIji) + (e=-+5ai,bj, Hzai,bj=(aj(bi)-(~blii)+(~“-t~)~=ni,bj
V@), dr=ihr;l,,
2@(rx
(27)
where iN = eL, and N is real and skew Hermitian. The CHF equations give Re(U$)H1,,,=O
+ i Im ( U&)Hzni,bj= 0 ,
(26)
(28)
3
Nb+Im(U&)Hm,bj=O
3
(29)
and so U is purely imaginary. The energy is expressed as ~=$-E~ml-fEnE,~~r+...,
(30)
where (19) (20)
for perturbations out of real orbitals. For a uniform electric field, the perturbation to the Hamiltonian ( =eE*r) is real, 9@=EA,
m,=-2iNi=0
(31)
is the magnetic moment (which vanishes within RHF), and &- --2Qf/‘++Im(U$)N$
(32) 273
Volume 2 17, number
3
is the magnetisability,
CHEMICAL
PHYSICS
14 January
LETTERS
where
1 e2 p(r,Mr2)
J[p] = -4ne, 2 s is the Coulomb energy is
These expressions are equivalent Stevens et al. [ 5 1.
Colwell et al. [ 1 ] and others have implemented coupled Kohn-Sham theory (which is exactly analogous to the above but using Kohn-Sham theory rather than RHF) and successfully calculated static polarisabilities and hyperpolarisabilities. However, a straightforward application of ordinary density functional theory to magnetic properties does not work because the exchange correlation part of the magnetic Hessian vanishes. Partly in order to overcome this kind of problem, Vignale, Rasolt and Geldart [ 41 have introduced current density functional theory where the exchange-correlation energy is a functional, not only of the density p(r) but also of the currentj, (r). They have also formulated the theory in terms of the gauge invariant vector y(r) , where
Ij,(r)/Ar)l
,
(34)
and we follow their theory in this section, before applying it in the next section to get new formulae for electric and magnetic properties. In the presence of an electromagnetic field the current density functional theory expression for the energy of the system previously.discussed (4) becomes
(35)
(37) 2
energy. The non-interacting
kinetic
where x’[p, j,] is a single Slater determinant of the M doubly occupied non-interacting orbitals which are the lowest eigenfunctions of the appropriate orbital equation. The generalised Hohenberg-Kohn theorem states that gvA is minimised when p and j, take the correct values for the potentials A and VeXt, and this minimum value is the true ground state energy for the system. The appropriate orbital equation for the non-interacting orbitals must have the form
(
&
where the current P=
V&r)
( -ifiV+eA,e)2+
>
@i=Ei@i,
(39)
and density are now given by
ye; 3
(41 and the effective potentials found. The energy is
+
Aeff and
Vcrr are to be
I w(r) dr+ I pV,,(r) dr+&W,l
which can be written
+ &
s
+&,[pJ,l 3
where u(r) is the nuclear 274
dr
’
,
(42)
as
(A’-A&)pdr+e2
j, * (A -4d
dr
+J[PI + s pv(r) dr+ s pK,,(r> dr
where
+Ecc[P,jpl
dr
to those given by
4. Current density functional theory
v(r)=Vx
Ic-r2l
1994
attraction
potential,
.
(43)
(36)
Taking the functional derivative of (43) with respect to p, (6&,/6p) lip, gives rise to the definition
and
of v,ff,
Volume 2 17, number 3
CHEMICAL PHYSICS LETTERS
I4 January 1994
be substituted back into the orbital equation (39) and the energy expression (43). + T/H(r) + ~x’x,(r)3
(44)
and the functional derivative with respect to j,,
(Wc/6jp)[P,j,l Ipgives 4ff(r) =A(r) +4,(r) ,
(45)
where
5. Energy derivatives in current Kohn-Sham
We now rewrite this in standard Kohn-Sham [ 6,7 ] language. The standard ground state energy expression is &=2(iIhIi)+2(iiIjj)+Ex,[p],
(46) (47)
theory
(53)
where the usual exchange term of HF theory is replaced by the exchange-correlation energy functional. Within current density functional theory this becomes, from eq. (43 ) , 8=2(iIh,,Ii)+2(iiIjj)+Ex,,[p,
v]
(48) (54)
The orbital equation becomes
where we define
- ~V2+v(r)+VH(r)+Vs(r) +
e*
~A*+Kx,-
ifie 2m
(A-V+ V-A)
>
+i=~igi,
- E
(49)
and, given E,, [p, j,] , this can in principle be solved to give the orbitals, and hence the density, current and energy. Vignale, Rasolt and Geldart [ 41 now recast the problem in terms of the gauge invariant vector y(r)=Vx v,(r)/p(r)]. It can be shown that Exc[~Jpl =E&,
~1 ,
(50)
i.e. the exchange-correlation energy is actually a functional of p and v and so is gauge invariant. Taking the functional derivatives of (50) with respect to p and v gives new expressions for V,,(r) and A,,(r), SE
V,,(r) = $
[p,v]l,-e’*,
(A,,~V+V*A,,)
=h+h,,
= i Vx $$
[P,Vllp,
(55)
The Kohn-Sham matrix is now
(56)
In the rest of this section we consider the derivatives of the energy with respect to a real external field .@, and show that a Hellmann-Feynman identity holds. We follow a similar procedure to the RHF case (section 2) with the same definition of 0 etc. Working within the LDA approximation, i.e. assuming there are no Vp or V-j, dependent terms in the potential, dE m =2(il $Fj li)+2(il
s
Ii)
(51) +2&[FY-(ilV,,Ir)]+
eA,,(r)
.
(52) -e
giving a gauge invariant expression for V, which can
s
dlip(r)~&l df”
s$$dr dr ’
(57)
275
Volume 2 17, number 3
CHEMICAL PHYSICS LETTERS
where
(58)
substituting Hermitian, CL!?
m
Note
(59) Then
in (57) and remembering one obtains
8 is skew
Ii).
(65)
Note that this is dh/d@ not dh,,/dP, i.e. we have verified that the Hellmann-Feynman identity holds for current density functional theory, as we would expect from the Kohn-Sham condition Fy ~0. Therefore first- and second-order electric and magnetic properties are given by equations which are formally the same as in the RHF case ( 15) and ( 16), but I i) now refers to @p, i.e. the unperturbed KohnSham orbitals. To evaluate these we need the U matrices, which are now obtained from the coupled-perturbed Kohn-Sham equations, described in the next section.
6. Coupled-perturbed current Kohn-Sham
but d& dp
=2(il$
14 January 1994
(61)
The coupled-perturbed Kohn-Sham equations obtained by differentiating the condition Fy=O
and
(62)
theory are
(66)
with respect to FA and evaluating entiating (56) gives
at 9~0.
Differ-
so
(63)
Then
dV
+@I xc dp”
14) )
and it can be seen that
276
(67)
14 January 1994
CHEMICAL PHYSICS LETTERS
Volume 2 17, number 3
choosing oij= oa:ob=0 as before
dl/,,
(69)
(” d3@ 6=. IO
where (77)
For perturbations out of real orbitals j, is first order in 4 and so is A,, and so the last term in (68) can be ignored when evaluating at T= 0. Also, from (62) (71)
dhxc _ li)=4iRe(&j)e
(al s
9-O x
k
d(Axca) =lp
X
dp
@bbQj
dr-4iIm(
obj>e
9F=o
J
k~ia~l,=,E~&~$).
(78)
So the coupled perturbed Kohn-Sham equations have the same form as the CPHF equations ( 18 ), but now Hlai,bj and Hzoi,bj, the electric and magnetic Hessians respectively, are given by
(72)
hence +4e
JjaiylF=?ejdr
(79)
(73)
(74)
If the orbitals are initially real, jj,: ikj, is pure imaginary, i.e. kjr= -kv is real. Hence di, dP
p=o
= -4Im(
Od)kj,
(75)
,??tF*[p, v] =EkFA[p]o+
and therefore dv, =-41m(Uti)VX d@ y=o
where the j matrices have been reintroduced in place of ik. Vignale et al. [ 41 give an expression for the LDA exchange-correlation energy in the case of a weak field (which may be the applicable case). It appears that it can be written
kjr 0
p
,
(76)
J
F(p”‘)
1v(r) I’dr
and sufficient data is given to enable F(p’/‘)
(81) to be 277
Volume 2 17, number 3
CHEMICAL PHYSICS LETTERS
fitted. Specifically for this form of the functional, setting F(Y)
=g(P) 3
(82)
(83)
E:,,(r) = ~xo(r> +&T’(P) I y(r) 12,
-L(r)
= j Vx
KWpPl
(84)
14 January 1994
where the j matrices have been reintroduced in place of ik. HI and Hz are now real, and so the same deductions can be made as in the RHF case. So for purely electric fields U is real and one recovers the usual CKS equations. For purely magnetic fields U is pure imaginary, and one obtains
hence, as v=O when Pc=O
&
=o.
dv
(85)
.%=O
As differentiations w.r.t. r and w.r.t. p or v do not commute, but differentiations w.r.t. r and w.r.t. 9 do, it is easier to obtain dA/dp and dA/dv by differentiating (84) w.r.t. 9 directly:
dJ4 --dP
dp
--+x(2gv)+~vxv
- dPt”
+
;v(2‘&)xv+
fvx(2g-&J,
1 (86)
hence
%I,_,= $+d$J =-41m(Ubj)
fVx[2gVx(4>].
(87)
Therefore, from ( 69 ) Ii)=-4iIm(Ubj)f?
-4ej
2g(Vx$)-(Vx$)dr],
(91)
enabling Ubjto be found, and hence the magnetisability to be evaluated using eq. (32) as all the other quantities are known. Although more work needs to be done before we can present any results from calculations on specific molecules, all the quantities in eq. (9 1) except g are calculated as a matter of routine in a standard DFT code. g is given in numerical form by Vignale et al. [ 41, and needs to be fitted to a simple function, but we do not envisage any difficulties with this. We believe that the theory as presented here is complete, and hence that it is appropriate to present it at this stage. Clearly there are questions about the numerical accuracy of this method in practice, and we can make no predictions about this yet. Our key interest is to extend this work to non-linear optics to derive working expressions for the determination of frequency-dependent polarisabilities. Now that we have obtained the form of the magnetic Hessian we anticipate that this will be more straightforward.
(a;;;TVx,2gVx(&$‘),dr References =-4iIm(LTb,)ej2g(VxF)-(VxF)dr.
(88) Hence now Hlai,bj=(~a-Ei)~=i,,+4(ailbj)
(89)
-4e/ 278
2g(Vx$)-(Vx%)dr,
(90)
[ 1 ] SM. Colwell, C.M. Murray, N.C. Handy and R.D. Amos, Chem. Phys. Letters 210 (1993) 261. [ 2 ] G. Senatore and K.R. Subbaswamy, Phys. Rev. A 34 ( 1986) 3619. [ 31 V.G. Malkin, O.L. Malkina and D.R. Salahub, Chem. Phys. Letters 204 (1993) 87. [4] G. Vignale, M. Rasolt and D.J.W. Geldart, Advan. Quantum Chem. 21 (1990) 235; G. Vignale and M. Rasolt, Phys. Rev. Letters 59 ( 1987) 2360; Phys. Rev. B 37 (1988) 10685. [ 51R.M. Stevens, R.M. Pitzer and W.N. Lipscomb, J. Chem. Phys. 38 (1963) 550. [ 61 P. Hohenberg and W. Kohn, Phys. Rev. 136 ( 1964) 864. [7] W. Kohn and L.J. Sham, Phys. Rev. A 140 (1965) 1133.