-
Notifications
You must be signed in to change notification settings - Fork 14
/
Copy pathgear.f
201 lines (165 loc) · 8.51 KB
/
gear.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
********************************************************************************
** FICHE F.2. GEAR 5-VALUE PREDICTOR-CORRECTOR ALGORITHM **
** This FORTRAN code is intended to illustrate points made in the text. **
** To our knowledge it works correctly. However it is the responsibility of **
** the user to test it, if it is to be used in a research application. **
********************************************************************************
C *******************************************************************
C ** GEAR 5-VALUE PREDICTOR-CORRECTOR ALGORITHM FOR TRANSLATION. **
C ** **
C ** REFERENCE: **
C ** **
C ** GEAR, NUMERICAL INITIAL VALUE PROBLEMS IN ORDINARY **
C ** DIFFERENTIAL EQUATIONS (PRENTICE-HALL,1971). **
C ** **
C ** SUPPLIED ROUTINES: **
C ** **
C ** SUBROUTINE PREDIC ( DT ) **
C ** PREDICTS THE NEW POSITIONS, VELOCITIES, ETC. **
C ** SUBROUTINE CORREC ( DT, M, K ) **
C ** CORRECTS THE POSITIONS, VELOCITIES ETC. USING GEAR METHOD **
C ** **
C ** PRINCIPAL VARIABLES: **
C ** **
C ** INTEGER N NUMBER OF MOLECULES **
C ** REAL DT TIMESTEP **
C ** REAL RX(N),RY(N),RZ(N) POSITIONS **
C ** REAL VX(N),VY(N),VZ(N) VELOCITIES **
C ** REAL AX(N),AY(N),AZ(N) ACCELERATIONS **
C ** REAL BX(N),BY(N),BZ(N) THIRD DERIVATIVES **
C ** REAL CX(N),CY(N),CZ(N) FOURTH DERIVATIVES **
C ** REAL FX(N),FY(N),FZ(N) FORCES **
C ** **
C ** USAGE: **
C ** **
C ** AT EACH TIMESTEP, CALL PREDIC, FORCE, CORREC IN ORDER **
C ** FOLLOWED BY ACCUMULATION OF THERMODYNAMIC QUANTITIES. **
C ** THE FORCE ROUTINE (NOT SUPPLIED HERE: SEE F.17) CALCULATES **
C ** POTENTIAL ENERGY AND FORCES ON ALL ATOMS. **
C *******************************************************************
SUBROUTINE PREDIC ( DT )
COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, AX, AY, AZ,
: BX, BY, BZ, CX, CY, CZ, FX, FY, FZ
C *******************************************************************
C ** PREDICTOR ROUTINE **
C ** **
C ** IN TIMESTEP-SCALED VARIABLES THE PREDICTOR IS THE PASCAL **
C ** TRIANGLE MATRIX. IN UNSCALED VARIABLES IT IS A TAYLOR SERIES **
C ** **
C ** USAGE: **
C ** **
C ** PREDIC IS CALLED TO ADVANCE THE COORDINATES, VELOCITIES ETC. **
C ** BY ONE TIMESTEP DT, PRIOR TO FORCE EVALUATION. **
C *******************************************************************
INTEGER N
PARAMETER ( N = 108 )
REAL DT
REAL RX(N), RY(N), RZ(N)
REAL VX(N), VY(N), VZ(N)
REAL AX(N), AY(N), AZ(N)
REAL BX(N), BY(N), BZ(N)
REAL CX(N), CY(N), CZ(N)
REAL FX(N), FY(N), FZ(N)
INTEGER I
REAL C1, C2, C3, C4
C *******************************************************************
C1 = DT
C2 = C1 * DT / 2.0
C3 = C2 * DT / 3.0
C4 = C3 * DT / 4.0
DO 100 I = 1, N
RX(I) = RX(I) + C1*VX(I) + C2*AX(I) + C3*BX(I) + C4*CX(I)
RY(I) = RY(I) + C1*VY(I) + C2*AY(I) + C3*BY(I) + C4*CY(I)
RZ(I) = RZ(I) + C1*VZ(I) + C2*AZ(I) + C3*BZ(I) + C4*CZ(I)
VX(I) = VX(I) + C1*AX(I) + C2*BX(I) + C3*CX(I)
VY(I) = VY(I) + C1*AY(I) + C2*BY(I) + C3*CY(I)
VZ(I) = VZ(I) + C1*AZ(I) + C2*BZ(I) + C3*CZ(I)
AX(I) = AX(I) + C1*BX(I) + C2*CX(I)
AY(I) = AY(I) + C1*BY(I) + C2*CY(I)
AZ(I) = AZ(I) + C1*BZ(I) + C2*CZ(I)
BX(I) = BX(I) + C1*CX(I)
BY(I) = BY(I) + C1*CY(I)
BZ(I) = BZ(I) + C1*CZ(I)
100 CONTINUE
RETURN
END
SUBROUTINE CORREC ( DT, M, K )
COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, AX, AY, AZ,
: BX, BY, BZ, CX, CY, CZ, FX, FY, FZ
C *******************************************************************
C ** CORRECTOR ROUTINE **
C ** **
C ** CORRECTS POSITIONS, VELOCITIES ETC. AFTER FORCE EVALUATION. **
C ** IN TIMESTEP-SCALED VARIABLES THE NUMERICAL COEFFICIENTS ARE **
C ** GIVEN BY GEAR (REF ABOVE): 19/120, 3/4, 1, 1/2, 1/12. **
C ** IN UNSCALED FORM THESE MUST BE MULTIPLIED BY FACTORS **
C ** INVOLVING THE TIMESTEP AS SHOWN HERE. **
C ** **
C ** PRINCIPAL VARIABLES: **
C ** **
C ** REAL M ATOMIC MASS **
C ** REAL K KINETIC ENERGY PER ATOM **
C ** REAL GEAR0,GEAR1,GEAR2,GEAR3 GEAR COEFFICIENTS **
C ** **
C ** USAGE: **
C ** **
C ** IT IS ASSUMED THAT INTERMOLECULAR FORCES HAVE BEEN CALCULATED **
C ** AND STORED IN FX,FY,FZ. CORREC SIMPLY APPLIES THE CORRECTOR **
C ** EQUATIONS BASED ON THE DIFFERENCES BETWEEN PREDICTED AND **
C ** EVALUATED ACCELERATIONS. IT ALSO CALCULATES KINETIC ENERGY. **
C *******************************************************************
INTEGER N
PARAMETER ( N = 108)
REAL DT
REAL RX(N), RY(N), RZ(N)
REAL VX(N), VY(N), VZ(N)
REAL AX(N), AY(N), AZ(N)
REAL BX(N), BY(N), BZ(N)
REAL CX(N), CY(N), CZ(N)
REAL FX(N), FY(N), FZ(N)
REAL M, K
INTEGER I
REAL AXI, AYI, AZI
REAL CORRX, CORRY, CORRZ
REAL C1, C2, C3, C4
REAL CR, CV, CB, CC
REAL GEAR0, GEAR1, GEAR3, GEAR4
PARAMETER ( GEAR0 = 19.0 / 120.0, GEAR1 = 3.0 / 4.0,
: GEAR3 = 1.0 / 2.0, GEAR4 = 1.0 / 12.0 )
C *******************************************************************
C1 = DT
C2 = C1 * DT / 2.0
C3 = C2 * DT / 3.0
C4 = C3 * DT / 4.0
CR = GEAR0 * C2
CV = GEAR1 * C2 / C1
CB = GEAR3 * C2 / C3
CC = GEAR4 * C2 / C4
K = 0.0
DO 100 I = 1, N
AXI = FX(I) / M
AYI = FY(I) / M
AZI = FZ(I) / M
CORRX = AXI - AX(I)
CORRY = AYI - AY(I)
CORRZ = AZI - AZ(I)
RX(I) = RX(I) + CR * CORRX
RY(I) = RY(I) + CR * CORRY
RZ(I) = RZ(I) + CR * CORRZ
VX(I) = VX(I) + CV * CORRX
VY(I) = VY(I) + CV * CORRY
VZ(I) = VZ(I) + CV * CORRZ
AX(I) = AXI
AY(I) = AYI
AZ(I) = AZI
BX(I) = BX(I) + CB * CORRX
BY(I) = BY(I) + CB * CORRY
BZ(I) = BZ(I) + CB * CORRZ
CX(I) = CX(I) + CC * CORRX
CY(I) = CY(I) + CC * CORRY
CZ(I) = CZ(I) + CC * CORRZ
K = K + VX(I) ** 2 + VY(I) ** 2 + VZ(I) ** 2
100 CONTINUE
K = 0.5 * M * K
RETURN
END