[Wien] setrmt_lapw #5

[John Rundgren]

Dear Professor Blaha,

I should be most grateful if you would criticize the following code
extending vcoul beyond RMT.

Attached is TiO2.vcoul for LM=0,0 corresponding to rmt=1.92,1.74
(setrmt) and rmt=2.0,1.6 (UG example). Encouraging observations:
  vcoul is smooth and differentiable;
  vcoul is almost unaffected by rmt change;
  vcoul max. value is similar for Ti and O (the same interstice).
But there can still be misunderstandings in the suggested code.
Modifications marked !JR or !PB concern all LM's:

1) Extended radii rxt > RMT and corresponding max. subscripts:
    rext=(/2.6d0,2.2d0/)
    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

2) The part of code where exp(iKr) is expanded in Bessel functions:
    CVOUT=(0.d0,0.d0)
    allocate(cvoutx(jri(jatom):jrx(jatom),ncom+3)) !JR
    CVOUTx=(0.d0,0.d0)                             !JR

    do jj=jri(jatom),jrx(jatom)                      !JR
      r(jj)=r0(jatom)*exp(dx(jatom)*dble(jj-1))      !JR
        
      DO J=2,NKK
        ARG=r(jj)*ABSK(J)                            !PB 
        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
       CVOUTx(jj,LM1)=CVOUTx(jj,LM1) + POTK(J)*BES(L)*YKA(LMMULT1) !JR
       if(jj==jri(jatom)) CVOUT(LM1)=CVOUTx(jj,LM1)                !JR
        ENDDO !LM1

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

      ENDDO !J
    enddo !jj

      IF(SWITCH.EQ.'COUL') then
        write(3,'("ATOM",i4,/,(2f15.8))') &
          jatom,( rhopw0r(j),sqfp*rhopw0(j),j=jrj-1,jrj+1)
      endif

      DO LM1=1,LLMM
        M=LM(2,LM1,JATOM)
        IF(M.NE.0) THEN
          IMAG1=(1.d0,0.d0)
          IF(LM(1,LM1,JATOM).LT.0) IMAG1=-IMAG
          IF(MOD(M,2).EQ.1) IMAG1=-IMAG1
          do jj=jri(jatom),jrx(jatom)                 !JR
            CVOUTx(jj,lm1)=CVOUTx(jj,lm1)/IMAG1*SQRT2 !JR
          enddo                                       !JR
          cvout(lm1)=cvoutx(jri(jatom),lm1)           !JR
        ENDIF
      ENDDO

4) Statements connected with SURFIN:
      DO J=1,jrx(jatom)            !JR
        RLDM1(J)=R(J)**L
      ENDDO
      DO J=1,JRJ
        RLDM2(J)=( 1.0D0-(R(J)/RMT(JATOM))**(2*L+1))/R(J)**(L+1)
      ENDDO
      rldm2(jrj+1:jrx(jatom))=0.d0 !JR

      DO 53 IRAD=1,jrx(jatom)      !JR
        rj1=rldm2(irad)
        DO J=1,irad-1

VALUE(J)=(CLM(J,LM1,JATOM)-delcc*r(j)*r(j)*pi4*LMOD)*RLDM1(J)*rj1
        ENDDO
        rj2=RLDM1(IRAD)
        DO J=irad,jrx(jatom)       !JR

VALUE(J)=(CLM(J,LM1,JATOM)-delcc*r(j)*r(j)*pi4*LMOD)*RJ2*RLDM2(J)
       ENDDO
!
!     integration in two parts, simpson + 3/8 formula at start(end)
!     for even points
!
       CALL CHARG2 (R,DX(JATOM),VALUE,1,IRAD,VLIN1)
       CALL CHARG3 (R,DX(JATOM),VALUE,IRAD,jrx(jatom),VLIN2) !JR
       ACORR=VALUE(1)*R(1)*0.5D0
       VOLINT=VLIN1+VLIN2+ACORR
       VOLINT =VOLINT*8.0*PI/dble(2*L+1)*VINF
       if(r(irad)<=rmt(jatom))then                             !JR
         SURFIN=CVOUT(LM1)*(-4.D0*PI)*(R(IRAD)/RMT(JATOM))**L* &
                IMAG**L*VOUTF
         ZSHIFT=(-2.D0*ZZ(JATOM)/R(IRAD)+ &
                2.D0*ZZ(JATOM)/RMT(JATOM))*SQFP
         V(IRAD,LM1,JATOM)=(VOLINT-SURFIN+LMOD*ZSHIFT)
       elseif(r(irad)>rmt(jatom))then                          !JR
         SURFIN=cvoutx(irad,LM1)*(-4.D0*PI)*IMAG**L*VOUTF      !PB
         zshift=0.d0                                           !JR
       endif                                                   !JR
       V(IRAD,LM1,JATOM)=(VOLINT-SURFIN+LMOD*ZSHIFT)

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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: VC.ps
Type: application/postscript
Size: 21309 bytes
Desc: not available
URL: <http://zeus.theochem.tuwien.ac.at/pipermail/wien/attachments/20121004/a71ecfb0/attachment-0001.ps>