[Wien] setrmt_lapw #4

[Peter Blaha]

I can see that now the potential is at least continuous, but the slopes
still do not match.

I don't have the equations with me at home, but as far as I can see:

I hope, you have redimensioned the arrays    r, jrx  properly.

 >        elseif(r(irad)>rmt(jatom))then
 >          SURFIN=CVOUT(irad,LM1)*(-4.D0*PI)*(RMT(JATOM)/R(IRAD))**L* &
 >            IMAG**L*VOUTF
 >        endif

I don't think this is ok. The modification of cvout by r(i)/r(RMT) is valid only
inside the spheres, where the coulomb potential depends on on both, the
charge inside the sphere and the (scaled) potential at RMT.

Outside the sphere, we only need the bessel expansion of the plane waves, i.e.

           SURFIN=CVOUT(irad,LM1)*(-4.D0*PI)*IMAG**L*VOUTF
----------------------------------------------------------------------------------
PS: We want the coefficients V_LM(r) of V(r) = SUM_LM (V_LM(r) Y_LM(r-hat))   (1)

In the interstitial the potential is given as:
                              V(r) = 1/vol sum_K (c_K exp (iKr) )        (2)


We expand the exp (iKr) = 4 pi sum_LM (i**L j_L(Kr) Y_LM(K) * Y_LM(r)    (3)

we set equ (1) = equ (2); put (3) into (2); and multiply this equality by Y_LM(r)
and integrate over the angles.
Then from (1) only one LM remains, i.e. only V_LM(r) remains (thats the coefficient you are looking for) and

V_LM(r) = 4 pi/vol sum_K [ c_K  i**L j_L(Kr) Y_LM(K) ]

This is (besides the constants)    CVOUT(jj,LM1)=CVOUT(jj,LM1) + POTK(J)*BES(L)*YKA(LMMULT1)

Depending how you later on "use" the calculated quantities, it may need an additional
r**2, or a factor of 1/sqrt(4pi), which is equal to Y_00 ! or similar....

and this is probably added in the SURFIN line.

Best regards
Peter


Am 30.09.2012 16:03, schrieb 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
>
>
>
> _______________________________________________
> Wien mailing list
> Wien at zeus.theochem.tuwien.ac.at
> http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
>

-- 
-----------------------------------------
Peter Blaha
Inst. Materials Chemistry, TU Vienna
Getreidemarkt 9, A-1060 Vienna, Austria
Tel: +43-1-5880115671
Fax: +43-1-5880115698
email: pblaha at theochem.tuwien.ac.at
-----------------------------------------