[Wien] mbj on Diamond
Kamil Klier
kk04 at Lehigh.EDU
Sat Nov 6 14:59:45 CET 2010
Dear Prof. Blaha,
Is the brj.f file referred to in this e-mail (attached as brj.f_new)
one that should replace any previous brj.f files (for example,
attached as brj.f_old)?
Regards,
Kamil Klier
Quoting Peter Blaha <pblaha at theochem.tuwien.ac.at>:
> Hi,
> After I got the struct file, I could verify the problem.
>
> As expected, it is again in the SRC_lapw0/brj.f subroutine,
> where one has to find the root of a function.
>
> For strange densities this is numerically non-trivial.
>
> The problem at the nucleus of heavy elements was solved before, but
> here is the problem in the interstitial, when rho is very small and
> also grad-rho, tau and laplace-rho are sufficiently small.
>
> Then the required functions are "nearly" zero (lt. 10^-6) for a range
> of x=10-30; but using x=30 produces a V-xc potential of -100 Ry,
> which is the reason for your "Eigenvalues below zero".
>
> When such problems occur again, please check also case.output0.
> The Fourier coefficients of Vxc must "converge", i.e. (0 0 0) should be
> order one, while (0 0 30) should be order 10^-5 .
>
> The attached subroutine brj.f should fix these problems (at least
> your case converges
> smoothly).
>
> Dear All,
>
>
>
> We are performing the mbj calculations for a carbon based compound.
> According to the usersguide there are three SCF cycles for mbj
> calculations: first a regular calculations within LDA/GGA (we use
> the PBE-GGA here), second one more cycle run_lapw ?NI ?i 1 , and
> third the mbj run after changing the potential energy functional
> indxc=28 in case.in0 and index=50 in case.in0_grr.
>
> Here we call the regular SCF cycles C1.scf, second one-more SCF
> cycle as C2.scf, and the third the mbj as cycle C3.scf.
>
> The first regular cycle and the second run_lapw ?NI ?i 1 are
> converged smoothly. However, the third mbj cycle is stopped at lapw2
> in its second iteration.
>
> We analyzed the problem to find the source of the error. The result
> is given below, where the C2.scf line refers to the last :ITE of the
> second one more SCF cycle, and the C3.scf refers to the first :ITE
> of the third mbj run:
>
>
>
> C2.scf::NTO033: TOTAL CHARGE IN SPHERE 1 = 3.9781366
>
> C3.scf::NTO033: TOTAL CHARGE IN SPHERE 1 = 2.4250427
>
>
>
> C2.scf::CTO033: TOTAL CHARGE IN SPHERE 1 = 3.9781254
>
> C3.scf::CTO033: TOTAL CHARGE IN SPHERE 1 = 3.9470631
>
>
>
> C2.scf::DIS : CHARGE DISTANCE ( 0.0000355 for atom 33 spin
> 1) 0.0000136
>
> C3.scf::DIS : CHARGE DISTANCE ( 1.8978668 for atom 25 spin
> 1) 1.5016586
>
>
>
> C2.scf::NEC01: NUCLEAR AND ELECTRONIC CHARGE 366.00000
> 365.98257 1.00005
>
> C3.scf::NEC01: NUCLEAR AND ELECTRONIC CHARGE 366.00000
> 365.98171 1.00005
>
>
>
> C2.scf::FER : F E R M I - ENERGY(TETRAH.M.)= 0.21390
>
> C3.scf::FER : F E R M I - ENERGY(TETRAH.M.)= -1.44751
>
>
>
> The result clearly shows that there is a jump in :NTO, :DIS, and
> :FER (But in :CTO) after changing the functional to index=28.
>
>
>
> --
>
> P.Blaha
> --------------------------------------------------------------------------
> Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
> Phone: +43-1-58801-15671 FAX: +43-1-58801-15698
> Email: blaha at theochem.tuwien.ac.at WWW: http://info.tuwien.ac.at/theochem/
> --------------------------------------------------------------------------
>
----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.
-------------- next part --------------
SUBROUTINE BRJ(RHO,GRHO,G2RHO,TAU,VXBRJ,ir)
!A. D. Becke and M. R. Roussel, Phys. Rev. A 39, 3761 (1989).
!A. D. Becke and E. R. Johnson, J. Chem. Phys. 124, 221101 (2006).
USE xcparam
IMPLICIT REAL*8(A-H,O-Z)
SAVE X,iint,isphere
DATA X/0D0/,iint/0/,isphere/0/
PI = 4D0*DATAN(1D0)
IF (RHO .GT. 1D-18) THEN
TOL = 1D-6
F = 10D0*TOL
NLOOP = 0
X1 = 0D0
X2 = 0D0
NIF = 0
DABSFOLD = -1D0
tautf = (3d0/10d0)*(3d0*pi**2)**(2d0/3d0)*(2d0*rho)**(5d0/3d0)
tauw = 0.125d0*grho*grho*2.d0/rho
if(tau.lt.tauw) then
tau_falsch=tau
tau=tauw
endif
if(tau.gt.(tautf+tauw)) then
tau_falsch=tau
tau=tautf+tauw
endif
if(tau.eq.tauw .and. rho.lt.10.d0.and.ir.lt.900.and.isphere.eq.0) then
print*,'sphere:rho,tauw,grho,g2rho',rho,tau,grho,g2rho,'tauwrong=',tau_falsch
isphere=1
endif
if(tau.eq.(tautf+tauw) .and. rho.lt.10.d0.and.ir.lt.900.and.isphere.eq.0) then
print*,'sphere:rho,tauw+tautf,grho,g2rho',rho,tau,grho,g2rho,'tauwrong=',tau_falsch
isphere=1
endif
if(tau.eq.tauw .and. ir.gt.900.and.iint.lt.10) then
print*,'int:rho,tauw,grho,g2rho',rho,tau,grho,g2rho,'tauwrong=',tau_falsch
iint=iint+1
endif
if(tau.eq.(tautf+tauw) .and. ir.gt.900.and.iint.lt.10) then
print*,'int:rho,tauw+tautf,grho,g2rho',rho,tau,grho,g2rho,'tauwrong=',tau_falsch
iint=iint+1
endif
D = TAU - 0.25D0*GRHO**2D0/RHO
Q = (1D0/6D0)*(G2RHO - 2D0*0.8D0*D)
10 DO WHILE (DABS(F) .GE. TOL)
F = X*DEXP(-2D0*X/3D0)*Q - &
(2D0/3D0)*PI**(2D0/3D0)*(X-2D0)*RHO**(5D0/3D0)
DF = (1D0-(2D0/3D0)*X)*DEXP(-2D0*X/3D0)*Q - &
(2D0/3D0)*PI**(2D0/3D0)*RHO**(5D0/3D0)
IF (DABS(DF) .GT. 1D-18) THEN
X = X - F/DF
ELSE
X = X - F/SIGN(1D-18,DF)
ENDIF
IF ((NLOOP .GE. 500) .OR. (DABS(X) .GE. 1D20) .OR. (DABS(F) .GE. 1D20)) THEN
IF (DABS(DABS(F)-DABSFOLD) .LT. 1D-8) THEN
NIF = NIF + 1
ENDIF
DABSFOLD = DABS(F)
IF (NIF .EQ. 5) THEN
TOL = 5D0*TOL
NIF = 0
ENDIF
X = X1
X1 = X1 + 1D0
F = 10D0*TOL
NLOOP = 0
GOTO 10
ENDIF
NLOOP = NLOOP + 1
ENDDO
if(dabs(df) .le. tol) then
tol=tol/10.d0
if(tol.gt.1.d-12) goto 10
endif
IF (X .LT. 0D0) THEN
X = X2
X2 = X2 + 1D0
F = 10D0*TOL
NLOOP = 0
GOTO 10
ENDIF
B = (X**3D0*DEXP(-X)/(8D0*PI*RHO))**(1D0/3D0)
IF (B .GT. 1D-18) THEN
VXBRJ = -(1D0 - DEXP(-X) - 0.5D0*X*DEXP(-X))/B
! vxbrj1=vxbrj
IF (TAU .GE. 0D0) THEN
VXBRJ = XCCONST*VXBRJ + (3D0*XCCONST-2D0)*DSQRT(5D0/12D0)/PI*DSQRT(TAU/RHO)
ELSE
VXBRJ = XCCONST*VXBRJ - (3D0*XCCONST-2D0)*DSQRT(5D0/12D0)/PI*DSQRT(ABS(TAU/RHO))
ENDIF
ELSE
VXBRJ = 0D0
ENDIF
ELSE
VXBRJ = 0D0
ENDIF
! if(ir.gt.900.and.iint.lt.100.and.abs(Vxbrj).gt.5.d0) print*,'brj:rho,tau,tauw,grho,g2rho',rho,tau,tauw,grho,g2rho,'VXBRJ',Vxbrj,d,q,b,x,DSQRT(TAU/RHO),vxbrj1,ir
RETURN
END
-------------- next part --------------
SUBROUTINE BRJ(RHO,GRHO,G2RHO,TAU,VXBRJ)
!A. D. Becke and M. R. Roussel, Phys. Rev. A 39, 3761 (1989).
!A. D. Becke and E. R. Johnson, J. Chem. Phys. 124, 221101 (2006).
USE xcparam
IMPLICIT REAL*8(A-H,O-Z)
SAVE X
DATA X/0D0/
PI = 4D0*DATAN(1D0)
IF (RHO .GT. 1D-18) THEN
TOL = 1D-4
F = 10D0*TOL
NLOOP = 0
X1 = 0D0
X2 = 0D0
tauw = 0.125d0*grho*grho*2.d0/rho
if(tau.lt.tauw) then
tau_falsch=tau
tau=tauw
endif
D = TAU - 0.25D0*GRHO**2D0/RHO
Q = (1D0/6D0)*(G2RHO - 2D0*0.8D0*D)
if(tau.eq.tauw .and. q.lt.-1.d9) then
q1=q
!print*,'rho,x,q,tau,grho,g2rho',rho,x,q1,tau,grho,g2rho,'tf=',tau_falsch
q=-1.d9
endif
10 DO WHILE (DABS(F) .GE. TOL)
F = X*DEXP(-2D0*X/3D0)*Q - &
(2D0/3D0)*PI**(2D0/3D0)*(X-2D0)*RHO**(5D0/3D0)
DF = (1D0-(2D0/3D0)*X)*DEXP(-2D0*X/3D0)*Q - &
(2D0/3D0)*PI**(2D0/3D0)*RHO**(5D0/3D0)
IF (DABS(DF) .GT. 1D-18) THEN
X = X - F/DF
ELSE
X = X - F/SIGN(1D-18,DF)
ENDIF
!print*,'f,df,x,nloop', f,df,x,nloop
IF ((NLOOP .GE. 400) .OR. (DABS(X) .GE. 1D10) .OR. (DABS(F) .GE. 1D10)) THEN
X = X1
X1 = X1 + 1D0
F = 10D0*TOL
NLOOP = 0
GOTO 10
ENDIF
NLOOP = NLOOP + 1
ENDDO
IF (X .LT. 0D0) THEN
X = X2
X2 = X2 + 1D0
F = 10D0*TOL
NLOOP = 0
GOTO 10
ENDIF
B = (X**3D0*DEXP(-X)/(8D0*PI*RHO))**(1D0/3D0)
IF (B .GT. 1D-18) THEN
VXBRJ = -(1D0 - DEXP(-X) - 0.5D0*X*DEXP(-X))/B
IF (TAU .GE. 0D0) THEN
VXBRJ = XCCONST*VXBRJ + (3D0*XCCONST-2D0)*DSQRT(5D0/12D0)/PI*DSQRT(TAU/RHO)
ELSE
VXBRJ = XCCONST*VXBRJ - (3D0*XCCONST-2D0)*DSQRT(5D0/12D0)/PI*DSQRT(ABS(TAU/RHO))
ENDIF
ELSE
VXBRJ = 0D0
ENDIF
ELSE
VXBRJ = 0D0
ENDIF
RETURN
END
More information about the Wien
mailing list