[Wien] setrmt_lapw #5
Gavin Abo
gsabo at crimson.ua.edu
Tue Oct 9 00:21:13 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?
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:
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
>>
>
More information about the Wien
mailing list