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