1 |
|
2 |
subroutine long_range_force ( check_exclude, & |
3 |
pair_i, pair_j, n_pairs, natoms, sigma, epslon, & |
4 |
rcut, rlist, box_x, box_y, box_z, & |
5 |
rx, ry, rz, fx, fy, fz, & |
6 |
v, & |
7 |
update, point, list, maxnab, & |
8 |
rx0, ry0, rz0, & |
9 |
ux, uy, uz, & |
10 |
tx, ty, tz, & |
11 |
ex, ey, ez, & |
12 |
mu, rrf, rtaper, & |
13 |
dielectric, & |
14 |
is_dipole, is_vdw, is_lj, is_ssd, & |
15 |
Axx, Axy, Axz, & |
16 |
Ayx, Ayy, Ayz, & |
17 |
Azx, Azy, Azz ) |
18 |
|
19 |
implicit none |
20 |
|
21 |
! This routine implements the Verlet neighbor list for all of |
22 |
! the long range force calculations. After a pair has been |
23 |
! determined to lie within the neighbor list, the pair is passed |
24 |
! to the appropriate force functions. |
25 |
|
26 |
|
27 |
|
28 |
! Passed Arguments |
29 |
|
30 |
integer :: n_pairs ! the number of excluded pairs |
31 |
integer :: natoms ! the number of atoms in the system |
32 |
integer :: maxnab ! the max number of neighbors for the verlet list |
33 |
|
34 |
logical(kind=2) :: check_exclude ! boolean to check for exclusion pairs |
35 |
logical(kind=2) :: update ! boolean to update the neighbor list |
36 |
|
37 |
double precision rcut ! the VDW/LJ cutoff radius |
38 |
double precision rlist ! Verlet list cutoff |
39 |
double precision box_x, box_y, box_z ! periodic boundry conditions |
40 |
double precision rrf ! dipole reaction field cutoff |
41 |
double precision rtaper ! the taper radius of the rxn field |
42 |
double precision dielectric ! the dielectric cnst. of the rxn field |
43 |
double precision v ! the potential energy |
44 |
|
45 |
|
46 |
! Passed Arrays |
47 |
|
48 |
integer, dimension(n_pairs):: pair_i, pair_j ! the excluded pairs i and j |
49 |
integer, dimension(natoms) :: point ! the Verlet list ptr array |
50 |
integer, dimension(maxnab) :: list ! the verlet list |
51 |
|
52 |
logical(kind=2), dimension(natoms) :: is_dipole ! dipole boolean array |
53 |
logical(kind=2), dimension(natoms) :: is_vdw ! VDW boolean array |
54 |
logical(kind=2), dimension(natoms) :: is_lj ! LJ boolean array |
55 |
logical(kind=2), dimension(natoms) :: is_ssd ! ssd boolean array |
56 |
|
57 |
double precision, dimension(natoms) :: sigma ! VDW/LJ distance prmtr. |
58 |
double precision, dimension(natoms) :: epslon ! VDW/LJ well depth prmtr. |
59 |
double precision, dimension(natoms) :: rx, ry, rz ! positions |
60 |
double precision, dimension(natoms) :: fx, fy, fz ! forces |
61 |
double precision, dimension(natoms) :: rx0, ry0, rz0 ! intial verlet positions |
62 |
double precision, dimension(natoms) :: ux, uy, uz ! dipole unit vectors |
63 |
double precision, dimension(natoms) :: tx, ty, tz ! torsions |
64 |
double precision, dimension(natoms) :: ex, ey, ez ! reacion field |
65 |
double precision, dimension(natoms) :: mu ! dipole moments |
66 |
double precision, dimension(natoms) :: Axx, Axy, Axz ! rotational array |
67 |
double precision, dimension(natoms) :: Ayx, Ayy, Ayz ! |
68 |
double precision, dimension(natoms) :: Azx, Azy, Azz ! |
69 |
|
70 |
! local variables |
71 |
|
72 |
double precision rxi, ryi, rzi |
73 |
double precision rcutsq, rlstsq, rrfsq, rijsq, rij |
74 |
double precision rxij, ryij, rzij |
75 |
double precision prerf ! a reaction field preterm |
76 |
|
77 |
double precision, dimension(9,natoms) :: A |
78 |
|
79 |
integer :: i, j, k, ix, jx, nlist |
80 |
integer :: jbeg, jend, jnab |
81 |
|
82 |
logical :: exclude_temp1, exclude_temp2, exclude |
83 |
|
84 |
|
85 |
!******************************************************************* |
86 |
|
87 |
do i=1, natoms |
88 |
A(1,i) = Axx(i) |
89 |
A(2,i) = Axy(i) |
90 |
A(3,i) = Axz(i) |
91 |
|
92 |
A(4,i) = Ayx(i) |
93 |
A(5,i) = Ayy(i) |
94 |
A(6,i) = Ayz(i) |
95 |
|
96 |
A(7,i) = Azx(i) |
97 |
A(8,i) = Azy(i) |
98 |
A(9,i) = Azz(i) |
99 |
end do |
100 |
|
101 |
|
102 |
rlstsq = rlist * rlist |
103 |
rcutsq = rcut * rcut |
104 |
rrfsq = rrf * rrf |
105 |
|
106 |
prerf = 14.39257d0 * 2.0d0 * ( dielectric - 1.0d0 ) / & |
107 |
( ( 2.0d0 * dielectric + 1.0d0 ) * rrfsq * rrf ) |
108 |
|
109 |
! ** zero forces ** |
110 |
|
111 |
do I = 1, natoms |
112 |
|
113 |
fx(i) = 0.0d0 |
114 |
fy(i) = 0.0d0 |
115 |
fz(i) = 0.0d0 |
116 |
|
117 |
tx(i) = 0.0d0 |
118 |
ty(i) = 0.0d0 |
119 |
tz(i) = 0.0d0 |
120 |
|
121 |
ex(i) = 0.0d0 |
122 |
ey(i) = 0.0d0 |
123 |
ez(i) = 0.0d0 |
124 |
|
125 |
end do |
126 |
|
127 |
V = 0.0d0 |
128 |
|
129 |
if ( update ) then |
130 |
|
131 |
|
132 |
! ** save current configuration, construct ** |
133 |
! ** neighbour list and calculate forces ** |
134 |
|
135 |
call long_range_save(natoms,rx,ry,rz,rx0,ry0,rz0) |
136 |
|
137 |
nlist = 0 |
138 |
|
139 |
do i = 1, natoms - 1 |
140 |
|
141 |
point(i) = nlist + 1 |
142 |
|
143 |
rxi = rx(i) |
144 |
ryi = ry(i) |
145 |
rzi = rz(i) |
146 |
|
147 |
do j = i + 1, natoms |
148 |
|
149 |
rxij = rx(j) - rxi |
150 |
ryij = ry(j) - ryi |
151 |
rzij = rz(j) - rzi |
152 |
|
153 |
exclude = .false. |
154 |
|
155 |
if( check_exclude ) then |
156 |
|
157 |
do k = 1, n_pairs |
158 |
|
159 |
! add 1 for c -> fortran indexing |
160 |
ix = pair_i(k) + 1 |
161 |
jx = pair_j(k) + 1 |
162 |
exclude_temp1 = (i.eq.ix).and.(j.eq.jx) |
163 |
exclude_temp2 = (i.eq.jx).and.(j.eq.ix) |
164 |
|
165 |
if(exclude_temp1.or.exclude_temp2) then |
166 |
|
167 |
exclude = .true. |
168 |
endif |
169 |
enddo |
170 |
endif |
171 |
|
172 |
if(.not.exclude) then |
173 |
|
174 |
rxij = rxij - box_x * dsign( 1.0d0, rxij ) & |
175 |
* int( (dabs( rxij / box_x ) + 0.5d0) ) |
176 |
ryij = ryij - box_y * dsign( 1.0d0, ryij ) & |
177 |
* int( (dabs( ryij / box_y ) + 0.5d0) ) |
178 |
rzij = rzij - box_z * dsign( 1.0d0, rzij ) & |
179 |
* int( (dabs( rzij / box_z ) + 0.5d0) ) |
180 |
|
181 |
rijsq = rxij * rxij + ryij * ryij + rzij * rzij |
182 |
|
183 |
if ( rijsq .lt. rlstsq ) then |
184 |
|
185 |
nlist = nlist + 1 |
186 |
list(nlist) = j |
187 |
|
188 |
! ** remove this check if maxnab is appropriate ** |
189 |
|
190 |
if ( nlist .eq. maxnab ) then |
191 |
write(*,*) i, j, nlist, maxnab, rijsq, rlstsq |
192 |
stop 'list too small' |
193 |
endif |
194 |
|
195 |
if ( rijsq .lt. rcutsq ) then |
196 |
|
197 |
if ( is_vdw(i) .and. is_vdw(j) ) then |
198 |
call force_VDW( i, j, rcutsq, rijsq, sigma, epslon, & |
199 |
v, fx, fy, fz, rxij, ryij, rzij, natoms ) |
200 |
endif |
201 |
|
202 |
if ( is_lj(i) .and. is_lj(j) ) then |
203 |
call force_lj( i, j, rcutsq, rijsq, sigma, epslon, v, & |
204 |
fx, fy, fz, rxij, ryij, rzij, natoms ) |
205 |
endif |
206 |
|
207 |
if( is_ssd(i) .and. is_ssd(j) ) then |
208 |
|
209 |
rij = dsqrt(rijsq); |
210 |
|
211 |
call ssd_forces( natoms, i, j, i, j, rxij, ryij, rzij, & |
212 |
rij, v, rijsq, fx, fy, fz, tx, ty, tz, A ) |
213 |
end if |
214 |
|
215 |
endif |
216 |
|
217 |
if( is_dipole(i) .and. is_dipole(j) ) then |
218 |
if( rijsq .lt. rrfsq ) then |
219 |
|
220 |
call dipole_rf( i, j, rrfsq, rijsq, v, rxij, ryij, & |
221 |
rzij, rrf, dielectric, rtaper, natoms, & |
222 |
ux, uy, uz, mu, fx, fy, fz, tx, ty, tz, & |
223 |
ex, ey, ez, prerf ) |
224 |
endif |
225 |
end if |
226 |
endif |
227 |
endif |
228 |
|
229 |
enddo |
230 |
|
231 |
end do |
232 |
|
233 |
point(natoms) = nlist + 1 |
234 |
|
235 |
else |
236 |
|
237 |
|
238 |
!** use the list to find the neighbours ** |
239 |
|
240 |
do i = 1, natoms - 1 |
241 |
|
242 |
jbeg = point(i) |
243 |
jend = point(i+1) - 1 |
244 |
|
245 |
!** check that atom i has neighbours ** |
246 |
|
247 |
if( jbeg .le. jend ) then |
248 |
|
249 |
rxi = rx(i) |
250 |
ryi = ry(i) |
251 |
rzi = rz(i) |
252 |
|
253 |
do jnab = jbeg, jend |
254 |
|
255 |
j = list(jnab) |
256 |
|
257 |
rxij = rx(j) - rxi |
258 |
ryij = ry(j) - ryi |
259 |
rzij = rz(j) - rzi |
260 |
|
261 |
rxij = rxij - box_x * dsign( 1.0d0, rxij ) & |
262 |
* int( dabs( rxij / box_x ) + 0.5d0 ) |
263 |
ryij = ryij - box_y * dsign( 1.0d0, ryij ) & |
264 |
* int( dabs( ryij / box_y ) + 0.5d0 ) |
265 |
rzij = rzij - box_z * dsign( 1.0d0, rzij ) & |
266 |
* int( dabs( rzij / box_z ) + 0.5d0 ) |
267 |
|
268 |
rijsq = rxij * rxij + ryij * ryij + rzij * rzij |
269 |
|
270 |
if ( rijsq .lt. rcutsq ) then |
271 |
|
272 |
if ( is_vdw(i) .and. is_vdw(j) ) then |
273 |
call force_VDW( i, j, rcutsq, rijsq, sigma, epslon, & |
274 |
v, fx, fy, fz, rxij, ryij, rzij, natoms ) |
275 |
endif |
276 |
|
277 |
if ( is_lj(i) .and. is_lj(j) ) then |
278 |
call force_lj( i, j, rcutsq, rijsq, sigma, epslon, v, & |
279 |
fx, fy, fz, rxij, ryij, rzij, natoms ) |
280 |
endif |
281 |
|
282 |
if( is_ssd(i) .and. is_ssd(j) ) then |
283 |
|
284 |
rij = dsqrt(rijsq); |
285 |
|
286 |
call ssd_forces( natoms, i, j, i, j, rxij, ryij, rzij, & |
287 |
rij, v, rijsq, fx, fy, fz, tx, ty, tz, A ) |
288 |
end if |
289 |
|
290 |
endif |
291 |
|
292 |
|
293 |
if( is_dipole(i) .and. is_dipole(j) ) then |
294 |
if( rijsq .lt. rrfsq ) then |
295 |
|
296 |
call dipole_rf( i, j, rrfsq, rijsq, v, rxij, ryij, & |
297 |
rzij, rrf, dielectric, rtaper, natoms, & |
298 |
ux, uy, uz, mu, fx, fy, fz, tx, ty, tz, & |
299 |
ex, ey, ez, prerf ) |
300 |
endif |
301 |
endif |
302 |
|
303 |
end do |
304 |
endif |
305 |
end do |
306 |
endif |
307 |
|
308 |
! calculate the reaction field effect |
309 |
|
310 |
!!$ call reaction_field( v, prerf, natoms, mu, ux, uy, uz, tx, ty, tz, & |
311 |
!!$ ex, ey, ez, is_dipole ) |
312 |
|
313 |
|
314 |
end subroutine long_range_force |
315 |
|
316 |
|
317 |
|
318 |
subroutine long_range_check ( rcut, rlist, update, natoms, & |
319 |
rx,ry,rz,rx0,ry0,rz0) |
320 |
|
321 |
|
322 |
!******************************************************************* |
323 |
!** decides whether the list needs to be reconstructed. ** |
324 |
!** ** |
325 |
!** principal variables: ** |
326 |
!** ** |
327 |
!** real rx(n),ry(n),rz(n) atom positions ** |
328 |
!** real rx0(n),ry0(n),rz0(n) coordinates at last update ** |
329 |
!** real rlist radius of verlet list ** |
330 |
!** real rcut cutoff distance for forces ** |
331 |
!** real dispmx largest displacement ** |
332 |
!** integer n number of atoms ** |
333 |
!** logical update if true the list is updated ** |
334 |
!** ** |
335 |
!** usage: ** |
336 |
!** ** |
337 |
!** check is called to set update before every call to force. ** |
338 |
!******************************************************************* |
339 |
|
340 |
|
341 |
! Passed variables |
342 |
|
343 |
integer :: natoms |
344 |
|
345 |
logical(kind=2) :: update |
346 |
|
347 |
double precision rcut, rlist |
348 |
|
349 |
! Passed arrays |
350 |
double precision, dimension(natoms) :: rx, ry, rz |
351 |
double precision, dimension(natoms) :: rx0, ry0, rz0 |
352 |
|
353 |
! local variables |
354 |
|
355 |
double precision dispmx |
356 |
|
357 |
integer :: i |
358 |
|
359 |
!******************************************************************* |
360 |
|
361 |
!** calculate maximum displacement since last update ** |
362 |
|
363 |
dispmx = 0.0d0 |
364 |
|
365 |
do i = 1, natoms |
366 |
|
367 |
!write(*,*) 'calling displacement', i |
368 |
|
369 |
dispmx = max ( dabs ( rx(i) - rx0(i) ), dispmx ) |
370 |
dispmx = max ( dabs ( ry(i) - ry0(i) ), dispmx ) |
371 |
dispmx = max ( dabs ( rz(i) - rz0(i) ), dispmx ) |
372 |
|
373 |
!write(*,*) 'called displacement', i |
374 |
|
375 |
end do |
376 |
|
377 |
!** a conservative test of the list skin crossing ** |
378 |
|
379 |
dispmx = 2.0d0 * dsqrt ( 3.0d0 * dispmx ** 2 ) |
380 |
|
381 |
update = ( dispmx .gt. ( rlist - rcut ) ) |
382 |
|
383 |
|
384 |
return |
385 |
end subroutine long_range_check |
386 |
|
387 |
|
388 |
|
389 |
subroutine long_range_save(natoms,rx,ry,rz,rx0,ry0,rz0) |
390 |
|
391 |
|
392 |
!******************************************************************* |
393 |
!** saves current configuration for future checking. ** |
394 |
!** ** |
395 |
!** principal variables: ** |
396 |
!** ** |
397 |
!** real rx(n),ry(n),rz(n) atom positions ** |
398 |
!** real rx0(n),ry0(n),rz0(n) storage locations for save ** |
399 |
!** integer n number of atoms ** |
400 |
!** ** |
401 |
!** usage: ** |
402 |
!** ** |
403 |
!** save is called whenever the new verlet list is constructed. ** |
404 |
!******************************************************************* |
405 |
|
406 |
integer :: i |
407 |
|
408 |
integer :: natoms |
409 |
|
410 |
double precision, dimension(natoms) :: rx, ry, rz |
411 |
double precision, dimension(natoms) :: rx0, ry0, rz0 |
412 |
|
413 |
!******************************************************************* |
414 |
|
415 |
do i = 1, natoms |
416 |
|
417 |
rx0(i) = rx(i) |
418 |
ry0(i) = ry(i) |
419 |
rz0(i) = rz(i) |
420 |
|
421 |
end do |
422 |
|
423 |
return |
424 |
end subroutine long_range_save |
425 |
|