Geometrical derivatives and magnetic properties in atomic-orbital density-based Hartree–Fock theory

Geometrical derivatives and magnetic properties in atomic-orbital density-based Hartree–Fock theory

Helena Larsen^a) and Trygve Helgaker

Department of Chemistry, University of Oslo, P.O. Box 1033 Blindern, N-0315 Oslo, Norway Jeppe Olsen and Poul Jørgensen

Department of Chemistry, University of Aarhus, Langelandsgade 140, DK-8000 A˚ rhus C, Denmark.

共Received 18 July 2001; accepted 11 September 2001兲

A reformulation of Hartree–Fock theory for time-independent molecular properties with perturbation-dependent basis sets and which refers strictly to the atomic-orbital basis is presented.

The formalism is based on a recently proposed exponential parametrization of the one-electron atomic-orbital density matrix. In the presented formulation, only multiplications and additions of sparse matrices are needed. Linear scaling with system size is therefore obtainable, making this formulation ideally suited to large molecular systems. The paper contains general formulas for molecular energy derivatives up to fourth order, with special attention given to molecular gradients, molecular Hessians, magnetizabilities, and nuclear magnetic shieldings. © 2001 American Institute of Physics. 关DOI: 10.1063/1.1415082兴

I. INTRODUCTION

A major challenge of quantum chemistry is to reduce the scaling of the standard electronic-structure methods so as to treat large molecules. Ideally, for a sufficiently large system, the cost of the calculation should scale linearly with the size of the system. With conventional techniques, systems con- taining more than a few hundred atoms are beyond the pow- ers of present ab initio electronic-structure models. However, in order for quantum chemistry to contribute to many areas of biological and medical interest, it must be able to handle much larger systems such as proteins and nucleic acids. In recent years, therefore, much effort has been put into the development of such techniques.

In particular, within the Hartree–Fock and Kohn–Sham self-consistent-field theories, which formally scale as the fourth power in the size of the system, significant progress has been made. While these methods are important in them- selves in that they provide good compromises between cost and accuracy, the achievement of linear scaling for these methods is also a necessary first step towards the develop- ment of linear scaling for the more complex, correlated ab initio methods such as Møller–Plesset and coupled-cluster theories.

By using the fast multipole method共FMM兲or the tree- code method combined with prescreening, the complexity of the long-range Coulomb problem in Hartree–Fock and Kohn–Sham theories 共i.e., the construction of the Fock and Kohn–Sham matrices兲has been reduced to linear.^{1– 6}How- ever, another bottleneck in these theories is the optimization of the energy, which is traditionally achieved by diagonaliza- tion of the Fock/Kohn–Sham matrices, a step that scales cu- bically with system size.

To achieve linear scaling, several alternatives to diago- nalization have been suggested, based either on the use of localized orbitals^7–9or on the direct optimization of the den- sity matrix in the atomic-orbital 共AO兲 basis, using the ele- ments of the density matrix as variational parameters.^10–14In the latter approach, the sparsity of the Fock/Kohn–Sham, overlap, and density matrices is exploited to obtain linear scaling. Thus, it has now become possible to compute ener- gies and to optimize geometries for systems containing thou- sands of atoms using minimal basis sets. These ideas are also applicable for semiempirical methods. Some studies involv- ing correlated methods also exist.^15–18For molecular proper- ties, only a few studies have appeared, including studies of force constants and forces.¹⁴

In Refs. 19 and 20, we proposed an exponential param- etrization of the AO density matrix, discussing how it can be used to optimize the density matrix and thus the energy. In Ref. 21, this parametrization was used to develop a time- dependent Hartree–Fock/Kohn–Sham response theory that refers solely to the AO basis and thus is suitable for linear scaling. We here extend these studies to molecular properties where the AOs depend on the perturbation. While general expressions are given for derivatives to fourth order, molecu- lar gradients, molecular Hessians, magnetizabilities, and nuclear magnetic shieldings receive special attention and have 共except for the shieldings兲been implemented in a de- velopment version of the Dalton program. As will become apparent from this paper, the density-based derivative theory is simpler than the corresponding molecular-orbital 共MO兲 based theory, involving a smaller set of AO matrices and mathematical operations.

The calculation of derivatives of the Hartree–Fock en- ergy for perturbation-dependent basis sets has a long history.

The expression for the Hartree–Fock gradient was derived by Bratozˇ in 1958²² and the first practical implementation was reported by Pulay in 1969.²³ Molecular Hessians were

a兲Permanent address: Department of Chemistry, University of Aarhus, Langelandsgade 140, DK-8000 A˚ rhus C, Denmark.

10344

0021-9606/2001/115(22)/10344/9/$18.00 © 2001 American Institute of Physics

first implemented by Thomsen and Swanstrøm in 1973,²⁴ based on the work of Bratozˇ²² and Gerratt and Mills.²⁵The first practical implementation was reported by Pople et al. in 1979.²⁶Third derivative expressions were derived by Moccia in 1970²⁷and implemented by Gaw and Handy in 1984.²⁸

