(git:34ef472)
helium_sampling.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Methods for sampling helium variables
10 !> \author Lukasz Walewski
11 !> \date 2009-06-10
12 ! **************************************************************************************************
13 MODULE helium_sampling
14 
15  USE cp_external_control, ONLY: external_control
16  USE cp_log_handling, ONLY: cp_add_default_logger,&
17  cp_get_default_logger,&
18  cp_logger_create,&
19  cp_logger_release,&
20  cp_logger_type,&
21  cp_rm_default_logger
22  USE cp_output_handling, ONLY: cp_add_iter_level,&
23  cp_iterate,&
24  cp_rm_iter_level
25  USE global_types, ONLY: global_environment_type
26  USE helium_common, ONLY: &
27  helium_boxmean_3d, helium_calc_plength, helium_calc_rdf, helium_calc_rho, helium_com, &
28  helium_eval_expansion, helium_pbc, helium_rotate, helium_set_rdf_coord_system, &
29  helium_spline, helium_total_moment_of_inertia, helium_total_projected_area, &
30  helium_total_winding_number, helium_update_transition_matrix
31  USE helium_interactions, ONLY: helium_bead_solute_e_f,&
32  helium_calc_energy,&
33  helium_solute_e_f
34  USE helium_io, ONLY: &
35  helium_print_accepts, helium_print_action, helium_print_coordinates, helium_print_energy, &
36  helium_print_force, helium_print_perm, helium_print_plength, helium_print_rdf, &
37  helium_print_rho, helium_print_vector, helium_write_line
38  USE helium_types, ONLY: e_id_total,&
39  helium_solvent_p_type,&
40  helium_solvent_type
41  USE helium_worm, ONLY: helium_sample_worm
42  USE input_constants, ONLY: &
43  helium_forces_average, helium_forces_last, helium_mdist_exponential, &
44  helium_mdist_gaussian, helium_mdist_linear, helium_mdist_quadratic, helium_mdist_singlev, &
45  helium_mdist_uniform, helium_sampling_ceperley, helium_sampling_worm
46  USE input_cp2k_restarts, ONLY: write_restart
47  USE kinds, ONLY: default_string_length,&
48  dp
49  USE machine, ONLY: m_walltime
50  USE physcon, ONLY: angstrom
51  USE pint_public, ONLY: pint_com_pos
52  USE pint_types, ONLY: pint_env_type
53  USE splines_types, ONLY: spline_data_type
54 #include "../base/base_uses.f90"
55 
56  IMPLICIT NONE
57 
58  PRIVATE
59 
60  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
61  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_sampling'
62 
63  PUBLIC :: helium_do_run
64  PUBLIC :: helium_sample
65  PUBLIC :: helium_step
66 
67 CONTAINS
68 
69 ! ***************************************************************************
70 !> \brief Performs MC simulation for helium (only)
71 !> \param helium_env ...
72 !> \param globenv ...
73 !> \date 2009-07-14
74 !> \par History
75 !> 2016-07-14 Modified to work with independent helium_env [cschran]
76 !> \author Lukasz Walewski
77 !> \note This routine gets called only when HELIUM_ONLY is set to .TRUE.,
78 !> so do not put any property calculations here!
79 ! **************************************************************************************************
80  SUBROUTINE helium_do_run(helium_env, globenv)
81  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
82  TYPE(global_environment_type), POINTER :: globenv
83 
84  INTEGER :: k, num_steps, step, tot_steps
85  LOGICAL :: should_stop
86  TYPE(cp_logger_type), POINTER :: logger
87  TYPE(pint_env_type) :: pint_env
88 
89  NULLIFY (logger)
90  logger => cp_get_default_logger()
91 
92  num_steps = 0
93 
94  IF (ASSOCIATED(helium_env)) THEN
95  DO k = 1, SIZE(helium_env)
96 
97  NULLIFY (helium_env(k)%helium%logger)
98  helium_env(k)%helium%logger => cp_get_default_logger()
99 
100  ! create iteration level
101  ! Although the Helium MC is not really an md it is most of the times
102  ! embedded in the pint code so the same step counter applies.
103  IF (k .EQ. 1) THEN
104  CALL cp_add_iter_level(helium_env(k)%helium%logger%iter_info, "PINT")
105  CALL cp_iterate(helium_env(k)%helium%logger%iter_info, &
106  iter_nr=helium_env(k)%helium%first_step)
107  END IF
108  tot_steps = helium_env(k)%helium%first_step
109  num_steps = helium_env(k)%helium%num_steps
110 
111  ! set the properties accumulated over the whole MC process to 0
112  helium_env(k)%helium%proarea%accu(:) = 0.0_dp
113  helium_env(k)%helium%prarea2%accu(:) = 0.0_dp
114  helium_env(k)%helium%wnmber2%accu(:) = 0.0_dp
115  helium_env(k)%helium%mominer%accu(:) = 0.0_dp
116  IF (helium_env(k)%helium%rho_present) THEN
117  helium_env(k)%helium%rho_accu(:, :, :, :) = 0.0_dp
118  END IF
119  IF (helium_env(k)%helium%rdf_present) THEN
120  helium_env(k)%helium%rdf_accu(:, :) = 0.0_dp
121  END IF
122  END DO
123  END IF
124 
125  ! Distribute steps for processors without helium_env
126  CALL logger%para_env%bcast(num_steps)
127  CALL logger%para_env%bcast(tot_steps)
128 
129  DO step = 1, num_steps
130 
131  tot_steps = tot_steps + 1
132 
133  IF (ASSOCIATED(helium_env)) THEN
134  DO k = 1, SIZE(helium_env)
135  IF (k .EQ. 1) THEN
136  CALL cp_iterate(helium_env(k)%helium%logger%iter_info, &
137  last=(step == helium_env(k)%helium%num_steps), &
138  iter_nr=tot_steps)
139  END IF
140  helium_env(k)%helium%current_step = tot_steps
141  END DO
142  END IF
143 
144  CALL helium_step(helium_env, pint_env)
145 
146  ! call write_restart here to avoid interference with PINT write_restart
147  IF (ASSOCIATED(helium_env)) THEN
148  CALL write_restart(root_section=helium_env(1)%helium%input, helium_env=helium_env)
149  END IF
150 
151  ! exit from the main loop if soft exit has been requested
152  CALL external_control(should_stop, "PINT", globenv=globenv)
153  IF (should_stop) EXIT
154 
155  END DO
156 
157  ! remove iteration level
158  IF (ASSOCIATED(helium_env)) THEN
159  CALL cp_rm_iter_level(helium_env(1)%helium%logger%iter_info, "PINT")
160  END IF
161 
162  RETURN
163  END SUBROUTINE helium_do_run
164 
165 ! ***************************************************************************
166 !> \brief Sample the helium phase space
167 !> \param helium_env ...
168 !> \param pint_env ...
169 !> \date 2009-10-27
170 !> \par History
171 !> 2016-07-14 Modified to work with independent helium_env [cschran]
172 !> \author Lukasz Walewski
173 !> \note Samples helium variable space according to multilevel Metropolis
174 !> MC scheme, calculates the forces exerted by helium solvent on the
175 !> solute and stores them in helium%force_avrg array. The forces are
176 !> averaged over outer MC loop.
177 !> \note The implicit assumption is that we simulate solute _with_ solvent
178 !> most of the time, so for performance reasons I do not check that
179 !> everywhere I should. This leads to some redundancy in the case of
180 !> helium-only calculations.
181 ! **************************************************************************************************
182  SUBROUTINE helium_sample(helium_env, pint_env)
183 
184  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
185  TYPE(pint_env_type), INTENT(IN) :: pint_env
186 
187  CHARACTER(len=default_string_length) :: msg_str
188  INTEGER :: i, irot, iweight, k, nslices, nsteps, &
189  num_env, offset, sel_mp_source
190  REAL(kind=dp) :: inv_num_env, inv_xn, rnd, rtmp, rweight
191  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work_2d
192  TYPE(cp_logger_type), POINTER :: logger
193 
194  NULLIFY (logger)
195  IF (SIZE(helium_env) < 1) THEN
196  logger => cp_get_default_logger()
197  ELSE
198  CALL cp_logger_create(logger, para_env=helium_env(1)%comm, template_logger=cp_get_default_logger())
199  CALL cp_add_default_logger(logger)
200  END IF
201 
202  DO k = 1, SIZE(helium_env)
203 
204  IF (helium_env(k)%helium%solute_present) THEN
205  helium_env(k)%helium%force_avrg(:, :) = 0.0_dp
206  helium_env(k)%helium%current_step = pint_env%iter
207  helium_env(k)%helium%origin = pint_com_pos(pint_env)
208  END IF
209 
210  helium_env(k)%helium%energy_avrg(:) = 0.0_dp
211  helium_env(k)%helium%plength_avrg(:) = 0.0_dp
212  helium_env(k)%helium%num_accepted(:, :) = 0.0_dp
213 
214  ! helium parallelization: each processor gets different RN stream and
215  ! runs independent helium simulation, the properties and forces are
216  ! averaged over parallel helium environments once per step.
217  inv_xn = 0.0_dp
218  SELECT CASE (helium_env(k)%helium%sampling_method)
219 
220  CASE (helium_sampling_worm)
221 
222  CALL helium_sample_worm(helium_env(k)%helium, pint_env)
223 
224  inv_xn = 1.0_dp/real(helium_env(k)%helium%worm_nstat, dp)
225 
226  CASE (helium_sampling_ceperley)
227  ! helium sampling (outer MC loop)
228  DO irot = 1, helium_env(k)%helium%iter_rot
229 
230  ! rotate helium beads in imaginary time at random, store current
231  ! 'rotation state' in helium%relrot which is within (0, helium%beads-1)
232  ! (this is needed to sample different fragments of the permutation
233  ! paths in try_permutations)
234  rnd = helium_env(k)%helium%rng_stream_uniform%next()
235  nslices = int(rnd*helium_env(k)%helium%beads)
236  CALL helium_rotate(helium_env(k)%helium, nslices)
237 
238  CALL helium_try_permutations(helium_env(k)%helium, pint_env)
239 
240  ! calculate & accumulate instantaneous properties for averaging
241  IF (helium_env(k)%helium%solute_present) THEN
242  IF (helium_env(k)%helium%get_helium_forces == helium_forces_average) THEN
243  CALL helium_solute_e_f(pint_env, helium_env(k)%helium, rtmp)
244  helium_env(k)%helium%force_avrg(:, :) = helium_env(k)%helium%force_avrg(:, :) + &
245  helium_env(k)%helium%force_inst(:, :)
246  END IF
247  END IF
248  CALL helium_calc_energy(helium_env(k)%helium, pint_env)
249 
250  helium_env(k)%helium%energy_avrg(:) = helium_env(k)%helium%energy_avrg(:) + helium_env(k)%helium%energy_inst(:)
251  CALL helium_calc_plength(helium_env(k)%helium)
252  helium_env(k)%helium%plength_avrg(:) = helium_env(k)%helium%plength_avrg(:) + helium_env(k)%helium%plength_inst(:)
253 
254  ! instantaneous force output according to HELIUM%PRINT%FORCES_INST
255  ! Warning: file I/O here may cost A LOT of cpu time!
256  ! switched off here to save cpu
257  !CALL helium_print_force_inst( helium_env(k)%helium, error )
258 
259  END DO ! outer MC loop
260 
261  ! If only the last forces are taken, calculate them for the last outer MC loop configuration
262  IF ((helium_env(k)%helium%solute_present) .AND. (helium_env(k)%helium%get_helium_forces == helium_forces_last)) THEN
263  CALL helium_solute_e_f(pint_env, helium_env(k)%helium, rtmp)
264  helium_env(k)%helium%force_avrg(:, :) = helium_env(k)%helium%force_inst(:, :)
265  END IF
266 
267  ! restore the original alignment of beads in imaginary time
268  ! (this is useful to make a single bead's position easy to follow
269  ! in the trajectory, otherwise its index would change every step)
270  CALL helium_rotate(helium_env(k)%helium, -helium_env(k)%helium%relrot)
271 
272  inv_xn = 1.0_dp/real(helium_env(k)%helium%iter_rot, dp)
273 
274  CASE DEFAULT
275  WRITE (msg_str, *) helium_env(k)%helium%sampling_method
276  msg_str = "Sampling method ("//trim(adjustl(msg_str))//") not supported."
277  cpabort(msg_str)
278 
279  END SELECT
280 
281  ! actually average the properties over the outer MC loop
282  IF ((helium_env(k)%helium%solute_present) .AND. (helium_env(k)%helium%get_helium_forces == helium_forces_average)) THEN
283  helium_env(k)%helium%force_avrg(:, :) = helium_env(k)%helium%force_avrg(:, :)*inv_xn
284  END IF
285  helium_env(k)%helium%energy_avrg(:) = helium_env(k)%helium%energy_avrg(:)*inv_xn
286  helium_env(k)%helium%plength_avrg(:) = helium_env(k)%helium%plength_avrg(:)*inv_xn
287 
288  ! instantaneous properties
289  helium_env(k)%helium%proarea%inst(:) = helium_total_projected_area(helium_env(k)%helium)
290  helium_env(k)%helium%prarea2%inst(:) = helium_env(k)%helium%proarea%inst(:)*helium_env(k)%helium%proarea%inst(:)
291  helium_env(k)%helium%wnumber%inst(:) = helium_total_winding_number(helium_env(k)%helium)
292  helium_env(k)%helium%wnmber2%inst(:) = helium_env(k)%helium%wnumber%inst(:)*helium_env(k)%helium%wnumber%inst(:)
293  helium_env(k)%helium%mominer%inst(:) = helium_total_moment_of_inertia(helium_env(k)%helium)
294 
295  ! properties accumulated over the whole MC process
296  helium_env(k)%helium%proarea%accu(:) = helium_env(k)%helium%proarea%accu(:) + helium_env(k)%helium%proarea%inst(:)
297  helium_env(k)%helium%prarea2%accu(:) = helium_env(k)%helium%prarea2%accu(:) + helium_env(k)%helium%prarea2%inst(:)
298  helium_env(k)%helium%wnmber2%accu(:) = helium_env(k)%helium%wnmber2%accu(:) + helium_env(k)%helium%wnmber2%inst(:)
299  helium_env(k)%helium%mominer%accu(:) = helium_env(k)%helium%mominer%accu(:) + helium_env(k)%helium%mominer%inst(:)
300  IF (helium_env(k)%helium%rho_present) THEN
301  CALL helium_calc_rho(helium_env(k)%helium)
302  helium_env(k)%helium%rho_accu(:, :, :, :) = helium_env(k)%helium%rho_accu(:, :, :, :) + &
303  helium_env(k)%helium%rho_inst(:, :, :, :)
304  END IF
305  IF (helium_env(k)%helium%rdf_present) THEN
306  CALL helium_set_rdf_coord_system(helium_env(k)%helium, pint_env)
307  CALL helium_calc_rdf(helium_env(k)%helium)
308  helium_env(k)%helium%rdf_accu(:, :) = helium_env(k)%helium%rdf_accu(:, :) + helium_env(k)%helium%rdf_inst(:, :)
309  END IF
310 
311  ! running averages (restart-aware)
312  nsteps = helium_env(k)%helium%current_step - helium_env(k)%helium%first_step
313  iweight = helium_env(k)%helium%averages_iweight
314  rweight = real(iweight, dp)
315  rtmp = 1.0_dp/(real(max(1, nsteps + iweight), dp))
316  helium_env(k)%helium%proarea%ravr(:) = (helium_env(k)%helium%proarea%accu(:) + &
317  rweight*helium_env(k)%helium%proarea%rstr(:))*rtmp
318  helium_env(k)%helium%prarea2%ravr(:) = (helium_env(k)%helium%prarea2%accu(:) + &
319  rweight*helium_env(k)%helium%prarea2%rstr(:))*rtmp
320  helium_env(k)%helium%wnmber2%ravr(:) = (helium_env(k)%helium%wnmber2%accu(:) + &
321  rweight*helium_env(k)%helium%wnmber2%rstr(:))*rtmp
322  helium_env(k)%helium%mominer%ravr(:) = (helium_env(k)%helium%mominer%accu(:) + &
323  rweight*helium_env(k)%helium%mominer%rstr(:))*rtmp
324 
325  END DO
326 
327  ! average over helium environments sitting at different processors
328 !TODO the following involves message passing, maybe move it to i/o routines?
329  num_env = helium_env(1)%helium%num_env
330  inv_num_env = 1.0_dp/real(num_env, dp)
331 
332  !energy_avrg:
333  DO k = 2, SIZE(helium_env)
334  helium_env(1)%helium%energy_avrg(:) = helium_env(1)%helium%energy_avrg(:) + &
335  helium_env(k)%helium%energy_avrg(:)
336  END DO
337  CALL helium_env(1)%comm%sum(helium_env(1)%helium%energy_avrg(:))
338  helium_env(1)%helium%energy_avrg(:) = helium_env(1)%helium%energy_avrg(:)*inv_num_env
339  DO k = 2, SIZE(helium_env)
340  helium_env(k)%helium%energy_avrg(:) = helium_env(1)%helium%energy_avrg(:)
341  END DO
342 
343  !plength_avrg:
344  DO k = 2, SIZE(helium_env)
345  helium_env(1)%helium%plength_avrg(:) = helium_env(1)%helium%plength_avrg(:) + &
346  helium_env(k)%helium%plength_avrg(:)
347  END DO
348  CALL helium_env(1)%comm%sum(helium_env(1)%helium%plength_avrg(:))
349  helium_env(1)%helium%plength_avrg(:) = helium_env(1)%helium%plength_avrg(:)*inv_num_env
350  DO k = 2, SIZE(helium_env)
351  helium_env(k)%helium%plength_avrg(:) = helium_env(1)%helium%plength_avrg(:)
352  END DO
353 
354  !num_accepted:
355  DO k = 2, SIZE(helium_env)
356  helium_env(1)%helium%num_accepted(:, :) = helium_env(1)%helium%num_accepted(:, :) + &
357  helium_env(k)%helium%num_accepted(:, :)
358  END DO
359  CALL helium_env(1)%comm%sum(helium_env(1)%helium%num_accepted(:, :))
360  helium_env(1)%helium%num_accepted(:, :) = helium_env(1)%helium%num_accepted(:, :)*inv_num_env
361  DO k = 2, SIZE(helium_env)
362  helium_env(k)%helium%num_accepted(:, :) = helium_env(1)%helium%num_accepted(:, :)
363  END DO
364 
365  !force_avrg:
366  IF (helium_env(1)%helium%solute_present) THEN
367  IF (helium_env(1)%helium%get_helium_forces == helium_forces_average) THEN
368  DO k = 2, SIZE(helium_env)
369  helium_env(1)%helium%force_avrg(:, :) = helium_env(1)%helium%force_avrg(:, :) + &
370  helium_env(k)%helium%force_avrg(:, :)
371  END DO
372  CALL helium_env(1)%comm%sum(helium_env(1)%helium%force_avrg(:, :))
373  helium_env(1)%helium%force_avrg(:, :) = helium_env(1)%helium%force_avrg(:, :)*inv_num_env
374  DO k = 2, SIZE(helium_env)
375  helium_env(k)%helium%force_avrg(:, :) = helium_env(1)%helium%force_avrg(:, :)
376  END DO
377  ELSE IF (helium_env(1)%helium%get_helium_forces == helium_forces_last) THEN
378  IF (logger%para_env%is_source()) THEN
379  sel_mp_source = int(helium_env(1)%helium%rng_stream_uniform%next() &
380  *real(helium_env(1)%helium%num_env, dp))
381  END IF
382  CALL helium_env(1)%comm%bcast(sel_mp_source, logger%para_env%source)
383 
384  offset = 0
385  DO i = 1, logger%para_env%mepos
386  offset = offset + helium_env(1)%env_all(i)
387  END DO
388  ALLOCATE (work_2d(SIZE(helium_env(1)%helium%force_avrg, 1), &
389  SIZE(helium_env(1)%helium%force_avrg, 2)))
390  work_2d(:, :) = 0.0_dp
391  IF (sel_mp_source .GE. offset .AND. sel_mp_source .LT. offset + SIZE(helium_env)) THEN
392  work_2d(:, :) = helium_env(sel_mp_source - offset + 1)%helium%force_avrg(:, :)
393  END IF
394  CALL helium_env(1)%comm%sum(work_2d(:, :))
395  DO k = 1, SIZE(helium_env)
396  helium_env(k)%helium%force_avrg(:, :) = work_2d(:, :)
397  END DO
398  DEALLOCATE (work_2d)
399  END IF
400  END IF
401 
402  IF (SIZE(helium_env) > 0) THEN
403  CALL cp_rm_default_logger()
404  CALL cp_logger_release(logger)
405  END IF
406 
407  RETURN
408  END SUBROUTINE helium_sample
409 
410 ! ***************************************************************************
411 !> \brief Perform MC step for helium
412 !> \param helium_env ...
413 !> \param pint_env ...
414 !> \date 2009-06-12
415 !> \par History
416 !> 2016-07-14 Modified to work with independent helium_env [cschran]
417 !> \author Lukasz Walewski
418 ! **************************************************************************************************
419  SUBROUTINE helium_step(helium_env, pint_env)
420 
421  TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
422  TYPE(pint_env_type), INTENT(IN) :: pint_env
423 
424  CHARACTER(len=*), PARAMETER :: routinen = 'helium_step'
425 
426  CHARACTER(len=default_string_length) :: msgstr, stmp, time_unit
427  INTEGER :: handle, k
428  REAL(kind=dp) :: time_start, time_stop, time_used
429  REAL(kind=dp), DIMENSION(:), POINTER :: DATA
430 
431  CALL timeset(routinen, handle)
432  time_start = m_walltime()
433 
434  IF (ASSOCIATED(helium_env)) THEN
435  ! perform the actual phase space sampling
436  CALL helium_sample(helium_env, pint_env)
437 
438  ! print properties
439  CALL helium_print_energy(helium_env)
440  CALL helium_print_plength(helium_env)
441  CALL helium_print_accepts(helium_env)
442  CALL helium_print_force(helium_env)
443  CALL helium_print_perm(helium_env)
444  CALL helium_print_coordinates(helium_env)
445  CALL helium_print_action(pint_env, helium_env)
446 
447  IF (helium_env(1)%helium%rdf_present) CALL helium_print_rdf(helium_env)
448  IF (helium_env(1)%helium%rho_present) CALL helium_print_rho(helium_env)
449 
450  NULLIFY (data)
451  ALLOCATE (DATA(3*SIZE(helium_env)))
452 
453  WRITE (stmp, *) helium_env(1)%helium%apref
454  DATA(:) = 0.0_dp
455  DO k = 1, SIZE(helium_env)
456  DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%proarea%inst(:)
457  END DO
458  CALL helium_print_vector(helium_env, &
459  "MOTION%PINT%HELIUM%PRINT%PROJECTED_AREA", &
460  DATA, &
461  angstrom*angstrom, &
462  (/"A_x", "A_y", "A_z"/), &
463  "prefactor = "//trim(adjustl(stmp))//" [Angstrom^-2]", &
464  "helium-pa")
465 
466  DATA(:) = 0.0_dp
467  DO k = 1, SIZE(helium_env)
468  DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%prarea2%ravr(:)
469  END DO
470  CALL helium_print_vector(helium_env, &
471  "MOTION%PINT%HELIUM%PRINT%PROJECTED_AREA_2_AVG", &
472  DATA, &
473  angstrom*angstrom*angstrom*angstrom, &
474  (/"<A_x^2>", "<A_y^2>", "<A_z^2>"/), &
475  "prefactor = "//trim(adjustl(stmp))//" [Angstrom^-2]", &
476  "helium-pa-avg", &
477  "REWIND", &
478  .true.)
479 
480  DATA(:) = 0.0_dp
481  DO k = 1, SIZE(helium_env)
482  DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%mominer%inst(:)
483  END DO
484  CALL helium_print_vector(helium_env, &
485  "MOTION%PINT%HELIUM%PRINT%MOMENT_OF_INERTIA", &
486  DATA, &
487  angstrom*angstrom, &
488  (/"I_x/m", "I_y/m", "I_z/m"/), &
489  "prefactor = "//trim(adjustl(stmp))//" [Angstrom^-2]", &
490  "helium-mi")
491 
492  DATA(:) = 0.0_dp
493  DO k = 1, SIZE(helium_env)
494  DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%mominer%ravr
495  END DO
496  CALL helium_print_vector(helium_env, &
497  "MOTION%PINT%HELIUM%PRINT%MOMENT_OF_INERTIA_AVG", &
498  DATA, &
499  angstrom*angstrom, &
500  (/"<I_x/m>", "<I_y/m>", "<I_z/m>"/), &
501  "prefactor = "//trim(adjustl(stmp))//" [Angstrom^-2]", &
502  "helium-mi-avg", &
503  "REWIND", &
504  .true.)
505 
506  DATA(:) = 0.0_dp
507  DO k = 1, SIZE(helium_env)
508  DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%wnumber%inst
509  END DO
510  WRITE (stmp, *) helium_env(1)%helium%wpref
511  CALL helium_print_vector(helium_env, &
512  "MOTION%PINT%HELIUM%PRINT%WINDING_NUMBER", &
513  DATA, &
514  angstrom, &
515  (/"W_x", "W_y", "W_z"/), &
516  "prefactor = "//trim(adjustl(stmp))//" [Angstrom^-2]", &
517  "helium-wn")
518 
519  DATA(:) = 0.0_dp
520  DO k = 1, SIZE(helium_env)
521  DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%wnmber2%ravr
522  END DO
523  CALL helium_print_vector(helium_env, &
524  "MOTION%PINT%HELIUM%PRINT%WINDING_NUMBER_2_AVG", &
525  DATA, &
526  angstrom*angstrom, &
527  (/"<W_x^2>", "<W_y^2>", "<W_z^2>"/), &
528  "prefactor = "//trim(adjustl(stmp))//" [Angstrom^-2]", &
529  "helium-wn-avg", &
530  "REWIND", &
531  .true.)
532  DEALLOCATE (data)
533 
534  time_stop = m_walltime()
535  time_used = time_stop - time_start
536  time_unit = "sec"
537  IF (time_used .GE. 60.0_dp) THEN
538  time_used = time_used/60.0_dp
539  time_unit = "min"
540  IF (time_used .GE. 60.0_dp) THEN
541  time_used = time_used/60.0_dp
542  time_unit = "hours"
543  END IF
544  END IF
545  msgstr = "MC step"
546  stmp = ""
547  WRITE (stmp, *) helium_env(1)%helium%current_step
548  msgstr = trim(adjustl(msgstr))//" "//trim(adjustl(stmp))//" of"
549  stmp = ""
550  WRITE (stmp, *) helium_env(1)%helium%last_step
551  msgstr = trim(adjustl(msgstr))//" "//trim(adjustl(stmp))//" in"
552  stmp = ""
553  WRITE (stmp, '(F20.1)') time_used
554  msgstr = trim(adjustl(msgstr))//" "//trim(adjustl(stmp))
555  msgstr = trim(adjustl(msgstr))//" "//trim(adjustl(time_unit))//"."
556  CALL helium_write_line(trim(msgstr))
557 
558  ! print out the total energy - for regtest evaluation
559  stmp = ""
560  WRITE (stmp, *) helium_env(1)%helium%energy_avrg(e_id_total)
561  msgstr = "Total energy = "//trim(adjustl(stmp))
562  CALL helium_write_line(trim(msgstr))
563  END IF
564 
565  CALL timestop(handle)
566 
567  RETURN
568  END SUBROUTINE helium_step
569 
570 ! ***************************************************************************
571 !> \brief ...
572 !> \param helium ...
573 !> \param pint_env path integral environment
574 !> \par History
575 !> 2010-06-17 added different distributions for m-sampling [lwalewski]
576 !> 2010-06-17 ratio for m-value added (m-sampling related) [lwalewski]
577 !> \author hforbert
578 ! **************************************************************************************************
579  SUBROUTINE helium_try_permutations(helium, pint_env)
580  TYPE(helium_solvent_type), POINTER :: helium
581  TYPE(pint_env_type), INTENT(IN) :: pint_env
582 
583  CHARACTER(len=default_string_length) :: err_str, stmp
584  INTEGER :: cyclen, ni
585  LOGICAL :: accepted, selected
586  REAL(kind=dp) :: r, rnd, x, y, z
587 
588  IF (helium%maxcycle > 1) CALL helium_update_transition_matrix(helium)
589 
590  ! consecutive calls to helium_slice_metro_cyclic require that
591  ! ALL(helium%pos.EQ.helium%work) - see a comment there,
592  ! besides there is a magic regarding helium%permutation
593  ! you have been warned
594  helium%work(:, :, :) = helium%pos(:, :, :)
595 
596  ! the inner MC loop (without rotation in imaginary time)
597  DO ni = 1, helium%iter_norot
598 
599  ! set the probability threshold for m_value: 1/(1+(m-1)/helium%m_ratio)
600  r = 1.0_dp/(1.0_dp + (helium%maxcycle - 1)/helium%m_ratio)
601 
602  ! draw permutation length for this trial from the distribution of choice
603  !
604  SELECT CASE (helium%m_dist_type)
605 
606  CASE (helium_mdist_singlev)
607  x = helium%rng_stream_uniform%next()
608  IF (x .LT. r) THEN
609  cyclen = 1
610  ELSE
611  cyclen = helium%m_value
612  END IF
613 
614  CASE (helium_mdist_uniform)
615  x = helium%rng_stream_uniform%next()
616  IF (x .LT. r) THEN
617  cyclen = helium%m_value
618  ELSE
619  DO
620  x = helium%rng_stream_uniform%next()
621  cyclen = int(helium%maxcycle*x) + 1
622  IF (cyclen .NE. helium%m_value) EXIT
623  END DO
624  END IF
625 
626  CASE (helium_mdist_linear)
627  x = helium%rng_stream_uniform%next()
628  IF (x .LT. r) THEN
629  cyclen = helium%m_value
630  ELSE
631  DO
632  x = helium%rng_stream_uniform%next()
633  y = sqrt(2.0_dp*x)
634  cyclen = int(helium%maxcycle*y/sqrt(2.0_dp)) + 1
635  IF (cyclen .NE. helium%m_value) EXIT
636  END DO
637  END IF
638 
639  CASE (helium_mdist_quadratic)
640  x = helium%rng_stream_uniform%next()
641  IF (x .LT. r) THEN
642  cyclen = helium%m_value
643  ELSE
644  DO
645  x = helium%rng_stream_uniform%next()
646  y = (3.0_dp*x)**(1.0_dp/3.0_dp)
647  z = 3.0_dp**(1.0_dp/3.0_dp)
648  cyclen = int(helium%maxcycle*y/z) + 1
649  IF (cyclen .NE. helium%m_value) EXIT
650  END DO
651  END IF
652 
653  CASE (helium_mdist_exponential)
654  x = helium%rng_stream_uniform%next()
655  IF (x .LT. r) THEN
656  cyclen = helium%m_value
657  ELSE
658  DO
659  DO
660  x = helium%rng_stream_uniform%next()
661  IF (x .GE. 0.01_dp) EXIT
662  END DO
663  z = -log(0.01_dp)
664  y = log(x)/z + 1.0_dp;
665  cyclen = int(helium%maxcycle*y) + 1
666  IF (cyclen .NE. helium%m_value) EXIT
667  END DO
668  END IF
669 
670  CASE (helium_mdist_gaussian)
671  x = helium%rng_stream_uniform%next()
672  IF (x .LT. r) THEN
673  cyclen = 1
674  ELSE
675  DO
676  x = helium%rng_stream_gaussian%next()
677  cyclen = int(x*0.75_dp + helium%m_value - 0.5_dp) + 1
678  IF (cyclen .NE. 1) EXIT
679  END DO
680  END IF
681 
682  CASE DEFAULT
683  WRITE (stmp, *) helium%m_dist_type
684  err_str = "M distribution type unknown ("//trim(adjustl(stmp))//")"
685  cpabort(err_str)
686  cyclen = 1
687 
688  END SELECT
689 
690  IF (cyclen < 1) cyclen = 1
691  IF (cyclen > helium%maxcycle) cyclen = helium%maxcycle
692  helium%num_accepted(1, cyclen) = helium%num_accepted(1, cyclen) + 1
693 
694  ! check, if permutation of this length can be constructed
695  IF (cyclen == 1) THEN
696  rnd = helium%rng_stream_uniform%next()
697  helium%ptable(1) = 1 + int(rnd*helium%atoms)
698  helium%ptable(2) = -1
699  helium%pweight = 0.0_dp
700  selected = .true.
701  ELSE
702  selected = helium_select_permutation(helium, cyclen)
703  END IF
704 
705  IF (selected) THEN
706  ! permutation was successfully selected - actually sample this permutation
707  accepted = helium_slice_metro_cyclic(helium, pint_env, cyclen)
708  ELSE
709  accepted = .false.
710  END IF
711 
712 !if (any(helium%pos .ne. helium%work)) then
713 ! print *, ""
714 ! print *, "pos and work are different!"
715 ! print *, ""
716 !end if
717 
718  IF (accepted) THEN
719  ! positions changed
720  IF (helium%solute_present) THEN
721  ! use COM of the solute, but it did not move here so do nothing to save cpu
722  ELSE
723  IF (helium%periodic) THEN
724  ! use center of the cell, but it did not move here so do nothing to save cpu
725  ELSE
726  ! update center of mass
727  helium%center(:) = helium_com(helium)
728  END IF
729  END IF
730  END IF
731 
732  END DO
733 
734  RETURN
735  END SUBROUTINE helium_try_permutations
736 
737 ! ***************************************************************************
738 !> \brief Sample path variables for permutation of length <cyclen>
739 !> \param helium ...
740 !> \param pint_env ...
741 !> \param cyclen ...
742 !> \return ...
743 !> \author hforbert
744 ! **************************************************************************************************
745  FUNCTION helium_slice_metro_cyclic(helium, pint_env, cyclen) RESULT(res)
746  TYPE(helium_solvent_type), POINTER :: helium
747  TYPE(pint_env_type), INTENT(IN) :: pint_env
748  INTEGER, INTENT(IN) :: cyclen
749  LOGICAL :: res
750 
751  CHARACTER(len=default_string_length) :: err_str, stmp
752  INTEGER :: c, ia, ib, ic, ifix, j, k, l, level, &
753  pk1, pk2, stride
754  INTEGER, DIMENSION(:), POINTER :: p, perm
755  LOGICAL :: nperiodic, should_reject
756  REAL(kind=dp) :: cell_size, ds, dtk, e1, e2, pds, &
757  prev_ds, r, rtmp, rtmpo, sigma, x
758  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work3
759  REAL(kind=dp), DIMENSION(3) :: bis, biso, new_com, rm1, rm2, tmp1, tmp2
760  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pos, work
761  TYPE(spline_data_type), POINTER :: u0
762 
763 ! trial permutation implicit in p
764 ! since we (momentarily) only use cyclic permutations:
765 ! cyclen = 1 : no permutation, sample p[0] anew
766 ! cyclen = 2 : p[0] -> p[1] -> p[0]
767 ! cyclen = 3 : p[0] -> p[1] -> p[2] -> p[0]
768 ! cyclen = 4 : p[0] -> p[1] -> p[2] -> p[3] -> p[0]
769 
770  p => helium%ptable
771  prev_ds = helium%pweight
772 
773  helium%num_accepted(2, cyclen) = helium%num_accepted(2, cyclen) + 1
774  level = 1
775  res = .false.
776 
777  ! nothing to be done
778  IF (helium%bisection == 0) RETURN
779 
780  ! On entry helium_slice_metro_cyclic ASSUMES that work is a copy of pos
781  ! and constructs the trial move there. If the move is accepted, the
782  ! changed parts are copied to pos, but if it fails, only the CHANGED parts
783  ! of work are overwritten by the corresponding parts of pos. So on exit work
784  ! will AGAIN be a copy of pos (either way accept/reject). This is done here
785  ! so we do not have to copy around the whole pos array in the caller, and
786  ! the caller does not need to know which parts of work might have changed.
787  pos => helium%pos
788  work => helium%work
789  perm => helium%permutation
790  u0 => helium%u0
791  cell_size = (0.5_dp*helium%cell_size)**2
792  nperiodic = .NOT. helium%periodic
793 
794  pds = prev_ds
795  ifix = helium%beads - helium%bisection + 1
796 
797  ! sanity checks
798  !
799  IF (ifix < 1) THEN
800  WRITE (stmp, *) ifix
801  err_str = "ifix<1 test failed (ifix="//trim(adjustl(stmp))//")"
802  cpabort(err_str)
803  END IF
804  !
805  j = 1
806  k = helium%bisection
807  DO
808  IF (k < 2) EXIT
809  j = j*2
810  k = k/2
811  END DO
812  !
813  IF (j /= helium%bisection) THEN
814  WRITE (stmp, *) helium%bisection
815  err_str = "helium%bisection not a power of 2 (helium%bisection="//trim(adjustl(stmp))//")"
816  cpabort(err_str)
817  END IF
818  !
819  IF (helium%bisection < 2) THEN
820  WRITE (stmp, *) helium%bisection
821  err_str = "helium%bisection less than 2 (helium%bisection="//trim(adjustl(stmp))//")"
822  cpabort(err_str)
823  END IF
824 
825  stride = helium%bisection
826  DO
827  IF (stride <= 2) EXIT
828  ! calc new trial positions
829  ! trial action: 0.5*stride*endpointapprox
830  sigma = sqrt(0.5_dp*helium%hb2m*(stride/2)*helium%tau)
831  dtk = 0.0_dp
832  ds = 0.0_dp
833 
834  j = ifix + stride/2
835  DO
836  IF (j > helium%beads - stride/2) EXIT
837  pk1 = j - stride/2
838  pk2 = j + stride/2
839  ! calculate log(T(s)):
840  DO k = 1, cyclen
841  CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, p(k), pk2), bis)
842  tmp1(:) = bis(:) - pos(:, p(k), j)
843  CALL helium_pbc(helium, tmp1)
844  tmp1(:) = tmp1(:)/sigma
845  dtk = dtk - 0.5_dp*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
846  END DO
847  ! calculate log(T(sprime)) and sprime itself
848  DO k = 1, cyclen
849  CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, p(k), pk2), tmp1)
850  DO c = 1, 3
851  x = helium%rng_stream_gaussian%next(variance=1.0_dp)
852  x = sigma*x
853  tmp1(c) = tmp1(c) + x
854  tmp2(c) = x
855  END DO
856  CALL helium_pbc(helium, tmp1)
857  CALL helium_pbc(helium, tmp2)
858  work(:, p(k), j) = tmp1(:)
859  tmp2(:) = tmp2(:)/sigma
860  dtk = dtk + 0.5_dp*(tmp2(1)*tmp2(1) + tmp2(2)*tmp2(2) + tmp2(3)*tmp2(3))
861  END DO
862  j = j + stride
863  END DO
864 
865  j = helium%beads - stride/2 + 1
866  pk1 = j - stride/2
867  DO k = 1, cyclen
868  CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, perm(p(k)), 1), bis)
869  tmp1(:) = bis(:) - pos(:, p(k), j)
870  CALL helium_pbc(helium, tmp1)
871  tmp1(:) = tmp1(:)/sigma
872  dtk = dtk - 0.5_dp*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
873  END DO
874  DO k = 1, cyclen
875  CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, perm(p(1 + mod(k, cyclen))), 1), tmp1)
876  DO c = 1, 3
877  x = helium%rng_stream_gaussian%next(variance=1.0_dp)
878  x = sigma*x
879  tmp1(c) = tmp1(c) + x
880  tmp2(c) = x
881  END DO
882  CALL helium_pbc(helium, tmp1)
883  CALL helium_pbc(helium, tmp2)
884  work(:, p(k), j) = tmp1(:)
885  tmp2(:) = tmp2(:)/sigma
886  dtk = dtk + 0.5_dp*(tmp2(1)*tmp2(1) + tmp2(2)*tmp2(2) + tmp2(3)*tmp2(3))
887  END DO
888  ! ok we got the new positions
889  ! calculate action_k(s)-action_k(sprime)
890  x = 1.0_dp/(helium%tau*helium%hb2m*stride)
891  j = ifix
892  DO
893  IF (j > helium%beads - stride/2) EXIT
894  pk1 = j + stride/2
895  DO k = 1, cyclen
896  tmp1(:) = pos(:, p(k), j) - pos(:, p(k), pk1)
897  CALL helium_pbc(helium, tmp1)
898  ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
899  tmp1(:) = work(:, p(k), j) - work(:, p(k), pk1)
900  CALL helium_pbc(helium, tmp1)
901  ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
902  ! interaction change
903  IF (helium%solute_present) THEN
904  CALL helium_bead_solute_e_f(pint_env, helium, p(k), pk1, energy=e1)
905  CALL helium_bead_solute_e_f(pint_env, helium, p(k), pk1, work(:, p(k), pk1), e2)
906  ds = ds + (stride/2)*(e1 - e2)*helium%tau
907  END IF
908  DO l = 1, helium%atoms
909  IF (l /= p(k)) THEN
910  tmp1(:) = pos(:, p(k), pk1) - pos(:, l, pk1)
911  CALL helium_pbc(helium, tmp1)
912  r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
913  IF ((r < cell_size) .OR. nperiodic) THEN
914  r = sqrt(r)
915  ds = ds + real(stride/2, dp)*helium_spline(u0, r)
916  END IF
917  tmp1(:) = work(:, p(k), pk1) - work(:, l, pk1)
918  CALL helium_pbc(helium, tmp1)
919  r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
920  IF ((r < cell_size) .OR. nperiodic) THEN
921  r = sqrt(r)
922  ds = ds - real(stride/2, dp)*helium_spline(u0, r)
923  END IF
924  END IF
925  END DO
926  ! counted p[k], p[m] twice. subtract those again
927  IF (k < cyclen) THEN
928  DO l = k + 1, cyclen
929  tmp1(:) = pos(:, p(k), pk1) - pos(:, p(l), pk1)
930  CALL helium_pbc(helium, tmp1)
931  r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
932  IF ((r < cell_size) .OR. nperiodic) THEN
933  r = sqrt(r)
934  ds = ds - real(stride/2, dp)*helium_spline(u0, r)
935  END IF
936  tmp1(:) = work(:, p(k), pk1) - work(:, p(l), pk1)
937  CALL helium_pbc(helium, tmp1)
938  r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
939  IF ((r < cell_size) .OR. nperiodic) THEN
940  r = sqrt(r)
941  ds = ds + real(stride/2, dp)*helium_spline(u0, r)
942  END IF
943  END DO
944  END IF
945  END DO
946  j = j + stride/2
947  END DO
948  ! last link
949  pk1 = helium%beads - stride/2 + 1
950  DO k = 1, cyclen
951  tmp1(:) = pos(:, p(k), pk1) - pos(:, perm(p(k)), 1)
952  CALL helium_pbc(helium, tmp1)
953  ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
954  tmp1(:) = work(:, p(k), pk1) - work(:, perm(p(1 + mod(k, cyclen))), 1)
955  CALL helium_pbc(helium, tmp1)
956  ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
957  END DO
958  ! ok now accept or reject:
959  rtmp = helium%rng_stream_uniform%next()
960 ! IF ((dtk+ds-pds < 0.0_dp).AND.(EXP(dtk+ds-pds)<rtmp)) THEN
961  IF (dtk + ds - pds < 0.0_dp) THEN
962  IF (exp(dtk + ds - pds) < rtmp) THEN
963  DO c = ifix, helium%beads
964  DO k = 1, cyclen
965  work(:, p(k), c) = pos(:, p(k), c)
966  END DO
967  END DO
968  RETURN
969  END IF
970  END IF
971  ! accepted. go on to the next level
972  helium%num_accepted(level + 2, cyclen) = helium%num_accepted(level + 2, cyclen) + 1
973  level = level + 1
974  pds = ds
975  stride = stride/2
976  END DO
977  ! we are on the lowest level now
978  ! calc new trial positions
979  ! trial action: endpointapprox for T, full action otherwise
980  sigma = sqrt(0.5_dp*helium%hb2m*helium%tau)
981  dtk = 0.0_dp
982  ds = 0.0_dp
983  DO j = ifix + 1, helium%beads - 1, 2
984  pk1 = j - 1
985  pk2 = j + 1
986  ! calculate log(T(s)):
987  DO k = 1, cyclen
988  CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, p(k), pk2), bis)
989  tmp1(:) = bis(:) - pos(:, p(k), j)
990  CALL helium_pbc(helium, tmp1)
991  tmp1(:) = tmp1(:)/sigma
992  dtk = dtk - 0.5_dp*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
993  END DO
994  ! calculate log(T(sprime)) and sprime itself
995  DO k = 1, cyclen
996  CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, p(k), pk2), tmp1)
997  DO c = 1, 3
998  x = helium%rng_stream_gaussian%next(variance=1.0_dp)
999  x = sigma*x
1000  tmp1(c) = tmp1(c) + x
1001  tmp2(c) = x
1002  END DO
1003  CALL helium_pbc(helium, tmp1)
1004  CALL helium_pbc(helium, tmp2)
1005  work(:, p(k), j) = tmp1(:)
1006  tmp2(:) = tmp2(:)/sigma
1007  dtk = dtk + 0.5_dp*(tmp2(1)*tmp2(1) + tmp2(2)*tmp2(2) + tmp2(3)*tmp2(3))
1008  END DO
1009  END DO
1010  j = helium%beads
1011  pk1 = j - 1
1012  DO k = 1, cyclen
1013  CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, perm(p(k)), 1), bis)
1014  tmp1(:) = bis(:) - pos(:, p(k), j)
1015  CALL helium_pbc(helium, tmp1)
1016  tmp1(:) = tmp1(:)/sigma
1017  dtk = dtk - 0.5_dp*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
1018  END DO
1019  DO k = 1, cyclen
1020  CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, perm(p(1 + mod(k, cyclen))), 1), tmp1)
1021  DO c = 1, 3
1022  x = helium%rng_stream_gaussian%next(variance=1.0_dp)
1023  x = sigma*x
1024  tmp1(c) = tmp1(c) + x
1025  tmp2(c) = x
1026  END DO
1027  CALL helium_pbc(helium, tmp1)
1028  CALL helium_pbc(helium, tmp2)
1029  work(:, p(k), j) = tmp1(:)
1030  tmp2 = tmp2/sigma
1031  dtk = dtk + 0.5_dp*(tmp2(1)*tmp2(1) + tmp2(2)*tmp2(2) + tmp2(3)*tmp2(3))
1032  END DO
1033  ! ok we got the new positions.
1034  ! calculate action_k(s)-action_k(sprime)
1035  ! interaction change
1036 !TODO interaction ONLY here? or some simplified 12-6 in the upper part?
1037  IF (helium%solute_present) THEN
1038  DO j = ifix, helium%beads
1039  DO k = 1, cyclen
1040  CALL helium_bead_solute_e_f(pint_env, helium, p(k), j, energy=e1)
1041  CALL helium_bead_solute_e_f(pint_env, helium, p(k), j, work(:, p(k), j), e2)
1042  ds = ds + (e1 - e2)*helium%tau
1043  END DO
1044  END DO
1045  END IF
1046  ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
1047  x = 1.0_dp/(helium%tau*helium%hb2m*stride)
1048  DO j = ifix, helium%beads - 1
1049  pk1 = j + 1
1050  DO k = 1, cyclen
1051  tmp1(:) = pos(:, p(k), j) - pos(:, p(k), pk1)
1052  CALL helium_pbc(helium, tmp1)
1053  ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
1054  tmp1(:) = work(:, p(k), j) - work(:, p(k), pk1)
1055  CALL helium_pbc(helium, tmp1)
1056  ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
1057  DO l = 1, helium%atoms
1058  IF (l /= p(k)) THEN
1059  rm1(:) = pos(:, p(k), j) - pos(:, l, j)
1060  rm2(:) = pos(:, p(k), pk1) - pos(:, l, pk1)
1061  ds = ds + helium_eval_expansion(helium, rm1, rm2, work3)
1062  rm1(:) = work(:, p(k), j) - work(:, l, j)
1063  rm2(:) = work(:, p(k), pk1) - work(:, l, pk1)
1064  ds = ds - helium_eval_expansion(helium, rm1, rm2, work3)
1065  END IF
1066  END DO
1067  ! counted p[k], p[m] twice. subtract those again
1068  IF (k < cyclen) THEN
1069  DO l = k + 1, cyclen
1070  rm1(:) = pos(:, p(k), j) - pos(:, p(l), j)
1071  rm2(:) = pos(:, p(k), pk1) - pos(:, p(l), pk1)
1072  ds = ds - helium_eval_expansion(helium, rm1, rm2, work3)
1073  rm1(:) = work(:, p(k), j) - work(:, p(l), j)
1074  rm2(:) = work(:, p(k), pk1) - work(:, p(l), pk1)
1075  ds = ds + helium_eval_expansion(helium, rm1, rm2, work3)
1076  END DO
1077  END IF
1078  END DO
1079  END DO
1080  j = helium%beads
1081  pk1 = 1
1082  DO k = 1, cyclen
1083  tmp1(:) = pos(:, p(k), j) - pos(:, perm(p(k)), 1)
1084  CALL helium_pbc(helium, tmp1)
1085  ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
1086  tmp1(:) = work(:, p(k), j) - work(:, perm(p(1 + mod(k, cyclen))), 1)
1087  CALL helium_pbc(helium, tmp1)
1088  ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
1089  DO l = 1, helium%atoms
1090  IF (l /= p(k)) THEN
1091  rm1(:) = pos(:, p(k), j) - pos(:, l, j)
1092  rm2(:) = pos(:, perm(p(k)), 1) - pos(:, perm(l), 1)
1093  ds = ds + helium_eval_expansion(helium, rm1, rm2, work3)
1094  END IF
1095  END DO
1096  ! counted p[k], p[m] twice. subtract those again
1097  IF (k < cyclen) THEN
1098  DO l = k + 1, cyclen
1099  rm1(:) = pos(:, p(k), j) - pos(:, p(l), j)
1100  rm2(:) = pos(:, perm(p(k)), pk1) - pos(:, perm(p(l)), pk1)
1101  ds = ds - helium_eval_expansion(helium, rm1, rm2, work3)
1102  END DO
1103  END IF
1104  END DO
1105  IF (cyclen > 1) THEN
1106  !k,c,l
1107  c = perm(p(1))
1108  DO k = 1, cyclen - 1
1109  perm(p(k)) = perm(p(k + 1))
1110  END DO
1111  perm(p(cyclen)) = c
1112  END IF
1113  DO k = 1, cyclen
1114  DO l = 1, helium%atoms
1115  IF (l /= p(k)) THEN
1116  rm1(:) = work(:, p(k), j) - work(:, l, j)
1117  rm2(:) = work(:, perm(p(k)), 1) - work(:, perm(l), 1)
1118  ds = ds - helium_eval_expansion(helium, rm1, rm2, work3)
1119  END IF
1120  END DO
1121  ! counted p[k], p[m] twice. subtract those again
1122  IF (k < cyclen) THEN
1123  DO l = k + 1, cyclen
1124  rm1(:) = work(:, p(k), j) - work(:, p(l), j)
1125  rm2(:) = work(:, perm(p(k)), 1) - work(:, perm(p(l)), 1)
1126  ds = ds + helium_eval_expansion(helium, rm1, rm2, work3)
1127  END DO
1128  END IF
1129  END DO
1130  DEALLOCATE (work3)
1131  ! ok now accept or reject:
1132  rtmp = helium%rng_stream_uniform%next()
1133 ! IF ((dtk+ds-pds<0.0_dp).AND.(EXP(dtk+ds-pds)<rtmp)) THEN
1134  IF (dtk + ds - pds < 0.0_dp) THEN
1135  IF (exp(dtk + ds - pds) < rtmp) THEN
1136  DO c = ifix, helium%beads
1137  DO k = 1, cyclen
1138  work(:, p(k), c) = pos(:, p(k), c)
1139  END DO
1140  END DO
1141  IF (cyclen > 1) THEN
1142  c = perm(p(cyclen))
1143  DO k = cyclen - 1, 1, -1
1144  perm(p(k + 1)) = perm(p(k))
1145  END DO
1146  perm(p(1)) = c
1147  END IF
1148  RETURN
1149  END IF
1150  END IF
1151  ! accepted.
1152 
1153  ! rejection of the whole move if any centroid is farther away
1154  ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
1155  ! AND ist not moving towards the center
1156  IF (.NOT. helium%periodic) THEN
1157  IF (helium%solute_present) THEN
1158  new_com(:) = helium%center(:)
1159  ELSE
1160  new_com(:) = 0.0_dp
1161  DO k = 1, helium%atoms
1162  DO l = 1, helium%beads
1163  new_com(:) = new_com(:) + helium%work(:, k, l)
1164  END DO
1165  END DO
1166  new_com(:) = new_com(:)/helium%atoms/helium%beads
1167  END IF
1168  ! Cycle through atoms (ignore connectivity)
1169  should_reject = .false.
1170  outer: DO ia = 1, helium%atoms
1171  bis(:) = 0.0_dp
1172  DO ib = 1, helium%beads
1173  DO ic = 1, 3
1174  bis(ic) = bis(ic) + work(ic, ia, ib) - new_com(ic)
1175  END DO
1176  END DO
1177  bis(:) = bis(:)/helium%beads
1178  rtmp = dot_product(bis, bis)
1179  IF (rtmp .GE. helium%droplet_radius**2) THEN
1180  biso(:) = 0.0_dp
1181  DO ib = 1, helium%beads
1182  DO ic = 1, 3
1183  biso(ic) = biso(ic) + pos(ic, ia, ib) - helium%center(ic)
1184  END DO
1185  END DO
1186  biso(:) = biso(:)/helium%beads
1187  rtmpo = dot_product(biso, biso)
1188  ! only reject if it moves away from COM outside the droplet radius
1189  IF (rtmpo < rtmp) THEN
1190  ! found - this one does not fit within R from the center
1191  should_reject = .true.
1192  EXIT outer
1193  END IF
1194  END IF
1195  END DO outer
1196  IF (should_reject) THEN
1197  ! restore work and perm, then return
1198  DO c = ifix, helium%beads
1199  DO k = 1, cyclen
1200  work(:, p(k), c) = pos(:, p(k), c)
1201  END DO
1202  END DO
1203  IF (cyclen > 1) THEN
1204  c = perm(p(cyclen))
1205  DO k = cyclen - 1, 1, -1
1206  perm(p(k + 1)) = perm(p(k))
1207  END DO
1208  perm(p(1)) = c
1209  END IF
1210  RETURN
1211  END IF
1212  END IF
1213  ! accepted.
1214 
1215  ! copy trial over to the real thing
1216  DO c = ifix, helium%beads
1217  DO k = 1, cyclen
1218  pos(:, p(k), c) = work(:, p(k), c)
1219  END DO
1220  END DO
1221  DO k = 1, cyclen
1222  helium%iperm(perm(p(k))) = p(k)
1223  END DO
1224  helium%num_accepted(level + 2, cyclen) = helium%num_accepted(level + 2, cyclen) + 1
1225  res = .true.
1226 
1227  RETURN
1228  END FUNCTION helium_slice_metro_cyclic
1229 
1230 ! ***************************************************************************
1231 !> \brief ...
1232 !> \param helium ...
1233 !> \param len ...
1234 !> \return ...
1235 !> \author hforbert
1236 ! **************************************************************************************************
1237  FUNCTION helium_select_permutation(helium, len) RESULT(res)
1238  TYPE(helium_solvent_type), POINTER :: helium
1239  INTEGER, INTENT(IN) :: len
1240  LOGICAL :: res
1241 
1242  INTEGER :: i, j, k, n
1243  INTEGER, DIMENSION(:), POINTER :: iperm, p, perm
1244  INTEGER, DIMENSION(:, :), POINTER :: nmatrix
1245  REAL(kind=dp) :: rnd, s1, s2, t
1246  REAL(kind=dp), DIMENSION(:, :), POINTER :: ipmatrix, pmatrix, tmatrix
1247 
1248  s1 = 0.0_dp
1249  s2 = 0.0_dp
1250  res = .false.
1251  n = helium%atoms
1252  tmatrix => helium%tmatrix
1253  pmatrix => helium%pmatrix
1254  ipmatrix => helium%ipmatrix
1255  perm => helium%permutation
1256  iperm => helium%iperm
1257  p => helium%ptable
1258  nmatrix => helium%nmatrix
1259 
1260  p(len + 1) = -1
1261  rnd = helium%rng_stream_uniform%next()
1262  p(1) = int(n*rnd) + 1
1263  DO k = 1, len - 1
1264  t = helium%rng_stream_uniform%next()
1265  ! find the corresponding path to connect to
1266  ! using the precalculated optimal decision tree:
1267  i = n - 1
1268  DO
1269  IF (tmatrix(p(k), i) > t) THEN
1270  i = nmatrix(p(k), 2*i - 1)
1271  ELSE
1272  i = nmatrix(p(k), 2*i)
1273  END IF
1274  IF (i < 0) EXIT
1275  END DO
1276  i = -i
1277  ! which particle was it previously connected to?
1278  p(k + 1) = iperm(i)
1279  ! is it unique? quit if it was already part of the permutation
1280  DO j = 1, k
1281  IF (p(j) == p(k + 1)) RETURN
1282  END DO
1283  ! acummulate the needed values for the final
1284  ! accept/reject step:
1285  s1 = s1 + ipmatrix(p(k), i)
1286  s2 = s2 + ipmatrix(p(k), perm(p(k)))
1287  END DO
1288  ! close the permutation loop:
1289  s1 = s1 + ipmatrix(p(len), perm(p(1)))
1290  s2 = s2 + ipmatrix(p(len), perm(p(len)))
1291  ! final accept/reject:
1292  rnd = helium%rng_stream_uniform%next()
1293  t = s1*rnd
1294  IF (t > s2) RETURN
1295  ! ok, we have accepted the permutation
1296  ! calculate the action bias for the subsequent resampling
1297  ! of the paths:
1298  s1 = pmatrix(p(len), perm(p(1))) - pmatrix(p(len), perm(p(len)))
1299  DO k = 1, len - 1
1300  s1 = s1 + pmatrix(p(k), perm(p(k + 1))) - pmatrix(p(k), perm(p(k)))
1301  END DO
1302  helium%pweight = s1
1303  res = .true.
1304  RETURN
1305  END FUNCTION helium_select_permutation
1306 
1307 END MODULE helium_sampling
cp_external_control
Routines to handle the external control of CP2K.
Definition: cp_external_control.F:15
cp_external_control::external_control
subroutine, public external_control(should_stop, flag, globenv, target_time, start_time, force_check)
External manipulations during a run : when the <PROJECT_NAME>.EXIT_$runtype command is sent the progr...
Definition: cp_external_control.F:90
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_rm_default_logger
subroutine, public cp_rm_default_logger()
the cousin of cp_add_default_logger, decrements the stack, so that the default logger is what it has ...
Definition: cp_log_handling.F:213
cp_log_handling::cp_logger_release
subroutine, public cp_logger_release(logger)
releases this logger
Definition: cp_log_handling.F:423
cp_log_handling::cp_logger_create
subroutine, public cp_logger_create(logger, para_env, print_level, default_global_unit_nr, default_local_unit_nr, global_filename, local_filename, close_global_unit_on_dealloc, iter_info, close_local_unit_on_dealloc, suffix, template_logger)
initializes a logger
Definition: cp_log_handling.F:285
cp_log_handling::cp_add_default_logger
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
Definition: cp_log_handling.F:188
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_iterate
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
Definition: cp_output_handling.F:574
cp_output_handling::cp_rm_iter_level
subroutine, public cp_rm_iter_level(iteration_info, level_name, n_rlevel_att)
Removes an iteration level.
Definition: cp_output_handling.F:655
cp_output_handling::cp_add_iter_level
subroutine, public cp_add_iter_level(iteration_info, level_name, n_rlevel_new)
Adds an iteration level.
Definition: cp_output_handling.F:612
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition: global_types.F:21
helium_common
Independent helium subroutines shared with other modules.
Definition: helium_common.F:14
helium_common::helium_total_moment_of_inertia
real(kind=dp) function, dimension(3), public helium_total_moment_of_inertia(helium)
Return total moment of inertia divided by m_He.
Definition: helium_common.F:2122
helium_common::helium_com
pure real(kind=dp) function, dimension(3), public helium_com(helium)
Calculate center of mass.
Definition: helium_common.F:2387
helium_common::helium_rotate
pure subroutine, public helium_rotate(helium, nslices)
Rotate helium particles in imaginary time by nslices.
Definition: helium_common.F:796
helium_common::helium_pbc
subroutine, public helium_pbc(helium, r, enforce)
General PBC routine for helium.
Definition: helium_common.F:81
helium_common::helium_boxmean_3d
subroutine, public helium_boxmean_3d(helium, a, b, c)
Calculate the point equidistant from two other points a and b within the helium box - 3D version.
Definition: helium_common.F:205
helium_common::helium_calc_rho
subroutine, public helium_calc_rho(helium)
Calculate helium density distribution functions and store them in heliumrho_inst.
Definition: helium_common.F:370
helium_common::helium_eval_expansion
real(kind=dp) function, public helium_eval_expansion(helium, r, rp, work, action)
Calculate the pair-product action or energy by evaluating the power series expansion according to Eq....
Definition: helium_common.F:832
helium_common::helium_calc_plength
pure subroutine, public helium_calc_plength(helium)
Calculate probability distribution of the permutation lengths.
Definition: helium_common.F:763
helium_common::helium_update_transition_matrix
subroutine, public helium_update_transition_matrix(helium)
...
Definition: helium_common.F:1058
helium_common::helium_calc_rdf
subroutine, public helium_calc_rdf(helium)
Calculate helium radial distribution functions wrt positions given by <centers> using the atomic weig...
Definition: helium_common.F:630
helium_common::helium_total_projected_area
real(kind=dp) function, dimension(3), public helium_total_projected_area(helium)
Return total projected area.
Definition: helium_common.F:1863
helium_common::helium_spline
real(kind=dp) function, public helium_spline(spl, xx)
...
Definition: helium_common.F:1423
helium_common::helium_total_winding_number
real(kind=dp) function, dimension(3), public helium_total_winding_number(helium)
Return total winding number.
Definition: helium_common.F:1685
helium_common::helium_set_rdf_coord_system
pure subroutine, public helium_set_rdf_coord_system(helium, pint_env)
Set coordinate system for RDF calculations.
Definition: helium_common.F:2363
helium_interactions
Methods that handle helium-solvent and helium-helium interactions.
Definition: helium_interactions.F:13
helium_interactions::helium_solute_e_f
subroutine, public helium_solute_e_f(pint_env, helium, energy)
Calculate total helium-solute interaction energy and forces.
Definition: helium_interactions.F:448
helium_interactions::helium_bead_solute_e_f
subroutine, public helium_bead_solute_e_f(pint_env, helium, helium_part_index, helium_slice_index, helium_r_opt, energy, force)
Calculate general helium-solute interaction energy (and forces) between one helium bead and the corre...
Definition: helium_interactions.F:352
helium_interactions::helium_calc_energy
subroutine, public helium_calc_energy(helium, pint_env)
Calculate the helium energy (including helium-solute interaction)
Definition: helium_interactions.F:74
helium_io
I/O subroutines for helium.
Definition: helium_io.F:13
helium_io::helium_print_vector
subroutine, public helium_print_vector(helium_env, pkey, DATA, uconv, col_label, cmmnt, fname, fpos, avg)
Print a 3D real vector according to printkey <pkey>
Definition: helium_io.F:585
helium_io::helium_print_action
subroutine, public helium_print_action(pint_env, helium_env)
Print helium action file to HELIUMPRINTACTION.
Definition: helium_io.F:988
helium_io::helium_print_rdf
subroutine, public helium_print_rdf(helium_env)
Print radial distribution functions according to HELIUMPRINTRDF.
Definition: helium_io.F:1403
helium_io::helium_print_force
subroutine, public helium_print_force(helium_env)
Print helium force according to HELIUMPRINTFORCE.
Definition: helium_io.F:1785
helium_io::helium_print_accepts
subroutine, public helium_print_accepts(helium_env)
Print acceptance counts according to HELIUMPRINTACCEPTS.
Definition: helium_io.F:726
helium_io::helium_print_energy
subroutine, public helium_print_energy(helium_env)
Print energies according to HELIUMPRINTENERGY.
Definition: helium_io.F:475
helium_io::helium_print_perm
subroutine, public helium_print_perm(helium_env)
Print permutation state according to HELIUMPRINTPERM.
Definition: helium_io.F:795
helium_io::helium_print_rho
subroutine, public helium_print_rho(helium_env)
Print densities according to HELIUMPRINTRHO.
Definition: helium_io.F:1534
helium_io::helium_print_plength
subroutine, public helium_print_plength(helium_env)
Print permutation length according to HELIUMPRINTPLENGTH.
Definition: helium_io.F:1726
helium_io::helium_write_line
subroutine, public helium_write_line(line)
Writes out a line of text to the default output unit.
Definition: helium_io.F:447
helium_io::helium_print_coordinates
subroutine, public helium_print_coordinates(helium_env)
Print coordinates according to HELIUMPRINTCOORDINATES.
Definition: helium_io.F:1107
helium_sampling
Methods for sampling helium variables.
Definition: helium_sampling.F:13
helium_sampling::helium_step
subroutine, public helium_step(helium_env, pint_env)
Perform MC step for helium.
Definition: helium_sampling.F:420
helium_sampling::helium_do_run
subroutine, public helium_do_run(helium_env, globenv)
Performs MC simulation for helium (only)
Definition: helium_sampling.F:81
helium_sampling::helium_sample
subroutine, public helium_sample(helium_env, pint_env)
Sample the helium phase space.
Definition: helium_sampling.F:183
helium_types
Data types representing superfluid helium.
Definition: helium_types.F:15
helium_types::e_id_total
integer, parameter, public e_id_total
Energy contributions - symbolic names for indexing energy arrays.
Definition: helium_types.F:38
helium_worm
Methods dealing with the canonical worm alogrithm.
Definition: helium_worm.F:13
helium_worm::helium_sample_worm
subroutine, public helium_sample_worm(helium, pint_env)
Main worm sampling routine.
Definition: helium_worm.F:50
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::helium_forces_average
integer, parameter, public helium_forces_average
Definition: input_constants.F:330
input_constants::helium_mdist_exponential
integer, parameter, public helium_mdist_exponential
Definition: input_constants.F:351
input_constants::helium_sampling_ceperley
integer, parameter, public helium_sampling_ceperley
Definition: input_constants.F:347
input_constants::helium_forces_last
integer, parameter, public helium_forces_last
Definition: input_constants.F:330
input_constants::helium_mdist_gaussian
integer, parameter, public helium_mdist_gaussian
Definition: input_constants.F:351
input_constants::helium_mdist_quadratic
integer, parameter, public helium_mdist_quadratic
Definition: input_constants.F:351
input_constants::helium_mdist_uniform
integer, parameter, public helium_mdist_uniform
Definition: input_constants.F:351
input_constants::helium_sampling_worm
integer, parameter, public helium_sampling_worm
Definition: input_constants.F:347
input_constants::helium_mdist_singlev
integer, parameter, public helium_mdist_singlev
Definition: input_constants.F:351
input_constants::helium_mdist_linear
integer, parameter, public helium_mdist_linear
Definition: input_constants.F:351
input_cp2k_restarts
Set of routines to dump the restart file of CP2K.
Definition: input_cp2k_restarts.F:13
input_cp2k_restarts::write_restart
subroutine, public write_restart(md_env, force_env, root_section, coords, vels, pint_env, helium_env)
checks if a restart needs to be written and does so, updating all necessary fields in the input file....
Definition: input_cp2k_restarts.F:112
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition: machine.F:123
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144
pint_public
Public path integral routines that can be called from other modules.
Definition: pint_public.F:15
pint_public::pint_com_pos
pure real(kind=dp) function, dimension(3), public pint_com_pos(pint_env)
Return the center of mass of the PI system.
Definition: pint_public.F:45
pint_types
Definition: pint_types.F:8
splines_types
routines for handling splines_types
Definition: splines_types.F:14