Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 10 : 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 10 : NULLIFY (logger)
90 20 : logger => cp_get_default_logger()
91 :
92 10 : num_steps = 0
93 :
94 10 : IF (ASSOCIATED(helium_env)) THEN
95 20 : DO k = 1, SIZE(helium_env)
96 :
97 10 : NULLIFY (helium_env(k)%helium%logger)
98 10 : 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 10 : IF (k == 1) THEN
104 10 : 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 10 : iter_nr=helium_env(k)%helium%first_step)
107 : END IF
108 10 : tot_steps = helium_env(k)%helium%first_step
109 10 : num_steps = helium_env(k)%helium%num_steps
110 :
111 : ! set the properties accumulated over the whole MC process to 0
112 40 : helium_env(k)%helium%proarea%accu(:) = 0.0_dp
113 40 : helium_env(k)%helium%prarea2%accu(:) = 0.0_dp
114 40 : helium_env(k)%helium%wnmber2%accu(:) = 0.0_dp
115 40 : helium_env(k)%helium%mominer%accu(:) = 0.0_dp
116 10 : IF (helium_env(k)%helium%rho_present) THEN
117 4222 : helium_env(k)%helium%rho_accu(:, :, :, :) = 0.0_dp
118 : END IF
119 20 : IF (helium_env(k)%helium%rdf_present) THEN
120 1002 : 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 10 : CALL logger%para_env%bcast(num_steps)
127 10 : CALL logger%para_env%bcast(tot_steps)
128 :
129 48 : DO step = 1, num_steps
130 :
131 38 : tot_steps = tot_steps + 1
132 :
133 38 : IF (ASSOCIATED(helium_env)) THEN
134 76 : DO k = 1, SIZE(helium_env)
135 38 : IF (k == 1) THEN
136 : CALL cp_iterate(helium_env(k)%helium%logger%iter_info, &
137 : last=(step == helium_env(k)%helium%num_steps), &
138 38 : iter_nr=tot_steps)
139 : END IF
140 76 : helium_env(k)%helium%current_step = tot_steps
141 : END DO
142 : END IF
143 :
144 38 : CALL helium_step(helium_env, pint_env)
145 :
146 : ! call write_restart here to avoid interference with PINT write_restart
147 38 : IF (ASSOCIATED(helium_env)) THEN
148 38 : 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 38 : CALL external_control(should_stop, "PINT", globenv=globenv)
153 48 : IF (should_stop) EXIT
154 :
155 : END DO
156 :
157 : ! remove iteration level
158 10 : IF (ASSOCIATED(helium_env)) THEN
159 10 : CALL cp_rm_iter_level(helium_env(1)%helium%logger%iter_info, "PINT")
160 : END IF
161 :
162 10 : RETURN
163 320 : 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 117 : 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 117 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work_2d
192 : TYPE(cp_logger_type), POINTER :: logger
193 :
194 117 : NULLIFY (logger)
195 117 : IF (SIZE(helium_env) < 1) THEN
196 0 : logger => cp_get_default_logger()
197 : ELSE
198 117 : CALL cp_logger_create(logger, para_env=helium_env(1)%comm, template_logger=cp_get_default_logger())
199 117 : CALL cp_add_default_logger(logger)
200 : END IF
201 :
202 270 : DO k = 1, SIZE(helium_env)
203 :
204 153 : IF (helium_env(k)%helium%solute_present) THEN
205 7066 : helium_env(k)%helium%force_avrg(:, :) = 0.0_dp
206 109 : helium_env(k)%helium%current_step = pint_env%iter
207 436 : helium_env(k)%helium%origin = pint_com_pos(pint_env)
208 : END IF
209 :
210 1683 : helium_env(k)%helium%energy_avrg(:) = 0.0_dp
211 3369 : helium_env(k)%helium%plength_avrg(:) = 0.0_dp
212 2178 : 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 153 : inv_xn = 0.0_dp
218 180 : SELECT CASE (helium_env(k)%helium%sampling_method)
219 :
220 : CASE (helium_sampling_worm)
221 :
222 27 : CALL helium_sample_worm(helium_env(k)%helium, pint_env)
223 :
224 27 : 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 7146 : 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 7020 : rnd = helium_env(k)%helium%rng_stream_uniform%next()
235 7020 : nslices = INT(rnd*helium_env(k)%helium%beads)
236 7020 : CALL helium_rotate(helium_env(k)%helium, nslices)
237 :
238 7020 : CALL helium_try_permutations(helium_env(k)%helium, pint_env)
239 :
240 : ! calculate & accumulate instantaneous properties for averaging
241 7020 : IF (helium_env(k)%helium%solute_present) THEN
242 3620 : IF (helium_env(k)%helium%get_helium_forces == helium_forces_average) THEN
243 2620 : 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 120520 : helium_env(k)%helium%force_inst(:, :)
246 : END IF
247 : END IF
248 7020 : CALL helium_calc_energy(helium_env(k)%helium, pint_env)
249 :
250 77220 : helium_env(k)%helium%energy_avrg(:) = helium_env(k)%helium%energy_avrg(:) + helium_env(k)%helium%energy_inst(:)
251 7020 : CALL helium_calc_plength(helium_env(k)%helium)
252 220446 : 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 126 : IF ((helium_env(k)%helium%solute_present) .AND. (helium_env(k)%helium%get_helium_forces == helium_forces_last)) THEN
263 10 : CALL helium_solute_e_f(pint_env, helium_env(k)%helium, rtmp)
264 460 : 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 126 : CALL helium_rotate(helium_env(k)%helium, -helium_env(k)%helium%relrot)
271 :
272 126 : inv_xn = 1.0_dp/REAL(helium_env(k)%helium%iter_rot, dp)
273 :
274 : CASE DEFAULT
275 0 : WRITE (msg_str, *) helium_env(k)%helium%sampling_method
276 0 : msg_str = "Sampling method ("//TRIM(ADJUSTL(msg_str))//") not supported."
277 153 : CPABORT(msg_str)
278 :
279 : END SELECT
280 :
281 : ! actually average the properties over the outer MC loop
282 153 : IF ((helium_env(k)%helium%solute_present) .AND. (helium_env(k)%helium%get_helium_forces == helium_forces_average)) THEN
283 3772 : helium_env(k)%helium%force_avrg(:, :) = helium_env(k)%helium%force_avrg(:, :)*inv_xn
284 : END IF
285 1683 : helium_env(k)%helium%energy_avrg(:) = helium_env(k)%helium%energy_avrg(:)*inv_xn
286 3369 : helium_env(k)%helium%plength_avrg(:) = helium_env(k)%helium%plength_avrg(:)*inv_xn
287 :
288 : ! instantaneous properties
289 612 : helium_env(k)%helium%proarea%inst(:) = helium_total_projected_area(helium_env(k)%helium)
290 612 : helium_env(k)%helium%prarea2%inst(:) = helium_env(k)%helium%proarea%inst(:)*helium_env(k)%helium%proarea%inst(:)
291 612 : helium_env(k)%helium%wnumber%inst(:) = helium_total_winding_number(helium_env(k)%helium)
292 612 : helium_env(k)%helium%wnmber2%inst(:) = helium_env(k)%helium%wnumber%inst(:)*helium_env(k)%helium%wnumber%inst(:)
293 612 : 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 612 : helium_env(k)%helium%proarea%accu(:) = helium_env(k)%helium%proarea%accu(:) + helium_env(k)%helium%proarea%inst(:)
297 612 : helium_env(k)%helium%prarea2%accu(:) = helium_env(k)%helium%prarea2%accu(:) + helium_env(k)%helium%prarea2%inst(:)
298 612 : helium_env(k)%helium%wnmber2%accu(:) = helium_env(k)%helium%wnmber2%accu(:) + helium_env(k)%helium%wnmber2%inst(:)
299 612 : helium_env(k)%helium%mominer%accu(:) = helium_env(k)%helium%mominer%accu(:) + helium_env(k)%helium%mominer%inst(:)
300 153 : IF (helium_env(k)%helium%rho_present) THEN
301 12 : CALL helium_calc_rho(helium_env(k)%helium)
302 : helium_env(k)%helium%rho_accu(:, :, :, :) = helium_env(k)%helium%rho_accu(:, :, :, :) + &
303 25332 : helium_env(k)%helium%rho_inst(:, :, :, :)
304 : END IF
305 153 : IF (helium_env(k)%helium%rdf_present) THEN
306 12 : CALL helium_set_rdf_coord_system(helium_env(k)%helium, pint_env)
307 12 : CALL helium_calc_rdf(helium_env(k)%helium)
308 6012 : 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 153 : nsteps = helium_env(k)%helium%current_step - helium_env(k)%helium%first_step
313 153 : iweight = helium_env(k)%helium%averages_iweight
314 153 : rweight = REAL(iweight, dp)
315 153 : rtmp = 1.0_dp/(REAL(MAX(1, nsteps + iweight), dp))
316 : helium_env(k)%helium%proarea%ravr(:) = (helium_env(k)%helium%proarea%accu(:) + &
317 612 : rweight*helium_env(k)%helium%proarea%rstr(:))*rtmp
318 : helium_env(k)%helium%prarea2%ravr(:) = (helium_env(k)%helium%prarea2%accu(:) + &
319 612 : rweight*helium_env(k)%helium%prarea2%rstr(:))*rtmp
320 : helium_env(k)%helium%wnmber2%ravr(:) = (helium_env(k)%helium%wnmber2%accu(:) + &
321 612 : rweight*helium_env(k)%helium%wnmber2%rstr(:))*rtmp
322 : helium_env(k)%helium%mominer%ravr(:) = (helium_env(k)%helium%mominer%accu(:) + &
323 729 : 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 117 : num_env = helium_env(1)%helium%num_env
330 117 : inv_num_env = 1.0_dp/REAL(num_env, dp)
331 :
332 : !energy_avrg:
333 153 : DO k = 2, SIZE(helium_env)
334 : helium_env(1)%helium%energy_avrg(:) = helium_env(1)%helium%energy_avrg(:) + &
335 513 : helium_env(k)%helium%energy_avrg(:)
336 : END DO
337 117 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%energy_avrg(:))
338 1287 : helium_env(1)%helium%energy_avrg(:) = helium_env(1)%helium%energy_avrg(:)*inv_num_env
339 153 : DO k = 2, SIZE(helium_env)
340 513 : helium_env(k)%helium%energy_avrg(:) = helium_env(1)%helium%energy_avrg(:)
341 : END DO
342 :
343 : !plength_avrg:
344 153 : DO k = 2, SIZE(helium_env)
345 : helium_env(1)%helium%plength_avrg(:) = helium_env(1)%helium%plength_avrg(:) + &
346 333 : helium_env(k)%helium%plength_avrg(:)
347 : END DO
348 6189 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%plength_avrg(:))
349 3153 : helium_env(1)%helium%plength_avrg(:) = helium_env(1)%helium%plength_avrg(:)*inv_num_env
350 153 : DO k = 2, SIZE(helium_env)
351 333 : helium_env(k)%helium%plength_avrg(:) = helium_env(1)%helium%plength_avrg(:)
352 : END DO
353 :
354 : !num_accepted:
355 153 : DO k = 2, SIZE(helium_env)
356 : helium_env(1)%helium%num_accepted(:, :) = helium_env(1)%helium%num_accepted(:, :) + &
357 369 : helium_env(k)%helium%num_accepted(:, :)
358 : END DO
359 3735 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%num_accepted(:, :))
360 1926 : helium_env(1)%helium%num_accepted(:, :) = helium_env(1)%helium%num_accepted(:, :)*inv_num_env
361 153 : DO k = 2, SIZE(helium_env)
362 369 : helium_env(k)%helium%num_accepted(:, :) = helium_env(1)%helium%num_accepted(:, :)
363 : END DO
364 :
365 : !force_avrg:
366 117 : IF (helium_env(1)%helium%solute_present) THEN
367 73 : IF (helium_env(1)%helium%get_helium_forces == helium_forces_average) THEN
368 82 : DO k = 2, SIZE(helium_env)
369 : helium_env(1)%helium%force_avrg(:, :) = helium_env(1)%helium%force_avrg(:, :) + &
370 1702 : helium_env(k)%helium%force_avrg(:, :)
371 : END DO
372 4186 : CALL helium_env(1)%comm%sum(helium_env(1)%helium%force_avrg(:, :))
373 2116 : helium_env(1)%helium%force_avrg(:, :) = helium_env(1)%helium%force_avrg(:, :)*inv_num_env
374 82 : DO k = 2, SIZE(helium_env)
375 1702 : helium_env(k)%helium%force_avrg(:, :) = helium_env(1)%helium%force_avrg(:, :)
376 : END DO
377 27 : ELSE IF (helium_env(1)%helium%get_helium_forces == helium_forces_last) THEN
378 27 : IF (logger%para_env%is_source()) THEN
379 : sel_mp_source = INT(helium_env(1)%helium%rng_stream_uniform%next() &
380 16 : *REAL(helium_env(1)%helium%num_env, dp))
381 : END IF
382 27 : CALL helium_env(1)%comm%bcast(sel_mp_source, logger%para_env%source)
383 :
384 27 : offset = 0
385 38 : DO i = 1, logger%para_env%mepos
386 38 : offset = offset + helium_env(1)%env_all(i)
387 : END DO
388 : ALLOCATE (work_2d(SIZE(helium_env(1)%helium%force_avrg, 1), &
389 108 : SIZE(helium_env(1)%helium%force_avrg, 2)))
390 3294 : work_2d(:, :) = 0.0_dp
391 27 : IF (sel_mp_source >= offset .AND. sel_mp_source < offset + SIZE(helium_env)) THEN
392 2032 : work_2d(:, :) = helium_env(sel_mp_source - offset + 1)%helium%force_avrg(:, :)
393 : END IF
394 27 : CALL helium_env(1)%comm%sum(work_2d(:, :))
395 54 : DO k = 1, SIZE(helium_env)
396 3321 : helium_env(k)%helium%force_avrg(:, :) = work_2d(:, :)
397 : END DO
398 27 : DEALLOCATE (work_2d)
399 : END IF
400 : END IF
401 :
402 117 : IF (SIZE(helium_env) > 0) THEN
403 117 : CALL cp_rm_default_logger()
404 117 : CALL cp_logger_release(logger)
405 : END IF
406 :
407 117 : 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 110 : 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 110 : REAL(KIND=dp), DIMENSION(:), POINTER :: DATA
430 :
431 110 : CALL timeset(routineN, handle)
432 110 : time_start = m_walltime()
433 :
434 110 : IF (ASSOCIATED(helium_env)) THEN
435 : ! perform the actual phase space sampling
436 105 : CALL helium_sample(helium_env, pint_env)
437 :
438 : ! print properties
439 105 : CALL helium_print_energy(helium_env)
440 105 : CALL helium_print_plength(helium_env)
441 105 : CALL helium_print_accepts(helium_env)
442 105 : CALL helium_print_force(helium_env)
443 105 : CALL helium_print_perm(helium_env)
444 105 : CALL helium_print_coordinates(helium_env)
445 105 : CALL helium_print_action(pint_env, helium_env)
446 :
447 105 : IF (helium_env(1)%helium%rdf_present) CALL helium_print_rdf(helium_env)
448 105 : IF (helium_env(1)%helium%rho_present) CALL helium_print_rho(helium_env)
449 :
450 105 : NULLIFY (DATA)
451 315 : ALLOCATE (DATA(3*SIZE(helium_env)))
452 :
453 105 : WRITE (stmp, *) helium_env(1)%helium%apref
454 510 : DATA(:) = 0.0_dp
455 240 : DO k = 1, SIZE(helium_env)
456 645 : 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 420 : "helium-pa")
465 :
466 510 : DATA(:) = 0.0_dp
467 240 : DO k = 1, SIZE(helium_env)
468 645 : 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 420 : .TRUE.)
479 :
480 510 : DATA(:) = 0.0_dp
481 240 : DO k = 1, SIZE(helium_env)
482 645 : 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 420 : "helium-mi")
491 :
492 510 : DATA(:) = 0.0_dp
493 240 : DO k = 1, SIZE(helium_env)
494 645 : 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 420 : .TRUE.)
505 :
506 510 : DATA(:) = 0.0_dp
507 240 : DO k = 1, SIZE(helium_env)
508 645 : DATA((k - 1)*3 + 1:k*3) = helium_env(k)%helium%wnumber%inst
509 : END DO
510 105 : 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 420 : "helium-wn")
518 :
519 510 : DATA(:) = 0.0_dp
520 240 : DO k = 1, SIZE(helium_env)
521 645 : 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 420 : .TRUE.)
532 105 : DEALLOCATE (DATA)
533 :
534 105 : time_stop = m_walltime()
535 105 : time_used = time_stop - time_start
536 105 : time_unit = "sec"
537 105 : IF (time_used >= 60.0_dp) THEN
538 0 : time_used = time_used/60.0_dp
539 0 : time_unit = "min"
540 0 : IF (time_used >= 60.0_dp) THEN
541 0 : time_used = time_used/60.0_dp
542 0 : time_unit = "hours"
543 : END IF
544 : END IF
545 105 : msgstr = "MC step"
546 105 : stmp = ""
547 105 : WRITE (stmp, *) helium_env(1)%helium%current_step
548 105 : msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(stmp))//" of"
549 105 : stmp = ""
550 105 : WRITE (stmp, *) helium_env(1)%helium%last_step
551 105 : msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(stmp))//" in"
552 105 : stmp = ""
553 105 : WRITE (stmp, '(F20.1)') time_used
554 105 : msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(stmp))
555 105 : msgstr = TRIM(ADJUSTL(msgstr))//" "//TRIM(ADJUSTL(time_unit))//"."
556 105 : CALL helium_write_line(TRIM(msgstr))
557 :
558 : ! print out the total energy - for regtest evaluation
559 105 : stmp = ""
560 105 : WRITE (stmp, *) helium_env(1)%helium%energy_avrg(e_id_total)
561 105 : msgstr = "Total energy = "//TRIM(ADJUSTL(stmp))
562 105 : CALL helium_write_line(TRIM(msgstr))
563 : END IF
564 :
565 110 : CALL timestop(handle)
566 :
567 110 : RETURN
568 110 : 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 7020 : 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 7020 : 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==helium%work) - see a comment there,
592 : ! besides there is a magic regarding helium%permutation
593 : ! you have been warned
594 17717940 : helium%work(:, :, :) = helium%pos(:, :, :)
595 :
596 : ! the inner MC loop (without rotation in imaginary time)
597 1348020 : DO ni = 1, helium%iter_norot
598 :
599 : ! set the probability threshold for m_value: 1/(1+(m-1)/helium%m_ratio)
600 1341000 : 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 1341000 : SELECT CASE (helium%m_dist_type)
605 :
606 : CASE (helium_mdist_singlev)
607 0 : x = helium%rng_stream_uniform%next()
608 0 : IF (x < r) THEN
609 0 : cyclen = 1
610 : ELSE
611 0 : cyclen = helium%m_value
612 : END IF
613 :
614 : CASE (helium_mdist_uniform)
615 1341000 : x = helium%rng_stream_uniform%next()
616 1341000 : IF (x < r) THEN
617 350513 : cyclen = helium%m_value
618 : ELSE
619 : DO
620 1320591 : x = helium%rng_stream_uniform%next()
621 1320591 : cyclen = INT(helium%maxcycle*x) + 1
622 1320591 : IF (cyclen /= helium%m_value) EXIT
623 : END DO
624 : END IF
625 :
626 : CASE (helium_mdist_linear)
627 0 : x = helium%rng_stream_uniform%next()
628 0 : IF (x < r) THEN
629 0 : cyclen = helium%m_value
630 : ELSE
631 : DO
632 0 : x = helium%rng_stream_uniform%next()
633 0 : y = SQRT(2.0_dp*x)
634 0 : cyclen = INT(helium%maxcycle*y/SQRT(2.0_dp)) + 1
635 0 : IF (cyclen /= helium%m_value) EXIT
636 : END DO
637 : END IF
638 :
639 : CASE (helium_mdist_quadratic)
640 0 : x = helium%rng_stream_uniform%next()
641 0 : IF (x < r) THEN
642 0 : cyclen = helium%m_value
643 : ELSE
644 : DO
645 0 : x = helium%rng_stream_uniform%next()
646 0 : y = (3.0_dp*x)**(1.0_dp/3.0_dp)
647 0 : z = 3.0_dp**(1.0_dp/3.0_dp)
648 0 : cyclen = INT(helium%maxcycle*y/z) + 1
649 0 : IF (cyclen /= helium%m_value) EXIT
650 : END DO
651 : END IF
652 :
653 : CASE (helium_mdist_exponential)
654 0 : x = helium%rng_stream_uniform%next()
655 0 : IF (x < r) THEN
656 0 : cyclen = helium%m_value
657 : ELSE
658 : DO
659 : DO
660 0 : x = helium%rng_stream_uniform%next()
661 0 : IF (x >= 0.01_dp) EXIT
662 : END DO
663 0 : z = -LOG(0.01_dp)
664 0 : y = LOG(x)/z + 1.0_dp;
665 0 : cyclen = INT(helium%maxcycle*y) + 1
666 0 : IF (cyclen /= helium%m_value) EXIT
667 : END DO
668 : END IF
669 :
670 : CASE (helium_mdist_gaussian)
671 0 : x = helium%rng_stream_uniform%next()
672 0 : IF (x < r) THEN
673 0 : cyclen = 1
674 : ELSE
675 : DO
676 0 : x = helium%rng_stream_gaussian%next()
677 0 : cyclen = INT(x*0.75_dp + helium%m_value - 0.5_dp) + 1
678 0 : IF (cyclen /= 1) EXIT
679 : END DO
680 : END IF
681 :
682 : CASE DEFAULT
683 0 : WRITE (stmp, *) helium%m_dist_type
684 0 : err_str = "M distribution type unknown ("//TRIM(ADJUSTL(stmp))//")"
685 0 : CPABORT(err_str)
686 1341000 : cyclen = 1
687 :
688 : END SELECT
689 :
690 1341000 : IF (cyclen < 1) cyclen = 1
691 1341000 : IF (cyclen > helium%maxcycle) cyclen = helium%maxcycle
692 1341000 : helium%num_accepted(1, cyclen) = helium%num_accepted(1, cyclen) + 1
693 :
694 : ! check, if permutation of this length can be constructed
695 1341000 : IF (cyclen == 1) THEN
696 350513 : rnd = helium%rng_stream_uniform%next()
697 350513 : helium%ptable(1) = 1 + INT(rnd*helium%atoms)
698 350513 : helium%ptable(2) = -1
699 350513 : helium%pweight = 0.0_dp
700 : selected = .TRUE.
701 : ELSE
702 990487 : selected = helium_select_permutation(helium, cyclen)
703 : END IF
704 :
705 990487 : IF (selected) THEN
706 : ! permutation was successfully selected - actually sample this permutation
707 350548 : 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 357568 : IF (accepted) THEN
719 : ! positions changed
720 128676 : 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 60356 : 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 0 : helium%center(:) = helium_com(helium)
728 : END IF
729 : END IF
730 : END IF
731 :
732 : END DO
733 :
734 7020 : 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 350548 : 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 350548 : 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 350548 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work3
759 : REAL(KIND=dp), DIMENSION(3) :: bis, biso, new_com, rm1, rm2, tmp1, tmp2
760 350548 : 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 350548 : p => helium%ptable
771 350548 : prev_ds = helium%pweight
772 :
773 350548 : helium%num_accepted(2, cyclen) = helium%num_accepted(2, cyclen) + 1
774 350548 : level = 1
775 350548 : res = .FALSE.
776 :
777 : ! nothing to be done
778 350548 : 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 350548 : pos => helium%pos
788 350548 : work => helium%work
789 350548 : perm => helium%permutation
790 350548 : u0 => helium%u0
791 350548 : cell_size = (0.5_dp*helium%cell_size)**2
792 350548 : nperiodic = .NOT. helium%periodic
793 :
794 350548 : pds = prev_ds
795 350548 : ifix = helium%beads - helium%bisection + 1
796 :
797 : ! sanity checks
798 : !
799 350548 : IF (ifix < 1) THEN
800 0 : WRITE (stmp, *) ifix
801 0 : err_str = "ifix<1 test failed (ifix="//TRIM(ADJUSTL(stmp))//")"
802 0 : CPABORT(err_str)
803 : END IF
804 : !
805 350548 : j = 1
806 350548 : k = helium%bisection
807 1051644 : DO
808 1402192 : IF (k < 2) EXIT
809 1051644 : j = j*2
810 1051644 : k = k/2
811 : END DO
812 : !
813 350548 : IF (j /= helium%bisection) THEN
814 0 : WRITE (stmp, *) helium%bisection
815 0 : err_str = "helium%bisection not a power of 2 (helium%bisection="//TRIM(ADJUSTL(stmp))//")"
816 0 : CPABORT(err_str)
817 : END IF
818 : !
819 350548 : IF (helium%bisection < 2) THEN
820 0 : WRITE (stmp, *) helium%bisection
821 0 : err_str = "helium%bisection less than 2 (helium%bisection="//TRIM(ADJUSTL(stmp))//")"
822 0 : CPABORT(err_str)
823 : END IF
824 :
825 350548 : stride = helium%bisection
826 399882 : DO
827 750430 : IF (stride <= 2) EXIT
828 : ! calc new trial positions
829 : ! trial action: 0.5*stride*endpointapprox
830 585672 : sigma = SQRT(0.5_dp*helium%hb2m*(stride/2)*helium%tau)
831 585672 : dtk = 0.0_dp
832 585672 : ds = 0.0_dp
833 :
834 585672 : j = ifix + stride/2
835 235124 : DO
836 820796 : IF (j > helium%beads - stride/2) EXIT
837 235124 : pk1 = j - stride/2
838 235124 : pk2 = j + stride/2
839 : ! calculate log(T(s)):
840 470250 : DO k = 1, cyclen
841 235126 : CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, p(k), pk2), bis)
842 940504 : tmp1(:) = bis(:) - pos(:, p(k), j)
843 235126 : CALL helium_pbc(helium, tmp1)
844 940504 : tmp1(:) = tmp1(:)/sigma
845 470250 : 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 470250 : DO k = 1, cyclen
849 235126 : CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, p(k), pk2), tmp1)
850 940504 : DO c = 1, 3
851 705378 : x = helium%rng_stream_gaussian%next(variance=1.0_dp)
852 705378 : x = sigma*x
853 705378 : tmp1(c) = tmp1(c) + x
854 940504 : tmp2(c) = x
855 : END DO
856 235126 : CALL helium_pbc(helium, tmp1)
857 235126 : CALL helium_pbc(helium, tmp2)
858 940504 : work(:, p(k), j) = tmp1(:)
859 940504 : tmp2(:) = tmp2(:)/sigma
860 470250 : dtk = dtk + 0.5_dp*(tmp2(1)*tmp2(1) + tmp2(2)*tmp2(2) + tmp2(3)*tmp2(3))
861 : END DO
862 235124 : j = j + stride
863 : END DO
864 :
865 585672 : j = helium%beads - stride/2 + 1
866 585672 : pk1 = j - stride/2
867 1171381 : DO k = 1, cyclen
868 585709 : CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, perm(p(k)), 1), bis)
869 2342836 : tmp1(:) = bis(:) - pos(:, p(k), j)
870 585709 : CALL helium_pbc(helium, tmp1)
871 2342836 : tmp1(:) = tmp1(:)/sigma
872 1171381 : dtk = dtk - 0.5_dp*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
873 : END DO
874 1171381 : DO k = 1, cyclen
875 585709 : CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, perm(p(1 + MOD(k, cyclen))), 1), tmp1)
876 2342836 : DO c = 1, 3
877 1757127 : x = helium%rng_stream_gaussian%next(variance=1.0_dp)
878 1757127 : x = sigma*x
879 1757127 : tmp1(c) = tmp1(c) + x
880 2342836 : tmp2(c) = x
881 : END DO
882 585709 : CALL helium_pbc(helium, tmp1)
883 585709 : CALL helium_pbc(helium, tmp2)
884 2342836 : work(:, p(k), j) = tmp1(:)
885 2342836 : tmp2(:) = tmp2(:)/sigma
886 1171381 : 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 585672 : x = 1.0_dp/(helium%tau*helium%hb2m*stride)
891 585672 : j = ifix
892 1055920 : DO
893 1641592 : IF (j > helium%beads - stride/2) EXIT
894 1055920 : pk1 = j + stride/2
895 2111881 : DO k = 1, cyclen
896 4223844 : tmp1(:) = pos(:, p(k), j) - pos(:, p(k), pk1)
897 1055961 : CALL helium_pbc(helium, tmp1)
898 1055961 : ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
899 4223844 : tmp1(:) = work(:, p(k), j) - work(:, p(k), pk1)
900 1055961 : CALL helium_pbc(helium, tmp1)
901 1055961 : ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
902 : ! interaction change
903 1055961 : IF (helium%solute_present) THEN
904 543746 : CALL helium_bead_solute_e_f(pint_env, helium, p(k), pk1, energy=e1)
905 543746 : CALL helium_bead_solute_e_f(pint_env, helium, p(k), pk1, work(:, p(k), pk1), e2)
906 543746 : ds = ds + (stride/2)*(e1 - e2)*helium%tau
907 : END IF
908 32647563 : DO l = 1, helium%atoms
909 32647563 : IF (l /= p(k)) THEN
910 122142564 : tmp1(:) = pos(:, p(k), pk1) - pos(:, l, pk1)
911 30535641 : CALL helium_pbc(helium, tmp1)
912 30535641 : r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
913 30535641 : IF ((r < cell_size) .OR. nperiodic) THEN
914 28604255 : r = SQRT(r)
915 28604255 : ds = ds + REAL(stride/2, dp)*helium_spline(u0, r)
916 : END IF
917 122142564 : tmp1(:) = work(:, p(k), pk1) - work(:, l, pk1)
918 30535641 : CALL helium_pbc(helium, tmp1)
919 30535641 : r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
920 30535641 : IF ((r < cell_size) .OR. nperiodic) THEN
921 28603810 : r = SQRT(r)
922 28603810 : 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 2111881 : IF (k < cyclen) THEN
928 82 : DO l = k + 1, cyclen
929 164 : tmp1(:) = pos(:, p(k), pk1) - pos(:, p(l), pk1)
930 41 : CALL helium_pbc(helium, tmp1)
931 41 : r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
932 41 : IF ((r < cell_size) .OR. nperiodic) THEN
933 41 : r = SQRT(r)
934 41 : ds = ds - REAL(stride/2, dp)*helium_spline(u0, r)
935 : END IF
936 164 : tmp1(:) = work(:, p(k), pk1) - work(:, p(l), pk1)
937 41 : CALL helium_pbc(helium, tmp1)
938 41 : r = tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3)
939 82 : IF ((r < cell_size) .OR. nperiodic) THEN
940 41 : r = SQRT(r)
941 41 : ds = ds + REAL(stride/2, dp)*helium_spline(u0, r)
942 : END IF
943 : END DO
944 : END IF
945 : END DO
946 1055920 : j = j + stride/2
947 : END DO
948 : ! last link
949 585672 : pk1 = helium%beads - stride/2 + 1
950 1171381 : DO k = 1, cyclen
951 2342836 : tmp1(:) = pos(:, p(k), pk1) - pos(:, perm(p(k)), 1)
952 585709 : CALL helium_pbc(helium, tmp1)
953 585709 : ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
954 2342836 : tmp1(:) = work(:, p(k), pk1) - work(:, perm(p(1 + MOD(k, cyclen))), 1)
955 585709 : CALL helium_pbc(helium, tmp1)
956 1171381 : 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 585672 : rtmp = helium%rng_stream_uniform%next()
960 : ! IF ((dtk+ds-pds < 0.0_dp).AND.(EXP(dtk+ds-pds)<rtmp)) THEN
961 585672 : IF (dtk + ds - pds < 0.0_dp) THEN
962 376906 : IF (EXP(dtk + ds - pds) < rtmp) THEN
963 1672110 : DO c = ifix, helium%beads
964 3158710 : DO k = 1, cyclen
965 11892520 : 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 399882 : helium%num_accepted(level + 2, cyclen) = helium%num_accepted(level + 2, cyclen) + 1
973 399882 : level = level + 1
974 399882 : pds = ds
975 399882 : 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 164758 : sigma = SQRT(0.5_dp*helium%hb2m*helium%tau)
981 164758 : dtk = 0.0_dp
982 164758 : ds = 0.0_dp
983 659032 : DO j = ifix + 1, helium%beads - 1, 2
984 494274 : pk1 = j - 1
985 494274 : pk2 = j + 1
986 : ! calculate log(T(s)):
987 988548 : DO k = 1, cyclen
988 494274 : CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, p(k), pk2), bis)
989 1977096 : tmp1(:) = bis(:) - pos(:, p(k), j)
990 494274 : CALL helium_pbc(helium, tmp1)
991 1977096 : tmp1(:) = tmp1(:)/sigma
992 988548 : 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 1153306 : DO k = 1, cyclen
996 494274 : CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, p(k), pk2), tmp1)
997 1977096 : DO c = 1, 3
998 1482822 : x = helium%rng_stream_gaussian%next(variance=1.0_dp)
999 1482822 : x = sigma*x
1000 1482822 : tmp1(c) = tmp1(c) + x
1001 1977096 : tmp2(c) = x
1002 : END DO
1003 494274 : CALL helium_pbc(helium, tmp1)
1004 494274 : CALL helium_pbc(helium, tmp2)
1005 1977096 : work(:, p(k), j) = tmp1(:)
1006 1977096 : tmp2(:) = tmp2(:)/sigma
1007 988548 : dtk = dtk + 0.5_dp*(tmp2(1)*tmp2(1) + tmp2(2)*tmp2(2) + tmp2(3)*tmp2(3))
1008 : END DO
1009 : END DO
1010 164758 : j = helium%beads
1011 164758 : pk1 = j - 1
1012 329516 : DO k = 1, cyclen
1013 164758 : CALL helium_boxmean_3d(helium, pos(:, p(k), pk1), pos(:, perm(p(k)), 1), bis)
1014 659032 : tmp1(:) = bis(:) - pos(:, p(k), j)
1015 164758 : CALL helium_pbc(helium, tmp1)
1016 659032 : tmp1(:) = tmp1(:)/sigma
1017 329516 : dtk = dtk - 0.5_dp*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
1018 : END DO
1019 329516 : DO k = 1, cyclen
1020 164758 : CALL helium_boxmean_3d(helium, work(:, p(k), pk1), work(:, perm(p(1 + MOD(k, cyclen))), 1), tmp1)
1021 659032 : DO c = 1, 3
1022 494274 : x = helium%rng_stream_gaussian%next(variance=1.0_dp)
1023 494274 : x = sigma*x
1024 494274 : tmp1(c) = tmp1(c) + x
1025 659032 : tmp2(c) = x
1026 : END DO
1027 164758 : CALL helium_pbc(helium, tmp1)
1028 164758 : CALL helium_pbc(helium, tmp2)
1029 659032 : work(:, p(k), j) = tmp1(:)
1030 659032 : tmp2 = tmp2/sigma
1031 329516 : 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 164758 : IF (helium%solute_present) THEN
1038 772362 : DO j = ifix, helium%beads
1039 1458906 : DO k = 1, cyclen
1040 686544 : CALL helium_bead_solute_e_f(pint_env, helium, p(k), j, energy=e1)
1041 686544 : CALL helium_bead_solute_e_f(pint_env, helium, p(k), j, work(:, p(k), j), e2)
1042 2059632 : ds = ds + (e1 - e2)*helium%tau
1043 : END DO
1044 : END DO
1045 : END IF
1046 494274 : ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
1047 164758 : x = 1.0_dp/(helium%tau*helium%hb2m*stride)
1048 1318064 : DO j = ifix, helium%beads - 1
1049 1153306 : pk1 = j + 1
1050 2471370 : DO k = 1, cyclen
1051 4613224 : tmp1(:) = pos(:, p(k), j) - pos(:, p(k), pk1)
1052 1153306 : CALL helium_pbc(helium, tmp1)
1053 1153306 : ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
1054 4613224 : tmp1(:) = work(:, p(k), j) - work(:, p(k), pk1)
1055 1153306 : CALL helium_pbc(helium, tmp1)
1056 1153306 : ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
1057 34347327 : DO l = 1, helium%atoms
1058 34347327 : IF (l /= p(k)) THEN
1059 128162860 : rm1(:) = pos(:, p(k), j) - pos(:, l, j)
1060 128162860 : rm2(:) = pos(:, p(k), pk1) - pos(:, l, pk1)
1061 32040715 : ds = ds + helium_eval_expansion(helium, rm1, rm2, work3)
1062 128162860 : rm1(:) = work(:, p(k), j) - work(:, l, j)
1063 128162860 : rm2(:) = work(:, p(k), pk1) - work(:, l, pk1)
1064 32040715 : 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 2306612 : IF (k < cyclen) THEN
1069 0 : DO l = k + 1, cyclen
1070 0 : rm1(:) = pos(:, p(k), j) - pos(:, p(l), j)
1071 0 : rm2(:) = pos(:, p(k), pk1) - pos(:, p(l), pk1)
1072 0 : ds = ds - helium_eval_expansion(helium, rm1, rm2, work3)
1073 0 : rm1(:) = work(:, p(k), j) - work(:, p(l), j)
1074 0 : rm2(:) = work(:, p(k), pk1) - work(:, p(l), pk1)
1075 0 : ds = ds + helium_eval_expansion(helium, rm1, rm2, work3)
1076 : END DO
1077 : END IF
1078 : END DO
1079 : END DO
1080 164758 : j = helium%beads
1081 164758 : pk1 = 1
1082 329516 : DO k = 1, cyclen
1083 659032 : tmp1(:) = pos(:, p(k), j) - pos(:, perm(p(k)), 1)
1084 164758 : CALL helium_pbc(helium, tmp1)
1085 164758 : ds = ds + x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
1086 659032 : tmp1(:) = work(:, p(k), j) - work(:, perm(p(1 + MOD(k, cyclen))), 1)
1087 164758 : CALL helium_pbc(helium, tmp1)
1088 164758 : ds = ds - x*(tmp1(1)*tmp1(1) + tmp1(2)*tmp1(2) + tmp1(3)*tmp1(3))
1089 4906761 : DO l = 1, helium%atoms
1090 4906761 : IF (l /= p(k)) THEN
1091 18308980 : rm1(:) = pos(:, p(k), j) - pos(:, l, j)
1092 18308980 : rm2(:) = pos(:, perm(p(k)), 1) - pos(:, perm(l), 1)
1093 4577245 : 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 329516 : IF (k < cyclen) THEN
1098 0 : DO l = k + 1, cyclen
1099 0 : rm1(:) = pos(:, p(k), j) - pos(:, p(l), j)
1100 0 : rm2(:) = pos(:, perm(p(k)), pk1) - pos(:, perm(p(l)), pk1)
1101 0 : ds = ds - helium_eval_expansion(helium, rm1, rm2, work3)
1102 : END DO
1103 : END IF
1104 : END DO
1105 164758 : IF (cyclen > 1) THEN
1106 : !k,c,l
1107 0 : c = perm(p(1))
1108 0 : DO k = 1, cyclen - 1
1109 0 : perm(p(k)) = perm(p(k + 1))
1110 : END DO
1111 0 : perm(p(cyclen)) = c
1112 : END IF
1113 329516 : DO k = 1, cyclen
1114 4906761 : DO l = 1, helium%atoms
1115 4906761 : IF (l /= p(k)) THEN
1116 18308980 : rm1(:) = work(:, p(k), j) - work(:, l, j)
1117 18308980 : rm2(:) = work(:, perm(p(k)), 1) - work(:, perm(l), 1)
1118 4577245 : 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 329516 : IF (k < cyclen) THEN
1123 0 : DO l = k + 1, cyclen
1124 0 : rm1(:) = work(:, p(k), j) - work(:, p(l), j)
1125 0 : rm2(:) = work(:, perm(p(k)), 1) - work(:, perm(p(l)), 1)
1126 0 : ds = ds + helium_eval_expansion(helium, rm1, rm2, work3)
1127 : END DO
1128 : END IF
1129 : END DO
1130 164758 : DEALLOCATE (work3)
1131 : ! ok now accept or reject:
1132 164758 : rtmp = helium%rng_stream_uniform%next()
1133 : ! IF ((dtk+ds-pds<0.0_dp).AND.(EXP(dtk+ds-pds)<rtmp)) THEN
1134 164758 : IF (dtk + ds - pds < 0.0_dp) THEN
1135 100337 : IF (EXP(dtk + ds - pds) < rtmp) THEN
1136 320841 : DO c = ifix, helium%beads
1137 606033 : DO k = 1, cyclen
1138 2281536 : work(:, p(k), c) = pos(:, p(k), c)
1139 : END DO
1140 : END DO
1141 35649 : IF (cyclen > 1) THEN
1142 0 : c = perm(p(cyclen))
1143 0 : DO k = cyclen - 1, 1, -1
1144 0 : perm(p(k + 1)) = perm(p(k))
1145 : END DO
1146 0 : perm(p(1)) = c
1147 : END IF
1148 35649 : 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 129109 : IF (.NOT. helium%periodic) THEN
1157 19322 : IF (helium%solute_present) THEN
1158 77288 : new_com(:) = helium%center(:)
1159 : ELSE
1160 0 : new_com(:) = 0.0_dp
1161 0 : DO k = 1, helium%atoms
1162 0 : DO l = 1, helium%beads
1163 0 : new_com(:) = new_com(:) + helium%work(:, k, l)
1164 : END DO
1165 : END DO
1166 0 : new_com(:) = new_com(:)/helium%atoms/helium%beads
1167 : END IF
1168 : ! Cycle through atoms (ignore connectivity)
1169 19322 : should_reject = .FALSE.
1170 114428 : outer: DO ia = 1, helium%atoms
1171 95539 : bis(:) = 0.0_dp
1172 1624163 : DO ib = 1, helium%beads
1173 6210035 : DO ic = 1, 3
1174 6114496 : bis(ic) = bis(ic) + work(ic, ia, ib) - new_com(ic)
1175 : END DO
1176 : END DO
1177 382156 : bis(:) = bis(:)/helium%beads
1178 382156 : rtmp = DOT_PRODUCT(bis, bis)
1179 114428 : IF (rtmp >= helium%droplet_radius**2) THEN
1180 433 : biso(:) = 0.0_dp
1181 7361 : DO ib = 1, helium%beads
1182 28145 : DO ic = 1, 3
1183 27712 : biso(ic) = biso(ic) + pos(ic, ia, ib) - helium%center(ic)
1184 : END DO
1185 : END DO
1186 1732 : biso(:) = biso(:)/helium%beads
1187 1732 : rtmpo = DOT_PRODUCT(biso, biso)
1188 : ! only reject if it moves away from COM outside the droplet radius
1189 433 : 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 19322 : IF (should_reject) THEN
1197 : ! restore work and perm, then return
1198 3897 : DO c = ifix, helium%beads
1199 7361 : DO k = 1, cyclen
1200 27712 : work(:, p(k), c) = pos(:, p(k), c)
1201 : END DO
1202 : END DO
1203 433 : IF (cyclen > 1) THEN
1204 0 : c = perm(p(cyclen))
1205 0 : DO k = cyclen - 1, 1, -1
1206 0 : perm(p(k + 1)) = perm(p(k))
1207 : END DO
1208 0 : perm(p(1)) = c
1209 : END IF
1210 433 : RETURN
1211 : END IF
1212 : END IF
1213 : ! accepted.
1214 :
1215 : ! copy trial over to the real thing
1216 1158084 : DO c = ifix, helium%beads
1217 2187492 : DO k = 1, cyclen
1218 8235264 : pos(:, p(k), c) = work(:, p(k), c)
1219 : END DO
1220 : END DO
1221 257352 : DO k = 1, cyclen
1222 257352 : helium%iperm(perm(p(k))) = p(k)
1223 : END DO
1224 128676 : helium%num_accepted(level + 2, cyclen) = helium%num_accepted(level + 2, cyclen) + 1
1225 128676 : res = .TRUE.
1226 :
1227 128676 : RETURN
1228 701096 : END FUNCTION helium_slice_metro_cyclic
1229 :
1230 : ! ***************************************************************************
1231 : !> \brief ...
1232 : !> \param helium ...
1233 : !> \param len ...
1234 : !> \return ...
1235 : !> \author hforbert
1236 : ! **************************************************************************************************
1237 990487 : 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 990487 : INTEGER, DIMENSION(:), POINTER :: iperm, p, perm
1244 990487 : INTEGER, DIMENSION(:, :), POINTER :: nmatrix
1245 : REAL(KIND=dp) :: rnd, s1, s2, t
1246 990487 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ipmatrix, pmatrix, tmatrix
1247 :
1248 990487 : s1 = 0.0_dp
1249 990487 : s2 = 0.0_dp
1250 990487 : res = .FALSE.
1251 990487 : n = helium%atoms
1252 990487 : tmatrix => helium%tmatrix
1253 990487 : pmatrix => helium%pmatrix
1254 990487 : ipmatrix => helium%ipmatrix
1255 990487 : perm => helium%permutation
1256 990487 : iperm => helium%iperm
1257 990487 : p => helium%ptable
1258 990487 : nmatrix => helium%nmatrix
1259 :
1260 990487 : p(len + 1) = -1
1261 1980974 : rnd = helium%rng_stream_uniform%next()
1262 990487 : p(1) = INT(n*rnd) + 1
1263 1028377 : DO k = 1, len - 1
1264 1015528 : t = helium%rng_stream_uniform%next()
1265 : ! find the corresponding path to connect to
1266 : ! using the precalculated optimal decision tree:
1267 1015528 : i = n - 1
1268 : DO
1269 1067754 : IF (tmatrix(p(k), i) > t) THEN
1270 51667 : i = nmatrix(p(k), 2*i - 1)
1271 : ELSE
1272 1016087 : i = nmatrix(p(k), 2*i)
1273 : END IF
1274 1067754 : IF (i < 0) EXIT
1275 : END DO
1276 1015528 : i = -i
1277 : ! which particle was it previously connected to?
1278 1015528 : p(k + 1) = iperm(i)
1279 : ! is it unique? quit if it was already part of the permutation
1280 1079057 : DO j = 1, k
1281 1079057 : IF (p(j) == p(k + 1)) RETURN
1282 : END DO
1283 : ! acummulate the needed values for the final
1284 : ! accept/reject step:
1285 37890 : s1 = s1 + ipmatrix(p(k), i)
1286 50739 : s2 = s2 + ipmatrix(p(k), perm(p(k)))
1287 : END DO
1288 : ! close the permutation loop:
1289 12849 : s1 = s1 + ipmatrix(p(len), perm(p(1)))
1290 12849 : s2 = s2 + ipmatrix(p(len), perm(p(len)))
1291 : ! final accept/reject:
1292 12849 : rnd = helium%rng_stream_uniform%next()
1293 12849 : t = s1*rnd
1294 12849 : 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 35 : s1 = pmatrix(p(len), perm(p(1))) - pmatrix(p(len), perm(p(len)))
1299 70 : DO k = 1, len - 1
1300 70 : s1 = s1 + pmatrix(p(k), perm(p(k + 1))) - pmatrix(p(k), perm(p(k)))
1301 : END DO
1302 35 : helium%pweight = s1
1303 35 : res = .TRUE.
1304 35 : RETURN
1305 990487 : END FUNCTION helium_select_permutation
1306 :
1307 : END MODULE helium_sampling
|