The above formulations and implementations of deriva- tives were all given in the MO basis. In 1977, Dodds, McWeeny, and Sadlej presented a density-matrix formulation of first and second derivatives in the AO basis.²⁹ However, their formulation differs from ours in that it requires the di- agonalization of the Fock matrix to determine the perturbed density matrix. In our formulation, no such diagonalization is required—the derivatives are obtained directly from the共dif- ferentiated兲Hamiltonian and overlap integrals in the AO ba- sis and from the 共differentiated兲 density matrix in the AO basis. As such, the present formulation is closer to the recent work of Ochsenfeld and Head-Gordon,¹⁴ who have derived the molecular gradient and Hessian in a density-based for- mulation of Hartree–Fock theory. The main difference be- tween our derivation and that of Ochsenfeld and Head- Gordon is that we use an explicit parametrization of the AO density matrix, where redundancies have been identified and removed. Without the proper removal of redundancies, the solution of the response equations and thus the calculation of second- and higher-order molecular properties may become difficult or even impossible for large systems.

After a discussion of the Hartree–Fock energy and the parametrization of the density matrix in Sec. II, we consider general expressions for the derivatives to fourth order in Sec.

III. The derivatives of the AO density matrix are then dis- cussed in Sec. IV. In Sec. V, we consider in greater detail the evaluation of a few selected properties: the molecular gradi- ent, the molecular Hessian, and the magnetizability and nuclear shielding constants of closed-shell molecules.

II. PARAMETRIZATION OF THE HARTREE–FOCK ENERGY

In the present section, we review the density-based for- mulation of the Hartree–Fock energy. After a brief discus- sion of the Hartree–Fock energy in Sec. II A, we consider the exponential parametrization of the AO density matrix in Sec. II B and the Hartree–Fock variational conditions in Sec.

II C.

A. The Hartree–Fock energy

In the atomic spin–orbital basis, the Hartree–Fock en- ergy may be written in the form

E⫽Tr Dh⫹ ¹2Tr DG共D兲⫹h_nuc. 共1兲 Here D is the one-electron density matrix, h is the one- electron Hamiltonian matrix with elements

h_␮␯⫽

冕

冉

^兺

冊

and G(D) contains the matrix elements

G_␮␯共D兲⫽

兺

where

g_␮␯␳␴⫽

冕 冕

dx₁dx₂. 共4兲 Both spin and space coordinates are included in x. The last term in Eq.共1兲contains the nuclear–nuclear contribution to the total energy; for the calculation of magnetic properties, it also contains purely nuclear magnetic interactions, such as the Zeeman interaction between the nuclear magnetic dipole moments and the external magnetic field.

For Eq.共1兲to represent a valid Hartree–Fock energy, the AO density matrix must satisfy certain conditions. In particu- lar, in the spin–orbital representation, the density matrix must fulfill the symmetry, rank, and idempotency conditions:

D^†⫽D, 共5兲

Tr DS⫽N_e, 共6兲

DSD⫽D, 共7兲

where N_e is the number of electrons. When the Hartree–

Fock energy is optimized or perturbed in any manner, the modifications to the density matrix must comply with these conditions. Although this may be achieved by a constrained optimization of the energy, it is more convenient to work with an unconstrained, explicit parametrization of the AO density matrix. Such a parametrization was presented in Refs. 19 and 20.

B. Parametrization of the AO density matrix

Let D be any matrix that satisfies the symmetry, rank, and idempotency conditions Eqs. 共5兲–共7兲. From this refer- ence density matrix, we may generate any valid Hartree–

Fock density matrix by carrying out the transformation D共X兲⫽exp共⫺XS兲D exp共SX兲. 共8兲 Here S is the AO overlap matrix

S_␮␯⫽

冕

while X is an anti-Hermitian matrix, which may be decom- posed into an antisymmetric real part X^R and a symmetric imaginary part X^I

X⫽X^R⫹iX^I. 共10兲

The need to introduce imaginary rotations arises since we shall here consider not only real perturbations 共nuclear dis- placements兲, but also perturbations associated with imagi- nary operators 共magnetic perturbations兲.

Since D(X) is a valid Hartree–Fock density matrix pro- vided X is anti-Hermitian, we may now optimize the Hartree–Fock energy in an unconstrained manner by varying the independent elements of X^R and X^I. However, it should be noted that, although the parametrization Eq.共8兲is uncon- strained, it is also redundant in the sense that not all nonzero anti-Hermitian matrices X generate a modified density ma- trix. To avoid such redundant rotations

D共X_red兲⫽D共0兲⫽D, 共11兲

the anti-Hermitian matrix X should comply with the projec- tion relation

X⫽P共X兲⬅PXQ^†⫹QXP^†, 共12兲 where

