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