(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
13MODULE 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
67CONTAINS
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
1307END 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_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_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_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
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
global_types::global_environment_type
contains the initially parsed file and the initial parallel environment
Definition global_types.F:66
helium_types::helium_solvent_p_type
data structure for array of solvent helium environments
Definition helium_types.F:300
helium_types::helium_solvent_type
data structure for solvent helium
Definition helium_types.F:93
pint_types::pint_env_type
environment for a path integral run
Definition pint_types.F:112
splines_types::spline_data_type
Data-structure that holds all needed information about a specific spline interpolation.
Definition splines_types.F:40