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