C^^^^^^^^^^^^^^^^LOCAL ZZZ,EKK^^^^^^^^^^^^^^^^^^^^^^^ 5/1, 2009
SUBROUTINE MSD(IREC8,IREC9,IMD,ICAR,ISTR,ZVCOU)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
DIMENSION ZDEVC(KNG1,KEG),EMAS(KEG),ZVCOU(KNG)
C-----ARRAYS FOR MFFT--------------------------------------------
DIMENSION ZEV3C(IFX2-1,IFY2,IFZ2),ZV3D(IFX2-1,IFY2,IFZ2)
DIMENSION XL(KNG1),YL1(KNG1),YL2(KNG1),YL3(KNG1)
C DIMENSION ZV11D(KFFT),ZC11D(KFFT)
DIMENSION IWL(8*IFX2+28),IWM(8*IFY2+28),IWN(8*IFZ2+28)
& ,IWORK(2*IFX2)
DIMENSION ZZZ(KNG1,KEG,KNV3),EKK(KEG,KNV3)
C EQUIVALENCE (ZAJ,ZZZ),(EKO,EKK)
C & ,(SNL,SSS),(ZFBB,ZFC),(OCCUP,OCCUU),(RAK,RAA)
C DIMENSION SSS(KNG1,KNV3,KTYP,10),OCCUU(KEG,KNV3)
CCCCC!$OMP THREADPRIVATE( /WORK8/ ,/RECIP/,/NONLC/, /XY123/ )
CCCCC!$OMP THREADPRIVATE( /ZAJEKO/, /WORK4/, /PSSNL/, /WORK8/,
CCCCC!$OMP& /RECIP/, /NONLC/ )
C!$OMP THREADPRIVATE( /ZAJEKO/, /WORK4/, /PSSNL/ )
C================================================================
C EQUIVALENCE (ZV1D(1),ZV3D(1,1,1)),(ZC1D(1),ZEV3C(1,1,1))
C EQUIVALENCE (ZC11D(1),ZEV3C(1,1,1))
C==== ATTN KX1 OR KX11 ==========================================
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
ZVCOU(I)= PAI4*ZCHGO(I)*RGG(I)
ZVP(I) = ZVXC(I) + ZPSCC(I) + ZVCOU(I)
2001 CONTINUE
C*****---- SET UP C3FFT (FAST FOURIER TRANSFORMATION) -----
C CALL C3FFT(ZEV3C,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
C & ,0,0,1,IWORK,IERR)
C CALL C3FFT_DXML(ZEV3C,KFT1,KFT1-1,KFT2,KFT3,0)
C DO 234 I=1,KSUM
C ZV1D(I)=DCMPLX(0.0D0,0.0D0)
C 234 CONTINUE
DO 7731 KK=1,KFT3
DO 7732 JJ=1,KFT2
DO 7733 II=1,KFT1-1
ZV3D(II,JJ,KK)=DCMPLX(0.0D0,0.0D0)
7733 CONTINUE
7732 CONTINUE
7731 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_MKL2(ZV3D,KFT1,KFT1-1,KFT2,KFT3,1,IERR)
C CALL C3FFT(ZV3D ,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
C & ,0,1,1,IWORK,IERR2)
C CALL C3FFT_DXML(ZV3D,KFT1,KFT1-1,KFT2,KFT3,1)
C
C DO 3501 NNN = 1,KV3
C DO 3502 IBAN = NBD1,NBD2
C EKK(IBAN,NNN) = EKO(IBAN,NNN)
C DO 3503 I= 1,IBA(NNN)
C ZZZ(I,IBAN,NNN) = ZAJ(I,IBAN,NNN)
C 3503 CONTINUE
C 3502 CONTINUE
C 3501 CONTINUE
C
C-----K-POINT LOOP---------------------------------------------
C
C OEPNMP 4/9, 2009 FOR SGI-SUPERCOMPUTER
C
C X,Y1,Y2,Y3, (MSD)
!$OMP PARALLEL DEFAULT(NONE)
!$OMP& FIRSTPRIVATE( KFT1,KFT2,KFT3,KSUM,KVOL, ZV3D, DTIM )
C!$OMP& COPYIN( ZAJ, EKO, ZFBB, SNL, CWL )
CC!$OMP& /ZAJEKO/, /RECIP/, /NONLC/, /WORK8/ )
!$OMP& PRIVATE(NNN,IBAN,VKIN,AKX,AKY,AKZ,I,J,K,L,M,
!$OMP& IIBA,ITY,CD,CS,CP,LNUM,TMPP,ZEV3C,
!$OMP& I1,L1,L2,L3,ierr1,IA,ZTMP,
!$OMP& ierr,ZPGG,WDI,ZROF,ZDEVC,EMAS,ESHI,ZNRM,
!$OMP& XL,YL1,YL2,YL3,
!$OMP& DENOM,JBAN,EE,ZTV,IV )
!$OMP& SHARED( KV3,NBD1,NBD2,
!$OMP& VX,VY,VZ,ZI,IBA,IGF1,IGF2,IGF3,
!$OMP& UNIVOL,WS,WP,WD,IATOM,
!$OMP& KFTYPE, ZFM2,NBASE, ZZZ,EKK, ZAJ,EKO,ZFBB,SNL,CWL,
!$OMP& GX,GY,GZ,NLSPD,
!$OMP& ZVP )
!$OMP DO
DO 100 NNN=1,KV3
AKX = VX(NNN)
AKY = VY(NNN)
AKZ = VZ(NNN)
IIBA = IBA(NNN)
C WRITE (6,*) 'CHECK POINT 1 IN MSD -- OPENMP, NNN = ',NNN
DO 3 IBAN=NBD1,NBD2
EMAS(IBAN)=0.0D0
EKK(IBAN,NNN) = EKO(IBAN,NNN)
DO 3503 I= 1,KNG1
ZZZ(I,IBAN,NNN) = ZAJ(I,IBAN,NNN)
3503 CONTINUE
3 CONTINUE
DO 2011 J=1,KNG1
XL(J)=0.0D0
YL1(J)=0.0D0
YL2(J)=0.0D0
YL3(J)=0.0D0
2011 CONTINUE
C
C
C WRITE (6,*) 'CHECK POINT 2 IN MSD -- OPENMP'
DO 2002 J=1,IIBA
C WRITE (6,*) 'CHECK POINT 22 IN MSD -- OPENMP',J,IIBA
XL(J) = 0.0D0
VKIN =
& ( (AKX+GX(NBASE(J,NNN)))**2
& + (AKY+GY(NBASE(J,NNN)))**2
& + (AKZ+GZ(NBASE(J,NNN)))**2 )/2.D0
YL3(J) = VKIN + DREAL(ZVP(1))
C & Y3(J),ZVP(1)
2002 CONTINUE
C WRITE (6,*) 'CHECK POINT 3 IN MSD -- OPENMP'
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)*SNL(I,NNN,ITY,L)**2
YL3(I) = YL3(I) + TMPP*DFLOAT(IATOM(ITY))
XL(I) = XL(I) + TMPP*DFLOAT(IATOM(ITY))
7010 CONTINUE
7111 CONTINUE
7000 CONTINUE
C WRITE (6,*) 'CHECK POINT 4 IN MSD -- OPENMP'
C-----BAND LOOP----------------------------------------------
DO 200 IBAN=NBD1,NBD2
C WRITE (6,*) 'CHECK POINT 5 IN MSD -- OPENMP, IBAND = ',IBAN
C DO 3231 I=1,KSUM
C ZC11D(I)=CMPLX(0.0D0,0.0D0)
C 3231 CONTINUE
DO 7131 KK=1,KFT3
DO 7132 JJ=1,KFT2
DO 7133 II=1,KFT1-1
ZEV3C(II,JJ,KK)=DCMPLX(0.0D0,0.0D0)
7133 CONTINUE
7132 CONTINUE
7131 CONTINUE
C WRITE (6,*) 'CHECK POINT 6 IN MSD -- OPENMP, IBAND = ',IBAN
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)
3232 CONTINUE
C WRITE (6,*) 'CHECK POINT 7 IN MSD -- OPENMP, IBAND = ',IBAN
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
C
C CALL C3FFT_MKL(ZRB,KFT1,KFT1-1,KFT2,KFT3,IFFT,IERR)
C
CALL C3FFT_MKL2(ZEV3C,KFT1,KFT1-1,KFT2,KFT3,1,IERR1)
C WRITE (6,*) 'CHECK POINT 9 IN MSD -- OPENMP, IBAND = ',IBAN
C DO 71 I=1,KNG1
DO 71 I=1,IIBA
YL1(I) = - DREAL(ZZZ(I,IBAN,NNN)*XL(I))
YL2(I) = - DIMAG(ZZZ(I,IBAN,NNN)*XL(I))
71 CONTINUE
C WRITE (6,*) 'CHECK POINT A IN MSD -- OPENMP, IBAND = ',IBAN
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 = SNL(I,NNN,KFTYPE(IA),L)
& * ZFBB(IBAN,NNN,IA,L)
YL1(I) = YL1(I) + DREAL(ZTMP*ZFM2(NBASE(I,NNN),IA))
YL2(I) = YL2(I) + DIMAG(ZTMP*ZFM2(NBASE(I,NNN),IA))
400 CONTINUE
411 CONTINUE
401 CONTINUE
C DO 240 I=1,KSUM
C ZC11D(I)=ZV11D(I)*ZC11D(I)
C 240 CONTINUE
DO 7231 KK=1,KFT3
DO 7232 JJ=1,KFT2
DO 7233 II=1,KFT1-1
ZEV3C(II,JJ,KK)=ZV3D(II,JJ,KK)*ZEV3C(II,JJ,KK)
7233 CONTINUE
7232 CONTINUE
7231 CONTINUE
C*****---- FAST FOURIER TRANSFORMATION -----
C
C CALL C3FFT_MKL(ZRB,KFT1,KFT1-1,KFT2,KFT3,IFFT,IERR)
C
CALL C3FFT_MKL2(ZEV3C,KFT1,KFT1-1,KFT2,KFT3,-1,IERR)
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 =(YL3(I)-EKK(IBAN,NNN))
ZROF =(ZPGG+YL1(I)+ZI*YL2(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 WRITE (6,*) 'CHECK POINT B IN MSD -- OPENMP END 200-LOOP'
C >>>>> AFTER C(T+DT)=C(T)+DC(T) <<<<<
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C GRAM-SCHMIDT ORTHOGONALIZATION
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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******************************************************************
C WRITE (6,*) 'CHECK POINT C IN MSD -- OPENMP END GRAM=SCHMIDT'
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)
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 WRITE (6,*) 'CHECK POINT D IN MSD -- OPENMP BOTTOM OF 100-LOOP'
100 CONTINUE
!$OMP END DO
!$OMP END PARALLEL
C
DO 3511 NNN = 1,KV3
DO 3512 IBAN = NBD1,NBD2
EKO(IBAN,NNN) = EKK(IBAN,NNN)
DO 3513 I= 1,KNG1
ZAJ(I,IBAN,NNN) = ZZZ(I,IBAN,NNN)
3513 CONTINUE
3512 CONTINUE
3511 CONTINUE
C
RETURN
END
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
C Program COMPLEX_3D_DOUBLE_EX1
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE C3FFT_MKL2(ZRB,KFT1,KFT1_1,KFT2,KFT3,IFFT,IERR)
Use MKL_DFTI
INCLUDE 'PACVPP'
integer II, JJ, KK
integer IFFT, IERR, KFT1,KFT2,KFT3,KFT1_1, INDX
complex(8) :: X_IN ((IFX2-1)*IFY2*IFZ2)
complex(8) :: ZRB(IFX2-1,IFY2,IFZ2)
type(DFTI_DESCRIPTOR), POINTER :: Desc_Handle
integer status
real(8) Scale
integer lengths(3)
IERR = 0
DO 111 KK=1,KFT3
DO 112 JJ=1,KFT2
DO 113 II=1,KFT1-1
INDX=(KFT1-1)*KFT2*(KK-1)+(KFT1-1)*(JJ-1)+II
X_IN(INDX)=ZRB(II,JJ,KK)
113 CONTINUE
112 CONTINUE
111 CONTINUE
!
! Put transform parameters
!
lengths(1) = KFT1_1
lengths(2) = KFT2
lengths(3) = KFT3
!
! Create Dfti descriptor
!
Status = DftiCreateDescriptor( Desc_Handle,
& DFTI_DOUBLE, DFTI_COMPLEX, 3, lengths)
C
IF (IFFT .EQ. -1) THEN
C
Status = DftiCommitDescriptor( Desc_Handle )
!
! Compute Forward transform
!
Status = DftiComputeForward( Desc_Handle, X_IN)
C
ELSE IF (IFFT .EQ. 1) THEN
!
! Set Scale number for Backward transform
!
Scale = 1.0
Status = DftiSetValue(Desc_Handle, DFTI_BACKWARD_SCALE, Scale)
!
! Commit Dfti descriptor
!
Status = DftiCommitDescriptor( Desc_Handle )
!
! Compute Backward transform
!
Status = DftiComputeBackward( Desc_Handle, X_IN)
C
END IF
C
DO 211 KK=1,KFT3
DO 212 JJ=1,KFT2
DO 213 II=1,KFT1-1
INDX=(KFT1-1)*KFT2*(KK-1)+(KFT1-1)*(JJ-1)+II
ZRB(II,JJ,KK)=X_IN(INDX)
213 CONTINUE
212 CONTINUE
211 CONTINUE
C
Status = DftiFreeDescriptor(Desc_Handle)
C
end