C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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