******************************************************************************** * * * 2DHF version 1-2003 * * Copyright (C) 1996 Jacek Kobus, Leif Laaksonen, Dage Sundholm * * * * This software may be used and distributed according to the terms * * of the GNU General Public License, see README and COPYING. * * * ******************************************************************************** c ### totenp ### c c Calculate the total energy in the case of exchange c potentials being kept on a disk subroutine totenp (iorb,psi,pot,excp,e,f0,wgt1,wgt2, & wk0,wk1,wk2,wk3) implicit integer*4 (i-n) implicit real*8 (a-h,o-z) include 'commons8.inc' dimension psi(*),pot(*),excp(*),e(*),f0(*),wgt1(*),wgt2(*), & wk0(*),wk1(*),wk2(*),wk3(*) data eps /1.d-7/ if (iorb.eq.0) then ictot=0 vkt =0.d0 vnt =0.d0 evt =0.d0 epott=0.d0 return endif if (islat.ne.0) then c the Slater exchange approximation is used call totensl (psi,pot,excp,e,f0,wgt1,wgt2,wk0,wk1,wk2,wk3) return endif engt(1)=eng(1) if (nel.lt.2) return ictot=ictot+1 woneel=0.d0 wdcoul=0.d0 engt(1)=0.d0 c calculate first contributions from one particle operators and c coulomb potential contributions within the same shell if (idbg(70).ne.0) then write(*,*) 'orb t v(n) t+v(n) ', & ' v(c-within shell)' endif norb2=norb*norb iborb=i1b(iorb) isiorb=i1si(iorb) nmut=i1mu(iorb) ibpot=i2b(iorb) isipot=i2si(iorb) oc=occ(iorb) ipe=mgx(6,iorb) ipex=ipe-2*(ipe/2) c calculate derivatives over mu and ni variables by means of matrix c multiplication if (ipex.eq.0) then isym= 1 else isym=-1 endif call fill (nni,nmut,isym,psi(iborb),wk3) call difni (nmut,wk3,wk0,wk1,wk2) call refill(nni,nmut,wk1,wk0) call difmu (nmut,wk3,wk2) call refill(nni,nmut,wk0,wk2) c add derivatives over mu and ni call add (isiorb,wk0,wk1) c add contribution from phi part of laplasian c e enters the expression with minus sign which is already c incorporated in e if (ipe.ne.0) then w=dble(ipe*ipe) call dcopy (isiorb,e,1,wk0,1) if (ipe.ne.1) then call dscal (isiorb,w,wk0,1) endif call prod (isiorb,psi(iborb),wk0) call add (isiorb,wk0,wk1) endif call dcopy (isiorb,wk1,1,wk2,1) call prod (isiorb,psi(iborb),wk2) w=ddot(isiorb,wgt1,1,wk2,1) vk(iorb)=w call dcopy (isiorb,f0,1,wk0,1) call prod (isiorb,psi(iborb),wk0) call prod2 (isiorb,psi(iborb),wk0,wk2) w=ddot(isiorb,wgt1,1,wk2,1) vn(iorb)=w call add (isiorb,wk0,wk1) call prod (isiorb,psi(iborb),wk1) w = ddot(isiorb,wgt1,1,wk1,1) woneel=oc*w c contribution from coulomb interaction within the same shell if (oc.gt.1.d0) then if (isipot.ne.isiorb) then write(iout6,*) 'coulomb potentials and orbitals have to ', & 'defined on the same number of subgrids' stop 'toten' endif call dcopy (isipot,pot(ibpot),1,wk2,1) call prod (isiorb,psi(iborb),wk2) call prod (isiorb,psi(iborb),wk2) w=ddot(isiorb,wgt2,1,wk2,1) c wdcoul=wdcoul+oc*(oc-1.d0)/2.d0*w wdcoul= oc*(oc-1.d0)/2.d0*w endif vkt=vkt+oc*vk(iorb) vnt=vnt+oc*vn(iorb) if (idbg(70).ne.0) then write(*,6000) iorb, vk(iorb),vn(iorb),vk(iorb)+vn(iorb),w endif 10 continue 6000 format(' ',i4,2x,4d20.9) 6010 format(' ',i4,2x,4d20.9) 6012 format(' ',2i4,2x,4d20.9) c contribution from coulomb and exchange interaction between shells wndc =0.d0 wex1 =0.d0 wex2 =0.d0 iorb1=iorb iborb1=i1b(iorb1) isiorb1=i1si(iorb1) nmut1=i1mu(iorb1) ibpot1=i2b(iorb1) isipot1=i2si(iorb1) oc1=occ(iorb1) ipe1=mgx(6,iorb1) iex1 =iorb1+iorb1*(iorb1-1)/2 call dcopy (isiorb1,psi(iborb1),1,wk0,1) call prod (isiorb1,psi(iborb1),wk0) c calculate exchange interaction within pi, delta, etc. open shell if (ipe1.gt.0.and.abs(oc1-1.d0).gt.eps) then call exint (iorb1,ocx1) ibex=i3b(iex1) isiex1=i3si(iex1) ngrid=min(isiorb1,isiex1) call prod2 (ngrid,wk0,excp(ibex),wk1) w=ddot (ngrid,wgt2,1,wk1,1) wex1=wex1+ocx1*w if (idbg(70).ne.0) then write(*,6020) iorb1,w 6020 format(' ','orb v(x-within shell) ',i4,2x,4d15.5) endif endif do iorb2=iorb1+1,norb iborb2=i1b(iorb2) isiorb2=i1si(iorb2) nmut2=i1mu(iorb2) ibpot2=i2b(iorb2) isipot2=i2si(iorb2) oc2=occ(iorb2) c coulomb interaction between shells ngrid=min(isipot2,isiorb1) call prod2 (ngrid,pot(ibpot2),wk0,wk1) w=ddot(ngrid,wgt2,1,wk1,1) wndc=wndc+oc1*oc2*w if (idbg(70).ne.0) then write(*,6022) iorb1,iorb2,w 6022 format(' ','orb1 orb2 v(c-inter shell) ',2i4,2x,4d15.5) endif ipe2=mgx(6,iorb2) iex=iorb1+iorb2*(iorb2-1)/2 ibex=i3b(iex) isiex=i3si(iex) c exchange interaction between shells (same lambda) call excont (iorb1,iorb2,ocx1,ocx2) ngrid=min(isiorb1,isiorb2,isiex) call prod2 (ngrid,excp(ibex),psi(iborb1),wk1) call prod (ngrid,psi(iborb2),wk1) w=ddot(ngrid,wgt2,1,wk1,1) wex1=wex1+ocx1*w if (idbg(70).ne.0) then write(*,6024) iorb1,iorb2,w 6024 format(' ','orb1 orb2 v(x-inter shell)a',2i4,2x,4d15.5) endif if (ilc(iex).gt.1) then call prod2 (ngrid,excp(ibex+ngrid),psi(iborb1),wk1) call prod (ngrid,psi(iborb2),wk1) w=ddot (ngrid,wgt2,1,wk1,1) wex2=wex2+ocx2*w if (idbg(70).ne.0) then write(*,6026) iorb1,iorb2,w 6026 format(' ','orb1 orb2 v(x-inter shell)b',2i4,2x,4d15.5) endif endif enddo evt=evt+woneel+wdcoul+wndc-wex1-wex2 epott=epott+oc*vn(iorb)+wdcoul+wndc-wex1-wex2 virrat=(epott+z1*z2/r)/vkt if (ictot.eq.norb) then iready=ictot etot=evt+z1*z2/r engt(1)=evt ictot=0 vkt =0.d0 vnt =0.d0 evt =0.d0 epott=0.d0 endif if (idbg(70).ne.0) then write(*,*) 'total energy: ', evt write(*,*) 'contributions: one-electron coulomb', & ' exchange ' write(*,*) woneel,wdcoul+wndc,wex1+wex2 endif return end