Sample program list
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE PSEUDO(TOTCH,ETOT1)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DIMENSION VD(KIZA),VDNL(KIZA)
C DATA FOR ATOMS
TOTCH = 0.D0
ETOT1 = 0.D0
UNIRV = 1.0D0/UNIVOL
RG(1) =0.0D0
RGG(1)=0.0D0
DO 2000 I=2,KG
RG(I) =1.0D0/(GR(I))
RGG(I)=1.0D0/(GR(I)*GR(I))
2000 CONTINUE
C
REWIND 16
DO 500 II=1,KTYP
ETOT1 = ETOT1 + FLOAT(IATOM(II))*PAI*ACHG(II)
& * ( AC(II,1)/BC(II,1) + AC(II,2)/BC(II,2) )/UNIVOL
TOTCH=TOTCH + FLOAT(IATOM(II))*ACHG(II)
C
IF (ITER.EQ.0) THEN
IF (NLSPD(II).EQ.1) THEN
READ(16,*) MESHR,NMES,DX,RAD,VD,VDNL <--- Read local ps data!
DO 1210 N=1,MESHR
VDD(N,II)=VD(N)
VDDNL(N,II)=VDNL(N)
1210 CONTINUE
ELSE
DO 1212 N=1,MESHR
VDD(N,II) =0.0D0
VDDNL(N,II)=0.0D0
1212 CONTINUE
END IF
ELSE
DO 1211 N=1,MESHR
VD(N)=VDD(N,II)
VDNL(N)=VDDNL(N,II)
1211 CONTINUE
END IF
C
S=0.0D0
DO 1200 N=1,MESHR
S=S + OMO(N)*VDNL(N)*(RAD(N)**3)
1200 CONTINUE
PSC(1,II)=0.0D0
DSC(1,II)=0.0D0
ETOT1 =ETOT1+FLOAT(IATOM(II))*4.0D0*PAI*DX*S/UNIVOL
C CORE-PART NCP MATRIX ELEMENT; PSCO
C1=AC(II,1)*4.D0*PAI*ACHG(II)
C2=AC(II,2)*4.D0*PAI*ACHG(II)
B1=-0.25D0/BC(II,1)
B2=-0.25D0/BC(II,2)
DO 82 I=2,KG
GG=GR(I)**2
PSC(I,II)=-( C1*EXP(B1*GG) + C2*EXP(B2*GG) )/(GG*UNIVOL)
DSC(I,II)=-( C1*(B1-RGG(I))*EXP(B1*GG)
& +C2*(B2-RGG(I))*EXP(B2*GG) )*RGG(I)*UNIRV
82 CONTINUE
C
IF (NLSPD(II).EQ.1) THEN
C
DO 1201 N=1,MESHR
R=RAD(N)
FAC1=OMO(N)*VDNL(N)*RAD(N)**3*4.0D0*PAI*DX/UNIVOL
FAC2=OMO(N)*VDNL(N)*RAD(N)**4*2.0D0*PAI*DX
DO 1110 I=1,KG
XX(I)=R*GR(I)
1110 CONTINUE
CALL DSJNV(0,KG,XX,Y11)
CALL DSJNV(1,KG,XX,Y22)
DO 1120 I=2,KG
PSC(I,II)=PSC(I,II)+FAC1*Y11(I)
DSC(I,II)=DSC(I,II)-FAC2*Y22(I)*RG(I)*UNIRV
1120 CONTINUE
1201 CONTINUE
END IF
500 CONTINUE
C ETOT1 -> TOTAL ENERGY CORRECTION FROM THE PSEUDOPOTENTIAL
C ETOT2 -> TOTAL ENERGY CONTRIBUTION FROM THE EWALD SUMMATION
WRITE(6,*) 'TOTAL VAL CHG PER UNIT-CELL = ',TOTCH
WRITE(6,*) 'TOTAL ENERGY CORRECTION FROM PS= ',ETOT1
WRITE(6,*) ' UNIVOL= ',UNIVOL
DO 75 I=1,KG
ZPSCC(I) = DCMPLX(0.0D0,0.0D0)
ZDSCC(I) = DCMPLX(0.0D0,0.0D0)
75 CONTINUE
DO 73 IT=1,KTYP
DO 72 I=1,KG
ZPSCC(I) = ZPSCC(I) + PSC(I,IT )*ZFM3(I,IT)
ZDSCC(I) = ZDSCC(I) + DSC(I,IT )*ZFM3(I,IT)
72 CONTINUE
73 CONTINUE
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE KBMAT
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
C VPP-PARALLEL
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
DIMENSION SSS(KNG1,KNV3,KTYP,10),RAA(KNG1,KNV3)
!XOCL LOCAL SSS(:,/IP,:,:),RAA(:,/IP)
EQUIVALENCE (SNL,SSS),(RAK,RAA)
C K-B TYPE ADDITION
DIMENSION RKS(KIZA),RKP(KIZA),RKD(KIZA),DVS(KIZA)
& ,DVP(KIZA),DVD(KIZA)
NSEK=2
SC=DSQRT(PAI4)
PC=DSQRT(3.0D0*PAI4)
DC=DSQRT(5.0D0*PAI4)
CO1=DSQRT(3.0D0)
REWIND 15
C INTEGRATION OF K-B TYPE POTENTIAL
DO 100 ITY = 1,KTYP
IF (NLSPD(ITY).EQ.1) THEN
READ (15,*) MESHR,NMES,DX,WS(ITY),WP(ITY) <--- Read non-local ps data!
& ,RAD,DVS,DVP
WD(ITY)=0.0D0
DO 111 N=1,MESHR
WVS(N,ITY)=DVS(N)
WVP(N,ITY)=DVP(N)
WVD(N,ITY)=0.0D0
111 CONTINUE
ELSE
READ (15,*) MESHR,NMES,DX,WS(ITY),WP(ITY),WD(ITY)
& ,RAD,DVS,DVP,DVD
DO 110 N=1,MESHR
WVS(N,ITY)=DVS(N)
WVP(N,ITY)=DVP(N)
WVD(N,ITY)=DVD(N)
110 CONTINUE
END IF
100 CONTINUE
C
!XOCL SPREAD DO /IP
DO 1003 IK=1,KV3
AKX = VX(IK)
AKY = VY(IK)
AKZ = VZ(IK)
C IIBA = IBA(IK)
DO 1005 I=1,IBA2(IK)
RK=SQRT((AKX+GX(NBMAT(I,IK)))**2
& +(AKY+GY(NBMAT(I,IK)))**2 + (AKZ+GZ(NBMAT(I,IK)))**2 )
*VOCL STMT,IF(90)
IF (RK.NE.0.0D0) THEN
RAA(I,IK)=1.0D0/RK
ELSE
RAA(I,IK)=0.0D0
END IF
1005 CONTINUE
1003 CONTINUE
!XOCL END SPREAD
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 200 NNN=1,KV3
AKX = VX(NNN)
AKY = VY(NNN)
AKZ = VZ(NNN)
IIBA = IBA2(NNN)
DO 1000 J =1,IIBA
QX(J) = AKX + GX(NBMAT(J,NNN))
QY(J) = AKY + GY(NBMAT(J,NNN))
QZ(J) = AKZ + GZ(NBMAT(J,NNN))
QA1(J) = 0.5D0*(2.0D0*QZ(J)*QZ(J)-QX(J)*QX(J)-QY(J)*QY(J))
QA2(J) = 0.5D0*CO1*(QX(J)*QX(J)-QY(J)*QY(J))
QA3(J) = CO1*QX(J)*QY(J)
QA4(J) = CO1*QX(J)*QZ(J)
QA5(J) = CO1*QY(J)*QZ(J)
AK2(J) = SQRT( QX(J)**2 + QY(J)**2 + QZ(J)**2 )
1000 CONTINUE
DO 1100 ITY=1,KTYP
DO 1111 KN=1,KNG11
SSS(KN,NNN,ITY,1)=0.0D0
SSS(KN,NNN,ITY,2)=0.0D0
SSS(KN,NNN,ITY,3)=0.0D0
SSS(KN,NNN,ITY,4)=0.0D0
1111 CONTINUE
C
DO 7320 N=1,MESHR
R=RAD(N)
FACS= OMO(N)*WVS(N,ITY)*R**NSEK*SC*DX
FACP= OMO(N)*WVP(N,ITY)*R**NSEK*PC*DX
DO 7777 KN=1,IIBA
X(KN)=AK2(KN)*R
7777 CONTINUE
CALL DSJNV(0,IIBA,X,Y1)
CALL DSJNV(1,IIBA,X,Y2)
DO 7330 KN=1,IIBA
C RKS(N)=$DMSJM(0,AK2(KN)*RAD(N))
SSS(KN,NNN,ITY,1) =SSS(KN,NNN,ITY,1) +FACS*Y1(KN)
SSS(KN,NNN,ITY,2) =SSS(KN,NNN,ITY,2)
& +FACP*Y2(KN)*QX(KN)*RAA(KN,NNN)
SSS(KN,NNN,ITY,3) =SSS(KN,NNN,ITY,3)
& +FACP*Y2(KN)*QY(KN)*RAA(KN,NNN)
SSS(KN,NNN,ITY,4) =SSS(KN,NNN,ITY,4)
& +FACP*Y2(KN)*QZ(KN)*RAA(KN,NNN)
7330 CONTINUE
7320 CONTINUE
C
IF (NLSPD(ITY).EQ.2) THEN
DO 1112 KN=1,KNG11
SSS(KN,NNN,ITY,5)=0.0D0
SSS(KN,NNN,ITY,6)=0.0D0
SSS(KN,NNN,ITY,7)=0.0D0
SSS(KN,NNN,ITY,8)=0.0D0
SSS(KN,NNN,ITY,9)=0.0D0
1112 CONTINUE
C
DO 7321 N=1,MESHR
R=RAD(N)
FACD= OMO(N)*WVD(N,ITY)*R**NSEK*DC*DX
DO 7778 KN=1,IIBA
X(KN)=AK2(KN)*R
7778 CONTINUE
CALL DSJNV(2,IIBA,X,Y3)
DO 7331 KN=1,IIBA
C RKS(N)=$DMSJM(0,AK2(KN)*RAD(N))
SSS(KN,NNN,ITY,5) =SSS(KN,NNN,ITY,5)
& +FACD*Y3(KN)*QA1(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SSS(KN,NNN,ITY,6) =SSS(KN,NNN,ITY,6)
& +FACD*Y3(KN)*QA2(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SSS(KN,NNN,ITY,7) =SSS(KN,NNN,ITY,7)
& +FACD*Y3(KN)*QA3(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SSS(KN,NNN,ITY,8) =SSS(KN,NNN,ITY,8)
& +FACD*Y3(KN)*QA4(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SSS(KN,NNN,ITY,9) =SSS(KN,NNN,ITY,9)
& +FACD*Y3(KN)*QA5(KN)*RAA(KN,NNN)*RAA(KN,NNN)
7331 CONTINUE
7321 CONTINUE
END IF
1100 CONTINUE
200 CONTINUE
!XOCL END SPREAD
C!XOCL END PARALLEL
RETURN
END
Description
Subroutine PSEUDO
A Subroutine "PSEUDO" is to transform from real space local pseudopotential data to reciprocal space data by using the Fourier transformation.
In this transformation, an integration of radial direction is performed as follows,
V_loc(G) = 4*pai*Int{j_0(Gr)*V_loc(r)*r^3}dr/UNIVOL.
- int{} :Integration in a real space(radial part)
- V_loc(G):Local pseudopotential in a reciprocal space
- V_loc(r):Local pseudopotential in a real space(radial part)
- j_0(Gr) :Bessel function(l=0)
- UNIVOL :Volume of a unit cell
r^2 transforms to r^3 in a logarithmic radial space. V_loc(r) consists of two parts as follows,
V_loc(r) = {V_d(r)-Vcore(r)} + Vcore(r)
Vcore(r) is a core potential defined in BHS paper, V_d(r) is a d pseudopotential. It is possible to consider two cases. One case is s,p and d non-local, the other is s,p(only) non-local. In the s,p and d non-local case, {V_d(r)-Vcore(r)} is treated as a non-local part. Therefore, only Vcore(r) is integrated in the radial space, where it is possible to integrate analytically as for Vcore(r).
In the later case, {V_d(r)-Vcore(r)} is integrated numerically and the integration of Vcore(r) is also calculated analytically.
Subroutine KBMAT
A Subroutine "KBMAT" is to transform from real space non-local pseudopotential data to reciprocal space data by using the Fourier transformation.
Non-local pseudopotentials are treated as separable forms(see from Kleinman and Bylander[Phys. Rev. Lett., 48, 1425(1982)]).
Also, in this routine, two cases are considered that one is s,p and d non-local, the other is s and p(only) non-local.
The way of transformation(and integration) is almost same of the local pseudopotential case.
A formulation of integration in order to obtain V_nl(G)[s,p,d and f parts] for a radial direction is shown in [V_nl(G) <-- FFT -- V_nl(r)](with png image[about 30kb]) .
[To Guide][Top]