您所在位置：主题：求大神知道Fortran编程解决单晶铜分子动力学模拟问题急死了
程序跟子程序如下，建立的平面模型，解决单晶铜梁分子弯曲的分子动力学模拟。提取应力应变问题上不知道如何下手了，还有温度控制是否合理，有大神能帮我看一下这个程序么，研究生毕业成问题了，急急急急急急急急急急急急！！！！！
主程序：
DIMENSION VX(3810,2),VY(3810,2),NFORM(3810,40)
DIMENSION LENG(3810,2),F(3810,3),FF(3810,3)
DIMENSION RIJ(3810,40),ROI(3810,2),XYL(3810,3)
DOUBLE PRECISION XYL,VX,VY,VM,T0,FF,PE,RIJ,ROI,F,H,E,EE
DOUBLE PRECISION ZPVX2,ZPVY2,BETAX,BETAY
DOUBLE PRECISION kb,xis1,xis2
OPEN(3,FILE='VX.DAT',STATUS='old')
OPEN(4,FILE='VY.DAT',STATUS='old')
OPEN(10,FILE='XYL.DAT',STATUS='old')
OPEN(11,FILE='RESULT.DAT',STATUS='old')
OPEN(12,FILE='AA.DAT',STATUS='old')
OPEN(15,FILE='F.DAT',STATUS='old')
OPEN(16,FILE='H.DAT',STATUS='old')
OPEN(17,FILE='E.DAT',STATUS='old')
npart=3810
cum=1.06E-25
kb=1.380650524E-23
ev=1.602E-19
T=293
D=0.000000000255
xis1=3.917489092446246E-6
xis2=1.847872213418041E-5
CALL CRADOM(VX,VY,npart)
WRITE(3,102)((VX(I,J),J=1,2),I=1,npart)
WRITE(4,102)((VY(I,J),J=1,2),I=1,npart)
102 FORMAT(1X,f8.1,F14.6)
CALL BCONSTRUCT(XYL)
WRITE(12,*)((XYL(I,J),J=1,3),I=1,npart)
CALL DLISTFORM(XYL,NFORM,LENG,npart)
VM=0.0
DO 702 I=1,npart
VM=VM+VX(I,2)**2+VY(I,2)**2
702 CONTINUE
EK=0.5*cum*VM !新动能J
T0=EK/(npart-1)/kb
EK=EK/ev
CALL EPENERGY(LENG,NFORM,XYL,PE,ROI,RIJ,npart)
WRITE(*,307)PE,EK,T0
307FORMAT(1X,'NUM=0',6X,'PE=',F18.4,2X,'EK=',F10.4,2X,'T0=',F10.4)
CALL FORCE(LENG,XYL,RIJ,ROI,NFORM,F,npart)
NUMB=1 !更新邻域列表及循环控制用
1000 DO 603 I=1,3810
FF(I,1)=I
FF(I,2)=F(I,2)
FF(I,3)=F(I,3)
603CONTINUE
DT=1.0
IF(NUMB.EQ.INT(NUMB/20000)*20000) THEN
DO 506 I=1,59
XYL(I,2)=XYL(I,2)
XYL(I,3)=XYL(I,3)
506CONTINUE
DO 507 I=60,3751
XYL(I,2)=XYL(I,2)+xis1*VX(I,2)*DT+xis2*FF(I,2)*DT**2
XYL(I,3)=XYL(I,3)+xis1*VY(I,2)*DT+xis2*FF(I,3)*DT**2
507CONTINUE
DO 505 I=3752,3810
XYL(I,2)=XYL(I,2)
XYL(I,3)=XYL(I,3)+0.5
505CONTINUE
else
DO 502 I=1,59
XYL(I,2)=XYL(I,2)
XYL(I,3)=XYL(I,3)
502CONTINUE
DO 516 I=3752,3810
XYL(I,2)=XYL(I,2)
XYL(I,3)=XYL(I,3)
516CONTINUE
DO 508 I=60,3751
XYL(I,2)=XYL(I,2)+xis1*VX(I,2)*DT+xis2*FF(I,2)*DT**2
XYL(I,3)=XYL(I,3)+xis1*VY(I,2)*DT+xis2*FF(I,3)*DT**2
508CONTINUE
end if
IF(NUMB.EQ.INT(NUMB/50)*50) THEN
CALL DLISTFORM(XYL,NFORM,LENG,npart)
END IF
CALL EPENERGY(LENG,NFORM,XYL,PE,ROI,RIJ,npart)
CALL FORCE(LENG,XYL,RIJ,ROI,NFORM,F,npart)
DO 708 I=1,npart
VX(I,2)=VX(I,2)+10.0/2.12*(F(I,2)+FF(I,2))*DT
708 CONTINUE
DO 701 I=21,3790
VY(I,2)=VY(I,2)+10.0/2.12*(F(I,3)+FF(I,3))*DT
701 CONTINUE
VM=0.0
DO 703 I=1,npart
VM=VM+VX(I,2)**2+VY(I,2)**2
703 CONTINUE
EK=0.5*cum*VM !新动能J
T0=EK/(npart-1)/kb
EK=EK/ev
IF(NUMB.EQ.INT(NUMB/10)*10) THEN
ZPVX2=0.0
ZPVY2=0.0
DO 705 I=1,npart
ZPVX2=ZPVX2+VX(I,2)**2
ZPVY2=ZPVY2+VY(I,2)**2
705 CONTINUE
BETAX=DSQRT((npart-1)*KB*T/(cum*ZPVX2))
BETAY=DSQRT((npart-1)*KB*T/(cum*ZPVY2))
DO 711 I=1,npart
VX(I,2)=BETAX*VX(I,2)
VY(I,2)=BETAY*VY(I,2)
711CONTINUE
END IF
WRITE(*,706)NUMB,BETAX,BETAY
706 FORMAT(1X,'NUMB=',I8,2X,'BETAX='F10.6,2X,'BETAY='F10.6)
IF(NUMB.EQ.INT(NUMB/10)*10) THEN
WRITE(11,999)NUMB,PE,EK,T0
999 FORMAT(1X,'NUMB=',2X,I8,2X,'PE=',2X,F18.4,2X,'EK=',2X,F16.4,2X,'T0=',2X,F16.4)
end if
IF(NUMB.EQ.INT(NUMB/10000)*10000) THEN
WRITE(15,9990)(FF(I,2),FF(I,3),I=1880,1900)
9990 FORMAT(2X,'FF(I,2)=',2X,F16.4,2X,'FF(I,3)=',1X,F18.6)
END IF
IF(NUMB.EQ.INT(NUMB/10)*10) THEN
H=XYL(3800,3)-173.205077648163
END IF
WRITE(16,9977)NUMB,H
9977 FORMAT(1X,'NUMB=',I8,2X,F10.6)
DO 123 I=1763,2099,14
E=E+FF(I,2)/0.065025
EE=E/20
123 CONTINUE
WRITE(17,9987)NUMB,EE
9987 FORMAT(I8,2X,F40.20)
NUMB=NUMB+1
IF(NUMB.LE.300000) GOTO 1000
WRITE(10,997)((XYL(I,J),J=1,3),I=1,npart)
997 FORMAT(1X,F10.1,2X,2F16.4)
CLOSE(3)
CLOSE(4)
CLOSE(10)
CLOSE(11)
CLOSE(12)
CLOSE(15)
CLOSE(16)
CLOSE(17)
END
建立初始构型的程序：
SUBROUTINE BCONSTRUCT(XYLZ)
DIMENSION XY(201,20,2),NP(3920),XYN(3920,3),XYLZ(3810,3)
DOUBLE PRECISION XY,XYN,XYLZ
DO 15 I=1,201,2
DO 10 J=1,20
K=(I-1)*20-INT(I/2)+J
10NP(K)=K
15CONTINUE
DO 17 I=2,200,2
DO 16 J=1,19
K=(I-1)*19+INT(I/2)+J
16NP(K)=K
17CONTINUE
DO 11 I=1,3920
XYN(I,1)=NP(I)
11CONTINUE
DO 23 I=1,201,2
DO 22 J=1,20
K=(I-1)*20-INT(I/2)+J
XYN(K,2)=J-1
22XYN(K,3)=(I-1)*SQRT(3.0)/2.0
23CONTINUE
DO 25 I=2,200,2
DO 24 J=1,19
K=(I-1)*19+INT(I/2)+J
XYN(K,2)=1.0/2+(J-1)
24XYN(K,3)=(I-1)*SQRT(3.0)/2.0
25CONTINUE
DO 26 I=1,1762
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I,2)
26XYLZ(I,3)=XYN(I,3)
DO 27 I=1763,1776
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+6,2)
27XYLZ(I,3)=XYN(I+6,3)
DO 28 I=1777,1790
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+11,2)
28XYLZ(I,3)=XYN(I+11,3)
DO 29 I=1791,1804
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+17,2)
29XYLZ(I,3)=XYN(I+17,3)
DO 30 I=1805,1818
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+22,2)
30XYLZ(I,3)=XYN(I+22,3)
DO 31 I=1819,1832
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+28,2)
31XYLZ(I,3)=XYN(I+28,3)
DO 32 I=1833,1846
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+33,2)
32XYLZ(I,3)=XYN(I+33,3)
DO 33 I=1847,1860
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+39,2)
33XYLZ(I,3)=XYN(I+39,3)
DO 34 I=1861,1874
XYLZ(I,1)=XYN(I+44,1)
XYLZ(I,2)=XYN(I+44,2)
34 XYLZ(I,3)=XYN(I+44,3)
DO 35 I=1875,1888
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+50,2)
35XYLZ(I,3)=XYN(I+50,3)
DO 36 I=1889,1902
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+55,2)
36XYLZ(I,3)=XYN(I+55,3)
DO 37 I=1903,1916
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+61,2)
37XYLZ(I,3)=XYN(I+61,3)
DO 38 I=1917,1930
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+66,2)
38XYLZ(I,3)=XYN(I+66,3)
DO 39 I=1931,1944
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+72,2)
39XYLZ(I,3)=XYN(I+72,3)
DO 40 I=1945,1958
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+77,2)
40XYLZ(I,3)=XYN(I+77,3)
DO 41 I=1959,1972
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+83,2)
41XYLZ(I,3)=XYN(I+83,3)
DO 42 I=1973,1986
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+88,2)
42XYLZ(I,3)=XYN(I+88,3)
DO 46 I=1987,2000
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+94,2)
46XYLZ(I,3)=XYN(I+94,3)
DO 43 I=2001,2014
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+99,2)
43XYLZ(I,3)=XYN(I+99,3)
DO 44 I=2015,2028
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+105,2)
44XYLZ(I,3)=XYN(I+105,3)
DO 45 I=2029,3810
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+110,2)
45XYLZ(I,3)=XYN(I+110,3)
END
给粒子随即速度的子程序：
SUBROUTINE CRADOM(VXZ,VYZ,npart)
DIMENSION RAN(npart,2),VXZ(npart,2),VYZ(npart,2)
DOUBLE PRECISION SUMVXZ,SUMVYZ,RAN,VXZ,VYZ,PVXZ,PVYZ
CALL RANDOM_SEED()
CALL RANDOM_NUMBER(RAN)
DO 101 I=1,npart
DO 101 J=1,2
VXZ(I,1)=I
VXZ(I,2)=RAN(I,1)
VYZ(I,1)=I
VYZ(I,2)=RAN(I,2)
101CONTINUE
SUMVXZ=0.0
SUMVYZ=0.0
DO 300 I=1,npart
SUMVXZ=SUMVXZ+VXZ(I,2)
SUMVYZ=SUMVYZ+VYZ(I,2)
300 CONTINUE
PVXZ=SUMVXZ/npart
PVYZ=SUMVYZ/npart
DO 301 I=1,npart
VXZ(I,2)=VXZ(I,2)-PVXZ
VYZ(I,2)=VYZ(I,2)-PVYZ
301 CONTINUE
end
邻域列表的程序：
SUBROUTINE DLISTFORM(XYLZ,NFORMZ,LENG,npart)
DIMENSION XYLZ(npart,3),NFORMZ(npart,40),LENG(npart,2)
DOUBLE PRECISION XYLZ,QQ
DO 220 I=1,npart
DO 220 J=1,2
LENG(I,J)=0
220CONTINUE
DO 200 I=1,npart
DO 200 J=1,40
NFORMZ(I,J)=0
200 CONTINUE
DO 204 I=1,npart
NUM=1
NFORMZ(I,1)=I
LENG(I,1)=I
DO 201 J=1,npart
IF(I.NE.J)THEN
QQ=(XYLZ(I,2)-XYLZ(J,2))**2+(XYLZ(I,3)-XYLZ(J,3))**2
ELSE
GOTO 201
ENDIF
IF(QQ.LE.(5.0625))THEN
NUM=NUM+1
NFORMZ(I,NUM)=J
LENG(I,2)=NUM
ELSE
GOTO 201
ENDIF
201 CONTINUE
204 CONTINUE
END
求构型势能的程序：
SUBROUTINE EPENERGY(LENG,NFORM,XYL,PE,ROI,RIJ,npart)
DIMENSION XYL(npart,3),NFORM(npart,40),LENG(npart,2)
DIMENSION RIJ(npart,40),Q(npart,40),ROI(npart,2),PHIJ(npart,2)
DOUBLE PRECISION XYL,PE,A1,A2,C1,C2,D
DOUBLE PRECISION RIJ,FIJ,PHIJ,PPHIJ,ROI,Q
rc1=1.65
rc2=1.95
A1=8289.5
A2=0.018325
C1=10.727
C2=0.31976
D=13.079
PE=0.0
DO 306 I=1,npart
DO 306 J=1,40
Q(I,J)=0.0
RIJ(I,J)=0.0
306CONTINUE
DO 310 I=1,npart
DO 310 J=1,2
ROI(I,J)=0.0
PHIJ(I,J)=0.0
310 CONTINUE
DO 304 I=1,npart
PHIJ(I,1)=I
Q(I,1)=I
RIJ(I,1)=I
ROI(I,1)=I
DO 305 J=2,LENG(I,2)
NN=NFORM(I,J)
Q(I,J)=(XYL(I,2)-XYL(NN,2))**2+(XYL(I,3)-XYL(NN,3))**2
RIJ(I,J)=DSQRT(Q(I,J))
IF(RIJ(I,J).LT.RC2)THEN
FIJ=A2*(RC2-RIJ(I,J))**2*DEXP(-C2*RIJ(I,J))
ELSE
FIJ=0.0
END IF
ROI(I,2)=ROI(I,2)+FIJ
305CONTINUE
DO 307 J=2,LENG(I,2)
IF(RIJ(I,J).LT.RC1)THEN
PPHIJ=A1*(RC1-RIJ(I,J))**2*DEXP(-C1*RIJ(I,J))
ELSE
PPHIJ=0.0
END IF
PHIJ(I,2)=PHIJ(I,2)+PPHIJ
307CONTINUE
PE=PE+(D*ROI(I,2)*DLOG(ROI(I,2))+0.5*PHIJ(I,2))
304 CONTINUE
END
求构型分子间作用力的程序：
SUBROUTINE FORCE(LENG,XYL,RIJ,ROI,NFORM,F,npart)
DIMENSION XYL(npart,3),NFORM(npart,40),RIJ(npart,40),ROI(npart,2)
DIMENSION FX(npart,40),FY(npart,40),F(npart,3),FBETWEEN(npart,40)
DIMENSION FBETWEENA(npart,40),FBETWEENB(npart,40),LENG(npart,2)
DIMENSION AX(npart,40),AY(npart,40)
DOUBLE PRECISION XYL,RIJ,ROI,FX,FY,F,A1,A2,C1,C2,D
DOUBLE PRECISION AX,AY,FBETWEEN,FBETWEENA,FBETWEENB
rc1=1.65
rc2=1.95
A1=8289.5
A2=0.018325
C1=10.727
C2=0.31976
D=13.079
DO 400 I=1,npart
DO 400 J=1,40
FX(I,J)=0.0
FY(I,J)=0.0
FBETWEEN(I,J)=0.0
FBETWEENA(I,J)=0.0
FBETWEENB(I,J)=0.0
400 CONTINUE
DO 403 I=1,npart
DO 403 J=1,3
F(I,J)=0.0
403 CONTINUE
DO 401 I=1,npart
FX(I,1)=I
FY(I,1)=I
F(I,1)=I
FBETWEEN(I,1)=I
FBETWEENA(I,1)=I
FBETWEENB(I,1)=I
DO 402 J=2,LENG(I,2)
IF(RIJ(I,J).LT.RC2)THEN
M=NFORM(I,J)
FBETWEENA(I,J)=D*(DLOG(ROI(I,2))+DLOG(ROI(M,2))+2)*A2* DEXP(-C2*RIJ(I,J))*(2*(RC2-RIJ(I,J))+C2*(RC2-RIJ(I,J))**2)
ELSE
FBETWEENA(I,J)=0.0
END IF
IF(RIJ(I,J).LT.RC1)THEN
FBETWEENB(I,J)=A1*DEXP(-C1*RIJ(I,J))*(2*(RC1-RIJ(I,J))+C1*(RC1-RIJ(I,J))**2)
ELSE
FBETWEENB(I,J)=0.0
END IF
FBETWEEN(I,J)=FBETWEENA(I,J)+FBETWEENB(I,J)
NIJ=NFORM(I,J)
AX(I,1)=I
AX(I,J)=XYL(NIJ,2)-XYL(I,2)
AY(I,1)=I
AY(I,J)=XYL(NIJ,3)-XYL(I,3)
FX(I,J)=-0.6275817586098889*FBETWEEN(I,J)*AX(I,J)/RIJ(I,J)
FY(I,J)=-0.6275817586098889*FBETWEEN(I,J)*AY(I,J)/RIJ(I,J)
402CONTINUE
DO 404 N=2,LENG(I,2)
F(I,2)=F(I,2)+FX(I,N)
F(I,3)=F(I,3)+FY(I,N)
404 CONTINUE
401 CONTINUE
END
主程序：
DIMENSION VX(3810,2),VY(3810,2),NFORM(3810,40)
DIMENSION LENG(3810,2),F(3810,3),FF(3810,3)
DIMENSION RIJ(3810,40),ROI(3810,2),XYL(3810,3)
DOUBLE PRECISION XYL,VX,VY,VM,T0,FF,PE,RIJ,ROI,F,H,E,EE
DOUBLE PRECISION ZPVX2,ZPVY2,BETAX,BETAY
DOUBLE PRECISION kb,xis1,xis2
OPEN(3,FILE='VX.DAT',STATUS='old')
OPEN(4,FILE='VY.DAT',STATUS='old')
OPEN(10,FILE='XYL.DAT',STATUS='old')
OPEN(11,FILE='RESULT.DAT',STATUS='old')
OPEN(12,FILE='AA.DAT',STATUS='old')
OPEN(15,FILE='F.DAT',STATUS='old')
OPEN(16,FILE='H.DAT',STATUS='old')
OPEN(17,FILE='E.DAT',STATUS='old')
npart=3810
cum=1.06E-25
kb=1.380650524E-23
ev=1.602E-19
T=293
D=0.000000000255
xis1=3.917489092446246E-6
xis2=1.847872213418041E-5
CALL CRADOM(VX,VY,npart)
WRITE(3,102)((VX(I,J),J=1,2),I=1,npart)
WRITE(4,102)((VY(I,J),J=1,2),I=1,npart)
102 FORMAT(1X,f8.1,F14.6)
CALL BCONSTRUCT(XYL)
WRITE(12,*)((XYL(I,J),J=1,3),I=1,npart)
CALL DLISTFORM(XYL,NFORM,LENG,npart)
VM=0.0
DO 702 I=1,npart
VM=VM+VX(I,2)**2+VY(I,2)**2
702 CONTINUE
EK=0.5*cum*VM !新动能J
T0=EK/(npart-1)/kb
EK=EK/ev
CALL EPENERGY(LENG,NFORM,XYL,PE,ROI,RIJ,npart)
WRITE(*,307)PE,EK,T0
307FORMAT(1X,'NUM=0',6X,'PE=',F18.4,2X,'EK=',F10.4,2X,'T0=',F10.4)
CALL FORCE(LENG,XYL,RIJ,ROI,NFORM,F,npart)
NUMB=1 !更新邻域列表及循环控制用
1000 DO 603 I=1,3810
FF(I,1)=I
FF(I,2)=F(I,2)
FF(I,3)=F(I,3)
603CONTINUE
DT=1.0
IF(NUMB.EQ.INT(NUMB/20000)*20000) THEN
DO 506 I=1,59
XYL(I,2)=XYL(I,2)
XYL(I,3)=XYL(I,3)
506CONTINUE
DO 507 I=60,3751
XYL(I,2)=XYL(I,2)+xis1*VX(I,2)*DT+xis2*FF(I,2)*DT**2
XYL(I,3)=XYL(I,3)+xis1*VY(I,2)*DT+xis2*FF(I,3)*DT**2
507CONTINUE
DO 505 I=3752,3810
XYL(I,2)=XYL(I,2)
XYL(I,3)=XYL(I,3)+0.5
505CONTINUE
else
DO 502 I=1,59
XYL(I,2)=XYL(I,2)
XYL(I,3)=XYL(I,3)
502CONTINUE
DO 516 I=3752,3810
XYL(I,2)=XYL(I,2)
XYL(I,3)=XYL(I,3)
516CONTINUE
DO 508 I=60,3751
XYL(I,2)=XYL(I,2)+xis1*VX(I,2)*DT+xis2*FF(I,2)*DT**2
XYL(I,3)=XYL(I,3)+xis1*VY(I,2)*DT+xis2*FF(I,3)*DT**2
508CONTINUE
end if
IF(NUMB.EQ.INT(NUMB/50)*50) THEN
CALL DLISTFORM(XYL,NFORM,LENG,npart)
END IF
CALL EPENERGY(LENG,NFORM,XYL,PE,ROI,RIJ,npart)
CALL FORCE(LENG,XYL,RIJ,ROI,NFORM,F,npart)
DO 708 I=1,npart
VX(I,2)=VX(I,2)+10.0/2.12*(F(I,2)+FF(I,2))*DT
708 CONTINUE
DO 701 I=21,3790
VY(I,2)=VY(I,2)+10.0/2.12*(F(I,3)+FF(I,3))*DT
701 CONTINUE
VM=0.0
DO 703 I=1,npart
VM=VM+VX(I,2)**2+VY(I,2)**2
703 CONTINUE
EK=0.5*cum*VM !新动能J
T0=EK/(npart-1)/kb
EK=EK/ev
IF(NUMB.EQ.INT(NUMB/10)*10) THEN
ZPVX2=0.0
ZPVY2=0.0
DO 705 I=1,npart
ZPVX2=ZPVX2+VX(I,2)**2
ZPVY2=ZPVY2+VY(I,2)**2
705 CONTINUE
BETAX=DSQRT((npart-1)*KB*T/(cum*ZPVX2))
BETAY=DSQRT((npart-1)*KB*T/(cum*ZPVY2))
DO 711 I=1,npart
VX(I,2)=BETAX*VX(I,2)
VY(I,2)=BETAY*VY(I,2)
711CONTINUE
END IF
WRITE(*,706)NUMB,BETAX,BETAY
706 FORMAT(1X,'NUMB=',I8,2X,'BETAX='F10.6,2X,'BETAY='F10.6)
IF(NUMB.EQ.INT(NUMB/10)*10) THEN
WRITE(11,999)NUMB,PE,EK,T0
999 FORMAT(1X,'NUMB=',2X,I8,2X,'PE=',2X,F18.4,2X,'EK=',2X,F16.4,2X,'T0=',2X,F16.4)
end if
IF(NUMB.EQ.INT(NUMB/10000)*10000) THEN
WRITE(15,9990)(FF(I,2),FF(I,3),I=1880,1900)
9990 FORMAT(2X,'FF(I,2)=',2X,F16.4,2X,'FF(I,3)=',1X,F18.6)
END IF
IF(NUMB.EQ.INT(NUMB/10)*10) THEN
H=XYL(3800,3)-173.205077648163
END IF
WRITE(16,9977)NUMB,H
9977 FORMAT(1X,'NUMB=',I8,2X,F10.6)
DO 123 I=1763,2099,14
E=E+FF(I,2)/0.065025
EE=E/20
123 CONTINUE
WRITE(17,9987)NUMB,EE
9987 FORMAT(I8,2X,F40.20)
NUMB=NUMB+1
IF(NUMB.LE.300000) GOTO 1000
WRITE(10,997)((XYL(I,J),J=1,3),I=1,npart)
997 FORMAT(1X,F10.1,2X,2F16.4)
CLOSE(3)
CLOSE(4)
CLOSE(10)
CLOSE(11)
CLOSE(12)
CLOSE(15)
CLOSE(16)
CLOSE(17)
END
建立初始构型的程序：
SUBROUTINE BCONSTRUCT(XYLZ)
DIMENSION XY(201,20,2),NP(3920),XYN(3920,3),XYLZ(3810,3)
DOUBLE PRECISION XY,XYN,XYLZ
DO 15 I=1,201,2
DO 10 J=1,20
K=(I-1)*20-INT(I/2)+J
10NP(K)=K
15CONTINUE
DO 17 I=2,200,2
DO 16 J=1,19
K=(I-1)*19+INT(I/2)+J
16NP(K)=K
17CONTINUE
DO 11 I=1,3920
XYN(I,1)=NP(I)
11CONTINUE
DO 23 I=1,201,2
DO 22 J=1,20
K=(I-1)*20-INT(I/2)+J
XYN(K,2)=J-1
22XYN(K,3)=(I-1)*SQRT(3.0)/2.0
23CONTINUE
DO 25 I=2,200,2
DO 24 J=1,19
K=(I-1)*19+INT(I/2)+J
XYN(K,2)=1.0/2+(J-1)
24XYN(K,3)=(I-1)*SQRT(3.0)/2.0
25CONTINUE
DO 26 I=1,1762
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I,2)
26XYLZ(I,3)=XYN(I,3)
DO 27 I=1763,1776
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+6,2)
27XYLZ(I,3)=XYN(I+6,3)
DO 28 I=1777,1790
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+11,2)
28XYLZ(I,3)=XYN(I+11,3)
DO 29 I=1791,1804
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+17,2)
29XYLZ(I,3)=XYN(I+17,3)
DO 30 I=1805,1818
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+22,2)
30XYLZ(I,3)=XYN(I+22,3)
DO 31 I=1819,1832
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+28,2)
31XYLZ(I,3)=XYN(I+28,3)
DO 32 I=1833,1846
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+33,2)
32XYLZ(I,3)=XYN(I+33,3)
DO 33 I=1847,1860
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+39,2)
33XYLZ(I,3)=XYN(I+39,3)
DO 34 I=1861,1874
XYLZ(I,1)=XYN(I+44,1)
XYLZ(I,2)=XYN(I+44,2)
34 XYLZ(I,3)=XYN(I+44,3)
DO 35 I=1875,1888
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+50,2)
35XYLZ(I,3)=XYN(I+50,3)
DO 36 I=1889,1902
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+55,2)
36XYLZ(I,3)=XYN(I+55,3)
DO 37 I=1903,1916
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+61,2)
37XYLZ(I,3)=XYN(I+61,3)
DO 38 I=1917,1930
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+66,2)
38XYLZ(I,3)=XYN(I+66,3)
DO 39 I=1931,1944
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+72,2)
39XYLZ(I,3)=XYN(I+72,3)
DO 40 I=1945,1958
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+77,2)
40XYLZ(I,3)=XYN(I+77,3)
DO 41 I=1959,1972
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+83,2)
41XYLZ(I,3)=XYN(I+83,3)
DO 42 I=1973,1986
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+88,2)
42XYLZ(I,3)=XYN(I+88,3)
DO 46 I=1987,2000
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+94,2)
46XYLZ(I,3)=XYN(I+94,3)
DO 43 I=2001,2014
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+99,2)
43XYLZ(I,3)=XYN(I+99,3)
DO 44 I=2015,2028
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+105,2)
44XYLZ(I,3)=XYN(I+105,3)
DO 45 I=2029,3810
XYLZ(I,1)=XYN(I,1)
XYLZ(I,2)=XYN(I+110,2)
45XYLZ(I,3)=XYN(I+110,3)
END
给粒子随即速度的子程序：
SUBROUTINE CRADOM(VXZ,VYZ,npart)
DIMENSION RAN(npart,2),VXZ(npart,2),VYZ(npart,2)
DOUBLE PRECISION SUMVXZ,SUMVYZ,RAN,VXZ,VYZ,PVXZ,PVYZ
CALL RANDOM_SEED()
CALL RANDOM_NUMBER(RAN)
DO 101 I=1,npart
DO 101 J=1,2
VXZ(I,1)=I
VXZ(I,2)=RAN(I,1)
VYZ(I,1)=I
VYZ(I,2)=RAN(I,2)
101CONTINUE
SUMVXZ=0.0
SUMVYZ=0.0
DO 300 I=1,npart
SUMVXZ=SUMVXZ+VXZ(I,2)
SUMVYZ=SUMVYZ+VYZ(I,2)
300 CONTINUE
PVXZ=SUMVXZ/npart
PVYZ=SUMVYZ/npart
DO 301 I=1,npart
VXZ(I,2)=VXZ(I,2)-PVXZ
VYZ(I,2)=VYZ(I,2)-PVYZ
301 CONTINUE
end
邻域列表的程序：
SUBROUTINE DLISTFORM(XYLZ,NFORMZ,LENG,npart)
DIMENSION XYLZ(npart,3),NFORMZ(npart,40),LENG(npart,2)
DOUBLE PRECISION XYLZ,QQ
DO 220 I=1,npart
DO 220 J=1,2
LENG(I,J)=0
220CONTINUE
DO 200 I=1,npart
DO 200 J=1,40
NFORMZ(I,J)=0
200 CONTINUE
DO 204 I=1,npart
NUM=1
NFORMZ(I,1)=I
LENG(I,1)=I
DO 201 J=1,npart
IF(I.NE.J)THEN
QQ=(XYLZ(I,2)-XYLZ(J,2))**2+(XYLZ(I,3)-XYLZ(J,3))**2
ELSE
GOTO 201
ENDIF
IF(QQ.LE.(5.0625))THEN
NUM=NUM+1
NFORMZ(I,NUM)=J
LENG(I,2)=NUM
ELSE
GOTO 201
ENDIF
201 CONTINUE
204 CONTINUE
END
求构型势能的程序：
SUBROUTINE EPENERGY(LENG,NFORM,XYL,PE,ROI,RIJ,npart)
DIMENSION XYL(npart,3),NFORM(npart,40),LENG(npart,2)
DIMENSION RIJ(npart,40),Q(npart,40),ROI(npart,2),PHIJ(npart,2)
DOUBLE PRECISION XYL,PE,A1,A2,C1,C2,D
DOUBLE PRECISION RIJ,FIJ,PHIJ,PPHIJ,ROI,Q
rc1=1.65
rc2=1.95
A1=8289.5
A2=0.018325
C1=10.727
C2=0.31976
D=13.079
PE=0.0
DO 306 I=1,npart
DO 306 J=1,40
Q(I,J)=0.0
RIJ(I,J)=0.0
306CONTINUE
DO 310 I=1,npart
DO 310 J=1,2
ROI(I,J)=0.0
PHIJ(I,J)=0.0
310 CONTINUE
DO 304 I=1,npart
PHIJ(I,1)=I
Q(I,1)=I
RIJ(I,1)=I
ROI(I,1)=I
DO 305 J=2,LENG(I,2)
NN=NFORM(I,J)
Q(I,J)=(XYL(I,2)-XYL(NN,2))**2+(XYL(I,3)-XYL(NN,3))**2
RIJ(I,J)=DSQRT(Q(I,J))
IF(RIJ(I,J).LT.RC2)THEN
FIJ=A2*(RC2-RIJ(I,J))**2*DEXP(-C2*RIJ(I,J))
ELSE
FIJ=0.0
END IF
ROI(I,2)=ROI(I,2)+FIJ
305CONTINUE
DO 307 J=2,LENG(I,2)
IF(RIJ(I,J).LT.RC1)THEN
PPHIJ=A1*(RC1-RIJ(I,J))**2*DEXP(-C1*RIJ(I,J))
ELSE
PPHIJ=0.0
END IF
PHIJ(I,2)=PHIJ(I,2)+PPHIJ
307CONTINUE
PE=PE+(D*ROI(I,2)*DLOG(ROI(I,2))+0.5*PHIJ(I,2))
304 CONTINUE
END
求构型分子间作用力的程序：
SUBROUTINE FORCE(LENG,XYL,RIJ,ROI,NFORM,F,npart)
DIMENSION XYL(npart,3),NFORM(npart,40),RIJ(npart,40),ROI(npart,2)
DIMENSION FX(npart,40),FY(npart,40),F(npart,3),FBETWEEN(npart,40)
DIMENSION FBETWEENA(npart,40),FBETWEENB(npart,40),LENG(npart,2)
DIMENSION AX(npart,40),AY(npart,40)
DOUBLE PRECISION XYL,RIJ,ROI,FX,FY,F,A1,A2,C1,C2,D
DOUBLE PRECISION AX,AY,FBETWEEN,FBETWEENA,FBETWEENB
rc1=1.65
rc2=1.95
A1=8289.5
A2=0.018325
C1=10.727
C2=0.31976
D=13.079
DO 400 I=1,npart
DO 400 J=1,40
FX(I,J)=0.0
FY(I,J)=0.0
FBETWEEN(I,J)=0.0
FBETWEENA(I,J)=0.0
FBETWEENB(I,J)=0.0
400 CONTINUE
DO 403 I=1,npart
DO 403 J=1,3
F(I,J)=0.0
403 CONTINUE
DO 401 I=1,npart
FX(I,1)=I
FY(I,1)=I
F(I,1)=I
FBETWEEN(I,1)=I
FBETWEENA(I,1)=I
FBETWEENB(I,1)=I
DO 402 J=2,LENG(I,2)
IF(RIJ(I,J).LT.RC2)THEN
M=NFORM(I,J)
FBETWEENA(I,J)=D*(DLOG(ROI(I,2))+DLOG(ROI(M,2))+2)*A2* DEXP(-C2*RIJ(I,J))*(2*(RC2-RIJ(I,J))+C2*(RC2-RIJ(I,J))**2)
ELSE
FBETWEENA(I,J)=0.0
END IF
IF(RIJ(I,J).LT.RC1)THEN
FBETWEENB(I,J)=A1*DEXP(-C1*RIJ(I,J))*(2*(RC1-RIJ(I,J))+C1*(RC1-RIJ(I,J))**2)
ELSE
FBETWEENB(I,J)=0.0
END IF
FBETWEEN(I,J)=FBETWEENA(I,J)+FBETWEENB(I,J)
NIJ=NFORM(I,J)
AX(I,1)=I
AX(I,J)=XYL(NIJ,2)-XYL(I,2)
AY(I,1)=I
AY(I,J)=XYL(NIJ,3)-XYL(I,3)
FX(I,J)=-0.6275817586098889*FBETWEEN(I,J)*AX(I,J)/RIJ(I,J)
FY(I,J)=-0.6275817586098889*FBETWEEN(I,J)*AY(I,J)/RIJ(I,J)
402CONTINUE
DO 404 N=2,LENG(I,2)
F(I,2)=F(I,2)+FX(I,N)
F(I,3)=F(I,3)+FY(I,N)
404 CONTINUE
401 CONTINUE
END