[Wien] setrmt_lapw #5

[Peter Blaha]

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
>

-- 

                                       P.Blaha
--------------------------------------------------------------------------
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-165300             FAX: +43-1-58801-165982
Email: blaha at theochem.tuwien.ac.at    WWW: http://info.tuwien.ac.at/theochem/
--------------------------------------------------------------------------