P⫽DS, 共13兲

Q⫽1⫺DS, 共14兲

are nonorthogonal projectors onto the occupied and virtual orbital spaces, respectively.

C. Variational conditions

The transformed AO density matrix D(X) may be evalu- ated using the asymmetric Baker–Campbell–Hausdorff 共BCH兲expansion¹⁹

D共X兲⫽D⫹关D,X兴S⫹¹2关关D,X兴S,D兴S⫹•••, 共15兲 where we have introduced the S commutator

关D,X兴S⫽DSX⫺XSD. 共16兲

Expanding the Hartree–Fock energy Eq.共1兲to first order in X, we obtain the variational condition

␦^E⫽Tr F␦^D共X兲⫽Tr F关D,␦^X兴S

⫽Tr共FDS⫺SDF兲␦^X

⫽0, 共17兲

where we have introduced the Fock matrix

F⫽h⫹G共D兲. 共18兲

Since Eq.共17兲should hold for arbitrary X, we may write the variational conditions in the form

FDS⫽SDF. 共19兲

Together with the symmetry, rank, and idempotency condi- tions Eqs. 共5兲–共7兲, this equation uniquely determines the Hartree–Fock AO density matrix and thus the optimized Hartree–Fock energy given in Eq.共1兲.

III. DERIVATIVES OF THE HARTREE–FOCK ENERGY Having considered the representation of the Hartree–

Fock energy and density matrix in the AO basis, let us now turn our attention to the derivatives of the Hartree–Fock en- ergy E(x) with respect to a set of external parameters x.

共Previously, we have used x for the combined spatial and spin electronic coordinates; in the following, x will be used only for the external parameters.兲 For much of our discus- sion, the nature of these parameters is unimportant. However, in the final sections of this paper, we shall consider in greater detail two special cases: geometrical distortions and external magnetic fields. These perturbations are of special interest since, for their accurate description, both require the use of perturbation-dependent basis sets: atom-fixed AOs for geo- metrical distortions and field-dependent London AOs for ex- ternal magnetic fields. The atom-fixed London AOs may be written in the general form

␹␮共r_M兲⫽exp关⫺¹2iB⫻共M⫺O兲•r兴S_lm共r_M兲exp共⫺␣^rM 2 兲, 共20兲 where r_M⫽r⫺M is the position of the electron relative to the AO center M; O is the origin of the vector potential A

⫽ ¹2B⫻(r⫺O) chosen to represent B in the Hamiltonian;

␣⬎0 is the orbital exponent; and S_lm(r_M) is a real solid- harmonic function. Clearly, these orbitals depend both on the molecular geometry and on the external magnetic field. The use of such perturbation-dependent AOs imply that the over- lap matrix S changes as we apply the perturbation. As we shall see in Sec. IV, this fact has important consequences for the evaluation of derivatives.

By straightforward differentiation of Eq.共1兲with respect to x, we obtain the following expressions for the derivatives of the Hartree–Fock energy to fourth order:

E^a⫽Tr Dh^a⫹¹2Tr DG^a共D兲⫹Tr D^aF⫹h_nuc^a , 共21兲 E^ab⫽Tr Dh^ab⫹ ¹2TrDG^ab共D兲⫹Tr D^abF

⫹2 P_abTr D^a关h^b⫹G^b共D兲⫹ ¹2G共D^b兲兴⫹h_nuc^ab , 共22兲

E^abc⫽Tr Dh^abc⫹ ¹2Tr DG^abc共D兲⫹Tr D^abcF

⫹3 P_abcTr D^a关h^bc⫹G^bc共D兲⫹G^b共D^c兲兴

⫹3 P_abcTr D^ab关h^c⫹G^c共D兲⫹G共D^c兲兴⫹h_nuc^abc, 共23兲

E^abcd⫽Tr Dh^abcd⫹ ¹2Tr DG^abcd共D兲⫹Tr D^abcdF

⫹4 P_abcdTr D^a关h^bcd⫹G^bcd共D兲⫹³2G^bc共D^d兲兴

⫹6 P_abcdTr D^ab关h^cd⫹G^cd共D兲⫹2G^c共D^d兲

⫹¹2G共D^cd兲兴⫹4 P_abcdTr D^abc关h^d⫹G^d共D兲⫹G共D^d兲兴

⫹h_nuc^abcd. 共24兲 We have here used the short-hand notation

f^a⫽d f共x兲

dx_a 共25兲

for the first derivative and likewise for the higher derivatives.

Also, we have for brevity introduced the symmetrization op- erator

P_abA^aB^b⫽¹2共A^aB^b⫹A^bB^a兲 共26兲 and similarly for P_abc and P_abcd. Finally, in deriving Eqs.

共21兲–共24兲, we have made frequent use of the identity

Tr AG共B兲⫽Tr BG共A兲, 共27兲

which follows from the particle symmetry of the two- electron integrals.

