Line data Source code
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 106955608 : 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 106955608 : IF (.NOT. (helium%periodic .OR. PRESENT(enforce))) RETURN
90 :
91 101673513 : cell_size = helium%cell_size
92 101673513 : cell_size_inv = helium%cell_size_inv
93 :
94 101673513 : rx = r(1)*cell_size_inv
95 101673513 : IF (rx > 0.5_dp) THEN
96 4098339 : rx = rx - REAL(INT(rx + 0.5_dp), dp)
97 97575174 : ELSEIF (rx < -0.5_dp) THEN
98 4159268 : rx = rx - REAL(INT(rx - 0.5_dp), dp)
99 : END IF
100 :
101 101673513 : ry = r(2)*cell_size_inv
102 101673513 : IF (ry > 0.5_dp) THEN
103 3960982 : ry = ry - REAL(INT(ry + 0.5_dp), dp)
104 97712531 : ELSEIF (ry < -0.5_dp) THEN
105 8203072 : ry = ry - REAL(INT(ry - 0.5_dp), dp)
106 : END IF
107 :
108 101673513 : rz = r(3)*cell_size_inv
109 101673513 : IF (rz > 0.5_dp) THEN
110 4236780 : rz = rz - REAL(INT(rz + 0.5_dp), dp)
111 97436733 : ELSEIF (rz < -0.5_dp) THEN
112 10562463 : rz = rz - REAL(INT(rz - 0.5_dp), dp)
113 : END IF
114 :
115 101673513 : IF (helium%cell_shape == helium_cell_shape_octahedron) THEN
116 101673513 : corr = 0.0_dp
117 101673513 : IF (rx > 0.0_dp) THEN
118 50015863 : corr = corr + rx
119 50015863 : sx = 0.5_dp
120 : ELSE
121 51657650 : corr = corr - rx
122 51657650 : sx = -0.5_dp
123 : END IF
124 101673513 : IF (ry > 0.0_dp) THEN
125 48620754 : corr = corr + ry
126 48620754 : sy = 0.5_dp
127 : ELSE
128 53052759 : corr = corr - ry
129 53052759 : sy = -0.5_dp
130 : END IF
131 101673513 : IF (rz > 0.0_dp) THEN
132 50109946 : corr = corr + rz
133 50109946 : sz = 0.5_dp
134 : ELSE
135 51563567 : corr = corr - rz
136 51563567 : sz = -0.5_dp
137 : END IF
138 101673513 : IF (corr > 0.75_dp) THEN
139 36034635 : rx = rx - sx
140 36034635 : ry = ry - sy
141 36034635 : rz = rz - sz
142 : END IF
143 : END IF
144 :
145 101673513 : r(1) = rx*cell_size
146 101673513 : r(2) = ry*cell_size
147 101673513 : 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 367358409 : 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 367358409 : IF (helium%periodic) THEN
172 348193214 : cell_size_inv = helium%cell_size_inv
173 348193214 : rx = ABS(r(1)*cell_size_inv)
174 348193214 : rx = ABS(rx - REAL(INT(rx + 0.5_dp), dp))
175 348193214 : ry = ABS(r(2)*cell_size_inv)
176 348193214 : ry = ABS(ry - REAL(INT(ry + 0.5_dp), dp))
177 348193214 : rz = ABS(r(3)*cell_size_inv)
178 348193214 : rz = ABS(rz - REAL(INT(rz + 0.5_dp), dp))
179 348193214 : IF (helium%cell_shape == helium_cell_shape_octahedron) THEN
180 348193214 : corr = rx + ry + rz - 0.75_dp
181 348193214 : IF (corr < 0.0_dp) corr = 0.0_dp
182 348193214 : helium_pbc_r2 = helium%cell_size**2*(rx*rx + ry*ry + rz*rz - corr)
183 : ELSE
184 0 : helium_pbc_r2 = helium%cell_size**2*(rx*rx + ry*ry + rz*rz)
185 : END IF
186 : ELSE
187 19165195 : helium_pbc_r2 = (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
188 : END IF
189 367358409 : 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 2959734 : 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 11838936 : c(:) = b(:) - a(:)
211 2959734 : CALL helium_pbc(helium, c)
212 11838936 : c(:) = a(:) + 0.5_dp*c(:)
213 2959734 : CALL helium_pbc(helium, c)
214 2959734 : 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 12 : 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 12 : INTEGER, DIMENSION(:), POINTER :: atoms_in_cycle
233 : LOGICAL :: atoms_left, path_end_reached
234 12 : LOGICAL, DIMENSION(:), POINTER :: atom_was_used
235 :
236 12 : CALL timeset(routineN, handle)
237 :
238 12 : NULLIFY (atoms_in_cycle)
239 12 : atoms_in_cycle => helium%itmp_atoms_1d
240 396 : atoms_in_cycle(:) = 0
241 :
242 12 : NULLIFY (atom_was_used)
243 12 : atom_was_used => helium%ltmp_atoms_1d
244 396 : atom_was_used(:) = .FALSE.
245 :
246 396 : helium%atom_plength(:) = 0
247 :
248 : start_idx = 1
249 : DO
250 384 : clen = 0
251 384 : path_end_reached = .FALSE.
252 384 : curr_idx = start_idx
253 384 : DO ia = 1, helium%atoms
254 384 : clen = clen + 1
255 384 : atoms_in_cycle(clen) = curr_idx
256 384 : atom_was_used(curr_idx) = .TRUE.
257 :
258 : ! follow to the next atom in the cycle
259 384 : curr_idx = helium%permutation(curr_idx)
260 384 : IF (curr_idx .EQ. start_idx) THEN
261 : path_end_reached = .TRUE.
262 : EXIT
263 : END IF
264 : END DO
265 384 : err_str = "Permutation path is not cyclic."
266 384 : IF (.NOT. path_end_reached) THEN
267 0 : CPABORT(err_str)
268 : END IF
269 :
270 : ! assign the cycle length to the collected atoms
271 768 : DO ia = 1, clen
272 768 : helium%atom_plength(atoms_in_cycle(ia)) = clen
273 : END DO
274 :
275 : ! go to the next "unused" atom
276 384 : atoms_left = .FALSE.
277 6720 : DO ia = 1, helium%atoms
278 6720 : 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 384 : IF (.NOT. atoms_left) EXIT
286 : END DO
287 :
288 396 : atoms_in_cycle(:) = 0
289 396 : NULLIFY (atoms_in_cycle)
290 396 : atom_was_used(:) = .FALSE.
291 12 : NULLIFY (atom_was_used)
292 :
293 12 : CALL timestop(handle)
294 :
295 12 : 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 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: cur_cycle, used_indices
315 0 : TYPE(int_arr_ptr), DIMENSION(:), POINTER :: my_cycles
316 :
317 0 : nsize = SIZE(permutation)
318 :
319 : ! most pesimistic scenario: no cycles longer than 1
320 0 : ALLOCATE (my_cycles(nsize))
321 :
322 : ! go over the permutation and identify cycles
323 0 : curat = 1
324 0 : nused = 0
325 0 : ncycl = 0
326 0 : DO WHILE (curat .LE. nsize)
327 :
328 : ! get the permutation cycle the current atom belongs to
329 0 : cur_cycle => helium_cycle_of(curat, permutation)
330 :
331 : ! include the current cycle in the pool of "used" indices
332 0 : nused = nused + SIZE(cur_cycle)
333 0 : CALL reallocate(used_indices, 1, nused)
334 0 : used_indices(nused - SIZE(cur_cycle) + 1:nused) = cur_cycle(1:SIZE(cur_cycle))
335 :
336 : ! store the pointer to the current cycle
337 0 : ncycl = ncycl + 1
338 0 : my_cycles(ncycl)%iap => cur_cycle
339 :
340 : ! free the local pointer
341 0 : NULLIFY (cur_cycle)
342 :
343 : ! try to increment the current atom index
344 0 : DO WHILE (ANY(used_indices .EQ. curat))
345 0 : curat = curat + 1
346 : END DO
347 :
348 : END DO
349 :
350 0 : DEALLOCATE (used_indices)
351 : NULLIFY (used_indices)
352 :
353 : ! assign the result
354 0 : ALLOCATE (cycles(ncycl))
355 0 : cycles(1:ncycl) = my_cycles(1:ncycl)
356 :
357 0 : DEALLOCATE (my_cycles)
358 0 : NULLIFY (my_cycles)
359 :
360 0 : 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 12 : 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 12 : 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 12 : TYPE(int_arr_ptr), DIMENSION(:), POINTER :: cycles
384 :
385 12 : CALL timeset(routineN, handle)
386 :
387 48 : maxr_half(:) = helium%rho_maxr/2.0_dp
388 12 : invdr = 1.0_dp/helium%rho_delr
389 12 : invd3r = invdr*invdr*invdr
390 12 : invp = 1.0_dp/REAL(helium%beads, dp)
391 12 : nbin = helium%rho_nbin
392 : ! assign the cycle lengths to all the atoms
393 12 : CALL helium_calc_atom_cycle_length(helium)
394 : NULLIFY (cycl_len)
395 12 : cycl_len => helium%atom_plength
396 72 : DO ir = 1, rho_num ! loop over densities ---
397 :
398 60 : 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 12 : SELECT CASE (ir)
404 :
405 : CASE (rho_atom_number)
406 12 : ii = helium%rho_property(ir)%component_index(1)
407 9912 : helium%rho_incr(ii, :, :) = invp
408 :
409 : CASE (rho_projected_area)
410 0 : vtotal(:) = helium_total_projected_area(helium)
411 0 : DO ia = 1, helium%atoms
412 0 : DO ib = 1, helium%beads
413 0 : vlink(:) = helium_link_projected_area(helium, ia, ib)
414 0 : DO ic = 1, 3
415 0 : ii = helium%rho_property(ir)%component_index(ic)
416 0 : 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 0 : vtotal(:) = helium_total_winding_number(helium)
435 0 : DO id = 1, 3
436 0 : ii = helium%rho_property(ir)%component_index(id)
437 0 : helium%rho_incr(ii, :, :) = 0.0_dp
438 : END DO
439 0 : NULLIFY (cycles)
440 0 : cycles => helium_calc_cycles(helium%permutation)
441 0 : DO ic = 1, SIZE(cycles)
442 0 : wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos)
443 0 : DO ia = 1, SIZE(cycles(ic)%iap)
444 0 : aa = cycles(ic)%iap(ia)
445 0 : DO ib = 1, helium%beads
446 0 : vlink(:) = helium_link_winding_number(helium, aa, ib)
447 0 : DO id = 1, 3
448 0 : IF (ABS(wn(id)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
449 0 : ii = helium%rho_property(ir)%component_index(id)
450 0 : 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 0 : DEALLOCATE (cycles)
457 :
458 : CASE (rho_winding_cycle)
459 0 : vtotal(:) = helium_total_winding_number(helium)
460 0 : DO id = 1, 3
461 0 : ii = helium%rho_property(ir)%component_index(id)
462 0 : helium%rho_incr(ii, :, :) = 0.0_dp
463 : END DO
464 0 : NULLIFY (cycles)
465 0 : cycles => helium_calc_cycles(helium%permutation)
466 : ! compute number of atoms in all winding cycles
467 0 : nw(:) = 0
468 0 : DO ic = 1, SIZE(cycles)
469 0 : wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos)
470 0 : DO id = 1, 3
471 0 : IF (ABS(wn(id)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
472 0 : 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 0 : DO ic = 1, SIZE(cycles)
478 0 : wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos)
479 0 : DO id = 1, 3
480 0 : IF (ABS(wn(id)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
481 0 : DO ia = 1, SIZE(cycles(ic)%iap)
482 0 : aa = cycles(ic)%iap(ia)
483 0 : DO ib = 1, helium%beads
484 0 : IF (nw(id) .GT. 0) THEN ! this test should always get passed
485 0 : ii = helium%rho_property(ir)%component_index(id)
486 0 : rtmp = invp/REAL(nw(id), dp)
487 0 : 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 0 : DEALLOCATE (cycles)
495 :
496 : CASE (rho_moment_of_inertia)
497 0 : vtotal(:) = helium_total_moment_of_inertia(helium)
498 12 : DO ia = 1, helium%atoms
499 0 : DO ib = 1, helium%beads
500 0 : vlink(:) = helium_link_moment_of_inertia(helium, ia, ib)
501 0 : DO ic = 1, 3
502 0 : ii = helium%rho_property(ir)%component_index(ic)
503 0 : 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 12 : n_out_of_range = 0
516 25332 : helium%rho_inst(:, :, :, :) = 0.0_dp
517 396 : DO ia = 1, helium%atoms
518 : ! bin the bead positions of the current atom using the increments set above
519 9996 : DO ib = 1, helium%beads
520 : ! map the current bead position to the corresponding voxel
521 38400 : r(:) = helium%pos(:, ia, ib) - helium%center(:)
522 : ! enforce PBC even if this is a non-periodic calc to avoid density leakage
523 9600 : CALL helium_pbc(helium, r, enforce=.TRUE.)
524 : ! set up bin indices (translate by L/2 to avoid non-positive array indices)
525 9600 : bx = INT((r(1) + maxr_half(1))*invdr) + 1
526 9600 : by = INT((r(2) + maxr_half(2))*invdr) + 1
527 9600 : bz = INT((r(3) + maxr_half(3))*invdr) + 1
528 : ! check that the resulting bin numbers are within array bounds
529 9600 : ltmp1 = (0 .LT. bx) .AND. (bx .LE. nbin)
530 9600 : ltmp2 = (0 .LT. by) .AND. (by .LE. nbin)
531 9600 : ltmp3 = (0 .LT. bz) .AND. (bz .LE. nbin)
532 9984 : 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 19200 : DO ic = 1, helium%rho_num_act
536 19200 : 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 0 : n_out_of_range = n_out_of_range + 1
540 : END IF
541 : END DO
542 : END DO
543 : ! scale by volume element
544 25332 : 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 12 : WRITE (msgstr, *) n_out_of_range
549 12 : msgstr = "Number of bead positions out of range: "//TRIM(ADJUSTL(msgstr))
550 12 : IF (n_out_of_range .GT. 0) THEN
551 0 : CPABORT(msgstr)
552 : END IF
553 :
554 12 : CALL timestop(handle)
555 :
556 12 : 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 12 : 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 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pref
640 : REAL(KIND=dp), DIMENSION(3) :: r, r0
641 :
642 12 : CALL timeset(routineN, handle)
643 :
644 12 : invdr = 1.0_dp/helium%rdf_delr
645 12 : nbin = helium%rdf_nbin
646 12 : n_out_of_range = 0
647 6012 : helium%rdf_inst(:, :) = 0.0_dp
648 :
649 12 : ind_hehe = 0
650 : ! Calculate He-He RDF
651 12 : IF (helium%rdf_he_he) THEN
652 12 : ind_hehe = 1
653 312 : DO ib = 1, helium%beads
654 9912 : DO ia = 1, helium%atoms
655 38400 : r0(:) = helium%pos(:, ia, ib)
656 :
657 317100 : DO ic = 1, helium%atoms
658 307200 : IF (ia == ic) CYCLE
659 1190400 : r(:) = helium%pos(:, ic, ib) - r0(:)
660 297600 : CALL helium_pbc(helium, r)
661 297600 : ri = SQRT(r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
662 297600 : bin = INT(ri*invdr) + 1
663 307200 : IF ((0 .LT. bin) .AND. (bin .LE. nbin)) THEN
664 : ! increment the RDF value for He atoms inside the r_6 sphere
665 194384 : helium%rdf_inst(ind_hehe, bin) = helium%rdf_inst(ind_hehe, bin) + 1.0_dp
666 : ELSE
667 103216 : 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 12 : IF (helium%solute_present .AND. helium%rdf_sol_he) THEN
676 0 : DO ib = 1, helium%beads
677 0 : DO ia = 1, helium%solute_atoms
678 0 : r0(:) = helium%rdf_centers(ib, 3*(ia - 1) + 1:3*(ia - 1) + 3)
679 :
680 0 : DO ic = 1, helium%atoms
681 0 : r(:) = helium%pos(:, ic, ib) - r0(:)
682 0 : CALL helium_pbc(helium, r)
683 0 : ri = SQRT(r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
684 0 : bin = INT(ri*invdr) + 1
685 0 : IF ((0 .LT. bin) .AND. (bin .LE. nbin)) THEN
686 : ! increment the RDF value for He atoms inside the r_6 sphere
687 0 : helium%rdf_inst(ind_hehe + ia, bin) = helium%rdf_inst(ind_hehe + ia, bin) + 1.0_dp
688 : ELSE
689 0 : 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 12 : IF (.NOT. helium%periodic) THEN
700 0 : IF (n_out_of_range .GT. 0) THEN
701 0 : WRITE (msgstr, *) n_out_of_range
702 0 : msgstr = "Number of bead positions out of range: "//TRIM(ADJUSTL(msgstr))
703 0 : CPWARN(msgstr)
704 : END IF
705 : END IF
706 :
707 : ! normalize the histogram to get g(r)
708 36 : ALLOCATE (pref(helium%rdf_num))
709 24 : pref(:) = 0.0_dp
710 12 : IF (helium%periodic) THEN
711 : ! Use helium density for normalization
712 24 : pref(:) = helium%density*helium%beads*helium%atoms
713 : ! Correct for He-He-RDF
714 12 : IF (helium%rdf_he_he) THEN
715 12 : 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 0 : pref(:) = 0.5_dp*helium%rdf_inst(:, 1)
721 0 : DO bin = 2, helium%rdf_nbin - 1
722 0 : pref(:) = pref(:) + helium%rdf_inst(:, bin)
723 : END DO
724 0 : pref(:) = pref(:) + 0.5_dp*helium%rdf_inst(:, helium%rdf_nbin)
725 :
726 : !set integral of histogram to number of atoms:
727 0 : pref(:) = pref(:)/helium%atoms
728 : ! Correct for He-He-RDF
729 0 : IF (helium%rdf_he_he) THEN
730 0 : pref(1) = pref(1)*helium%atoms/MAX(helium%atoms - 1, 1)
731 : END IF
732 : END IF
733 : ! Volume integral first:
734 3012 : DO bin = 1, helium%rdf_nbin
735 3000 : rlower = REAL(bin - 1, dp)*helium%rdf_delr
736 3000 : rupper = rlower + helium%rdf_delr
737 3000 : nideal = (rupper**3 - rlower**3)*4.0_dp*pi/3.0_dp
738 6012 : helium%rdf_inst(:, bin) = helium%rdf_inst(:, bin)/nideal
739 : END DO
740 : ! No normalization for density
741 24 : pref(:) = 1.0_dp/pref(:)
742 24 : DO ia = 1, helium%rdf_num
743 3024 : helium%rdf_inst(ia, :) = helium%rdf_inst(ia, :)*pref(ia)
744 : END DO
745 :
746 12 : DEALLOCATE (pref)
747 :
748 12 : CALL timestop(handle)
749 :
750 12 : 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 7047 : PURE SUBROUTINE helium_calc_plength(helium)
763 :
764 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
765 :
766 : INTEGER :: i, j, k
767 :
768 221151 : helium%plength_inst(:) = 0.0_dp
769 221151 : DO i = 1, helium%atoms
770 214104 : j = helium%permutation(i)
771 214104 : k = 1
772 0 : DO
773 214104 : IF (j == i) EXIT
774 0 : k = k + 1
775 0 : j = helium%permutation(j)
776 : END DO
777 221151 : helium%plength_inst(k) = helium%plength_inst(k) + 1
778 : END DO
779 221151 : helium%plength_inst(:) = helium%plength_inst(:)/helium%atoms
780 :
781 7047 : 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 7146 : 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 7146 : b = helium%beads
802 7146 : n = helium%atoms
803 7146 : i = MOD(nslices, b)
804 7146 : IF (i < 0) i = i + b
805 7146 : IF ((i >= b) .OR. (i < 1)) RETURN
806 6750 : helium%relrot = MOD(helium%relrot + i, b)
807 75215 : DO k = 1, i
808 8470671 : helium%work(:, :, k) = helium%pos(:, :, k)
809 : END DO
810 75850 : DO k = i + 1, b
811 8514570 : helium%pos(:, :, k - i) = helium%pos(:, :, k)
812 : END DO
813 75215 : DO k = 1, i
814 2174079 : DO j = 1, n
815 8463921 : 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 73517397 : 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 73517397 : POINTER :: offdiagsplines
845 : TYPE(spline_data_type), POINTER :: endpspline
846 :
847 73517397 : endpspline => helium%u0
848 73517397 : offdiagsplines => helium%uoffdiag
849 73517397 : nocut = .TRUE.
850 73517397 : IF (PRESENT(action)) THEN
851 0 : IF (.NOT. action) THEN
852 0 : endpspline => helium%e0
853 0 : offdiagsplines => helium%eoffdiag
854 0 : nocut = .FALSE.
855 : END IF
856 : END IF
857 :
858 73517397 : ar = SQRT(helium_pbc_r2(helium, r))
859 73517397 : arp = SQRT(helium_pbc_r2(helium, rp))
860 :
861 73517397 : IF (helium%periodic .AND. ((ar > 0.5_dp*helium%cell_size) &
862 : .OR. (arp > 0.5_dp*helium%cell_size))) THEN
863 6754137 : v = 0.0_dp
864 6754137 : IF (arp > 0.5_dp*helium%cell_size) THEN
865 4649159 : IF (nocut) v = v + helium_spline(endpspline, 0.5_dp*helium%cell_size)
866 : ELSE
867 2104978 : v = v + helium_spline(endpspline, arp)
868 : END IF
869 6754137 : IF (ar > 0.5_dp*helium%cell_size) THEN
870 4649772 : IF (nocut) v = v + helium_spline(endpspline, 0.5_dp*helium%cell_size)
871 : ELSE
872 2104365 : v = v + helium_spline(endpspline, ar)
873 : END IF
874 6754137 : res = 0.5_dp*v
875 : ELSE
876 : ! end-point action (first term):
877 66763260 : v = 0.5_dp*(helium_spline(endpspline, ar) + helium_spline(endpspline, arp))
878 267053040 : s = helium_pbc_r2(helium, r - rp)
879 66763260 : q = 0.5_dp*(ar + arp)
880 66763260 : z = (ar - arp)**2
881 66763260 : nsp = ((helium%pdx + 2)*(helium%pdx + 1))/2 - 1
882 66763260 : a = endpspline%x1
883 66763260 : b = endpspline%invh
884 66763260 : IF (q < a) THEN
885 0 : b = b*(q - a)
886 0 : a = 1.0_dp - b
887 : work(1:nsp) = a*offdiagsplines(1:nsp, 1, 1) + b*(offdiagsplines(1:nsp, 1, 2) - &
888 0 : offdiagsplines(1:nsp, 2, 2)*endpspline%h26)
889 66763260 : ELSE IF (q > endpspline%xn) THEN
890 764446 : b = b*(q - endpspline%xn) + 1.0_dp
891 764446 : a = 1.0_dp - b
892 764446 : tline = endpspline%n
893 : work(1:nsp) = b*offdiagsplines(1:nsp, 1, tline) + a*(offdiagsplines(1:nsp, 1, tline - 1) - &
894 4586676 : offdiagsplines(1:nsp, 2, tline - 1)*endpspline%h26)
895 : ELSE
896 65998814 : a = (a - q)*b
897 65998814 : tline = INT(1.0_dp - a)
898 65998814 : a = a + REAL(tline, kind=dp)
899 65998814 : b = 1.0_dp - a
900 65998814 : 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 395992884 : offdiagsplines(1:nsp, 2, tline + 1))
904 : END IF
905 66763260 : work(nsp + 1) = v
906 66763260 : tline = 1
907 66763260 : v = work(1)
908 200289780 : DO i = 1, helium%pdx
909 133526520 : v = v*z
910 133526520 : tline = tline + 1
911 133526520 : b = work(tline)
912 333816300 : DO j = 1, i
913 200289780 : tline = tline + 1
914 333816300 : b = b*s + work(tline)
915 : END DO
916 200289780 : v = v + b
917 : END DO
918 : res = v
919 : END IF
920 73517397 : 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 5057417 : 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 5057417 : POINTER :: offdiagsplines
951 : TYPE(spline_data_type), POINTER :: endpspline
952 :
953 5057417 : endpspline => helium%u0
954 5057417 : offdiagsplines => helium%uoffdiag
955 5057417 : nocut = .TRUE.
956 5057417 : IF (PRESENT(energy)) THEN
957 3289662 : IF (energy) THEN
958 3289662 : endpspline => helium%e0
959 3289662 : offdiagsplines => helium%eoffdiag
960 3289662 : nocut = .FALSE.
961 : END IF
962 : END IF
963 :
964 5057417 : ndim = helium%pdx
965 5057417 : nsp = ((ndim + 2)*(ndim + 1))/2 - 1
966 35401919 : vcoef(1:nsp + 1) = 0.0_dp
967 84366303 : totalv = 0.0_dp
968 84366303 : DO i = 1, nchain
969 84366303 : aij(i) = SQRT(helium_pbc_r2(helium, rchain(:, i)))
970 : END DO
971 79308886 : DO i = 2, nchain
972 297005876 : rchain(:, i - 1) = rchain(:, i - 1) - rchain(:, i)
973 79308886 : rchain(1, i - 1) = helium_pbc_r2(helium, rchain(:, i - 1))
974 : END DO
975 : !
976 5057417 : IF (helium%periodic) THEN
977 3319852 : ch = 0.5_dp*helium%cell_size
978 3319852 : IF (nocut) THEN
979 159247 : DO i = 2, nchain - 1
980 159247 : totalv = totalv + helium_spline(endpspline, MIN(aij(i), ch))
981 : END DO
982 41292 : totalv = totalv + 0.5_dp*helium_spline(endpspline, MIN(aij(1), ch))
983 41292 : totalv = totalv + 0.5_dp*helium_spline(endpspline, MIN(aij(nchain), ch))
984 : ELSE
985 67679200 : DO i = 2, nchain - 1
986 64400640 : ar = aij(i)
987 67679200 : IF (ar < ch) totalv = totalv + helium_spline(endpspline, ar)
988 : END DO
989 3278560 : ar = aij(1)
990 3278560 : IF (ar < ch) totalv = totalv + 0.5_dp*helium_spline(endpspline, ar)
991 3278560 : ar = aij(nchain)
992 3278560 : IF (ar < ch) totalv = totalv + 0.5_dp*helium_spline(endpspline, ar)
993 : END IF
994 : ELSE
995 6413022 : DO i = 2, nchain - 1
996 6413022 : totalv = totalv + helium_spline(endpspline, aij(i))
997 : END DO
998 1737565 : totalv = totalv + 0.5_dp*helium_spline(endpspline, aij(1))
999 1737565 : totalv = totalv + 0.5_dp*helium_spline(endpspline, aij(nchain))
1000 : END IF
1001 :
1002 79308886 : DO chainpos = 1, nchain - 1
1003 74251469 : ar = aij(chainpos)
1004 74251469 : arp = aij(chainpos + 1)
1005 74251469 : IF (helium%periodic .AND. ((ar > 0.5_dp*helium%cell_size) &
1006 5057417 : .OR. (arp > 0.5_dp*helium%cell_size))) THEN
1007 : CYCLE
1008 : ELSE
1009 67951527 : q = 0.5_dp*(ar + arp)
1010 67951527 : s = rchain(1, chainpos)
1011 67951527 : z = (ar - arp)**2
1012 67951527 : a = endpspline%x1
1013 67951527 : b = endpspline%invh
1014 67951527 : IF (q < a) THEN
1015 678 : b = b*(q - a)
1016 678 : a = 1.0_dp - b
1017 : vcoef(1:nsp) = a*offdiagsplines(1:nsp, 1, 1) + b*(offdiagsplines(1:nsp, 1, 2) - &
1018 4068 : offdiagsplines(1:nsp, 2, 2)*endpspline%h26)
1019 67950849 : ELSE IF (q > endpspline%xn) THEN
1020 5555970 : b = b*(q - endpspline%xn) + 1.0_dp
1021 5555970 : a = 1.0_dp - b
1022 5555970 : tline = endpspline%n
1023 : vcoef(1:nsp) = b*offdiagsplines(1:nsp, 1, tline) + a*(offdiagsplines(1:nsp, 1, tline - 1) - &
1024 33335820 : offdiagsplines(1:nsp, 2, tline - 1)*endpspline%h26)
1025 : ELSE
1026 62394879 : a = (a - q)*b
1027 62394879 : tline = INT(1.0_dp - a)
1028 62394879 : a = a + REAL(tline, kind=dp)
1029 62394879 : b = 1.0_dp - a
1030 62394879 : 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 374369274 : offdiagsplines(1:nsp, 2, tline + 1))
1034 : END IF
1035 67951527 : tline = 1
1036 67951527 : v = vcoef(1)
1037 203854581 : DO i = 1, ndim
1038 135903054 : tline = tline + 1
1039 135903054 : b = vcoef(tline)
1040 339757635 : DO j = 1, i
1041 203854581 : tline = tline + 1
1042 339757635 : b = b*s + vcoef(tline)
1043 : END DO
1044 203854581 : v = v*z + b
1045 : END DO
1046 67951527 : totalv = totalv + v
1047 : END IF
1048 : END DO
1049 5057417 : res = totalv
1050 5057417 : END FUNCTION helium_eval_chain
1051 :
1052 : ! **************************************************************************************************
1053 : !> \brief ...
1054 : !> \param helium ...
1055 : !> \author Harald Forbert
1056 : ! **************************************************************************************************
1057 6600 : 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 6600 : INTEGER, ALLOCATABLE, DIMENSION(:) :: lens, order
1063 6600 : INTEGER, DIMENSION(:), POINTER :: perm
1064 6600 : INTEGER, DIMENSION(:, :), POINTER :: nmatrix
1065 : REAL(KIND=dp) :: f, q, t, v
1066 6600 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: p
1067 : REAL(KIND=dp), DIMENSION(3) :: r
1068 6600 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ipmatrix, pmatrix, tmatrix
1069 6600 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pos
1070 :
1071 6600 : nb = helium%atoms
1072 : !TODO: check allocation status
1073 19800 : ALLOCATE (p(2*nb))
1074 19800 : ALLOCATE (order(nb))
1075 19800 : ALLOCATE (lens(2*nb))
1076 6600 : b = helium%beads - helium%bisection + 1
1077 6600 : f = -0.5_dp/(helium%hb2m*helium%tau*helium%bisection)
1078 6600 : tmatrix => helium%tmatrix
1079 6600 : pmatrix => helium%pmatrix
1080 6600 : ipmatrix => helium%ipmatrix
1081 6600 : nmatrix => helium%nmatrix
1082 6600 : perm => helium%permutation
1083 6600 : pos => helium%pos
1084 217800 : DO i = 1, nb
1085 6969600 : DO j = 1, nb
1086 6758400 : v = 0.0_dp
1087 27033600 : r(:) = pos(:, i, b) - pos(:, j, 1)
1088 6758400 : CALL helium_pbc(helium, r)
1089 6758400 : v = v + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
1090 6969600 : pmatrix(i, j) = f*v
1091 : END DO
1092 211200 : t = pmatrix(i, perm(i)) ! just some reference
1093 211200 : v = 0.0_dp
1094 6969600 : DO j = 1, nb
1095 6758400 : tmatrix(i, j) = EXP(pmatrix(i, j) - t)
1096 6969600 : v = v + tmatrix(i, j)
1097 : END DO
1098 : ! normalize
1099 : q = t + LOG(v)
1100 211200 : t = 1.0_dp/v
1101 6969600 : DO j = 1, nb
1102 6758400 : tmatrix(i, j) = tmatrix(i, j)*t
1103 6969600 : 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 6969600 : DO j = 1, nb
1146 6758400 : order(j) = j
1147 6969600 : p(j) = tmatrix(i, j)
1148 : END DO
1149 211200 : IF (nb > 1) THEN ! if nb = 1 it is already sorted.
1150 211200 : k = nb/2 + 1
1151 211200 : c = nb
1152 9715200 : DO
1153 9926400 : IF (k > 1) THEN
1154 : ! building up the heap:
1155 3379200 : k = k - 1
1156 3379200 : n = order(k)
1157 3379200 : v = p(k)
1158 : ELSE
1159 : ! removing the top of the heap
1160 6547200 : n = order(c)
1161 6547200 : v = p(c)
1162 6547200 : order(c) = order(1)
1163 6547200 : p(c) = p(1)
1164 6547200 : c = c - 1
1165 6547200 : IF (c == 1) THEN
1166 211200 : order(1) = n
1167 211200 : p(1) = v
1168 : EXIT
1169 : END IF
1170 : END IF
1171 9715200 : m = k
1172 9715200 : j = 2*k
1173 : ! restoring heap order between k and c
1174 21788037 : DO
1175 31503237 : IF (j > c) EXIT
1176 24171708 : IF (j < c) THEN
1177 22911836 : IF (p(j) < p(j + 1)) j = j + 1
1178 : END IF
1179 24171708 : IF (v >= p(j)) EXIT
1180 21788037 : order(m) = order(j)
1181 21788037 : p(m) = p(j)
1182 21788037 : m = j
1183 24171708 : j = 2*j
1184 : END DO
1185 9715200 : order(m) = n
1186 9715200 : 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 6969600 : DO j = nb + 1, 2*nb
1200 6969600 : 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 6758400 : DO m = nb + 1, 2*nb - 1
1209 : ! get next smallest element
1210 6547200 : IF (p(j) < p(c + 1)) THEN
1211 2319018 : v = p(j)
1212 2319018 : lens(j) = m
1213 2319018 : j = j + 1
1214 : ELSE
1215 4228182 : v = p(c + 1)
1216 4228182 : lens(c) = m
1217 4228182 : c = c + 1
1218 : END IF
1219 : ! get the second next smallest element
1220 6758400 : IF (p(j) < p(c + 1)) THEN
1221 4439382 : p(m + 1) = v + p(j)
1222 4439382 : lens(j) = m
1223 4439382 : j = j + 1
1224 : ELSE
1225 2107818 : p(m + 1) = v + p(c + 1)
1226 2107818 : lens(c) = m
1227 2107818 : 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 211200 : lens(2*nb - 1) = 0
1236 13305600 : DO m = 2*nb - 2, 1, -1
1237 13305600 : 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 6969600 : DO j = 1, nb
1253 6758400 : lens(j + nb) = 0
1254 6969600 : 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 13305600 : DO k = 1, 2*nb - 2
1265 13094400 : j = nb + 1 + (k - 1)/2
1266 13094400 : IF (lens(c) > lens(m + 1)) THEN
1267 6758400 : nmatrix(i, k) = -order(c)
1268 6758400 : lens(j + 1) = lens(c) - 1
1269 6758400 : v = p(c)
1270 6758400 : c = c + 1
1271 : ELSE
1272 6336000 : nmatrix(i, k) = m - nb
1273 6336000 : lens(j + 1) = lens(m + 1) - 1
1274 6336000 : v = p(m)
1275 6336000 : m = m + 1
1276 : END IF
1277 13094400 : p(j) = p(j) + v
1278 13305600 : 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 211200 : p(nb - 1) = 0.0_dp
1294 6765000 : DO m = nb - 1, 1, -1
1295 : ! if right child is a node, set its offset and
1296 : ! change its decision value
1297 6547200 : IF (nmatrix(i, 2*m) > 0) THEN
1298 564170 : p(nmatrix(i, 2*m)) = tmatrix(i, m)
1299 564170 : 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 6758400 : IF (nmatrix(i, 2*m - 1) > 0) THEN
1304 5771830 : p(nmatrix(i, 2*m - 1)) = p(m)
1305 5771830 : 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 6600 : helium%pweight = 0.0_dp
1399 6600 : t = helium%rng_stream_uniform%next()
1400 6600 : helium%ptable(1) = 1 + INT(t*nb)
1401 6600 : helium%ptable(2) = -1
1402 :
1403 : ! recalculate inverse permutation table (just in case)
1404 217800 : DO i = 1, nb
1405 217800 : helium%iperm(perm(i)) = i
1406 : END DO
1407 :
1408 : ! clean up:
1409 6600 : DEALLOCATE (lens)
1410 6600 : DEALLOCATE (order)
1411 6600 : DEALLOCATE (p)
1412 :
1413 6600 : END SUBROUTINE helium_update_transition_matrix
1414 :
1415 : ! **************************************************************************************************
1416 : !> \brief ...
1417 : !> \param spl ...
1418 : !> \param xx ...
1419 : !> \return ...
1420 : !> \author Harald Forbert
1421 : ! **************************************************************************************************
1422 342551603 : 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 342551603 : IF (xx < spl%x1) THEN
1430 1995 : b = spl%invh*(xx - spl%x1)
1431 1995 : a = 1.0_dp - b
1432 1995 : res = a*spl%y(1) + b*(spl%y(2) - spl%y2(2)*spl%h26)
1433 342549608 : ELSE IF (xx > spl%xn) THEN
1434 8988646 : b = spl%invh*(xx - spl%xn) + 1.0_dp
1435 8988646 : a = 1.0_dp - b
1436 8988646 : res = b*spl%y(spl%n) + a*(spl%y(spl%n - 1) - spl%y2(spl%n - 1)*spl%h26)
1437 : ELSE
1438 333560962 : res = spline_value(spl, xx)
1439 : END IF
1440 342551603 : 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 0 : 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 0 : dx = helium%pos(1, ia, ib) - helium%pos(1, ja, jb)
1463 0 : dy = helium%pos(2, ia, ib) - helium%pos(2, ja, jb)
1464 0 : dz = helium%pos(3, ia, ib) - helium%pos(3, ja, jb)
1465 0 : rij = SQRT(dx*dx + dy*dy + dz*dz)
1466 :
1467 0 : 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 250 : 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 250 : INTEGER, DIMENSION(:), POINTER :: cycle_index
1494 : LOGICAL :: found, new_cycle
1495 :
1496 250 : NULLIFY (cycle_index)
1497 250 : cycle_index => helium%itmp_atoms_1d
1498 1500 : cycle_index(:) = 0
1499 :
1500 250 : num_cycles = 0
1501 250 : found = .FALSE.
1502 250 : cycle_num = -1
1503 1000 : 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 950 : IF (found) THEN
1509 : EXIT
1510 : END IF
1511 :
1512 : ! initialize current cycle index with the current atom
1513 750 : cycle_idx = ia
1514 :
1515 750 : atom_idx = ia
1516 1000 : 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 750 : atom_idx = permutation(atom_idx)
1522 :
1523 750 : 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 1750 : DO ic = 1, num_cycles
1529 1750 : IF (cycle_index(ic) .EQ. cycle_idx) THEN
1530 0 : new_cycle = .FALSE.
1531 : END IF
1532 : END DO
1533 :
1534 750 : IF (new_cycle) THEN
1535 : ! increase number of cycles and update the current cycle's index
1536 750 : num_cycles = num_cycles + 1
1537 750 : cycle_index(num_cycles) = cycle_idx
1538 : END IF
1539 :
1540 : ! if this was the atom in question
1541 750 : IF (ia .EQ. atom_number) THEN
1542 : ! save the cycle index it belongs to
1543 250 : cycle_num = cycle_idx
1544 :
1545 : ! exit the loop over atoms, we've found what we've been looking for
1546 250 : 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 0 : 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 250 : IF (.NOT. found) THEN
1564 0 : 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 250 : cycle_number = 0
1573 750 : DO ic = 1, num_cycles
1574 750 : IF (cycle_index(ic) .EQ. cycle_num) THEN
1575 : cycle_number = ic
1576 : EXIT
1577 : END IF
1578 : END DO
1579 :
1580 250 : NULLIFY (cycle_index)
1581 :
1582 250 : 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 250 : 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 250 : atom_idx = atom_number
1605 250 : path_length = 0
1606 250 : path_end_reached = .FALSE.
1607 250 : DO ia = 1, helium%atoms
1608 250 : path_length = path_length + 1
1609 250 : atom_idx = permutation(atom_idx)
1610 250 : IF (atom_idx .EQ. atom_number) THEN
1611 : path_end_reached = .TRUE.
1612 : EXIT
1613 : END IF
1614 : END DO
1615 :
1616 250 : IF (.NOT. path_end_reached) THEN
1617 0 : path_length = -1
1618 : END IF
1619 :
1620 250 : 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 250 : 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 250 : INTEGER, DIMENSION(:), POINTER :: my_cycle
1640 : LOGICAL :: cycle_end_reached
1641 :
1642 250 : nsize = SIZE(permutation)
1643 :
1644 : ! maximum possible cycle length is the number of atoms
1645 250 : NULLIFY (my_cycle)
1646 750 : ALLOCATE (my_cycle(nsize))
1647 :
1648 : ! traverse the permutation table
1649 250 : len = 0
1650 250 : icur = element
1651 250 : cycle_end_reached = .FALSE.
1652 250 : DO ia = 1, nsize
1653 250 : len = len + 1
1654 250 : my_cycle(len) = icur
1655 250 : icur = permutation(icur)
1656 250 : IF (icur .EQ. element) THEN
1657 : cycle_end_reached = .TRUE.
1658 : EXIT
1659 : END IF
1660 : END DO
1661 :
1662 250 : IF (.NOT. cycle_end_reached) THEN
1663 : ! return null
1664 0 : NULLIFY (CYCLE)
1665 : ELSE
1666 : ! assign the result
1667 750 : ALLOCATE (CYCLE(len))
1668 1000 : CYCLE(1:len) = my_cycle(1:len)
1669 : END IF
1670 :
1671 : ! clean up
1672 250 : DEALLOCATE (my_cycle)
1673 250 : NULLIFY (my_cycle)
1674 :
1675 250 : 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 153 : 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 153 : REAL(KIND=dp), DIMENSION(:), POINTER :: ri, rj
1692 :
1693 612 : wnum(:) = 0.0_dp
1694 3369 : DO ia = 1, helium%atoms
1695 : ! sum of contributions from the rest of bead pairs
1696 64128 : DO ib = 1, helium%beads - 1
1697 60912 : ri => helium%pos(:, ia, ib)
1698 60912 : rj => helium%pos(:, ia, ib + 1)
1699 243648 : rcur(:) = ri(:) - rj(:)
1700 60912 : CALL helium_pbc(helium, rcur)
1701 246864 : wnum(:) = wnum(:) + rcur(:)
1702 : END DO
1703 : ! contribution from the last and the first bead
1704 3216 : ri => helium%pos(:, ia, helium%beads)
1705 3216 : rj => helium%pos(:, helium%permutation(ia), 1)
1706 12864 : rcur(:) = ri(:) - rj(:)
1707 3216 : CALL helium_pbc(helium, rcur)
1708 13017 : wnum(:) = wnum(:) + rcur(:)
1709 : END DO
1710 :
1711 153 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: r1, r2
1730 :
1731 0 : IF (ib .EQ. helium%beads) THEN
1732 0 : ja1 = ia
1733 0 : ja2 = helium%permutation(ia)
1734 0 : jb1 = ib
1735 0 : jb2 = 1
1736 : ELSE
1737 0 : ja1 = ia
1738 0 : ja2 = ia
1739 0 : jb1 = ib
1740 0 : jb2 = ib + 1
1741 : END IF
1742 0 : r1 => helium%pos(:, ja1, jb1)
1743 0 : r2 => helium%pos(:, ja2, jb2)
1744 0 : wnum(:) = r1(:) - r2(:)
1745 0 : CALL helium_pbc(helium, wnum)
1746 :
1747 0 : 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 0 : 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 0 : wnum(:) = 0.0_dp
1765 0 : DO ia = 1, helium%atoms
1766 0 : DO ib = 1, helium%beads
1767 0 : wnum(:) = wnum(:) + helium_link_winding_number(helium, ia, ib)
1768 : END DO
1769 : END DO
1770 :
1771 0 : 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 500 : 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 250 : REAL(KIND=dp), DIMENSION(:), POINTER :: ri, rj
1792 :
1793 250 : nsize = SIZE(CYCLE)
1794 :
1795 : ! traverse the path
1796 1000 : wnum(:) = 0.0_dp
1797 500 : DO ia = 1, nsize
1798 : ! contributions from all bead pairs of the current atom
1799 4000 : DO ib = 1, helium%beads - 1
1800 3750 : ri => pos(:, CYCLE(ia), ib)
1801 3750 : rj => pos(:, CYCLE(ia), ib + 1)
1802 15000 : rcur(:) = ri(:) - rj(:)
1803 3750 : CALL helium_pbc(helium, rcur)
1804 15250 : 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 250 : i1 = CYCLE(ia)
1809 250 : IF (ia .EQ. nsize) THEN
1810 250 : i2 = CYCLE(1)
1811 : ELSE
1812 0 : i2 = CYCLE(ia + 1)
1813 : END IF
1814 250 : ri => pos(:, i1, helium%beads)
1815 250 : rj => pos(:, i2, 1)
1816 1000 : rcur(:) = ri(:) - rj(:)
1817 250 : CALL helium_pbc(helium, rcur)
1818 1250 : wnum(:) = wnum(:) + rcur(:)
1819 : END DO
1820 :
1821 250 : 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 0 : 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 0 : TYPE(int_arr_ptr), DIMENSION(:), POINTER :: cycles
1839 :
1840 : ! decompose the current permutation state into permutation cycles
1841 :
1842 0 : NULLIFY (cycles)
1843 0 : cycles => helium_calc_cycles(helium%permutation)
1844 :
1845 0 : wnum(:) = 0.0_dp
1846 0 : DO ic = 1, SIZE(cycles)
1847 0 : wn = helium_cycle_winding_number(helium, cycles(ic)%iap, helium%pos)
1848 0 : wnum(:) = wnum(:) + wn(:)
1849 : END DO
1850 :
1851 0 : DEALLOCATE (cycles)
1852 :
1853 0 : 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 153 : 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 612 : area(:) = 0.0_dp
1871 3369 : DO ia = 1, helium%atoms
1872 : ! contributions from all links of the current atom
1873 64128 : DO ib = 1, helium%beads - 1
1874 243648 : r1(:) = helium%pos(:, ia, ib)
1875 243648 : r2(:) = helium%pos(:, ia, ib + 1)
1876 : ! comment out for non-PBC version -->
1877 243648 : r12(:) = r2(:) - r1(:)
1878 60912 : CALL helium_pbc(helium, r1)
1879 60912 : CALL helium_pbc(helium, r12)
1880 243648 : r2(:) = r1(:) + r12(:)
1881 : ! comment out for non-PBC version <--
1882 60912 : rcur(1) = r1(2)*r2(3) - r1(3)*r2(2)
1883 60912 : rcur(2) = r1(3)*r2(1) - r1(1)*r2(3)
1884 60912 : rcur(3) = r1(1)*r2(2) - r1(2)*r2(1)
1885 246864 : 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 12864 : r1(:) = helium%pos(:, ia, helium%beads)
1890 12864 : r2(:) = helium%pos(:, helium%permutation(ia), 1)
1891 : ! comment out for non-PBC version -->
1892 12864 : r12(:) = r2(:) - r1(:)
1893 3216 : CALL helium_pbc(helium, r1)
1894 3216 : CALL helium_pbc(helium, r12)
1895 12864 : r2(:) = r1(:) + r12(:)
1896 : ! comment out for non-PBC version <--
1897 3216 : rcur(1) = r1(2)*r2(3) - r1(3)*r2(2)
1898 3216 : rcur(2) = r1(3)*r2(1) - r1(1)*r2(3)
1899 3216 : rcur(3) = r1(1)*r2(2) - r1(2)*r2(1)
1900 13017 : area(:) = area(:) + rcur(:)
1901 : END DO
1902 612 : area(:) = 0.5_dp*area(:)
1903 :
1904 153 : 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 0 : 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 0 : IF (ib .EQ. helium%beads) THEN
1925 0 : ja1 = ia
1926 0 : ja2 = helium%permutation(ia)
1927 0 : jb1 = ib
1928 0 : jb2 = 1
1929 : ELSE
1930 0 : ja1 = ia
1931 0 : ja2 = ia
1932 0 : jb1 = ib
1933 0 : jb2 = ib + 1
1934 : END IF
1935 0 : r1(:) = helium%pos(:, ja1, jb1)
1936 0 : r2(:) = helium%pos(:, ja2, jb2)
1937 : ! comment out for non-PBC version -->
1938 0 : r12(:) = r2(:) - r1(:)
1939 0 : CALL helium_pbc(helium, r1)
1940 0 : CALL helium_pbc(helium, r12)
1941 0 : r2(:) = r1(:) + r12(:)
1942 : ! comment out for non-PBC version <--
1943 0 : area(1) = r1(2)*r2(3) - r1(3)*r2(2)
1944 0 : area(2) = r1(3)*r2(1) - r1(1)*r2(3)
1945 0 : area(3) = r1(1)*r2(2) - r1(2)*r2(1)
1946 0 : area(:) = 0.5_dp*area(:)
1947 :
1948 0 : 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 0 : 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 0 : area(:) = 0.0_dp
1966 0 : DO ia = 1, helium%atoms
1967 0 : DO ib = 1, helium%beads
1968 0 : area(:) = area(:) + helium_link_projected_area(helium, ia, ib)
1969 : END DO
1970 : END DO
1971 :
1972 0 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: ri, rj
1991 :
1992 0 : nsize = SIZE(CYCLE)
1993 :
1994 : ! traverse the path
1995 0 : rsum(:) = 0.0_dp
1996 0 : DO ia = 1, nsize
1997 : ! contributions from all bead pairs of the current atom
1998 0 : DO ib = 1, helium%beads - 1
1999 0 : ri => helium%pos(:, CYCLE(ia), ib)
2000 0 : rj => helium%pos(:, CYCLE(ia), ib + 1)
2001 0 : rcur(1) = ri(2)*rj(3) - ri(3)*rj(2)
2002 0 : rcur(2) = ri(3)*rj(1) - ri(1)*rj(3)
2003 0 : rcur(3) = ri(1)*rj(2) - ri(2)*rj(1)
2004 0 : 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 0 : i1 = CYCLE(ia)
2009 0 : IF (ia .EQ. nsize) THEN
2010 0 : i2 = CYCLE(1)
2011 : ELSE
2012 0 : i2 = CYCLE(ia + 1)
2013 : END IF
2014 0 : ri => helium%pos(:, i1, helium%beads)
2015 0 : rj => helium%pos(:, i2, 1)
2016 0 : rcur(1) = ri(2)*rj(3) - ri(3)*rj(2)
2017 0 : rcur(2) = ri(3)*rj(1) - ri(1)*rj(3)
2018 0 : rcur(3) = ri(1)*rj(2) - ri(2)*rj(1)
2019 0 : rsum(:) = rsum(:) + rcur(:)
2020 : END DO
2021 0 : area(:) = 0.5_dp*rsum(:)
2022 :
2023 0 : 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 0 : 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 0 : nsize = SIZE(CYCLE)
2044 :
2045 : ! traverse the path
2046 0 : area(:) = 0.0_dp
2047 0 : DO ia = 1, nsize
2048 : ! contributions from all bead pairs of the current atom
2049 0 : DO ib = 1, helium%beads - 1
2050 0 : r1(:) = helium%pos(:, CYCLE(ia), ib)
2051 0 : r2(:) = helium%pos(:, CYCLE(ia), ib + 1)
2052 0 : r12(:) = r2(:) - r1(:)
2053 0 : CALL helium_pbc(helium, r1)
2054 0 : CALL helium_pbc(helium, r12)
2055 0 : r2(:) = r1(:) + r12(:)
2056 0 : rcur(1) = r1(2)*r2(3) - r1(3)*r2(2)
2057 0 : rcur(2) = r1(3)*r2(1) - r1(1)*r2(3)
2058 0 : rcur(3) = r1(1)*r2(2) - r1(2)*r2(1)
2059 0 : 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 0 : i1 = CYCLE(ia)
2064 0 : IF (ia .EQ. nsize) THEN
2065 0 : i2 = CYCLE(1)
2066 : ELSE
2067 0 : i2 = CYCLE(ia + 1)
2068 : END IF
2069 0 : r1(:) = helium%pos(:, i1, helium%beads)
2070 0 : r2(:) = helium%pos(:, i2, 1)
2071 0 : r12(:) = r2(:) - r1(:)
2072 0 : CALL helium_pbc(helium, r1)
2073 0 : CALL helium_pbc(helium, r12)
2074 0 : r2(:) = r1(:) + r12(:)
2075 0 : rcur(1) = r1(2)*r2(3) - r1(3)*r2(2)
2076 0 : rcur(2) = r1(3)*r2(1) - r1(1)*r2(3)
2077 0 : rcur(3) = r1(1)*r2(2) - r1(2)*r2(1)
2078 0 : area(:) = area(:) + rcur(:)
2079 : END DO
2080 0 : area(:) = 0.5_dp*area(:)
2081 :
2082 0 : 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 0 : 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 0 : TYPE(int_arr_ptr), DIMENSION(:), POINTER :: cycles
2100 :
2101 : ! decompose the current permutation state into permutation cycles
2102 :
2103 0 : NULLIFY (cycles)
2104 0 : cycles => helium_calc_cycles(helium%permutation)
2105 :
2106 0 : area(:) = 0.0_dp
2107 0 : DO ic = 1, SIZE(cycles)
2108 0 : pa = helium_cycle_projected_area(helium, cycles(ic)%iap)
2109 0 : area(:) = area(:) + pa(:)
2110 : END DO
2111 :
2112 0 : 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 153 : 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 612 : com(:) = helium%center(:)
2130 :
2131 612 : moit(:) = 0.0_dp
2132 3369 : DO ia = 1, helium%atoms
2133 : ! contributions from all the links of the current atom
2134 64128 : DO ib = 1, helium%beads - 1
2135 243648 : r1(:) = helium%pos(:, ia, ib) - com(:)
2136 243648 : r2(:) = helium%pos(:, ia, ib + 1) - com(:)
2137 : ! comment out for non-PBC version -->
2138 243648 : r12(:) = r2(:) - r1(:)
2139 60912 : CALL helium_pbc(helium, r1)
2140 60912 : CALL helium_pbc(helium, r12)
2141 243648 : r2(:) = r1(:) + r12(:)
2142 : ! comment out for non-PBC version <--
2143 60912 : rcur(1) = r1(2)*r2(2) + r1(3)*r2(3)
2144 60912 : rcur(2) = r1(3)*r2(3) + r1(1)*r2(1)
2145 60912 : rcur(3) = r1(1)*r2(1) + r1(2)*r2(2)
2146 246864 : 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 12864 : r1(:) = helium%pos(:, ia, helium%beads) - com(:)
2151 12864 : r2(:) = helium%pos(:, helium%permutation(ia), 1) - com(:)
2152 : ! comment out for non-PBC version -->
2153 12864 : r12(:) = r2(:) - r1(:)
2154 3216 : CALL helium_pbc(helium, r1)
2155 3216 : CALL helium_pbc(helium, r12)
2156 12864 : r2(:) = r1(:) + r12(:)
2157 : ! comment out for non-PBC version <--
2158 3216 : rcur(1) = r1(2)*r2(2) + r1(3)*r2(3)
2159 3216 : rcur(2) = r1(3)*r2(3) + r1(1)*r2(1)
2160 3216 : rcur(3) = r1(1)*r2(1) + r1(2)*r2(2)
2161 13017 : moit(:) = moit(:) + rcur(:)
2162 : END DO
2163 612 : moit(:) = moit(:)/helium%beads
2164 :
2165 153 : 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 0 : 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 0 : com(:) = helium%center(:)
2186 :
2187 0 : IF (ib .EQ. helium%beads) THEN
2188 0 : ja1 = ia
2189 0 : ja2 = helium%permutation(ia)
2190 0 : jb1 = ib
2191 0 : jb2 = 1
2192 : ELSE
2193 0 : ja1 = ia
2194 0 : ja2 = ia
2195 0 : jb1 = ib
2196 0 : jb2 = ib + 1
2197 : END IF
2198 0 : r1(:) = helium%pos(:, ja1, jb1) - com(:)
2199 0 : r2(:) = helium%pos(:, ja2, jb2) - com(:)
2200 : ! comment out for non-PBC version -->
2201 0 : r12(:) = r2(:) - r1(:)
2202 0 : CALL helium_pbc(helium, r1)
2203 0 : CALL helium_pbc(helium, r12)
2204 0 : r2(:) = r1(:) + r12(:)
2205 : ! comment out for non-PBC version <--
2206 0 : moit(1) = r1(2)*r2(2) + r1(3)*r2(3)
2207 0 : moit(2) = r1(3)*r2(3) + r1(1)*r2(1)
2208 0 : moit(3) = r1(1)*r2(1) + r1(2)*r2(2)
2209 0 : moit(:) = moit(:)/helium%beads
2210 :
2211 0 : 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 0 : 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 0 : moit(:) = 0.0_dp
2229 0 : DO ia = 1, helium%atoms
2230 0 : DO ib = 1, helium%beads
2231 0 : moit(:) = moit(:) + helium_link_moment_of_inertia(helium, ia, ib)
2232 : END DO
2233 : END DO
2234 :
2235 0 : 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 0 : 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 0 : nsize = SIZE(CYCLE)
2257 :
2258 : ! traverse the path
2259 0 : moit(:) = 0.0_dp
2260 0 : com(:) = helium_com(helium)
2261 0 : DO ia = 1, nsize
2262 : ! contributions from all bead pairs of the current atom
2263 0 : DO ib = 1, helium%beads - 1
2264 0 : ri = pos(:, CYCLE(ia), ib) - com(:)
2265 0 : rj = pos(:, CYCLE(ia), ib + 1) - com(:)
2266 0 : rcur(1) = ri(2)*rj(2) + ri(3)*rj(3)
2267 0 : rcur(2) = ri(3)*rj(3) + ri(1)*rj(1)
2268 0 : rcur(3) = ri(1)*rj(1) + ri(2)*rj(2)
2269 0 : 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 0 : i1 = CYCLE(ia)
2274 0 : IF (ia .EQ. nsize) THEN
2275 0 : i2 = CYCLE(1)
2276 : ELSE
2277 0 : i2 = CYCLE(ia + 1)
2278 : END IF
2279 : ! rotation invariant bead index
2280 0 : ri = pos(:, i1, helium%beads) - com(:)
2281 0 : rj = pos(:, i2, 1) - com(:)
2282 0 : rcur(1) = ri(2)*rj(2) + ri(3)*rj(3)
2283 0 : rcur(2) = ri(3)*rj(3) + ri(1)*rj(1)
2284 0 : rcur(3) = ri(1)*rj(1) + ri(2)*rj(2)
2285 0 : moit(:) = moit(:) + rcur(:)
2286 : END DO
2287 0 : moit(:) = moit(:)/helium%beads
2288 :
2289 0 : 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 0 : 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 0 : TYPE(int_arr_ptr), DIMENSION(:), POINTER :: cycles
2307 :
2308 : ! decompose the current permutation state into permutation cycles
2309 :
2310 0 : NULLIFY (cycles)
2311 0 : cycles => helium_calc_cycles(helium%permutation)
2312 :
2313 0 : moit(:) = 0.0_dp
2314 0 : DO ic = 1, SIZE(cycles)
2315 0 : pa = helium_cycle_moment_of_inertia(helium, cycles(ic)%iap, helium%pos)
2316 0 : moit(:) = moit(:) + pa(:)
2317 : END DO
2318 :
2319 0 : DEALLOCATE (cycles)
2320 :
2321 0 : 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 27 : 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 27 : IF (helium%solute_present) THEN
2340 68 : helium%center(:) = pint_com_pos(pint_env)
2341 : ELSE
2342 10 : IF (helium%periodic) THEN
2343 40 : helium%center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/)
2344 : ELSE
2345 0 : helium%center(:) = helium_com(helium)
2346 : END IF
2347 : END IF
2348 :
2349 27 : 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 12 : 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 12 : IF (helium%solute_present .AND. helium%rdf_sol_he) THEN
2370 : ! Account for unequal number of beads for solute and helium
2371 0 : DO i = 1, helium%beads
2372 0 : j = ((i - 1)*helium%solute_beads)/helium%beads + 1
2373 0 : helium%rdf_centers(i, :) = pint_env%x(j, :)
2374 : END DO
2375 : END IF
2376 :
2377 12 : 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 0 : 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 0 : com(:) = 0.0_dp
2394 0 : DO ia = 1, helium%atoms
2395 0 : DO ib = 1, helium%beads
2396 0 : com(:) = com(:) + helium%pos(:, ia, ib)
2397 : END DO
2398 : END DO
2399 0 : com(:) = com(:)/helium%atoms/helium%beads
2400 :
2401 0 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: r1, r2
2420 :
2421 0 : IF (ib .EQ. helium%beads) THEN
2422 0 : ia1 = ia
2423 0 : ia2 = helium%permutation(ia)
2424 0 : ib1 = ib
2425 0 : ib2 = 1
2426 : ELSE
2427 0 : ia1 = ia
2428 0 : ia2 = ia
2429 0 : ib1 = ib
2430 0 : ib2 = ib + 1
2431 : END IF
2432 0 : r1 => helium%pos(:, ia1, ib1)
2433 0 : r2 => helium%pos(:, ia2, ib2)
2434 0 : lvec(:) = r2(:) - r1(:)
2435 0 : CALL helium_pbc(helium, lvec)
2436 :
2437 0 : END FUNCTION helium_link_vector
2438 :
2439 : ! **************************************************************************************************
2440 : !> \brief ...
2441 : !> \param helium ...
2442 : !> \param ia ...
2443 : !> \param ib ...
2444 : !> \return ...
2445 : ! **************************************************************************************************
2446 0 : 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 0 : rperp(1) = SQRT(ri(2)*ri(2) + ri(3)*ri(3))
2454 0 : rperp(2) = SQRT(ri(3)*ri(3) + ri(1)*ri(1))
2455 0 : rperp(3) = SQRT(ri(1)*ri(1) + ri(2)*ri(2))
2456 : END ASSOCIATE
2457 :
2458 0 : 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 250 : 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 250 : 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 1000 : inum(:) = NINT(wcur(:))
2479 :
2480 250 : 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 250 : 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 250 : is_winding = .FALSE.
2509 : NULLIFY (CYCLE)
2510 250 : CYCLE => helium_cycle_of(atmidx, permutation)
2511 1250 : wnum(:) = bohr*helium_cycle_winding_number(helium, CYCLE, pos)
2512 250 : inum(:) = helium_wnumber_to_integer(helium, wnum)
2513 1000 : DO ic = 1, 3
2514 1000 : IF (ABS(inum(ic)) .GT. 0) THEN
2515 : is_winding = .TRUE.
2516 : EXIT
2517 : END IF
2518 : END DO
2519 250 : DEALLOCATE (CYCLE)
2520 :
2521 250 : END FUNCTION helium_is_winding
2522 :
2523 : END MODULE helium_common
|
