| 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 |
character(len=80) :: junk |
| 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, INFO=la_info) |
| 329 |
if (la_info /= 0) then |
| 330 |
write(junk,*) 'LA_GETRF error in LU-Decomposition ',la_info |
| 331 |
call error('REMOVE_ANGULAR_MOMENTUM', & |
| 332 |
junk) |
| 333 |
|
| 334 |
end if |
| 335 |
|
| 336 |
call la_getri(inertia, ipiv, INFO=la_info) |
| 337 |
if (la_info /= 0) then |
| 338 |
call error('REMOVE_ANGULAR_MOMENTUM', & |
| 339 |
'LA_GETRI error in Inverting Moment of Inertia tensor') |
| 340 |
end if |
| 341 |
|
| 342 |
inertia_inverse = inertia |
| 343 |
|
| 344 |
! calculate the angular velocities: omega = I^-1 . L |
| 345 |
|
| 346 |
omega = matmul(inertia_inverse, angmom) |
| 347 |
|
| 348 |
! subtract out center of mass velocity and angular momentum from |
| 349 |
! particle velocities |
| 350 |
|
| 351 |
do i = 1, nlocal |
| 352 |
this_q(1:3) = q_adj_com(1:3,i) |
| 353 |
v(1:3,i) = v_adj_com(1:3,i) - cross_product(omega, this_q) |
| 354 |
enddo |
| 355 |
|
| 356 |
return |
| 357 |
end subroutine remove_angular_momentum |
| 358 |
|
| 359 |
! ====================================================================================== |
| 360 |
!. functions |
| 361 |
function get_vbar(mass, this_temp) result(vbar) |
| 362 |
|
| 363 |
real( kind = DP ) :: vbar |
| 364 |
real( kind = DP ),intent(in) :: this_temp, mass |
| 365 |
real( kind = DP ), parameter :: kb = 1.381E-23_DP |
| 366 |
real( kind = DP ) :: kebar, vbar2 |
| 367 |
|
| 368 |
|
| 369 |
kebar = 3.0E0_DP * kb * this_temp / 2.0E0_DP |
| 370 |
|
| 371 |
vbar2 = 2.0E0_DP * kebar / (mass * 1.661E-27_DP) |
| 372 |
|
| 373 |
vbar = dsqrt(vbar2/3.0E0_DP) * 1.0E-5_DP |
| 374 |
|
| 375 |
return |
| 376 |
end function get_vbar |
| 377 |
|
| 378 |
subroutine calc_temp(this_temp,ke) |
| 379 |
use constants |
| 380 |
integer i, dim |
| 381 |
real( kind = DP ), intent(out) :: this_temp |
| 382 |
real( kind = DP ), intent(out), optional :: ke |
| 383 |
real( kind = DP ) :: ktmp, converter, kebar, kb |
| 384 |
real( kind = DP ) :: ke_local |
| 385 |
real( kind = DP ) :: ke_global |
| 386 |
logical :: do_ke |
| 387 |
|
| 388 |
do_ke = .false. |
| 389 |
if (present(ke)) do_ke = .true. |
| 390 |
|
| 391 |
converter = 1.0E0_DP/2.390664d3 |
| 392 |
|
| 393 |
if (do_ke) ke = 0.0E0_DP |
| 394 |
ke_local = 0.0E0_DP |
| 395 |
this_temp = 0.0E0_DP |
| 396 |
ke_global = 0.0E0_DP |
| 397 |
|
| 398 |
! write(*,*) 'ndim: ',ndim |
| 399 |
! write(*,*) 'nlocal: ',nlocal |
| 400 |
|
| 401 |
do i = 1, nlocal |
| 402 |
ktmp = 0.0E0_DP |
| 403 |
do dim = 1, ndim |
| 404 |
ktmp = ktmp + v(dim,i)**2 |
| 405 |
end do |
| 406 |
ke_local = ke_local + 0.5E0_DP * mass(i) * ktmp |
| 407 |
end do |
| 408 |
|
| 409 |
|
| 410 |
|
| 411 |
#ifdef MPI |
| 412 |
call mpi_allreduce(ke_local,ke_global,1,mpi_double_precision, & |
| 413 |
mpi_sum,mpi_comm_world,mpi_err) |
| 414 |
#else |
| 415 |
ke_global = ke_local |
| 416 |
#endif |
| 417 |
|
| 418 |
ke_global = ke_global * AMUA2FS2_TO_KCALMOL |
| 419 |
|
| 420 |
! ke_global = ke_global/converter |
| 421 |
|
| 422 |
kebar = kcal_to_kj*1000.0E0_DP * ke_global / dble(natoms) / NAVOGADRO |
| 423 |
! kebar = 4.184E3_DP * ke_global / dble(natoms) / NAVOGADRO |
| 424 |
this_temp = 2.0E0_DP * kebar / (3.0E0_DP * K_BOLTZ) |
| 425 |
|
| 426 |
if (do_ke) ke = ke_global |
| 427 |
! write(*,*) 'ke, kb, temp = ', ke, kb, temp |
| 428 |
|
| 429 |
return |
| 430 |
end subroutine calc_temp |
| 431 |
|
| 432 |
function get_com() result(com) |
| 433 |
real( kind = DP ), dimension( ndim ) :: com_local |
| 434 |
real( kind = DP ), dimension( ndim ) :: com |
| 435 |
integer :: i |
| 436 |
|
| 437 |
com = 0.0E0_DP |
| 438 |
com_local = 0.0E0_DP |
| 439 |
|
| 440 |
do i = 1, nlocal |
| 441 |
com_local(1:ndim) = com_local(1:ndim) + mass(i)*q(1:ndim,i) |
| 442 |
enddo |
| 443 |
|
| 444 |
#ifdef MPI |
| 445 |
call mpi_allreduce(com_local,com,ndim,mpi_double_precision, & |
| 446 |
mpi_sum,mpi_comm_world,mpi_err) |
| 447 |
if (mpi_err /= 0) then |
| 448 |
call error("PRINT_ANGULAR_MOMENTUM", "MPI allreduce com failed") |
| 449 |
endif |
| 450 |
#else |
| 451 |
com = com_local |
| 452 |
#endif |
| 453 |
|
| 454 |
com(1:ndim) = com(1:ndim)/total_mass |
| 455 |
|
| 456 |
|
| 457 |
return |
| 458 |
|
| 459 |
end function get_com |
| 460 |
|
| 461 |
|
| 462 |
! . ============================================================================= |
| 463 |
function get_com_v() result(com_v) |
| 464 |
!. ============================================================================= |
| 465 |
!. returns center of mass velocity |
| 466 |
|
| 467 |
real( kind = DP ), dimension( ndim ) :: com_p_local |
| 468 |
real( kind = DP ), dimension( ndim ) :: com_p |
| 469 |
real( kind = DP ), dimension( ndim ) :: com_v |
| 470 |
integer :: i |
| 471 |
|
| 472 |
com_p = 0.0E0_DP |
| 473 |
com_p_local = 0.0E0_DP |
| 474 |
|
| 475 |
do i = 1, nlocal |
| 476 |
com_p_local(1:ndim) = com_p_local(1:ndim) + mass(i)*v(1:ndim,i) |
| 477 |
enddo |
| 478 |
|
| 479 |
#ifdef MPI |
| 480 |
call mpi_allreduce(com_p_local,com_p,ndim,mpi_double_precision, & |
| 481 |
mpi_sum,mpi_comm_world,mpi_err) |
| 482 |
if (mpi_err /= 0) then |
| 483 |
call error("PRINT_ANGULAR_MOMENTUM", "MPI allreduce com failed") |
| 484 |
endif |
| 485 |
#else |
| 486 |
com_p = com_p_local |
| 487 |
#endif |
| 488 |
|
| 489 |
do i = 1, ndim |
| 490 |
com_v(i) = com_p(i)/total_mass |
| 491 |
enddo |
| 492 |
|
| 493 |
return |
| 494 |
|
| 495 |
end function get_com_v |
| 496 |
|
| 497 |
|
| 498 |
!Maybe a subroutine to get com and com_v without looping twice |
| 499 |
|
| 500 |
subroutine get_com_all(com, com_v) |
| 501 |
|
| 502 |
real( kind = DP ), dimension( ndim ) :: com_local |
| 503 |
real( kind = DP ), dimension( ndim ),intent(out) :: com |
| 504 |
real( kind = DP ), dimension( ndim ) :: com_p_local |
| 505 |
real( kind = DP ), dimension( ndim ) :: com_p |
| 506 |
real( kind = DP ), dimension( ndim ),intent(out) :: com_v |
| 507 |
integer :: i |
| 508 |
|
| 509 |
com = 0.0E0_DP |
| 510 |
com_local = 0.0E0_DP |
| 511 |
com_p = 0.0E0_DP |
| 512 |
com_p_local = 0.0E0_DP |
| 513 |
|
| 514 |
do i = 1, nlocal |
| 515 |
com_local(1:ndim) = com_local(1:ndim) + mass(i)*q(1:ndim,i) |
| 516 |
com_p_local(1:ndim) = com_p_local(1:ndim) + mass(i)*v(1:ndim,i) |
| 517 |
enddo |
| 518 |
|
| 519 |
#ifdef MPI |
| 520 |
call mpi_allreduce(com_local,com,ndim,mpi_double_precision, & |
| 521 |
mpi_sum,mpi_comm_world,mpi_err) |
| 522 |
if (mpi_err /= 0) then |
| 523 |
call error("PRINT_ANGULAR_MOMENTUM", "MPI allreduce com failed") |
| 524 |
endif |
| 525 |
|
| 526 |
call mpi_allreduce(com_p_local,com_p,ndim,mpi_double_precision, & |
| 527 |
mpi_sum,mpi_comm_world,mpi_err) |
| 528 |
if (mpi_err /= 0) then |
| 529 |
call error("PRINT_ANGULAR_MOMENTUM", "MPI allreduce com failed") |
| 530 |
endif |
| 531 |
#else |
| 532 |
com = com_local |
| 533 |
com_p = com_p_local |
| 534 |
#endif |
| 535 |
|
| 536 |
com(1:ndim) = com(1:ndim)/total_mass |
| 537 |
com_v(1:ndim) = com_p(1:ndim)/total_mass |
| 538 |
|
| 539 |
return |
| 540 |
|
| 541 |
end subroutine get_com_all |
| 542 |
|
| 543 |
end module velocity |
