[Wien] setrmt_lapw #5

Peter Blaha pblaha at theochem.tuwien.ac.at
Tue Oct 9 07:47:39 CEST 2012


> I could be wrong, but something doesn't seem right.  Can someone help clarify how the equations relate to the code?
>
> It was previously given that:
>
> V(r) = 1/vol sum_K (c_K exp (iKr) )        (2)
> exp (iKr) = 4 pi sum_LM (i**L j_L(Kr) Y_LM(K) * Y_LM(r)    (3)
>
> Combining (2) and (3),
>
> V(r) = 4 pi/vol i**L sum_K {sum_LM [c_K j_L(Kr) Y_LM(K) * Y_LM(r)] } (4)
>
> A part of the code is:
>
> CVOUTx(jj,LM1)=CVOUTx(jj,LM1) + POTK(J)*BES(L)*YKA(LMMULT1)
>
> Is c_K, j_L(Kr), and Y_LM(K) in the equation the code YKA(LMMULT1), BES(L), and POTK(J)? Where is Y_LM(r) in the code?

We are calculating the expansion coefficients V_LM  of     V(r)= sum_LM (V_LM(r) Y_LM(r)

Put this in the left side of equ. (4), multiply both sides with Y_L'M'(r) and integrate over angles.
Then Y_LM(r) vanishes and Y_L'M'(r) remains.

Of course, to get the spherically averaged potential, you have to multiply with 1/sqrt(4pi)=Y_00



>
> Another part of the code is:
>
>         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)
>
> Regarding the code, VOLINT and VOUTF, where is this in the equation?  I believe V(IRAD,LM1,JATOM) is the potential V(r) or is it just the coefficient V_LM(r)?  The code
> 1/vol seems to be missing in the SURFIN line.  However, in lapw0:

VOLINT should be zero for V(r>rmt) and VOUTF is some normalization factor (may contain volume?, but in particular
some constants for cubic harmonics,...)

Of course, the proper addition of some constants, Y_OO factors or r**L factors is the "art" and has to be checked carefully.
But the latest plots look absolutely correct.

>
> POTK(J)=(TMP*POTK(J)*INST(J)-RHOK(J))*TMP1/(ABSK(J)*ABSK(J))
>
> where
>
>        TMP=4.D0*PI/VOL
>        TMP1=-8.D0*PI
>
> Is the 4 pi/vol needed in SURFIN or is it already included in POTK(J)?
>
> The "CVOUTx=(0.d0,0.d0)" may need to be put inside the "do jj=jri(jatom),jrx(jatom)" so that sum starts at zero for each R.
>
> Kind regards.
>
> On 10/8/2012 2:09 AM, Peter Blaha wrote:
>> The ps file looks now as one would expect it and also the code looks reasonable,
>> although this is difficult to judge from a quick glance without running it.
>>
>> But when the plots look as they do now, I'd say this is ok.
>>
>> Am 04.10.2012 14:06, schrieb 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
>>>
>>>
>>>
>>> _______________________________________________
>>> Wien mailing list
>>> Wien at zeus.theochem.tuwien.ac.at
>>> http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
>>>
>>
>
> _______________________________________________
> 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
-----------------------------------------


More information about the Wien mailing list