C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C @@ BAND STRUCTURE AND TOTAL ENERGY CALCULATION @@
C @@ (ORIGINAL FROM DR. H. ISHIDA) @@
C @@ 1989 5/4 (START) BY KOBAYASHI KAZUAKI @@
C @@ MODIFIED BY Y. MORIKAWA(SYMMETRY)@@
C @@ 1992 7/14 MODIFIED TO STRESS CALCULATION @@
C @@ 1993 3/16 MODIFIED TO CONSTANT PRESSURE @@
C @@ AUTOMATIC OPT (ANDERSEN ARG) @@
C @@ 1994 10/25 MODIFIED TO INTRODUCTION OF PCC @@
C @@ 1994 11/22 MODIFIED TO CONSTANT PRESSURE @@
C @@ AUTOMATIC OPT (PARRINELLO-RAHMAN)@@
C @@ 1994 12/16 MODIFIED TO INPUT/OUTPUT OF RECAL@@
C @@ 1995 1/2 MODIFIED TO VPP-PARALEEL VERSION @@
C @@ 1995 2/13 VPP-PARALLEL VER 1.01@@
C @@ 1995 6/2 VPP-PARALLEL VER 1.1 @@
C @@ 1995 7/17 MODIFIED TO ADD NON-LOCAL D 1.2 @@
C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
C VPP-PARALLEL
C*****FOR PCC
DIMENSION RHPC(KIZA,KTYP)
INTEGER KTPCC(KTYP),KMDTP(KATM)
C*****FOR KPMIRB
DIMENSION QWGT(KNV3)
C*****FOR SYMM
DIMENSION FORC(KATM,3),FORCW(KATM,3)
DIMENSION EPP(3,3)
DIMENSION RX(125),RY(125),RZ(125),RR(125)
INTEGER IRX(125),IRY(125),IRZ(125)
!XOCL PARALLEL REGION
KOPR=KO
KIMG= 1
C-----DUMMY KIMG--------------------------------------ATTN]]] 8*KG1
ITIME1=0
WRITE (6,*) 'START OF PROGRAM (REVPE_D.F)'
IREC8=0
IREC9=0
KREC8=16*KNG1
KREC9=16*KFFT
KR8=KREC8/32760
KR9=KREC9/32760
IF (KR8.NE.0) THEN
IREC8=KR8+1
IF (MOD(KREC8,IREC8).EQ.0) THEN
KREC8=KREC8/IREC8
ELSE
KREC8=KREC8/IREC8+1
END IF
END IF
IF (KR9.NE.0) THEN
IREC9=KR9+1
IF (MOD(KREC9,IREC9).EQ.0) THEN
KREC9=KREC9/IREC9
ELSE
KREC9=KREC9/IREC9+1
END IF
END IF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1400 CONTINUE
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
INCLUDE 'REINIT'
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C FOR MOVING ATOMS KMDTP=1 : FIXED ATOMS KMDTP=0
WRITE(6,*) 'MD CONDITION FOR EACH ATOMS'
DO 101 I=1,KATM
KMDTP(I) = 1
101 CONTINUE
C DO 103 I= 9,12
C KMDTP(I) = 0
C 103 CONTINUE
WRITE(6,102) (I,KMDTP(I),I=1,KATM)
102 FORMAT(20I3)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL SIMP
CALL TIME(ITIME1,'SIMP ',2)
1300 CONTINUE
C^^^^ READ UNIT = 1 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL INPUT1(ALFA)
CALL TIME(ITIME1,'INPUT1',2)
IF (ICONT.EQ.1.AND.ICAF.EQ.1) IMD=0
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL INFOUT(ICAR,ICAF,IPRE,IPCC,IMD,IOVE,IRES,IORW,ISSS
& ,ISTMD,ICOST,IOUT,IMDI,NBAND1,NBAND2
& ,IREC8,IREC9,KREC8,KREC9,KBZTYP
& ,WIDTH,CMIX,DTIMEL,DTIMIO,DTIMUC,KTPCC)
CALL TIME(ITIME1,'INFOUT',2)
C^^^^ UNIT CELL CHANGE ^^^^^^^^^^^^^^^^^^^^^^^^^^^
ITER=0
1200 CONTINUE
IF (ICAR.NE.2) THEN
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL LATTIC
< (RX,RY,RZ,RR,IRX,IRY,IRZ)
CALL TIME(ITIME1,'LATTIC',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL GSTEP1(IPRE)
CALL TIME(ITIME1,'GSTEP1',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL HPSORT(KG,IG1,IG2,IG3,GX,GY,GZ,GR)
CALL TIME(ITIME1,'BUBBLG',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IF(KBZTYP.GE.2) THEN
CALL TIME(ITIME1,' ',1)
CALL SYMM(KOPR,NOPR,KBZTYP,125,6,RX,RY,RZ)
CALL TIME(ITIME1,'SYMM ',2)
END IF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL GSTEPF
CALL TIME(ITIME1,'GSTEPF',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL KSTEP
> (KNV3,KBZTYP,ALTV,RLTV,
> KNVX,KNVY,KNVZ,KNVX2,KNVY2,KNVZ2,6,
< KV3,VX,VY,VZ,QWGT)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL BASNUM(QWGT)
CALL TIME(ITIME1,'BASNUM',2)
C^^^^ PRE-CALCULATION IPRE=1 ^^^^^^^^^^^^^^^^^^^^^
IF (IPRE.EQ.1) THEN
WRITE (6,*) 'THIS IS PRE-CALCULATION] PROGRAM STOPS]'
STOP
END IF
ELSE IF (ICAR.EQ.2 .AND. ISTR.EQ.1) THEN
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL LATSCA
< (RX,RY,RZ,RR,IRX,IRY,IRZ)
WRITE (6,*) ALTV(1,1),ALTV(1,2),ALTV(1,3)
WRITE (6,*) ALTV(2,1),ALTV(2,2),ALTV(2,3)
WRITE (6,*) ALTV(3,1),ALTV(3,2),ALTV(3,3)
C WRITE (6,*) EPP(1,1),EPP(1,2),EPP(1,3)
C WRITE (6,*) EPP(2,1),EPP(2,2),EPP(2,3)
C WRITE (6,*) EPP(3,1),EPP(3,2),EPP(3,3)
CALL TIME(ITIME1,'LATTIC',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C WRITE (6,*) 'CALL START GSTSCA'
CALL GSTSCA(IPRE,EPP)
C WRITE (6,*) 'CALL END GSTSCA'
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C CALL HPSORT(KG,IG1,IG2,IG3,GX,GY,GZ,GR)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IF(KBZTYP.GE.2) THEN
CALL SYMSCA(KOPR,NOPR,KBZTYP,NEIBRD,NFOUT
& ,RX,RY,RZ,EPP)
END IF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL GSFSCA
CALL TIME(ITIME1,'GSTEPF',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL KSTEP
> (KNV3,KBZTYP,ALTV,RLTV,
> KNVX,KNVY,KNVZ,KNVX2,KNVY2,KNVZ2,6,
< KV3,VX,VY,VZ,QWGT)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
C CALL BASNUM(QWGT)
CALL TIME(ITIME1,'BASNUM',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
END IF
1100 CONTINUE
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL FORM
CALL TIME(ITIME1,'FORM ',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IF (IPCC.EQ.1) THEN
C CALL TIME(ITIME1,' ',1)
CALL PCC(KTPCC,RHPC,
< SCHGPC,EPC,PCM)
C CALL TIME(ITIME1,'PCC ',2)
ELSE IF (ITER.EQ.0) THEN
DO 1010 I=1,KNG
ZRHPC(I)=DCMPLX(0.D0,0.D0)
ZRRPC(I)=DCMPLX(0.D0,0.D0)
EPC = 0.D0
1010 CONTINUE
END IF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
IF (ITER.EQ.0 .OR. ISTR.EQ.1) THEN
WRITE (6,*) 'CALL PSEUDO'
CALL PSEUDO(TOTCH,ETOT1)
ELSE
WRITE (6,*) 'CALL PSEUMD'
CALL PSELMD
END IF
CALL TIME(ITIME1,'PSEUDO',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
IF (ICAR.EQ.0 .OR. ICAR.EQ.2) CALL EWVEC(ETOT2)
IF (ICAR.EQ.1) CALL EWVMD(ETOT2)
CALL TIME(ITIME1,'EWALDV',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IF (ITER.EQ.0.AND.ICONT.EQ.0) THEN
CALL TIME(ITIME1,' ',1)
CALL INTCHG(ALFA)
CALL TIME(ITIME1,'INTCHG',2)
END IF
C^^^^ LOOP ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1000 CONTINUE
ITER=ITER+1
IF (ICONT.EQ.1 .AND. IORW.EQ.0) THEN
IORW=1
CALL EVIN(IREC8,IREC9,EPC)
WRITE (6,*) 'END OF EVIN'
END IF
IF (ICAF.NE.0.AND.ITER.GT.IMD) ICAR=ICAF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL XCFFT(1,SCHGPC)
C^^^^^^^^^^^^^^^^^^^^^^^^^DEBUG^^^^^^^^^^^^^^^^^^^
IF (ITER.EQ.1) THEN
CALL TIME(ITIME1,' ',1)
CALL KBMAT
CALL DIAGON(IREC8)
CALL KBINT
CALL TIME(ITIME1,'KBMSD ',2)
ELSE
CALL TIME(ITIME1,' ',1)
CALL MSD(IREC8,IREC9,IMD,ICAR,ISTR)
CALL TIME(ITIME1,'MSD ',2)
END IF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IF (ICAR.EQ.2) THEN
IF (MOD(ITER,IOVE).EQ.IRES) THEN
ISTR=1
ELSE
ISTR=0
END IF
END IF
C
CALL TIME(ITIME1,' ',1)
CALL FERMI(TOTCH,QWGT,WIDTH)
CALL TIME(ITIME1,'FERMI ',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
IF (MOD(ITER,IMDI).EQ.IRES) THEN
IF (ITER.EQ.1 .OR. ITER.GT.IMD) THEN
WRITE (6,*) 'CALL FORCE'
CALL FORCE(IREC8)
ELSE
CALL FORZFB
END IF
ELSE
WRITE (6,*) 'CALL FORZFB'
CALL FORZFB
END IF
CALL TIME(ITIME1,'FORCE ',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL CHAVER(IREC8,IREC9,KOPR,NOPR,KBZTYP)
CALL TIME(ITIME1,'CHAVER',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
IF (IPCC.EQ.1) CALL XCFFT(3,SCHGPC)
CALL TIME(ITIME1,'XCFFT ',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IF (ITER.EQ.IMD+1) ISTR=1
C STRESS
C IF (ITER.EQ.100 .OR. ITER.EQ.300) THEN
IF (ISTR.EQ.1) THEN
CALL TIME(ITIME1,' ',1)
CALL STRESS(IPCC,SCHGPC
& ,ETOT1,TOTCH,EPC,PCM,KOPR,KBZTYP)
CALL TIME(ITIME1,'STRESS',2)
END IF
CALL TIME(ITIME1,'FORLOC',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IF (ITER.GT.IMD .AND. ICAR.GE.1) THEN
DTIO=0.0D0
IF (MOD(ITER,IMDI).EQ.IRES .OR. ISTR.EQ.1) THEN
WRITE (6,*) 'CALL FORLOC AND MD'
CALL TIME(ITIME1,' ',1)
CALL FORLOC(KTPCC)
IF(KBZTYP.GE.2) THEN
DO 2500 IA=1,KATM
FORC(IA,1)=REAL(ZFORC2(IA,1))
FORC(IA,2)=REAL(ZFORC2(IA,2))
FORC(IA,3)=REAL(ZFORC2(IA,3))
2500 CONTINUE
CALL FORCES(KO,NOPR,FORC,FORCW)
DO 2501 IA=1,KATM
ZFORC2(IA,1)=CMPLX(FORC(IA,1),0.0D0)
ZFORC2(IA,2)=CMPLX(FORC(IA,2),0.0D0)
ZFORC2(IA,3)=CMPLX(FORC(IA,3),0.0D0)
2501 CONTINUE
END IF
IF (MOD(ITER,10).EQ.1) THEN
WRITE (17,*) 'ITER = ',ITER
WRITE (17,3500)(I,ZFORC2(I,1),ZFORC2(I,2),ZFORC2(I,3),I=1,KATM)
END IF
3500 FORMAT(I3,6F11.7)
C
DTIO=DTIMIO
CALL MD(IMD,KMDTP,MDFLAG,ICAR,ITERMD,ITERMX,ISTR,ISSS
& ,DTUC,EPP)
END IF
END IF
IF (ITER.EQ.ISTMD) DTIM=DTIMEL
IF (ITER.EQ.ISTMD) ICAR=0
IF (ITER.EQ.ISTMD) ICAF=0
C^^^^^^^^^^^^^^^^^^^^^^^^^DEBUG^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL ENERGY(TOTCH,ETOT1,ETOT2,ETONEW,SCHGPC,EPC)
CALL TIME(ITIME1,'ENERGY',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL TIME(ITIME1,' ',1)
CALL CONV2(ETONEW,ETOOLD,IMD,ITERMD,CMIX)
CALL TIME(ITIME1,'CONV2 ',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IF (MOD(ITER,ICOST).EQ.IRES) THEN
WRITE (6,*) 'READ 2 --- REINIT OR NOT'
REWIND 2
READ (2,*) IREINT
WRITE (6,*) 'READ 2 --- REINIT OR NOT',IREINT
IF (IREINT.GE.1) THEN
WRITE (6,*) 'READ 2 --- REINIT'
READ (2,*) DTIMEL,DTIMIO,DTIMUC
READ (2,*) CMIX,WIDTH,TEMPSB
C SET A EXTERNAL PRESSURE
READ (2,*) PEX(1)
READ (2,*) PEX(2)
READ (2,*) PEX(3)
READ (2,*) PEX(4)
READ (2,*) PEX(5)
READ (2,*) PEX(6)
READ (2,*) IPRE,ITRA1,ITRA2,IPCC,KBZTYP
READ (2,*) IMDC,IMDIC,IOVEC,IRESC,ISTOP,ISTEL,ICOST,IOUT
DO 1238 ITYP=1,KTYP
READ (2,*) KTPCC(ITYP)
1238 CONTINUE
ISTMD = ISTOP - ISTEL
DTIM = DTIMEL
DTIO = DTIMIO
DTUC = DTIMUC
C.................................................
WRITE (6,*) DTIMEL,DTIMIO,DTIMUC
WRITE (6,*) CMIX,WIDTH,TEMPSB
WRITE (6,*) PEX(1)
WRITE (6,*) PEX(2)
WRITE (6,*) PEX(3)
WRITE (6,*) PEX(4)
WRITE (6,*) PEX(5)
WRITE (6,*) PEX(6)
WRITE (6,*) IPRE,ICAF,ICAR,ITRA1,ITRA2,IPCC,KBZTYP
WRITE (6,*) IMDC,IMDIC,IOVEC,IRESC,ISTOP,ISTEL,ICOST,IOUT
DO 1336 ITYP=1,KTYP
WRITE (6,*) KTPCC(ITYP)
1336 CONTINUE
END IF
C
WRITE (6,*) 'EVOUT ITER = ',ITER
IF (IOUT.EQ.1) THEN
CALL EVOUT(IREC8,IREC9,ALFA,EPC)
END IF
IF (IREINT.EQ.2) GO TO 1300
IF (IREINT.EQ.3) GO TO 1400
END IF
C
IF (ITER.GE. ISTOP ) THEN
CALL EVOU2(ALFA)
GOTO 4100
END IF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IF (ITER.GT.IMD) THEN
IF (ISTR.EQ.1) GO TO 1200
IF (MOD(ITER,IMDI).EQ.IRES) GO TO 1100
END IF
GO TO 1000
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
4100 CONTINUE
!XOCL END PARALLEL
C WRITE(34,4200) (ZCHG(I),I=1,KG)
C*****DO 4110 IBAN = 38,44,2
C NNN=1
C WRITE(12,4210) IBAN,NNN
C* IREC = KV3*(IBAN-1)+IWRT(NNN)
C IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1
C READ(80,REC=IREC) ZV1
C* WRITE(12,4200) ( ZV1(I),I=1,IBA(NNN) )
C WRITE(12,4200) ( ZAJ(I,IBAN,NNN),I=1,IBA(NNN) )
C NNN=5
C WRITE(12,4210) IBAN,NNN
C* IREC = KV3*(IBAN-1)+IWRT(NNN)
C IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1
C READ(80,REC=IREC) ZV1
C* WRITE(12,4200) ( ZV1(I),I=1,IBA(NNN) )
C WRITE(12,4200) ( ZAJ(I,IBAN,NNN),I=1,IBA(NNN) )
C NNN=17
C WRITE(12,4210) IBAN,NNN
C* IREC = KV3*(IBAN-1)+IWRT(NNN)
C IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1
C READ(80,REC=IREC) ZV1
C* WRITE(12,4200) ( ZV1(I),I=1,IBA(NNN) )
C WRITE(12,4200) ( ZAJ(I,IBAN,NNN),I=1,IBA(NNN) )
C NNN=21
C WRITE(12,4210) IBAN,NNN
C* IREC = KV3*(IBAN-1)+IWRT(NNN)
C IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1
C READ(80,REC=IREC) ZV1
C* WRITE(12,4200) ( ZV1(I),I=1,IBA(NNN) )
C WRITE(12,4200) ( ZAJ(I,IBAN,NNN),I=1,IBA(NNN) )
C4110 CONTINUE
C WRITE(14,4220) KG,KX1,KY1,KZ1
C WRITE(14,4230) (IG1(I),I=1,KG)
C WRITE(14,4230) (IG2(I),I=1,KG)
C WRITE(14,4230) (IG3(I),I=1,KG)
C4200 FORMAT(4D18.10)
C4210 FORMAT(' ',I4,3X,I4)
C4220 FORMAT(4I6)
C4230 FORMAT(20I4)
C CLOSE(UNIT=90,STATUS='KEEP')
STOP
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE INFOUT(ICAR,ICAF,IPRE,IPCC,IMD,IOVE,IRES,IORW,ISSS
& ,ISTMD,ICOST,IOUT,IMDI,NBAND1,NBAND2
& ,IREC8,IREC9,KREC8,KREC9,KBZTYP
& ,WIDTH,CMIX,DTIMEL,DTIMIO,DTIMUC,KTPCC)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
INTEGER KTPCC(KTYP),KMDTP(KATM)
C
WRITE (6,*) 'IPRE = ',IPRE
WRITE (6,*) 'ICAR = ',ICAR
WRITE (6,*) 'ICAF = ',ICAF
WRITE (6,*) 'IPRI = ',IPRI
WRITE (6,*) 'ICONT = ',ICONT
WRITE (6,*) 'IPCC = ',IPCC
WRITE (6,*) 'KBZTYP = ',KBZTYP
WRITE (6,*) 'IMD = ',IMD
WRITE (6,*) 'IOVE = ',IOVE
WRITE (6,*) 'IMDI = ',IMDI
WRITE (6,*) 'IRES = ',IRES
WRITE (6,*) 'IORW = ',IORW
WRITE (6,*) 'ISSS = ',ISSS
WRITE (6,*) 'ISTOP = ',ISTOP
WRITE (6,*) 'ISTMD = ',ISTMD
WRITE (6,*) 'ICOST = ',ICOST
WRITE (6,*) 'IOUT = ',IOUT
WRITE (6,*) 'IREC8 = ',IREC8
WRITE (6,*) 'IREC9 = ',IREC9
WRITE (6,*) 'KREC8 = ',KREC8
WRITE (6,*) 'KREC9 = ',KREC9
WRITE (6,*) 'UP-BAND = ',NBAND1
WRITE (6,*) 'DOWN-BAND = ',NBAND2
C OUTPUT INCLUDE FILE INFORMATION
WRITE (6,*) 'OUTPUT INCLUDE FILE INFORMATION'
WRITE (6,*) 'KNX,Y,Z =',KNX,KNY,KNZ
WRITE (6,*) 'KNG,1,11 =',KNG,KNG1,KNG11
WRITE (6,*) 'KNBMX =',KNBMX
WRITE (6,*) 'KNVX,Y,Z =',KNVX,KNVY,KNVZ
WRITE (6,*) 'KVX2,Y2,Z2=',KNVX2,KNVY2,KNVZ2
WRITE (6,*) 'KEG,KBD1,2=',KEG,KBD1,KBD2
WRITE (6,*) 'KTYP,KATM =',KTYP,KATM
WRITE (6,*) 'EWALD =',KLX,KLY,KLZ,KMX,KMY,KMZ
WRITE (6,*) 'DTIMEL =',DTIMEL
WRITE (6,*) 'DTIMIO =',DTIMIO
WRITE (6,*) 'DTIMUC =',DTIMUC
WRITE (6,*) 'WIDTH =',WIDTH
WRITE (6,*) 'CMIX =',CMIX
WRITE (6,*) 'TEMPSB =',TEMPSB
DO 1000 ITYP=1,KTYP
WRITE (6,*) 'KTPCC =',KTPCC(ITYP)
1000 CONTINUE
DO 1010 I=1,6
WRITE (6,*) 'PEX = ',I,PEX(I)
1010 CONTINUE
DO 2000 I=1,KATM
WRITE (6,*) 'AMION = ',AMION(I)
2000 CONTINUE
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE SIMP
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
PARAMETER(KIZA=421)
COMMON/OM/ OMO(KIZA)
OMO(1) = 1.0D0
DO 1101 N = 2,KIZA-1,2
OMO(N) = 4.0D0
OMO(N+1) = 2.0D0
1101 CONTINUE
OMO(KIZA) = 1.0D0
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE INPUT1(ALFA)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
C'''''READING THE DATA (UNIT=1)''''''''''''''''''''''''''''''''''
REWIND 1
READ(1,*) ITEMAX,PPMIX,CONV,GMAX,ICONT
WRITE(6,99) ITEMAX,PPMIX,CONV,GMAX,ICONT
99 FORMAT(' ',I6,F8.4,D12.4,F8.4,I4)
DO 100 I=1,3
READ(1,*) ALTV(1,I),ALTV(2,I),ALTV(3,I)
WRITE(6,200) ALTV(1,I),ALTV(2,I),ALTV(3,I)
100 CONTINUE
200 FORMAT(3(F20.10))
READ(1,*) KCOTYP
WRITE(6,*)KCOTYP,' COORDINATES 0:NORMALIZED 1:CARTESIAN '
IF (KCOTYP.EQ.0) THEN
DO 300 IA=1,KATM
READ(1,*) PPOS(IA),QPOS(IA),RPOS(IA)
WRITE(6,200) PPOS(IA),QPOS(IA),RPOS(IA)
300 CONTINUE
DO 400 IA=1,KATM
CATX(IA)=PPOS(IA)*ALTV(1,1)+QPOS(IA)*ALTV(1,2)
& +RPOS(IA)*ALTV(1,3)
CATY(IA)=PPOS(IA)*ALTV(2,1)+QPOS(IA)*ALTV(2,2)
& +RPOS(IA)*ALTV(2,3)
CATZ(IA)=PPOS(IA)*ALTV(3,1)+QPOS(IA)*ALTV(3,2)
& +RPOS(IA)*ALTV(3,3)
400 CONTINUE
ELSE IF(KCOTYP.EQ.1) THEN
DO 500 IA=1,KATM
READ(1,*) CATX(IA),CATY(IA),CATZ(IA)
WRITE(6,200) CATX(IA),CATY(IA),CATZ(IA)
500 CONTINUE
ELSE
WRITE(6,*) ' KCOTYP = ',KCOTYP
STOP
END IF
NNATM=0
DO 600 IT=1,KTYP
READ(1,*) IATOM(IT),NLSPD(IT)
WRITE(6,*) 'IATOM,NLSPD=',IATOM(IT),NLSPD(IT)
READ(1,*) AICH,ALFA
WRITE(6,700) IT,AICH,ALFA
READ(1,*) ACHG(IT),AC(IT,1),AC(IT,2),BC(IT,1),BC(IT,2)
WRITE(6,701) ACHG(IT),AC(IT,1),AC(IT,2),BC(IT,1),BC(IT,2)
DO 800 IA=NNATM+1,NNATM+IATOM(IT)
AICHG(IA) = AICH
KFTYPE(IA) = IT
800 CONTINUE
NNATM=NNATM+IATOM(IT)
600 CONTINUE
700 FORMAT(' ','TYPE=',I3,' CHARGE=',F8.4,'GAUSS COEFF.=',F8.4)
701 FORMAT(' ','AC,BC=',5F12.6)
IF (NNATM.NE.KATM) THEN
WRITE(6,*) ' NNATM= ',NNATM
STOP
END IF
KV3 = KNV3
C'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
PINIT=PPMIX
DO 1 I=1,ITEMAX
PMIX(I)=PPMIX
1 CONTINUE
DO 2 K=1,KV3
DO 3 I=1,KEG
OCCUP(I,K)=1.0D0
3 CONTINUE
2 CONTINUE
DO 4 I=1,KATM
VEL(I,1) =0.0D0
VEL(I,2) =0.0D0
VEL(I,3) =0.0D0
4 CONTINUE
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE LATTIC
< (RX,RY,RZ,RR,IRX,IRY,IRZ)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
REAL RX(125),RY(125),RZ(125),RR(125)
INTEGER IRX(125),IRY(125),IRZ(125)
NEIBR = 2
NEIBRD = 125
C
ALEN1=SQRT(ALTV(1,1)**2+ALTV(2,1)**2+ALTV(3,1)**2)
ALEN2=SQRT(ALTV(1,2)**2+ALTV(2,2)**2+ALTV(3,2)**2)
ALEN3=SQRT(ALTV(1,3)**2+ALTV(2,3)**2+ALTV(3,3)**2)
ALMAX=MAX(ALEN1,ALEN2,ALEN3)
ALF=1.0D0/ALMAX
WRITE(6,*) ' <<>>'
FFF= ALTV(1,1)*(ALTV(2,2)*ALTV(3,3)-ALTV(3,2)*ALTV(2,3))
& + ALTV(2,1)*(ALTV(3,2)*ALTV(1,3)-ALTV(1,2)*ALTV(3,3))
& + ALTV(3,1)*(ALTV(1,2)*ALTV(2,3)-ALTV(2,2)*ALTV(1,3))
UNIVOL=ABS( FFF )
WRITE(6,708) UNIVOL
708 FORMAT(1H ,'VOLUME OF A UNIT CELL=',F15.6)
FFF=2.D0*PAI/FFF
RLTV(1,1)=(ALTV(2,2)*ALTV(3,3)-ALTV(3,2)*ALTV(2,3))*FFF
RLTV(2,1)=(ALTV(3,2)*ALTV(1,3)-ALTV(1,2)*ALTV(3,3))*FFF
RLTV(3,1)=(ALTV(1,2)*ALTV(2,3)-ALTV(2,2)*ALTV(1,3))*FFF
RLTV(1,2)=(ALTV(2,3)*ALTV(3,1)-ALTV(3,3)*ALTV(2,1))*FFF
RLTV(2,2)=(ALTV(3,3)*ALTV(1,1)-ALTV(1,3)*ALTV(3,1))*FFF
RLTV(3,2)=(ALTV(1,3)*ALTV(2,1)-ALTV(2,3)*ALTV(1,1))*FFF
RLTV(1,3)=(ALTV(2,1)*ALTV(3,2)-ALTV(3,1)*ALTV(2,2))*FFF
RLTV(2,3)=(ALTV(3,1)*ALTV(1,2)-ALTV(1,1)*ALTV(3,2))*FFF
RLTV(3,3)=(ALTV(1,1)*ALTV(2,2)-ALTV(2,1)*ALTV(1,2))*FFF
WRITE(6,707) ((RLTV(I,J),I=1,3),J=1,3)
707 FORMAT((1H ,3(F10.6,5X)))
DO 800 I=1,3
DO 810 J=1,3
ALINV(I,J)=RLTV(J,I)/(2.D0*PAI)
810 CONTINUE
800 CONTINUE
DO 703 I=1,3
DO 703 J=1,3
FFF = ALTV(I,1)*ALINV(1,J) + ALTV(I,2)*ALINV(2,J)
& + ALTV(I,3)*ALINV(3,J)
WRITE(6,702) I,J,FFF
702 FORMAT(1H ,'ALTV(',I3,')*ALINV(',I3,')=',F20.15)
703 CONTINUE
C--*--GENERATE TRANSLATIONAL VECTOR
MM=0
DO 500 I=-NEIBR,NEIBR
DO 510 J=-NEIBR,NEIBR
DO 520 K=-NEIBR,NEIBR
MM = MM + 1
IRX(MM) = I
IRY(MM) = I
IRZ(MM) = I
F1 = DFLOAT(I)
F2 = DFLOAT(J)
F3 = DFLOAT(K)
RX(MM) = ALTV(1,1)*F1 + ALTV(1,2)*F2 + ALTV(1,3)*F3
RY(MM) = ALTV(2,1)*F1 + ALTV(2,2)*F2 + ALTV(2,3)*F3
RZ(MM) = ALTV(3,1)*F1 + ALTV(3,2)*F2 + ALTV(3,3)*F3
RR(MM) = SQRT( RX(MM)**2 + RY(MM)**2 + RZ(MM)**2 )
520 CONTINUE
510 CONTINUE
500 CONTINUE
IF(MM.NE.NEIBRD) THEN
WRITE(6,*) ' MM.NE.NEIBRD MM,NEIBRD =',MM,NEIBRD
END IF
CALL HPSORT(NEIBRD,IRX,IRY,IRZ,RX,RY,RZ,RR)
IF(IPRI.GE.2) THEN
WRITE(6,*) ' TRANSLATIONAL VECTOR'
WRITE(6,600) (I,IRX(I),IRY(I),IRZ(I),
& RX(I) ,RY(I) ,RZ(I) ,RR(I),I=1,NEIBRD)
END IF
600 FORMAT((' ',4I5,4F12.6))
WRITE (6,*) 'UNIVOL(IN LATTIC) = ',UNIVOL
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE LATSCA
< (RX,RY,RZ,RR,IRX,IRY,IRZ)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
REAL RX(125),RY(125),RZ(125),RR(125)
INTEGER IRX(125),IRY(125),IRZ(125)
NEIBR = 2
NEIBRD = 125
C
ALEN1=SQRT(ALTV(1,1)**2+ALTV(2,1)**2+ALTV(3,1)**2)
ALEN2=SQRT(ALTV(1,2)**2+ALTV(2,2)**2+ALTV(3,2)**2)
ALEN3=SQRT(ALTV(1,3)**2+ALTV(2,3)**2+ALTV(3,3)**2)
ALMAX=MAX(ALEN1,ALEN2,ALEN3)
ALF=1.0D0/ALMAX
C WRITE(6,*) ' <<>>'
FFF= ALTV(1,1)*(ALTV(2,2)*ALTV(3,3)-ALTV(3,2)*ALTV(2,3))
& + ALTV(2,1)*(ALTV(3,2)*ALTV(1,3)-ALTV(1,2)*ALTV(3,3))
& + ALTV(3,1)*(ALTV(1,2)*ALTV(2,3)-ALTV(2,2)*ALTV(1,3))
UNIVOL=ABS( FFF )
WRITE(6,708) UNIVOL
708 FORMAT(1H ,'VOLUME AT LATSCA CELL=',F15.6)
FFF=2.D0*PAI/FFF
RLTV(1,1)=(ALTV(2,2)*ALTV(3,3)-ALTV(3,2)*ALTV(2,3))*FFF
RLTV(2,1)=(ALTV(3,2)*ALTV(1,3)-ALTV(1,2)*ALTV(3,3))*FFF
RLTV(3,1)=(ALTV(1,2)*ALTV(2,3)-ALTV(2,2)*ALTV(1,3))*FFF
RLTV(1,2)=(ALTV(2,3)*ALTV(3,1)-ALTV(3,3)*ALTV(2,1))*FFF
RLTV(2,2)=(ALTV(3,3)*ALTV(1,1)-ALTV(1,3)*ALTV(3,1))*FFF
RLTV(3,2)=(ALTV(1,3)*ALTV(2,1)-ALTV(2,3)*ALTV(1,1))*FFF
RLTV(1,3)=(ALTV(2,1)*ALTV(3,2)-ALTV(3,1)*ALTV(2,2))*FFF
RLTV(2,3)=(ALTV(3,1)*ALTV(1,2)-ALTV(1,1)*ALTV(3,2))*FFF
RLTV(3,3)=(ALTV(1,1)*ALTV(2,2)-ALTV(2,1)*ALTV(1,2))*FFF
WRITE(6,707) ((RLTV(I,J),I=1,3),J=1,3)
707 FORMAT((1H ,3(F10.6,5X)))
DO 800 I=1,3
DO 810 J=1,3
ALINV(I,J)=RLTV(J,I)/(2.D0*PAI)
810 CONTINUE
800 CONTINUE
DO 703 I=1,3
DO 703 J=1,3
FFF = ALTV(I,1)*ALINV(1,J) + ALTV(I,2)*ALINV(2,J)
& + ALTV(I,3)*ALINV(3,J)
C WRITE(6,702) I,J,FFF
C 702 FORMAT(1H ,'ALTV(',I3,')*ALINV(',I3,')=',F20.15)
703 CONTINUE
C--*--GENERATE TRANSLATIONAL VECTOR
MM=0
DO 500 I=-NEIBR,NEIBR
DO 510 J=-NEIBR,NEIBR
DO 520 K=-NEIBR,NEIBR
MM = MM + 1
IRX(MM) = I
IRY(MM) = I
IRZ(MM) = I
F1 = DFLOAT(I)
F2 = DFLOAT(J)
F3 = DFLOAT(K)
RX(MM) = ALTV(1,1)*F1 + ALTV(1,2)*F2 + ALTV(1,3)*F3
RY(MM) = ALTV(2,1)*F1 + ALTV(2,2)*F2 + ALTV(2,3)*F3
RZ(MM) = ALTV(3,1)*F1 + ALTV(3,2)*F2 + ALTV(3,3)*F3
RR(MM) = SQRT( RX(MM)**2 + RY(MM)**2 + RZ(MM)**2 )
520 CONTINUE
510 CONTINUE
500 CONTINUE
IF(MM.NE.NEIBRD) THEN
WRITE(6,*) ' MM.NE.NEIBRD MM,NEIBRD =',MM,NEIBRD
END IF
C CALL HPSORT(NEIBRD,IRX,IRY,IRZ,RX,RY,RZ,RR) <------- STRESS
IF(IPRI.GE.2) THEN
WRITE(6,*) ' TRANSLATIONAL VECTOR'
WRITE(6,600) (I,IRX(I),IRY(I),IRZ(I),
& RX(I) ,RY(I) ,RZ(I) ,RR(I),I=1,NEIBRD)
END IF
600 FORMAT((' ',4I5,4F12.6))
C WRITE (6,*) 'UNIVOL(IN LATTIC) = ',UNIVOL
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE GSTEP1(IPRE)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
KMAX=KNX+1
KMAY=KNY+1
KMAZ=KNZ+1
IF (IPRE.NE.1) THEN
KMAX=KNX-1
KMAY=KNY-1
KMAZ=KNZ-1
END IF
KMAX2=2*KMAX+1
KMAY2=2*KMAY+1
KMAZ2=2*KMAZ+1
C ## LARGEST WAVE VECTOR; CUT OFF WAVE NUMBER FOR WAVE FUNCTION ##
GMAX2=GMAX*2.0D0
C //////////////////////////////////////////////
C // GENERATION OF RECIPROCAL LATTICE VECTORS //
C // STEP 1 //
C //////////////////////////////////////////////
B1X = RLTV(1,1)
B1Y = RLTV(2,1)
B1Z = RLTV(3,1)
B2X = RLTV(1,2)
B2Y = RLTV(2,2)
B2Z = RLTV(3,2)
B3X = RLTV(1,3)
B3Y = RLTV(2,3)
B3Z = RLTV(3,3)
MM=0
DO 20 I1=1,KMAX2
F1=FLOAT(I1-KMAX-1)
DO 20 I2=1,KMAY2
F2=FLOAT(I2-KMAY-1)
DO 20 I3=1,KMAZ2
F3=FLOAT(I3-KMAZ-1)
FX= F1*B1X + F2*B2X + F3*B3X
FY= F1*B1Y + F2*B2Y + F3*B3Y
FZ= F1*B1Z + F2*B2Z + F3*B3Z
FR=SQRT(FX*FX+FY*FY+FZ*FZ)
IF(FR.LE.GMAX2) THEN
MM=MM+1
IG1(MM)=I1-KMAX-1
IG2(MM)=I2-KMAY-1
IG3(MM)=I3-KMAZ-1
GX(MM)=FX
GY(MM)=FY
GZ(MM)=FZ
GR(MM)=FR
END IF
20 CONTINUE
WRITE(6,*) 'NUMBER OF GENERATED VECTORS=',MM
IF(MM.NE.KG) WRITE(6,*) '**WARNING] KG MUST BE EQUAL TO ',MM
C-----------------------------------------------------------------------
KG=MM
IXMA=0
IYMA=0
IZMA=0
DO 643 I=1,KG
IXMA= MAX(IXMA,ABS(IG1(I)))
IYMA= MAX(IYMA,ABS(IG2(I)))
IZMA= MAX(IZMA,ABS(IG3(I)))
643 CONTINUE
IXMA = IXMA +1
IYMA = IYMA +1
IZMA = IZMA +1
WRITE(6,*) '**KX1 SHOULD BE =',IXMA
WRITE(6,*) '**KY1 SHOULD BE =',IYMA
WRITE(6,*) '**KZ1 SHOULD BE =',IZMA
C
KX1=IXMA
KY1=IYMA
KZ1=IZMA
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE GSTSCA(IPRE,EPP)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
REAL EPP(3,3),ERE(3,3),EPCI(3,3)
C SXX, SXY, SXZ
C SYX, SYY, SYZ
C SZX, SZY, SZZ
C X Y Z
C 1 ALTV(1,1),ALTV(2,1),ALTV(3,1) AX1,AY1,AZ1
C 2 ALTV(1,2),ALTV(2,2),ALTV(3,2) AX2,AY2,AZ2
C 3 ALTV(1,3),ALTV(2,3),ALTV(3,3) AX3,AY3,AZ3
C
C CALCULATE INVERSE MATRIX OF STRESS
SXX=EPP(1,1)
SXY=EPP(1,2)
SXZ=EPP(1,3)
SYY=EPP(2,2)
SYZ=EPP(2,3)
SZZ=EPP(3,3)
C
ADJA=SXX*SYY*SZZ+2.0D0*SXY*SXZ*SYZ-SXZ*SXZ*SYY-SYZ*SYZ*SXX
& -SXY*SXY*SZZ
ADJAI=1.0D0/ADJA
C
ERE(1,1)=ADJAI*(SYY*SZZ-SYZ*SYZ)
ERE(1,2)=ADJAI*(SYZ*SXZ-SZZ*SXY)
ERE(1,3)=ADJAI*(SXY*SYZ-SXZ*SYY)
ERE(2,1)=ERE(1,2)
ERE(2,2)=ADJAI*(SZZ*SXX-SXZ*SXZ)
ERE(2,3)=ADJAI*(SXZ*SXY-SXX*SYZ)
ERE(3,1)=ERE(1,3)
ERE(3,2)=ERE(2,3)
ERE(3,3)=ADJAI*(SXX*SYY-SXY*SXY)
C CHECK!
DO 2000 I=1,3
DO 2100 J=1,3
EPCI(I,J)=EPP(I,1)*ERE(1,J)+EPP(I,2)*ERE(2,J)
& +EPP(I,3)*ERE(3,J)
2100 CONTINUE
2000 CONTINUE
C
WRITE (6,*) 'CHECK! EPP*ERE'
WRITE (6,*) EPCI(1,1),EPCI(1,2),EPCI(1,3)
WRITE (6,*) EPCI(2,1),EPCI(2,2),EPCI(2,3)
WRITE (6,*) EPCI(3,1),EPCI(3,2),EPCI(3,3)
C
DO 1000 I=1,KNG
GXX =ERE(1,1)*GX(I)+ERE(1,2)*GY(I)+ERE(1,3)*GZ(I)
GYY =ERE(2,1)*GX(I)+ERE(2,2)*GY(I)+ERE(2,3)*GZ(I)
GZZ =ERE(3,1)*GX(I)+ERE(3,2)*GY(I)+ERE(3,3)*GZ(I)
GX(I)=GXX
GY(I)=GYY
GZ(I)=GZZ
GR(I)=SQRT(GX(I)*GX(I)+GY(I)*GY(I)+GZ(I)*GZ(I))
1000 CONTINUE
C
C ERE(1,1)=1.0D0/EPP(1,1)
C ERE(2,2)=1.0D0/EPP(2,2)
C ERE(3,3)=1.0D0/EPP(3,3)
RETURN
END
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE SYMM(KOPR,NOPR,NBZTYP,NEIBRD,NFOUT,RX,RY,RZ)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DIMENSION NPPP(KNG,KO),NAAA(KATM,KO)
DIMENSION OO(3,3,KO),TAA(3,KO)
DIMENSION RX(NEIBRD),RY(NEIBRD),RZ(NEIBRD)
!XOCL SUBPROCESSOR PT(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IQ=(PT,INDEX=1:KO,PART=BAND)
!XOCL LOCAL NPPP(:,/IQ),NAAA(:,/IQ),OO(:,:,/IQ),TAA(:,/IQ)
EQUIVALENCE (NGPT,NPPP),(NAPT,NAAA),(OP,OO),(TAU,TAA)
C----------------------------------------------------------------------
IF (NBZTYP.EQ.5) THEN
A = ABS(ALTV(1,1))*2.D0
ELSE IF(NBZTYP.EQ.6.OR.NBZTYP.EQ.7) THEN
A = ABS(ALTV(1,1))
B = ABS(ALTV(2,2))
ELSE IF(NBZTYP.EQ.8 .OR. NBZTYP.EQ.10) THEN
A = ABS(ALTV(1,1))
B = ABS(ALTV(3,3))
ELSE IF(NBZTYP.EQ.9.OR.NBZTYP.GE.11) THEN
A = ABS(ALTV(1,1))
B = ABS(ALTV(2,2))
C = ABS(ALTV(3,3))
END IF
C---------------------------------------------
CALL OPGR
> (NBZTYP,KOPR,A,B,C,NFOUT,
< NOPR)
C--*--SYMMETRY PAIRS FOR G-POINTS ------------
DO 100 I = 1,KG
PX=GX(I)
PY=GY(I)
PZ=GZ(I)
!XOCL SPREAD DO /IQ
DO 110 NO=1,NOPR
FX=OO(1,1,NO)*PX +OO(1,2,NO)*PY +OO(1,3,NO)*PZ
FY=OO(2,1,NO)*PX +OO(2,2,NO)*PY +OO(2,3,NO)*PZ
FZ=OO(3,1,NO)*PX +OO(3,2,NO)*PY +OO(3,3,NO)*PZ
DO 120 J = 1,KG
FF1 = ABS(GX(J)-FX)+ABS(GY(J)-FY)+ABS(GZ(J)-FZ)
IF (FF1.LE.1.D-5) THEN
NPPP(I,NO) = J
GOTO 110
END IF
120 CONTINUE
WRITE(NFOUT,130) I,NO
STOP
110 CONTINUE
!XOCL END SPREAD
100 CONTINUE
130 FORMAT(' ','THERE IS NO PAIR FOR NG,NOP=',2I8)
C--*--SYMMETRY PAIRS FOR ATOMS
DDD = 1.0D-5
DO 200 I = 1,KATM
PX=CATX(I)
PY=CATY(I)
PZ=CATZ(I)
!XOCL SPREAD DO /IQ
DO 210 NO=1,NOPR
FX=OO(1,1,NO)*PX+OO(1,2,NO)*PY+OO(1,3,NO)*PZ+TAA(1,NO)
FY=OO(2,1,NO)*PX+OO(2,2,NO)*PY+OO(2,3,NO)*PZ+TAA(2,NO)
FZ=OO(3,1,NO)*PX+OO(3,2,NO)*PY+OO(3,3,NO)*PZ+TAA(3,NO)
DO 220 J = 1,KATM
DO 230 K = 1,NEIBRD
FFX = ABS( FX - CATX(J) - RX(K) )
FFY = ABS( FY - CATY(J) - RY(K) )
FFZ = ABS( FZ - CATZ(J) - RZ(K) )
IF( (FFX.LE.DDD).AND.
& (FFY.LE.DDD).AND.
& (FFZ.LE.DDD) ) THEN
NAAA(I,NO) = J
GOTO 210
END IF
230 CONTINUE
220 CONTINUE
WRITE(NFOUT,*) ' NO PAIR I, NO ',I,NO
STOP
210 CONTINUE
!XOCL END SPREAD
200 CONTINUE
C--*--OUTPUT
IF(IPRI.GE.2) THEN
WRITE(NFOUT,*) ' NGPT '
DO 310 I=1,20
WRITE(NFOUT,*) ' NG =',I
WRITE(NFOUT,300) (NGPT(I,J),J=1,NOPR)
310 CONTINUE
DO 320 I=KG-20,KG
WRITE(NFOUT,*) ' NG =',I
WRITE(NFOUT,300) (NGPT(I,J),J=1,NOPR)
320 CONTINUE
WRITE(NFOUT,*) ' NAPT '
DO 330 I=1,KATM
WRITE(NFOUT,*) ' NA =',I
WRITE(NFOUT,300) (NAPT(I,J),J=1,NOPR)
330 CONTINUE
300 FORMAT((8I8))
END IF
RETURN
END
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE SYMSCA(KOPR,NOPR,NBZTYP,NEIBRD,NFOUT,RX,RY,RZ,EPP)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DIMENSION TAA(3,KO)
!XOCL SUBPROCESSOR PT(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IQ=(PT,INDEX=1:KO,PART=BAND)
!XOCL LOCAL TAA(:,/IQ)
EQUIVALENCE (TAU,TAA)
DIMENSION EPP(3,3)
C----------------------------------------------------------------------
!XOCL SPREAD DO /IQ
DO 1000 I=1,NOPR
C TAU(1,I)=EPP(1,1)*TAU(1,I)+EPP(1,2)*TAU(2,I)+EPP(1,3)*TAU(3,I)
C TAU(2,I)=EPP(1,2)*TAU(1,I)+EPP(2,2)*TAU(2,I)+EPP(2,3)*TAU(3,I)
C TAU(3,I)=EPP(1,3)*TAU(1,I)+EPP(2,3)*TAU(2,I)+EPP(3,3)*TAU(3,I)
TAA(1,I)=EPP(1,1)*TAA(1,I)
TAA(2,I)=EPP(2,2)*TAA(2,I)
TAA(3,I)=EPP(3,3)*TAA(3,I)
1000 CONTINUE
!XOCL END SPREAD
C DO 1100 I=1,NOPR
C WRITE (6,600) TAU(1,I),TAU(2,I),TAU(3,I)
C 600 FORMAT(1H ,'TAU(1,2,3) = ',3F12.6)
C1100 CONTINUE
RETURN
END
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE FORCES
> (KOPR,NOPR,FORC,FORCW)
C FORCE SYMMETRIZATION
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DIMENSION NAAA(KATM,KO)
DIMENSION OO(3,3,KO)
!XOCL SUBPROCESSOR PT(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IQ=(PT,INDEX=1:KO,PART=BAND)
!XOCL LOCAL NAAA(:,/IQ),OO(:,:,/IQ)
EQUIVALENCE (NAPT,NAAA),(OP,OO)
DIMENSION FORC(KATM,3),FORCW(KATM,3)
C----------------------------------------------------------------------
DENOM = 1.0D0/FLOAT(NOPR)
DO 100 IA = 1,KATM
FORCW(IA,1) = 0.0D0
FORCW(IA,2) = 0.0D0
FORCW(IA,3) = 0.0D0
100 CONTINUE
!XOCL SPREAD DO /IQ
DO 210 IOP = 1,NOPR
DO 200 IA= 1,KATM
IAA= NAAA(IA,IOP)
F1 = FORC(IAA,1)
F2 = FORC(IAA,2)
F3 = FORC(IAA,3)
F4 = OO(1,1,IOP)*F1 +OO(2,1,IOP)*F2 +OO(3,1,IOP)*F3
F5 = OO(1,2,IOP)*F1 +OO(2,2,IOP)*F2 +OO(3,2,IOP)*F3
F6 = OO(1,3,IOP)*F1 +OO(2,3,IOP)*F2 +OO(3,3,IOP)*F3
FORCW(IA,1) = FORCW(IA,1) + F4
FORCW(IA,2) = FORCW(IA,2) + F5
FORCW(IA,3) = FORCW(IA,3) + F6
200 CONTINUE
210 CONTINUE
!XOCL END SPREAD SUM(FORCW)
DO 300 IA = 1,KATM
FORC(IA,1) = FORCW(IA,1) * DENOM
FORC(IA,2) = FORCW(IA,2) * DENOM
FORC(IA,3) = FORCW(IA,3) * DENOM
300 CONTINUE
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE GSTEPF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DO 30 I=1,KG
IF (IG1(I).LE.-1) IGF1(I)=IG1(I)+2*IFX-1
IF (IG1(I).GE. 0) IGF1(I)=IG1(I)+1
IF (IG2(I).LE.-1) IGF2(I)=IG2(I)+2*IFY-1
IF (IG2(I).GE. 0) IGF2(I)=IG2(I)+1
IF (IG3(I).LE.-1) IGF3(I)=IG3(I)+2*IFZ-1
IF (IG3(I).GE. 0) IGF3(I)=IG3(I)+1
30 CONTINUE
C --------------------TETRAHEDRON VOLUME-------------------------
FFF = RLTV(1,1)*(RLTV(2,2)*RLTV(3,3)-RLTV(3,2)*RLTV(2,3))
& + RLTV(2,1)*(RLTV(3,2)*RLTV(1,3)-RLTV(1,2)*RLTV(3,3))
& + RLTV(3,1)*(RLTV(1,2)*RLTV(2,3)-RLTV(2,2)*RLTV(1,3))
RVOL = ABS(FFF)
DO 1000 K=1,KNZ2
DO 1000 J=1,KNY2
DO 1000 I=1,KNX2
IGPO(I,J,K) = 0
1000 CONTINUE
C CALCULATION OF POINTER IGPO
DO 1100 I=1,KG
IGPO(IG1(I)+KX1,IG2(I)+KY1,IG3(I)+KZ1)= I
1100 CONTINUE
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE GSFSCA
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
C --------------------TETRAHEDRON VOLUME-------------------------
FFF = RLTV(1,1)*(RLTV(2,2)*RLTV(3,3)-RLTV(3,2)*RLTV(2,3))
& + RLTV(2,1)*(RLTV(3,2)*RLTV(1,3)-RLTV(1,2)*RLTV(3,3))
& + RLTV(3,1)*(RLTV(1,2)*RLTV(2,3)-RLTV(2,2)*RLTV(1,3))
RVOL = ABS(FFF)
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE BASNUM(QWGT)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DIMENSION QWGT(KNV3)
C //////////////////////////////
C // CUT OFF FOR WAVEFUNCTION //
C //////////////////////////////
JGB=0
JG2=0
JJT=0
GRVT=0.D0
DO 1000 I=1,KV3
DX = VX(I)
DY = VY(I)
DZ = VZ(I)
JJ = 0
KK = 0
GRVMX = 0.D0
DO 1100 J=1,KG
GRVV = SQRT((DX+GX(J))**2+(DY+GY(J))**2+(DZ+GZ(J))**2)
IF(GRVV.LE.GMAX) THEN
JJ = JJ + 1
NBASE(JJ,I)=J
IF(GRVMX.LT.GRVV) GRVMX = GRVV
END IF
IF(GRVV.LE.GMAX/2.0D0) THEN
KK = KK + 1
NBMAT(KK,I)=J
END IF
1100 CONTINUE
IBA(I) = JJ
JJT = JJT+JJ
GRVT = GRVT+GRVMX*QWGT(I)
IBA2(I)= KK
WRITE(6,1200) I,JJ,KK,GRVT
1200 FORMAT(' ',' K = ',I4,' JJ = ',I5,' KK = ',I5,' GRV = ',F12.8)
IF(JJ.GT.JGB) JGB = JJ
IF(KK.GT.JG2) JG2 = KK
1000 CONTINUE
WRITE(6,*) 'MAXIMUM NUMBER OF BASES FOR WAVE FUNCTION (KG1)=',JGB
WRITE(6,*) ' JJT = ',JJT,' MEAN GRV = ',GRVT
WRITE(6,*) 'MAXIMUM NUMBER OF BASES FOR (GMAX/2) (KG2)=',JG2
C <----- ATTN] KG1 WO SAITEIGI SHITE IRU]
KG1=JGB
NBMX = 0
DO 2000 IK = 1,KV3
IIBA = IBA(IK)
DO 2010 I= 1,IIBA
NBMX=MAX(NBASE(I,IK),NBMX)
2010 CONTINUE
2000 CONTINUE
WRITE(6,2020) NBMX
2020 FORMAT(' ','MAX NBMX = ',I6)
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE FORM
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DO 1000 IT = 1,KTYP
DO 1010 N = 1,KG
ZFM3(N,IT) = DCMPLX(0.0D0,0.0D0)
1010 CONTINUE
1000 CONTINUE
DO 2000 IA=1,KATM
DO 2010 I=1,KG
FX = GX(I)
FY = GY(I)
FZ = GZ(I)
ZFORT =CDEXP(-ZI*(CATX(IA)*FX+CATY(IA)*FY+CATZ(IA)*FZ))
*VOCL STMT,IF(90)
IF(I.LE.KNBMX) ZFM2(I,IA) = ZFORT
ZFM3(I,KFTYPE(IA)) = ZFM3(I,KFTYPE(IA)) + ZFORT
2010 CONTINUE
2000 CONTINUE
RETURN
END
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE PCC(KTPCC,RHPC,
< SCHGPC,EPC,PCM)
C SUBROUTINE FOR PARTIAL CORE CORRECTION
C MODIFIED FROM ATOMIC SEP.1990 MORIKAWA
C REF. S.G.LOUIE,S.FROYEN,AND M.L.COHEN,
C PHYS.REV.B26,1738(1982)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
INTEGER I,IA,ITY,N,KTPCC(KTYP)
REAL EPC,PCM,RS
REAL XH,H,H3,S,SCHGPC
REAL RHPC(KIZA,KTYP),CHGPC(KTYP),EPCP(KTYP),EMCP(KTYP)
DATA XH/32.0D0/
H = 1.D0/XH
H3 = H/3.D0
CO3=4.D0/3.D0
C WRITE (6,*) 'UNIVOL IN PCC =',UNIVOL
C///////////// FOURIER TRANSFORMATION OF PARTIAL CHARGE DENSITY /////
IF (ITER.EQ.0) THEN
DO 1110 ITY=1,KTYP
DO 1111 N=1,KIZA
RHPC(N,ITY)=0.0D0
1111 CONTINUE
1110 CONTINUE
C----------------------------------------
REWIND 27
READ(27,1100) (I,RAD(N),N=1,KIZA)
1100 FORMAT(3(I4,D22.14))
IF (I.NE.KIZA) THEN
WRITE (6,*) ' READ IN PCC '
STOP
END IF
WRITE(6,*) I,' RAD(KIZA)= ',RAD(KIZA)
DO 1000 ITY=1,KTYP
IF (KTPCC(ITY).EQ.0) GOTO 1000
READ(27,*) CHGPC(ITY)
READ(27,1100) (I,RHPC(N,ITY),N=1,KIZA)
IF (I.NE.KIZA) THEN
WRITE (6,*) ' READ IN PCC ITY= ',ITY
STOP
END IF
1000 CONTINUE
END IF
C--------------------------------
DO 1001 ITY=1,KTYP
WRITE(6,*) ITY,' RHPC(KIZA,ITY)= ',RHPC(KIZA,ITY)
S = 0.D0
EPCP(ITY) = 0.D0
EMCP(ITY) = 0.D0
DO 1150 N=1,KIZA
IF (RHPC(N,ITY).GE.1.D-40) THEN
RS = (3.D0/RHPC(N,ITY))**(1.D0/3.D0)*RAD(N)
EPCP(ITY) = EPCP(ITY)-(0.458D0/RS+0.44D0/(RS+7.8D0))
& *RHPC(N,ITY)*OMO(N)
EMCP(ITY) = EMCP(ITY)-(0.458D0*CO3/RS+(0.44D0*CO3*RS +
& 3.432D0)/(RS+7.8D0)**2)*RHPC(N,ITY)*OMO(N)
END IF
S = S + RHPC(N,ITY)*OMO(N)
1150 CONTINUE
S = S *H3
EPCP(ITY) = EPCP(ITY)*H3
EMCP(ITY) = EMCP(ITY)*H3
WRITE (6,*) ' CHGPC(ITY) = ',CHGPC(ITY)
IF (ABS(S-CHGPC(ITY)).GE.1.D-7) THEN
WRITE (6,*) ITY,'PARTIAL CORE CHARGE S-CHGPC='
/ ,S-CHGPC(ITY)
STOP
END IF
C------------------------------------
DO 1300 I=1,KG
RHPCG(I,ITY)=0.0D0
RRPCG(I,ITY)=0.0D0
1300 CONTINUE
DO 1200 N=1,KIZA
R=RAD(N)
FAC1=OMO(N)*RHPC(N,ITY)*H3/UNIVOL
FAC2=OMO(N)*RHPC(N,ITY)*R*H3/UNIVOL
DO 1210 I=1,KG
XX(I)=R*GR(I)
1210 CONTINUE
CALL DSJNV(0,KG,XX,Y11)
CALL DSJNV(1,KG,XX,Y22)
DO 1220 I=1,KG
RHPCG(I,ITY)=RHPCG(I,ITY)+FAC1*Y11(I)
RRPCG(I,ITY)=RRPCG(I,ITY)-FAC2*Y22(I)
1220 CONTINUE
1200 CONTINUE
1001 CONTINUE
SCHGPC = 0.D0
EPC = 0.D0
PCM = 0.0D0
DO 1410 IA = 1,KATM
IF (KTPCC(KFTYPE(IA)).EQ.0) GOTO 1410
SCHGPC = SCHGPC+CHGPC(KFTYPE(IA))
EPC = EPC +EPCP (KFTYPE(IA))
PCM = PCM +EMCP (KFTYPE(IA))
1410 CONTINUE
WRITE (6,*) 'NO-OVERLAP P.C.(USED) EPC = ',EPC
WRITE (6,*) 'NO-OVERLAP P.C.(USED) PCM = ',PCM
WRITE (6,*) 'TOTAL PARTIAL CHARGE SCHGPC = ',SCHGPC
C///////////////////////////////////////////////////////////////////
ZRHPC(1)=DCMPLX(0.D0,0.D0)
ZRRPC(1)=DCMPLX(0.D0,0.D0)
RG(1)=1.0D0
DO 1500 I=2,KG
ZRHPC(I)=DCMPLX(0.D0,0.D0)
ZRRPC(I)=DCMPLX(0.D0,0.D0)
RG(I)=1.0D0/GR(I)
1500 CONTINUE
DO 2000 IT=1,KTYP
IF (KTPCC(IT).EQ.0) GOTO 2000
DO 2100 I=1,KG
ZRHPC(I) = ZRHPC(I) + ZFM3(I,IT)*RHPCG(I,IT)
ZRRPC(I) = ZRRPC(I) + ZFM3(I,IT)*RRPCG(I,IT)*RG(I)
2100 CONTINUE
2000 CONTINUE
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE INTCHG(ALFA)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DO 128 I=1,KG
ZCHG(I) = DCMPLX(0.0D0,0.0D0)
128 CONTINUE
DO 15 IT=1,KTYP
DO 10 I=1,KG
ZCHG(I) = ZCHG(I) + ZFM3(I,IT)*
& ACHG(IT)*EXP(-GR(I)*GR(I)*0.25D0/ALFA)/UNIVOL
10 CONTINUE
15 CONTINUE
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE FORLOC(KTPCC)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DIMENSION ZFPC(6)
INTEGER KTPCC(KTYP)
C FORCE2 AND PRESSURE ADDITION SSN(KNG1)
DO 161 IA=1,KATM
DO 162 I=1,KG
FX = GX(I)
FY = GY(I)
FZ = GZ(I)
ZFORT =CDEXP(-ZI*(CATX(IA)*FX+CATY(IA)*FY+CATZ(IA)*FZ))
ZFTMP =ZFORT *PSC(I,KFTYPE(IA))*DCONJG(ZCHG(I))
ZFORC2(IA,1)=ZFORC2(IA,1)+GX(I)*ZFTMP
ZFORC2(IA,2)=ZFORC2(IA,2)+GY(I)*ZFTMP
ZFORC2(IA,3)=ZFORC2(IA,3)+GZ(I)*ZFTMP
162 CONTINUE
IF (KTPCC(KFTYPE(IA)).EQ.1) THEN
DO 1000 N=1,6
ZFPC(N)=DCMPLX(0.0D0,0.0D0)
1000 CONTINUE
DO 200 I=1,KG
C 4/20, 1999, MODIFIED SERIOUS BUG
FX = GX(I)
FY = GY(I)
FZ = GZ(I)
ZFORT =CDEXP(-ZI*(CATX(IA)*FX+CATY(IA)*FY+CATZ(IA)*FZ))
ZFTMP =ZFORT *RHPCG(I,KFTYPE(IA))*DCONJG(ZVXC(I))
ZFPC(1)=ZFPC(1)+GX(I)*ZFTMP
ZFPC(2)=ZFPC(2)+GY(I)*ZFTMP
ZFPC(3)=ZFPC(3)+GZ(I)*ZFTMP
C [COMMENT OUT](in Japanese, 10/5, 1999)
C ZFTMP =ZFORT *RHPCG(I,KFTYPE(IA))*DCONJG(ZVXCPC(I))
C ZFPC(4)=ZFPC(4)+GX(I)*ZFTMP
C ZFPC(5)=ZFPC(5)+GY(I)*ZFTMP
C ZFPC(6)=ZFPC(6)+GZ(I)*ZFTMP
200 CONTINUE
ZFORC2(IA,1)=ZFORC2(IA,1)+ZFPC(1)-ZFPC(4)
ZFORC2(IA,2)=ZFORC2(IA,2)+ZFPC(2)-ZFPC(5)
ZFORC2(IA,3)=ZFORC2(IA,3)+ZFPC(3)-ZFPC(6)
C IF(MOD(ITER,50).EQ.0) THEN
C WRITE(6,1100) IA,ZFPC(1),ZFPC(2),ZFPC(3),
C & ZFPC(4),ZFPC(5),ZFPC(6)
C 1100 FORMAT(I4,6D12.4,/,4X,6D12.4)
C END IF
END IF
C COMBINE THE FORCE1 AND FORCE2.
161 CONTINUE
DO 261 IA=1,KATM
ZFORC2(IA,1)=UNIVOL*ZI*ZFORC2(IA,1) + DCMPLX(FFF1(IA,1),0.0D0)
ZFORC2(IA,2)=UNIVOL*ZI*ZFORC2(IA,2) + DCMPLX(FFF1(IA,2),0.0D0)
ZFORC2(IA,3)=UNIVOL*ZI*ZFORC2(IA,3) + DCMPLX(FFF1(IA,3),0.0D0)
261 CONTINUE
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE CONV2(ETONEW,ETOOLD,IMD,ITERMD,CMIX)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
IF (ITER.EQ.IMD+1) WRITE (6,*) 'MIXING RATE = ',CMIX
IF (ITER.GE.IMD+1) PMIX(ITER)=CMIX
ZFFF=DCMPLX(0.0D0,0.0D0)
PPMIX = PMIX(ITER)
QQMIX = 1.D0 - PPMIX
DO 1100 I=1,KG
ZUU1=ZCHG(I)-ZCHGO(I)
ZFFF=ZFFF+DCONJG(ZUU1)*ZUU1
ZCHO1 = ZCHGO(I)
ZCHGO(I) = ZCHG(I)
ZCHG(I) = ZCHO1*PPMIX + ZCHGO(I)*QQMIX
1100 CONTINUE
DCHG=DSQRT(DREAL(ZFFF))
5000 CONTINUE
IF(ITER.NE.1) THEN
EPP2=ABS(ETONEW-ETOOLD)
IF(ETOOLD.LT.ETONEW) THEN
C-------------------------------WRITE(6,*) '**********************'
WRITE(6,*) '>> ETOOLD.LT.ETONEW <<'
C-------------------------------WRITE(6,*) '**********************'
WRITE(6,*) 'ETONEW-ETTOOLD=',
& ETONEW-ETOOLD
EPP2=1.0D0
END IF
ETOOLD=ETONEW
ELSE
EPP2 =1.0D5
ETOOLD=ETONEW
END IF
EDEV=27200.0D0*EPP2
WRITE(6,610) ITER,EPP2,EDEV,DCHG
610 FORMAT(1H ,' ITER=',I3,' ET(H)=',D15.7,' ET(M)=',F10.6,' DC='
& ,D15.7)
IF( ITER.LE.ITEMAX ) THEN
C IF(EPP2.LE.CONV) GOTO 2000
IF(ITER.NE.ITEMAX) THEN
IFLAG=1
GO TO 9999
END IF
END IF
WRITE(6,*) '////////////////////////////////////////////'
WRITE(6,*) '// UNLUCKY! CONVERGENCE WAS NOT ACHIEVED! //'
WRITE(6,*) '////////////////////////////////////////////'
GO TO 2001
2000 CONTINUE
WRITE(6,*) '///////////////////////////////////////////////'
WRITE(6,*) '// CONGRATULATION! CONVERGENCE WAS ACHIEVED! //'
WRITE(6,*) '///////////////////////////////////////////////'
2001 CONTINUE
IFLAG=0
9999 CONTINUE
RETURN
END
C================================================================
SUBROUTINE EVIN(IREC8,IREC9,EPC)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DIMENSION DVS(KIZA),DVP(KIZA),DVD(KIZA)
REWIND 28
REWIND 29
REWIND 30
REWIND 31
REWIND 32
REWIND 33
REWIND 34
REWIND 35
REWIND 36
REWIND 37
READ (28,*) ITER,EPC,WS,WP,WD
READ (29) SNL
READ (30) SNL2
READ (31) SNL3
READ (32) ZAJ
READ (33) ZFBB,ZFBB2
READ (34) EKO
READ (35) ZCHG
READ (36) VEL
READ (37) ZRHPC
WRITE (6,*) 'EVIN ITER = ',ITER,EPC
DO 1000 I=NBD1,NBD2
DO 1010 K=1,KNV3
WRITE (6,600) EKO(I,K)
600 FORMAT(1H ,4F12.6)
1010 CONTINUE
1000 CONTINUE
C NSEK = 2 BECAUSE DVS,DVP INCLUDES RAD(N).
NSEK=2
REWIND 15
DO 100 ITY = 1,KTYP
IF (NLSPD(ITY).EQ.1) THEN
READ (15,*) MESHR,NMES,DX,WS(ITY),WP(ITY)
& ,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
RETURN
END
C================================================================
SUBROUTINE EVOUT(IREC8,IREC9,ALFA,EPC)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
REWIND 11
REWIND 20
REWIND 28
REWIND 29
REWIND 30
REWIND 31
REWIND 32
REWIND 33
REWIND 34
REWIND 35
REWIND 36
REWIND 37
WRITE (20,300) EF
300 FORMAT(1H ,'E-K CURVE EF =',F12.6)
DO 100 NNN=1,KV3
WRITE(20,303) VX(NNN),VY(NNN),VZ(NNN)
303 FORMAT((1H ,3(F10.6,2X)))
WRITE(20,302) (EKO(I,NNN),I=1,KEG)
302 FORMAT((1H ,5(F10.6,2X)))
100 CONTINUE
C
WRITE (28,*) ITER,EPC,WS,WP,WD
WRITE (29) SNL
WRITE (30) SNL2
WRITE (31) SNL3
WRITE (32) ZAJ
WRITE (33) ZFBB,ZFBB2
WRITE (34) EKO
WRITE (35) ZCHG
WRITE (36) VEL
WRITE (37) ZRHPC
C'''''WRITE 11!'''''''''''''''''''''''''''''''''''''''''''''''''''
WRITE(11,399) ITEMAX,PINIT,CONV,GMAX,ICONT
399 FORMAT(' ',I6,F8.4,D12.4,F8.4,I4)
DO 400 I=1,3
WRITE(11,500) ALTV(1,I),ALTV(2,I),ALTV(3,I)
400 CONTINUE
WRITE(11,*) KCOTYP,' COORDINATES 0:NORMALIZED 1:CARTESIAN '
IF(KCOTYP.EQ.0) THEN
DO 410 IA=1,KATM
PPOS(IA) = ALINV(1,1)*CATX(IA)
& + ALINV(1,2)*CATY(IA) + ALINV(1,3)*CATZ(IA)
QPOS(IA) = ALINV(2,1)*CATX(IA)
& + ALINV(2,2)*CATY(IA) + ALINV(2,3)*CATZ(IA)
RPOS(IA) = ALINV(3,1)*CATX(IA)
& + ALINV(3,2)*CATY(IA) + ALINV(3,3)*CATZ(IA)
410 CONTINUE
WRITE(11,500)(PPOS(IA),QPOS(IA),RPOS(IA),IA=1,KATM)
ELSE
WRITE(11,500)(CATX(IA),CATY(IA),CATZ(IA),IA=1,KATM)
END IF
NNATM=0
DO 420 IT=1,KTYP
NNATM=NNATM+IATOM(IT)
WRITE(11,*) IATOM(IT),NLSPD(IT)
WRITE(11,502) AICHG(NNATM),ALFA
WRITE(11,502) ACHG(IT),AC(IT,1),AC(IT,2),BC(IT,1),BC(IT,2)
420 CONTINUE
500 FORMAT(3(F20.10))
502 FORMAT(5(F15.8))
C'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
RETURN
END
C----------------------------------------------------------------
SUBROUTINE EVOU2(ALFA)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
REWIND 11
REWIND 20
WRITE (20,300) EF
300 FORMAT(1H ,'E-K CURVE EF =',F12.6)
DO 100 NNN=1,KV3
WRITE(20,303) VX(NNN),VY(NNN),VZ(NNN)
303 FORMAT((1H ,3(F10.6,2X)))
WRITE(20,302) (EKO(I,NNN),I=1,KEG)
302 FORMAT((1H ,5(F10.6,2X)))
100 CONTINUE
C'''''WRITE 11!'''''''''''''''''''''''''''''''''''''''''''''''''''
WRITE(11,399) ITEMAX,PINIT,CONV,GMAX,ICONT
399 FORMAT(' ',I6,F8.4,D12.4,F8.4,I4)
DO 400 I=1,3
WRITE(11,500) ALTV(1,I),ALTV(2,I),ALTV(3,I)
400 CONTINUE
WRITE(11,*) KCOTYP,' COORDINATES 0:NORMALIZED 1:CARTESIAN '
IF(KCOTYP.EQ.0) THEN
DO 410 IA=1,KATM
PPOS(IA) = ALINV(1,1)*CATX(IA)
& + ALINV(1,2)*CATY(IA) + ALINV(1,3)*CATZ(IA)
QPOS(IA) = ALINV(2,1)*CATX(IA)
& + ALINV(2,2)*CATY(IA) + ALINV(2,3)*CATZ(IA)
RPOS(IA) = ALINV(3,1)*CATX(IA)
& + ALINV(3,2)*CATY(IA) + ALINV(3,3)*CATZ(IA)
410 CONTINUE
WRITE(11,500)(PPOS(IA),QPOS(IA),RPOS(IA),IA=1,KATM)
ELSE
WRITE(11,500)(CATX(IA),CATY(IA),CATZ(IA),IA=1,KATM)
END IF
NNATM=0
DO 420 IT=1,KTYP
NNATM=NNATM+IATOM(IT)
WRITE(11,*) IATOM(IT),NLSPD(IT)
WRITE(11,502) AICHG(NNATM),ALFA
WRITE(11,502) ACHG(IT),AC(IT,1),AC(IT,2),BC(IT,1),BC(IT,2)
420 CONTINUE
500 FORMAT(3(F20.10))
502 FORMAT(5(F15.8))
C'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
RETURN
END
C------------------------------------------------------------
SUBROUTINE XCFFT(III,SCHGPC)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
C-----ARRAYS FOR MFFT--------------------------------------------
DIMENSION IWL(8*IFX2+28),IWM(8*IFY2+28),IWN(8*IFZ2+28)
& ,IWORK(2*IFX2)
DIMENSION ZRB(IFX2,IFY2,IFZ2)
C============================================================
EQUIVALENCE (ZV1D(1),ZRB(1,1,1))
C============================================================
KFT1 = IFX2
KFT2 = IFY2
KFT3 = IFZ2
KSUM=KFT1*KFT2*KFT3
KVOL=(KFT1-1)*KFT2*KFT3
THIRD=1.D0/3.D0
ALFF=0.7D0
CO1=1.5D0*ALFF*(3.D0/PAI)**THIRD
CO2=(0.75D0/PAI)**THIRD
CO3=4.D0/3.D0
VCEL=ABS(UNIVOL/FLOAT(KVOL))
C*****---- SET UP C3FFT (FAST FOURIER TRANSFORMATION) -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,0,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(SET UP)] IERR = ',IERR
STOP
END IF
DO 111 I=1,KSUM
ZV1D(I)=DCMPLX(0.0D0,0.0D0)
111 CONTINUE
IF (III.EQ.3) THEN
DO 110 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZRB(L1,L2,L3)= ZRHPC(I)
110 CONTINUE
ELSE
DO 210 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZRB(L1,L2,L3)=ZCHG(I)+ZRHPC(I)
210 CONTINUE
END IF
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERRI)
IF (IERRI.NE.0) THEN
WRITE (6,*) ' C3FFT(INVERSE)] IERR = ',IERRI
STOP
END IF
DO 1000 I=1,KSUM
CHGB1(I) = DREAL(ZV1D(I))
1000 CONTINUE
DO 20 I=1,KSUM
IF( CHGB1(I).LE.0.0D0 ) THEN
IF ( MOD(I,IFX2).NE.0 ) THEN
C IF ( CHGB1(I).LE.-1.0D-5 ) THEN
C WRITE(6,610) I,CHGB1(I)
C 610 FORMAT(1H ,'**WARNING CHG.DEN<0.0 AT ',
C & I5,2D15.7,'***')
C END IF
CHGB1(I) = 1.D-40
ELSE
CHGB1(I) = 1.D-40
END IF
END IF
20 CONTINUE
S = 0.D0
DO 1010 I=1,KSUM
S = S+CHGB1(I)
1010 CONTINUE
S=S*VCEL - SCHGPC
IF (MOD(ITER,10).EQ.0) THEN
WRITE (6,*) 'REAL TOTAL CHARGE = ',S,' IN XCFFT'
END IF
C*****<<< WIGNER INTERPOLATION >>>
IF((III.EQ.1).OR.(III.EQ.3)) THEN
DO 1020 I=1,KSUM
RS=CO2*( 1.0D0/CHGB1(I) )**THIRD
RHO= -0.458D0*CO3/RS-(0.44D0*CO3*RS + 3.432D0)/
& (RS+7.8D0)**2
ZV1D(I)=DCMPLX(RHO,0.0D0)
1020 CONTINUE
ELSE
DO 1030 I=1,KSUM
RS=CO2*( 1.0D0/CHGB1(I) )**THIRD
RHO= -0.458D0/RS-0.44D0/(RS+7.8D0)
ZV1D(I)=DCMPLX(RHO,0.0D0)
1030 CONTINUE
END IF
C*****---- FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,-1,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(FFT)] IERR = ',IERR
STOP
END IF
DENOM=1.0D0/DBLE(KVOL)
C----------------------------------------------------
IF (III.EQ.1) THEN
DO 120 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZVXC(I)=ZRB(L1,L2,L3)*DENOM
120 CONTINUE
ELSE IF (III.EQ.2) THEN
DO 121 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZEXC(I)=ZRB(L1,L2,L3)*DENOM
121 CONTINUE
ELSE IF (III.EQ.3) THEN
DO 122 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZVXCPC(I)=ZRB(L1,L2,L3)*DENOM
122 CONTINUE
END IF
RETURN
END
C------------------------------------------------------------
SUBROUTINE XSTPC(III,SCHGPC)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
C STRESS FOR PARTIAL CORE CORRECTION
C-----ARRAYS FOR MFFT--------------------------------------------
DIMENSION IWL(8*IFX2+28),IWM(8*IFY2+28),IWN(8*IFZ2+28)
& ,IWORK(2*IFX2)
DIMENSION ZRB(IFX2,IFY2,IFZ2),ZRX(KFFT)
& ,ZRY(KFFT),ZRZ(KFFT)
C============================================================
EQUIVALENCE (ZV1D(1),ZRB(1,1,1))
C============================================================
ICOUNT=0
KFT1 = IFX2
KFT2 = IFY2
KFT3 = IFZ2
KSUM=KFT1*KFT2*KFT3
KVOL=(KFT1-1)*KFT2*KFT3
THIRD=1.D0/3.D0
ALFF=0.7D0
CO1=1.5D0*ALFF*(3.D0/PAI)**THIRD
CO2=(0.75D0/PAI)**THIRD
CO3=4.D0/3.D0
VCEL=ABS(UNIVOL/FLOAT(KVOL))
C*****---- SET UP C3FFT (FAST FOURIER TRANSFORMATION) -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,0,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(SET UP)] IERR = ',IERR
STOP
END IF
C----- X ---------------------------
DO 111 I=1,KSUM
ZV1D(I)=DCMPLX(0.0D0,0.0D0)
111 CONTINUE
DO 110 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZRB(L1,L2,L3)= ZRRPC(I)*GX(I)*GX(I)
110 CONTINUE
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERRI)
IF (IERRI.NE.0) THEN
WRITE (6,*) ' C3FFT(INVERSE)] IERR = ',IERRI
STOP
END IF
DO 1005 I=1,KSUM
ZRX(I) = ZV1D(I)
1005 CONTINUE
C----- Y ---------------------------
DO 113 I=1,KSUM
ZV1D(I)=DCMPLX(0.0D0,0.0D0)
113 CONTINUE
DO 134 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZRB(L1,L2,L3)= ZRRPC(I)*GY(I)*GY(I)
134 CONTINUE
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERRI)
IF (IERRI.NE.0) THEN
WRITE (6,*) ' C3FFT(INVERSE)] IERR = ',IERRI
STOP
END IF
DO 1006 I=1,KSUM
ZRY(I) = ZV1D(I)
1006 CONTINUE
C----- Z ---------------------------
DO 114 I=1,KSUM
ZV1D(I)=DCMPLX(0.0D0,0.0D0)
114 CONTINUE
DO 124 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZRB(L1,L2,L3)= ZRRPC(I)*GZ(I)*GZ(I)
124 CONTINUE
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERRI)
IF (IERRI.NE.0) THEN
WRITE (6,*) ' C3FFT(INVERSE)] IERR = ',IERRI
STOP
END IF
DO 1007 I=1,KSUM
ZRZ(I) = ZV1D(I)
1007 CONTINUE
C-------------------------------
DO 112 I=1,KSUM
ZV1D(I)=DCMPLX(0.0D0,0.0D0)
112 CONTINUE
DO 210 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZRB(L1,L2,L3)=ZCHG(I)+ZRHPC(I)
210 CONTINUE
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERRI)
IF (IERRI.NE.0) THEN
WRITE (6,*) ' C3FFT(INVERSE)] IERR = ',IERRI
STOP
END IF
DO 1000 I=1,KSUM
CHGB1(I) = DREAL(ZV1D(I))
1000 CONTINUE
DO 20 I=1,KSUM
IF( CHGB1(I).LE.0.0D0 ) THEN
IF ( MOD(I,IFX2).NE.0 ) THEN
IF ( CHGB1(I).LE.-1.0D-5 ) THEN
ICOUNT=ICOUNT+1
C WRITE(6,610) I,CHGB1(I)
C 610 FORMAT(1H ,'**WARNING CHG.DEN<0.0 AT ',
C & I5,2D15.7,'***')
END IF
CHGB1(I) = 1.D-40
ELSE
CHGB1(I) = 1.D-40
END IF
END IF
20 CONTINUE
S = 0.D0
DO 1010 I=1,KSUM
S = S+CHGB1(I)
1010 CONTINUE
S=S*VCEL - SCHGPC
WRITE (6,*) 'REAL TOTAL CHARGE = ',S,' IN XSTPC'
C -- X ------------
DO 1020 I=1,KSUM
RS=CO2*( 1.0D0/CHGB1(I) )**THIRD
RH1 = -0.458D0*CO3/RS-(0.44D0*CO3*RS + 3.432D0)/
& (RS+7.8D0)**2
RH2 = -0.458D0/RS-0.44D0/(RS+7.8D0)
C
RHV(I)=RH1 - RH2
ZV1D(I)=ZRX(I)*(RHV(I))/CHGB1(I)
1020 CONTINUE
C*****---- FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,-1,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(FFT)] IERR = ',IERR
STOP
END IF
DENOM=1.0D0/DBLE(KVOL)
DO 30 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZXXX(I)=ZRB(L1,L2,L3)*DENOM
30 CONTINUE
C -- Y -------------
DO 1021 I=1,KSUM
ZV1D(I)=ZRY(I)*(RHV(I))/CHGB1(I)
1021 CONTINUE
C*****---- FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,-1,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(FFT)] IERR = ',IERR
STOP
END IF
DO 31 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZYYY(I)=ZRB(L1,L2,L3)*DENOM
31 CONTINUE
C -- Z -------------
DO 1022 I=1,KSUM
ZV1D(I)=ZRZ(I)*(RHV(I))/CHGB1(I)
1022 CONTINUE
C*****---- FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,-1,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(FFT)] IERR = ',IERR
STOP
END IF
DO 32 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZZZZ(I)=ZRB(L1,L2,L3)*DENOM
32 CONTINUE
WRITE (6,*) 'ICOUNT = ',ICOUNT,' XSTPC END'
RETURN
END
C /////////////////////////////////////
C // CALCULATION OF THE TOTAL ENERGY //
C /////////////////////////////////////
SUBROUTINE ENERGY(TOTCH,ETOT1,ETOT2,ETOTAL,SCHGPC,EPC)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
DIMENSION EKK(KEG,KNV3),OCCUU(KEG,KNV3)
!XOCL LOCAL EKK(:,/IP),OCCUU(:,/IP)
EQUIVALENCE (EKO,EKK),(OCCUP,OCCUU)
C /////////////////////////////////////////////////
C // PSEUDOPOTENTIAL-CORRECTION AND EWALD ENERGY //
C /////////////////////////////////////////////////
ETOTAL = ETOT1*TOTCH + ETOT2
C ///////////////////////////////////////////////////////////
C // BAND ENERGY(1) <- "VXC" AND "VHR" WITH "OLD CHARGE" //
C ///////////////////////////////////////////////////////////
CALL XCFFT(2,SCHGPC)
ETOTAL = ETOTAL - UNIVOL*DREAL(DCONJG(ZVXC(1))*ZCHG(1))
& + UNIVOL*DREAL(DCONJG(ZEXC(1))*(ZCHG(1)+ZRHPC(1)))
DO 10 I=2,KG
ETOTAL = ETOTAL - UNIVOL*
& ( DREAL(DCONJG(ZVXC(I))*ZCHG(I)) +
& PAI4*REAL( (DCONJG(ZCHGO(I)))*ZCHG(I) )*RGG(I) )
& + UNIVOL*
& ( DREAL(DCONJG(ZEXC(I))*(ZCHG(I)+ZRHPC(I))) +
& PAI2*REAL( (DCONJG(ZCHG(I)))*ZCHG(I) )*RGG(I) )
10 CONTINUE
C //////////////////////////////////////////////
C // BAND ENERGY(3); BRILLOUIN ZONE SUMMATION //
C //////////////////////////////////////////////
FFF = FLOAT(KV3)
TTT = 0.D0
!XOCL SPREAD DO /IP
DO 2 I=1,KV3
DO 100 IBAN=NBD1,NBD2
TTT = TTT + OCCUU(IBAN,I)*
& EKK(IBAN,I)
100 CONTINUE
2 CONTINUE
!XOCL END SPREAD SUM(TTT)
ETOTAL = ETOTAL + 2.D0*TTT/FFF
C FOR PARTIAL CORE CORRECTION 90.11.3
ETOTAL = ETOTAL-EPC <-- [Attn!]
WRITE(6,1104) ITER,ETOTAL,ETOTAL
WRITE(19,1104) ITER,ETOTAL,ETOTAL
1104 FORMAT(1H ,'TOTAL ENERGY FOR',I4,'-TH ITERATION=',
& F12.7,5X,D15.7)
RETURN
END
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
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 PSELMD
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
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 DIAGON(IREC8)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
PARAMETER(KNMT=KNG11*(KNG11+1)/2,KNM9=9*KNG11)
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
DIMENSION ZZ2(KNG1,KEG,KNV3),SSS(KNG1,KNV3,KTYP,10)
DIMENSION EKK(KEG,KNV3),OCCUU(KEG,KNV3)
!XOCL LOCAL ZZ2(:,:,/IP),SSS(:,/IP,:,:)
!XOCL LOCAL EKK(:,/IP),OCCUU(:,/IP)
EQUIVALENCE (ZAJ,ZZ2),(SNL,SSS),(EKO,EKK)
& ,(OCCUP,OCCUU)
C EIGEN-VALUE PROBLEM
DIMENSION ZAAA(KNMT)
& ,WWW(3*KNG11),ZWW(5*KNG11),ZAWORK(KNG1)
& ,IFLG(KEG)
& ,ZVN(KNG11,KEG),ZZZ(KNG11,KNG11),EG(KEG)
C & ,EVR(KNG11,KEG),EVI(KNG11,KEG),EG(KEG)
C & ,WK1(KNG11,6),WK2(8*KEG),ZWK3(KNG11,2),IWK(6*KEG+KNG11)
C & ,IFLG(KEG)
C & ,ZAWORK(KNG1)
C & ,WOK(KNM9),CCC(KNG11,KNG11)
C
DIMENSION PPM(KNG11),IJG(KNG11)
C
IF (KNG1.NE.KNG11) THEN
WRITE (6,*) 'WORNING] KNG1 IS NOT EQUAL TO KNG11'
END IF
C
DO 4 I=1,KG
ZCHGO(I) = ZCHG(I)
ZCHG(I) = DCMPLX(0.0D0,0.0D0)
4 CONTINUE
C ////////////////////////////
C // DIAGONALIZE AT K-POINT //
C ////////////////////////////
C VPP-PARALLEL START
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 100 NNN=1,KNV3
C IWRT(NNN) =NNN
AKX = VX(NNN)
AKY = VY(NNN)
AKZ = VZ(NNN)
IIKB = IBA(NNN)
IIBA = IBA2(NNN)
C PSEUDOPOTENTIAL NON-LOCAL PART MATRIX ELEMENT
DO 71 I=1,KNG11*(KNG11+1)/2
ZAAA(I)=DCMPLX(0.0D0,0.0D0)
71 CONTINUE
DO 272 I=1,KNG11
DO 274 K=1,KEG
C EVR(I,K)=0.0D0
C EVI(I,K)=0.0D0
ZVN(I,K)=CMPLX(0.0,0.0)
274 CONTINUE
DO 273 J=1,KNG11
C CCC(J,I)=0.0D0
ZZZ(J,I)=CMPLX(0.0,0.0)
273 CONTINUE
272 CONTINUE
C
DO 6 IT=1,KTYP
CS=1.0D0/(WS(IT)*UNIVOL)
CP=1.0D0/(WP(IT)*UNIVOL)
CWL(1)=CS
CWL(2)=CP
CWL(3)=CP
CWL(4)=CP
IF (NLSPD(IT).EQ.2) THEN
CD=1.0D0/(WD(IT)*UNIVOL)
CWL(5)=CD
CWL(6)=CD
CWL(7)=CD
CWL(8)=CD
CWL(9)=CD
END IF
IF (NLSPD(IT).EQ.1) THEN
LNUM = 4
ELSE
LNUM = 9
END IF
DO 5 I=1,IIBA
I1 = NBMAT(I,NNN)
L1=IG1(I1)+KX1
L2=IG2(I1)+KY1
L3=IG3(I1)+KZ1
II=I*(I-1)/2
DO 8700 L=1,LNUM
DO 7600 J =1,I
PPMT = CWL(L)*SSS(I,NNN,IT,L)*SSS(J,NNN,IT,L)
ZAAA(II+J) = ZAAA(II+J) + PPMT *
& ZFM3( IGPO( L1-IG1(NBMAT(J,NNN)),
& L2-IG2(NBMAT(J,NNN)),
& L3-IG3(NBMAT(J,NNN)) ) ,IT)
7600 CONTINUE
8700 CONTINUE
5 CONTINUE
6 CONTINUE
C MATRIX ELEMENT AND KINETIC ENERGY
DO 230 I=1,IIBA
PPM(I) = ( (AKX+GX(NBMAT(I,NNN)))**2
& + (AKY+GY(NBMAT(I,NNN)))**2
& + (AKZ+GZ(NBMAT(I,NNN)))**2 )/2.D0
230 CONTINUE
DO 222 I=1,IIBA
I1=NBMAT(I,NNN)
I2=I*(I-1)/2
L1=IG1(I1)+KX1
L2=IG2(I1)+KY1
L3=IG3(I1)+KZ1
ZAAA(I*(I+1)/2) = ZAAA(I*(I+1)/2) + PPM(I)
DO 222 J=1,I
C HARTREE, EXCHANGE AND CORE POTENTIAL
IJG(J) = ( IGPO( L1-IG1(NBMAT(J,NNN)), L2-IG2(NBMAT(J,NNN)),
& L3-IG3(NBMAT(J,NNN)) ) )
IF(IJG(J).NE.1) THEN
ZAAA(I2+J) = ZAAA(I2+J) + PAI4*ZCHGO(IJG(J))/GR(IJG(J))**2
END IF
ZAAA(I2+J) = ZAAA(I2+J) + ZVXC(IJG(J)) + ZPSCC(IJG(J))
222 CONTINUE
C=====DIAGONALIZATION SSL2
DO 7701 I=1,IIBA
MM = I*(I-1)/2
DO 7702 J=1,I
ZZZ(I,J)=ZAAA(MM+J)
IF (I.NE.J) ZZZ(J,I)=CONJG(ZZZ(I,J))
C RMAT = DREAL(ZAAA(MM+J))
C CMAT = DIMAG(ZAAA(MM+J))
C CCC(I,J) = RMAT
C IF (I.NE.J) CCC(J,I) = CMAT
7702 CONTINUE
7701 CONTINUE
C""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
KBAS = KNG11
KGB = KEG
IOPT = 2
EPS = 1.0D-15
C CALL HZES1M(ZAAA,IIBA,KGB,KGB,EPS,IOPT,EG,EVR,KBAS,EVI,IFLG
C & ,WK1,WK2,ZWK3,IWK,ICON)
C WRITE(6,*) KBAS,IIBA,KGB
C CALL DHEIG2(CCC,KBAS,IIBA,-KGB,EG,EVR,EVI,WOK,ICON)
CALL CHOBSD(ZZZ,KBAS,IIBA,EG,-KGB,ZVN,KGB,EPS,WWW,ZWW,ICON)
C""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
IF(ICON.NE.0) THEN
WRITE(6,*) '*****!! AT HZES1M, ICON=',ICON
END IF
C
DO 7164 M=1,KEG
EKK(M,NNN)=EG(M)
C WRITE(6,*) 'NNN = ',NNN,'IBAN = ',M
C WRITE(6,*) 'ENERGY = ',EKK(M,NNN)
7164 CONTINUE
C OUTPUT EIGENVECTORS ON FILE 80 (DIRECT-ACCESS) (KNG11 <---> IIBA)
C WRITE(6,*) 'ZAJ = '
DO 231 IBAN = NBD1,NBD2
DO 33 J=1,IIBA
C ZAJ(J,IBAN,NNN) = EVR(J,IBAN) + ZI*EVI(J,IBAN)
ZZ2(J,IBAN,NNN) = ZVN(J,IBAN)
C*********ZV1(J) = EVR(J,IBAN) + ZI*EVI(J,IBAN)
33 CONTINUE
DO 97 M=1,IIKB
ZAWORK(M)=0.D0
97 CONTINUE
DO 91 J=1,IIBA
DO 92 M=1,IIKB
IF (NBASE(M,NNN) .EQ. NBMAT(J,NNN)) THEN
ZAWORK(M)=ZZ2(J,IBAN,NNN)
ENDIF
92 CONTINUE
91 CONTINUE
DO 93 M=1,IIKB
ZZ2(M,IBAN,NNN)=ZAWORK(M)
93 CONTINUE
C**************DO 34 J=IIBA+1,KNG1
C ZV1(J) = DCMPLX(0.0D0,0.0D0)
C 34 CONTINUE
C IREC = KV3*(IBAN-1)+IWRT(NNN)
C IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1
C**************WRITE(80,REC=IREC) ZV1
231 CONTINUE
100 CONTINUE
!XOCL END SPREAD
C!XOCL END PARALLEL
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)
& ,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
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE KBINT
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
C K-B TYPE ADDITION
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),SS2(KNG1,KNV3,KTYP,9)
& ,SS3(KNG1,KNV3,KTYP,4)
& ,RAA(KNG1,KNV3)
!XOCL LOCAL SSS(:,/IP,:,:),SS2(:,/IP,:,:),SS3(:,/IP,:)
!XOCL LOCAL RAA(:,/IP)
EQUIVALENCE (SNL,SSS),(SNL2,SS2),(SNL3,SS3)
& ,(RAK,RAA)
C
C K-B TYPE ADDITION
C INTEGRATION OF K-B TYPE POTENTIAL
C MODIFIED TO STRESS CALCULATION FOR KB-SEPARABLE FORM
NSEK=2
CO0=DSQRT(2.0D0)
CO1=DSQRT(3.0D0)
SC=DSQRT(PAI4)
PC=DSQRT(3.0D0*PAI4)
DC=DSQRT(5.0D0*PAI4)
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,IBA(IK)
RK=SQRT((AKX+GX(NBASE(I,IK)))**2
& +(AKY+GY(NBASE(I,IK)))**2 + (AKZ+GZ(NBASE(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 1006 NNN=1,KV3
AKX = VX(NNN)
AKY = VY(NNN)
AKZ = VZ(NNN)
IIBA = IBA(NNN)
C
DO 7500 J =1,IIBA
QX(J)=AKX+GX(NBASE(J,NNN))
QY(J)=AKY+GY(NBASE(J,NNN))
QZ(J)=AKZ+GZ(NBASE(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 )
7500 CONTINUE
C REWIND 15
DO 7300 ITY=1,KTYP
C ATTN] IL=1,4 (S,P ONLY)
DO 7311 KN=1,KNG1
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
SS2(KN,NNN,ITY,1)=0.0D0
SS2(KN,NNN,ITY,2)=0.0D0
SS2(KN,NNN,ITY,3)=0.0D0
SS2(KN,NNN,ITY,4)=0.0D0
SS3(KN,NNN,ITY,1)=0.0D0
7311 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
FAC1=-OMO(N)*WVS(N,ITY)*R**(NSEK+1)*SC*DX
FAC2= OMO(N)*WVP(N,ITY)*R**(NSEK+1)*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)
CALL DSJNV(2,IIBA,X,Y3)
DO 7878 KN=1,IIBA
Y3(KN)=(Y1(KN)-2.0D0*Y3(KN))/3.0D0
7878 CONTINUE
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)
SS2(KN,NNN,ITY,1) =SS2(KN,NNN,ITY,1)+FAC1*Y2(KN)
SS2(KN,NNN,ITY,2) =SS2(KN,NNN,ITY,2)
& +FAC2*Y3(KN)*QX(KN)*RAA(KN,NNN)
SS2(KN,NNN,ITY,3) =SS2(KN,NNN,ITY,3)
& +FAC2*Y3(KN)*QY(KN)*RAA(KN,NNN)
SS2(KN,NNN,ITY,4) =SS2(KN,NNN,ITY,4)
& +FAC2*Y3(KN)*QZ(KN)*RAA(KN,NNN)
SS3(KN,NNN,ITY,1) =SS3(KN,NNN,ITY,1)+FACP*Y2(KN)
7330 CONTINUE
7320 CONTINUE
C
IF (NLSPD(ITY).EQ.2) THEN
C
DO 7314 KN=1,KNG1
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
SSS(KN,NNN,ITY,10)=0.0D0
SS2(KN,NNN,ITY,5)=0.0D0
SS2(KN,NNN,ITY,6)=0.0D0
SS2(KN,NNN,ITY,7)=0.0D0
SS2(KN,NNN,ITY,8)=0.0D0
SS2(KN,NNN,ITY,9)=0.0D0
SS3(KN,NNN,ITY,2)=0.0D0
SS3(KN,NNN,ITY,3)=0.0D0
SS3(KN,NNN,ITY,4)=0.0D0
7314 CONTINUE
C
DO 7321 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
FACD= OMO(N)*WVD(N,ITY)*R**NSEK*DC*DX
FAC1=-OMO(N)*WVS(N,ITY)*R**(NSEK+1)*SC*DX
FAC2= OMO(N)*WVP(N,ITY)*R**(NSEK+1)*PC*DX
FAC3= OMO(N)*WVD(N,ITY)*R**(NSEK+1)*DC*DX
DO 8778 KN=1,IIBA
X(KN)=AK2(KN)*R
8778 CONTINUE
CALL DSJNV(0,IIBA,X,Y1)
CALL DSJNV(1,IIBA,X,Y2)
CALL DSJNV(2,IIBA,X,Y3)
CALL DSJNV(3,IIBA,X,Y4)
DO 7788 KN=1,IIBA
YD(KN)= Y3(KN)
Y3(KN)=(Y1(KN)-2.0D0*Y3(KN))/3.0D0
Y4(KN)=(2.0D0*Y2(KN)-3.0D0*Y4(KN))/5.0D0
7788 CONTINUE
DO 7331 KN=1,IIBA
C ATTN] FOR SSL2 (IN ISSP)
SSS(KN,NNN,ITY,5) =SSS(KN,NNN,ITY,5)
& +FACD*YD(KN)*QA1(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SSS(KN,NNN,ITY,6) =SSS(KN,NNN,ITY,6)
& +FACD*YD(KN)*QA2(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SSS(KN,NNN,ITY,7) =SSS(KN,NNN,ITY,7)
& +FACD*YD(KN)*QA3(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SSS(KN,NNN,ITY,8) =SSS(KN,NNN,ITY,8)
& +FACD*YD(KN)*QA4(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SSS(KN,NNN,ITY,9) =SSS(KN,NNN,ITY,9)
& +FACD*YD(KN)*QA5(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SSS(KN,NNN,ITY,10)=SSS(KN,NNN,ITY,10)+FACD*YD(KN)
SS2(KN,NNN,ITY,5)=SS2(KN,NNN,ITY,5)
& +FAC3*Y4(KN)*QA1(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SS2(KN,NNN,ITY,6)=SS2(KN,NNN,ITY,6)
& +FAC3*Y4(KN)*QA2(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SS2(KN,NNN,ITY,7)=SS2(KN,NNN,ITY,7)
& +FAC3*Y4(KN)*QA3(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SS2(KN,NNN,ITY,8)=SS2(KN,NNN,ITY,8)
& +FAC3*Y4(KN)*QA4(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SS2(KN,NNN,ITY,9)=SS2(KN,NNN,ITY,9)
& +FAC3*Y4(KN)*QA5(KN)*RAA(KN,NNN)*RAA(KN,NNN)
SS3(KN,NNN,ITY,2)=SS3(KN,NNN,ITY,2)+CO1*FACD*QX(KN)*YD(KN)
& *RAA(KN,NNN)
SS3(KN,NNN,ITY,3)=SS3(KN,NNN,ITY,3)+CO1*FACD*QY(KN)*YD(KN)
& *RAA(KN,NNN)
SS3(KN,NNN,ITY,4)=SS3(KN,NNN,ITY,4)+CO1*FACD*QZ(KN)*YD(KN)
& *RAA(KN,NNN)
7331 CONTINUE
7321 CONTINUE
END IF
7300 CONTINUE
1006 CONTINUE
!XOCL END SPREAD
C!XOCL END PARALLEL
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE FORCE(IREC8)
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)
& ,ZZZ(KNG1,KEG,KNV3),ZFC(KEG,KNV3,KATM,10)
& ,OCCUU(KEG,KNV3)
!XOCL LOCAL SSS(:,/IP,:,:),ZZZ(:,:,/IP),ZFC(:,/IP,:,:)
!XOCL LOCAL OCCUU(:,/IP)
EQUIVALENCE (SNL,SSS),(ZAJ,ZZZ),(ZFBB,ZFC)
& ,(OCCUP,OCCUU)
DIMENSION ZFX1(KEG),ZFX2(KEG),ZFX3(KEG),ZFX4(KEG)
DIMENSION ZFY1(KEG),ZFY2(KEG),ZFY3(KEG),ZFY4(KEG)
DIMENSION ZFZ1(KEG),ZFZ2(KEG),ZFZ3(KEG),ZFZ4(KEG)
DIMENSION ZFB1(KEG),ZFB2(KEG),ZFB3(KEG),ZFB4(KEG)
C
DIMENSION ZFX5(KEG),ZFX6(KEG),ZFX7(KEG),ZFX8(KEG),ZFX9(KEG)
DIMENSION ZFY5(KEG),ZFY6(KEG),ZFY7(KEG),ZFY8(KEG),ZFY9(KEG)
DIMENSION ZFZ5(KEG),ZFZ6(KEG),ZFZ7(KEG),ZFZ8(KEG),ZFZ9(KEG)
DIMENSION ZFB5(KEG),ZFB6(KEG),ZFB7(KEG),ZFB8(KEG),ZFB9(KEG)
DIMENSION ZFBA(KEG)
C
DO 1000 IA=1,KATM
ZFORC2(IA,1)=DCMPLX(0.0D0,0.0D0)
ZFORC2(IA,2)=DCMPLX(0.0D0,0.0D0)
ZFORC2(IA,3)=DCMPLX(0.0D0,0.0D0)
1000 CONTINUE
C VPP-PAEALLEL START
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 2100 IK=1,KV3
IIBA = IBA(IK)
C*******DO 1 IBAN=NBD1,NBD2
C IREC = KV3*(IBAN-1)+IK
C IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1
C READ(80,REC=IREC) ZV1
C DO 2 I=1,IIBA
C ZEVC(I,IBAN) = ZV1(I)
C 2 CONTINUE
C***1 CONTINUE
DO 2110 IA=1,KATM
CS=1.0D0/(WS(KFTYPE(IA))*UNIVOL)
CP=1.0D0/(WP(KFTYPE(IA))*UNIVOL)
DO 2115 IBAN=NBD1,NBD2
ZFX1(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFX2(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFX3(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFX4(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFY1(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFY2(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFY3(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFY4(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFZ1(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFZ2(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFZ3(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFZ4(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFB1(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFB2(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFB3(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFB4(IBAN)=DCMPLX(0.0D0,0.0D0)
2115 CONTINUE
DO 2200 IBAN=NBD1,NBD2
DO 1510 I=1,IIBA
I1=NBASE(I,IK)
ZTMP=ZZZ(I,IBAN,IK)*DCONJG( ZFM2( I1 ,IA ) )
ZTMP1=ZTMP*SSS(I,IK,KFTYPE(IA),1)
ZTMP2=ZTMP*SSS(I,IK,KFTYPE(IA),2)
ZTMP3=ZTMP*SSS(I,IK,KFTYPE(IA),3)
ZTMP4=ZTMP*SSS(I,IK,KFTYPE(IA),4)
ZFB1(IBAN)=ZFB1(IBAN)+ZTMP1
ZFB2(IBAN)=ZFB2(IBAN)+ZTMP2
ZFB3(IBAN)=ZFB3(IBAN)+ZTMP3
ZFB4(IBAN)=ZFB4(IBAN)+ZTMP4
ZFX1(IBAN)=ZFX1(IBAN)+GX(I1)*ZTMP1
ZFX2(IBAN)=ZFX2(IBAN)+GX(I1)*ZTMP2
ZFX3(IBAN)=ZFX3(IBAN)+GX(I1)*ZTMP3
ZFX4(IBAN)=ZFX4(IBAN)+GX(I1)*ZTMP4
ZFY1(IBAN)=ZFY1(IBAN)+GY(I1)*ZTMP1
ZFY2(IBAN)=ZFY2(IBAN)+GY(I1)*ZTMP2
ZFY3(IBAN)=ZFY3(IBAN)+GY(I1)*ZTMP3
ZFY4(IBAN)=ZFY4(IBAN)+GY(I1)*ZTMP4
ZFZ1(IBAN)=ZFZ1(IBAN)+GZ(I1)*ZTMP1
ZFZ2(IBAN)=ZFZ2(IBAN)+GZ(I1)*ZTMP2
ZFZ3(IBAN)=ZFZ3(IBAN)+GZ(I1)*ZTMP3
ZFZ4(IBAN)=ZFZ4(IBAN)+GZ(I1)*ZTMP4
1510 CONTINUE
2200 CONTINUE
DO 2205 IBAN=NBD1,NBD2
ZFC(IBAN,IK,IA,1)=CS*ZFB1(IBAN)
ZFC(IBAN,IK,IA,2)=CP*ZFB2(IBAN)
ZFC(IBAN,IK,IA,3)=CP*ZFB3(IBAN)
ZFC(IBAN,IK,IA,4)=CP*ZFB4(IBAN)
ZFORC2(IA,1)=ZFORC2(IA,1)+OCCUU(IBAN,IK)*
& ( DCONJG(ZFX1(IBAN))*ZFC(IBAN,IK,IA,1)
& +DCONJG(ZFX2(IBAN))*ZFC(IBAN,IK,IA,2)
& +DCONJG(ZFX3(IBAN))*ZFC(IBAN,IK,IA,3)
& +DCONJG(ZFX4(IBAN))*ZFC(IBAN,IK,IA,4) )
ZFORC2(IA,2)=ZFORC2(IA,2)+OCCUU(IBAN,IK)*
& ( DCONJG(ZFY1(IBAN))*ZFC(IBAN,IK,IA,1)
& +DCONJG(ZFY2(IBAN))*ZFC(IBAN,IK,IA,2)
& +DCONJG(ZFY3(IBAN))*ZFC(IBAN,IK,IA,3)
& +DCONJG(ZFY4(IBAN))*ZFC(IBAN,IK,IA,4) )
ZFORC2(IA,3)=ZFORC2(IA,3)+OCCUU(IBAN,IK)*
& ( DCONJG(ZFZ1(IBAN))*ZFC(IBAN,IK,IA,1)
& +DCONJG(ZFZ2(IBAN))*ZFC(IBAN,IK,IA,2)
& +DCONJG(ZFZ3(IBAN))*ZFC(IBAN,IK,IA,3)
& +DCONJG(ZFZ4(IBAN))*ZFC(IBAN,IK,IA,4) )
2205 CONTINUE
C
IF (NLSPD(KFTYPE(IA)).EQ.2) THEN
C
CD=1.0D0/(WD(KFTYPE(IA))*UNIVOL)
DO 2135 IBAN=NBD1,NBD2
ZFX5(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFX6(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFX7(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFX8(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFX9(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFY5(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFY6(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFY7(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFY8(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFY9(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFZ5(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFZ6(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFZ7(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFZ8(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFZ9(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFB5(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFB6(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFB7(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFB8(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFB9(IBAN)=DCMPLX(0.0D0,0.0D0)
ZFBA(IBAN)=DCMPLX(0.0D0,0.0D0)
2135 CONTINUE
DO 2230 IBAN=NBD1,NBD2
DO 1513 I=1,IIBA
I1=NBASE(I,IK)
ZTMP=ZZZ(I,IBAN,IK)*DCONJG( ZFM2( I1 ,IA ) )
ZTMP5=ZTMP*SSS(I,IK,KFTYPE(IA),5)
ZTMP6=ZTMP*SSS(I,IK,KFTYPE(IA),6)
ZTMP7=ZTMP*SSS(I,IK,KFTYPE(IA),7)
ZTMP8=ZTMP*SSS(I,IK,KFTYPE(IA),8)
ZTMP9=ZTMP*SSS(I,IK,KFTYPE(IA),9)
ZTMPA=ZTMP*SSS(I,IK,KFTYPE(IA),10)
ZFB5(IBAN)=ZFB5(IBAN)+ZTMP5
ZFB6(IBAN)=ZFB6(IBAN)+ZTMP6
ZFB7(IBAN)=ZFB7(IBAN)+ZTMP7
ZFB8(IBAN)=ZFB8(IBAN)+ZTMP8
ZFB9(IBAN)=ZFB9(IBAN)+ZTMP9
ZFBA(IBAN)=ZFBA(IBAN)+ZTMPA
ZFX5(IBAN)=ZFX5(IBAN)+GX(I1)*ZTMP5
ZFX6(IBAN)=ZFX6(IBAN)+GX(I1)*ZTMP6
ZFX7(IBAN)=ZFX7(IBAN)+GX(I1)*ZTMP7
ZFX8(IBAN)=ZFX8(IBAN)+GX(I1)*ZTMP8
ZFX9(IBAN)=ZFX9(IBAN)+GX(I1)*ZTMP9
ZFY5(IBAN)=ZFY5(IBAN)+GY(I1)*ZTMP5
ZFY6(IBAN)=ZFY6(IBAN)+GY(I1)*ZTMP6
ZFY7(IBAN)=ZFY7(IBAN)+GY(I1)*ZTMP7
ZFY8(IBAN)=ZFY8(IBAN)+GY(I1)*ZTMP8
ZFY9(IBAN)=ZFY9(IBAN)+GY(I1)*ZTMP9
ZFZ5(IBAN)=ZFZ5(IBAN)+GZ(I1)*ZTMP5
ZFZ6(IBAN)=ZFZ6(IBAN)+GZ(I1)*ZTMP6
ZFZ7(IBAN)=ZFZ7(IBAN)+GZ(I1)*ZTMP7
ZFZ8(IBAN)=ZFZ8(IBAN)+GZ(I1)*ZTMP8
ZFZ9(IBAN)=ZFZ9(IBAN)+GZ(I1)*ZTMP9
1513 CONTINUE
2230 CONTINUE
DO 2235 IBAN=NBD1,NBD2
ZFC(IBAN,IK,IA,5)=CD*ZFB5(IBAN)
ZFC(IBAN,IK,IA,6)=CD*ZFB6(IBAN)
ZFC(IBAN,IK,IA,7)=CD*ZFB7(IBAN)
ZFC(IBAN,IK,IA,8)=CD*ZFB8(IBAN)
ZFC(IBAN,IK,IA,9)=CD*ZFB9(IBAN)
ZFC(IBAN,IK,IA,10)=CD*ZFBA(IBAN)
C 4/20, 1999, MODIFIED ZFX,Y,Z1 ---> ZFX,Y,Z5
ZFORC2(IA,1)=ZFORC2(IA,1)+OCCUU(IBAN,IK)*
& ( DCONJG(ZFX5(IBAN))*ZFC(IBAN,IK,IA,5)
& +DCONJG(ZFX6(IBAN))*ZFC(IBAN,IK,IA,6)
& +DCONJG(ZFX7(IBAN))*ZFC(IBAN,IK,IA,7)
& +DCONJG(ZFX8(IBAN))*ZFC(IBAN,IK,IA,8)
& +DCONJG(ZFX9(IBAN))*ZFC(IBAN,IK,IA,9) )
ZFORC2(IA,2)=ZFORC2(IA,2)+OCCUU(IBAN,IK)*
& ( DCONJG(ZFY5(IBAN))*ZFC(IBAN,IK,IA,5)
& +DCONJG(ZFY6(IBAN))*ZFC(IBAN,IK,IA,6)
& +DCONJG(ZFY7(IBAN))*ZFC(IBAN,IK,IA,7)
& +DCONJG(ZFY8(IBAN))*ZFC(IBAN,IK,IA,8)
& +DCONJG(ZFY9(IBAN))*ZFC(IBAN,IK,IA,9) )
ZFORC2(IA,3)=ZFORC2(IA,3)+OCCUU(IBAN,IK)*
& ( DCONJG(ZFZ5(IBAN))*ZFC(IBAN,IK,IA,5)
& +DCONJG(ZFZ6(IBAN))*ZFC(IBAN,IK,IA,6)
& +DCONJG(ZFZ7(IBAN))*ZFC(IBAN,IK,IA,7)
& +DCONJG(ZFZ8(IBAN))*ZFC(IBAN,IK,IA,8)
& +DCONJG(ZFZ9(IBAN))*ZFC(IBAN,IK,IA,9) )
2235 CONTINUE
END IF
C
2110 CONTINUE
2100 CONTINUE
!XOCL END SPREAD SUM(ZFORC2)
C!XOCL END PARALLEL
ZCCC = 4.D0*ZI*RVOL/(2.D0*PAI)**3/FLOAT(KV3)
C---------------------------------------------------------
DO 1501 IA=1,KATM
ZFORC2(IA,1)=ZCCC*DIMAG( ZFORC2(IA,1) )
ZFORC2(IA,2)=ZCCC*DIMAG( ZFORC2(IA,2) )
ZFORC2(IA,3)=ZCCC*DIMAG( ZFORC2(IA,3) )
1501 CONTINUE
3000 FORMAT(I4,6F12.8)
RETURN
END
C^^^^^1992 1/7 FOR STRESS ^^^^^^^^^^^
SUBROUTINE FORZFB
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)
& ,ZZZ(KNG1,KEG,KNV3),ZFC(KEG,KNV3,KATM,10)
!XOCL LOCAL SSS(:,/IP,:,:),ZZZ(:,:,/IP),ZFC(:,/IP,:,:)
EQUIVALENCE (SNL,SSS),(ZAJ,ZZZ),(ZFBB,ZFC)
C
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 1000 IA=1,KATM
IF (NLSPD(KFTYPE(IA)).EQ.1) THEN
LNUM = 4
ELSE
LNUM = 9
END IF
DO 1003 L=1,LNUM
DO 1001 IK=1,KV3
DO 1002 IBAN=NBD1,NBD2
ZFC(IBAN,IK,IA,L) =DCMPLX(0.0D0,0.0D0)
1002 CONTINUE
1001 CONTINUE
1003 CONTINUE
1000 CONTINUE
!XOCL END SPREAD
C!XOCL END PARALLEL
C VPP-PARALLEL START
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 2100 IK=1,KV3
DO 2110 IA=1,KATM
CS=1.0D0/(WS(KFTYPE(IA))*UNIVOL)
CP=1.0D0/(WP(KFTYPE(IA))*UNIVOL)
CWL(1)=CS
CWL(2)=CP
CWL(3)=CP
CWL(4)=CP
C
IF (NLSPD(KFTYPE(IA)).EQ.1) THEN
LNUM = 4
ELSE
LNUM = 9
CD=1.0D0/(WD(KFTYPE(IA))*UNIVOL)
CWL(5)=CD
CWL(6)=CD
CWL(7)=CD
CWL(8)=CD
CWL(9)=CD
END IF
DO 3222 L=1,LNUM
DO 3200 IBAN=NBD1,NBD2
C 1991 11/28 I ---> I1
DO 3510 I=1,IBA(IK)
I1 = NBASE(I,IK)
L1 = IG1(I1)+KX1
L2 = IG2(I1)+KY1
L3 = IG3(I1)+KZ1
ZTMP=ZZZ(I,IBAN,IK)*DCONJG( ZFM2( I1 ,IA ) )
ZFC(IBAN,IK,IA,L)=ZFC(IBAN,IK,IA,L)+
& ZTMP*SSS(I,IK,KFTYPE(IA),L)
3510 CONTINUE
ZFC(IBAN,IK,IA,L)=CWL(L)*ZFC(IBAN,IK,IA,L)
3200 CONTINUE
3222 CONTINUE
C
2110 CONTINUE
2100 CONTINUE
!XOCL END SPREAD
C!XOCL END PARALLEL
C
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE CHAVER(IREC8,IREC9,KOPR,NOPR,KBZTYP)
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 ZZZ(KNG1,KEG,KNV3),OCCUU(KEG,KNV3)
!XOCL LOCAL ZZZ(:,:,/IP),OCCUU(:,/IP)
EQUIVALENCE (ZAJ,ZZZ),(OCCUP,OCCUU)
C-----ARRAYS FOR MFFT--------------------------------------------
DIMENSION ZC3D1(IFX2,IFY2,IFZ2)
DIMENSION IWL(8*IFX2+28),IWM(8*IFY2+28),IWN(8*IFZ2+28)
& ,IWORK(2*IFX2)
C================================================================
EQUIVALENCE (ZC1D(1),ZC3D1(1,1,1))
C================================================================
C REWIND 90
KFT1 = IFX2
KFT2 = IFY2
KFT3 = IFZ2
KSUM=KFT1*KFT2*KFT3
KVOL=(KFT1-1)*KFT2*KFT3
CCC = 2.D0*RVOL/(2.D0*PAI)**3/FLOAT(KV3)
KIMG=1
C*****---- SET UP C3FFT (FAST FOURIER TRANSFORMATION) -----
CALL C3FFT(ZC3D1,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,0,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(SET UP)] IERR = ',IERR
STOP
END IF
DO 7701 I=1,KSUM
CHGB1(I)=0.0D0
7701 CONTINUE
C VPP-PARALLEL START
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 7820 IK=1,KV3
DO 7810 IBAN=NBD1,NBD2
C******* BUILD UP CHARGE DENSITY FOR IBAN-TH BAND, IK-TH K-POINT
C*********NAC = (IBAN-1)*KV3 + IWRT(IK)
C IF (IREC8.NE.0) NAC=IREC8*NAC-IREC8+1
C*********READ(80,REC=NAC) ZV1
IIBA=IBA(IK)
DO 7830 I=1,KSUM
ZC1D(I)=DCMPLX(0.0D0,0.0D0)
7830 CONTINUE
DO 232 I=1,IIBA
I1=NBASE(I,IK)
L1=IGF1(I1)
L2=IGF2(I1)
L3=IGF3(I1)
ZC3D1(L1,L2,L3) = ZZZ(I,IBAN,IK)
C***********ZC3D1(L1,L2,L3) = ZV1(I)
232 CONTINUE
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZC3D1,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERR1)
IF (IERR1.NE.0) THEN
WRITE (6,*) ' C3FFT(IFFT C)] IERR1 = ',IERR1
STOP
END IF
C---------IREC = KV3*(IBAN-1)+IWRT(IK)
C IF (IREC.LE.IRECMX) THEN
C IF (IREC9.NE.0) IREC=IREC9*IREC-IREC9+1
C WRITE(90,REC=IREC) ZC1D
C---------END IF
DO 240 I=1,KSUM
CHGB1(I)=CHGB1(I)+OCCUU(IBAN,IK)*
& DCONJG(ZC1D(I))*ZC1D(I)
240 CONTINUE
7810 CONTINUE
7820 CONTINUE
!XOCL END SPREAD SUM(CHGB1)
C!XOCL END PARALLEL
DO 7900 I=1,KSUM
ZC1D(I)=CCC*CHGB1(I)
7900 CONTINUE
C*****---- FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZC3D1,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,-1,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(FFT)] IERR = ',IERR
STOP
END IF
DENOM=1.0D0/DBLE(KVOL)
DO 30 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZCHG(I)=ZC3D1(L1,L2,L3)*DENOM
30 CONTINUE
C---------------KNGP ---> KNG KPO ---> KOPR -----------------------
IF (KBZTYP.GE.2) THEN
CALL CHGAVR(KOPR,NOPR,KIMG)
END IF
C---------------------------------------------------------------
RETURN
END
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE CHGAVR(KOPR,NOPR,KIMG)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DIMENSION NPPP(KNG,KO)
DIMENSION TAA(3,KO)
!XOCL SUBPROCESSOR PT(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IQ=(PT,INDEX=1:KO,PART=BAND)
!XOCL LOCAL NPPP(:,/IQ),TAA(:,/IQ)
EQUIVALENCE (NGPT,NPPP),(TAU,TAA)
FI = 1.D0/DFLOAT(NOPR)
ZI = CMPLX(0.0D0,1.0D0)
C--*-- CASE FOR REAL CHARGE
DO 200 NG = 1,KG
ZXXX(NG) = CMPLX(0.D0,0.0D0)
200 CONTINUE
!XOCL SPREAD DO /IQ
DO 210 NO = 1,NOPR
TX = TAA(1,NO)
TY = TAA(2,NO)
TZ = TAA(3,NO)
DO 220 NG = 1,KG
NGP= NPPP(NG,NO)
FX = GX(NGP)
FY = GY(NGP)
FZ = GZ(NGP)
FP = FX*TX + FY*TY + FZ*TZ
FC = EXP(ZI*FP)
ZCR= ZCHG(NGP)
ZXXX(NG) = ZXXX(NG) + FC*ZCR
220 CONTINUE
210 CONTINUE
!XOCL END SPREAD SUM(ZXXX)
DO 230 NG = 1,KG
ZCHG(NG) = ZXXX(NG) *FI
230 CONTINUE
RETURN
END
C ////////////////////////////////////////////////////////////
C // EWALD SUMMATION FOR TOTAL ENERGY , FORCE (VEC VERSION) //
C // 7 OCTOBER, 1986 ---> 3/15,1990 //
C ////////////////////////////////////////////////////////////
SUBROUTINE EWVEC(ETOT2)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
PARAMETER(KEWA=KATM*(KATM+1)/2)
DIMENSION CR(KLM),CI(KLM)
DIMENSION CRX(KLM),CRY(KLM),CRZ(KLM),CLX(KLM),CLY(KLM),CLZ(KLM)
DIMENSION SSS(KLM),SSI(KLM)
DIMENSION ITYP(KATM)
DIMENSION FORC1(KATM,KATM,3)
C
CON=4.D0*PAI/UNIVOL
CC=2.D0/SQRT(PAI)
ETOT2=0.D0
IRL=KLX2*KLY2*KLZ2
IKM=KMX2*KMY2*KMZ2
MM=0
DO 2 II=1,KTYP
I1=MM+1
I2=MM+IATOM(II)
MM=I2
DO 1 JJ=I1,I2
C <------ ATTN]]] KFTYPE <------- ]]]]]
ITYP(JJ)=II
1 CONTINUE
2 CONTINUE
IF(MM.NE.KATM)
& WRITE(6,*) '****WARNING "MM.NE.KATM" IN EWALD DO-LOOP 1****'
C ************************************************************
C ** THE RESULTS SHOULD BE INDEPENDENT OF THE VALUE OF ALF. **
C ************************************************************
A1X=ALTV(1,1)
A1Y=ALTV(2,1)
A1Z=ALTV(3,1)
A2X=ALTV(1,2)
A2Y=ALTV(2,2)
A2Z=ALTV(3,2)
A3X=ALTV(1,3)
A3Y=ALTV(2,3)
A3Z=ALTV(3,3)
B1X=RLTV(1,1)
B1Y=RLTV(2,1)
B1Z=RLTV(3,1)
B2X=RLTV(1,2)
B2Y=RLTV(2,2)
B2Z=RLTV(3,2)
B3X=RLTV(1,3)
B3Y=RLTV(2,3)
B3Z=RLTV(3,3)
C
DO 110 J=1,KATM
DO 111 I=1,KATM
FORC1(I,J,1)=0.0D0
FORC1(I,J,2)=0.0D0
FORC1(I,J,3)=0.0D0
111 CONTINUE
110 CONTINUE
C
MM=0
DO 10 I1=1,KLX2
F1=FLOAT(I1-KLX-1)
DO 10 I2=1,KLY2
F2=FLOAT(I2-KLY-1)
DO 10 I3=1,KLZ2
F3=FLOAT(I3-KLZ-1)
MM=MM+1
CX(MM)= F1*A1X + F2*A2X + F3*A3X
CY(MM)= F1*A1Y + F2*A2Y + F3*A3Y
CZ(MM)= F1*A1Z + F2*A2Z + F3*A3Z
10 CONTINUE
C
MK=0
DO 20 I1=1,KMX2
F1=FLOAT(I1-KMX-1)
DO 20 I2=1,KMY2
F2=FLOAT(I2-KMY-1)
DO 20 I3=1,KMZ2
F3=FLOAT(I3-KMZ-1)
MK=MK+1
CKX(MK)= F1*B1X + F2*B2X + F3*B3X
CKY(MK)= F1*B1Y + F2*B2Y + F3*B3Y
CKZ(MK)= F1*B1Z + F2*B2Z + F3*B3Z
CKR(MK)= SQRT( CKX(MK)**2 + CKY(MK)**2 + CKZ(MK)**2 )
20 CONTINUE
C
IF (IRL.NE.MM.OR.IKM.NE.MK)
& WRITE (6,*) 'WORNING] IRL .NE. MM OR IKM .NE. MK]'
C STRESS 1
DO 1000 IS=1,6
SIGEWA(IS)=0.0D0
1000 CONTINUE
C
DO 100 I=1,KATM
DO 101 J=I,KATM
ITE=1
ALFF=ALF
77 CONTINUE
AL1=SQRT(ALFF)
AL2=0.25D0/ALFF
SUM =0.0D0
SIG1 =0.0D0
SIG2 =0.0D0
SIG3 =0.0D0
SIG4 =0.0D0
SIG5 =0.0D0
SIG6 =0.0D0
SF1X =0.0D0
SF1Y =0.0D0
SF1Z =0.0D0
SF1XI=0.0D0
SF1YI=0.0D0
SF1ZI=0.0D0
C
R1=CATX(I)-CATX(J)
R2=CATY(I)-CATY(J)
R3=CATZ(I)-CATZ(J)
C ///////////////////////////////////
C // STEP 1; REAL-SPACE SUMMATION //
C ///////////////////////////////////
DO 11 L=1,MM
SSS(L)=0.0D0
SSI(L)=0.0D0
CRX(L)=CX(L)+R1
CRY(L)=CY(L)+R2
CRZ(L)=CZ(L)+R3
CLX(L)=CX(L)-R1
CLY(L)=CY(L)-R2
CLZ(L)=CZ(L)-R3
CR(L)=SQRT((CX(L)+R1)**2 + (CY(L)+R2)**2 + (CZ(L)+R3)**2 )
CI(L)=SQRT((CX(L)-R1)**2 + (CY(L)-R2)**2 + (CZ(L)-R3)**2 )
11 CONTINUE
C
DO 12 L=1,MM
IF(CR(L).EQ.0.0D0) THEN
SUM=SUM-CC*AL1
ELSE
SUM=SUM+ (1.D0-ERF(AL1*CR(L)))/CR(L)
C-----STRESS CALCULATION
STMP = - (1.D0-ERF(AL1*CR(L)))/(CR(L)*CR(L)*CR(L))
& - CC*AL1*EXP(-ALFF*CR(L)*CR(L))/(CR(L)*CR(L))
SIG1=SIG1+STMP*CRX(L)*CRX(L)
SIG2=SIG2+STMP*CRX(L)*CRY(L)
SIG3=SIG3+STMP*CRX(L)*CRZ(L)
SIG4=SIG4+STMP*CRY(L)*CRY(L)
SIG5=SIG5+STMP*CRY(L)*CRZ(L)
SIG6=SIG6+STMP*CRZ(L)*CRZ(L)
END IF
12 CONTINUE
C-----FORCE CALCULATION (R & K-SPACE) ---------------------------------
IF (I.NE.J) THEN
DO 14 L=1,MM
IF (CR(L).NE.0.0D0)
& SSS(L)=(1.0D0-ERF(AL1*CR(L)))/(CR(L)**3)+(CC*AL1
& *EXP(-ALFF*CR(L)*CR(L)))/(CR(L)*CR(L))
IF (CI(L).NE.0.0D0)
& SSI(L)=(1.0D0-ERF(AL1*CI(L)))/(CI(L)**3)+(CC*AL1
& *EXP(-ALFF*CI(L)*CI(L)))/(CI(L)*CI(L))
14 CONTINUE
C
DO 15 L=1,MM
SF1X =SF1X +CRX(L)*SSS(L)
SF1Y =SF1Y +CRY(L)*SSS(L)
SF1Z =SF1Z +CRZ(L)*SSS(L)
SF1XI=SF1XI+CLX(L)*SSI(L)
SF1YI=SF1YI+CLY(L)*SSI(L)
SF1ZI=SF1ZI+CLZ(L)*SSI(L)
15 CONTINUE
C FORCE CALCULATION (K-SPACE) FOR OPTICAL VIBRATION
DO 24 L=1,MK
IF (CKR(L).EQ.0.0D0) GOTO 24
SSK =(CON*SIN(CKX(L)*R1+CKY(L)*R2+CKZ(L)*R3)
& *EXP(-AL2*CKR(L)*CKR(L)))/(CKR(L)*CKR(L))
SF1X =SF1X +CKX(L)*SSK
SF1Y =SF1Y +CKY(L)*SSK
SF1Z =SF1Z +CKZ(L)*SSK
SF1XI=SF1XI-CKX(L)*SSK
SF1YI=SF1YI-CKY(L)*SSK
SF1ZI=SF1ZI-CKZ(L)*SSK
24 CONTINUE
C
FORC1(I,J,1)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1X
FORC1(I,J,2)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1Y
FORC1(I,J,3)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1Z
FORC1(J,I,1)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1XI
FORC1(J,I,2)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1YI
FORC1(J,I,3)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1ZI
END IF
C /////////////////////////////////////////
C // STEP 2; RECIPROCAL-SPACE SUMMATION //
C /////////////////////////////////////////
DO 21 L=1,MK
C CX = CY = CZ =0
IF (CKR(L).EQ.0.0D0) GOTO 21
STMP=CON*EXP(-AL2*CKR(L)*CKR(L))
& *COS(CKX(L)*R1+CKY(L)*R2+CKZ(L)*R3)/(CKR(L)*CKR(L))
SUM=SUM+ STMP
SIG1=SIG1+STMP*( 2.0D0*CKX(L)*CKX(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) - 1.0D0 )
SIG2=SIG2+STMP*( 2.0D0*CKX(L)*CKY(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) )
SIG3=SIG3+STMP*( 2.0D0*CKX(L)*CKZ(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) )
SIG4=SIG4+STMP*( 2.0D0*CKY(L)*CKY(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) - 1.0D0 )
SIG5=SIG5+STMP*( 2.0D0*CKY(L)*CKZ(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) )
SIG6=SIG6+STMP*( 2.0D0*CKZ(L)*CKZ(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) - 1.0D0 )
21 CONTINUE
SUM= SUM -PAI/ALFF/UNIVOL
SIG1=SIG1 +PAI/ALFF/UNIVOL
SIG4=SIG4 +PAI/ALFF/UNIVOL
SIG6=SIG6 +PAI/ALFF/UNIVOL
C ////////////
C // STEP 3 //
C ////////////
IF(ITE.EQ.1) THEN
SUM1 =SUM
SI21 =SIG1
C------------SF2X =SF1X
C SF2Y =SF1Y
C SF2Z =SF1Z
C SF2XI=SF1XI
C SF2YI=SF1YI
C------------SF2ZI=SF1ZI
ALFF=2.D0*ALFF
ITE=2
C
C IF (MOD(ITER,ILAM).EQ.1) GOTO 77
GOTO 77
ELSE
IF(ABS(SUM-SUM1).GT.1.D-5) WRITE(6,*)
& '******WARNING SUM.NE.SUM1***********'
END IF
C WRITE (6,22) I,J,SUM1 ,SUM
C WRITE (6,22) I,J,SI21 ,SIG1
C-----WRITE (6,22) I,J,SF2X ,SF1X
C WRITE (6,22) I,J,SF2Y ,SF1Y
C WRITE (6,22) I,J,SF2Z ,SF1Z
C WRITE (6,22) I,J,SF2XI,SF1XI
C WRITE (6,22) I,J,SF2YI,SF1YI
C WRITE (6,22) I,J,SF2ZI,SF1ZI
FAC=1.D0
IF(I.EQ.J) FAC=0.5D0
ETOT2= ETOT2 + FAC*SUM*ACHG(ITYP(I))*ACHG(ITYP(J))
SIGEWA(1)=SIGEWA(1)+FAC*SIG1*ACHG(ITYP(I))*ACHG(ITYP(J))
SIGEWA(2)=SIGEWA(2)+FAC*SIG2*ACHG(ITYP(I))*ACHG(ITYP(J))
SIGEWA(3)=SIGEWA(3)+FAC*SIG3*ACHG(ITYP(I))*ACHG(ITYP(J))
SIGEWA(4)=SIGEWA(4)+FAC*SIG4*ACHG(ITYP(I))*ACHG(ITYP(J))
SIGEWA(5)=SIGEWA(5)+FAC*SIG5*ACHG(ITYP(I))*ACHG(ITYP(J))
SIGEWA(6)=SIGEWA(6)+FAC*SIG6*ACHG(ITYP(I))*ACHG(ITYP(J))
22 FORMAT(1H ,2(I5,1X),' EWALD-SUMM=',D16.8,2X,D16.8)
101 CONTINUE
100 CONTINUE
SIGEWA(1)=SIGEWA(1)/UNIVOL
SIGEWA(2)=SIGEWA(2)/UNIVOL
SIGEWA(3)=SIGEWA(3)/UNIVOL
SIGEWA(4)=SIGEWA(4)/UNIVOL
SIGEWA(5)=SIGEWA(5)/UNIVOL
SIGEWA(6)=SIGEWA(6)/UNIVOL
C
DO 150 I=1,KATM
FFF1(I,1)=0.0D0
FFF1(I,2)=0.0D0
FFF1(I,3)=0.0D0
150 CONTINUE
C
DO 200 I=1,KATM
DO 201 J=1,KATM
FFF1(I,1)=FFF1(I,1)+FORC1(I,J,1)
FFF1(I,2)=FFF1(I,2)+FORC1(I,J,2)
FFF1(I,3)=FFF1(I,3)+FORC1(I,J,3)
201 CONTINUE
200 CONTINUE
C
IF (ITER.EQ.0) THEN
WRITE(6,*) 'TOTAL ENERGY CONTRIBUTION FROM THE EWALD SUMM.....'
WRITE(6,*) ETOT2
END IF
RETURN
END
C ////////////////////////////////////////////////////////////
C // EWALD SUMMATION FOR TOTAL ENERGY , FORCE (VEC VERSION) //
C // 7 OCTOBER, 1986 ---> 3/15,1990 //
C ////////////////////////////////////////////////////////////
SUBROUTINE EWVMD(ETOT2)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
PARAMETER(KEWA=KATM*(KATM+1)/2)
DIMENSION CR(KLM),CI(KLM)
DIMENSION CRX(KLM),CRY(KLM),CRZ(KLM),CLX(KLM),CLY(KLM),CLZ(KLM)
DIMENSION SSS(KLM),SSI(KLM)
DIMENSION ITYP(KATM)
DIMENSION FORC1(KATM,KATM,3)
C
CON=4.D0*PAI/UNIVOL
CC=2.D0/SQRT(PAI)
ETOT2=0.D0
IRL=KLX2*KLY2*KLZ2
IKM=KMX2*KMY2*KMZ2
MM=0
DO 2 II=1,KTYP
I1=MM+1
I2=MM+IATOM(II)
MM=I2
DO 1 JJ=I1,I2
C <------ ATTN]]] KFTYPE <------- ]]]]]
ITYP(JJ)=II
1 CONTINUE
2 CONTINUE
IF(MM.NE.KATM)
& WRITE(6,*) '****WARNING "MM.NE.KATM" IN EWALD DO-LOOP 1****'
C
DO 110 J=1,KATM
DO 111 I=1,KATM
FORC1(I,J,1)=0.0D0
FORC1(I,J,2)=0.0D0
FORC1(I,J,3)=0.0D0
111 CONTINUE
110 CONTINUE
C
DO 100 I=1,KATM
DO 101 J=I,KATM
ITE=1
ALFF=ALF
77 CONTINUE
AL1=SQRT(ALFF)
AL2=0.25D0/ALFF
SUM =0.0D0
SF1X =0.0D0
SF1Y =0.0D0
SF1Z =0.0D0
SF1XI=0.0D0
SF1YI=0.0D0
SF1ZI=0.0D0
C
R1=CATX(I)-CATX(J)
R2=CATY(I)-CATY(J)
R3=CATZ(I)-CATZ(J)
C ///////////////////////////////////
C // STEP 1; REAL-SPACE SUMMATION //
C ///////////////////////////////////
DO 11 L=1,IRL
SSS(L)=0.0D0
SSI(L)=0.0D0
CRX(L)=CX(L)+R1
CRY(L)=CY(L)+R2
CRZ(L)=CZ(L)+R3
CLX(L)=CX(L)-R1
CLY(L)=CY(L)-R2
CLZ(L)=CZ(L)-R3
CR(L)=SQRT((CX(L)+R1)**2 + (CY(L)+R2)**2 + (CZ(L)+R3)**2 )
CI(L)=SQRT((CX(L)-R1)**2 + (CY(L)-R2)**2 + (CZ(L)-R3)**2 )
11 CONTINUE
C
DO 12 L=1,IRL
IF(CR(L).EQ.0.0D0) THEN
SUM=SUM-CC*AL1
ELSE
SUM=SUM+ (1.D0-ERF(AL1*CR(L)))/CR(L)
END IF
12 CONTINUE
C-----FORCE CALCULATION (R-SPACE) ---------------------------------
IF (I.NE.J) THEN
DO 14 L=1,IRL
IF (CR(L).NE.0.0D0)
& SSS(L)=(1.0D0-ERF(AL1*CR(L)))/(CR(L)**3)+(CC*AL1
& *EXP(-ALFF*CR(L)*CR(L)))/(CR(L)*CR(L))
IF (CI(L).NE.0.0D0)
& SSI(L)=(1.0D0-ERF(AL1*CI(L)))/(CI(L)**3)+(CC*AL1
& *EXP(-ALFF*CI(L)*CI(L)))/(CI(L)*CI(L))
14 CONTINUE
C
DO 15 L=1,IRL
SF1X =SF1X +CRX(L)*SSS(L)
SF1Y =SF1Y +CRY(L)*SSS(L)
SF1Z =SF1Z +CRZ(L)*SSS(L)
SF1XI=SF1XI+CLX(L)*SSI(L)
SF1YI=SF1YI+CLY(L)*SSI(L)
SF1ZI=SF1ZI+CLZ(L)*SSI(L)
15 CONTINUE
C FORCE CALCULATION (K-SPACE) FOR OPTICAL VIBRATION
DO 24 L=1,IKM
IF (CKR(L).EQ.0.0D0) GOTO 24
SSK =(CON*SIN(CKX(L)*R1+CKY(L)*R2+CKZ(L)*R3)
& *EXP(-AL2*CKR(L)*CKR(L)))/(CKR(L)*CKR(L))
SF1X =SF1X +CKX(L)*SSK
SF1Y =SF1Y +CKY(L)*SSK
SF1Z =SF1Z +CKZ(L)*SSK
SF1XI=SF1XI-CKX(L)*SSK
SF1YI=SF1YI-CKY(L)*SSK
SF1ZI=SF1ZI-CKZ(L)*SSK
24 CONTINUE
C
FORC1(I,J,1)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1X
FORC1(I,J,2)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1Y
FORC1(I,J,3)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1Z
FORC1(J,I,1)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1XI
FORC1(J,I,2)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1YI
FORC1(J,I,3)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1ZI
END IF
C /////////////////////////////////////////
C // STEP 2; RECIPROCAL-SPACE SUMMATION //
C /////////////////////////////////////////
DO 21 L=1,IKM
C CX = CY = CZ =0
IF (CKR(L).EQ.0.0D0) GOTO 21
SUM=SUM+ CON*EXP(-AL2*CKR(L)*CKR(L))
& *COS(CKX(L)*R1+CKY(L)*R2+CKZ(L)*R3)/(CKR(L)*CKR(L))
21 CONTINUE
SUM= SUM -PAI/ALFF/UNIVOL
C ////////////
C // STEP 3 //
C ////////////
IF(ITE.EQ.1) THEN
SUM1 =SUM
C------------SF2X =SF1X
C SF2Y =SF1Y
C SF2Z =SF1Z
C SF2XI=SF1XI
C SF2YI=SF1YI
C------------SF2ZI=SF1ZI
ALFF=2.D0*ALFF
ITE=2
ELSE
IF(ABS(SUM-SUM1).GT.1.D-5) WRITE(6,*)
& '******WARNING SUM.NE.SUM1***********'
END IF
C-----WRITE (6,22) I,J,SUM1 ,SUM
C WRITE (6,22) I,J,SF2X ,SF1X
C WRITE (6,22) I,J,SF2Y ,SF1Y
C WRITE (6,22) I,J,SF2Z ,SF1Z
C WRITE (6,22) I,J,SF2XI,SF1XI
C WRITE (6,22) I,J,SF2YI,SF1YI
C WRITE (6,22) I,J,SF2ZI,SF1ZI
FAC=1.D0
IF(I.EQ.J) FAC=0.5D0
ETOT2= ETOT2 + FAC*SUM*ACHG(ITYP(I))*ACHG(ITYP(J))
22 FORMAT(1H ,2(I5,1X),' EWALD-SUMM=',D16.8,2X,D16.8)
101 CONTINUE
100 CONTINUE
C
DO 150 I=1,KATM
FFF1(I,1)=0.0D0
FFF1(I,2)=0.0D0
FFF1(I,3)=0.0D0
150 CONTINUE
C
DO 200 I=1,KATM
DO 201 J=1,KATM
FFF1(I,1)=FFF1(I,1)+FORC1(I,J,1)
FFF1(I,2)=FFF1(I,2)+FORC1(I,J,2)
FFF1(I,3)=FFF1(I,3)+FORC1(I,J,3)
201 CONTINUE
200 CONTINUE
C
WRITE(6,*) 'ETOT2(EWALD) = ',ETOT2
RETURN
END
C-----------------------------------------------------------------------
SUBROUTINE MD(IMD,KMDTP,MDFLAG,ICAR,ITERMD,ITERMX,ISTR,ISSS
& ,DTUC,EPP)
C-----------------------------------------------------------------------
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
INTEGER KMDTP(KATM)
C MOLECULAR DYNAMICS PART (SIMPLE VERLET ALGORITHM)
DIMENSION VEPSI(6),EPP(3,3),VAVB(6),EPSIR(6)
DIMENSION ALT(3,3)
C WRITE (6,*) 'IN MD ICAR,ISTR = ',ICAR,ISTR,ISSS
TKIN=0.0D0
FCMX=0.0D0
FDEL =-5.0D-12
DO 10 I=1,KATM
IF (KMDTP(I).EQ.1) THEN
FCX = DREAL(ZFORC2(I,1))
FCY = DREAL(ZFORC2(I,2))
FCZ = DREAL(ZFORC2(I,3))
FCA = SQRT(FCX**2+FCY**2+FCZ**2)
IF (FCA.GT.FCMX) FCMX=FCA
END IF
10 CONTINUE
C WRITE (6,*) 'IN MD MDFLAG = ',MDFLAG
IF(MDFLAG.EQ.0) THEN
DO 100 I=1,KATM
IF (KMDTP(I).EQ.1) THEN
FCX = DREAL(ZFORC2(I,1))
FCY = DREAL(ZFORC2(I,2))
FCZ = DREAL(ZFORC2(I,3))
C** FCA = SQRT(FCX**2+FCY**2+FCZ**2)
C** IF (FCA.GT.FCMX) FCMX=FCA
FVX=FCX*VEL(I,1)
IF (FVX.LT.0.0D0) VEL(I,1)=0.0D0
FVY=FCY*VEL(I,2)
IF (FVY.LT.0.0D0) VEL(I,2)=0.0D0
FVZ=FCZ*VEL(I,3)
IF (FVZ.LT.0.0D0) VEL(I,3)=0.0D0
C
VEL(I,1) = VEL(I,1) + DTIO*FCX/AMION(I)
VEL(I,2) = VEL(I,2) + DTIO*FCY/AMION(I)
VEL(I,3) = VEL(I,3) + DTIO*FCZ/AMION(I)
CATX(I) = CATX(I) + DTIO*VEL(I,1)
CATY(I) = CATY(I) + DTIO*VEL(I,2)
CATZ(I) = CATZ(I) + DTIO*VEL(I,3)
END IF
100 CONTINUE
C CONSIDER BOUNDARY CONDITION
C*****IF (MOD(ITER,100).EQ.0) THEN
DO 200 I=1,KATM
P = ALINV(1,1)*CATX(I)+ALINV(1,2)*CATY(I)+ALINV(1,3)*CATZ(I)
Q = ALINV(2,1)*CATX(I)+ALINV(2,2)*CATY(I)+ALINV(2,3)*CATZ(I)
R = ALINV(3,1)*CATX(I)+ALINV(3,2)*CATY(I)+ALINV(3,3)*CATZ(I)
C IF (P.GE.1.0D0) P = P-1.0D0
C IF (P.LT.0.0D0) P = P+1.0D0
C IF (Q.GE.1.0D0) Q = Q-1.0D0
C IF (Q.LT.0.0D0) Q = Q+1.0D0
C IF (R.GE.1.0D0) R = R-1.0D0
C IF (R.LT.0.0D0) R = R+1.0D0
C* RECALCULATE X,Y,Z
CATX(I)=P*ALTV(1,1)+Q*ALTV(1,2)+R*ALTV(1,3)
CATY(I)=P*ALTV(2,1)+Q*ALTV(2,2)+R*ALTV(2,3)
CATZ(I)=P*ALTV(3,1)+Q*ALTV(3,2)+R*ALTV(3,3)
PPOS(I)=P
QPOS(I)=Q
RPOS(I)=R
200 CONTINUE
C WRITE (6,*) 'NEXT IS EPP SET]'
C------------------------------------------------------------
C CHANGE A SHAPE OF AN UNIT CELL BY USING STRESSES -
C 1993 4/9 -
C------------------------------------------------------------
IF (ICAR.EQ.2 .AND. ISTR.EQ.1) THEN
WRITE (6,*) 'NOW MD SET EPP',ITER
DO 6005 I=1,6
VAVB(I)=0.0D0
6005 CONTINUE
DO 6010 IA=1,KATM
C WRITE (6,*) 'VEL(IA) = ',IA,VEL(IA,1),VEL(IA,2),VEL(IA,3)
VAVB(1)=VAVB(1)+VEL(IA,1)*VEL(IA,1)/UNIVOL
VAVB(2)=VAVB(2)+VEL(IA,1)*VEL(IA,2)/UNIVOL
VAVB(3)=VAVB(3)+VEL(IA,1)*VEL(IA,3)/UNIVOL
VAVB(4)=VAVB(4)+VEL(IA,2)*VEL(IA,2)/UNIVOL
VAVB(5)=VAVB(5)+VEL(IA,2)*VEL(IA,3)/UNIVOL
VAVB(6)=VAVB(6)+VEL(IA,3)*VEL(IA,3)/UNIVOL
6010 CONTINUE
IF (ISSS.EQ.0) THEN
DO 6000 I=1,6
C VEPSI(I)=-DTUC*(SIGSTR(I)+VAVB(I)-PEX(I))
EPSIR(I)=-DTUC*(SIGSTR(I)+VAVB(I)-PEX(I))
C EPSIR(I)= DTUC*VEPSI(I)
6000 CONTINUE
ISSS=1
END IF
DO 6100 I=1,6
C CHKVF=VEPSI(I)*(SIGSTR(I)+VAVB(I)-PEX(I))
CHKVF=EPSIR(I)*(SIGSTR(I)+VAVB(I)-PEX(I))
IF (CHKVF.GE.0.0D0) THEN
EPSIR(I)=0.0D0
C VEPSI(I)=0.0D0
END IF
C VEPSI(I)=VEPSI(I) -DTUC*(SIGSTR(I)+VAVB(I)-PEX(I))
EPSIR(I)= -DTUC*(SIGSTR(I)+VAVB(I)-PEX(I))
C EPSIR(I)=EPSIR(I) +DTUC*VEPSI(I)
6100 CONTINUE
EPP(1,1)=EPSIR(1) + 1.0D0
EPP(1,2)=EPSIR(2)
EPP(1,3)=EPSIR(3)
EPP(2,1)=EPSIR(2)
EPP(2,2)=EPSIR(4) + 1.0D0
EPP(2,3)=EPSIR(5)
EPP(3,1)=EPSIR(3)
EPP(3,2)=EPSIR(5)
EPP(3,3)=EPSIR(6) + 1.0D0
WRITE (6,*) 'EPP(1,1) = ',EPP(1,1),VAVB(1)
WRITE (6,*) 'EPP(1,2) = ',EPP(1,2),VAVB(2)
WRITE (6,*) 'EPP(1,3) = ',EPP(1,3),VAVB(3)
WRITE (6,*) 'EPP(2,1) = ',EPP(2,1),VAVB(2)
WRITE (6,*) 'EPP(2,2) = ',EPP(2,2),VAVB(4)
WRITE (6,*) 'EPP(2,3) = ',EPP(2,3),VAVB(5)
WRITE (6,*) 'EPP(3,1) = ',EPP(3,1),VAVB(3)
WRITE (6,*) 'EPP(3,2) = ',EPP(3,2),VAVB(5)
WRITE (6,*) 'EPP(3,3) = ',EPP(3,3),VAVB(6)
5000 CONTINUE
C CHANGE AXES OF A UNIT CELL
DO 6200 I=1,3
ALT(I,1)=EPP(I,1)*ALTV(1,1)+EPP(I,2)*ALTV(2,1)
& +EPP(I,3)*ALTV(3,1)
ALT(I,2)=EPP(I,1)*ALTV(1,2)+EPP(I,2)*ALTV(2,2)
& +EPP(I,3)*ALTV(3,2)
ALT(I,3)=EPP(I,1)*ALTV(1,3)+EPP(I,2)*ALTV(2,3)
& +EPP(I,3)*ALTV(3,3)
6200 CONTINUE
C
DO 6210 I=1,3
ALTV(I,1)=ALT(I,1)
ALTV(I,2)=ALT(I,2)
ALTV(I,3)=ALT(I,3)
6210 CONTINUE
C
DO 6500 IA=1,KATM
P=PPOS(IA)
Q=QPOS(IA)
R=RPOS(IA)
CATX(IA)=P*ALTV(1,1)+Q*ALTV(1,2)+R*ALTV(1,3)
CATY(IA)=P*ALTV(2,1)+Q*ALTV(2,2)+R*ALTV(2,3)
CATZ(IA)=P*ALTV(3,1)+Q*ALTV(3,2)+R*ALTV(3,3)
6500 CONTINUE
C
END IF
C
C IF(MOD(ITER,10).EQ.0) THEN
WRITE (18,*) 'ITER = ',ITER
WRITE (18,300)(I,CATX(I),CATY(I),CATZ(I),
& PPOS(I),QPOS(I),RPOS(I) ,I=1,KATM)
C END IF
IF((FCMX.LT.FDEL*0.5D0).OR.(ITER.GE.ITERMX)) THEN
MDFLAG = 1
ICAR = 0
ITERMD = ITER
WRITE(6,*) ' ITERMD= ',ITERMD
END IF
END IF
IF((FCMX.GT.FDEL).AND.(ITER.LE.ITERMX)) THEN
MDFLAG = 0
ICAR = 2
ITERMD = 5000
END IF
300 FORMAT(' ',I3,3F15.10,3F10.6)
302 FORMAT(' ',I3,3F20.10)
WRITE(21,305) ITER,FCMX,MDFLAG
305 FORMAT(' ',I5,F12.6,3X,I3,'ITER AND FCMX')
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE MSD(IREC8,IREC9,IMD,ICAR,ISTR)
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),ZZZ(KNG1,KEG,KNV3)
& ,ZFC(KEG,KNV3,KATM,10)
DIMENSION EKK(KEG,KNV3),OCCUU(KEG,KNV3),RAA(KNG1,KNV3)
!XOCL LOCAL SSS(:,/IP,:,:),ZZZ(:,:,/IP),ZFC(:,/IP,:,:)
!XOCL LOCAL EKK(:,/IP),OCCUU(:,/IP),RAA(:,/IP)
EQUIVALENCE (SNL,SSS),(ZAJ,ZZZ),(ZFBB,ZFC)
& ,(EKO,EKK),(OCCUP,OCCUU),(RAK,RAA)
C EIGEN-VALUE PROBLEM
DIMENSION ZDEVC(KNG1,KEG)
& ,EMAS(KEG)
C-----ARRAYS FOR MFFT--------------------------------------------
DIMENSION ZEV3C(IFX2,IFY2,IFZ2),ZV3D(IFX2,IFY2,IFZ2)
DIMENSION IWL(8*IFX2+28),IWM(8*IFY2+28),IWN(8*IFZ2+28)
& ,IWORK(2*IFX2)
C================================================================
EQUIVALENCE (ZV1D(1),ZV3D(1,1,1)),(ZC1D(1),ZEV3C(1,1,1))
C==== ATTN KX1 OR KX11 ==========================================
C REWIND 90
KFT1 = IFX2
KFT2 = IFY2
KFT3 = IFZ2
KSUM=KFT1*KFT2*KFT3
KVOL=(KFT1-1)*KFT2*KFT3
IF (ITER.EQ.2) WRITE (6,*) 'DTIME = ',DTIM
C---------- MASS OF ELECTRON -----
DO 4 I=1,KG
ZCHGO(I) = ZCHG(I)
ZCHG(I) = DCMPLX(0.0D0,0.0D0)
4 CONTINUE
C FOR STRESS (CALL KBINT & FORZFB 1992 1/7)
IF(ISTR.EQ.1) THEN
WRITE (6,*) 'KBINT ITER>IMD IN MSD'
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL KBINT
C^^^^^ STRESS ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL FORZFB
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
END IF
C ////////////////////////////////////
C // CALCULATE THE VERLET ALGORITHM //
C ////////////////////////////////////
ZVP(1) = ZVXC(1) + ZPSCC(1)
DO 2001 I=2,KG
ZVP(I) = ZVXC(I) + ZPSCC(I) + PAI4*ZCHGO(I)*RGG(I)
2001 CONTINUE
C*****---- SET UP C3FFT (FAST FOURIER TRANSFORMATION) -----
CALL C3FFT(ZEV3C,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,0,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(SET UP)] IERR = ',IERR
STOP
END IF
DO 234 I=1,KSUM
ZV1D(I)=DCMPLX(0.0D0,0.0D0)
234 CONTINUE
DO 233 I=1,KG
LF1=IGF1(I)
LF2=IGF2(I)
LF3=IGF3(I)
ZV3D(LF1,LF2,LF3) = ZVP(I)
233 CONTINUE
CALL TIME(ITIME1,' ',1)
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZV3D ,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERR2)
CALL TIME(ITIME1,'C3FFT ',2)
IF (IERR2.NE.0) THEN
WRITE (6,*) ' C3FFT(IFFT V)] IERR2 = ',IERR2
STOP
END IF
C WRITE(6,*) 'ZV1D = '
C DO 27 I=1,20
C WRITE(6,*) ZV1D(I)
C 27 CONTINUE
C-----K-POINT LOOP---------------------------------------------
C VPP-PARALLEL START 2
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 100 NNN=1,KV3
C IWRT(NNN) =NNN
AKX = VX(NNN)
AKY = VY(NNN)
AKZ = VZ(NNN)
IIBA = IBA(NNN)
DO 3 IBAN=NBD1,NBD2
EMAS(IBAN)=0.0D0
3 CONTINUE
DO 2002 J=1,IIBA
X(J) = 0.0D0
VKIN =
& ( (AKX+GX(NBASE(J,NNN)))**2
& + (AKY+GY(NBASE(J,NNN)))**2
& + (AKZ+GZ(NBASE(J,NNN)))**2 )/2.D0
Y3(J) = VKIN + DREAL(ZVP(1))
2002 CONTINUE
DO 7000 ITY=1,KTYP
CS=1.0D0/(WS(ITY)*UNIVOL)
CP=1.0D0/(WP(ITY)*UNIVOL)
CWL(1)=CS
CWL(2)=CP
CWL(3)=CP
CWL(4)=CP
IF (NLSPD(ITY).EQ.1) THEN
LNUM = 4
ELSE
LNUM = 9
CD=1.0D0/(WD(ITY)*UNIVOL)
CWL(5)=CD
CWL(6)=CD
CWL(7)=CD
CWL(8)=CD
CWL(9)=CD
END IF
C
DO 7111 L=1,LNUM
DO 7010 I=1,IIBA
TMPP = CWL(L)*SSS(I,NNN,ITY,L)**2
Y3(I) = Y3(I) + TMPP*DFLOAT(IATOM(ITY))
X(I) = X(I) + TMPP*DFLOAT(IATOM(ITY))
7010 CONTINUE
7111 CONTINUE
7000 CONTINUE
C-----BAND LOOP----------------------------------------------
DO 200 IBAN=NBD1,NBD2
C*********IREC = KV3*(IBAN-1)+IWRT(NNN)
C IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1
C READ(80,REC=IREC) ZV1
C IREC = KV3*(IBAN-1)+IWRT(NNN)
C IF (IREC.LE.1000) THEN
C IF (IREC.LE.IRECMX) THEN
C IF (IREC9.NE.0) IREC=IREC9*IREC-IREC9+1
C READ(90,REC=IREC) ZRC
C ELSE
C READ(90,REC=IREC) ZRC
C---------------------------------------90.1.22 Y.M.----
DO 3231 I=1,KSUM
ZC1D(I)=CMPLX(0.0D0,0.0D0)
3231 CONTINUE
DO 3232 I=1,IIBA
I1=NBASE(I,NNN)
L1=IGF1(I1)
L2=IGF2(I1)
L3=IGF3(I1)
ZEV3C(L1,L2,L3) = ZZZ(I,IBAN,NNN)
C*******ZEV3C(L1,L2,L3) = ZV1(I)
3232 CONTINUE
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZEV3C,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERR1)
IF (IERR1.NE.0) THEN
WRITE (6,*) ' C3FFT(IFFT C)] IERR1 = ',IERR1
STOP
END IF
CXX DO 3240 I=1,KSUM
CXX ZRC(I)=ZC1D(I)
C3240 CONTINUE
C END IF
C ----- CALCULATE THE MATRIX ELEMENTS -----
DO 71 I=1,KNG1
C ZPNL(I) = - ZZZ(I,IBAN,NNN)*X(I)
Y1(I) = - DREAL(ZZZ(I,IBAN,NNN)*X(I))
Y2(I) = - DIMAG(ZZZ(I,IBAN,NNN)*X(I))
71 CONTINUE
DO 401 IA=1,KATM
C
IF (NLSPD(KFTYPE(IA)).EQ.1) THEN
LNUM = 4
ELSE
LNUM = 9
END IF
C
DO 411 L =1,LNUM
DO 400 I =1,IIBA
ZTMP = SSS(I,NNN,KFTYPE(IA),L)
& * ZFC(IBAN,NNN,IA,L)
C ZPNL(I) = ZPNL(I) + ZTMP*ZFM2(NBASE(I,NNN),IA)
Y1(I) = Y1(I) + DREAL(ZTMP*ZFM2(NBASE(I,NNN),IA))
Y2(I) = Y2(I) + DIMAG(ZTMP*ZFM2(NBASE(I,NNN),IA))
400 CONTINUE
411 CONTINUE
401 CONTINUE
DO 240 I=1,KSUM
C ZC1D(I)=ZV1D(I)*ZRC(I)
ZC1D(I)=ZV1D(I)*ZC1D(I)
240 CONTINUE
C*****---- FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZEV3C,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,-1,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(FFT)] IERR = ',IERR
STOP
END IF
DENOM=1.0D0/DBLE(KVOL)
DO 30 I=1,IIBA
I1=NBASE(I,NNN)
L1=IGF1(I1)
L2=IGF2(I1)
L3=IGF3(I1)
ZPGG=ZEV3C(L1,L2,L3)*DENOM
WDI =(Y3(I)-EKK(IBAN,NNN))
C ZROF =(ZPGG+ZPNL(I)-ZVP(1)*ZZZ(I,IBAN,NNN))
ZROF =(ZPGG+Y1(I)+ZI*Y2(I)-ZVP(1)*ZZZ(I,IBAN,NNN))
ZDEVC(I,IBAN)=-(WDI*ZZZ(I,IBAN,NNN)+ZROF)
EMAS(IBAN)=EMAS(IBAN)+DCONJG(ZZZ(I,IBAN,NNN))*ZDEVC(I,IBAN)
ESHI=DEXP(-WDI*DTIM)
ZZZ(I,IBAN,NNN)=ZZZ(I,IBAN,NNN)*ESHI+(ESHI-1.0D0)*ZROF/WDI
30 CONTINUE
C-------------------------------------------------------------------
200 CONTINUE
C >>>>> AFTER C(T+DT)=C(T)+DC(T) <<<<<
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C GRAM-SCHMIDT ORTHOGONALIZATION
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C IF (NNN.EQ.1) CALL TIME(ITIME1,' ',1)
C CALL GRAMMD(NNN,IIBA)
C IF (NNN.EQ.1) CALL TIME(ITIME1,'GRAM ',2)
DO 1000 IBAN=NBD1,NBD2-1
ZNRM=DCMPLX(0.0D0,0.0D0)
DO 1010 I=1,IIBA
ZNRM=ZNRM+DCONJG(ZZZ(I,IBAN,NNN))*ZZZ(I,IBAN,NNN)
1010 CONTINUE
DENOM=1.0D0/DSQRT(DREAL(ZNRM))
DO 1020 I=1,IIBA
ZZZ(I,IBAN,NNN)=ZZZ(I,IBAN,NNN)*DENOM
1020 CONTINUE
DO 1030 JBAN=IBAN+1,NBD2
ZNRM=DCMPLX(0.0D0,0.0D0)
DO 1040 I=1,IIBA
ZNRM=ZNRM+DCONJG(ZZZ(I,IBAN,NNN))*ZZZ(I,JBAN,NNN)
1040 CONTINUE
DO 1050 I=1,IIBA
ZZZ(I,JBAN,NNN)=ZZZ(I,JBAN,NNN)-ZNRM*ZZZ(I,IBAN,NNN)
1050 CONTINUE
1030 CONTINUE
1000 CONTINUE
ZNRM=DCMPLX(0.0D0,0.0D0)
DO 1100 I=1,IIBA
ZNRM=ZNRM+DCONJG(ZZZ(I,NBD2,NNN))*ZZZ(I,NBD2,NNN)
1100 CONTINUE
DENOM=1.0D0/DSQRT(DREAL(ZNRM))
DO 1110 I=1,IIBA
ZZZ(I,NBD2,NNN)=ZZZ(I,NBD2,NNN)*DENOM
1110 CONTINUE
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C CHECK
C**************************************************************
C IF (ITER.GT.40) THEN
C DO 700 IBAN=NBD1,NBD2
C DO 720 JBAN=NBD1,IBAN
C ZNORM=DCMPLX(0.0D0,0.0D0)
C DO 7001 I=1,IIBA
C ZNORM=ZNORM+DCONJG(ZAJ(I,IBAN,NNN))*ZAJ(I,JBAN,NNN)
C7001 CONTINUE
C RNOR=DSQRT(CDABS(ZNORM))
C IF (IBAN.EQ.JBAN) THEN
C IF(DABS(1.0D0-RNOR).GT.1.0D-7) THEN
C WRITE(6,*) 'NRM.NE.1.0 I= ',IBAN,'NRM=',RNOR
C END IF
C ELSE
C IF(RNOR.GT.1.0D-7) THEN
C WRITE(6,*) 'NRM.NE.0.0 I J= ',IBAN,JBAN
C & ,'NRM= ',RNOR
C END IF
C IFLAG=1
C END IF
C 720 CONTINUE
C 700 CONTINUE
C END IF
C******************************************************************
DO 261 IBAN=NBD1,NBD2
DO 260 I=1,IIBA
EKK(IBAN,NNN)=EKK(IBAN,NNN)
& -2.0D0*DCONJG(ZZZ(I,IBAN,NNN))*ZDEVC(I,IBAN)
260 CONTINUE
EKK(IBAN,NNN)=EKK(IBAN,NNN)+EMAS(IBAN)
C EG(IBAN)=EKK(IBAN,NNN)
261 CONTINUE
DO 262 IBAN=NBD1,NBD2-1
DO 262 JBAN=IBAN+1,NBD2
IF(EKK(JBAN,NNN).LT.EKK(IBAN,NNN)) THEN
EE =EKK(IBAN,NNN)
EKK(IBAN,NNN)=EKK(JBAN,NNN)
EKK(JBAN,NNN)=EE
DO 270 IV=1,IIBA
ZTV = ZZZ(IV,IBAN,NNN)
ZZZ(IV,IBAN,NNN) = ZZZ(IV,JBAN,NNN)
ZZZ(IV,JBAN,NNN) = ZTV
270 CONTINUE
END IF
262 CONTINUE
C DO 263 IBAN=NBD1,NBD2
C EKK(IBAN,NNN)=EG(IBAN)
C WRITE(6,*) 'NNN = ',NNN,'IBAN = ',IBAN
C WRITE(6,*) EKO(IBAN,NNN)
C 263 CONTINUE
C OUTPUT EIGENVECTORS ON FILE 80 (DIRECT-ACCESS) (KG1 <---> IIBA)
DO 231 IBAN = NBD1,NBD2
C**************IREC = KV3*(IBAN-1)+IWRT(NNN)
C IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1
C**************WRITE(80,REC=IREC) ZV1
231 CONTINUE
100 CONTINUE
!XOCL END SPREAD
C!XOCL END PARALLEL
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE STRESS(IPCC,SCHGPC,ETOT1,TOTCH,EPC,PCM
& ,KOPR,NBZTYP)
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 ZZZ(KNG1,KEG,KNV3),OCCUU(KEG,KNV3)
!XOCL LOCAL ZZZ(:,:,/IP),OCCUU(:,/IP)
EQUIVALENCE (ZAJ,ZZZ),(OCCUP,OCCUU)
DIMENSION OO(3,3,KO)
!XOCL INDEX PARTITION IQ=(PS,INDEX=1:KO,PART=BAND)
!XOCL LOCAL OO(:,:,/IQ)
EQUIVALENCE (OP,OO)
C
DIMENSION SSS(3,3),TTT(3,3)
C STRESS CALCULATION FOR KINETIC, HARTREE, EX-COR, LOCAL AND ETOT1
CCC =2.D0*RVOL/(2.D0*PAI)**3/FLOAT(KV3)
DO 1000 IS=1,6
SIGSTR(IS)=0.0D0
1000 CONTINUE
RGG(1)=0.0D0
ZXXX(1)=DCMPLX(0.0D0,0.0D0)
ZYYY(1)=DCMPLX(0.0D0,0.0D0)
ZZZZ(1)=DCMPLX(0.0D0,0.0D0)
DO 1001 I=2,KG
ZXXX(I)=DCMPLX(0.0D0,0.0D0)
ZYYY(I)=DCMPLX(0.0D0,0.0D0)
ZZZZ(I)=DCMPLX(0.0D0,0.0D0)
RGG(I)=1.0D0/(GR(I)*GR(I))
1001 CONTINUE
C KINETIC PART
C VPP-PARALLEL START
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 1100 IK=1,KV3
C IWRT(IK) =IK
AKX = VX(IK)
AKY = VY(IK)
AKZ = VZ(IK)
DO 1110 IBAN=NBD1,NBD2
CW=CCC*OCCUU(IBAN,IK)
DO 1120 I=1,IBA(IK)
I1 = NBASE(I,IK)
STMP=CW*DCONJG(ZZZ(I,IBAN,IK))*ZZZ(I,IBAN,IK)
SIGSTR(1)=SIGSTR(1)-STMP*(AKX+GX(I1))*(AKX+GX(I1))
SIGSTR(2)=SIGSTR(2)-STMP*(AKX+GX(I1))*(AKY+GY(I1))
SIGSTR(3)=SIGSTR(3)-STMP*(AKX+GX(I1))*(AKZ+GZ(I1))
SIGSTR(4)=SIGSTR(4)-STMP*(AKY+GY(I1))*(AKY+GY(I1))
SIGSTR(5)=SIGSTR(5)-STMP*(AKY+GY(I1))*(AKZ+GZ(I1))
SIGSTR(6)=SIGSTR(6)-STMP*(AKZ+GZ(I1))*(AKZ+GZ(I1))
1120 CONTINUE
1110 CONTINUE
1100 CONTINUE
!XOCL END SPREAD SUM(SIGSTR)
C!XOCL END PARALLEL
SSS1=SIGSTR(1)
SSS2=SIGSTR(2)
SSS3=SIGSTR(3)
SSS4=SIGSTR(4)
SSS5=SIGSTR(5)
SSS6=SIGSTR(6)
WRITE (6,*) 'TOTAL STR 1 1 = ',SIGSTR(1)
WRITE (6,*) ' 2 = ',SIGSTR(2)
WRITE (6,*) ' 3 = ',SIGSTR(3)
WRITE (6,*) ' 4 = ',SIGSTR(4)
WRITE (6,*) ' 5 = ',SIGSTR(5)
WRITE (6,*) ' 6 = ',SIGSTR(6)
C HARTREE PART I=1 -----> GR(1)=0.0D0 : RGG(1)=INFINITE
C EXCHANGE-CORRELATION + LOCAL PS PARTS
DO 1250 I=1,KG
ZVP(I)=ZVXC(I)
1250 CONTINUE
CALL XCFFT(1,SCHGPC)
CALL XCFFT(2,SCHGPC)
IF (IPCC.EQ.1) CALL XSTPC(1,SCHGPC)
STMQ=DCONJG(ZEXC(1))*(ZCHG(1)+ZRHPC(1))-DCONJG(ZVXC(1))*ZCHG(1)
& -DCONJG(ZVXC(1))*ZRHPC(1)
SS1=-DCONJG(ZEXC(1))*ZRRPC(1)
SSX=-DCONJG(ZXXX(1))*(ZCHG(1)+ZRHPC(1))
SSY=-DCONJG(ZYYY(1))*(ZCHG(1)+ZRHPC(1))
SSZ=-DCONJG(ZZZZ(1))*(ZCHG(1)+ZRHPC(1))
SIGSTR(1)=SIGSTR(1)+STMQ+SS1*GX(1)*GX(1)+SSX
SIGSTR(4)=SIGSTR(4)+STMQ+SS1*GY(1)*GY(1)+SSY
SIGSTR(6)=SIGSTR(6)+STMQ+SS1*GZ(1)*GZ(I)+SSZ
TTT1=SIGSTR(1)-SSS1
TTT2=SIGSTR(2)-SSS2
TTT3=SIGSTR(3)-SSS3
TTT4=SIGSTR(4)-SSS4
TTT5=SIGSTR(5)-SSS5
TTT6=SIGSTR(6)-SSS6
SSS1=SIGSTR(1)
SSS2=SIGSTR(2)
SSS3=SIGSTR(3)
SSS4=SIGSTR(4)
SSS5=SIGSTR(5)
SSS6=SIGSTR(6)
WRITE (6,*) 'TOTAL STR 2 1 = ',TTT1
WRITE (6,*) ' 2 = ',TTT2
WRITE (6,*) ' 3 = ',TTT3
WRITE (6,*) ' 4 = ',TTT4
WRITE (6,*) ' 5 = ',TTT5
WRITE (6,*) ' 6 = ',TTT6
STX=0.0D0
STY=0.0D0
STZ=0.0D0
STT=0.0D0
STU=0.0D0
STV=0.0D0
STW=0.0D0
DO 1200 I=2,KG
STMP=0.5D0*PAI4*DCONJG(ZCHG(I))*ZCHG(I)*RGG(I)
STMQ=DCONJG(ZEXC(I))*(ZCHG(I)+ZRHPC(I))-DCONJG(ZVXC(I))*ZCHG(I)
& -DCONJG(ZVXC(I))*ZRHPC(I)
STM1= DCONJG(ZEXC(I))*ZRRPC(I)
STMX= DCONJG(ZXXX(I))*(ZCHG(I)+ZRHPC(I))
STMY= DCONJG(ZYYY(I))*(ZCHG(I)+ZRHPC(I))
STMZ= DCONJG(ZZZZ(I))*(ZCHG(I)+ZRHPC(I))
C STX=STX+DREAL(DCONJG(ZVXC(I))*ZRRPC(I))*GX(I)*GX(I)
C STY=STY+DREAL(DCONJG(ZEXC(I))*ZRRPC(I))*GX(I)*GX(I)
C STZ=STZ+STMX
C STT=STT+STMY
STX=STX+STMP*(2.0D0*GX(I)*GX(I)*RGG(I)-1.0D0)+STMQ
& -STMS*2.0D0*GX(I)*GX(I)-STMR-STM1*GX(I)*GX(I)-STMX
STY=STY+STMP*(2.0D0*GX(I)*GX(I)*RGG(I)-1.0D0)
STZ=STZ+STMQ
STT=STT-STMS*2.0D0*GX(I)*GX(I)
STU=STU-STMR
STV=STV-STM1*GX(I)*GX(I)
STW=STW-STMX
STMR=ZPSCC(I)*DCONJG(ZCHG(I))
STMS=ZDSCC(I)*DCONJG(ZCHG(I))
SIGSTR(1)=SIGSTR(1)+STMP*(2.0D0*GX(I)*GX(I)*RGG(I)-1.0D0)+STMQ
& -STMS*2.0D0*GX(I)*GX(I)-STMR-STM1*GX(I)*GX(I)-STMX
SIGSTR(2)=SIGSTR(2)+STMP*(2.0D0*GX(I)*GY(I)*RGG(I))
& -STMS*2.0D0*GX(I)*GY(I)
SIGSTR(3)=SIGSTR(3)+STMP*(2.0D0*GX(I)*GZ(I)*RGG(I))
& -STMS*2.0D0*GX(I)*GZ(I)
SIGSTR(4)=SIGSTR(4)+STMP*(2.0D0*GY(I)*GY(I)*RGG(I)-1.0D0)+STMQ
& -STMS*2.0D0*GY(I)*GY(I)-STMR-STM1*GY(I)*GY(I)-STMY
SIGSTR(5)=SIGSTR(5)+STMP*(2.0D0*GY(I)*GZ(I)*RGG(I))
& -STMS*2.0D0*GY(I)*GZ(I)
SIGSTR(6)=SIGSTR(6)+STMP*(2.0D0*GZ(I)*GZ(I)*RGG(I)-1.0D0)+STMQ
& -STMS*2.0D0*GZ(I)*GZ(I)-STMR-STM1*GZ(I)*GZ(I)-STMZ
1200 CONTINUE
WRITE (6,*) 'ZVX,ZEX = ',STX,STY
WRITE (6,*) 'STM1(XY)= ',STZ,STT
WRITE (6,*) 'STU,V,W = ',STU,STV,STW
STOT=STY+STZ+STT+STU+STV+STW
WRITE (6,*) 'STOT = ',STOT
TTT1=SIGSTR(1)-SSS1
TTT2=SIGSTR(2)-SSS2
TTT3=SIGSTR(3)-SSS3
TTT4=SIGSTR(4)-SSS4
TTT5=SIGSTR(5)-SSS5
TTT6=SIGSTR(6)-SSS6
SSS1=SIGSTR(1)
SSS2=SIGSTR(2)
SSS3=SIGSTR(3)
SSS4=SIGSTR(4)
SSS5=SIGSTR(5)
SSS6=SIGSTR(6)
WRITE (6,*) 'TOTAL STR 3 1 = ',TTT1
WRITE (6,*) ' 2 = ',TTT2
WRITE (6,*) ' 3 = ',TTT3
WRITE (6,*) ' 4 = ',TTT4
WRITE (6,*) ' 5 = ',TTT5
WRITE (6,*) ' 6 = ',TTT6
DO 1251 I=1,KNG
ZVXC(I)=ZVP(I)
1251 CONTINUE
C CALL NON-LOCAL PART SUBROUTINE
CALL STRNL(ETOT1)
C NON-LOCAL PS AND EWALD PARTS
DO 1400 IS=1,6
SIGSTR(IS)=SIGSTR(IS)+SIGNL(IS)+SIGEWA(IS)
1400 CONTINUE
WRITE (6,*) 'SIGEWA 1 = ',SIGEWA(1)
WRITE (6,*) 'SIGEWA 2 = ',SIGEWA(2)
WRITE (6,*) 'SIGEWA 3 = ',SIGEWA(3)
WRITE (6,*) 'SIGEWA 4 = ',SIGEWA(4)
WRITE (6,*) 'SIGEWA 5 = ',SIGEWA(5)
WRITE (6,*) 'SIGEWA 6 = ',SIGEWA(6)
WRITE (6,*) 'SIGNL 1 = ',SIGNL(1)
WRITE (6,*) 'SIGNL 2 = ',SIGNL(2)
WRITE (6,*) 'SIGNL 3 = ',SIGNL(3)
WRITE (6,*) 'SIGNL 4 = ',SIGNL(4)
WRITE (6,*) 'SIGNL 5 = ',SIGNL(5)
WRITE (6,*) 'SIGNL 6 = ',SIGNL(6)
TTT1=SIGSTR(1)-SSS1
TTT2=SIGSTR(2)-SSS2
TTT3=SIGSTR(3)-SSS3
TTT4=SIGSTR(4)-SSS4
TTT5=SIGSTR(5)-SSS5
TTT6=SIGSTR(6)-SSS6
SSS1=SIGSTR(1)
SSS2=SIGSTR(2)
SSS3=SIGSTR(3)
SSS4=SIGSTR(4)
SSS5=SIGSTR(5)
SSS6=SIGSTR(6)
WRITE (6,*) 'TOTAL STR 4 1 = ',TTT1
WRITE (6,*) ' 2 = ',TTT2
WRITE (6,*) ' 3 = ',TTT3
WRITE (6,*) ' 4 = ',TTT4
WRITE (6,*) ' 5 = ',TTT5
WRITE (6,*) ' 6 = ',TTT6
SIGSTR(1)=SIGSTR(1)-ETOT1*TOTCH/UNIVOL
SIGSTR(4)=SIGSTR(4)-ETOT1*TOTCH/UNIVOL
SIGSTR(6)=SIGSTR(6)-ETOT1*TOTCH/UNIVOL
TTT1=SIGSTR(1)-SSS1
TTT2=SIGSTR(2)-SSS2
TTT3=SIGSTR(3)-SSS3
TTT4=SIGSTR(4)-SSS4
TTT5=SIGSTR(5)-SSS5
TTT6=SIGSTR(6)-SSS6
SSS1=SIGSTR(1)
SSS2=SIGSTR(2)
SSS3=SIGSTR(3)
SSS4=SIGSTR(4)
SSS5=SIGSTR(5)
SSS6=SIGSTR(6)
WRITE (6,*) 'TOTAL STR 5 1 = ',TTT1
WRITE (6,*) ' 2 = ',TTT2
WRITE (6,*) ' 3 = ',TTT3
WRITE (6,*) ' 4 = ',TTT4
WRITE (6,*) ' 5 = ',TTT5
WRITE (6,*) ' 6 = ',TTT6
WRITE (6,*) 'TOTAL SUMM 1 = ',SSS1
WRITE (6,*) ' 2 = ',SSS2
WRITE (6,*) ' 3 = ',SSS3
WRITE (6,*) ' 4 = ',SSS4
WRITE (6,*) ' 5 = ',SSS5
WRITE (6,*) ' 6 = ',SSS6
IF (NBZTYP.LE.1) GO TO 9000
DENOM = 1.0D0/FLOAT(KOPR)
SSS(1,1)=SIGSTR(1)
SSS(2,2)=SIGSTR(4)
SSS(3,3)=SIGSTR(6)
SSS(1,2)=SIGSTR(2)
SSS(1,3)=SIGSTR(3)
SSS(2,3)=SIGSTR(5)
SSS(2,1)=SSS(1,2)
SSS(3,1)=SSS(1,3)
SSS(3,2)=SSS(2,3)
DO 1210 I=1,3
DO 1220 J=1,3
TTT(I,J)=0.0D0
1220 CONTINUE
1210 CONTINUE
!XOCL SPREAD DO /IQ
DO 2400 IOP = 1,KOPR
DO 2200 I=1,3
DO 2300 J=1,3
SXX=SSS(1,1)
SXY=SSS(1,2)
SXZ=SSS(1,3)
SYY=SSS(2,2)
SYZ=SSS(2,3)
SZZ=SSS(3,3)
TXX=OO(1,I,IOP)*SSS(I,J)*OO(1,J,IOP)
TXY=OO(1,I,IOP)*SSS(I,J)*OO(2,J,IOP)
TXZ=OO(1,I,IOP)*SSS(I,J)*OO(3,J,IOP)
TYY=OO(2,I,IOP)*SSS(I,J)*OO(2,J,IOP)
TYZ=OO(2,I,IOP)*SSS(I,J)*OO(3,J,IOP)
TZZ=OO(3,I,IOP)*SSS(I,J)*OO(3,J,IOP)
TTT(1,1)=TTT(1,1) + OO(1,I,IOP)*SSS(I,J)*OO(1,J,IOP)
TTT(1,2)=TTT(1,2) + OO(1,I,IOP)*SSS(I,J)*OO(2,J,IOP)
TTT(1,3)=TTT(1,3) + OO(1,I,IOP)*SSS(I,J)*OO(3,J,IOP)
TTT(2,2)=TTT(2,2) + OO(2,I,IOP)*SSS(I,J)*OO(2,J,IOP)
TTT(2,3)=TTT(2,3) + OO(2,I,IOP)*SSS(I,J)*OO(3,J,IOP)
TTT(3,3)=TTT(3,3) + OO(3,I,IOP)*SSS(I,J)*OO(3,J,IOP)
2300 CONTINUE
2200 CONTINUE
2400 CONTINUE
!XOCL END SPREAD SUM(TTT)
SIGSTR(1)=TTT(1,1)*DENOM
SIGSTR(2)=TTT(1,2)*DENOM
SIGSTR(3)=TTT(1,3)*DENOM
SIGSTR(4)=TTT(2,2)*DENOM
SIGSTR(5)=TTT(2,3)*DENOM
SIGSTR(6)=TTT(3,3)*DENOM
SSS1=SIGSTR(1)
SSS2=SIGSTR(2)
SSS3=SIGSTR(3)
SSS4=SIGSTR(4)
SSS5=SIGSTR(5)
SSS6=SIGSTR(6)
WRITE (6,*) 'TOTAL SUMM OP1 = ',SSS1
WRITE (6,*) ' 2 = ',SSS2
WRITE (6,*) ' 3 = ',SSS3
WRITE (6,*) ' 4 = ',SSS4
WRITE (6,*) ' 5 = ',SSS5
WRITE (6,*) ' 6 = ',SSS6
9000 CONTINUE
RETURN
END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE STRNL(ETOT1)
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),SS2(KNG1,KNV3,KTYP,9)
& ,SS3(KNG1,KNV3,KTYP,4),ZZZ(KNG1,KEG,KNV3)
& ,ZFC(KEG,KNV3,KATM,10),ZFC2(KEG,KNV3,KATM,6)
& ,OCCUU(KEG,KNV3),RAA(KNG1,KNV3)
!XOCL LOCAL SSS(:,/IP,:,:),SS2(:,/IP,:,:),SS3(:,/IP,:)
!XOCL LOCAL ZZZ(:,:,/IP),ZFC(:,/IP,:,:),ZFC2(:,/IP,:,:)
!XOCL LOCAL OCCUU(:,/IP),RAA(:,/IP)
EQUIVALENCE (SNL,SSS),(SNL2,SS2),(SNL3,SS3),(ZAJ,ZZZ)
& ,(ZFBB,ZFC),(ZFBB2,ZFC2),(OCCUP,OCCUU)
& ,(RAK,RAA)
C STRESS CALCULATION FOR NON-LOCAL PSEUDOPOTENTIAL PART
KSTR=KEG*KNV3*KATM*6
CCC =2.D0*RVOL/(2.D0*PAI)**3/FLOAT(KV3)
RUNI=2.0D0*CCC
DO 1003 IS=1,6
SIGNL(IS)=0.0D0
1003 CONTINUE
C VPP-PARALLEL START 1
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 1901 IK=1,KV3
DO 1002 IBAN=NBD1,NBD2
DO 1902 IA=1,KATM
ZFC2(IBAN,IK,IA,1)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,2)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,3)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,4)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,5)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,6)=DCMPLX(0.0D0,0.0D0)
1902 CONTINUE
1002 CONTINUE
1901 CONTINUE
!XOCL END SPREAD
C!XOCL END PARALLEL
C
C VPP-PARALLEL 2
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 1010 IK=1,KV3
DO 1000 IA=1,KATM
RWS =WS(KFTYPE(IA))*UNIVOL
RWP =WP(KFTYPE(IA))*UNIVOL
CWL(1)= RWS
CWL(2)= RWP
CWL(3)= RWP
CWL(4)= RWP
IF (NLSPD(KFTYPE(IA)).EQ.1) THEN
LNUM = 4
ELSE
LNUM = 9
RWD =WD(KFTYPE(IA))*UNIVOL
CWL(5)= RWD
CWL(6)= RWD
CWL(7)= RWD
CWL(8)= RWD
CWL(9)= RWD
END IF
C
DO 1022 L=1,LNUM
DO 1020 IBAN=NBD1,NBD2
CW=0.5D0*RUNI*OCCUU(IBAN,IK)
STMP=CW*CWL(L)*DREAL(DCONJG(ZFC(IBAN,IK,IA,L))
& *ZFC(IBAN,IK,IA,L))
SIGNL(1)=SIGNL(1)-STMP
SIGNL(4)=SIGNL(4)-STMP
SIGNL(6)=SIGNL(6)-STMP
1020 CONTINUE
1022 CONTINUE
1000 CONTINUE
1010 CONTINUE
!XOCL END SPREAD SUM(SIGNL)
C!XOCL END PARALLEL
SSS1=SIGNL(1)
SSS2=SIGNL(2)
SSS3=SIGNL(3)
SSS4=SIGNL(4)
SSS5=SIGNL(5)
SSS6=SIGNL(6)
WRITE (6,*) 'SIGNL 1 1 = ',SSS1
WRITE (6,*) ' 2 = ',SSS2
WRITE (6,*) ' 3 = ',SSS3
WRITE (6,*) ' 4 = ',SSS4
WRITE (6,*) ' 5 = ',SSS5
WRITE (6,*) ' 6 = ',SSS6
C
SK1=0.0D0
SK2=0.0D0
SK3=0.0D0
SK4=0.0D0
SK5=0.0D0
SK6=0.0D0
C VPP-PARALLEL START 3
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 1110 IK=1,KV3
AKX = VX(IK)
AKY = VY(IK)
AKZ = VZ(IK)
DO 1100 IA=1,KATM
DO 1120 IBAN=NBD1,NBD2
CW=RUNI*OCCUU(IBAN,IK)
DO 1130 I=1,IBA(IK)
I1 = NBASE(I,IK)
ZTMP=ZZZ(I,IBAN,IK)*DCONJG(ZFM2(I1,IA))
& *SS2(I,IK,KFTYPE(IA),1)*RAA(I,IK)
ZFC2(IBAN,IK,IA,1)=ZFC2(IBAN,IK,IA,1)
& -ZTMP*(AKX+GX(I1))*(AKX+GX(I1))
ZFC2(IBAN,IK,IA,2)=ZFC2(IBAN,IK,IA,2)
& -ZTMP*(AKX+GX(I1))*(AKY+GY(I1))
ZFC2(IBAN,IK,IA,3)=ZFC2(IBAN,IK,IA,3)
& -ZTMP*(AKX+GX(I1))*(AKZ+GZ(I1))
ZFC2(IBAN,IK,IA,4)=ZFC2(IBAN,IK,IA,4)
& -ZTMP*(AKY+GY(I1))*(AKY+GY(I1))
ZFC2(IBAN,IK,IA,5)=ZFC2(IBAN,IK,IA,5)
& -ZTMP*(AKY+GY(I1))*(AKZ+GZ(I1))
ZFC2(IBAN,IK,IA,6)=ZFC2(IBAN,IK,IA,6)
& -ZTMP*(AKZ+GZ(I1))*(AKZ+GZ(I1))
1130 CONTINUE
SIGNL(1)=SIGNL(1)+CW *DREAL(
& ZFC2(IBAN,IK,IA,1)*DCONJG(ZFC(IBAN,IK,IA,1)) )
SIGNL(2)=SIGNL(2)+CW *DREAL(
& ZFC2(IBAN,IK,IA,2)*DCONJG(ZFC(IBAN,IK,IA,1)) )
SIGNL(3)=SIGNL(3)+CW *DREAL(
& ZFC2(IBAN,IK,IA,3)*DCONJG(ZFC(IBAN,IK,IA,1)) )
SIGNL(4)=SIGNL(4)+CW *DREAL(
& ZFC2(IBAN,IK,IA,4)*DCONJG(ZFC(IBAN,IK,IA,1)) )
SIGNL(5)=SIGNL(5)+CW *DREAL(
& ZFC2(IBAN,IK,IA,5)*DCONJG(ZFC(IBAN,IK,IA,1)) )
SIGNL(6)=SIGNL(6)+CW *DREAL(
& ZFC2(IBAN,IK,IA,6)*DCONJG(ZFC(IBAN,IK,IA,1)) )
1120 CONTINUE
1100 CONTINUE
1110 CONTINUE
!XOCL END SPREAD SUM(SIGNL)
C!XOCL END PARALLEL
TTT1=SIGNL(1)-SSS1
TTT2=SIGNL(2)-SSS2
TTT3=SIGNL(3)-SSS3
TTT4=SIGNL(4)-SSS4
TTT5=SIGNL(5)-SSS5
TTT6=SIGNL(6)-SSS6
SSS1=SIGNL(1)
SSS2=SIGNL(2)
SSS3=SIGNL(3)
SSS4=SIGNL(4)
SSS5=SIGNL(5)
SSS6=SIGNL(6)
WRITE (6,*) 'SIGNL 2 1 = ',TTT1
WRITE (6,*) ' 2 = ',TTT2
WRITE (6,*) ' 3 = ',TTT3
WRITE (6,*) ' 4 = ',TTT4
WRITE (6,*) ' 5 = ',TTT5
WRITE (6,*) ' 6 = ',TTT6
C 4
DO 1300 L =2,9
C VPP-PARALLEL START 4
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 1912 IK=1,KV3
DO 1913 IA=1,KATM
IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.4)
& GO TO 1913
DO 1911 IBAN=NBD1,NBD2
ZFC2(IBAN,IK,IA,1)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,2)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,3)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,4)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,5)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,6)=DCMPLX(0.0D0,0.0D0)
1911 CONTINUE
1913 CONTINUE
1912 CONTINUE
!XOCL END SPREAD
C!XOCL END PARALLEL
C
C VPP-PARALLEL START 5
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 1320 IK=1,KV3
AKX = VX(IK)
AKY = VY(IK)
AKZ = VZ(IK)
DO 1310 IA=1,KATM
IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.4)
& GO TO 1310
DO 1330 IBAN=NBD1,NBD2
CW=RUNI*OCCUU(IBAN,IK)
DO 1340 I=1,IBA(IK)
I1 = NBASE(I,IK)
ZTMP=ZZZ(I,IBAN,IK)*DCONJG(ZFM2(I1,IA))
& *( SS2(I,IK,KFTYPE(IA),L) )*RAA(I,IK)
ZFC2(IBAN,IK,IA,1)=ZFC2(IBAN,IK,IA,1)-ZTMP
& *(AKX+GX(I1))*(AKX+GX(I1))
ZFC2(IBAN,IK,IA,2)=ZFC2(IBAN,IK,IA,2)-ZTMP
& *(AKX+GX(I1))*(AKY+GY(I1))
ZFC2(IBAN,IK,IA,3)=ZFC2(IBAN,IK,IA,3)-ZTMP
& *(AKX+GX(I1))*(AKZ+GZ(I1))
ZFC2(IBAN,IK,IA,4)=ZFC2(IBAN,IK,IA,4)-ZTMP
& *(AKY+GY(I1))*(AKY+GY(I1))
ZFC2(IBAN,IK,IA,5)=ZFC2(IBAN,IK,IA,5)-ZTMP
& *(AKY+GY(I1))*(AKZ+GZ(I1))
ZFC2(IBAN,IK,IA,6)=ZFC2(IBAN,IK,IA,6)-ZTMP
& *(AKZ+GZ(I1))*(AKZ+GZ(I1))
1340 CONTINUE
SIGNL(1)=SIGNL(1)+CW *DREAL(
& ZFC2(IBAN,IK,IA,1)*DCONJG( ZFC(IBAN,IK,IA,L)) )
SIGNL(2)=SIGNL(2)+CW *DREAL(
& ZFC2(IBAN,IK,IA,2)*DCONJG( ZFC(IBAN,IK,IA,L)) )
SIGNL(3)=SIGNL(3)+CW *DREAL(
& ZFC2(IBAN,IK,IA,3)*DCONJG( ZFC(IBAN,IK,IA,L)) )
SIGNL(4)=SIGNL(4)+CW *DREAL(
& ZFC2(IBAN,IK,IA,4)*DCONJG( ZFC(IBAN,IK,IA,L)) )
SIGNL(5)=SIGNL(5)+CW *DREAL(
& ZFC2(IBAN,IK,IA,5)*DCONJG( ZFC(IBAN,IK,IA,L)) )
SIGNL(6)=SIGNL(6)+CW *DREAL(
& ZFC2(IBAN,IK,IA,6)*DCONJG( ZFC(IBAN,IK,IA,L)) )
1330 CONTINUE
1310 CONTINUE
1320 CONTINUE
!XOCL END SPREAD SUM(SIGNL)
C!XOCL END PARALLEL
1300 CONTINUE
TTT1=SIGNL(1)-SSS1
TTT2=SIGNL(2)-SSS2
TTT3=SIGNL(3)-SSS3
TTT4=SIGNL(4)-SSS4
TTT5=SIGNL(5)-SSS5
TTT6=SIGNL(6)-SSS6
SSS1=SIGNL(1)
SSS2=SIGNL(2)
SSS3=SIGNL(3)
SSS4=SIGNL(4)
SSS5=SIGNL(5)
SSS6=SIGNL(6)
WRITE (6,*) 'SIGNL 3 1 = ',TTT1
WRITE (6,*) ' 2 = ',TTT2
WRITE (6,*) ' 3 = ',TTT3
WRITE (6,*) ' 4 = ',TTT4
WRITE (6,*) ' 5 = ',TTT5
WRITE (6,*) ' 6 = ',TTT6
CWL(1)=1.0D0
CWL(2)=1.0D0
CWL(3)=1.0D0
CWL(4)=1.0D0
CWL(5)=2.0D0
CWL(6)=2.0D0
CWL(7)=2.0D0
CWL(8)=2.0D0
CWL(9)=2.0D0
CWL(10)=1.0D0
C 4
DO 1400 L=2,10
C VPP-PARALLEL START 6
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 1922 IK=1,KV3
DO 1923 IA=1,KATM
IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.4)
& GO TO 1923
DO 1921 IBAN=NBD1,NBD2
ZFC2(IBAN,IK,IA,1)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,2)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,3)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,4)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,5)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,6)=DCMPLX(0.0D0,0.0D0)
1921 CONTINUE
1923 CONTINUE
1922 CONTINUE
!XOCL END SPREAD
C!XOCL END PARALLEL
C
C VPP-PARALLEL START 7
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 1420 IK=1,KV3
AKX = VX(IK)
AKY = VY(IK)
AKZ = VZ(IK)
DO 1410 IA=1,KATM
IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.4)
& GO TO 1410
DO 1430 IBAN=NBD1,NBD2
CW=CWL(L)*RUNI*OCCUU(IBAN,IK)
DO 1440 I=1,IBA(IK)
I1 = NBASE(I,IK)
ZTMP=ZZZ(I,IBAN,IK)*DCONJG(ZFM2(I1,IA))
& *SSS(I,IK,KFTYPE(IA),L)*RAA(I,IK)*RAA(I,IK)
ZFC2(IBAN,IK,IA,1)=ZFC2(IBAN,IK,IA,1)+ZTMP
& *(AKX+GX(I1))*(AKX+GX(I1))
ZFC2(IBAN,IK,IA,2)=ZFC2(IBAN,IK,IA,2)+ZTMP
& *(AKX+GX(I1))*(AKY+GY(I1))
ZFC2(IBAN,IK,IA,3)=ZFC2(IBAN,IK,IA,3)+ZTMP
& *(AKX+GX(I1))*(AKZ+GZ(I1))
ZFC2(IBAN,IK,IA,4)=ZFC2(IBAN,IK,IA,4)+ZTMP
& *(AKY+GY(I1))*(AKY+GY(I1))
ZFC2(IBAN,IK,IA,5)=ZFC2(IBAN,IK,IA,5)+ZTMP
& *(AKY+GY(I1))*(AKZ+GZ(I1))
ZFC2(IBAN,IK,IA,6)=ZFC2(IBAN,IK,IA,6)+ZTMP
& *(AKZ+GZ(I1))*(AKZ+GZ(I1))
1440 CONTINUE
SIGNL(1)=SIGNL(1)+CW *DREAL(
& ZFC2(IBAN,IK,IA,1)*DCONJG(ZFC(IBAN,IK,IA,L)) )
SIGNL(2)=SIGNL(2)+CW *DREAL(
& ZFC2(IBAN,IK,IA,2)*DCONJG(ZFC(IBAN,IK,IA,L)) )
SIGNL(3)=SIGNL(3)+CW *DREAL(
& ZFC2(IBAN,IK,IA,3)*DCONJG(ZFC(IBAN,IK,IA,L)) )
SIGNL(4)=SIGNL(4)+CW *DREAL(
& ZFC2(IBAN,IK,IA,4)*DCONJG(ZFC(IBAN,IK,IA,L)) )
SIGNL(5)=SIGNL(5)+CW *DREAL(
& ZFC2(IBAN,IK,IA,5)*DCONJG(ZFC(IBAN,IK,IA,L)) )
SIGNL(6)=SIGNL(6)+CW *DREAL(
& ZFC2(IBAN,IK,IA,6)*DCONJG(ZFC(IBAN,IK,IA,L)) )
1430 CONTINUE
1410 CONTINUE
1420 CONTINUE
!XOCL END SPREAD SUM(SIGNL)
C!XOCL END PARALLEL
1400 CONTINUE
TTT1=SIGNL(1)-SSS1
TTT2=SIGNL(2)-SSS2
TTT3=SIGNL(3)-SSS3
TTT4=SIGNL(4)-SSS4
TTT5=SIGNL(5)-SSS5
TTT6=SIGNL(6)-SSS6
SSS1=SIGNL(1)
SSS2=SIGNL(2)
SSS3=SIGNL(3)
SSS4=SIGNL(4)
SSS5=SIGNL(5)
SSS6=SIGNL(6)
WRITE (6,*) 'SIGNL 4 1 = ',TTT1
WRITE (6,*) ' 2 = ',TTT2
WRITE (6,*) ' 3 = ',TTT3
WRITE (6,*) ' 4 = ',TTT4
WRITE (6,*) ' 5 = ',TTT5
WRITE (6,*) ' 6 = ',TTT6
C 1
DO 2331 L=1,4
C VPP-PARALLEL START 8
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 1942 IK=1,KV3
DO 1943 IA=1,KATM
IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.1)
& GO TO 1943
DO 1941 IBAN=NBD1,NBD2
ZFC2(IBAN,IK,IA,1)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,2)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,3)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,4)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,5)=DCMPLX(0.0D0,0.0D0)
ZFC2(IBAN,IK,IA,6)=DCMPLX(0.0D0,0.0D0)
1941 CONTINUE
1943 CONTINUE
1942 CONTINUE
!XOCL END SPREAD
C!XOCL END PARALLEL
C
C VPP-PARALLEL START 9
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
DO 1520 IK=1,KV3
C IWRT(IK) =IK
AKX = VX(IK)
AKY = VY(IK)
AKZ = VZ(IK)
DO 1510 IA=1,KATM
IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.1)
& GO TO 1510
IF (L.EQ.1) THEN
RWC=1.0D0/(WP(KFTYPE(IA))*UNIVOL)
ELSE
RWC=1.0D0/(WD(KFTYPE(IA))*UNIVOL)
END IF
DO 1530 IBAN=NBD1,NBD2
CW=RUNI*OCCUU(IBAN,IK)
DO 1540 I=1,IBA(IK)
I1 = NBASE(I,IK)
ZTMP=DCONJG(ZZZ(I,IBAN,IK))*ZFM2(I1,IA)
& *SS3(I,IK,KFTYPE(IA),L)*RAA(I,IK)
ZFC2(IBAN,IK,IA,1)=ZFC2(IBAN,IK,IA,1)-ZTMP
& *(AKX+GX(I1))
ZFC2(IBAN,IK,IA,2)=ZFC2(IBAN,IK,IA,2)-ZTMP
& *(AKY+GY(I1))
ZFC2(IBAN,IK,IA,3)=ZFC2(IBAN,IK,IA,3)-ZTMP
& *(AKZ+GZ(I1))
1540 CONTINUE
SIGNL(1)=SIGNL(1)-RWC*CW *DREAL(
& ZFC2(IBAN,IK,IA,1)*DCONJG(ZFC2(IBAN,IK,IA,1)) )
SIGNL(2)=SIGNL(2)-RWC*CW *DREAL(
& ZFC2(IBAN,IK,IA,1)*DCONJG(ZFC2(IBAN,IK,IA,2)) )
SIGNL(3)=SIGNL(3)-RWC*CW *DREAL(
& ZFC2(IBAN,IK,IA,1)*DCONJG(ZFC2(IBAN,IK,IA,3)) )
SIGNL(4)=SIGNL(4)-RWC*CW *DREAL(
& ZFC2(IBAN,IK,IA,2)*DCONJG(ZFC2(IBAN,IK,IA,2)) )
SIGNL(5)=SIGNL(5)-RWC*CW *DREAL(
& ZFC2(IBAN,IK,IA,2)*DCONJG(ZFC2(IBAN,IK,IA,3)) )
SIGNL(6)=SIGNL(6)-RWC*CW *DREAL(
& ZFC2(IBAN,IK,IA,3)*DCONJG(ZFC2(IBAN,IK,IA,3)) )
1530 CONTINUE
1510 CONTINUE
1520 CONTINUE
!XOCL END SPREAD SUM(SIGNL)
C!XOCL END PARALLEL
2331 CONTINUE
TTT1=SIGNL(1)-SSS1
TTT2=SIGNL(2)-SSS2
TTT3=SIGNL(3)-SSS3
TTT4=SIGNL(4)-SSS4
TTT5=SIGNL(5)-SSS5
TTT6=SIGNL(6)-SSS6
SSS1=SIGNL(1)
SSS2=SIGNL(2)
SSS3=SIGNL(3)
SSS4=SIGNL(4)
SSS5=SIGNL(5)
SSS6=SIGNL(6)
WRITE (6,*) 'SIGNL 5 1 = ',TTT1
WRITE (6,*) ' 2 = ',TTT2
WRITE (6,*) ' 3 = ',TTT3
WRITE (6,*) ' 4 = ',TTT4
WRITE (6,*) ' 5 = ',TTT5
WRITE (6,*) ' 6 = ',TTT6
RETURN
END
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE KSTEP
> (KNV3,NBZTYP,ALTV,RLTV,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,NFOUT,
< KV3,VX,VY,VZ,QWGT)
IMPLICIT LOGICAL(A-Z)
INTEGER KNV3,NBZTYP,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,KV3
INTEGER INDX(0:20,0:20,0:20),IK,NFOUT
REAL ALTV(3,3),RLTV(3,3),A,B,C,PAI2,A1,B1,C1
REAL VX(KNV3),VY(KNV3),VZ(KNV3),QWGT(KNV3)
LOGICAL SPINP
PAI2 = 8.D0*ATAN(1.D0)
SPINP = .FALSE.
C*** IF((NBZTYP.EQ.0).OR.(NBZTYP.EQ.1)) THEN
IF((NBZTYP.EQ.0).OR.(NBZTYP.EQ.1).OR.(NBZTYP.EQ.9)
> .OR.(NBZTYP.GE.11)) THEN
IF(NBZTYP.EQ.0) THEN
C----------------------------------------------
WRITE(NFOUT,*) ' SURFACE B.Z. '
CALL KPMSF
> (KNV3,RLTV,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,NFOUT,
< KV3,VX,VY,VZ)
C----------------------------------------------
C**** ELSE IF(NBZTYP.EQ.1) THEN
ELSE IF(NBZTYP.EQ.1.OR.NBZTYP.EQ.9.OR.NBZTYP.GE.11) THEN
C----------------------------------------------
WRITE(NFOUT,*) ' WHOLE B.Z. '
CALL KPMWBZ
> (KNV3,RLTV,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,NFOUT,
< KV3,VX,VY,VZ)
C----------------------------------------------
END IF
DO 100 IK=1,KNV3
QWGT(IK) = 1.D0/DFLOAT(KV3)
100 CONTINUE
ELSE IF((NBZTYP.GE.2).AND.(NBZTYP.LE.5)) THEN
IF(NBZTYP.EQ.2) THEN
C----------------------------------------------------------
WRITE(NFOUT,*) ' SIMPLE CUBIC LATTICE '
CALL SCCM
> (KNV3,NKX,SPINP,NFOUT,
< KV3,VX,VY,VZ,QWGT,INDX)
C----------------------------------------------------------
ELSE IF(NBZTYP.EQ.3) THEN
C----------------------------------------------------------
WRITE(NFOUT,*) ' BCC LATTICE '
CALL BCCM
> (KNV3,NKX,SPINP,NFOUT,
< KV3,VX,VY,VZ,QWGT,INDX)
C----------------------------------------------------------
ELSE IF((NBZTYP.EQ.4).OR.(NBZTYP.EQ.5)) THEN
C----------------------------------------------------------
WRITE(NFOUT,*) ' FCC LATTICE '
CALL FCCM
> (KNV3,NKX,SPINP,NFOUT,
< KV3,VX,VY,VZ,QWGT,INDX)
C----------------------------------------------------------
END IF
A = ABS(ALTV(1,1)) * 2.D0
B = PAI2/A
DO 200 IK = 1,KV3
VX(IK) = B*VX(IK)
VY(IK) = B*VY(IK)
VZ(IK) = B*VZ(IK)
200 CONTINUE
ELSE IF(NBZTYP.EQ.6) THEN
C----------------------------------------------------------
WRITE(NFOUT,*) ' HEX LATTICE '
CALL HEXM
> (KNV3,NKX,NKY,RLTV,SPINP,NFOUT,
< KV3,VX,VY,VZ,QWGT,INDX)
C----------------------------------------------------------
ELSE IF (NBZTYP.EQ.8 .OR. NBZTYP.EQ.10) THEN
WRITE (NFOUT,*) ' TETRAGONAL LATTICE '
CALL TETRAH
> (KNV3,NKX,NKZ,RLTV,SPINP,NFOUT,
< KV3,VX,VY,VZ,QWGT,INDX)
A = ABS(ALTV(1,1))*2.0D0
A1 = PAI2/A
B = ABS(ALTV(3,3))*2.0D0
B1 = PAI2/B
WRITE (6,*) A,B,A1,B1,ALTV(1,1),ALTV(3,3)
DO 201 IK = 1,KV3
VX(IK) = A1*VX(IK)
VY(IK) = A1*VY(IK)
VZ(IK) = B1*VZ(IK)
201 CONTINUE
ELSE IF (NBZTYP.EQ.9.OR.NBZTYP.GE.11) THEN
WRITE (NFOUT,*) ' ORTHORHONBIC LATTICE '
CALL APBO2
> (KNV3,NKX,NKY,NKZ,RLTV,SPINP,NFOUT,
< KV3,VX,VY,VZ,QWGT,INDX)
A = ABS(ALTV(1,1))*2.0D0
A1 = PAI2/A
B = ABS(ALTV(2,2))*2.0D0
B1 = PAI2/B
C = ABS(ALTV(3,3))*2.0D0
C1 = PAI2/C
DO 202 IK = 1,KV3
VX(IK) = A1*VX(IK)
VY(IK) = B1*VY(IK)
VZ(IK) = C1*VZ(IK)
202 CONTINUE
ELSE
WRITE(NFOUT,*) ' NBZTYP ERR IN KSTEP NBUB = ',NBZTYP
STOP
END IF
C KV3=1
C A = ABS(ALTV(1,1))
C A1 = PAI2/A
C B = ABS(ALTV(3,3))
C B1 = PAI2/B
C WRITE (6,*) A,B,A1,B1,ALTV(1,1),ALTV(3,3)
C DO 202 IK = 1,KV3
C VX(IK) = A1*VX(IK)
C VY(IK) = A1*VY(IK)
C VZ(IK) = B1*VZ(IK)
C 202 CONTINUE
WRITE(NFOUT,300) (IK,VX(IK),VY(IK),VZ(IK),QWGT(IK),IK=1,KV3)
300 FORMAT((' ','IK = ',I3,3F10.6,3X,'QWGT = ',F10.6))
RETURN
END
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE KPMWBZ
> (KNV3,RLTV,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,NFOUT,
< KV3,VX,VY,VZ)
C************************************************************
IMPLICIT LOGICAL(A-Z)
INTEGER IKP,KNV3,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,KV3,MM,I1,I2,I3,NFOUT
REAL RLTV(3,3),VX(KNV3),VY(KNV3),VZ(KNV3),FVX,FVY,FVZ
REAL DGVX,DGVY,DGVZ,FV1,FV2,FV3
FVX=DFLOAT(NKX*2)
FVY=DFLOAT(NKY*2)
FVZ=DFLOAT(NKZ*2)
IKP=0
IF (IKP.EQ.0) THEN
MM=0
C DGVX=0.5D0*(RLTV(1,1)/FVX+RLTV(1,2)/FVY+RLTV(1,3)/FVZ)
C DGVY=0.5D0*(RLTV(2,1)/FVX+RLTV(2,2)/FVY+RLTV(2,3)/FVZ)
C DGVZ=0.5D0*(RLTV(3,1)/FVX+RLTV(3,2)/FVY+RLTV(3,3)/FVZ)
DGVX=0.0D0
DGVY=0.0D0
DGVZ=0.0D0
DO 100 I1=1,NKX2
FV1=DFLOAT(I1-NKX-1)/FVX
DO 100 I2=1,NKY2
FV2=DFLOAT(I2-NKY-1)/FVY
DO 100 I3=1,NKZ2
FV3=DFLOAT(I3-NKZ-1)/FVZ
MM=MM+1
VX(MM)=RLTV(1,1)*FV1 + RLTV(1,2)*FV2 + RLTV(1,3)*FV3 +DGVX
VY(MM)=RLTV(2,1)*FV1 + RLTV(2,2)*FV2 + RLTV(2,3)*FV3 +DGVY
VZ(MM)=RLTV(3,1)*FV1 + RLTV(3,2)*FV2 + RLTV(3,3)*FV3 +DGVZ
100 CONTINUE
WRITE(NFOUT,*) 'NUMBER OF GENERATED K-POINTS=',MM
IF(MM.NE.KNV3) WRITE(NFOUT,*) '**WARNING MM SHOULD BE KNV3**'
KV3=MM
ELSE IF (IKP.EQ.1) THEN
VX(1)=0.0D0
VY(1)=0.0D0
VZ(1)=0.0D0
KV3=1
END IF
RETURN
END
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE KPMSF
> (KNV3,RLTV,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,NFOUT,
< KV3,VX,VY,VZ)
IMPLICIT LOGICAL(A-Z)
INTEGER KNV3,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,KV3,MM,I2,I3,NFOUT
REAL RLTV(3,3),VX(KNV3),VY(KNV3),VZ(KNV3),FVY,FVZ
REAL DGVY,DGVZ,FV2,FV3
FVY=DFLOAT(NKY*2)
FVZ=DFLOAT(NKZ*2)
MM=0
DGVY=0.5D0*(RLTV(2,2)/FVY+RLTV(2,3)/FVZ)
DGVZ=0.5D0*(RLTV(3,2)/FVY+RLTV(3,3)/FVZ)
DO 100 I2=1,NKY2
FV2=DFLOAT(I2-NKY-1)/FVY
DO 100 I3=1,NKZ2
FV3=DFLOAT(I3-NKZ-1)/FVZ
MM=MM+1
VX(MM) = 0.D0
VY(MM) = RLTV(2,2)*FV2 + RLTV(2,3)*FV3 + DGVY
VZ(MM) = RLTV(3,2)*FV2 + RLTV(3,3)*FV3 + DGVZ
100 CONTINUE
WRITE(NFOUT,*) 'NUMBER OF GENERATED K-POINTS=',MM
IF(MM.NE.KNV3) WRITE(NFOUT,*) '**WARN MM SHOULD BE KNV3**'
KV3=MM
RETURN
END