Please find enclosed brj.f file.

Sincerely
------------------------------------------------
Dr. K. C. Bhamu
(UGC-Dr. D. S. Kothari Postdoc Fellow)
Department of Physics
Goa University, Goa-403 206
India
Mob. No.  +91-9975238952

On Tue, Oct 25, 2016 at 5:04 PM, <t...@theochem.tuwien.ac.at> wrote:

> In this case I don't understand why lapw0 is hanging. Can you
> send the subroutine brj.f?
>
> On Tuesday 2016-10-25 12:18, Dr. K. C. Bhamu wrote:
>
> Date: Tue, 25 Oct 2016 12:18:41
>> From: Dr. K. C. Bhamu <kcbham...@gmail.com>
>> Reply-To: A Mailing list for WIEN2k users <w...@zeus.theochem.tuwien.ac.
>> at>
>> To: A Mailing list for WIEN2k users <wien@zeus.theochem.tuwien.ac.at>
>> Subject: Re: [Wien] regarding hanging of job and mBJ
>>
>>
>>      Probably your are using an old version of WIEN2k with the old
>>      scheme
>>      to find the solution of the nonlinear equation in brj.f. You
>>      should get
>>      the new version brj.f which avoids such problems.
>>
>>
>> No, I am using latest version of Wien2k with mkl+ifort.
>> It is not for all case. It happens for few cases and it doesn't matter
>> that
>> which case I am running.
>>
>>      _______________________________________________
>>      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
>>
>>
>>
>>
> _______________________________________________
> 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
>
>
      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

Reply via email to