1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Independent helium subroutines shared with other modules 8!> \author Lukasz Walewski 9!> \date 2009-07-14 10!> \note Avoiding circular deps: do not USE any other helium_* modules here. 11! ************************************************************************************************** 12MODULE helium_common 13 14 USE helium_types, ONLY: helium_solvent_type,& 15 int_arr_ptr,& 16 rho_atom_number,& 17 rho_moment_of_inertia,& 18 rho_num,& 19 rho_projected_area,& 20 rho_winding_cycle,& 21 rho_winding_number 22 USE input_constants, ONLY: denominator_inertia,& 23 denominator_natoms,& 24 denominator_rperp2,& 25 helium_cell_shape_octahedron 26 USE kinds, ONLY: default_string_length,& 27 dp 28 USE mathconstants, ONLY: pi 29 USE memory_utilities, ONLY: reallocate 30 USE parallel_rng_types, ONLY: next_random_number 31 USE physcon, ONLY: angstrom,& 32 bohr 33 USE pint_public, ONLY: pint_com_pos 34 USE pint_types, ONLY: pint_env_type 35 USE splines_methods, ONLY: spline_value 36 USE splines_types, ONLY: spline_data_type 37#include "../base/base_uses.f90" 38 39 IMPLICIT NONE 40 41 PRIVATE 42 43 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 44 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_common' 45 46 PUBLIC :: helium_pbc 47 PUBLIC :: helium_boxmean_3d 48 PUBLIC :: helium_calc_rdf 49 PUBLIC :: helium_calc_rho 50 PUBLIC :: helium_calc_plength 51 PUBLIC :: helium_rotate 52 PUBLIC :: helium_eval_expansion 53 PUBLIC :: helium_eval_chain 54 PUBLIC :: helium_update_transition_matrix 55 PUBLIC :: helium_update_coord_system 56 PUBLIC :: helium_spline 57 PUBLIC :: helium_cycle_number 58 PUBLIC :: helium_path_length 59 PUBLIC :: helium_is_winding 60 PUBLIC :: helium_cycle_of 61 PUBLIC :: helium_total_winding_number 62 PUBLIC :: helium_total_projected_area 63 PUBLIC :: helium_total_moment_of_inertia 64 PUBLIC :: helium_com 65 PUBLIC :: helium_set_rdf_coord_system 66 67CONTAINS 68 69! *************************************************************************** 70!> \brief General PBC routine for helium. 71!> \param helium ... 72!> \param r ... 73!> \param enforce ... 74!> \date 2019-12-09 75!> \author Harald Forbert 76!> \note Apply periodic boundary conditions as needed 77!> \note Combines several older routines into a single simpler/faster routine 78! ************************************************************************************************** 79 SUBROUTINE helium_pbc(helium, r, enforce) 80 81 TYPE(helium_solvent_type), POINTER :: helium 82 REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: r 83 LOGICAL, OPTIONAL :: enforce 84 85 REAL(KIND=dp) :: cell_size, cell_size_inv, corr, rx, ry, & 86 rz, sx, sy, sz 87 88 IF (.NOT. (helium%periodic .OR. PRESENT(enforce))) RETURN 89 90 cell_size = helium%cell_size 91 cell_size_inv = helium%cell_size_inv 92 93 rx = r(1)*cell_size_inv 94 IF (rx > 0.5_dp) THEN 95 rx = rx - REAL(INT(rx + 0.5_dp), dp) 96 ELSEIF (rx < -0.5_dp) THEN 97 rx = rx - REAL(INT(rx - 0.5_dp), dp) 98 END IF 99 100 ry = r(2)*cell_size_inv 101 IF (ry > 0.5_dp) THEN 102 ry = ry - REAL(INT(ry + 0.5_dp), dp) 103 ELSEIF (ry < -0.5_dp) THEN 104 ry = ry - REAL(INT(ry - 0.5_dp), dp) 105 END IF 106 107 rz = r(3)*cell_size_inv 108 IF (rz > 0.5_dp) THEN 109 rz = rz - REAL(INT(rz + 0.5_dp), dp) 110 ELSEIF (rz < -0.5_dp) THEN 111 rz = rz - REAL(INT(rz - 0.5_dp), dp) 112 END IF 113 114 IF (helium%cell_shape == helium_cell_shape_octahedron) THEN 115 corr = 0.0_dp 116 IF (rx > 0.0_dp) THEN 117 corr = corr + rx 118 sx = 0.5_dp 119 ELSE 120 corr = corr - rx 121 sx = -0.5_dp 122 END IF 123 IF (ry > 0.0_dp) THEN 124 corr = corr + ry 125 sy = 0.5_dp 126 ELSE 127 corr = corr - ry 128 sy = -0.5_dp 129 END IF 130 IF (rz > 0.0_dp) THEN 131 corr = corr + rz 132 sz = 0.5_dp 133 ELSE 134 corr = corr - rz 135 sz = -0.5_dp 136 END IF 137 IF (corr > 0.75_dp) THEN 138 rx = rx - sx 139 ry = ry - sy 140 rz = rz - sz 141 END IF 142 END IF 143 144 r(1) = rx*cell_size 145 r(2) = ry*cell_size 146 r(3) = rz*cell_size 147 148 RETURN 149 END SUBROUTINE helium_pbc 150 151! *************************************************************************** 152!> \brief find distance squared of nearest image 153!> \param helium ... 154!> \param r ... 155!> \return ... 156!> \date 2019-12-09 157!> \author Harald Forbert 158!> \note Given a distance vector r, return the distance squared 159!> of the nearest image. Cell information is taken from the 160!> helium environment argument. 161! ************************************************************************************************** 162 163 PURE FUNCTION helium_pbc_r2(helium, r) 164 165 TYPE(helium_solvent_type), POINTER :: helium 166 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r 167 REAL(KIND=dp) :: helium_pbc_r2 168 169 REAL(KIND=dp) :: cell_size_inv, corr, rx, ry, rz 170 171 IF (helium%periodic) THEN 172 cell_size_inv = helium%cell_size_inv 173 rx = ABS(r(1)*cell_size_inv) 174 rx = ABS(rx - REAL(INT(rx + 0.5_dp), dp)) 175 ry = ABS(r(2)*cell_size_inv) 176 ry = ABS(ry - REAL(INT(ry + 0.5_dp), dp)) 177 rz = ABS(r(3)*cell_size_inv) 178 rz = ABS(rz - REAL(INT(rz + 0.5_dp), dp)) 179 IF (helium%cell_shape == helium_cell_shape_octahedron) THEN 180 corr = rx + ry + rz - 0.75_dp 181 IF (corr < 0.0_dp) corr = 0.0_dp 182 helium_pbc_r2 = helium%cell_size**2*(rx*rx + ry*ry + rz*rz - corr) 183 ELSE 184 helium_pbc_r2 = helium%cell_size**2*(rx*rx + ry*ry + rz*rz) 185 END IF 186 ELSE 187 helium_pbc_r2 = (r(1)*r(1) + r(2)*r(2) + r(3)*r(3)) 188 END IF 189 RETURN 190 END FUNCTION helium_pbc_r2 191 192! *************************************************************************** 193!> \brief Calculate the point equidistant from two other points a and b 194!> within the helium box - 3D version 195!> \param helium helium environment for which 196!> \param a vectors for which to find the mean within the He box 197!> \param b vectors for which to find the mean within the He box 198!> \param c ... 199!> \par History 200!> 2009-10-02 renamed, originally was helium_boxmean [lwalewski] 201!> 2009-10-02 redesigned so it is now called as a subroutine [lwalewski] 202!> 2009-10-02 redesigned so it now gets/returns a 3D vectors [lwalewski] 203!> \author hforbert 204! ************************************************************************************************** 205 SUBROUTINE helium_boxmean_3d(helium, a, b, c) 206 207 TYPE(helium_solvent_type), POINTER :: helium 208 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a, b 209 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: c 210 211 c(:) = b(:) - a(:) 212 CALL helium_pbc(helium, c) 213 c(:) = a(:) + 0.5_dp*c(:) 214 CALL helium_pbc(helium, c) 215 RETURN 216 END SUBROUTINE helium_boxmean_3d 217 218! *************************************************************************** 219!> \brief Given the permutation state assign cycle lengths to all He atoms. 220!> \param helium ... 221!> \date 2011-07-06 222!> \author Lukasz Walewski 223!> \note The helium%atom_plength array is filled with cycle lengths, 224!> each atom gets the length of the permutation cycle it belongs to. 225! ************************************************************************************************** 226 SUBROUTINE helium_calc_atom_cycle_length(helium) 227 228 TYPE(helium_solvent_type), POINTER :: helium 229 230 CHARACTER(len=*), PARAMETER :: routineN = 'helium_calc_atom_cycle_length', & 231 routineP = moduleN//':'//routineN 232 233 CHARACTER(len=default_string_length) :: err_str 234 INTEGER :: clen, curr_idx, handle, ia, start_idx 235 INTEGER, DIMENSION(:), POINTER :: atoms_in_cycle 236 LOGICAL :: atoms_left, path_end_reached 237 LOGICAL, DIMENSION(:), POINTER :: atom_was_used 238 239 CALL timeset(routineN, handle) 240 241 NULLIFY (atoms_in_cycle) 242 atoms_in_cycle => helium%itmp_atoms_1d 243 atoms_in_cycle(:) = 0 244 245 NULLIFY (atom_was_used) 246 atom_was_used => helium%ltmp_atoms_1d 247 atom_was_used(:) = .FALSE. 248 249 helium%atom_plength(:) = 0 250 251 start_idx = 1 252 DO 253 clen = 0 254 path_end_reached = .FALSE. 255 curr_idx = start_idx 256 DO ia = 1, helium%atoms 257 clen = clen + 1 258 atoms_in_cycle(clen) = curr_idx 259 atom_was_used(curr_idx) = .TRUE. 260 261 ! follow to the next atom in the cycle 262 curr_idx = helium%permutation(curr_idx) 263 IF (curr_idx .EQ. start_idx) THEN 264 path_end_reached = .TRUE. 265 EXIT 266 END IF 267 END DO 268 err_str = "Permutation path is not cyclic." 269 IF (.NOT. path_end_reached) THEN 270 CPABORT(err_str) 271 END IF 272 273 ! assign the cycle length to the collected atoms 274 DO ia = 1, clen 275 helium%atom_plength(atoms_in_cycle(ia)) = clen 276 END DO 277 278 ! go to the next "unused" atom 279 atoms_left = .FALSE. 280 DO ia = 1, helium%atoms 281 IF (.NOT. atom_was_used(ia)) THEN 282 start_idx = ia 283 atoms_left = .TRUE. 284 EXIT 285 END IF 286 END DO 287 288 IF (.NOT. atoms_left) EXIT 289 END DO 290 291 atoms_in_cycle(:) = 0 292 NULLIFY (atoms_in_cycle) 293 atom_was_used(:) = .FALSE. 294 NULLIFY (atom_was_used) 295 296 CALL timestop(handle) 297 298 RETURN 299 END SUBROUTINE helium_calc_atom_cycle_length 300 301! *************************************************************************** 302!> \brief Decompose a permutation into cycles 303!> \param permutation ... 304!> \return ... 305!> \date 2013-12-11 306!> \author Lukasz Walewski 307!> \note Given a permutation return a pointer to an array of pointers, 308!> with each element pointing to a cycle of this permutation. 309!> \note This function allocates memory and returns a pointer to it, 310!> do not forget to deallocate once finished with the results. 311! ************************************************************************************************** 312 FUNCTION helium_calc_cycles(permutation) RESULT(cycles) 313 314 INTEGER, DIMENSION(:), POINTER :: permutation 315 TYPE(int_arr_ptr), DIMENSION(:), POINTER :: cycles 316 317 CHARACTER(len=*), PARAMETER :: routineN = 'helium_calc_cycles', & 318 routineP = moduleN//':'//routineN 319 320 INTEGER :: curat, ncycl, nsize, nused 321 INTEGER, DIMENSION(:), POINTER :: cur_cycle, used_indices 322 TYPE(int_arr_ptr), DIMENSION(:), POINTER :: my_cycles 323 324 nsize = SIZE(permutation) 325 326 ! most pesimistic scenario: no cycles longer than 1 327 ALLOCATE (my_cycles(nsize)) 328 329 ! go over the permutation and identify cycles 330 curat = 1 331 nused = 0 332 ncycl = 0 333 DO WHILE (curat .LE. nsize) 334 335 ! get the permutation cycle the current atom belongs to 336 cur_cycle => helium_cycle_of(curat, permutation) 337 338 ! include the current cycle in the pool of "used" indices 339 nused = nused + SIZE(cur_cycle) 340 CALL reallocate(used_indices, 1, nused) 341 used_indices(nused - SIZE(cur_cycle) + 1:nused) = cur_cycle(1:SIZE(cur_cycle)) 342 343 ! store the pointer to the current cycle 344 ncycl = ncycl + 1 345 my_cycles(ncycl)%iap => cur_cycle 346 347 ! free the local pointer 348 NULLIFY (cur_cycle) 349 350 ! try to increment the current atom index 351 DO WHILE (ANY(used_indices .EQ. curat)) 352 curat = curat + 1 353 END DO 354 355 END DO 356 357 DEALLOCATE (used_indices) 358 NULLIFY (used_indices) 359 360 ! assign the result 361 ALLOCATE (cycles(ncycl)) 362 cycles(1:ncycl) = my_cycles(1:ncycl) 363 364 DEALLOCATE (my_cycles) 365 NULLIFY (my_cycles) 366 367 RETURN 368 END FUNCTION helium_calc_cycles 369 370! *************************************************************************** 371!> \brief Calculate helium density distribution functions and store them 372!> in helium%rho_inst 373!> \param helium ... 374!> \date 2011-06-14 375!> \author Lukasz Walewski 376! ************************************************************************************************** 377 SUBROUTINE helium_calc_rho(helium) 378 379 TYPE(helium_solvent_type), POINTER :: helium 380 381 CHARACTER(len=*), PARAMETER :: routineN = 'helium_calc_rho', & 382 routineP = moduleN//':'//routineN 383 384 CHARACTER(len=default_string_length) :: msgstr 385 INTEGER :: aa, bx, by, bz, handle, ia, ib, ic, id, & 386 ii, ir, n_out_of_range, nbin 387 INTEGER, DIMENSION(3) :: nw 388 INTEGER, DIMENSION(:), POINTER :: cycl_len 389 LOGICAL :: ltmp1, ltmp2, ltmp3 390 REAL(KIND=dp) :: invd3r, invdr, invp, rtmp 391 REAL(KIND=dp), DIMENSION(3) :: maxr_half, r, vlink, vtotal, wn 392 TYPE(int_arr_ptr), DIMENSION(:), POINTER :: cycles 393 394 CALL timeset(routineN, handle) 395 396 maxr_half(:) = helium%rho_maxr/2.0_dp 397 invdr = 1.0_dp/helium%rho_delr 398 invd3r = invdr*invdr*invdr 399 invp = 1.0_dp/REAL(helium%beads, dp) 400 nbin = helium%rho_nbin 401 ! assign the cycle lengths to all the atoms 402 CALL helium_calc_atom_cycle_length(helium) 403 NULLIFY (cycl_len) 404 cycl_len => helium%atom_plength 405 DO ir = 1, rho_num ! loop over densities --- 406 407 IF (.NOT. helium%rho_property(ir)%is_calculated) THEN 408 ! skip densities that are not requested by the user 409 CYCLE 410 END IF 411 412 SELECT CASE (ir) 413 414 CASE (rho_atom_number) 415 ii = helium%rho_property(ir)%component_index(1) 416 helium%rho_incr(ii, :, :) = invp 417 418 CASE (rho_projected_area) 419 vtotal(:) = helium_total_projected_area(helium) 420 DO ia = 1, helium%atoms 421 DO ib = 1, helium%beads 422 vlink(:) = helium_link_projected_area(helium, ia, ib) 423 DO ic = 1, 3 424 ii = helium%rho_property(ir)%component_index(ic) 425 helium%rho_incr(ii, ia, ib) = vtotal(ic)*vlink(ic)*angstrom*angstrom*angstrom*angstrom 426 END DO 427 END DO 428 END DO 429 430! CASE (rho_winding_number) 431! vtotal(:) = helium_total_winding_number(helium) 432! DO ia = 1, helium%atoms 433! DO ib = 1, helium%beads 434! vlink(:) = helium_link_winding_number(helium,ia,ib) 435! DO ic = 1, 3 436! ii = helium%rho_property(ir)%component_index(ic) 437! helium%rho_incr(ii,ia,ib) = vtotal(ic)*vlink(ic)*angstrom*angstrom 438! END DO 439! END DO 440! END DO 441 442 CASE (rho_winding_number) 443 vtotal(:) = helium_total_winding_number(helium) 444 DO id = 1, 3 445 ii = helium%rho_property(ir)%component_index(id) 446 helium%rho_incr(ii, :, :) = 0.0_dp 447 END DO 448 NULLIFY (cycles) 449 cycles => helium_calc_cycles(helium%permutation) 450 DO ic = 1, SIZE(cycles) 451 wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos) 452 DO ia = 1, SIZE(cycles(ic)%iap) 453 aa = cycles(ic)%iap(ia) 454 DO ib = 1, helium%beads 455 vlink(:) = helium_link_winding_number(helium, aa, ib) 456 DO id = 1, 3 457 IF (ABS(wn(id)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN 458 ii = helium%rho_property(ir)%component_index(id) 459 helium%rho_incr(ii, aa, ib) = vtotal(id)*vlink(id)*angstrom*angstrom 460 END IF 461 END DO 462 END DO 463 END DO 464 END DO 465 DEALLOCATE (cycles) 466 467 CASE (rho_winding_cycle) 468 vtotal(:) = helium_total_winding_number(helium) 469 DO id = 1, 3 470 ii = helium%rho_property(ir)%component_index(id) 471 helium%rho_incr(ii, :, :) = 0.0_dp 472 END DO 473 NULLIFY (cycles) 474 cycles => helium_calc_cycles(helium%permutation) 475 ! compute number of atoms in all winding cycles 476 nw(:) = 0 477 DO ic = 1, SIZE(cycles) 478 wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos) 479 DO id = 1, 3 480 IF (ABS(wn(id)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN 481 nw(id) = nw(id) + SIZE(cycles(ic)%iap) 482 END IF 483 END DO 484 END DO 485 ! assign contributions to all beads of winding cycles only 486 DO ic = 1, SIZE(cycles) 487 wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos) 488 DO id = 1, 3 489 IF (ABS(wn(id)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN 490 DO ia = 1, SIZE(cycles(ic)%iap) 491 aa = cycles(ic)%iap(ia) 492 DO ib = 1, helium%beads 493 IF (nw(id) .GT. 0) THEN ! this test should always get passed 494 ii = helium%rho_property(ir)%component_index(id) 495 rtmp = invp/REAL(nw(id), dp) 496 helium%rho_incr(ii, aa, ib) = rtmp*vtotal(id)*vtotal(id)*angstrom*angstrom 497 END IF 498 END DO 499 END DO 500 END IF 501 END DO 502 END DO 503 DEALLOCATE (cycles) 504 505 CASE (rho_moment_of_inertia) 506 vtotal(:) = helium_total_moment_of_inertia(helium) 507 DO ia = 1, helium%atoms 508 DO ib = 1, helium%beads 509 vlink(:) = helium_link_moment_of_inertia(helium, ia, ib) 510 DO ic = 1, 3 511 ii = helium%rho_property(ir)%component_index(ic) 512 helium%rho_incr(ii, ia, ib) = vlink(ic)*angstrom*angstrom 513 END DO 514 END DO 515 END DO 516 517 CASE DEFAULT 518 ! do nothing 519 520 END SELECT 521 522 END DO ! loop over densities --- 523 524 n_out_of_range = 0 525 helium%rho_inst(:, :, :, :) = 0.0_dp 526 DO ia = 1, helium%atoms 527 ! bin the bead positions of the current atom using the increments set above 528 DO ib = 1, helium%beads 529 ! map the current bead position to the corresponding voxel 530 r(:) = helium%pos(:, ia, ib) - helium%center(:) 531 ! enforce PBC even if this is a non-periodic calc to avoid density leakage 532 CALL helium_pbc(helium, r, enforce=.TRUE.) 533 ! set up bin indices (translate by L/2 to avoid non-positive array indices) 534 bx = INT((r(1) + maxr_half(1))*invdr) + 1 535 by = INT((r(2) + maxr_half(2))*invdr) + 1 536 bz = INT((r(3) + maxr_half(3))*invdr) + 1 537 ! check that the resulting bin numbers are within array bounds 538 ltmp1 = (0 .LT. bx) .AND. (bx .LE. nbin) 539 ltmp2 = (0 .LT. by) .AND. (by .LE. nbin) 540 ltmp3 = (0 .LT. bz) .AND. (bz .LE. nbin) 541 IF (ltmp1 .AND. ltmp2 .AND. ltmp3) THEN 542 ! increment all the estimators (those that the current atom does not 543 ! contribute to have increment incr(ic)==0) 544 DO ic = 1, helium%rho_num_act 545 helium%rho_inst(ic, bx, by, bz) = helium%rho_inst(ic, bx, by, bz) + helium%rho_incr(ic, ia, ib) 546 END DO 547 ELSE 548 n_out_of_range = n_out_of_range + 1 549 END IF 550 END DO 551 END DO 552 ! scale by volume element 553 helium%rho_inst(:, :, :, :) = helium%rho_inst(:, :, :, :)*invd3r 554 555 ! stop if any bead fell out of the range 556 ! since enforced PBC should have taken care of such leaks 557 WRITE (msgstr, *) n_out_of_range 558 msgstr = "Number of bead positions out of range: "//TRIM(ADJUSTL(msgstr)) 559 IF (n_out_of_range .GT. 0) THEN 560 CPABORT(msgstr) 561 END IF 562 563 CALL timestop(handle) 564 565 RETURN 566 END SUBROUTINE helium_calc_rho 567 568#if 0 569! *************************************************************************** 570!> \brief Normalize superfluid densities according to the input keyword 571!> HELIUM%SUPERFLUID_ESTIMATOR%DENOMINATOR 572!> \param helium ... 573!> \param rho ... 574!> \date 2014-06-24 575!> \author Lukasz Walewski 576! ************************************************************************************************** 577 SUBROUTINE helium_norm_rho(helium, rho) 578 579 TYPE(helium_solvent_type), POINTER :: helium 580 REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: rho 581 582 CHARACTER(len=*), PARAMETER :: routineN = 'helium_norm_rho', & 583 routineP = moduleN//':'//routineN 584 585 INTEGER :: ix, iy, iz, ndim 586 REAL(KIND=dp) :: invatoms, rx, ry, rz 587 REAL(KIND=dp), DIMENSION(3) :: invmoit, invrperp, ro 588 589 SELECT CASE (helium%supest_denominator) 590 591 CASE (denominator_natoms) 592 invatoms = 1.0_dp/REAL(helium%atoms, dp) 593 rho(2, :, :, :) = rho(2, :, :, :)*invatoms 594 rho(3, :, :, :) = rho(3, :, :, :)*invatoms 595 rho(4, :, :, :) = rho(4, :, :, :)*invatoms 596 597 CASE (denominator_inertia) 598 invmoit(:) = REAL(helium%atoms, dp)/helium%mominer%ravr(:) 599 rho(2, :, :, :) = rho(2, :, :, :)*invmoit(1) 600 rho(3, :, :, :) = rho(3, :, :, :)*invmoit(2) 601 rho(4, :, :, :) = rho(4, :, :, :)*invmoit(3) 602 603 CASE (denominator_rperp2) 604 ndim = helium%rho_nbin 605 ro(:) = helium%center(:) - 0.5_dp*(helium%rho_maxr - helium%rho_delr) 606 DO ix = 1, ndim 607 DO iy = 1, ndim 608 DO iz = 1, ndim 609 rx = ro(1) + REAL(ix - 1, dp)*helium%rho_delr 610 ry = ro(2) + REAL(iy - 1, dp)*helium%rho_delr 611 rz = ro(3) + REAL(iz - 1, dp)*helium%rho_delr 612 invrperp(1) = 1.0_dp/(ry*ry + rz*rz) 613 invrperp(2) = 1.0_dp/(rz*rz + rx*rx) 614 invrperp(3) = 1.0_dp/(rx*rx + ry*ry) 615 rho(2, ix, iy, iz) = rho(2, ix, iy, iz)*invrperp(1) 616 rho(3, ix, iy, iz) = rho(3, ix, iy, iz)*invrperp(2) 617 rho(4, ix, iy, iz) = rho(4, ix, iy, iz)*invrperp(3) 618 END DO 619 END DO 620 END DO 621 622 CASE DEFAULT 623 ! do nothing 624 625 END SELECT 626 627 RETURN 628 END SUBROUTINE helium_norm_rho 629#endif 630 631! *************************************************************************** 632!> \brief Calculate helium radial distribution functions wrt positions given 633!> by <centers> using the atomic weights given by <weights>. 634!> \param helium ... 635!> \date 2009-07-22 636!> \par History 637!> 2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran] 638!> \author Lukasz Walewski 639! ************************************************************************************************** 640 SUBROUTINE helium_calc_rdf(helium) 641 642 TYPE(helium_solvent_type), POINTER :: helium 643 644 CHARACTER(len=*), PARAMETER :: routineN = 'helium_calc_rdf', & 645 routineP = moduleN//':'//routineN 646 647 CHARACTER(len=default_string_length) :: msgstr 648 INTEGER :: bin, handle, ia, ib, ic, ind_hehe, & 649 n_out_of_range, nbin 650 REAL(KIND=dp) :: invdr, nideal, ri, rlower, rupper 651 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pref 652 REAL(KIND=dp), DIMENSION(3) :: r, r0 653 654 CALL timeset(routineN, handle) 655 656 invdr = 1.0_dp/helium%rdf_delr 657 nbin = helium%rdf_nbin 658 n_out_of_range = 0 659 helium%rdf_inst(:, :) = 0.0_dp 660 661 ind_hehe = 0 662 ! Calculate He-He RDF 663 IF (helium%rdf_he_he) THEN 664 ind_hehe = 1 665 DO ib = 1, helium%beads 666 DO ia = 1, helium%atoms 667 r0(:) = helium%pos(:, ia, ib) 668 669 DO ic = 1, helium%atoms 670 IF (ia == ic) CYCLE 671 r(:) = helium%pos(:, ic, ib) - r0(:) 672 CALL helium_pbc(helium, r) 673 ri = SQRT(r(1)*r(1) + r(2)*r(2) + r(3)*r(3)) 674 bin = INT(ri*invdr) + 1 675 IF ((0 .LT. bin) .AND. (bin .LE. nbin)) THEN 676 ! increment the RDF value for He atoms inside the r_6 sphere 677 helium%rdf_inst(ind_hehe, bin) = helium%rdf_inst(ind_hehe, bin) + 1.0_dp 678 ELSE 679 n_out_of_range = n_out_of_range + 1 680 END IF 681 END DO 682 END DO 683 END DO 684 END IF 685 686 ! Calculate Solute-He RDF 687 IF (helium%solute_present .AND. helium%rdf_sol_he) THEN 688 DO ib = 1, helium%beads 689 DO ia = 1, helium%solute_atoms 690 r0(:) = helium%rdf_centers(ib, 3*(ia - 1) + 1:3*(ia - 1) + 3) 691 692 DO ic = 1, helium%atoms 693 r(:) = helium%pos(:, ic, ib) - r0(:) 694 CALL helium_pbc(helium, r) 695 ri = SQRT(r(1)*r(1) + r(2)*r(2) + r(3)*r(3)) 696 bin = INT(ri*invdr) + 1 697 IF ((0 .LT. bin) .AND. (bin .LE. nbin)) THEN 698 ! increment the RDF value for He atoms inside the r_6 sphere 699 helium%rdf_inst(ind_hehe + ia, bin) = helium%rdf_inst(ind_hehe + ia, bin) + 1.0_dp 700 ELSE 701 n_out_of_range = n_out_of_range + 1 702 END IF 703 END DO 704 END DO 705 END DO 706 707 END IF 708 709 ! for periodic case we intentionally truncate RDFs to spherical volume 710 ! so we skip atoms in the corners 711 IF (.NOT. helium%periodic) THEN 712 IF (n_out_of_range .GT. 0) THEN 713 WRITE (msgstr, *) n_out_of_range 714 msgstr = "Number of bead positions out of range: "//TRIM(ADJUSTL(msgstr)) 715 CPWARN(msgstr) 716 END IF 717 END IF 718 719 ! normalize the histogram to get g(r) 720 ALLOCATE (pref(helium%rdf_num)) 721 pref(:) = 0.0_dp 722 IF (helium%periodic) THEN 723 ! Use helium density for normalization 724 pref(:) = helium%density*helium%beads*helium%atoms 725 ! Correct for He-He-RDF 726 IF (helium%rdf_he_he) THEN 727 pref(1) = pref(1)/helium%atoms*(helium%atoms - 1) 728 END IF 729 ELSE 730 ! Non-periodic case has density of 0, use integral for normalzation 731 ! This leads to a unit of 1/volume of the RDF 732 pref(:) = 0.5_dp*helium%rdf_inst(:, 1) 733 DO bin = 2, helium%rdf_nbin - 1 734 pref(:) = pref(:) + helium%rdf_inst(:, bin) 735 END DO 736 pref(:) = pref(:) + 0.5_dp*helium%rdf_inst(:, helium%rdf_nbin) 737 738 !set integral of histogram to number of atoms: 739 pref(:) = pref(:)/helium%atoms 740 ! Correct for He-He-RDF 741 IF (helium%rdf_he_he) THEN 742 pref(1) = pref(1)*helium%atoms/(helium%atoms - 1) 743 END IF 744 END IF 745 ! Volume integral first: 746 DO bin = 1, helium%rdf_nbin 747 rlower = REAL(bin - 1, dp)*helium%rdf_delr 748 rupper = rlower + helium%rdf_delr 749 nideal = (rupper**3 - rlower**3)*4.0_dp*pi/3.0_dp 750 helium%rdf_inst(:, bin) = helium%rdf_inst(:, bin)/nideal 751 END DO 752 ! No normalization for density 753 pref(:) = 1.0_dp/pref(:) 754 DO ia = 1, helium%rdf_num 755 helium%rdf_inst(ia, :) = helium%rdf_inst(ia, :)*pref(ia) 756 END DO 757 758 DEALLOCATE (pref) 759 760 CALL timestop(handle) 761 762 RETURN 763 END SUBROUTINE helium_calc_rdf 764 765! *************************************************************************** 766!> \brief Calculate probability distribution of the permutation lengths 767!> \param helium ... 768!> \date 2010-06-07 769!> \author Lukasz Walewski 770!> \note Valid permutation path length is an integer (1, NATOMS), number 771!> of paths of a given length is calculated here and average over 772!> inner loop iterations and helium environments is done in 773!> helium_sample. 774! ************************************************************************************************** 775 SUBROUTINE helium_calc_plength(helium) 776 777 TYPE(helium_solvent_type), POINTER :: helium 778 779 CHARACTER(len=*), PARAMETER :: routineN = 'helium_calc_plength', & 780 routineP = moduleN//':'//routineN 781 782 INTEGER :: i, j, k 783 784 helium%plength_inst(:) = 0.0_dp 785 DO i = 1, helium%atoms 786 j = helium%permutation(i) 787 k = 1 788 DO 789 IF (j == i) EXIT 790 k = k + 1 791 j = helium%permutation(j) 792 END DO 793 helium%plength_inst(k) = helium%plength_inst(k) + 1 794 END DO 795 helium%plength_inst(:) = helium%plength_inst(:)/helium%atoms 796 797 RETURN 798 END SUBROUTINE helium_calc_plength 799 800! *************************************************************************** 801!> \brief Rotate helium particles in imaginary time by nslices 802!> \param helium ... 803!> \param nslices ... 804!> \author hforbert 805!> \note Positions of helium beads in helium%pos array are reorganized such 806!> that the indices are cyclically translated in a permutation-aware 807!> manner. helium%relrot is given a new value that represents the new 808!> 'angle' of the beads. This is done modulo helium%beads, so relrot 809!> should be always within 0 (no rotation) and helium%beads-1 (almost 810!> full rotation). [lwalewski] 811! ************************************************************************************************** 812 SUBROUTINE helium_rotate(helium, nslices) 813 TYPE(helium_solvent_type), POINTER :: helium 814 INTEGER, INTENT(IN) :: nslices 815 816 CHARACTER(len=*), PARAMETER :: routineN = 'helium_rotate', routineP = moduleN//':'//routineN 817 818 INTEGER :: b, i, j, k, n 819 820 b = helium%beads 821 n = helium%atoms 822 i = MOD(nslices, b) 823 IF (i < 0) i = i + b 824 IF ((i >= b) .OR. (i < 1)) RETURN 825 helium%relrot = MOD(helium%relrot + i, b) 826 DO k = 1, i 827 helium%work(:, :, k) = helium%pos(:, :, k) 828 END DO 829 DO k = i + 1, b 830 helium%pos(:, :, k - i) = helium%pos(:, :, k) 831 END DO 832 DO k = 1, i 833 DO j = 1, n 834 helium%pos(:, j, b - i + k) = helium%work(:, helium%permutation(j), k) 835 END DO 836 END DO 837 RETURN 838 END SUBROUTINE helium_rotate 839 840! *************************************************************************** 841!> \brief Calculate the pair-product action or energy by evaluating the 842!> power series expansion according to Eq. 4.46 in Ceperley 1995. 843!> \param helium ... 844!> \param r ... 845!> \param rp ... 846!> \param work ... 847!> \param action ... 848!> \return ... 849!> \author Harald Forbert 850! ************************************************************************************************** 851 FUNCTION helium_eval_expansion(helium, r, rp, work, action) RESULT(res) 852 853 TYPE(helium_solvent_type), POINTER :: helium 854 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r, rp 855 REAL(KIND=dp), CONTIGUOUS, DIMENSION(:), & 856 INTENT(INOUT) :: work 857 LOGICAL, INTENT(IN), OPTIONAL :: action 858 REAL(KIND=dp) :: res 859 860 INTEGER :: i, j, nsp, tline 861 LOGICAL :: nocut 862 REAL(KIND=dp) :: a, ar, arp, b, h26, q, s, v, z 863 REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), & 864 POINTER :: offdiagsplines 865 TYPE(spline_data_type), POINTER :: endpspline 866 867 endpspline => helium%u0 868 offdiagsplines => helium%uoffdiag 869 nocut = .TRUE. 870 IF (PRESENT(action)) THEN 871 IF (.NOT. action) THEN 872 endpspline => helium%e0 873 offdiagsplines => helium%eoffdiag 874 nocut = .FALSE. 875 END IF 876 END IF 877 878 ar = SQRT(helium_pbc_r2(helium, r)) 879 arp = SQRT(helium_pbc_r2(helium, rp)) 880 881 IF (helium%periodic .AND. ((ar > 0.5_dp*helium%cell_size) & 882 .OR. (arp > 0.5_dp*helium%cell_size))) THEN 883 v = 0.0_dp 884 IF (arp > 0.5_dp*helium%cell_size) THEN 885 IF (nocut) v = v + helium_spline(endpspline, 0.5_dp*helium%cell_size) 886 ELSE 887 v = v + helium_spline(endpspline, arp) 888 END IF 889 IF (ar > 0.5_dp*helium%cell_size) THEN 890 IF (nocut) v = v + helium_spline(endpspline, 0.5_dp*helium%cell_size) 891 ELSE 892 v = v + helium_spline(endpspline, ar) 893 END IF 894 res = 0.5_dp*v 895 ELSE 896 ! end-point action (first term): 897 v = 0.5_dp*(helium_spline(endpspline, ar) + helium_spline(endpspline, arp)) 898 s = helium_pbc_r2(helium, r - rp) 899 q = 0.5_dp*(ar + arp) 900 z = (ar - arp)**2 901 nsp = ((helium%pdx + 2)*(helium%pdx + 1))/2 - 1 902 a = endpspline%x1 903 b = endpspline%invh 904 IF (q < a) THEN 905 b = b*(q - a) 906 a = 1.0_dp - b 907 work(1:nsp) = a*offdiagsplines(1:nsp, 1, 1) + b*(offdiagsplines(1:nsp, 1, 2) - & 908 offdiagsplines(1:nsp, 2, 2)*endpspline%h26) 909 ELSE IF (q > endpspline%xn) THEN 910 b = b*(q - endpspline%xn) + 1.0_dp 911 a = 1.0_dp - b 912 tline = endpspline%n 913 work(1:nsp) = b*offdiagsplines(1:nsp, 1, tline) + a*(offdiagsplines(1:nsp, 1, tline - 1) - & 914 offdiagsplines(1:nsp, 2, tline - 1)*endpspline%h26) 915 ELSE 916 a = (a - q)*b 917 tline = INT(1.0_dp - a) 918 a = a + REAL(tline, kind=dp) 919 b = 1.0_dp - a 920 h26 = -a*b*endpspline%h26 921 work(1:nsp) = a*offdiagsplines(1:nsp, 1, tline) + b*offdiagsplines(1:nsp, 1, tline + 1) + & 922 h26*((a + 1.0_dp)*offdiagsplines(1:nsp, 2, tline) + (b + 1.0_dp)* & 923 offdiagsplines(1:nsp, 2, tline + 1)) 924 END IF 925 work(nsp + 1) = v 926 tline = 1 927 v = work(1) 928 DO i = 1, helium%pdx 929 v = v*z 930 tline = tline + 1 931 b = work(tline) 932 DO j = 1, i 933 tline = tline + 1 934 b = b*s + work(tline) 935 END DO 936 v = v + b 937 END DO 938 res = v 939 END IF 940 RETURN 941 END FUNCTION helium_eval_expansion 942 943! *************************************************************************** 944!> \brief Calculate the pair-product action or energy by evaluating the 945!> power series expansion according to Eq. 4.46 in Ceperley 1995. 946!> \param helium ... 947!> \param rchain ... 948!> \param nchain ... 949!> \param aij ... 950!> \param vcoef ... 951!> \param energy ... 952!> \return ... 953!> \author Harald Forbert 954!> \note This version calculates the action/energy for a chain segment 955! ************************************************************************************************** 956 FUNCTION helium_eval_chain(helium, rchain, nchain, aij, vcoef, energy) RESULT(res) 957 958 TYPE(helium_solvent_type), POINTER :: helium 959 INTEGER, INTENT(IN) :: nchain 960 REAL(KIND=dp), DIMENSION(3, nchain), INTENT(INOUT) :: rchain 961 REAL(KIND=dp), DIMENSION(nchain), INTENT(INOUT) :: aij 962 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: vcoef 963 LOGICAL, INTENT(IN), OPTIONAL :: energy 964 REAL(KIND=dp) :: res 965 966 INTEGER :: chainpos, i, j, ndim, nsp, tline 967 LOGICAL :: nocut 968 REAL(KIND=dp) :: a, ar, arp, b, ch, h26, q, s, totalv, v, & 969 z 970 REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), & 971 POINTER :: offdiagsplines 972 TYPE(spline_data_type), POINTER :: endpspline 973 974 endpspline => helium%u0 975 offdiagsplines => helium%uoffdiag 976 nocut = .TRUE. 977 IF (PRESENT(energy)) THEN 978 IF (energy) THEN 979 endpspline => helium%e0 980 offdiagsplines => helium%eoffdiag 981 nocut = .FALSE. 982 END IF 983 END IF 984 985 ndim = helium%pdx 986 nsp = ((ndim + 2)*(ndim + 1))/2 - 1 987 vcoef(1:nsp + 1) = 0.0_dp 988 totalv = 0.0_dp 989 DO i = 1, nchain 990 aij(i) = SQRT(helium_pbc_r2(helium, rchain(:, i))) 991 END DO 992 DO i = 2, nchain 993 rchain(:, i - 1) = rchain(:, i - 1) - rchain(:, i) 994 rchain(1, i - 1) = helium_pbc_r2(helium, rchain(:, i - 1)) 995 END DO 996 ! 997 IF (helium%periodic) THEN 998 ch = 0.5_dp*helium%cell_size 999 IF (nocut) THEN 1000 DO i = 2, nchain - 1 1001 totalv = totalv + helium_spline(endpspline, MIN(aij(i), ch)) 1002 END DO 1003 totalv = totalv + 0.5_dp*helium_spline(endpspline, MIN(aij(1), ch)) 1004 totalv = totalv + 0.5_dp*helium_spline(endpspline, MIN(aij(nchain), ch)) 1005 ELSE 1006 DO i = 2, nchain - 1 1007 ar = aij(i) 1008 IF (ar < ch) totalv = totalv + helium_spline(endpspline, ar) 1009 END DO 1010 ar = aij(1) 1011 IF (ar < ch) totalv = totalv + 0.5_dp*helium_spline(endpspline, ar) 1012 ar = aij(nchain) 1013 IF (ar < ch) totalv = totalv + 0.5_dp*helium_spline(endpspline, ar) 1014 END IF 1015 ELSE 1016 DO i = 2, nchain - 1 1017 totalv = totalv + helium_spline(endpspline, aij(i)) 1018 END DO 1019 totalv = totalv + 0.5_dp*helium_spline(endpspline, aij(1)) 1020 totalv = totalv + 0.5_dp*helium_spline(endpspline, aij(nchain)) 1021 END IF 1022 1023 DO chainpos = 1, nchain - 1 1024 ar = aij(chainpos) 1025 arp = aij(chainpos + 1) 1026 IF (helium%periodic .AND. ((ar > 0.5_dp*helium%cell_size) & 1027 .OR. (arp > 0.5_dp*helium%cell_size))) THEN 1028 CYCLE 1029 ELSE 1030 q = 0.5_dp*(ar + arp) 1031 s = rchain(1, chainpos) 1032 z = (ar - arp)**2 1033 a = endpspline%x1 1034 b = endpspline%invh 1035 IF (q < a) THEN 1036 b = b*(q - a) 1037 a = 1.0_dp - b 1038 vcoef(1:nsp) = a*offdiagsplines(1:nsp, 1, 1) + b*(offdiagsplines(1:nsp, 1, 2) - & 1039 offdiagsplines(1:nsp, 2, 2)*endpspline%h26) 1040 ELSE IF (q > endpspline%xn) THEN 1041 b = b*(q - endpspline%xn) + 1.0_dp 1042 a = 1.0_dp - b 1043 tline = endpspline%n 1044 vcoef(1:nsp) = b*offdiagsplines(1:nsp, 1, tline) + a*(offdiagsplines(1:nsp, 1, tline - 1) - & 1045 offdiagsplines(1:nsp, 2, tline - 1)*endpspline%h26) 1046 ELSE 1047 a = (a - q)*b 1048 tline = INT(1.0_dp - a) 1049 a = a + REAL(tline, kind=dp) 1050 b = 1.0_dp - a 1051 h26 = -a*b*endpspline%h26 1052 vcoef(1:nsp) = a*offdiagsplines(1:nsp, 1, tline) + b*offdiagsplines(1:nsp, 1, tline + 1) + & 1053 h26*((a + 1.0_dp)*offdiagsplines(1:nsp, 2, tline) + (b + 1.0_dp)* & 1054 offdiagsplines(1:nsp, 2, tline + 1)) 1055 END IF 1056 tline = 1 1057 v = vcoef(1) 1058 DO i = 1, ndim 1059 tline = tline + 1 1060 b = vcoef(tline) 1061 DO j = 1, i 1062 tline = tline + 1 1063 b = b*s + vcoef(tline) 1064 END DO 1065 v = v*z + b 1066 END DO 1067 totalv = totalv + v 1068 END IF 1069 END DO 1070 res = totalv 1071 RETURN 1072 END FUNCTION helium_eval_chain 1073 1074! ************************************************************************************************** 1075!> \brief ... 1076!> \param helium ... 1077!> \author Harald Forbert 1078! ************************************************************************************************** 1079 SUBROUTINE helium_update_transition_matrix(helium) 1080 1081 TYPE(helium_solvent_type), POINTER :: helium 1082 1083 INTEGER :: b, c, i, j, k, m, n, nb 1084 INTEGER, ALLOCATABLE, DIMENSION(:) :: lens, order 1085 INTEGER, DIMENSION(:), POINTER :: perm 1086 INTEGER, DIMENSION(:, :), POINTER :: nmatrix 1087 REAL(KIND=dp) :: f, q, t, v 1088 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: p 1089 REAL(KIND=dp), DIMENSION(3) :: r 1090 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ipmatrix, pmatrix, tmatrix 1091 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pos 1092 1093 nb = helium%atoms 1094 !TODO: check allocation status 1095 ALLOCATE (p(2*nb)) 1096 ALLOCATE (order(nb)) 1097 ALLOCATE (lens(2*nb)) 1098 b = helium%beads - helium%bisection + 1 1099 f = -0.5_dp/(helium%hb2m*helium%tau*helium%bisection) 1100 tmatrix => helium%tmatrix 1101 pmatrix => helium%pmatrix 1102 ipmatrix => helium%ipmatrix 1103 nmatrix => helium%nmatrix 1104 perm => helium%permutation 1105 pos => helium%pos 1106 DO i = 1, nb 1107 DO j = 1, nb 1108 v = 0.0_dp 1109 r(:) = pos(:, i, b) - pos(:, j, 1) 1110 CALL helium_pbc(helium, r) 1111 v = v + r(1)*r(1) + r(2)*r(2) + r(3)*r(3) 1112 pmatrix(i, j) = f*v 1113 END DO 1114 t = pmatrix(i, perm(i)) ! just some reference 1115 v = 0.0_dp 1116 DO j = 1, nb 1117 tmatrix(i, j) = EXP(pmatrix(i, j) - t) 1118 v = v + tmatrix(i, j) 1119 END DO 1120 ! normalize 1121 q = t + LOG(v) 1122 t = 1.0_dp/v 1123 DO j = 1, nb 1124 tmatrix(i, j) = tmatrix(i, j)*t 1125 ipmatrix(i, j) = 1.0_dp/tmatrix(i, j) 1126 END DO 1127 1128 ! at this point we have: 1129 ! tmatrix(i,j) = exp(-f*(r_i^b - r_j^1)**2) normalized such 1130 ! that sum_j tmatrix(i,j) = 1. 1131 ! ( tmatrix(k1,k2) = t_{k1,k2} / h_{k1} of ceperly. ) 1132 ! so tmatrix(i,j) is the probability to try to change a permutation 1133 ! with particle j (assuming particle i is already selected as well) 1134 ! ipmatrix(i,j) = 1.0/tmatrix(i,j) 1135 ! pmatrix(i,j) = log(tmatrix(i,j)) + some_offset(i) 1136 1137 ! generate optimal search tree so we can select which particle j 1138 ! belongs to a given random_number as fast as possible. 1139 ! (the traditional approach would be to generate a table 1140 ! of cumulative probabilities and to search that table) 1141 ! so for example if we have: 1142 ! tmatrix(i,:) = ( 0.1 , 0.4 , 0.2 , 0.3 ) 1143 ! traditionally we would build the running sum table: 1144 ! ( 0.1 , 0.5 , 0.7 , 1.0 ) and for a random number r 1145 ! would search this table for the lowest index larger than r 1146 ! (which would then be the particle index chosen by this random number) 1147 ! we build an optimal binary search tree instead, so here 1148 ! we would have: 1149 ! if ( r > 0.6 ) then take index 2, 1150 ! else if ( r > 0.3 ) then take index 4, 1151 ! else if ( r > 0.1 ) then take index 3 else index 1. 1152 ! the search tree is generated in tmatrix and nmatrix. 1153 ! tmatrix contains the decision values (0.6,0.3,0.1 in this case) 1154 ! and nmatrix contains the two branches (what to do if lower or higher) 1155 ! negative numbers in nmatrix mean take minus that index 1156 ! positive number means go down the tree to that next node, since we 1157 ! put the root of the tree at the end the arrays in the example would 1158 ! look like this: 1159 ! tmatrix(i,:) = ( 0.1 , 0.3 , 0.6 , arbitrary ) 1160 ! namtrix(i,:) = ( -1 , -3 , 1 , -4 , 2 , -2 , arb. , arb. ) 1161 ! 1162 ! the way to generate this tree may not be the best, but the 1163 ! tree generation itself shouldn't be needed quite that often: 1164 ! 1165 ! first sort values (with some variant of heap sort) 1166 1167 DO j = 1, nb 1168 order(j) = j 1169 p(j) = tmatrix(i, j) 1170 END DO 1171 IF (nb > 1) THEN ! if nb = 1 it is already sorted. 1172 k = nb/2 + 1 1173 c = nb 1174 DO 1175 IF (k > 1) THEN 1176 ! building up the heap: 1177 k = k - 1 1178 n = order(k) 1179 v = p(k) 1180 ELSE 1181 ! removing the top of the heap 1182 n = order(c) 1183 v = p(c) 1184 order(c) = order(1) 1185 p(c) = p(1) 1186 c = c - 1 1187 IF (c == 1) THEN 1188 order(1) = n 1189 p(1) = v 1190 EXIT 1191 END IF 1192 END IF 1193 m = k 1194 j = 2*k 1195 ! restoring heap order between k and c 1196 DO 1197 IF (j > c) EXIT 1198 IF (j < c) THEN 1199 IF (p(j) < p(j + 1)) j = j + 1 1200 END IF 1201 IF (v >= p(j)) EXIT 1202 order(m) = order(j) 1203 p(m) = p(j) 1204 m = j 1205 j = 2*j 1206 END DO 1207 order(m) = n 1208 p(m) = v 1209 END DO 1210 END IF 1211 1212 ! now: 1213 ! p(1:nb) : tmatrix(i,1:nb) sorted in ascending order 1214 ! order(1:nb) : corresponding index: p(j) == tmatrix(i,order(j)) 1215 ! for all j 1216 1217 ! merge sort with elements as we generate new interior search nodes 1218 ! by combining older elements/nodes 1219 1220 ! first fill unused part of array with guard values: 1221 DO j = nb + 1, 2*nb 1222 p(j) = 2.0_dp 1223 END DO 1224 1225 ! j - head of leaf queue 1226 ! c+1 - head of node queue in p (c in lens) 1227 ! m+1 - tail of node queue in p (m in lens) 1228 c = nb + 1 1229 j = 1 1230 DO m = nb + 1, 2*nb - 1 1231 ! get next smallest element 1232 IF (p(j) < p(c + 1)) THEN 1233 v = p(j) 1234 lens(j) = m 1235 j = j + 1 1236 ELSE 1237 v = p(c + 1) 1238 lens(c) = m 1239 c = c + 1 1240 END IF 1241 ! get the second next smallest element 1242 IF (p(j) < p(c + 1)) THEN 1243 p(m + 1) = v + p(j) 1244 lens(j) = m 1245 j = j + 1 1246 ELSE 1247 p(m + 1) = v + p(c + 1) 1248 lens(c) = m 1249 c = c + 1 1250 END IF 1251 END DO 1252 1253 ! lens(:) now has the tree with lens(j) pointing to its parent 1254 ! the root of the tree is at 2*nb-1 1255 ! calculate the depth of each node in the tree now: (root = 0) 1256 1257 lens(2*nb - 1) = 0 1258 DO m = 2*nb - 2, 1, -1 1259 lens(m) = lens(lens(m)) + 1 1260 END DO 1261 1262 ! lens(:) now has the depths of the nodes/leafs 1263 1264#if 0 1265 ! calculate average search depth (for information only) 1266 v = 0.0_dp 1267 DO j = 1, nb 1268 v = v + p(j)*lens(j) 1269 END DO 1270 PRINT *, "Expected number of comparisons with i=", i, v 1271#endif 1272 1273 ! reset the nodes, for the canonical tree we just need the leaf info 1274 DO j = 1, nb 1275 lens(j + nb) = 0 1276 p(j + nb) = 0.0_dp 1277 END DO 1278 1279 ! build the canonical tree (number of decisions on average are 1280 ! the same to the tree we build above, but it has better caching behavior 1281 1282 ! c head of leafs 1283 ! m head of interior nodes 1284 c = 1 1285 m = nb + 1 1286 DO k = 1, 2*nb - 2 1287 j = nb + 1 + (k - 1)/2 1288 IF (lens(c) > lens(m + 1)) THEN 1289 nmatrix(i, k) = -order(c) 1290 lens(j + 1) = lens(c) - 1 1291 v = p(c) 1292 c = c + 1 1293 ELSE 1294 nmatrix(i, k) = m - nb 1295 lens(j + 1) = lens(m + 1) - 1 1296 v = p(m) 1297 m = m + 1 1298 END IF 1299 p(j) = p(j) + v 1300 IF (MOD(k, 2) == 1) tmatrix(i, j - nb) = v 1301 END DO 1302 1303 ! now: 1304 ! nmatrix(i,2*j+1) left child of node j 1305 ! nmatrix(i,2*j+2) right child of node j 1306 ! children: 1307 ! negative : leaf with particle index == abs(value) 1308 ! positive : child node index 1309 ! p(j) weight of leaf j 1310 ! p(nb+j) weight of node j 1311 ! tmatrix(i,j) weight of left child of node j 1312 1313 ! fix offsets for decision tree: 1314 1315 p(nb - 1) = 0.0_dp 1316 DO m = nb - 1, 1, -1 1317 ! if right child is a node, set its offset and 1318 ! change its decision value 1319 IF (nmatrix(i, 2*m) > 0) THEN 1320 p(nmatrix(i, 2*m)) = tmatrix(i, m) 1321 tmatrix(i, nmatrix(i, 2*m)) = tmatrix(i, nmatrix(i, 2*m)) + tmatrix(i, m) 1322 END IF 1323 ! if left child is a node, set its offset and 1324 ! change its decision value 1325 IF (nmatrix(i, 2*m - 1) > 0) THEN 1326 p(nmatrix(i, 2*m - 1)) = p(m) 1327 tmatrix(i, nmatrix(i, 2*m - 1)) = tmatrix(i, nmatrix(i, 2*m - 1)) + p(m) 1328 END IF 1329 END DO 1330 1331 ! canonical optimal search tree done 1332 1333#if 0 1334 !some test code, to check if it gives the right distribution 1335 DO k = 1, nb 1336 p(k) = 1.0/ipmatrix(i, k) 1337 END DO 1338 lens(:) = 0 1339 ! number of random numbers to generate: 1340 c = 1000000000 1341 DO j = 1, c 1342 v = next_random_number(helium%rng_stream_uniform) 1343 ! walk down the search tree: 1344 k = nb - 1 1345 DO 1346 IF (tmatrix(i, k) > v) THEN 1347 k = nmatrix(i, 2*k - 1) 1348 ELSE 1349 k = nmatrix(i, 2*k) 1350 END IF 1351 IF (k < 0) EXIT 1352 END DO 1353 k = -k 1354 ! increment the counter for this particle index 1355 lens(k) = lens(k) + 1 1356 END DO 1357 ! search for maximum deviation from expectation value 1358 ! (relative to the expected variance) 1359 v = 0.0_dp 1360 k = -1 1361 DO j = 1, nb 1362 q = ABS((lens(j) - c*p(j))/SQRT(c*p(j))) 1363 !PRINT *,j,lens(j),c*p(j) 1364 IF (q > v) THEN 1365 v = q 1366 k = j 1367 END IF 1368 !PRINT *,lens(j),c*p(j),(lens(j)-c*p(j))/sqrt(c*p(j)) 1369 END DO 1370 PRINT *, "MAXDEV:", k, lens(k), c*p(k), v 1371 !PRINT *,"TMAT:",tmatrix(i,:) 1372 !PRINT *,"NMAT:",nmatrix(i,:) 1373 !STOP 1374#endif 1375#if 0 1376 !additional test code: 1377 p(:) = -1.0_dp 1378 p(nb - 1) = 0.0_dp 1379 p(2*nb - 1) = 1.0_dp 1380 DO j = nb - 1, 1, -1 1381 ! right child 1382 IF (nmatrix(i, 2*j) > 0) THEN 1383 c = nmatrix(i, 2*j) 1384 p(c) = tmatrix(i, j) 1385 p(c + nb) = p(j + nb) 1386 ELSE 1387 c = -nmatrix(i, 2*j) 1388 !PRINT *,c,1.0/ipmatrix(i,c),p(j+nb)-tmatrix(i,j) 1389 IF (ABS(1.0/ipmatrix(i, c) - (p(j + nb) - tmatrix(i, j))) > & 1390 10.0_dp*EPSILON(1.0_dp)) THEN 1391 PRINT *, "Probability mismatch for particle i->j", i, c 1392 PRINT *, "Got", p(j + nb) - tmatrix(i, j), "should be", 1.0/ipmatrix(i, c) 1393 STOP 1394 END IF 1395 END IF 1396 ! left child 1397 IF (nmatrix(i, 2*j - 1) > 0) THEN 1398 c = nmatrix(i, 2*j - 1) 1399 p(c + nb) = tmatrix(i, j) 1400 p(c) = p(j) 1401 ELSE 1402 c = -nmatrix(i, 2*j - 1) 1403 !PRINT *,c,1.0/ipmatrix(i,c),tmatrix(i,j)-p(j) 1404 IF (ABS(1.0/ipmatrix(i, c) - (tmatrix(i, j) - p(j))) > & 1405 10.0_dp*EPSILON(1.0_dp)) THEN 1406 PRINT *, "Probability mismatch for particle i->j", i, c 1407 PRINT *, "Got", tmatrix(i, j) - p(j), "should be", 1.0/ipmatrix(i, c) 1408 STOP 1409 END IF 1410 END IF 1411 END DO 1412 PRINT *, "Probabilities ok" 1413#endif 1414 1415 END DO 1416 1417 ! initialize trial permutation with some identity permutation 1418 ! (should not be taken, but just in case it does we have something valid) 1419 1420 helium%pweight = 0.0_dp 1421 t = next_random_number(helium%rng_stream_uniform) 1422 helium%ptable(1) = 1 + INT(t*nb) 1423 helium%ptable(2) = -1 1424 1425 ! recalculate inverse permutation table (just in case) 1426 DO i = 1, nb 1427 helium%iperm(perm(i)) = i 1428 END DO 1429 1430 ! clean up: 1431 DEALLOCATE (lens) 1432 DEALLOCATE (order) 1433 DEALLOCATE (p) 1434 1435 RETURN 1436 END SUBROUTINE helium_update_transition_matrix 1437 1438! ************************************************************************************************** 1439!> \brief ... 1440!> \param spl ... 1441!> \param xx ... 1442!> \return ... 1443!> \author Harald Forbert 1444! ************************************************************************************************** 1445 FUNCTION helium_spline(spl, xx) RESULT(res) 1446 TYPE(spline_data_type), POINTER :: spl 1447 REAL(KIND=dp), INTENT(IN) :: xx 1448 REAL(KIND=dp) :: res 1449 1450 REAL(KIND=dp) :: a, b 1451 1452 IF (xx < spl%x1) THEN 1453 b = spl%invh*(xx - spl%x1) 1454 a = 1.0_dp - b 1455 res = a*spl%y(1) + b*(spl%y(2) - spl%y2(2)*spl%h26) 1456 ELSE IF (xx > spl%xn) THEN 1457 b = spl%invh*(xx - spl%xn) + 1.0_dp 1458 a = 1.0_dp - b 1459 res = b*spl%y(spl%n) + a*(spl%y(spl%n - 1) - spl%y2(spl%n - 1)*spl%h26) 1460 ELSE 1461 res = spline_value(spl, xx) 1462 END IF 1463 RETURN 1464 END FUNCTION helium_spline 1465 1466! *************************************************************************** 1467!> \brief Return the distance <rij> between bead <ib> of atom <ia> 1468!> and bead <jb> of atom <ja>. 1469!> \param helium ... 1470!> \param ia ... 1471!> \param ib ... 1472!> \param ja ... 1473!> \param jb ... 1474!> \return ... 1475!> \date 2009-07-17 1476!> \author Lukasz Walewski 1477! ************************************************************************************************** 1478 FUNCTION helium_bead_rij(helium, ia, ib, ja, jb) RESULT(rij) 1479 1480 TYPE(helium_solvent_type), POINTER :: helium 1481 INTEGER, INTENT(IN) :: ia, ib, ja, jb 1482 REAL(KIND=dp) :: rij 1483 1484 REAL(KIND=dp) :: dx, dy, dz 1485 1486 dx = helium%pos(1, ia, ib) - helium%pos(1, ja, jb) 1487 dy = helium%pos(2, ia, ib) - helium%pos(2, ja, jb) 1488 dz = helium%pos(3, ia, ib) - helium%pos(3, ja, jb) 1489 rij = SQRT(dx*dx + dy*dy + dz*dz) 1490 1491 RETURN 1492 END FUNCTION helium_bead_rij 1493 1494! *************************************************************************** 1495!> \brief Given the atom number and permutation state return the cycle 1496!> number the atom belongs to. 1497!> \param helium ... 1498!> \param atom_number ... 1499!> \param permutation ... 1500!> \return ... 1501!> \date 2009-07-21 1502!> \author Lukasz Walewski 1503!> \note Cycles (or paths) are numbered from 1 to <num_cycles>, where 1504!> <num_cycles> is in the range of (1, <helium%atoms>). 1505!> if (num_cycles .EQ. 1) then all atoms belong to one cycle 1506!> if (num_cycles .EQ. helium%atoms) then there are no cycles of 1507!> length greater than 1 (i.e. no atoms are connected) 1508! ************************************************************************************************** 1509 FUNCTION helium_cycle_number(helium, atom_number, permutation) RESULT(cycle_number) 1510 1511 TYPE(helium_solvent_type), POINTER :: helium 1512 INTEGER, INTENT(IN) :: atom_number 1513 INTEGER, DIMENSION(:), POINTER :: permutation 1514 INTEGER :: cycle_number 1515 1516 INTEGER :: atom_idx, cycle_idx, cycle_num, ia, ib, & 1517 ic, num_cycles 1518 INTEGER, DIMENSION(:), POINTER :: cycle_index 1519 LOGICAL :: found, new_cycle 1520 1521 NULLIFY (cycle_index) 1522 cycle_index => helium%itmp_atoms_1d 1523 cycle_index(:) = 0 1524 1525 num_cycles = 0 1526 found = .FALSE. 1527 cycle_num = -1 1528 DO ia = 1, helium%atoms 1529 ! this loop reaches its maximum iteration count when atom in question 1530 ! is the last one (i.e. when atom_number .EQ. helium%atoms) 1531 1532 ! exit if we have found the cycle number for the atom in question 1533 IF (found) THEN 1534 EXIT 1535 END IF 1536 1537 ! initialize current cycle index with the current atom 1538 cycle_idx = ia 1539 1540 atom_idx = ia 1541 DO ib = 1, helium%atoms*helium%beads 1542 ! this loop reaches its maximum iteration count when all He atoms 1543 ! form one cycle (i.e. all beads belong to one path) 1544 1545 ! proceed along the path 1546 atom_idx = permutation(atom_idx) 1547 1548 IF (atom_idx .EQ. ia) THEN 1549 ! end of cycle detected (looped back to the first atom) 1550 1551 ! check if this is a new cycle 1552 new_cycle = .TRUE. 1553 DO ic = 1, num_cycles 1554 IF (cycle_index(ic) .EQ. cycle_idx) THEN 1555 new_cycle = .FALSE. 1556 END IF 1557 END DO 1558 1559 IF (new_cycle) THEN 1560 ! increase number of cycles and update the current cycle's index 1561 num_cycles = num_cycles + 1 1562 cycle_index(num_cycles) = cycle_idx 1563 END IF 1564 1565 ! if this was the atom in question 1566 IF (ia .EQ. atom_number) THEN 1567 ! save the cycle index it belongs to 1568 cycle_num = cycle_idx 1569 1570 ! exit the loop over atoms, we've found what we've been looking for 1571 found = .TRUE. 1572 END IF 1573 1574 ! exit the loop over beads, there are no more (new) beads in this cycle 1575 EXIT 1576 END IF 1577 1578 ! set the cycle index to the lowest atom index in this cycle 1579 IF (atom_idx .LT. cycle_idx) THEN 1580 cycle_idx = atom_idx 1581 END IF 1582 1583 END DO 1584 1585 END DO 1586 1587!TODO make it cp2k way 1588 IF (.NOT. found) THEN 1589 CPWARN("helium_cycle_number: we are going to return -1, problems ahead!") 1590 END IF 1591 1592 ! at this point we know the cycle index for atom <atom_number> 1593 ! but it is expressed as the atom number of the first atom in that cycle 1594 1595 ! renumber cycle indices, so that they form a range (1, <num_cycles>) 1596 ! (don't do it actually - just return the requested <cycle_number>) 1597 cycle_number = 0 1598 DO ic = 1, num_cycles 1599 IF (cycle_index(ic) .EQ. cycle_num) THEN 1600 cycle_number = ic 1601 EXIT 1602 END IF 1603 END DO 1604 1605 NULLIFY (cycle_index) 1606 1607 RETURN 1608 END FUNCTION helium_cycle_number 1609 1610! *************************************************************************** 1611!> \brief Given the atom number and permutation state return the length of 1612!> the path this atom belongs to. 1613!> \param helium ... 1614!> \param atom_number ... 1615!> \param permutation ... 1616!> \return ... 1617!> \date 2009-10-07 1618!> \author Lukasz Walewski 1619! ************************************************************************************************** 1620 FUNCTION helium_path_length(helium, atom_number, permutation) RESULT(path_length) 1621 1622 TYPE(helium_solvent_type), POINTER :: helium 1623 INTEGER, INTENT(IN) :: atom_number 1624 INTEGER, DIMENSION(:), POINTER :: permutation 1625 INTEGER :: path_length 1626 1627 INTEGER :: atom_idx, ia 1628 LOGICAL :: path_end_reached 1629 1630 atom_idx = atom_number 1631 path_length = 0 1632 path_end_reached = .FALSE. 1633 DO ia = 1, helium%atoms 1634 path_length = path_length + 1 1635 atom_idx = permutation(atom_idx) 1636 IF (atom_idx .EQ. atom_number) THEN 1637 path_end_reached = .TRUE. 1638 EXIT 1639 END IF 1640 END DO 1641 1642 IF (.NOT. path_end_reached) THEN 1643 path_length = -1 1644 END IF 1645 1646 RETURN 1647 END FUNCTION helium_path_length 1648 1649! *************************************************************************** 1650!> \brief Given an element of a permutation return the cycle it belongs to. 1651!> \param element ... 1652!> \param permutation ... 1653!> \return ... 1654!> \date 2013-12-10 1655!> \author Lukasz Walewski 1656!> \note This function allocates memory and returns a pointer to it, 1657!> do not forget to deallocate once finished with the results. 1658! ************************************************************************************************** 1659 FUNCTION helium_cycle_of(element, permutation) RESULT(CYCLE) 1660 1661 INTEGER, INTENT(IN) :: element 1662 INTEGER, DIMENSION(:), INTENT(IN), POINTER :: permutation 1663 INTEGER, DIMENSION(:), POINTER :: CYCLE 1664 1665 INTEGER :: ia, icur, len, nsize 1666 CHARACTER(len=*), PARAMETER :: routineN = 'helium_cycle_of', & 1667 routineP = moduleN//':'//routineN 1668 1669 INTEGER, DIMENSION(:), POINTER :: my_cycle 1670 LOGICAL :: cycle_end_reached 1671 1672 nsize = SIZE(permutation) 1673 1674 ! maximum possible cycle length is the number of atoms 1675 NULLIFY (my_cycle) 1676 ALLOCATE (my_cycle(nsize)) 1677 1678 ! traverse the permutation table 1679 len = 0 1680 icur = element 1681 cycle_end_reached = .FALSE. 1682 DO ia = 1, nsize 1683 len = len + 1 1684 my_cycle(len) = icur 1685 icur = permutation(icur) 1686 IF (icur .EQ. element) THEN 1687 cycle_end_reached = .TRUE. 1688 EXIT 1689 END IF 1690 END DO 1691 1692 IF (.NOT. cycle_end_reached) THEN 1693 ! return null 1694 NULLIFY (CYCLE) 1695 ELSE 1696 ! assign the result 1697 ALLOCATE (CYCLE(len)) 1698 CYCLE(1:len) = my_cycle(1:len) 1699 END IF 1700 1701 ! clean up 1702 DEALLOCATE (my_cycle) 1703 NULLIFY (my_cycle) 1704 1705 RETURN 1706 END FUNCTION helium_cycle_of 1707 1708! *************************************************************************** 1709!> \brief Return total winding number 1710!> \param helium ... 1711!> \return ... 1712!> \date 2014-04-24 1713!> \author Lukasz Walewski 1714! ************************************************************************************************** 1715 FUNCTION helium_total_winding_number(helium) RESULT(wnum) 1716 1717 TYPE(helium_solvent_type), POINTER :: helium 1718 REAL(KIND=dp), DIMENSION(3) :: wnum 1719 1720 CHARACTER(len=*), PARAMETER :: routineN = 'helium_total_winding_number', & 1721 routineP = moduleN//':'//routineN 1722 1723 INTEGER :: ia, ib 1724 REAL(KIND=dp), DIMENSION(3) :: rcur 1725 REAL(KIND=dp), DIMENSION(:), POINTER :: ri, rj 1726 1727 wnum(:) = 0.0_dp 1728 DO ia = 1, helium%atoms 1729 ! sum of contributions from the rest of bead pairs 1730 DO ib = 1, helium%beads - 1 1731 ri => helium%pos(:, ia, ib) 1732 rj => helium%pos(:, ia, ib + 1) 1733 rcur(:) = ri(:) - rj(:) 1734 CALL helium_pbc(helium, rcur) 1735 wnum(:) = wnum(:) + rcur(:) 1736 END DO 1737 ! contribution from the last and the first bead 1738 ri => helium%pos(:, ia, helium%beads) 1739 rj => helium%pos(:, helium%permutation(ia), 1) 1740 rcur(:) = ri(:) - rj(:) 1741 CALL helium_pbc(helium, rcur) 1742 wnum(:) = wnum(:) + rcur(:) 1743 END DO 1744 1745 END FUNCTION helium_total_winding_number 1746 1747! *************************************************************************** 1748!> \brief Return link winding number 1749!> \param helium ... 1750!> \param ia ... 1751!> \param ib ... 1752!> \return ... 1753!> \date 2014-04-24 1754!> \author Lukasz Walewski 1755! ************************************************************************************************** 1756 FUNCTION helium_link_winding_number(helium, ia, ib) RESULT(wnum) 1757 1758 TYPE(helium_solvent_type), POINTER :: helium 1759 INTEGER :: ia, ib 1760 REAL(KIND=dp), DIMENSION(3) :: wnum 1761 1762 CHARACTER(len=*), PARAMETER :: routineN = 'helium_link_winding_number', & 1763 routineP = moduleN//':'//routineN 1764 1765 INTEGER :: ja1, ja2, jb1, jb2 1766 REAL(KIND=dp), DIMENSION(:), POINTER :: r1, r2 1767 1768 IF (ib .EQ. helium%beads) THEN 1769 ja1 = ia 1770 ja2 = helium%permutation(ia) 1771 jb1 = ib 1772 jb2 = 1 1773 ELSE 1774 ja1 = ia 1775 ja2 = ia 1776 jb1 = ib 1777 jb2 = ib + 1 1778 END IF 1779 r1 => helium%pos(:, ja1, jb1) 1780 r2 => helium%pos(:, ja2, jb2) 1781 wnum(:) = r1(:) - r2(:) 1782 CALL helium_pbc(helium, wnum) 1783 1784 RETURN 1785 END FUNCTION helium_link_winding_number 1786 1787! *************************************************************************** 1788!> \brief Return total winding number (sum over all links) 1789!> \param helium ... 1790!> \return ... 1791!> \date 2014-04-24 1792!> \author Lukasz Walewski 1793!> \note mostly for sanity checks 1794! ************************************************************************************************** 1795 FUNCTION helium_total_winding_number_linkwise(helium) RESULT(wnum) 1796 1797 TYPE(helium_solvent_type), POINTER :: helium 1798 REAL(KIND=dp), DIMENSION(3) :: wnum 1799 1800 CHARACTER(len=*), PARAMETER :: routineN = 'helium_total_winding_number_linkwise', & 1801 routineP = moduleN//':'//routineN 1802 1803 INTEGER :: ia, ib 1804 1805 wnum(:) = 0.0_dp 1806 DO ia = 1, helium%atoms 1807 DO ib = 1, helium%beads 1808 wnum(:) = wnum(:) + helium_link_winding_number(helium, ia, ib) 1809 END DO 1810 END DO 1811 1812 RETURN 1813 END FUNCTION helium_total_winding_number_linkwise 1814 1815! *************************************************************************** 1816!> \brief Return cycle winding number 1817!> \param helium ... 1818!> \param CYCLE ... 1819!> \param pos ... 1820!> \return ... 1821!> \date 2014-04-24 1822!> \author Lukasz Walewski 1823! ************************************************************************************************** 1824 FUNCTION helium_cycle_winding_number(helium, CYCLE, pos) RESULT(wnum) 1825 1826 TYPE(helium_solvent_type), POINTER :: helium 1827 INTEGER, DIMENSION(:), POINTER :: CYCLE 1828 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pos 1829 REAL(KIND=dp), DIMENSION(3) :: wnum 1830 1831 CHARACTER(len=*), PARAMETER :: routineN = 'helium_cycle_winding_number', & 1832 routineP = moduleN//':'//routineN 1833 1834 INTEGER :: i1, i2, ia, ib, nsize 1835 REAL(KIND=dp), DIMENSION(3) :: rcur 1836 REAL(KIND=dp), DIMENSION(:), POINTER :: ri, rj 1837 1838 nsize = SIZE(CYCLE) 1839 1840 ! traverse the path 1841 wnum(:) = 0.0_dp 1842 DO ia = 1, nsize 1843 ! contributions from all bead pairs of the current atom 1844 DO ib = 1, helium%beads - 1 1845 ri => pos(:, CYCLE(ia), ib) 1846 rj => pos(:, CYCLE(ia), ib + 1) 1847 rcur(:) = ri(:) - rj(:) 1848 CALL helium_pbc(helium, rcur) 1849 wnum(:) = wnum(:) + rcur(:) 1850 END DO 1851 ! contribution from the last bead of the current atom 1852 ! and the first bead of the next atom 1853 i1 = CYCLE(ia) 1854 IF (ia .EQ. nsize) THEN 1855 i2 = CYCLE(1) 1856 ELSE 1857 i2 = CYCLE(ia + 1) 1858 END IF 1859 ri => pos(:, i1, helium%beads) 1860 rj => pos(:, i2, 1) 1861 rcur(:) = ri(:) - rj(:) 1862 CALL helium_pbc(helium, rcur) 1863 wnum(:) = wnum(:) + rcur(:) 1864 END DO 1865 1866 RETURN 1867 END FUNCTION helium_cycle_winding_number 1868 1869! *************************************************************************** 1870!> \brief Return total winding number (sum over all cycles) 1871!> \param helium ... 1872!> \return ... 1873!> \date 2014-04-24 1874!> \author Lukasz Walewski 1875!> \note mostly for sanity checks 1876! ************************************************************************************************** 1877 FUNCTION helium_total_winding_number_cyclewise(helium) RESULT(wnum) 1878 1879 TYPE(helium_solvent_type), POINTER :: helium 1880 REAL(KIND=dp), DIMENSION(3) :: wnum 1881 1882 CHARACTER(len=*), PARAMETER :: routineN = 'helium_total_winding_number_cyclewise', & 1883 routineP = moduleN//':'//routineN 1884 1885 INTEGER :: ic 1886 REAL(KIND=dp), DIMENSION(3) :: wn 1887 TYPE(int_arr_ptr), DIMENSION(:), POINTER :: cycles 1888 1889! decompose the current permutation state into permutation cycles 1890 1891 NULLIFY (cycles) 1892 cycles => helium_calc_cycles(helium%permutation) 1893 1894 wnum(:) = 0.0_dp 1895 DO ic = 1, SIZE(cycles) 1896 wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos) 1897 wnum(:) = wnum(:) + wn(:) 1898 END DO 1899 1900 DEALLOCATE (cycles) 1901 1902 RETURN 1903 END FUNCTION helium_total_winding_number_cyclewise 1904 1905! *************************************************************************** 1906!> \brief Return total projected area 1907!> \param helium ... 1908!> \return ... 1909!> \date 2014-04-24 1910!> \author Lukasz Walewski 1911! ************************************************************************************************** 1912 FUNCTION helium_total_projected_area(helium) RESULT(area) 1913 1914 TYPE(helium_solvent_type), POINTER :: helium 1915 REAL(KIND=dp), DIMENSION(3) :: area 1916 1917 CHARACTER(len=*), PARAMETER :: routineN = 'helium_total_projected_area', & 1918 routineP = moduleN//':'//routineN 1919 1920 INTEGER :: ia, ib 1921 REAL(KIND=dp), DIMENSION(3) :: r1, r12, r2, rcur 1922 1923 area(:) = 0.0_dp 1924 DO ia = 1, helium%atoms 1925 ! contributions from all links of the current atom 1926 DO ib = 1, helium%beads - 1 1927 r1(:) = helium%pos(:, ia, ib) 1928 r2(:) = helium%pos(:, ia, ib + 1) 1929 ! comment out for non-PBC version --> 1930 r12(:) = r2(:) - r1(:) 1931 CALL helium_pbc(helium, r1) 1932 CALL helium_pbc(helium, r12) 1933 r2(:) = r1(:) + r12(:) 1934 ! comment out for non-PBC version <-- 1935 rcur(1) = r1(2)*r2(3) - r1(3)*r2(2) 1936 rcur(2) = r1(3)*r2(1) - r1(1)*r2(3) 1937 rcur(3) = r1(1)*r2(2) - r1(2)*r2(1) 1938 area(:) = area(:) + rcur(:) 1939 END DO 1940 ! contribution from the last bead of the current atom 1941 ! and the first bead of the next atom 1942 r1(:) = helium%pos(:, ia, helium%beads) 1943 r2(:) = helium%pos(:, helium%permutation(ia), 1) 1944 ! comment out for non-PBC version --> 1945 r12(:) = r2(:) - r1(:) 1946 CALL helium_pbc(helium, r1) 1947 CALL helium_pbc(helium, r12) 1948 r2(:) = r1(:) + r12(:) 1949 ! comment out for non-PBC version <-- 1950 rcur(1) = r1(2)*r2(3) - r1(3)*r2(2) 1951 rcur(2) = r1(3)*r2(1) - r1(1)*r2(3) 1952 rcur(3) = r1(1)*r2(2) - r1(2)*r2(1) 1953 area(:) = area(:) + rcur(:) 1954 END DO 1955 area(:) = 0.5_dp*area(:) 1956 1957 RETURN 1958 END FUNCTION helium_total_projected_area 1959 1960! *************************************************************************** 1961!> \brief Return link projected area 1962!> \param helium ... 1963!> \param ia ... 1964!> \param ib ... 1965!> \return ... 1966!> \date 2014-04-24 1967!> \author Lukasz Walewski 1968! ************************************************************************************************** 1969 FUNCTION helium_link_projected_area(helium, ia, ib) RESULT(area) 1970 1971 TYPE(helium_solvent_type), POINTER :: helium 1972 INTEGER :: ia, ib 1973 REAL(KIND=dp), DIMENSION(3) :: area 1974 1975 CHARACTER(len=*), PARAMETER :: routineN = 'helium_link_projected_area', & 1976 routineP = moduleN//':'//routineN 1977 1978 INTEGER :: ja1, ja2, jb1, jb2 1979 REAL(KIND=dp), DIMENSION(3) :: r1, r12, r2 1980 1981 IF (ib .EQ. helium%beads) THEN 1982 ja1 = ia 1983 ja2 = helium%permutation(ia) 1984 jb1 = ib 1985 jb2 = 1 1986 ELSE 1987 ja1 = ia 1988 ja2 = ia 1989 jb1 = ib 1990 jb2 = ib + 1 1991 END IF 1992 r1(:) = helium%pos(:, ja1, jb1) 1993 r2(:) = helium%pos(:, ja2, jb2) 1994 ! comment out for non-PBC version --> 1995 r12(:) = r2(:) - r1(:) 1996 CALL helium_pbc(helium, r1) 1997 CALL helium_pbc(helium, r12) 1998 r2(:) = r1(:) + r12(:) 1999 ! comment out for non-PBC version <-- 2000 area(1) = r1(2)*r2(3) - r1(3)*r2(2) 2001 area(2) = r1(3)*r2(1) - r1(1)*r2(3) 2002 area(3) = r1(1)*r2(2) - r1(2)*r2(1) 2003 area(:) = 0.5_dp*area(:) 2004 2005 RETURN 2006 END FUNCTION helium_link_projected_area 2007 2008! *************************************************************************** 2009!> \brief Return total projected area (sum over all links) 2010!> \param helium ... 2011!> \return ... 2012!> \date 2014-04-24 2013!> \author Lukasz Walewski 2014!> \note mostly for sanity checks 2015! ************************************************************************************************** 2016 FUNCTION helium_total_projected_area_linkwise(helium) RESULT(area) 2017 2018 TYPE(helium_solvent_type), POINTER :: helium 2019 REAL(KIND=dp), DIMENSION(3) :: area 2020 2021 CHARACTER(len=*), PARAMETER :: routineN = 'helium_total_projected_area_linkwise', & 2022 routineP = moduleN//':'//routineN 2023 2024 INTEGER :: ia, ib 2025 2026 area(:) = 0.0_dp 2027 DO ia = 1, helium%atoms 2028 DO ib = 1, helium%beads 2029 area(:) = area(:) + helium_link_projected_area(helium, ia, ib) 2030 END DO 2031 END DO 2032 2033 END FUNCTION helium_total_projected_area_linkwise 2034 2035! *************************************************************************** 2036!> \brief Return cycle projected area 2037!> \param helium ... 2038!> \param CYCLE ... 2039!> \return ... 2040!> \date 2014-04-24 2041!> \author Lukasz Walewski 2042! ************************************************************************************************** 2043 FUNCTION helium_cycle_projected_area(helium, CYCLE) RESULT(area) 2044 2045 TYPE(helium_solvent_type), POINTER :: helium 2046 INTEGER, DIMENSION(:), POINTER :: CYCLE 2047 REAL(KIND=dp), DIMENSION(3) :: area 2048 2049 CHARACTER(len=*), PARAMETER :: routineN = 'helium_cycle_projected_area', & 2050 routineP = moduleN//':'//routineN 2051 2052 INTEGER :: i1, i2, ia, ib, nsize 2053 REAL(KIND=dp), DIMENSION(3) :: rcur, rsum 2054 REAL(KIND=dp), DIMENSION(:), POINTER :: ri, rj 2055 2056 nsize = SIZE(CYCLE) 2057 2058 ! traverse the path 2059 rsum(:) = 0.0_dp 2060 DO ia = 1, nsize 2061 ! contributions from all bead pairs of the current atom 2062 DO ib = 1, helium%beads - 1 2063 ri => helium%pos(:, CYCLE(ia), ib) 2064 rj => helium%pos(:, CYCLE(ia), ib + 1) 2065 rcur(1) = ri(2)*rj(3) - ri(3)*rj(2) 2066 rcur(2) = ri(3)*rj(1) - ri(1)*rj(3) 2067 rcur(3) = ri(1)*rj(2) - ri(2)*rj(1) 2068 rsum(:) = rsum(:) + rcur(:) 2069 END DO 2070 ! contribution from the last bead of the current atom 2071 ! and the first bead of the next atom 2072 i1 = CYCLE(ia) 2073 IF (ia .EQ. nsize) THEN 2074 i2 = CYCLE(1) 2075 ELSE 2076 i2 = CYCLE(ia + 1) 2077 END IF 2078 ri => helium%pos(:, i1, helium%beads) 2079 rj => helium%pos(:, i2, 1) 2080 rcur(1) = ri(2)*rj(3) - ri(3)*rj(2) 2081 rcur(2) = ri(3)*rj(1) - ri(1)*rj(3) 2082 rcur(3) = ri(1)*rj(2) - ri(2)*rj(1) 2083 rsum(:) = rsum(:) + rcur(:) 2084 END DO 2085 area(:) = 0.5_dp*rsum(:) 2086 2087 RETURN 2088 END FUNCTION helium_cycle_projected_area 2089 2090! *************************************************************************** 2091!> \brief Return cycle projected area (sum over all links) 2092!> \param helium ... 2093!> \param CYCLE ... 2094!> \return ... 2095!> \date 2014-04-24 2096!> \author Lukasz Walewski 2097!> \note mostly for sanity checks 2098! ************************************************************************************************** 2099 FUNCTION helium_cycle_projected_area_pbc(helium, CYCLE) RESULT(area) 2100 2101 TYPE(helium_solvent_type), POINTER :: helium 2102 INTEGER, DIMENSION(:), POINTER :: CYCLE 2103 REAL(KIND=dp), DIMENSION(3) :: area 2104 2105 CHARACTER(len=*), PARAMETER :: routineN = 'helium_cycle_projected_area_pbc', & 2106 routineP = moduleN//':'//routineN 2107 2108 INTEGER :: i1, i2, ia, ib, nsize 2109 REAL(KIND=dp), DIMENSION(3) :: r1, r12, r2, rcur 2110 2111 nsize = SIZE(CYCLE) 2112 2113 ! traverse the path 2114 area(:) = 0.0_dp 2115 DO ia = 1, nsize 2116 ! contributions from all bead pairs of the current atom 2117 DO ib = 1, helium%beads - 1 2118 r1(:) = helium%pos(:, CYCLE(ia), ib) 2119 r2(:) = helium%pos(:, CYCLE(ia), ib + 1) 2120 r12(:) = r2(:) - r1(:) 2121 CALL helium_pbc(helium, r1) 2122 CALL helium_pbc(helium, r12) 2123 r2(:) = r1(:) + r12(:) 2124 rcur(1) = r1(2)*r2(3) - r1(3)*r2(2) 2125 rcur(2) = r1(3)*r2(1) - r1(1)*r2(3) 2126 rcur(3) = r1(1)*r2(2) - r1(2)*r2(1) 2127 area(:) = area(:) + rcur(:) 2128 END DO 2129 ! contribution from the last bead of the current atom 2130 ! and the first bead of the next atom 2131 i1 = CYCLE(ia) 2132 IF (ia .EQ. nsize) THEN 2133 i2 = CYCLE(1) 2134 ELSE 2135 i2 = CYCLE(ia + 1) 2136 END IF 2137 r1(:) = helium%pos(:, i1, helium%beads) 2138 r2(:) = helium%pos(:, i2, 1) 2139 r12(:) = r2(:) - r1(:) 2140 CALL helium_pbc(helium, r1) 2141 CALL helium_pbc(helium, r12) 2142 r2(:) = r1(:) + r12(:) 2143 rcur(1) = r1(2)*r2(3) - r1(3)*r2(2) 2144 rcur(2) = r1(3)*r2(1) - r1(1)*r2(3) 2145 rcur(3) = r1(1)*r2(2) - r1(2)*r2(1) 2146 area(:) = area(:) + rcur(:) 2147 END DO 2148 area(:) = 0.5_dp*area(:) 2149 2150 RETURN 2151 END FUNCTION helium_cycle_projected_area_pbc 2152 2153! *************************************************************************** 2154!> \brief Return total projected area (sum over all cycles) 2155!> \param helium ... 2156!> \return ... 2157!> \date 2014-04-24 2158!> \author Lukasz Walewski 2159!> \note mostly for sanity checks 2160! ************************************************************************************************** 2161 FUNCTION helium_total_projected_area_cyclewise(helium) RESULT(area) 2162 2163 TYPE(helium_solvent_type), POINTER :: helium 2164 REAL(KIND=dp), DIMENSION(3) :: area 2165 2166 CHARACTER(len=*), PARAMETER :: routineN = 'helium_total_projected_area_cyclewise', & 2167 routineP = moduleN//':'//routineN 2168 2169 INTEGER :: ic 2170 REAL(KIND=dp), DIMENSION(3) :: pa 2171 TYPE(int_arr_ptr), DIMENSION(:), POINTER :: cycles 2172 2173! decompose the current permutation state into permutation cycles 2174 2175 NULLIFY (cycles) 2176 cycles => helium_calc_cycles(helium%permutation) 2177 2178 area(:) = 0.0_dp 2179 DO ic = 1, SIZE(cycles) 2180 pa = helium_cycle_projected_area(helium, cycles(ic)%iap) 2181 area(:) = area(:) + pa(:) 2182 END DO 2183 2184 RETURN 2185 END FUNCTION helium_total_projected_area_cyclewise 2186 2187! *************************************************************************** 2188!> \brief Return total moment of inertia divided by m_He 2189!> \param helium ... 2190!> \return ... 2191!> \date 2014-04-24 2192!> \author Lukasz Walewski 2193! ************************************************************************************************** 2194 FUNCTION helium_total_moment_of_inertia(helium) RESULT(moit) 2195 2196 TYPE(helium_solvent_type), POINTER :: helium 2197 REAL(KIND=dp), DIMENSION(3) :: moit 2198 2199 CHARACTER(len=*), PARAMETER :: routineN = 'helium_total_moment_of_inertia', & 2200 routineP = moduleN//':'//routineN 2201 2202 INTEGER :: ia, ib 2203 REAL(KIND=dp), DIMENSION(3) :: com, r1, r12, r2, rcur 2204 2205 com(:) = helium%center(:) 2206 2207 moit(:) = 0.0_dp 2208 DO ia = 1, helium%atoms 2209 ! contributions from all the links of the current atom 2210 DO ib = 1, helium%beads - 1 2211 r1(:) = helium%pos(:, ia, ib) - com(:) 2212 r2(:) = helium%pos(:, ia, ib + 1) - com(:) 2213 ! comment out for non-PBC version --> 2214 r12(:) = r2(:) - r1(:) 2215 CALL helium_pbc(helium, r1) 2216 CALL helium_pbc(helium, r12) 2217 r2(:) = r1(:) + r12(:) 2218 ! comment out for non-PBC version <-- 2219 rcur(1) = r1(2)*r2(2) + r1(3)*r2(3) 2220 rcur(2) = r1(3)*r2(3) + r1(1)*r2(1) 2221 rcur(3) = r1(1)*r2(1) + r1(2)*r2(2) 2222 moit(:) = moit(:) + rcur(:) 2223 END DO 2224 ! contribution from the last bead of the current atom 2225 ! and the first bead of the next atom 2226 r1(:) = helium%pos(:, ia, helium%beads) - com(:) 2227 r2(:) = helium%pos(:, helium%permutation(ia), 1) - com(:) 2228 ! comment out for non-PBC version --> 2229 r12(:) = r2(:) - r1(:) 2230 CALL helium_pbc(helium, r1) 2231 CALL helium_pbc(helium, r12) 2232 r2(:) = r1(:) + r12(:) 2233 ! comment out for non-PBC version <-- 2234 rcur(1) = r1(2)*r2(2) + r1(3)*r2(3) 2235 rcur(2) = r1(3)*r2(3) + r1(1)*r2(1) 2236 rcur(3) = r1(1)*r2(1) + r1(2)*r2(2) 2237 moit(:) = moit(:) + rcur(:) 2238 END DO 2239 moit(:) = moit(:)/helium%beads 2240 2241 RETURN 2242 END FUNCTION helium_total_moment_of_inertia 2243 2244! *************************************************************************** 2245!> \brief Return link moment of inertia divided by m_He 2246!> \param helium ... 2247!> \param ia ... 2248!> \param ib ... 2249!> \return ... 2250!> \date 2014-04-24 2251!> \author Lukasz Walewski 2252! ************************************************************************************************** 2253 FUNCTION helium_link_moment_of_inertia(helium, ia, ib) RESULT(moit) 2254 2255 TYPE(helium_solvent_type), POINTER :: helium 2256 INTEGER :: ia, ib 2257 REAL(KIND=dp), DIMENSION(3) :: moit 2258 2259 CHARACTER(len=*), PARAMETER :: routineN = 'helium_link_moment_of_inertia', & 2260 routineP = moduleN//':'//routineN 2261 2262 INTEGER :: ja1, ja2, jb1, jb2 2263 REAL(KIND=dp), DIMENSION(3) :: com, r1, r12, r2 2264 2265 com(:) = helium%center(:) 2266 2267 IF (ib .EQ. helium%beads) THEN 2268 ja1 = ia 2269 ja2 = helium%permutation(ia) 2270 jb1 = ib 2271 jb2 = 1 2272 ELSE 2273 ja1 = ia 2274 ja2 = ia 2275 jb1 = ib 2276 jb2 = ib + 1 2277 END IF 2278 r1(:) = helium%pos(:, ja1, jb1) - com(:) 2279 r2(:) = helium%pos(:, ja2, jb2) - com(:) 2280 ! comment out for non-PBC version --> 2281 r12(:) = r2(:) - r1(:) 2282 CALL helium_pbc(helium, r1) 2283 CALL helium_pbc(helium, r12) 2284 r2(:) = r1(:) + r12(:) 2285 ! comment out for non-PBC version <-- 2286 moit(1) = r1(2)*r2(2) + r1(3)*r2(3) 2287 moit(2) = r1(3)*r2(3) + r1(1)*r2(1) 2288 moit(3) = r1(1)*r2(1) + r1(2)*r2(2) 2289 moit(:) = moit(:)/helium%beads 2290 2291 RETURN 2292 END FUNCTION helium_link_moment_of_inertia 2293 2294! *************************************************************************** 2295!> \brief Return total moment of inertia (sum over all links) 2296!> \param helium ... 2297!> \return ... 2298!> \date 2014-04-24 2299!> \author Lukasz Walewski 2300!> \note mostly for sanity checks 2301! ************************************************************************************************** 2302 FUNCTION helium_total_moment_of_inertia_linkwise(helium) RESULT(moit) 2303 2304 TYPE(helium_solvent_type), POINTER :: helium 2305 REAL(KIND=dp), DIMENSION(3) :: moit 2306 2307 CHARACTER(len=*), PARAMETER :: routineN = 'helium_total_moment_of_inertia_linkwise', & 2308 routineP = moduleN//':'//routineN 2309 2310 INTEGER :: ia, ib 2311 2312 moit(:) = 0.0_dp 2313 DO ia = 1, helium%atoms 2314 DO ib = 1, helium%beads 2315 moit(:) = moit(:) + helium_link_moment_of_inertia(helium, ia, ib) 2316 END DO 2317 END DO 2318 2319 END FUNCTION helium_total_moment_of_inertia_linkwise 2320 2321! *************************************************************************** 2322!> \brief Return moment of inertia of a cycle divided by m_He 2323!> \param helium ... 2324!> \param CYCLE ... 2325!> \param pos ... 2326!> \return ... 2327!> \date 2014-04-24 2328!> \author Lukasz Walewski 2329! ************************************************************************************************** 2330 FUNCTION helium_cycle_moment_of_inertia(helium, CYCLE, pos) RESULT(moit) 2331 2332 TYPE(helium_solvent_type), POINTER :: helium 2333 INTEGER, DIMENSION(:), POINTER :: CYCLE 2334 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pos 2335 REAL(KIND=dp), DIMENSION(3) :: moit 2336 2337 CHARACTER(len=*), PARAMETER :: routineN = 'helium_cycle_moment_of_inertia', & 2338 routineP = moduleN//':'//routineN 2339 2340 INTEGER :: i1, i2, ia, ib, nsize 2341 REAL(KIND=dp), DIMENSION(3) :: com, rcur, ri, rj 2342 2343 nsize = SIZE(CYCLE) 2344 2345 ! traverse the path 2346 moit(:) = 0.0_dp 2347 com(:) = helium_com(helium) 2348 DO ia = 1, nsize 2349 ! contributions from all bead pairs of the current atom 2350 DO ib = 1, helium%beads - 1 2351 ri = pos(:, CYCLE(ia), ib) - com(:) 2352 rj = pos(:, CYCLE(ia), ib + 1) - com(:) 2353 rcur(1) = ri(2)*rj(2) + ri(3)*rj(3) 2354 rcur(2) = ri(3)*rj(3) + ri(1)*rj(1) 2355 rcur(3) = ri(1)*rj(1) + ri(2)*rj(2) 2356 moit(:) = moit(:) + rcur(:) 2357 END DO 2358 ! contribution from the last bead of the current atom 2359 ! and the first bead of the next atom 2360 i1 = CYCLE(ia) 2361 IF (ia .EQ. nsize) THEN 2362 i2 = CYCLE(1) 2363 ELSE 2364 i2 = CYCLE(ia + 1) 2365 END IF 2366 ! rotation invariant bead index 2367 ri = pos(:, i1, helium%beads) - com(:) 2368 rj = pos(:, i2, 1) - com(:) 2369 rcur(1) = ri(2)*rj(2) + ri(3)*rj(3) 2370 rcur(2) = ri(3)*rj(3) + ri(1)*rj(1) 2371 rcur(3) = ri(1)*rj(1) + ri(2)*rj(2) 2372 moit(:) = moit(:) + rcur(:) 2373 END DO 2374 moit(:) = moit(:)/helium%beads 2375 2376 RETURN 2377 END FUNCTION helium_cycle_moment_of_inertia 2378 2379! *************************************************************************** 2380!> \brief Return total moment of inertia (sum over all cycles) 2381!> \param helium ... 2382!> \return ... 2383!> \date 2014-04-24 2384!> \author Lukasz Walewski 2385!> \note mostly for sanity checks 2386! ************************************************************************************************** 2387 FUNCTION helium_total_moment_of_inertia_cyclewise(helium) RESULT(moit) 2388 2389 TYPE(helium_solvent_type), POINTER :: helium 2390 REAL(KIND=dp), DIMENSION(3) :: moit 2391 2392 CHARACTER(len=*), PARAMETER :: routineN = 'helium_total_moment_of_inertia_cyclewise', & 2393 routineP = moduleN//':'//routineN 2394 2395 INTEGER :: ic 2396 REAL(KIND=dp), DIMENSION(3) :: pa 2397 TYPE(int_arr_ptr), DIMENSION(:), POINTER :: cycles 2398 2399! decompose the current permutation state into permutation cycles 2400 2401 NULLIFY (cycles) 2402 cycles => helium_calc_cycles(helium%permutation) 2403 2404 moit(:) = 0.0_dp 2405 DO ic = 1, SIZE(cycles) 2406 pa = helium_cycle_moment_of_inertia(helium, cycles(ic)%iap, helium%pos) 2407 moit(:) = moit(:) + pa(:) 2408 END DO 2409 2410 DEALLOCATE (cycles) 2411 2412 RETURN 2413 END FUNCTION helium_total_moment_of_inertia_cyclewise 2414 2415! *************************************************************************** 2416!> \brief Set coordinate system, e.g. for RHO calculations 2417!> \param helium ... 2418!> \param pint_env ... 2419!> \date 2014-04-25 2420!> \author Lukasz Walewski 2421!> \note Sets the origin (center of the coordinate system) wrt which 2422!> spatial distribution functions are calculated. 2423!> \note It can be extended later to set the axes of the coordinate system 2424!> as well, e.g. for dynamic analysis with moving solute. 2425! ************************************************************************************************** 2426 SUBROUTINE helium_update_coord_system(helium, pint_env) 2427 2428 TYPE(helium_solvent_type), POINTER :: helium 2429 TYPE(pint_env_type), POINTER :: pint_env 2430 2431 CHARACTER(len=*), PARAMETER :: routineN = 'helium_update_coord_system', & 2432 routineP = moduleN//':'//routineN 2433 2434 IF (helium%solute_present) THEN 2435 helium%center(:) = pint_com_pos(pint_env) 2436 ELSE 2437 IF (helium%periodic) THEN 2438 helium%center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/) 2439 ELSE 2440 helium%center(:) = helium_com(helium) 2441 END IF 2442 END IF 2443 2444 RETURN 2445 END SUBROUTINE helium_update_coord_system 2446 2447! *************************************************************************** 2448!> \brief Set coordinate system for RDF calculations 2449!> \param helium ... 2450!> \param pint_env ... 2451!> \date 2014-04-25 2452!> \par History 2453!> 2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran] 2454!> \author Lukasz Walewski 2455!> \note Sets the centers wrt which radial distribution functions are 2456!> calculated. 2457! ************************************************************************************************** 2458 SUBROUTINE helium_set_rdf_coord_system(helium, pint_env) 2459 2460 TYPE(helium_solvent_type), POINTER :: helium 2461 TYPE(pint_env_type), POINTER :: pint_env 2462 2463 CHARACTER(len=*), PARAMETER :: routineN = 'helium_set_rdf_coord_system', & 2464 routineP = moduleN//':'//routineN 2465 2466 INTEGER :: i, j 2467 2468 IF (helium%solute_present .AND. helium%rdf_sol_he) THEN 2469 ! Account for unequal number of beads for solute and helium 2470 DO i = 1, helium%beads 2471 j = ((i - 1)*helium%solute_beads)/helium%beads + 1 2472 helium%rdf_centers(i, :) = pint_env%x(j, :) 2473 END DO 2474 END IF 2475 2476 RETURN 2477 END SUBROUTINE helium_set_rdf_coord_system 2478 2479! *************************************************************************** 2480!> \brief Calculate center of mass 2481!> \param helium ... 2482!> \return ... 2483!> \date 2014-04-09 2484!> \author Lukasz Walewski 2485! ************************************************************************************************** 2486 FUNCTION helium_com(helium) RESULT(com) 2487 2488 TYPE(helium_solvent_type), POINTER :: helium 2489 REAL(KIND=dp), DIMENSION(3) :: com 2490 2491 CHARACTER(len=*), PARAMETER :: routineN = 'helium_com', routineP = moduleN//':'//routineN 2492 2493 INTEGER :: ia, ib 2494 2495 com(:) = 0.0_dp 2496 DO ia = 1, helium%atoms 2497 DO ib = 1, helium%beads 2498 com(:) = com(:) + helium%pos(:, ia, ib) 2499 END DO 2500 END DO 2501 com(:) = com(:)/helium%atoms/helium%beads 2502 2503 END FUNCTION helium_com 2504 2505! *************************************************************************** 2506!> \brief Return link vector, i.e. displacement vector of two consecutive 2507!> beads along the path starting at bead ib of atom ia 2508!> \param helium ... 2509!> \param ia ... 2510!> \param ib ... 2511!> \return ... 2512!> \author Lukasz Walewski 2513! ************************************************************************************************** 2514 FUNCTION helium_link_vector(helium, ia, ib) RESULT(lvec) 2515 2516 TYPE(helium_solvent_type), POINTER :: helium 2517 INTEGER :: ia, ib 2518 REAL(KIND=dp), DIMENSION(3) :: lvec 2519 2520 CHARACTER(len=*), PARAMETER :: routineN = 'helium_link_vector', & 2521 routineP = moduleN//':'//routineN 2522 2523 INTEGER :: ia1, ia2, ib1, ib2 2524 REAL(KIND=dp), DIMENSION(:), POINTER :: r1, r2 2525 2526 IF (ib .EQ. helium%beads) THEN 2527 ia1 = ia 2528 ia2 = helium%permutation(ia) 2529 ib1 = ib 2530 ib2 = 1 2531 ELSE 2532 ia1 = ia 2533 ia2 = ia 2534 ib1 = ib 2535 ib2 = ib + 1 2536 END IF 2537 r1 => helium%pos(:, ia1, ib1) 2538 r2 => helium%pos(:, ia2, ib2) 2539 lvec(:) = r2(:) - r1(:) 2540 CALL helium_pbc(helium, lvec) 2541 2542 END FUNCTION helium_link_vector 2543 2544! ************************************************************************************************** 2545!> \brief ... 2546!> \param helium ... 2547!> \param ia ... 2548!> \param ib ... 2549!> \return ... 2550! ************************************************************************************************** 2551 FUNCTION helium_rperpendicular(helium, ia, ib) RESULT(rperp) 2552 2553 TYPE(helium_solvent_type), POINTER :: helium 2554 INTEGER :: ia, ib 2555 REAL(KIND=dp), DIMENSION(3) :: rperp 2556 2557 CHARACTER(len=*), PARAMETER :: routineN = 'helium_rperpendicular', & 2558 routineP = moduleN//':'//routineN 2559 2560 REAL(KIND=dp), DIMENSION(:), POINTER :: ri 2561 2562 ri => helium%pos(:, ia, ib) 2563 rperp(1) = SQRT(ri(2)*ri(2) + ri(3)*ri(3)) 2564 rperp(2) = SQRT(ri(3)*ri(3) + ri(1)*ri(1)) 2565 rperp(3) = SQRT(ri(1)*ri(1) + ri(2)*ri(2)) 2566 2567 RETURN 2568 END FUNCTION helium_rperpendicular 2569 2570! *************************************************************************** 2571!> \brief Convert a winding number vector from real number representation 2572!> (in internal units) to integer number representation (in cell 2573!> vector units) 2574!> \param helium ... 2575!> \param wnum ... 2576!> \return ... 2577!> \author Lukasz Walewski 2578! ************************************************************************************************** 2579 FUNCTION helium_wnumber_to_integer(helium, wnum) RESULT(inum) 2580 2581 TYPE(helium_solvent_type), POINTER :: helium 2582 REAL(KIND=dp), DIMENSION(3) :: wnum 2583 INTEGER, DIMENSION(3) :: inum 2584 2585 CHARACTER(len=*), PARAMETER :: routineN = 'helium_wnumber_to_integer', & 2586 routineP = moduleN//':'//routineN 2587 2588 REAL(KIND=dp), DIMENSION(3) :: wcur 2589 REAL(KIND=dp), DIMENSION(:, :), POINTER :: invcell 2590 2591 invcell => helium%cell_m_inv 2592 CALL DGEMV('N', 3, 3, 1.0_dp, invcell, SIZE(invcell, 1), wnum, 1, 0.0_dp, wcur, 1) 2593 inum(:) = NINT(wcur(:)) 2594 2595 RETURN 2596 END FUNCTION helium_wnumber_to_integer 2597 2598! *************************************************************************** 2599!> \brief Given the atom index and permutation state returns .TRUE. if the 2600!> atom belongs to a winding path, .FASLE. otherwise. 2601!> \param helium ... 2602!> \param atmidx ... 2603!> \param pos ... 2604!> \param permutation ... 2605!> \return ... 2606!> \date 2010-09-21 2607!> \author Lukasz Walewski 2608!> \note The path winds around the periodic box if any component of its 2609!> widing number vector differs from 0. 2610! ************************************************************************************************** 2611 FUNCTION helium_is_winding(helium, atmidx, pos, permutation) RESULT(is_winding) 2612 2613 TYPE(helium_solvent_type), POINTER :: helium 2614 INTEGER, INTENT(IN) :: atmidx 2615 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pos 2616 INTEGER, DIMENSION(:), POINTER :: permutation 2617 LOGICAL :: is_winding 2618 2619 CHARACTER(len=*), PARAMETER :: routineN = 'helium_is_winding', & 2620 routineP = moduleN//':'//routineN 2621 2622 INTEGER :: ic 2623 INTEGER, DIMENSION(3) :: inum 2624 INTEGER, DIMENSION(:), POINTER :: CYCLE 2625 REAL(KIND=dp), DIMENSION(3) :: wnum 2626 2627 is_winding = .FALSE. 2628 NULLIFY (CYCLE) 2629 CYCLE => helium_cycle_of(atmidx, permutation) 2630 wnum(:) = bohr*helium_cycle_winding_number(helium, CYCLE, pos) 2631 inum(:) = helium_wnumber_to_integer(helium, wnum) 2632 DO ic = 1, 3 2633 IF (ABS(inum(ic)) .GT. 0) THEN 2634 is_winding = .TRUE. 2635 EXIT 2636 END IF 2637 END DO 2638 DEALLOCATE (CYCLE) 2639 2640 RETURN 2641 END FUNCTION helium_is_winding 2642 2643END MODULE helium_common 2644