1 |
chuckv |
4 |
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 |