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 |