[Wien] setrmt_lapw #4

[John Rundgren]

Dear Professor Blaha,

I should be most grateful if you would criticize the following attempt
to extend vcoul beyond RMT. Below are my modifications following your
advice of 09/14/2012. TiO2 in UG is used for testing.

1) External radii > RMT are defined:
      rext=(/2.6d0,2.2d0/) !for TiO2.

      do JATOM=NSTART(myid),NSTOP(myid)
        do J=1,jri(jatom)+100
          R(J)=R0(JATOM)*EXP(dble(J-1)*DX(JATOM))
          if(r(j)>rext(jatom))then
            jrx(jatom)=j-1
            exit
          endif
        enddo
      enddo
      !cvout depends on jrx(jatom) below.
      allocate(cvout(maxval(jrx),ncom+3))

2) The part of code where exp(iKr) is expanded in Bessel functions.
      CVOUT=(0.d0,0.d0)

      do jj=jri(jatom),jrx(jatom)
        r(jj)=r0(jatom)*exp(dx(jatom)*dble(jj-1))
        
        DO J=2,NKK
          ARG=r(jj)*ABSK(J) !RMT(JATOM) replaced by r(jj).
          CALL SPHBES(lmax2+1,ARG,BES)
          CALL YKAV(J,JATOM,YKA,LMMTMX,LMMULT,LLMM)
          LMMULT1=0

          DO LM1=1,LLMM
            L=IABS( LM(1,LM1,JATOM) )
            M=LM(2,LM1,JATOM)
            LMMULT1=LMMULT1+1
            IF (LMMULT1.NE.1) THEN
              IF ((L.EQ.IABS(LMMULT(1,LMMULT1-1,JATOM))).AND. &
                (M.EQ.LMMULT(2,LMMULT1-1,JATOM))) THEN
                LMMULT1=LMMULT1-1
              endif
            endif
            CVOUT(jj,LM1)=CVOUT(jj,LM1) + POTK(J)*BES(L)*YKA(LMMULT1)
          ENDDO !LM1

          IF(SWITCH.EQ.'COUL')  call
rhopw(rhok(j),yka,ABSK(J),rhopw0r,jrj,rhopw0)

        ENDDO !J
      enddo !jj
      jrj=jrx(jatom)

      !cvout = constant in interval [R0,RMT] as in original code.
      do lm1=1,llmm
        do jj=1,jri(jatom)-1
          cvout(jj,lm1)=cvout(jri(jatom),lm1)
        enddo
      enddo

3) Statement
            CVOUT(lm1)=CVOUT(lm1)/IMAG1*SQRT2
   is replaced by
            CVOUT(1:jrx(jatom),lm1)=CVOUT(1:jrx(jatom),lm1)/IMAG1*SQRT2

4) Statement SURFIN=... is modified. Is **L due to a Legendre expansion?
      if(r(irad)<=rmt(jatom))then
        SURFIN=CVOUT(irad,LM1)*(-4.D0*PI)*(R(IRAD)/RMT(JATOM))**L* &
          IMAG**L*VOUTF
      elseif(r(irad)>rmt(jatom))then
        SURFIN=CVOUT(irad,LM1)*(-4.D0*PI)*(RMT(JATOM)/R(IRAD))**L* &
          IMAG**L*VOUTF
      endif

The modifications are meant to be valid for all LM's. Please, see the
Attachment depicting TiO2_R.vcoul in Rydberg units, LM=0,0. The test
indicates that "run SCF" is not perturbed by the modifications made.

I wonder what your decision will be about the modifications.
Best regards,
John Rundgren

-------------- next part --------------
A non-text attachment was scrubbed...
Name: TiO2_R.vcoul.ps
Type: application/postscript
Size: 17823 bytes
Desc: not available
URL: <http://zeus.theochem.tuwien.ac.at/pipermail/wien/attachments/20120930/981fc526/attachment.ps>