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