Analytical Hartree-Fock gradients with respect to the cell parameter systems periodic in on

a r X i v :p h y s i c s /0512172v 1 [p h y s i c s .c o m p -p h ] 19 D e c 2005

Analytical Hartree-Fock gradients with respect to the cell parameter:systems

periodic in one and two dimensions

K.Doll

Institut f¨u r Mathematische Physik,TU Braunschweig,Mendelssohnstra?e 3,D-38106Braunschweig,Germany

R.Dovesi,R.Orlando

Dipartimento di Chimica IFM,Universit`a di Torino,Via Giuria 5,I-10125Torino,Italy Analytical Hartree-Fock gradients with respect to the cell parameter have been implemented in

the electronic structure code CRYSTAL,for the case of one and two-dimensional periodicity.As in most molecular codes,Gaussian type orbitals are used to express the wavefunction.Examples demonstrate that the gradients have a good accuracy.

Keywords:Hartree-Fock,cell gradient,periodic systems,CRYSTAL

I.INTRODUCTION

Analytical gradients 1–7have become a standard tool in molecular quantum chemistry.They are indispensable for the optimization of structures,and many properties can be e?ciently computed with the help of analytical derivatives.The ?eld was pioneered by Pulay 8;the the-ory had already been derived earlier independently 9.The traditional quantum chemical methods are dif-?cult to apply to solids because of the large increase of the computational e?ort with the system size.Af-ter several decades of development,Hartree-Fock calcu-lations for solids can nowadays be routinely performed with the CRYSTAL code 10,11.Although Hartree-Fock calculations often have large errors due to the neglect of electronic correlation,a large interest has grown in the past few years due to the success of hybrid functionals which include an admixture of exact (Fock)exchange.Analytical gradients in the CRYSTAL code were ?rst implemented with respect to nuclear positions 12,13,and after the implementation of a scheme for geometry opti-mization,an e?cient structural optimization could be performed 14.In periodic systems,the cell parameter is another variable to be optimized.The ?rst gradi-ents with respect to the cell parameter,at the Hartree-Fock level,were for systems periodic in one dimension 15.Various groups have implemented these gradients in one dimension 16,17(see also the recent review article 18)or in two dimensions 19.For the general case,a strategy to compute cell parameter derivatives (and thus the stress tensor)was suggested with point charges 20,and an al-gorithm for structural optimization,based on redundant internal coordinates was proposed 21.Second analytical derivatives with respect to the cell parameter have also been implemented recently 22.

A ?rst big step of the corresponding implementa-tion in the CRYSTAL code were analytical Hartree-

Fock gradients with respect to the cell parameter in three dimensions23.It is important to note that the CRYSTAL code is based on the Ewald24,25method in three dimensions,so that computing analytical gradi-ents with respect to the cell parameter requires vari-ous additional derivatives:for example the reciprocal lattice vectors depend on the cell parameter,and var-ious others.This requires additional derivatives which were not yet available with the implementation of nu-clear gradients,and this has been documented in great detail23.The one and two-dimensional case are again dif-ferent because di?erent potentials are used:Parry’s po-tential in two dimensions26,27,and Saunders’potential in one dimension28.Parry’s potential is similar to Ewald’s potential,but modi?ed for the case of two dimensions. Saunders’potential relies on a real space approach. This article is intended to complement the?rst article on cell gradients23.Many parts have already been de-scribed in the?rst article,and therefore the main empha-sis is to delineate the di?erences due to the dimension-ality.The article consists thus of one section about the general di?erences to the three-dimensional case,one sec-tion about the two-dimensional case,one section about the one-dimensional case,and one section with examples. II.GENERAL DIFFERENCES WITH RESPECT TO THE THREE-DIMENSIONAL CASE

The main di?erence to the three-dimensional case is the way how the Coulomb energy is computed.The ex-pression to be evaluated is the Coulomb energy per cell:

1

E coul=

spin is ignored in the equations in this article.The ex-tension is straightforward,as was shown for the three-dimensional case23.Examples for spin-polarized calcula-tions are given in section V.

The potential function enters via the nuclear-nuclear repulsion(equation10in reference23),the nuclear at-traction integrals(equation34in reference23),and the ?eld integrals(equation43in reference23).Essentially, the derivatives are computed as described in the previous article23,there are only minor di?erences as described in section III and IV.

The derivatives of the other integrals(overlap,kinetic energy,multipoles,bielectronics)and the calculation of the energy-weighted density matrix is practically identi-cal to the three-dimensional case23.

Finally,the correction due to the spheropole(equation 47in reference23)is zero in one and two dimensions and thus does not have to be discussed.The spheropole is a correction which arises due to the Ewald method,when applied to the electronic charge distribution:the charge distribution is approximated by multipoles in the long range,and not approximated in the short range.The electrostatic potential is then computed as the sum of the Ewald potential of the multipoles and of the Coulomb potential of the charge distribution in the short range. Replacing the Ewald potential with the Coulomb poten-tial is correct,if the di?erence of multipolar charge dis-tribution and the exact charge distribution in the short range,has zero charge,dipole,quadrupole,and second spherical moment25.The second spherical moment can also be seen as the average electrostatic potential of a charge distribution(see the discussion in section3.2of reference25).Here,it corresponds to the average elec-trostatic potential of the di?erence of the exact and the approximated charge distribution.This term is?nite and in general non-zero,in the case of periodicity in three di-mensions.However,when the system has periodicity in less than three dimensions,the average electrostatic po-tential of a charge distribution with zero total charge, dipole and quadrupol,is zero.Therefore,there is no spheropole in less than three dimensions.

This can also seen from equation31in reference25. The average Coulomb potential is obtained as follows:

2π

Φ=?

the average Coulomb potential becomes zero,and the spheropole correction becomes zero.

It should be mentioned,that two-dimensional period-icity is implemented in the CRYSTAL code in such a way that there is only one slab which is not repeated in the third dimension.Still,the argument presented above holds in a similar way,and there is thus no spheropole correction in systems with less than three-dimensional periodicity.

The total energy is thus similar to the three-dimensional case 23,apart from the spheropole term which is zero:

E total =E kinetic +E NN +E coul ?nuc +E coul ?el +E exch ?el == g ,μ,ν

P ν g μ 0T μ 0ν g

+12

g ,μ,ν

P ν g μ 0

h,τ,σ

P σ hτ 0C μ 0ν g τ 0σ h ?

c

L l =0l m =?l

ηm l (ρc ; A c )M m lμ 0ν g c

?1

?a ij =

?

g ,μ,ν

P ν g μ 0

?T μ 0ν g

?a ij

?

g ,μ,ν

P ν g μ 0

?N μ 0ν g

2

g ,μ,ν

P ν g μ 0

τ,σ

P σ hτ 0

?C μ 0ν g τ 0σ h

?a ij

φτ( r ? A τ)φσ( r ? A σ? h )X m l ( r ? A c )d 3r M m lμ 0ν g c

+1

?a ij

? g ,μ,ν

?S μ 0ν g

III.THE TWO-DIMENSIONAL CASE

In the two-dimensional case,the primitive cell is given by two vectors,with two components: a1, a2.a ij are de?ned in such a way that a11=a1x is the x-component of a1,a12=a1y the y-component of a1,a21is the x-component of a2,and a22is the y-component of a2. a1 a2 = a1x a1y a2x a2y = a11a12a21a22 (7) A point g of the direct lattice is de?ned as g=n1 a1+ n2 a2,with n1,n2being integer numbers.The position of an atom c in a cell at the origin(i.e. g= 0)is given as

A c =f c,1 a1+f c,2 a2,and then in cell g the position will

be:

A c + g=(f c,1+n g,1) a1+(f c,2+n g,2) a2

We have used an additional index,i.e.n g,1means fac-tor n1of the lattice vector g.The cartesian t component (with t being x or y)of the vector A c+ g,indicated as A c,t+g t,is thus

A c,t+g t= 2m=1(f c,m+n g,m)a mt

As all the integrals depend on the position of the nuclei, the derivatives of the nuclear coordinates with respect to the cell parameters are required:

?A c,t+g t

γ| r? A a? h|)

2V| K| exp(2π| K|(z?A a,z)) 1?erf(√√

γ(z?A a,z)?π| K|γ

?2πγ(z?A a,z))?2√

V√

| r? A a? h|

is omitted from the sum.In the

reciprocal lattice series,the prime indicates that all terms with K= 0are included.

The error function erf is de?ned as in reference23, equation12.

Like the Ewald potential,Parry’s potential depends on the variables A c, h,V,γand K.The derivative with respect to the cell parameters thus requires deriva-tives with respect to these variables.For the deriva-tives with respect to A c and h this is like in the three-dimensional case.There are minor changes due to the two-dimensionality for the derivatives of the area V,of the K-vectors with respect to a ij and of the screening parameterγ.

1.Derivative of the area

The area is obtained as the magnitude of the cross product of the cell parameters:

V=| a1× a2|=|a1x a2y?a1y a2x|(10) If we assume that a1x a2y?a1y a2x is positive,then the derivatives?V

?a1x

=a2y(11)

?V

?a2x

=?a1y(13)

?V

V (a2y,?a2x); b2=

2π

?a1x

=?

?V

V

(17)

? b1

?a1y b 1

?a2x =?

?V

V

+

2π

?a2y =?

?V

V

+

2π

?a1x =?

?V

V

+

2π

?a1y =?

?V

V

+

2π

?a2x

=?

?V

V

(23)

? b2

?a2y

b

2

?a ij =

?γ

?a ij

=?

γ

?a ij

(25)

As a whole,Parry’s potential leeds to similar terms appearing in the derivatives as in the case of the Ewald potential.This is what was to be expected,as Parry’s potential is essentially obtained when Ewald’s approach to treat the Coulomb interaction is applied to a system with two-dimensional periodicity.

IV.THE ONE-DIMENSIONAL CASE

In the one-dimensional case,there is only one cell pa-rameter:a xx=a11=a=| a|.This case is somewhat dif-ferent from the two-and three-dimensional case because a pure real space approach is used in the CRYSTAL code for the potential to describe the Coulomb interaction28. The potential consists of a point charge+1,neutralized by a uniform charge distribution of length a,with charge density?1

Euler-MacLaurin summation rule.As a whole,the fol-lowing expression was obtained28:

Φ( r)=

M

n=?M′1a+ξ(M, r)+ξ(M,? r)(26)

The?rst term comprises the exact part,the next two (with the H function)the region due to the uniform charge density in the range of the exact sum(from?M a to M a),the remaining two terms(theξ-function)are the approximated part.The prime indicates that terms with | r?n a|=0are omitted.M is thus the number of cells, where the sum is performed exactly,and U=a(M+1

p2+α+p).ξ(M, r)and ξ(M,? r)are contributions from the long range part, which is approximated by the Euler-MacLaurin rectange rule summation formula.For more details,see reference 28.For the present purpose,it is important to note that the direct lattice vector a appears in the potential,but no screening parameterγand no reciprocal lattice vec-tors K as in the two-and three-dimensional case.This means that derivatives with respect to the nuclear co-ordinates A c and derivatives with respect to the direct lattice vectors n a appear,which are essentially given by the nuclear gradients,multiplied with the fractional co-ordinates.The derivatives with respect to a due to the H andξfunction are very lengthy,but still straightfor-ward.They are thus not discussed here,but formulas can be derived from Saunders’article28.

V.EXAMPLES

In this section,we give some numerical examples of the accuracy of the gradients.The tests considered are es-sentially identical or similar to the test cases distributed with the CRYSTAL code and with the ones from ref-erence14.Note that the fractional coordinates of the atoms were not optimized.

First,two systems with one-dimensional periodicity are considered.In table I,SN is periodically arranged. The analytical and numerical derivative agree well up to4 digits,and the minimum of the energy at a=4.42?A agrees with the place where the gradient changes its sign.In ta-ble II,such a comparison is done for polyglycine.The agreement of numerical and analytical gradients is simi-lar to SN,and again the vanishing of the gradient agrees with the minimum of the energy,to at least0.01?A.In table III,ferromagnetic NiO is studied at the level of unrestricted Hartree-Fock.The agreement of numerical and analytical gradient can be improved by increasing the ”ITOL”-parameters11,as described earlier12,23.Indeed, when increasing them from default values to higher ones, symmetric in ITOL4and ITOL5,then analytical and nu-merical gradient match better.Note that,when running at lower ITOL parameters,an inaccuracy is introduced

in the total energy expression and thus in the numerical gradients as well.The fact that numerical and analyti-cal gradients match less well at low ITOL values is thus a combination of an inaccuracy in the energy expression (which a?ects the numerical gradient)and an inaccuracy in the analytical gradient.Still,in all the tests performed so far,no severe error was found when using default val-ues for the ITOL https://www.360docs.net/doc/2e4035982.html,ing higher ITOL pa-rameters is mainly useful for tests of the correctness of the code.

Then,various systems with two-dimensional periodic-ity are considered.In table IV,3MgO layers are con-sidered.Numerical and analytical derivative agree to3 digits,and the minimum of the energy and the vanish-ing of the gradient agree also well.The same accuracy is found for Al2O3in table V,where a slab with6atomic layers is considered.In table VI,a Cr2O3slab was chosen as an example for unrestricted Hartree-Fock.The accu-racy is slightly worse when comparing the numerical and the analytical gradient.This can again be improved by increasing the”ITOL”-parameters.The minimum in the energy agrees already with default”ITOL”values to at least0.01?A.Finally,in table VII,LiF was arranged with two dimensional periodicity,without symmetry,in such a way that three components of the cell gradient(a1x, a1y,a2y)can be computed independently.This test thus demonstrates that these components are correctly com-puted.

In table VIII,the CPU times are displayed.The cal-culations were performed on a single CPU of a Compaq ES45,with a clock rate of1GHz.As in the three-dimensional case,we compare again the CPU time for the integrals with the time for the gradients.The CPU time for all the gradients(nuclear and cell gradients)is roughly?ve to ten times the CPU time for the integrals. This may become smaller in the future with further opti-mizations in the gradient code.Note that the CPU time for the self consistent?eld calculations is relatively high because a very low convergence threshold was chosen in order to ensure the accuracy of the succeeding gradient calculation(the gradient calculation is the more accurate, the more accurately the self consistent?eld equations are solved).

The CPU times thus indicate that analytical gradients can be computed at a relatively low https://www.360docs.net/doc/2e4035982.html,pared with numerical gradients,it appears that analytical gra-dients should usually be favorable,especially because nu-merical gradients will depend on the step size,and often it will be necessary to break a symmetry for a?nite dis-placement,to compute the numerical gradient.Numeri-cal gradients require at least one additional energy eval-uation for each coordinate to be optimized,which makes analytical gradients clearly favorable,if there is a large number of geometrical parameters.

VI.CONCLUSION

A formalism for the calculation of the analytical gra-dient of the Hartree-Fock energy,with respect to the cell parameter,has been presented and implemented in the code CRYSTAL,for the case of systems periodic in one and two dimensions.The implementation includes the cases of spin-restricted and unrestricted polarization.It was shown that a high accuracy can be achieved.

VII.ACKNOWLEDGMENT

The calculations were performed on a Compaq ES45 (computer center of the TU Braunschweig).

19M.Tobita,S.Hirata,and R.J.Bartlett,J.Chem.Phys.

118,5776(2003).

20K.N.Kudin and G.E.Scuseria,Phys.Rev.B61,5141 (2000).

21K.N.Kudin,G.E.Scuseria,and H.B.Schlegel,J.Phys. Chem.114,2919(2001).

22D.Jacquemin,J.-M.Andr′e,and B.Champagne,J.Chem. Phys.118,373(2003).

23K.Doll,R.Dovesi and R.Orlando,Theor.Chem.Acc.

112,394(2004).

24P.P.Ewald,Ann.Phys.(Leipzig)64,253(1921).

25V.R.Saunders,C.Freyria-Fava,R.Dovesi,L.Salasco,and C.Roetti,Mol.Phys.77,629(1992).

26D.E.Parry,Surf.Science49,433(1975);54,195(1976) (Erratum).

27D.M.Heyes,M.Barber,and J.H.R.Clarke,J.Chem. Soc.Faraday Trans.II,73,1485(1977).

28V.R.Saunders,C.Freyria-Fava,R.Dovesi,and C.Roetti, https://www.360docs.net/doc/2e4035982.html,m.84,156(1994).

TABLE I.SN,with one-dimensional periodicity.A comparison of analytical and numerical gradient is done for various unit cell lengths.A [3s 2p 1d ]basis set was used for S,and a [2s 1p ]basis set for N.a analytical derivative

numerical derivative

energy [?A ][E h /a 0][E h /a 0][E h ]4.300.041440.0414-893.8700814.410.003720.0037-893.8746394.420.000640.0006-893.8746804.43-0.00238-0.0024-893.8746634.500

-0.02208

-0.0221

-893.873013

TABLE II.Polyglycine.A comparison of analytical and numerical gradient is done for various unit cell lengths.Basis sets of the size [2s 1p ]were used for C,O,N and a [1s ]basis set for H.a analytical derivative

numerical derivative

energy [?A ][E h /a 0][E h /a 0][E h ]7.300.019560.0196-408.2201737.420.001160.0012-408.2224957.43-0.00030-0.0003-408.2225037.44-0.00175-0.0017-408.2224847.50

-0.01018

-0.0102

-408.221807

TABLE III.NiO,ferromagnetic,unrestricted Hartree-Fock.The gradient with respect to the cell parameter is computed for two di?erent values of the ITOL parameters.A [5s 4p 2d ]basis set for Ni was used,and a [4s 3p ]basis set for O.a analytical derivative

numerical derivative

energy [?A ][E h /a 0]

[E h /a 0]

[E h ]ITOL 666612(default)

5.00-0.10864-0.1074

-1581.454974ITOL 6661212

5.00

-0.10782

-0.1078

-1581.456358

TABLE IV.MgO surface,3atomic layers.The unit cell consists of 3Mg and 3O atoms,with a 1x =a 2y =a .Basis sets

of the size [3s 2p ]were used.The derivative with respect to

?

?a 1x ?a 1x ?a 2y ?a 2y ?a 1y

and ?

TABLE V.Al 2O 3,6atomic layers.The unit cell consists of 6Al and 4O atoms,with a 1x =

√

3/2?a .Basis

sets of the size [3s 2p 1d ]for Al and [2s 1p ]for O were chosen.The derivative with respect to

?

?a 1x

?a 1x

?a 2y

?a 2y

?a 1y

and

?

3/2?a 2y =

√

?a =??a

+?

?a

is displayed,

?

?a 2x

do not contribute.

a

analytical derivative

numerical derivative

energy [?A ][E h /a 0]

[E h /a 0]

[E h ]ITOL 666612(default)

4.700.134650.1379-4622.5897854.870.004260.0069-4622.6122784.88-0.002530.0001-4622.6123394.89-0.00921-0.0066-4622.612277

5.00-0.07676-0.0745

-4622.603638ITOL 6661212

4.88-0.00116-0.0011-4622.617935

5.00

-0.07539-0.0754

-4622.609006

TABLE VII.LiF,with a unit cell of a 1x =5?A ,a 2y =4?A ,and an angle of 60?,resulting in a 1y =2.5?A .The F atoms are

at (x=0.1,y=0(x and y in fractional units),z=0.1?A ),(x=0.5,y=0.5(x and y in fractional units),z=0.3?A ),the Li atoms at

(x=0.5,y=0(x and y in fractional units),z=0.2?

A ),and (x=0,y=0.5(x and y in fractional units),z=0.4?A ).A [2s 1p ]basis set was used for Li,a [4s 3p ]basis set for https://www.360docs.net/doc/2e4035982.html,ponent

analytical derivative

numerical derivative

[E h /a 0]

[E h /a 0]

?E

?a 1y -0.04415-0.0441

?E

TABLE VIII.CPU times for one single point calculation of the various systems.The calculations were performed on a Compaq ES45,using a single CPU(1GHz).The CPU times refer to the part for the integrals(all the integrals were written to disk),the self-consistent?eld(SCF)procedure,and to the calculation of all the gradients(i.e.nuclear gradients and cell gradients).

system CPU time,in seconds

integrals SCF gradients SN116

Polyglycine2417 NiO2149

MgO5352

Al2O381278

Cr2O327153176

LiF31820