(git:1f285aa)
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 ! **************************************************************************************************
14 MODULE 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 
68 CONTAINS
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 
2523 END MODULE helium_common
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