Equations共21兲–共24兲are general, expressing the deriva- tives of the Hartree–Fock energy in terms of derivatives of the AO integrals and the AO density matrix. For higher de- rivatives, a large number of terms arise, making the calcula- tion of such derivatives a daunting task. We note, however, that many of the derivative共Hamiltonian and density兲matri-

ces that contribute to a given order n in the differentiation also contribute to derivatives of orders k⬍n, thereby reduc- ing the overall cost of calculating all derivatives up to order n.

The derivatives of the one- and two-electron Hamil- tonian AO integrals that appear in Eqs. 共21兲–共24兲 may be calculated using any of the standard techniques of molecular integral evaluation, see, for example, Chap. 9 of Ref. 19. In the following, we shall assume that the AO integrals may be generated as needed and instead concentrate our attention on the calculation of the derivatives of the AO density matrix in Eqs. 共21兲–共24兲.

IV. DERIVATIVES OF THE AO DENSITY MATRIX Having considered the general structure of the Hartree–

Fock energy derivatives in Sec. III, we now turn our atten- tion to the derivatives of the AO density matrix

D共X兲⫽exp共⫺XS兲D₀exp共SX兲, 共28兲 where D₀ is the optimized AO density matrix of the unper- turbed system. Before taking derivatives, we note that all three matrices S, D₀, and X on the right-hand side of Eq.

共28兲are affected by the perturbation x:

共1兲 S changes because the AOs in Eq. 共20兲 depend on the perturbation.

共2兲 D₀changes to comply with the idempotency condition in Eq. 共6兲, which is affected by the perturbed overlap inte- gral.

共3兲 X changes to satisfy the Hartree–Fock condition in Eq.

共19兲, which is affected by the perturbed Hamiltonian and overlap integrals.

For perturbations that do not affect the AOs, S is un- changed and D₀ remains valid also after the perturbation. In such cases, only X changes as we perturb the system. In

general, however, S, D₀, and X all change with the pertur- bation, generating a large number of contributions to high- order derivatives of the AO density matrix Eq.共28兲, making it important to develop a good strategy for their evaluation.

The remainder of this section consists of three parts. In Sec. IV A, we discuss how to assemble the derivatives of the density matrix D(X) from the derivatives of D₀, S, and X, assuming that these are already available. Next, in Sec. IV B, we consider the evaluation of the derivatives of D₀ from the derivatives of the idempotency condition Eq. 共6兲. Finally, in Sec. IV C, we calculate the derivatives of X from the differ- entiated Hartree–Fock conditions Eq. 共19兲.

A. Derivatives of the BCH expansion

Carrying out the asymmetric BCH expansion of the ref- erence density matrix Eq. 共15兲 and differentiating at X⫽0, we obtain to fourth order:

D^a⫽D₀^a⫹关D₀,X兴S

a, 共29兲

D^ab⫽D₀^ab⫹关D₀,X兴S

ab⫹¹2关关D₀,X兴S,X兴S

ab, 共30兲

D^abc⫽D₀^abc⫹关D₀,X兴S

abc⫹ ¹2关关D₀,X兴S,X兴S abc

⫹¹6关关关D₀,X兴S,X兴S,X兴S

abc, 共31兲

D^abcd⫽D₀^abcd⫹关D₀,X兴S

abcd⫹¹2关关D₀,X兴S,X兴S abcd

⫹ ¹6关关关D₀,X兴S,X兴S,X兴S abcd

⫹ 24¹ 关关关关D₀,X兴S,X兴S, X兴S,X兴S

abcd. 共32兲 Since each S commutator is a linear combination of products of three matrices Eq. 共16兲, the differentiated nested commu- tators rapidly become complicated, with a large number of contributing terms. For completeness, we here list the full expressions for the derivatives to fourth order:

D^a⫽D₀^a⫹关D₀,X^a兴S, 共33兲

D^ab⫽P_ab共D₀^ab⫹关D₀,X^ab兴S⫹2关D₀^a,X^b兴S⫹2关D₀,X^a兴S^b⫹关关D₀,X^a兴S,X^b兴S兲, 共34兲

D^abc⫽P_abc共D₀^abc⫹关D₀,X^abc兴S⫹3关D₀^a,X^bc兴S⫹3关D₀,X^ab兴S^c⫹3关D₀^ab,X^c兴S⫹6关D₀^a,X^b兴S^c⫹3关D₀,X^a兴S^bc

⫹³2关关D₀,X^ab兴S,X^c兴S⫹³2关关D₀,X^a兴S,X^bc兴S⫹3关关D₀^a,X^b兴S,X^c兴S⫹3关关D₀,X^a兴S^b,X^c兴S⫹3关关D₀,X^a兴S,X^b兴S^c

⫹关关关D₀,X^a兴S,X^b兴S,X^c兴S兲, 共35兲

