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