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
cGEN
*pdir reserve = 2
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')
cGEN
*pdir release
STOP
END