[CPMD-list] LSD - KS diag.
Juerg Hutter
hutter at pci.unizh.ch
Fri May 23 10:53:16 CEST 2003
Hi Jochen
I had some troubles finding this bug but finaly I was
able to track it down. The fixed file vofrhob.F
is attached.
best wishes
Juerg
----------------------------------------------------------
Juerg Hutter Phone : ++41 1 635 4491
Physical Chemistry Institute FAX : ++41 1 635 6838
University of Zurich E-mail: hutter at pci.unizh.ch
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland
----------------------------------------------------------
On Mon, 19 May 2003, Jochen Blumberger wrote:
> Dear CPMD developers,
>
>
> I encountered a problem for the Kohn-Sham diagonlization with LSD in V3.7.1:
>
> The wavefucntion optimization for a sample of 32 waters gives
>
> ELECTRONIC GRADIENT:
> MAX. COMPONENT = 5.10300E-06 NORM = 7.29372E-08
> NUCLEAR GRADIENT:
> MAX. COMPONENT = 4.18221E-02 NORM = 1.19544E-02
>
>
> TOTAL INTEGRATED ELECTRONIC DENSITY
> IN G-SPACE = 256.000000
> IN R-SPACE = 256.000000
> TOTAL INTEGRATED SPIN DENSITY .000026
>
> (K+E1+L+N+X) TOTAL ENERGY = -548.80201537 A.U.
> (K) KINETIC ENERGY = 412.08735410 A.U.
> (E1=A-S+R) ELECTROSTATIC ENERGY = -367.57551500 A.U.
> (S) ESELF = 404.26149474 A.U.
> (R) ESR = 24.02130066 A.U.
> (L) LOCAL PSEUDOPOTENTIAL ENERGY = -513.72725987 A.U.
> (N) N-L PSEUDOPOTENTIAL ENERGY = 54.97660217 A.U.
> (X) EXCHANGE-CORRELATION ENERGY = -134.56319678 A.U.
> GRADIENT CORRECTION ENERGY = -7.53087044 A.U.
>
> ****************************************************************
>
>
> but the diagonalization using the Restart file of the optimization run
> gives:
>
> TOTAL INTEGRATED ELECTRONIC DENSITY
> IN G-SPACE = 256.000000
> IN R-SPACE = 256.000000
> TOTAL INTEGRATED SPIN DENSITY .000026
>
> (B+E2+X-V) TOTAL ENERGY = -549.88465526 A.U.
> (B) BAND ENERGY = -88.88410491 A.U.
> (E2=I-H-S+R) ELECTROSTATIC ENERGY = -499.36486171 A.U.
> (S) ESELF = 404.26149474 A.U.
> (R) ESR = 24.02130066 A.U.
> (X) EXCHANGE-CORRELATION ENERGY = -134.56319678 A.U.
> (V) EXCHANGE-CORRELATION POTEN. = -172.92750821 A.U.
> GRADIENT CORRECTION ENERGY = -7.53087044 A.U.
>
> ==------------------------------------------------------------==
> == NFI= 1 ETOT= -549.884655 TCPU= 14.14 ==
> == DRHOMAX= 3.550E-09 DETOT= 0.000E+00 THL= 0.000E+00 ==
> ==------------------------------------------------------------==
>
>
> In the LDA case wavefunction optimization and diagonlization gives the
> same total energy which also coincides with the total energy of the LSD
> wavefunction optimization.
> Here's the output for LDA diagonalization:
>
> TOTAL INTEGRATED ELECTRONIC DENSITY
> IN G-SPACE = 256.000000
> IN R-SPACE = 256.000000
>
> (B+E2+X-V) TOTAL ENERGY = -548.80201716 A.U.
> (B) BAND ENERGY = -88.88413317 A.U.
> (E2=I-H-S+R) ELECTROSTATIC ENERGY = -499.36483919 A.U.
> (S) ESELF = 404.26149474 A.U.
> (R) ESR = 24.02130066 A.U.
> (X) EXCHANGE-CORRELATION ENERGY = -134.56321558 A.U.
> (V) EXCHANGE-CORRELATION POTEN. = -174.01017083 A.U.
> GRADIENT CORRECTION ENERGY = -7.53087727 A.U.
>
> ==------------------------------------------------------------==
> == NFI= 1 ETOT= -548.802017 TCPU= 7.11 ==
> == DRHOMAX= 4.225E-09 DETOT= 0.000E+00 THL= 0.000E+00 ==
> ==------------------------------------------------------------==
>
>
> The KS energy levels of the LSD and LDA diagonalization are ~identical
> (as they should be) but the total energies are 549.88465526 for LSD and
> 548.80201716 for LDA. The difference comes from the
> EXCHANGE-CORRELATION POTENTIAL: -172.92750821 for LSD, -174.01017083 for
> LDA. Could there be a bug in the calculation of the total
> exch.corr.Potential in the LSD case ? (In previous versions LSD diag.
> worked well)
>
> many thanks
>
> best regards
>
> Jochen Blumberger
>
> _______________________________________________
> CPMD-list mailing list
> CPMD-list at cpmd.org
> http://www.cpmd.org/mailman/listinfo/cpmd-list
>
-------------- next part --------------
C ==================================================================
SUBROUTINE VOFRHOB(FION,RHOE,V,VTEMP,SCR,LSCR,TFOR,TSTRESS)
C ==--------------------------------------------------------------==
C == COMPUTES ==
C == THE ONE-PARTICLE POTENTIAL V IN REAL SPACE ==
C == THE TOTAL ENERGY ETOT ==
C == THE FORCES FION ACTING ON THE IONS ==
C ==--------------------------------------------------------------==
C == RHOE: in electronic density in real space ==
C == out total electronic potential in real space ==
C ==--------------------------------------------------------------==
IMPLICIT NONE
INCLUDE 'system.h'
INCLUDE 'spin.inc'
INCLUDE 'ions.inc'
INCLUDE 'ener.inc'
INCLUDE 'elct.inc'
INCLUDE 'cppt.inc'
INCLUDE 'pslo.inc'
INCLUDE 'simul.inc'
INCLUDE 'fft.inc'
INCLUDE 'geq0.inc'
INCLUDE 'cvan.inc'
INCLUDE 'nlps.inc'
INCLUDE 'sfac.inc'
INCLUDE 'nlcc.inc'
INCLUDE 'strs.inc'
INCLUDE 'str2.inc'
INCLUDE 'efld.inc'
INCLUDE 'extpot.inc'
C Arguments
INTEGER LSCR
REAL*8 FION(3,NAX,NSX),
& RHOE(NNR1,NLSD),SCR(LSCR)
COMPLEX*16 V(MAXFFT,NLSD),VTEMP(NHG,NLSD)
LOGICAL TFOR,TSTRESS
C Variables
DIMENSION F(N),DQG(*)
COMPLEX*16 VTMP(NHG,NLSX)
REAL*8 RHOVAL(NNR1,NLSD),EXTPOT(*),EXTF(*),GRAD(NNR1,NLSD*4)
POINTER (IP_VTMP,VTMP),(IP_GRAD,GRAD),(IP_RHOVAL,RHOVAL)
REAL*8 FNL(IMAGP,LDF1,NHXS,*),DFNL(*),SXC,SGCX,SGCC,VGC
INTEGER ISUB,IG,IR,NNRS,KK,
& IL_VTMP,IL_GRAD,IL_RHOVAL,LENGTH
LOGICAL DEBUG
C ==--------------------------------------------------------------==
CALL TISET(' VOFRHOB',ISUB)
C SCR partition.
DEBUG=.FALSE.
CALL SCRSET(SCR,LSCR,DEBUG)
IL_VTMP = MAX(NNR1,NHG*2)*NLSX !GCENER and COREC
IF(TGC) THEN
IL_GRAD=NNR1*NLSD*4 !GRADEN
ELSE
IL_GRAD=0
ENDIF
IF(TINLC) THEN
IL_RHOVAL=NNR1*NLSD !XCENER
ELSE
IL_RHOVAL=0
ENDIF
CALL SCRPTR(SCR,LSCR,IP_VTMP,IL_VTMP,'VTMP',DEBUG)
CALL SCRPTR(SCR,LSCR,IP_GRAD,IL_GRAD,'GRAD',DEBUG)
CALL SCRPTR(SCR,LSCR,IP_RHOVAL,IL_RHOVAL,'RHOVAL',DEBUG)
C ==--------------------------------------------------------------==
C == ADD CORE CHARGE TO RHOE ==
C ==--------------------------------------------------------------==
IF(TINLC) THEN
CALL DCOPY(NNR1*NLSD,RHOE(1,1),1,RHOVAL(1,1),1)
CALL COREC(RHOE,VTMP,V)
CALL DCOPY(2*NHG,VTEMP(1,1),1,VNLT(1),1)
ENDIF
C ==-------------------------------------------------------------==
C == VTEMP (Potential in G-Space) -FFT-> V(R) ==
C ==--------------------------------------------------------------==
CALL AZZERO(V,2*MAXFFT)
#ifdef __NEC
!CDIR NODEP
#endif
DO IG=1,NHG
V(INDZ(IG),1) = DCONJG(VTEMP(IG,1))
V(NZH(IG),1) = VTEMP(IG,1)
ENDDO
IF(GEQ0) V(NZH(1),1) = VTEMP(1,1)
CALL INVFFT(V)
C ==--------------------------------------------------------------==
C == ADD EXTERNAL POTENTIAL TO V ==
C ==--------------------------------------------------------------==
IF(TEXTFLD) THEN
C..potential from classical interface
DO IR=1,NNR1
V(IR,1)=V(IR,1)+DCMPLX(EXTF(IR),0.D0)
ENDDO
ENDIF
IF(TEXPOT) THEN
C..static external potential
DO IR=1,NNR1
V(IR,1)=V(IR,1)+DCMPLX(EXTPOT(IR),0.D0)
ENDDO
ENDIF
IF((TEXPOT.OR.TEXTFLD).AND.TINLC) THEN
CALL FWFFT(V)
DO IG=1,NHG
VNLT(IG) = V(NZH(IG),1)
ENDDO
CALL INVFFT(V)
ENDIF
C ==--------------------------------------------------------------==
C == COMPUTE EXCHANGE AND CORRELATION ENERGY (EXC) ==
C ==--------------------------------------------------------------==
IF(TINLC) THEN
CALL XCENER(SXC,VXC,RHOVAL,RHOE,V)
ELSE
CALL XCENER(SXC,VXC,RHOE,RHOE,V)
ENDIF
SGCX = 0.0D0
SGCC = 0.0D0
VGC = 0.0D0
IF(TGC) THEN
C ==------------------------------------------------------------==
C == CALCULATE THE GRADIENT OF THE DENSITY ==
C ==------------------------------------------------------------==
IF(TLSD) THEN
IF(TSTRESS.OR.TDIAG) THEN
CALL DCOPY(2*NNR1,V(1,1),1,DQG(1),1)
CALL DCOPY(2*NNR1,V(1,2),1,DQG(NNR1+1),1)
ENDIF
CALL FWFFT(V(1,1))
CALL ZGTHR(NHG,V(1,1),VTEMP(1,1),NZH)
CALL FWFFT(V(1,2))
CALL ZGTHR(NHG,V(1,2),VTEMP(1,2),NZH)
DO IR=1,NNR1
RHOE(IR,1)=RHOE(IR,1)-RHOE(IR,2)
ENDDO
CALL GRADEN(RHOE(1,1),V,GRAD(1,1),VTMP)
CALL GRADEN(RHOE(1,2),V,GRAD(1,5),VTMP)
ELSE
IF(TSTRESS.OR.TDIAG) CALL DCOPY(2*NNR1,V(1,1),1,DQG(1),1)
CALL FWFFT(V)
CALL ZGTHR(NHG,V,VTEMP,NZH)
CALL GRADEN(RHOE,V,GRAD,VTMP)
ENDIF
C ==--------------------------------------------------------------==
C == GRADIENT CORRECTION TO THE EXCHANGE ENERGY (EGCX) ==
C ==--------------------------------------------------------------==
IF(TLSD) THEN
IF(TSTRESS) CALL STOPGM('VOFRHOB','NO STRESS WITH LSD AND GC')
CALL GCLSD(SGCX,SGCC,RHOE,V,VTEMP,VTMP,GRAD)
IF(TSTRESS.OR.TDIAG) THEN
CALL DAXPY(2*NNR1,-1.0D0,V(1,1),1,DQG(1),1)
CALL DAXPY(2*NNR1,-1.0D0,V(1,2),1,DQG(NNR1+1),1)
CALL DSCAL(4*NNR1,-1.0D0,DQG(1),1)
VGC=0.0D0
IF(TINLC) THEN
DO IR=1,NNR1
VGC=VGC+(RHOVAL(IR,1)-RHOVAL(IR,2))*DREAL(DQG(IR))
VGC=VGC+RHOVAL(IR,2)*DREAL(DQG(NNR1+IR))
ENDDO
ELSE
DO IR=1,NNR1
VGC=VGC+RHOE(IR,1)*DREAL(DQG(IR))
VGC=VGC+RHOE(IR,2)*DREAL(DQG(NNR1+IR))
ENDDO
ENDIF
ENDIF
ELSE
CALL GCENER(SGCX,SGCC,RHOE,V,VTEMP,VTMP,GRAD,TSTRESS)
IF(TSTRESS.OR.TDIAG) THEN
CALL DAXPY(2*NNR1,-1.0D0,V(1,1),1,DQG(1),1)
CALL DSCAL(2*NNR1,-1.0D0,DQG(1),1)
VGC=0.0D0
IF(TINLC) THEN
DO IR=1,NNR1
VGC=VGC+RHOVAL(IR,1)*DREAL(DQG(IR))
ENDDO
ELSE
DO IR=1,NNR1
VGC=VGC+RHOE(IR,1)*DREAL(DQG(IR))
ENDDO
ENDIF
ENDIF
ENDIF
ENDIF
C ==--------------------------------------------------------------==
C == XC CONTRIBUTION TO STRESS (VANDERBILT CHARGE) ==
C ==--------------------------------------------------------------==
IF(TIVAN.AND.TSTRESS)
* CALL XCSTR(DRHOVG,RHOE,DQG,DQG(MAXFFT+1))
C ==--------------------------------------------------------------==
C == NLCC CONTRIBUTION TO STRESS ==
C ==--------------------------------------------------------------==
IF(TINLC.AND.TSTRESS) CALL NLCCSTR(DQG,DQG(NLSD*MAXFFT+1))
C ==--------------------------------------------------------------==
C == RELEASE TEMPORARY ARRAYS FROM SCRATCH SPACE ==
C ==--------------------------------------------------------------==
CALL SCRREL(SCR,LSCR,DEBUG)
C ==--------------------------------------------------------------==
C == V CONTAINS THE TOTAL POTENTIAL IN R-SPACE ==
C == MOVE IT TO RHOE ==
C ==--------------------------------------------------------------==
DO IR=1,NNR1
RHOE(IR,1)=DREAL(V(IR,1))
ENDDO
IF(TLSD) THEN
DO IR=1,NNR1
RHOE(IR,2)=DREAL(V(IR,2))
ENDDO
ENDIF
C ==--------------------------------------------------------------==
C == CORRECT ENERGY FOR CORE CHARGE ==
C ==--------------------------------------------------------------==
IF(TINLC) THEN
SXC=SXC-EXCCO
SGCX=SGCX-EGCXCO
SGCC=SGCC-EGCCCO
ENDIF
IF(TIVAN.OR.(TINLC.AND.TFOR)) THEN
C ==--------------------------------------------------------------==
C == TRANSF. TOTAL POTENTIAL IN G-SPACE ==
C == PUT THE TOTAL POTENTIAL IN G-SPACE INTO VTEMP ==
C ==--------------------------------------------------------------==
CALL FWFFT(V)
CALL ZGTHR(NHG,V,VTEMP,NZH)
IF(TLSD) THEN
CALL FWFFT(V(1,2))
CALL ZGTHR(NHG,V(1,2),VTEMP(1,2),NZH)
ENDIF
ENDIF
IF(TIVAN) THEN
IF(IMAGP.EQ.2)CALL STOPGM('VOFRHOB','K-POINT NOT IMPLEMENTED')
C ==--------------------------------------------------------------==
C == CALCULATE DEEQ FOR VANDERBILT PP ==
C == FORCE ON IONS DUE TO THE "VANDERBILT CHARGE" ==
C ==--------------------------------------------------------------==
IF(TLSD) THEN
CALL NEWD(FNL(1,1,1,1),DEEQ(1,1,1,1),F(1),
* VTEMP(1,1),FION,SCR,LSCR,NSUP,TFOR)
CALL NEWD(FNL(1,1,1,NSUP+1),DEEQ(1,1,1,2),F(NSUP+1),
& VTEMP(1,2),FION,SCR,LSCR,NSDOWN,TFOR)
ELSE
CALL NEWD(FNL,DEEQ,F,VTEMP,FION,SCR,LSCR,N,TFOR)
ENDIF
ENDIF
C ==--------------------------------------------------------------==
C == FORCE ON IONS DUE TO NLCC ==
C ==--------------------------------------------------------------==
IF(TINLC.AND.TFOR) CALL COFOR(FION,VTEMP)
C ==--------------------------------------------------------------==
C == SINGLE PIECES OF THE ENERGY: ==
C ==--------------------------------------------------------------==
NNRS = NR1S*NR2S*NR3S
EGC=(SGCC+SGCX)*OMEGA/DFLOAT(NNRS)
EXC=SXC*OMEGA/DFLOAT(NNRS) + EGC
VXC=(VXC + VGC)*OMEGA/DFLOAT(NNRS)
IF(TSTRESS) THEN
DO KK=1,6
DEXC(KK) = DEXC(KK) + (EXC-VXC)*DELTA(ALPHA(KK),BETA(KK))
ENDDO
ENDIF
CALL TIHALT(' VOFRHOB',ISUB)
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
SUBROUTINE GIVE_SCR_VOFRHOB(LVOFRHOB,TAG)
C ==--------------------------------------------------------------==
IMPLICIT NONE
INCLUDE 'system.h' !NHG TGC
INCLUDE 'spin.inc' !NLSD NLSX
INCLUDE 'pslo.inc' !TIVAN
INCLUDE 'nlcc.inc' !TINLC
C Arguments
INTEGER LVOFRHOB
CHARACTER TAG*30
C Variables
INTEGER LTGC,LTINLC
C ==--------------------------------------------------------------==
IF(TGC) THEN
LTGC=NNR1*NLSD*4 !GRADEN
ELSE
LTGC=0
ENDIF
IF(TINLC) THEN !XCENER
LTINLC=NNR1*NLSD
ELSE
LTINLC=0
ENDIF
LVOFRHOB = MAX(NNR1,NHG*2)*NLSX+LTGC+LTINLC !GCENER and COREC
IF(TIVAN) LVOFRHOB=MAX(LVOFRHOB,NHG*4+NAX) !NEWD
LVOFRHOB = LVOFRHOB + 100 !For boundary checks in SCRPTR
TAG='MAX(NNR1,NHG*2)*NLSX+NNR1*NLSD*4'
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
More information about the CPMD-list
mailing list