[Wien] mbj on Diamond
Kamil Klier
kk04 at Lehigh.EDU
Wed Nov 10 23:55:15 CET 2010
Dear Prof. Blaha,
Your attached brj.j triggers the following error message after a crash
in our CFN-BNL environment:
*******************************************************************************************************
LAPW0 END
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
lapw0 000000000040A79E brj_ 37 brj.f
lapw0 000000000047338D vxclm2_ 534 vxclm2.f
lapw0 000000000047B565 xcpot1_ 365 xcpot1.f
lapw0 00000000004393FA MAIN__ 1696 lapw0.F
lapw0 0000000000403A1C Unknown Unknown Unknown
libc.so.6 000000370E01D994 Unknown Unknown Unknown
lapw0 0000000000403929 Unknown Unknown Unknown
********************************************************************************************************
However, the file brj.j_current, which I am now attaching, does not
crash in the same runs, although the covergence is very slow.
This prompts the questions:
Q1: What's wrong with the attached brj.f_current ?
Q2: For my curiosity, what does the parameter "ir" do ? [I have
attempted to trace it in vxclm2.f,also attached, where it is
referenced, but the calls to BRJ do not include ",ir" as you suggest]
The runs involved are 128 atom supercells of ZnO with various dopants.
The 'regular scf' converges well, the slow convergence occurs in the
mBJ runs with option 28 in *in0 and 50 in *in0_grr, even though the
PRATT option in *inm is used.
I believe I have the latest version of WIEN2k, 10.2.
Regards,
Kamil Klier
Quoting Peter Blaha <pblaha at theochem.tuwien.ac.at>:
> Hi,
> No, use the attached one.
>
> Note: it has one more parameter (ir), thus the calling in vxclm2.f
> should also be modified and ",ir" should be added to all calls.
>
> Best regards
>
> Am 06.11.2010 14:59, schrieb Kamil Klier:
>> 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.
>>
>>
>>
>> _______________________________________________
>> Wien mailing list
>> Wien at zeus.theochem.tuwien.ac.at
>> http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
>
> --
>
> 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)
!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
-------------- next part --------------
! potential option 28: modified Becke-Johnson for the exchange potential and
! LDA for the correlation potential.
! LDA for the exchange-correlation energy functional.
!
else IF (IGRAD .EQ. 28) THEN
GDGGu=GXU*GGXU+GYU*GGYU+GZU*GGZU
GDGGd=GXD*GGXD+GYD*GGYD+GZD*GGZD
T=SQRT((GXU+GXD)**2+(GYU+GYD)**2+(GZU+GZD)**2)
call pbe2(fu,gmagu,gdggu,g2u,fd,gmagd,gdggd,g2d, &
t,gdggt,1,1, &
exupls,exdnls,vxupls,vxdnls,eclsd,vcupls,vcdnls, &
exuppw,exdnpw,vxuppw,vxdnpw,ecpw,vcuppw,vcdnpw, &
exuppb,exdnpb,vxuppb,vxdnpb,ecpbe,vcuppb,vcdnpb)
exu=2.d0*(exupls + eclsd)
exd=2.d0*(exdnls + eclsd)
tauup=(tauup+2.d0*g2u*0.25d0-vxuold*fu)*0.5d0
taudn=(taudn+2.d0*g2d*0.25d0-vxdold*fd)*0.5d0
XCCONST = XCA + XCB*GRR**XCP
CALL BRJ(FU,GMAGU,G2U,2D0*TAUUP,VXBRJU)
VXU = 2D0*(VXBRJU + VCUPLS)
CALL BRJ(FD,GMAGD,G2D,2D0*TAUDN,VXBRJD)
VXD = 2D0*(VXBRJD + VCDNLS)
! if(ir.ne.index.and.ir.gt.770.and.ir.lt.777) then
! write(6,'(i11,6e13.5)') ir, vxmix,grr,xcconst,fu,vxuold,vxu
! endif
! index=ir
VXU = vxmix*VXU + (1.d0-vxmix)*VXUOLD
VXD = vxmix*VXD + (1.d0-vxmix)*VXDOLD
if(fu.lt.tlev) then
exu=2.d0*(exupls+eclsd)
vxu=2.d0*(vxupls+vcupls)
endif
if(fd.lt.tlev) then
exd=2.d0*(exdnls+eclsd)
vxd=2.d0*(vxdnls+vcdnls)
endif
RETURN
More information about the Wien
mailing list