D^abcd⫽P_abcd共D₀^abcd⫹关D₀,X^abcd兴S⫹4关D₀^a,X^bcd兴S⫹4关D₀,X^abc兴S^d⫹6关D₀^ab,X^cd兴S⫹12关D₀^a,X^bc兴S^d⫹6关D₀,X^ab兴S^cd

⫹4关D₀^abc,X^d兴S⫹12关D₀^ab,X^c兴S^d⫹12关D₀^a,X^b兴S^cd⫹4关D₀,X^a兴S^bcd⫹2关关D₀,X^abc兴S,X^d兴S⫹2关关D₀,X^a兴S,X^bcd兴S

⫹3关关D₀,X^ab兴S,X^cd兴S⫹6关关D₀^a,X^bc兴S,X^d兴S⫹6关关D₀,X^ab兴S^c,X^d兴S⫹6关关D₀,X^ab兴S,X^c兴S^d⫹6关关D₀^a,X^b兴S,X^cd兴S

⫹6关关D₀,X^a兴S^b,X^cd兴S⫹6关关D₀,X^a兴S,X^bc兴S^d⫹6关关D₀^ab,X^c兴S,X^d兴S⫹6关关D₀,X^a兴S^bc,X^d兴S⫹6关关D₀,X^a兴S,X^b兴S^cd

⫹12关关D₀^a,X^b兴S^c,X^d兴S⫹12关关D₀^a,X^b兴S,X^c兴S^d⫹12关关D₀,X^a兴S^b,X^c兴S^d⫹2关关关D₀,X^ab兴S,X^c兴S,X^d兴S

⫹2关关关D₀,X^a兴S,X^bc兴S,X^d兴S⫹2关关关D₀,X^a兴S,X^b兴S,X^cd兴S⫹4关关关D₀^a,X^b兴S,X^c兴S,X^d兴S⫹4关关关D₀,X^a兴S^b,X^c兴S,X^d兴S

⫹4关关关D₀,X^a兴S,X^b兴S^c,X^d兴S⫹4关关关D₀,X^a兴S,X^b兴S,X^c兴S^d⫹关关关关D₀,X^a兴S,X^b兴S,X^c兴S,X^d兴S兲. 共36兲

Clearly, a brute-force calculation of these derivatives would be rather expensive. However, an inspection of the nested commutators that contribute to the nth derivative of D re- veals that the inner commutators are also needed for the de- rivatives of order k⬍n. The nth derivative of the AO density matrix may therefore be calculated in a recursive manner, as a sum of single commutators of matrices already available from lower-order derivatives.

To set up such a recursive scheme, we must determine the derivatives of single S commutators. To fourth order, the required derivatives may be written as

关A,X兴S

a⫽2 P_HASX^a, 共37兲

关A,X兴S

ab⫽2 P_HP_ab关ASX^ab⫹2共AS兲^aX^b兴, 共38兲 关A,X兴S

abc⫽2 P_HP_abc关ASX^abc⫹3共AS兲^aX^bc

⫹3共AS兲^abX^c兴, 共39兲 关A,X兴S

abcd⫽2 P_HP_abcd关ASX^abcd⫹4共AS兲^aX^bcd

⫹6共AS兲^abX^cd⫹4共AS兲^abcX^d兴, 共40兲 where the operator P_Hextracts the Hermitian part of a matrix A:

P_HA⫽¹2共A⫹A^†兲. 共41兲 From these expressions, the derivatives of the AO density matrix may be generated recursively, with extensive reuse of intermediates. Expressions for higher derivatives are easily set up in the same manner.

It should be noted that the complexity of the derivatives Eqs.共33兲–共36兲arises because the overlap matrix depends on the perturbation. For perturbations where the overlap matrix is unaffected by the perturbation, the derivatives take the simpler form:

D^a⫽关D₀,X^a兴S, 共42兲 D^ab⫽P_ab共关D₀,X^ab兴S⫹关关D₀,X^a兴S,X^b兴S兲, 共43兲 D^abc⫽P_abc共关D₀,X^abc兴S⫹ ³2关关D₀,X^ab兴S,X^c兴S

⫹ ³2关关D₀,X^a兴S,X^bc兴S

⫹关关关D₀,X^a兴S,X^b兴S,X^c兴S兲, 共44兲 D^abcd⫽P_abcd共关D₀,X^abcd兴S⫹2关关D₀,X^abc兴S,X^d兴S

⫹2关关D₀,X^a兴S,X^bcd兴S⫹3关关D₀,X^ab兴S,X^cd兴S

⫹2关关关D₀,X^ab兴S,X^c兴S,X^d兴S

⫹2关关关D₀,X^a兴S,X^bc兴S,X^d兴S

⫹2关关关D₀,X^a兴S,X^b兴S,X^cd兴S

⫹关关关关D₀,X^a兴S,X^b兴S,X^c兴S,X^d兴S兲. 共45兲 The only contributions are those that involve exclusively the derivatives of X.

