[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