Hi,
Most likely the problem is due to the Newton iteration procedure in brj.f
which does not stop. Replace this subroutine by the one that I attached
(it is based on an analytical representation of the solution).
By the way, it is nonsense to calculate the average of grad rho/rho when
there is vacuum. So, fix the value of c to the one obtained for the bulk.
F. Tran
On Fri, 18 Oct 2013, alpa dashora wrote:
Dear Prof. Blaha and Wien2k users,
I am trying to run mbj calculation for MoS2 10 layer slab. After few cycles
lapw0 hanged. I have seen the earlier posts discussed on the same topics. I
have used the vxclm2.f file as
provided by Prof. Blaha and removed case.in0_grr file. I have also increased
the Rmt*Kmax up to 9 but still I am not able to solve my problem.
Please suggest any other solution for this error or any other method to
calculate the correct band gap.
Thanks in advance.
With kind regards,
--
Alpa Dashora
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).
!E. Proynov, Z. Gan, and J. Kong, Chem. Phys. Lett. 455, 103 (2008).
use xcparam
implicit real*8(a-h,o-z)
real*8 :: a(1:3), b(0:5), c(0:5), d(0:5), e(0:5), yp(0:5)
save iint,isphere
data iint/0/,isphere/0/
pi = 4d0*atan(1d0)
vxbrj = 0d0
if (rho .gt. 1d-18) then
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.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.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
a(1) = 1.5255251812009530d0
a(2) = 0.4576575543602858d0
a(3) = 0.4292036732051034d0
c(0) = 0.7566445420735584d0
c(1) = -2.6363977871370960d0
c(2) = 5.4745159964232880d0
c(3) = -12.657308127108290d0
c(4) = 4.1250584725121360d0
c(5) = -30.425133957163840d0
b(0) = 0.4771976183772063d0
b(1) = -1.7799813494556270d0
b(2) = 3.8433841862302150d0
b(3) = -9.5912050880518490d0
b(4) = 2.1730180285916720d0
b(5) = -30.425133851603660d0
d(0) = 0.00004435009886795587d0
d(1) = 0.58128653604457910d0
d(2) = 66.742764515940610d0
d(3) = 434.26780897229770d0
d(4) = 824.7765766052239000d0
d(5) = 1657.9652731582120d0
e(0) = 0.00003347285060926091d0
e(1) = 0.47917931023971350d0
e(2) = 62.392268338574240d0
e(3) = 463.14816427938120d0
e(4) = 785.2360350104029000d0
e(5) = 1657.962968223273000000d0
dd = tau - 0.25d0*grho**2/rho
q = (g2rho - 1.6d0*dd)/6d0
if (abs(q) .gt. 1d-18) then
y = (2d0/3d0)*pi**(2d0/3d0)*rho**(5d0/3d0)/q
do i=0, 5
yp(i) = y**i
enddo
if (y .le. 0d0) then
g = -atan(a(1)*y + a(2)) + a(3)
p1 = sum(c(0:5)*yp(0:5))
p2 = sum(b(0:5)*yp(0:5))
elseif (y .gt. 0d0) then
z = 2.085749716493756d0*y
g = log(sqrt(1d0 + 1d0/z**2) + 1d0/z) + 2d0
p1 = sum(d(0:5)*yp(0:5))
p2 = sum(e(0:5)*yp(0:5))
endif
if (abs(p2) .gt. 1d-18) then
x = g*p1/p2
if (abs(x) .gt. 1d-18) then
vxbrj = -2d0*pi**(1d0/3d0)*rho**(1d0/3d0)*exp(x/3d0)/x*(1d0 - exp(-x) - 0.5d0*x*exp(-x))
if (tau .ge. 0d0) then
vxbrj = xcconst*vxbrj + (3d0*xcconst-2d0)*sqrt(5d0/12d0)/pi*sqrt(tau/rho)
else
vxbrj = xcconst*vxbrj - (3d0*xcconst-2d0)*sqrt(5d0/12d0)/pi*sqrt(abs(tau/rho))
endif
endif
endif
endif
endif
return
end
_______________________________________________
Wien mailing list
Wien@zeus.theochem.tuwien.ac.at
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at:
http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