In the evaluation of the density matrices discussed here, we have not taken into account the simplifications that fol- low from the 2n⫹1 rule. According to this rule, the deriva-

tives of X that occur only linearly in the derivatives Eq.

共21兲–共24兲 do not contribute to the energy derivatives and may consequently be omitted.³⁰ In practice, this means that the derivatives of X to order n determine the energy to order 2n⫹1. We shall return to this point later.

B. Derivatives of the reference density matrix

As the molecule is perturbed, we must make sure that the reference density D₀remains a valid density matrix. Oth- erwise, it can no longer be used as the generating density matrix of the exponential parametrization Eq. 共28兲. In par- ticular, it must satisfy the idempotency relation

D₀SD₀⫽D₀ 共46兲

as S changes. In the following, we shall use this relation to determine the derivatives of D₀.

Differentiating Eq.共46兲with respect to the perturbation, we obtain

D₀^a⫽D₀^aSD₀⫹D₀S^aD₀⫹D₀SD₀^a. 共47兲 Some simple manipulation involving the projectors Eqs.共13兲 and共14兲shows that Eq.共47兲uniquely determines the projec- tions of D₀^a onto the occupied–occupied and virtual–virtual orbital spaces

PD₀^aP^†⫽⫺D₀S^aD₀, 共48兲

QD₀^aQ^†⫽0. 共49兲

By contrast, PD₀^aQ^† and QD₀^aP^† are undetermined by the idempotency. We may therefore choose these parts of D₀^a freely without affecting the idempotency condition. A natural choice is to set both equal to zero,

PD₀^aQ^†⫽0, 共50兲

QD₀^aP^†⫽0. 共51兲

With these choices of PD₀^aQ^†and QD₀^aP^†, the first derivative of D₀ is uniquely defined as

D₀^a⫽⫺D₀S^aD₀. 共52兲 It is easily verified that Eq.共52兲satisfies not only the differ- entiated idempotency relation Eq.共7兲but also the differenti- ated symmetry and rank conditions in Eqs.共5兲and共6兲. It is therefore a valid derivative Hartree–Fock density matrix.

Once we have determined the first derivative of the ref- erence density matrix, its higher derivatives follow recur- sively by further differentiation of Eq.共52兲. In the following, we list the derivatives of the reference density matrix to fourth order:

D₀^a⫽⫺D₀S^aD₀, 共53兲 D₀^ab⫽⫺D₀S^abD₀⫹2 P_HP_abD₀S^aD₀S^bD₀, 共54兲 D₀^abc⫽⫺D₀S^abcD₀⫹P_HP_abc共6D₀S^abD₀S^cD₀

⫺6D₀S^aD₀S^bD₀S^cD₀兲, 共55兲

D₀^abcd⫽⫺D₀S^abcdD₀⫹P_HP_abcd共8D₀S^abcD₀S^dD₀

⫹6D₀S^abD₀S^cdD₀ 共56兲

⫺24D₀S^abD₀S^cD₀S^dD₀⫺12D₀S^aD₀S^bcD₀S^dD₀

⫹24D₀S^aD₀S^bD₀S^cD₀S^dD₀). 共57兲 The symmetrization operators are given by Eqs. 共26兲 and 共41兲. The corresponding expressions for higher derivatives are easily generated.

It should be understood that the choice of reference den- sity matrix is never unique—any valid density matrix will do. If the geometry is distorted, for example, we must choose one particular reference density matrix共among the infinitely many possible matrices兲 at each geometry. In this way, we establish a connection between density matrices at different geometries, in much the same way as connections are estab- lished between reference MOs when derivatives are taken in MO theories³¹—see also Refs. 32–34.

The choice of connection can sometimes have an impor- tant effect on the calculation of derivatives. Thus, if the con- nection is chosen such that D₀^abecomes very large compared with D^a, then X^a becomes large as well. Eventually, this situation may lead to numerical instabilities in the calcula- tion of derivatives, when the differentiated density is ob- tained from large contributions that nearly cancel. Our con- nection minimizes the norm of the first derivative Eq.共52兲by setting the occupied–virtual and virtual–occupied blocks equal to zero—see Eqs. 共50兲 and 共51兲. Conceptually, our choice of connection for the AO reference density matrix is similar to that of the natural connection in MO theory, in which the orthonormal orbitals of the perturbed system are required to have maximum overlap with the unperturbed set of MOs.³³

C. Response equations

To calculate the derivatives of the density matrix accord- ing to Eqs.共33兲–共36兲, it remains to discuss the evaluation of the derivatives of the rotation matrix X. Because of the 2n

⫹1 rule, we need only determine these derivatives to second order to calculate the energy to fourth order. Differentiating the variational conditions Eq.共19兲twice, we obtain

P_A共F^aDS⫹FD^aS⫹FDS^a兲⫽0, 共58兲 P_AP_ab共F^abDS⫹FD^abS⫹FDS^ab⫹2F^aD^bS⫹2F^aDS^b

