-
Notifications
You must be signed in to change notification settings - Fork 14
/
Copy pathverlet-list.f
366 lines (275 loc) · 13.7 KB
/
verlet-list.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
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
************************************************************************************************************************
** FICHE F.19. THE VERLET NEIGHBOUR LIST **
** 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. **
************************************************************************************************************************
** *******************************************************************
** ** Modified by MPA 22/5/91 correcting error in force subroutine **
C *******************************************************************
C ** FORCE ROUTINE USING A VERLET NEIGHBOUR LIST. **
C ** **
C ** PRINCIPAL VARIABLES: **
C ** **
C ** INTEGER N NUMBER OF ATOMS **
C ** REAL RCUT CUTOFF DISTANCE FOR THE FORCE **
C ** REAL RLIST OUTER RADIUS OF THE VERLET LIST **
C ** REAL SIGMA LENNARD JONES LENGTH PARAMETER **
C ** REAL V POTENTIAL ENERGY **
C ** REAL W VIRIAL **
C ** REAL RX(N),RY(N),RZ(N) ATOM POSITIONS **
C ** REAL FX(N),FY(N),FZ(N) FORCE ON AN ATOM **
C ** LOGICAL UPDATE IF TRUE THE LIST IS UPDATED **
C ** **
C ** ROUTINES SUPPLIED: **
C ** **
C ** SUBROUTINE FORCE ( RCUT, RLIST, SIGMA, UPDATE, V, W ) **
C ** CALCULATE FORCES USING THE VERLET LIST, UPDATES THE LIST **
C ** SUBROUTINE CHECK ( RCUT, RLIST, UPDATE ) **
C ** SETS UPDATE TO TRUE WHEN THE LIST NEEDS TO BE UPDATED **
C ** SUBROUTINE SAVE **
C ** SAVES CURRENT CONFIGURATION FOR FUTURE CHECKING. **
C ** **
C ** USAGE: **
C ** **
C ** AT THE START OF A RUN, FORCE IS CALLED WITH UPDATE SET TRUE. **
C ** THIS SETS UP THE INITIAL VERLET LIST AND SAVES THE POSITIONS. **
C ** THEREAFTER CHECK IS CALLED JUST BEFORE EACH CALL OF FORCE TO **
C ** DECIDE WHETHER OR NOT A LIST UPDATE IS NECESSARY. **
C ** THESE ROUTINES CAN BE USED IN A CONVENTIONAL MD PROGRAM AND **
C ** CAN BE EASILY ADAPTED FOR USE IN AN MC PROGRAM. FORCES IS **
C ** SPECIFIC TO A FLUID OF LENNARD JONES ATOMS. **
C *******************************************************************
SUBROUTINE FORCE ( RCUT, RLIST, SIGMA, UPDATE, V, W )
COMMON / BLOCK1 / RX, RY, RZ, FX, FY, FZ
COMMON / BLOCK2 / POINT, LIST
C *******************************************************************
C ** CALCULATES THE FORCE ON AN ATOM IN A LENNARD-JONES LIQUID **
C ** **
C ** PRINCIPAL VARIABLES: **
C ** **
C ** REAL RCUT CUTOFF DISTANCE FOR THE FORCE **
C ** REAL RLIST OUTER RADIUS OF THE VERLET LIST **
C ** REAL SIGMA LENNARD JONES LENGTH PARAMETER **
C ** REAL V POTENTIAL ENERGY **
C ** REAL W VIRIAL **
C ** REAL RX(N),RY(N),RZ(N) ATOM POSITIONS **
C ** REAL FX(N),FY(N),FZ(N) FORCE ON AN ATOM **
C ** INTEGER POINT(N) INDEX TO THE NEIGHBOUR LIST **
C ** INTEGER LIST(MAXNAB) VERLET NEIGHBOUR LIST **
C ** LOGICAL UPDATE IF TRUE THE LIST IS UPDATED **
C ** **
C ** USAGE: **
C ** **
C ** FORCE IS CALLED IN TWO MODES. IF UPDATE IS TRUE, THE VERLET **
C ** LIST IS UPDATED AND THE FORCES CALCULATED AT THE SAME TIME. **
C ** THE CURRENT CONFIGURATION IS SAVED FOR FUTURE CHECKING. **
C ** IF UPDATE IS FALSE, THE VERLET LIST IS USED TO FIND THE **
C ** NEIGHBOURS OF A GIVEN ATOM AND CALCULATE THE FORCES. **
C ** THE ATOMS ARE IN A BOX OF UNIT LENGTH, CENTRED AT THE ORIGIN. **
C *******************************************************************
INTEGER N, MAXNAB
PARAMETER ( N = 108, MAXNAB = N * 25 )
REAL RX(N), RY(N), RZ(N), FX(N), FY(N), FZ(N)
REAL RLIST, RCUT, SIGMA, V, W
INTEGER POINT(N), LIST(MAXNAB)
LOGICAL UPDATE
REAL RXI, RYI, RZI, FXI, FYI, FZI
REAL RIJSQ, SR2, SR6, VIJ, WIJ, FIJ
REAL SIGSQ, RCUTSQ, RLSTSQ
REAL RXIJ, RYIJ, RZIJ, FXIJ, FYIJ, FZIJ
INTEGER I, J, NLIST
INTEGER JBEG, JEND, JNAB
C *******************************************************************
RLSTSQ = RLIST * RLIST
SIGSQ = SIGMA * SIGMA
RCUTSQ = RCUT * RCUT
C ** ZERO FORCES **
DO 10 I = 1, N
FX(I) = 0.0
FY(I) = 0.0
FZ(I) = 0.0
10 CONTINUE
V = 0.0
W = 0.0
IF ( UPDATE ) THEN
C ** SAVE CURRENT CONFIGURATION, CONSTRUCT **
C ** NEIGHBOUR LIST AND CALCULATE FORCES **
CALL SAVE
NLIST = 0
DO 100 I = 1, N - 1
POINT(I) = NLIST + 1
RXI = RX(I)
RYI = RY(I)
RZI = RZ(I)
FXI = FX(I)
FYI = FY(I)
FZI = FZ(I)
DO 99 J = I + 1, N
RXIJ = RXI - RX(J)
RYIJ = RYI - RY(J)
RZIJ = RZI - RZ(J)
RXIJ = RXIJ - ANINT ( RXIJ )
RYIJ = RYIJ - ANINT ( RYIJ )
RZIJ = RZIJ - ANINT ( RZIJ )
RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ
IF ( RIJSQ .LT. RLSTSQ ) THEN
NLIST = NLIST + 1
LIST(NLIST) = J
C ** REMOVE THIS CHECK IF MAXNAB IS APPROPRIATE **
IF ( NLIST .EQ. MAXNAB ) STOP 'LIST TOO SMALL'
IF ( RIJSQ .LT. RCUTSQ ) THEN
SR2 = SIGSQ / RIJSQ
SR6 = SR2 * SR2 * SR2
VIJ = SR6 * ( SR6 - 1.0 )
WIJ = SR6 * ( SR6 - 0.5 )
V = V + VIJ
W = W + WIJ
FIJ = WIJ / RIJSQ
FXIJ = RXIJ * FIJ
FYIJ = RYIJ * FIJ
FZIJ = RZIJ * FIJ
FXI = FXI + FXIJ
FYI = FYI + FYIJ
FZI = FZI + FZIJ
FX(J) = FX(J) - FXIJ
FY(J) = FY(J) - FYIJ
FZ(J) = FZ(J) - FZIJ
ENDIF
ENDIF
99 CONTINUE
FX(I) = FXI
FY(I) = FYI
FZ(I) = FZI
100 CONTINUE
POINT(N) = NLIST + 1
ELSE
C ** USE THE LIST TO FIND THE NEIGHBOURS **
DO 200 I = 1, N - 1
C ** ERROR CORRECTED HERE 22/5/91 BY MPA **
C ** THANKS TO WERNER LOOSE FOR SPOTTING IT **
C ** AMAZING THAT NO-ONE POINTED IT OUT BEFORE **
C ** ERROR ** JBEG = LIST(I)
C ** ERROR ** JEND = LIST(I+1) - 1
JBEG = POINT(I)
JEND = POINT(I+1) - 1
C ** CHECK THAT ATOM I HAS NEIGHBOURS **
IF( JBEG .LE. JEND ) THEN
RXI = RX(I)
RYI = RY(I)
RZI = RZ(I)
FXI = FX(I)
FYI = FY(I)
FZI = FZ(I)
DO 199 JNAB = JBEG, JEND
J = LIST(JNAB)
RXIJ = RXI - RX(J)
RYIJ = RYI - RY(J)
RZIJ = RZI - RZ(J)
RXIJ = RXIJ - ANINT( RXIJ )
RYIJ = RYIJ - ANINT( RYIJ )
RZIJ = RZIJ - ANINT( RZIJ )
RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ
IF ( RIJSQ .LT. RCUTSQ ) THEN
SR2 = SIGSQ / RIJSQ
SR6 = SR2 * SR2 * SR2
VIJ = SR6 * ( SR6 - 1.0 )
WIJ = SR6 * ( SR6 - 0.5 )
V = V + VIJ
W = W + WIJ
FIJ = WIJ / RIJSQ
FXIJ = RXIJ * FIJ
FYIJ = RYIJ * FIJ
FZIJ = RZIJ * FIJ
FXI = FXI + FXIJ
FYI = FYI + FYIJ
FZI = FZI + FZIJ
FX(J) = FX(J) - FXIJ
FY(J) = FY(J) - FYIJ
FZ(J) = FZ(J) - FZIJ
ENDIF
199 CONTINUE
FX(I) = FXI
FY(I) = FYI
FZ(I) = FZI
ENDIF
200 CONTINUE
ENDIF
V = 4.0 * V
W = 48.0 * W / 3.0
DO 300 I = 1, N
FX(I) = 48.0 * FX(I)
FY(I) = 48.0 * FY(I)
FZ(I) = 48.0 * FZ(I)
300 CONTINUE
RETURN
END
SUBROUTINE CHECK ( RCUT, RLIST, UPDATE )
COMMON / BLOCK1 / RX, RY, RZ, FX, FY, FZ
COMMON / BLOCK3 / RX0, RY0, RZ0
C *******************************************************************
C ** DECIDES WHETHER THE LIST NEEDS TO BE RECONSTRUCTED. **
C ** **
C ** PRINCIPAL VARIABLES: **
C ** **
C ** REAL RX(N),RY(N),RZ(N) ATOM POSITIONS **
C ** REAL RX0(N),RY0(N),RZ0(N) COORDINATES AT LAST UPDATE **
C ** REAL RLIST RADIUS OF VERLET LIST **
C ** REAL RCUT CUTOFF DISTANCE FOR FORCES **
C ** REAL DISPMX LARGEST DISPLACEMENT **
C ** INTEGER N NUMBER OF ATOMS **
C ** LOGICAL UPDATE IF TRUE THE LIST IS UPDATED **
C ** **
C ** USAGE: **
C ** **
C ** CHECK IS CALLED TO SET UPDATE BEFORE EVERY CALL TO FORCE. **
C *******************************************************************
INTEGER N
PARAMETER ( N = 108 )
REAL RX(N), RY(N), RZ(N), FX(N), FY(N), FZ(N)
REAL RX0(N), RY0(N), RZ0(N)
REAL RCUT, RLIST
LOGICAL UPDATE
REAL DISPMX
INTEGER I
C *******************************************************************
C ** CALCULATE MAXIMUM DISPLACEMENT SINCE LAST UPDATE **
DISPMX = 0.0
DO 30 I = 1, N
DISPMX = MAX ( ABS ( RX(I) - RX0(I) ), DISPMX )
DISPMX = MAX ( ABS ( RY(I) - RY0(I) ), DISPMX )
DISPMX = MAX ( ABS ( RZ(I) - RZ0(I) ), DISPMX )
30 CONTINUE
C ** A CONSERVATIVE TEST OF THE LIST SKIN CROSSING **
DISPMX = 2.0 * SQRT ( 3.0 * DISPMX ** 2 )
UPDATE = ( DISPMX .GT. ( RLIST - RCUT ) )
RETURN
END
SUBROUTINE SAVE
COMMON / BLOCK1 / RX, RY, RZ, FX, FY, FZ
COMMON / BLOCK3 / RX0, RY0, RZ0
C *******************************************************************
C ** SAVES CURRENT CONFIGURATION FOR FUTURE CHECKING. **
C ** **
C ** PRINCIPAL VARIABLES: **
C ** **
C ** REAL RX(N),RY(N),RZ(N) ATOM POSITIONS **
C ** REAL RX0(N),RY0(N),RZ0(N) STORAGE LOCATIONS FOR SAVE **
C ** INTEGER N NUMBER OF ATOMS **
C ** **
C ** USAGE: **
C ** **
C ** SAVE IS CALLED WHENEVER THE NEW VERLET LIST IS CONSTRUCTED. **
C *******************************************************************
INTEGER N
PARAMETER ( N = 108 )
REAL RX(N), RY(N), RZ(N), FX(N), FY(N), FZ(N)
REAL RX0(N), RY0(N), RZ0(N)
INTEGER I
C *******************************************************************
DO 100 I = 1, N
RX0(I) = RX(I)
RY0(I) = RY(I)
RZ0(I) = RZ(I)
100 CONTINUE
RETURN
END