[Wien] setrmt_lapw #4

[John Rundgren]

Dear Professor Blaha,

Many thanks for the final remarks on the radial extension of vcoul. A
similar extension of clmsum is also needed (no spin). 

In lapw0.F a statement concerning clm is found, line 563,
    READ(8,2021) (CLM(J,LM1,JATOM), J=1,JRI(JATOM))
leading down to lines 1109:1141 where summations occur just as in the
Fourier extension of vcoul, namely, 
    DO J=2,NKK ...
    DO LM1=1,LLMM ...
    CQ(lm1)=CQ(lm1)+RHOKT*YKA(LM1)*BES(LL+1)

I should be most obliged for a piece of code in which CQ is transformed
into
    new_clm(j,lm1,jatom), j=1,jrx(jatom), jrx(jatom)>jri(jatom) 

Best regards,
John Rundgren


On Sun, 2012-09-30 at 20:35 +0200, Peter Blaha wrote:
> 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
> >
>