⫹2FD^aS^b兲⫽0. 共59兲 By analogy with Eq. 共41兲, we have here introduced the op- erator P_A, which extracts the anti-Hermitian part of a matrix A:

P_AA⫽¹2共A⫺A^†兲. 共60兲 Noting that the derivatives of the Fock matrix are given by

F^a⫽h^a⫹G^a共D兲⫹G共D^a兲, 共61兲 F^ab⫽h^ab⫹G^ab共D兲⫹G共D^ab兲⫹2 P_abG^a共D^b兲, 共62兲 and introducing the notation

K共A兲⫽P_A关G共A兲DS⫹FAS兴, 共63兲

we find that Eqs. 共58兲and共59兲may be written as

K共D^a兲⫽⫺P_A关h^aDS⫹G^a共D兲DS⫹FDS^a兴, 共64兲 K共D^ab兲⫽⫺P_A关h^abDS⫹G^ab共D兲DS⫹FDS^ab兴

⫺2 P_AP_ab关G^a共D^b兲DS⫹F^aD^bS

⫹F^aDS^b⫹FD^aS^b兴, 共65兲 where the unknown, perturbed densities D^aand D^aboccur on the left-hand side only. To comply with the symmetry, rank, and idempotency conditions, the perturbed densities must be written in the form of Eqs. 共33兲 and共34兲. We therefore re- write the derivatives of the density matrices as

D^a⫽关D₀,X^a兴S⫹D_2n^a _⫹1, 共66兲 D^ab⫽关D₀,X^ab兴S⫹D_2n^ab_⫹1 共67兲 with

D_2n^a _⫹1⫽D₀^a, 共68兲 D_2n^ab_⫹1⫽P_ab共D₀^ab⫹2关D₀^a,X^b兴S⫹2关D₀,X^a兴S^b

⫹关关D₀,X^a兴S,X^b兴S兲, 共69兲 and arrive at the following linear sets of equations:

K共关D₀,X^a兴S兲⫽⫺P_A关h^aDS⫹G^a共D兲DS⫹FDS^a

⫹K共D_2n^a _⫹1兲兴, 共70兲 K共关D₀,X^ab兴S兲⫽⫺P_A关h^abDS⫹G^ab共D兲DS⫹FDS^ab

⫹K共D_2n^ab_⫹1兲兴⫺2 P_AP_ab关G^a共D^b兲DS

⫹F^aD^bS⫹F^aDS^b⫹FD^aS^b兴. 共71兲 These are the Hartree–Fock density-matrix response equa- tions, equivalent to the first- and second-order coupled- perturbed Hartree–Fock equations of MO theory.

The left-hand sides of the first- and second-order re- sponse equations共70兲and共71兲are identical and the same for all perturbations; only the right-hand sides differ depending on the order of the differentiation and the nature of the per- turbation. Consequently, the same linear solver may be used in all cases. However, we note that, since the second-order equation 共70兲contains the first-order perturbed density ma- trix D^a, the first-order equations must be solved before the solution of the second-order equations is attempted.

The Hartree–Fock response equations may be solved by, for example, the conjugate-gradient method, which is well suited to the study of large systems. In solving the equations, redundancies must be properly treated 共by applying projec- tors兲and a suitable preconditioner must be found to speed up convergence.³⁵ The techniques are essentially the same as those needed for the optimization of the Hartree–Fock den- sity matrix. By exploiting the sparsity of the AO matrices, linear scaling is achieved for sufficiently large systems, as discussed in the next section.

V. SELECTED PROPERTIES

In the present section, we consider in more detail the evaluation of a few important molecular properties: the mo-

lecular gradient, the molecular Hessian, the magnetizability tensor, and the nuclear magnetic shielding tensors.

A. The molecular gradient

From Sec. III, we recall that the molecular gradient may be written in the general form

E^a⫽Tr Dh^a⫹ ¹2Tr DG^a共D兲⫹Tr D^aF⫹h_nuc^a 共72兲

⫽Tr Dh^a⫹ ¹2Tr DG^a共D兲⫹Tr共D₀^a⫹关D,X^a兴S兲F⫹h_nuc^a , 共73兲 where we have used Eq.共33兲.共Since, for the undifferentiated density matrix, there is no need to distinguish between D and D₀ when X⫽0, we here use D in place of D₀ except in derivatives.兲According to the 2n⫹1 rule, X^amakes no con- tribution to the gradient and may be omitted in the calcula- tion of the gradient, which therefore takes the form

E^a⫽Tr Dh^a⫹ ¹2Tr DG^a共D兲⫹Tr D₀^aF⫹h_nuc^a 共74兲

⫽Tr Dh^a⫹ ¹2Tr DG^a共D兲⫺Tr DS^aDF⫹h_nuc^a . 共75兲 To obtain the final, more explicit expression, we have used Eq. 共53兲.

