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