ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/nano_mpi/src/velocity_module.F90
Revision: 4
Committed: Mon Jun 10 17:18:36 2002 UTC (22 years, 1 month ago) by chuckv
File size: 13624 byte(s)
Log Message:
Import Root

File Contents

# Content
1 module velocity
2 use parameter
3 use simulation
4 use definitions, ONLY: DP, NDIM
5 use status, ONLY: error,info,warning
6 #ifdef MPI
7 use mpi_module
8 #endif
9 implicit none
10 PRIVATE
11
12
13
14 public :: calc_temp
15 public :: quench
16 public :: scale_velocities
17 public :: randomize_velocities
18 public :: remove_cm_momentum
19 public :: remove_angular_momentum
20 public :: get_com
21 public :: get_com_all
22
23 contains
24
25 subroutine quench()
26
27 v = 0.0E0_DP
28
29 call remove_cm_momentum()
30 call remove_angular_momentum()
31
32 return
33 end subroutine quench
34
35 subroutine scale_velocities(target_t)
36
37
38 real( kind = DP ), intent(in) :: target_t
39 integer :: i, dim
40 real( kind = DP ) :: temperature
41 real( kind = DP ) :: scaling
42
43
44 call calc_temp(temperature)
45
46 scaling = dsqrt(target_t / temperature)
47
48 do i = 1, nlocal
49 do dim = 1, ndim
50 v(dim,i) = v(dim,i) * scaling
51 end do
52 end do
53
54 call remove_cm_momentum()
55 call remove_angular_momentum()
56
57 end subroutine scale_velocities
58
59 !. ===============================================================================
60 subroutine randomize_velocities(temperature)
61 use model_module
62 use constants, ONLY : twopi,pi
63 use sprng_mod, ONLY : random_gauss
64 real( kind = DP ), intent(in) :: temperature
65
66 integer :: i, dim
67 integer :: ntmp
68
69 real( kind = DP ), allocatable, dimension(:) :: vbar
70
71 real( kind = DP ) :: mtmp,vbtmp
72 real( kind = DP ) :: vrand
73
74
75
76
77 allocate(vbar(natypes))
78
79
80 do i = 1, natypes
81 if (is_atype(i)) then
82 mtmp = atype_mass(i)
83 vbar(i) = get_vbar(mtmp,temperature)
84 endif
85 enddo
86
87
88
89 ! pick random velocity from gaussian distribution
90 do i = 1, nlocal
91 vbtmp = vbar(ident(i))
92 v(1:ndim,i) = vbtmp * random_gauss()
93 end do
94
95
96
97 call remove_cm_momentum()
98 call print_angular_momentum()
99 call remove_angular_momentum()
100 call print_angular_momentum()
101
102 deallocate(vbar)
103 end subroutine randomize_velocities
104
105 !.============================================================================================
106 ! . write total angular momentum
107 subroutine print_angular_momentum()
108 use vectors, ONLY : cross_product
109 ! . variables for 2-d angular momentum
110
111 ! . varibles for 3-d angular momentum
112 real( kind = DP ), dimension( ndim ) :: ltot
113 real( kind = DP ), dimension( ndim ) :: l
114 real( kind = DP ), dimension( ndim ) :: l_local
115 real( kind = DP ), dimension( ndim ) :: r
116 real( kind = DP ), dimension( ndim ) :: p
117
118 real( kind = DP ), dimension( ndim ) :: com
119
120
121 integer :: i
122 character(len=80) :: msg
123
124 ! . for single version natoms = nlocal
125
126 ! . zero variables
127 ltot = 0.0_DP
128 l_local = 0.0_DP
129 com = 0.0_DP
130
131 ! . get center of mass
132 com = get_com()
133
134 if (ndim /= 3) return
135
136 ! . find angular momentum
137
138 ! . handle 3-d case first
139 ! . result is a vector
140
141 do i = 1, nlocal
142 r(1:ndim) = q(1:ndim,i) - com(1:ndim)
143 p(1:ndim) = v(1:ndim,i)*mass(i)
144
145 l = cross_product(r,p)
146
147 l_local(1:ndim) = l_local(1:ndim) + l(1:ndim)
148 end do
149
150 #ifdef MPI
151 call mpi_allreduce(l_local,ltot,ndim,mpi_double_precision, &
152 mpi_sum,mpi_comm_world,mpi_err)
153 if (mpi_err /= 0) then
154 call error("PRINT_ANGULAR_MOMENTUM", "MPI allreduce lx_local failed")
155 endif
156 #else
157 ltot = l_local
158 #endif
159 write(msg,'(a6,3es12.3)') 'L = ', &
160 ltot(1), ltot(2), ltot(3)
161 call info('PRINT_ANGULAR_MOMENTUM',trim(msg))
162
163 return
164
165
166 end subroutine print_angular_momentum
167
168 !. =======================================================================================
169 ! . Removes center of mass momentum
170 subroutine remove_cm_momentum()
171 integer :: i
172 real( kind = dp ), dimension( ndim ) :: com_v
173
174
175
176 ! . initialize center of mass velocity
177 com_v = 0.0E0_DP
178
179 ! . call get_com_v to get center of mass velocity
180 com_v = get_com_v()
181
182 ! . subtract center of mass momentum from velocity
183 do i = 1, nlocal
184 v(1:ndim,i) = v(1:ndim,i) - com_v(1:ndim)
185 enddo
186
187 return
188 end subroutine remove_cm_momentum
189
190
191
192
193 subroutine remove_angular_momentum()
194 use f95_lapack, only: la_getrf,la_getri
195 use vectors, only: cross_product
196 integer :: i, j, k
197 integer :: ifail,la_info
198
199 real( kind = DP ) :: xx, yy, zz
200 real( kind = DP ) :: xy, xz, yz
201
202 real( kind = DP ),dimension(3,3) :: inertia, inertia_save, identity
203 real( kind = DP ),dimension(3,3) :: inertia_inverse
204 integer ,dimension(3) :: ipiv
205 real( kind = DP), dimension(3) :: angmom, omega, angmom_local
206 real( kind = DP ),dimension(3) :: this_q, this_v
207 real( kind = DP ),dimension(9) :: inertia_vec_local
208 real( kind = DP ),dimension(9) :: inertia_vec
209 real( kind = DP ),dimension(3) :: com, com_v
210 real( kind = DP ),dimension(3,nlocal) :: q_adj_com, v_adj_com
211
212 ! . if were not in 3-d, don't even bother
213
214 if (ndim /= 3) then
215 return
216 end if
217
218 !. initialize center of mass
219 !. and velocity center of mass
220 com = 0.0E0_DP
221 com_v = 0.0E0_DP
222
223 ! com = get_com()
224 ! com_v = get_com_v()
225
226 ! The get_com_all subroutine to avoid looping twice
227
228 call get_com_all(com, com_v)
229
230 ! . initialize components for inertia tensor
231 xx=0.0_DP
232 yy=0.0_DP
233 zz=0.0_DP
234 xy=0.0_DP
235 xz=0.0_DP
236 yz=0.0_DP
237
238 ! . build components of Inertia tensor
239 !
240 ! [ Ixx -Ixy -Ixz ]
241 ! J = | -Iyx Iyy -Iyz |
242 ! [ -Izx -Iyz Izz ]
243 ! See Fowles and Cassidy Chapter 9
244 ! or Goldstein Chapter 5
245
246 !
247 do i = 1, nlocal
248 ! get components for this atom i.
249 this_q(1:3) = q(1:3,i) - com(1:3)
250 this_v(1:3) = v(1:3,i) - com_v(1:3)
251
252 q_adj_com(1:3,i) = this_q(1:3)
253 v_adj_com(1:3,i) = this_v(1:3)
254
255 ! . compute moment of inertia coefficents
256 xx = xx + &
257 this_q(1)*this_q(1)*mass(i)
258 yy = yy + &
259 this_q(2)*this_q(2)*mass(i)
260 zz = zz + &
261 this_q(3)*this_q(3)*mass(i)
262 ! . compute products of inertia
263 xy = xy + &
264 this_q(1)*this_q(2)*mass(i)
265 xz = xz + &
266 this_q(1)*this_q(3)*mass(i)
267 yz = yz + &
268 this_q(2)*this_q(3)*mass(i)
269
270 angmom_local = angmom_local + &
271 cross_product(this_q,this_v) * mass(i)
272 enddo
273
274 ! . form the local inertia vector
275 ! . f90 mpi won't allow a rank > 1 array in allreduce
276 ! . so we make a vector length 9 and transfer it back
277 ! . to a tensor after the allreduce
278 inertia_vec_local(1)=yy+zz
279 inertia_vec_local(2)=-xy
280 inertia_vec_local(3)=-xz
281 inertia_vec_local(4)=-xy
282 inertia_vec_local(5)=xx+zz
283 inertia_vec_local(6)=-yz
284 inertia_vec_local(7)=-xz
285 inertia_vec_local(8)=-yz
286 inertia_vec_local(9)=xx+yy
287
288 ! . Sum and distribute inertia and angmom arrays
289 #ifdef MPI
290 call mpi_allreduce(inertia_vec_local,inertia_vec,9,mpi_double_precision, &
291 mpi_sum,mpi_comm_world,mpi_err)
292 if (mpi_err /= 0) then
293 call error("REMOVE_ANGULAR_MOMENTUM", "MPI allreduce inertia_vec_local failed")
294 endif
295
296 call mpi_allreduce(angmom_local,angmom,3,mpi_double_precision, &
297 mpi_sum,mpi_comm_world,mpi_err)
298 if (mpi_err /= 0) then
299 call error("REMOVE_ANGULAR_MOMENTUM", "MPI allreduce angmom_local failed")
300 endif
301
302 inertia(1,1) = inertia_vec(1)
303 inertia(2,1) = inertia_vec(2)
304 inertia(3,1) = inertia_vec(3)
305 inertia(1,2) = inertia_vec(4)
306 inertia(2,2) = inertia_vec(5)
307 inertia(3,2) = inertia_vec(6)
308 inertia(1,3) = inertia_vec(7)
309 inertia(2,3) = inertia_vec(8)
310 inertia(3,3) = inertia_vec(9)
311
312 #else
313 inertia(1,1) = inertia_vec_local(1)
314 inertia(2,1) = inertia_vec_local(2)
315 inertia(3,1) = inertia_vec_local(3)
316 inertia(1,2) = inertia_vec_local(4)
317 inertia(2,2) = inertia_vec_local(5)
318 inertia(3,2) = inertia_vec_local(6)
319 inertia(1,3) = inertia_vec_local(7)
320 inertia(2,3) = inertia_vec_local(8)
321 inertia(3,3) = inertia_vec_local(9)
322
323 angmom = angmom_local
324 #endif
325
326 ! invert the moment of inertia tensor by LU-decomposition / backsolving:
327
328 call la_getrf(inertia, IPIV=ipiv, INFO=la_info)
329 if (la_info /= 0) then
330 call error('REMOVE_ANGULAR_MOMENTUM', &
331 'LA_GETRF error in LU-Decomposition')
332 end if
333
334 call la_getri(inertia, ipiv, INFO=la_info)
335 if (la_info /= 0) then
336 call error('REMOVE_ANGULAR_MOMENTUM', &
337 'LA_GETRI error in Inverting Moment of Inertia tensor')
338 end if
339
340 inertia_inverse = inertia
341
342 ! calculate the angular velocities: omega = I^-1 . L
343
344 omega = matmul(inertia_inverse, angmom)
345
346 ! subtract out center of mass velocity and angular momentum from
347 ! particle velocities
348
349 do i = 1, nlocal
350 this_q(1:3) = q_adj_com(1:3,i)
351 v(1:3,i) = v_adj_com(1:3,i) - cross_product(omega, this_q)
352 enddo
353
354 return
355 end subroutine remove_angular_momentum
356
357 ! ======================================================================================
358 !. functions
359 function get_vbar(mass, this_temp) result(vbar)
360
361 real( kind = DP ) :: vbar
362 real( kind = DP ),intent(in) :: this_temp, mass
363 real( kind = DP ), parameter :: kb = 1.381E-23_DP
364 real( kind = DP ) :: kebar, vbar2
365
366
367 kebar = 3.0E0_DP * kb * this_temp / 2.0E0_DP
368
369 vbar2 = 2.0E0_DP * kebar / (mass * 1.661E-27_DP)
370
371 vbar = dsqrt(vbar2/3.0E0_DP) * 1.0E-5_DP
372
373 return
374 end function get_vbar
375
376 subroutine calc_temp(this_temp,ke)
377 use constants
378 integer i, dim
379 real( kind = DP ), intent(out) :: this_temp
380 real( kind = DP ), intent(out), optional :: ke
381 real( kind = DP ) :: ktmp, converter, kebar, kb
382 real( kind = DP ) :: ke_local
383 real( kind = DP ) :: ke_global
384 logical :: do_ke
385
386 do_ke = .false.
387 if (present(ke)) do_ke = .true.
388
389 converter = 1.0E0_DP/2.390664d3
390
391 if (do_ke) ke = 0.0E0_DP
392 ke_local = 0.0E0_DP
393 this_temp = 0.0E0_DP
394 ke_global = 0.0E0_DP
395
396 do i = 1, nlocal
397 ktmp = 0.0E0_DP
398 do dim = 1, ndim
399 ktmp = ktmp + v(dim,i)**2
400 end do
401 ke_local = ke_local + 0.5E0_DP * mass(i) * ktmp
402 end do
403
404
405
406 #ifdef MPI
407 call mpi_allreduce(ke_local,ke_global,1,mpi_double_precision, &
408 mpi_sum,mpi_comm_world,mpi_err)
409 #else
410 ke_global = ke_local
411 #endif
412
413 ke_global = ke_global * AMUA2FS2_TO_KCALMOL
414
415 ! ke_global = ke_global/converter
416
417 kebar = kcal_to_kj*1000.0E0_DP * ke_global / dble(natoms) / NAVOGADRO
418 ! kebar = 4.184E3_DP * ke_global / dble(natoms) / NAVOGADRO
419 this_temp = 2.0E0_DP * kebar / (3.0E0_DP * K_BOLTZ)
420
421 if (do_ke) ke = ke_global
422 ! write(*,*) 'ke, kb, temp = ', ke, kb, temp
423
424 return
425 end subroutine calc_temp
426
427 function get_com() result(com)
428 real( kind = DP ), dimension( ndim ) :: com_local
429 real( kind = DP ), dimension( ndim ) :: com
430 integer :: i
431
432 com = 0.0E0_DP
433 com_local = 0.0E0_DP
434
435 do i = 1, nlocal
436 com_local(1:ndim) = com_local(1:ndim) + mass(i)*q(1:ndim,i)
437 enddo
438
439 #ifdef MPI
440 call mpi_allreduce(com_local,com,ndim,mpi_double_precision, &
441 mpi_sum,mpi_comm_world,mpi_err)
442 if (mpi_err /= 0) then
443 call error("PRINT_ANGULAR_MOMENTUM", "MPI allreduce com failed")
444 endif
445 #else
446 com = com_local
447 #endif
448
449 com(1:ndim) = com(1:ndim)/total_mass
450
451
452 return
453
454 end function get_com
455
456
457 ! . =============================================================================
458 function get_com_v() result(com_v)
459 !. =============================================================================
460 !. returns center of mass velocity
461
462 real( kind = DP ), dimension( ndim ) :: com_p_local
463 real( kind = DP ), dimension( ndim ) :: com_p
464 real( kind = DP ), dimension( ndim ) :: com_v
465 integer :: i
466
467 com_p = 0.0E0_DP
468 com_p_local = 0.0E0_DP
469
470 do i = 1, nlocal
471 com_p_local(1:ndim) = com_p_local(1:ndim) + mass(i)*v(1:ndim,i)
472 enddo
473
474 #ifdef MPI
475 call mpi_allreduce(com_p_local,com_p,ndim,mpi_double_precision, &
476 mpi_sum,mpi_comm_world,mpi_err)
477 if (mpi_err /= 0) then
478 call error("PRINT_ANGULAR_MOMENTUM", "MPI allreduce com failed")
479 endif
480 #else
481 com_p = com_p_local
482 #endif
483
484 do i = 1, ndim
485 com_v(i) = com_p(i)/total_mass
486 enddo
487
488 return
489
490 end function get_com_v
491
492
493 !Maybe a subroutine to get com and com_v without looping twice
494
495 subroutine get_com_all(com, com_v)
496
497 real( kind = DP ), dimension( ndim ) :: com_local
498 real( kind = DP ), dimension( ndim ),intent(out) :: com
499 real( kind = DP ), dimension( ndim ) :: com_p_local
500 real( kind = DP ), dimension( ndim ) :: com_p
501 real( kind = DP ), dimension( ndim ),intent(out) :: com_v
502 integer :: i
503
504 com = 0.0E0_DP
505 com_local = 0.0E0_DP
506 com_p = 0.0E0_DP
507 com_p_local = 0.0E0_DP
508
509 do i = 1, nlocal
510 com_local(1:ndim) = com_local(1:ndim) + mass(i)*q(1:ndim,i)
511 com_p_local(1:ndim) = com_p_local(1:ndim) + mass(i)*v(1:ndim,i)
512 enddo
513
514 #ifdef MPI
515 call mpi_allreduce(com_local,com,ndim,mpi_double_precision, &
516 mpi_sum,mpi_comm_world,mpi_err)
517 if (mpi_err /= 0) then
518 call error("PRINT_ANGULAR_MOMENTUM", "MPI allreduce com failed")
519 endif
520
521 call mpi_allreduce(com_p_local,com_p,ndim,mpi_double_precision, &
522 mpi_sum,mpi_comm_world,mpi_err)
523 if (mpi_err /= 0) then
524 call error("PRINT_ANGULAR_MOMENTUM", "MPI allreduce com failed")
525 endif
526 #else
527 com = com_local
528 com_p = com_p_local
529 #endif
530
531 com(1:ndim) = com(1:ndim)/total_mass
532 com_v(1:ndim) = com_p(1:ndim)/total_mass
533
534 return
535
536 end subroutine get_com_all
537
538 end module velocity