Although the vanishing of X^a from the molecular gradi- ent follows from quite general considerations 共the 2n⫹1 rule兲, it is nevertheless instructive to demonstrate explicitly how this occurs for the gradient. Expanding the S commuta- tor, we obtain for the last contribution to the gradient Eq.

共73兲

Tr关D,X^a兴SF⫽Tr共DSX^aF⫺X^aSDF兲

⫽Tr X^a共FDS⫺SDF兲

⫽0, 共76兲

where we have invoked Eq.共19兲. Thus, X^a does not contrib- ute since it is multiplied by the Hartree–Fock variational conditions. Similar simplifications occur to higher orders, where the derivatives of X are multiplied either by the Hartree–Fock conditions or by their derivatives.

For sufficiently large molecules, the cost of evaluating the gradient Eq.共75兲scales asO( M ) where M is the size of the molecule—that is, linearly with the size of the system.

For the first two terms and for the last term in Eq.共75兲, linear scaling follows by noting that these terms may be calculated in the same manner as the energy in Eq. 共1兲, replacing un- differentiated overlap distributions

⍀_␮␯⫽␹_␮*␹_␯ 共77兲

by differentiated ones as appropriate. Since the number of derivatives arising from each overlap distribution is indepen- dent of the size of the system, the cost of evaluating the first, second, and fourth contributions to Eq. 共75兲 scales in the same manner as the cost of evaluating the total energy—that is, linearly共provided the Coulomb interactions are calculated by the FMM or an equivalent method兲. Linear scaling of the term containing S^ain Eq.共75兲follows by similar arguments.

Alternatively, we note that, for a given coordinate x_a, the number of nonzero elements in S^ascales asO(1); similarly,

the number of nonzero elements in DS^a scales asO(1), and so on. Since the number of coordinates x_a scales as O( M ), the cost of evaluating ⫺TrDS^aDF also scales as O( M ). In conclusion, the total molecular gradient may be calculated in a manner that scales linearly with the size of the system.

B. The molecular Hessian

Combining the general expression for the molecular Hessian Eq.共22兲with the 2n⫹1 rule, we obtain

E^ab⫽Tr Dh^ab⫹¹2Tr DG^ab共D兲⫹Tr D_2n^ab_⫹1F

⫹2 P_abTr D^a关h^b⫹G^b共D兲⫹ ¹2G共D^b兲兴⫹h_nuc^ab , 共78兲 where D_2n^ab_⫹1 is given by Eq. 共69兲. Thus, to calculate the molecular Hessian, we need only solve the first-order re- sponse equations Eq. 共70兲. Since the geometrical perturba- tions are real, the perturbed rotation matrix X^a is a real, antisymmetric matrix, containing only contributions from X^R in Eq.共10兲.

We also note that, since the expression for the molecular Hessian has been written in a form that is manifestly sym- metric in the indices a and b, it may be shown that the error in the Hessian is quadratic in the error in the perturbed den- sity matrix D^a—see, for example, Ref. 30. Thus, in order to calculate the Hessian to a numerical precision of 10^⫺6, we need only solve the linear equations to a precision of 10^⫺3. By invoking Eq.共70兲, we may write the molecular Hes- sian Eq. 共78兲in the alternative form

E^ab⫽Tr Dh^ab⫹¹2Tr DG^ab共D兲⫹Tr D₀^abF⫹Tr关D₀^b,X^a兴SF

⫹Tr关D,X^a兴S^bF⫹Tr D₀^b关h^a⫹G^a共D兲兴

⫹Tr D^a关h^b⫹G^b共D兲⫹G共D₀^b兲兴⫹h_nuc^ab , 共79兲 where we distinguish between D₀^aand D^a. In this manner, we may compute E^ab without having to calculate D^b. However, if this expression is used to compute the Hessian, the error in E^ab is no longer quadratic in the error in D^a. This is an important point in favor of Eq. 共78兲since, in the AO basis, the solution of the response equations is often difficult be- cause the electronic Hessian is not diagonally dominant.

Let us now consider the cost of evaluating the molecular Hessian for a large system of size M. As for the molecular gradient in Sec. V A, our discussion of the cost of evaluating the Hessian is qualitative rather than quantitative, its purpose being to indicate in general terms how linear scaling arises for large systems.

Let us consider, for a given coordinate x_aof atom A, the number of significant elements on the right-hand side of the first-order response equations in the form Eq. 共64兲. We first note that, for the derivative overlap matrix S^a, the number of nonzero elements S_␮␯^a scales asO(1) since only overlap dis- tributions⍀_␮␯ with AOs on A contribute. Next, we consider the scaling of the matrix

f^a⫽h^a⫹G^a共D兲. 共80兲 From Eqs.共2兲–共4兲, we see that f^ahas three parts: the kinetic, Coulomb, and exchange parts. For the kinetic and exchange parts, O(1) scaling of the number of significant elements