Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief tree nodes acceptance
10 : !> code is separated in 3 parts, first the acceptance criteria,
11 : !> second the tree node acceptance handling, searching etc. and
12 : !> than the acceptance probability handling
13 : !> \par History
14 : !> 11.2012 created [Mandes Schoenherr]
15 : !> \author Mandes
16 : ! **************************************************************************************************
17 :
18 : MODULE tmc_tree_acceptance
19 : USE cp_log_handling, ONLY: cp_to_string
20 : USE kinds, ONLY: dp
21 : USE physcon, ONLY: boltzmann,&
22 : joule
23 : USE tmc_calculations, ONLY: compute_estimated_prob,&
24 : get_scaled_cell
25 : USE tmc_dot_tree, ONLY: create_dot_color,&
26 : create_global_tree_dot_color
27 : USE tmc_file_io, ONLY: write_result_list_element
28 : USE tmc_move_handle, ONLY: prob_update
29 : USE tmc_move_types, ONLY: mv_type_MD,&
30 : mv_type_volume_move
31 : USE tmc_stati, ONLY: task_type_gaussian_adaptation
32 : USE tmc_tree_build, ONLY: remove_unused_g_tree
33 : USE tmc_tree_search, ONLY: get_subtree_elements_to_check,&
34 : search_canceling_elements,&
35 : search_next_gt_element_to_check,&
36 : search_parent_element
37 : USE tmc_tree_types, ONLY: &
38 : add_to_list, global_tree_type, gt_elem_list_type, status_accepted, status_accepted_result, &
39 : status_calc_approx_ener, status_calculate_MD, status_calculate_NMC_steps, &
40 : status_calculate_energy, status_calculated, status_cancel_ener, status_cancel_nmc, &
41 : status_canceled_ener, status_canceled_nmc, status_created, status_deleted, &
42 : status_deleted_result, status_rejected, status_rejected_result, tree_type
43 : USE tmc_types, ONLY: tmc_env_type,&
44 : tmc_param_type
45 : #include "../base/base_uses.f90"
46 :
47 : IMPLICIT NONE
48 :
49 : PRIVATE
50 :
51 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_tree_acceptance'
52 :
53 : PUBLIC :: acceptance_check
54 : PUBLIC :: check_acceptance_of_depending_subtree_nodes, &
55 : check_elements_for_acc_prob_update
56 : PUBLIC :: tree_update
57 :
58 : INTEGER, PARAMETER :: DEBUG = 0
59 :
60 : CONTAINS
61 :
62 : !============================================================================
63 : ! acceptance criteria calculations
64 : !============================================================================
65 : ! **************************************************************************************************
66 : !> \brief standard Monte Carlo and 2 potential acceptance check
67 : !> acceptance check of move from old(last accepted) to new configuration
68 : !> the sum of kinetic and potential energy is used
69 : !> acc(o->n)=min(1,exp( -beta*(H(n)-H(o)) ))
70 : !> \param tree_element new/actual configuration
71 : !> \param parent_element last accepted configuration
72 : !> \param tmc_params TMC global parameters
73 : !> \param temperature actual temperature configuration should be checked with
74 : !> \param diff_pot_check 2potential check or not?
75 : !> \param accept result (configuration accepted of rejected)
76 : !> \param rnd_nr random number for acceptance check
77 : !> \param approx_ener for NMC the approximated energies schould be used
78 : !> \author Mandes 12.2012
79 : ! **************************************************************************************************
80 8802 : SUBROUTINE acceptance_check(tree_element, parent_element, tmc_params, &
81 : temperature, diff_pot_check, accept, rnd_nr, &
82 : approx_ener)
83 : TYPE(tree_type), POINTER :: tree_element, parent_element
84 : TYPE(tmc_param_type), POINTER :: tmc_params
85 : REAL(KIND=dp) :: temperature
86 : LOGICAL :: diff_pot_check, accept
87 : REAL(KIND=dp) :: rnd_nr
88 : LOGICAL, OPTIONAL :: approx_ener
89 :
90 : CHARACTER(LEN=*), PARAMETER :: routineN = 'acceptance_check'
91 :
92 : INTEGER :: handle
93 : REAL(KIND=dp) :: ekin_last_acc, elem_ener, kB, parent_ener
94 :
95 4401 : CPASSERT(ASSOCIATED(tree_element))
96 4401 : CPASSERT(ASSOCIATED(parent_element))
97 4401 : CPASSERT(ASSOCIATED(tmc_params))
98 4401 : CPASSERT(temperature > 0.0_dp)
99 4401 : CPASSERT(rnd_nr >= 0.0_dp .AND. rnd_nr <= 1.0_dp)
100 :
101 4401 : kB = boltzmann/joule
102 :
103 : ! start the timing
104 4401 : CALL timeset(routineN, handle)
105 :
106 4401 : IF (tmc_params%task_type == task_type_gaussian_adaptation) THEN
107 0 : CPABORT("")
108 : !TODO CALL acc_check_gauss_adapt(f=tree_element%potential, ct=tree_element%ekin_before_md, acc=accept)
109 : !DO NOT DO RETURN
110 : END IF
111 :
112 : !-- using 2 different potentials, the energy difference between these potentials
113 : ! and the two configurations have to be regarded.
114 : ! The ensamble should be in equilibrium state of the approximate potential
115 4401 : IF (diff_pot_check .AND. (tmc_params%NMC_inp_file /= "")) THEN
116 66 : IF ((tree_element%potential == HUGE(tree_element%potential)) .OR. &
117 : tree_element%e_pot_approx == HUGE(tree_element%e_pot_approx)) THEN
118 : elem_ener = HUGE(elem_ener)
119 : ELSE
120 : !for different potentials we have to regard the differences in energy
121 : ! min(1,e^{-\beta*[(E_{exact}(n)-E_{approx}(n))-(E_{exact}(o)-E_{approx}(o))]})
122 : elem_ener = 1.0_dp/(kB*temperature)*tree_element%potential &
123 : - 1.0_dp/(kB*tmc_params%Temp(tree_element%temp_created)) &
124 66 : *tree_element%e_pot_approx
125 : END IF
126 : parent_ener = 1.0_dp/(kB*temperature)*parent_element%potential &
127 : - 1.0_dp/(kB*tmc_params%Temp(tree_element%temp_created)) &
128 66 : *parent_element%e_pot_approx
129 :
130 : !-- always accepted if new energy is smaller than old energy
131 66 : IF (elem_ener <= parent_ener) THEN
132 27 : accept = .TRUE.
133 : ELSE
134 : !-- gaussian distributed acceptance if new energy is greater than old energy
135 39 : IF (rnd_nr < &
136 : EXP(-(elem_ener - parent_ener))) THEN
137 8 : accept = .TRUE.
138 : ELSE
139 31 : accept = .FALSE.
140 : END IF
141 : END IF
142 : ELSE
143 : !-- standard MC check, but also using the kinetic energy for Hybrid Monte Carlo
144 4335 : IF (tree_element%move_type == mv_type_MD) THEN
145 0 : ekin_last_acc = tree_element%ekin_before_md
146 : ! using the kinetic energy before velocity change
147 : ELSE
148 4335 : ekin_last_acc = parent_element%ekin
149 : END IF
150 : ! comparing aproximated energies
151 4335 : IF (PRESENT(approx_ener)) THEN
152 : elem_ener = tree_element%e_pot_approx &
153 126 : + tree_element%ekin
154 : parent_ener = parent_element%e_pot_approx &
155 126 : + ekin_last_acc
156 : ELSE
157 : elem_ener = tree_element%potential &
158 4209 : + tree_element%ekin
159 : parent_ener = parent_element%potential &
160 4209 : + ekin_last_acc
161 : END IF
162 :
163 : !-- always accepted if new energy is smaller than old energy
164 4335 : IF (elem_ener <= parent_ener) THEN
165 418 : accept = .TRUE.
166 : ELSE
167 : !-- gaussian distributed acceptance if new energy is greater than old energy
168 3917 : IF (rnd_nr < &
169 : EXP(-1.0_dp/(kB*temperature)*(elem_ener - parent_ener))) THEN
170 229 : accept = .TRUE.
171 : ELSE
172 3688 : accept = .FALSE.
173 : END IF
174 : END IF
175 : END IF
176 :
177 : ! update the estimated energy acceptance probability distribution
178 4401 : IF (diff_pot_check) THEN
179 4275 : CPASSERT(ASSOCIATED(tmc_params%prior_NMC_acc))
180 4275 : tmc_params%prior_NMC_acc%counter = tmc_params%prior_NMC_acc%counter + 1
181 : tmc_params%prior_NMC_acc%aver = (tmc_params%prior_NMC_acc%aver*(tmc_params%prior_NMC_acc%counter - 1) + &
182 4275 : ((elem_ener - parent_ener)))/REAL(tmc_params%prior_NMC_acc%counter, KIND=dp)
183 : tmc_params%prior_NMC_acc%aver_2 = (tmc_params%prior_NMC_acc%aver_2*(tmc_params%prior_NMC_acc%counter - 1) + &
184 4275 : (elem_ener - parent_ener)**2)/REAL(tmc_params%prior_NMC_acc%counter, KIND=dp)
185 : END IF
186 :
187 : ! end the timing
188 4401 : CALL timestop(handle)
189 4401 : END SUBROUTINE acceptance_check
190 :
191 : ! **************************************************************************************************
192 : !> \brief parallel tempering swap acceptance check,
193 : !> check swapped configurations of different temperatures
194 : !> using sum of potential and kinetic energy
195 : !> always accepted if energy of configuration with higher temperature
196 : !> is smaller than energy of configuration with lower temperature
197 : !> acc(o->n)=min(1,exp((beta_i-beta_j)*(U_i-U_j)))
198 : !> \param tree_elem actual global tree element
199 : !> \param conf1 sub tree element of lower temperature
200 : !> \param conf2 sub tree element of higher temperature
201 : !> \param tmc_params TMC environment parameters
202 : !> \param accept acceptance of configurational change
203 : !> \author Mandes 12.2012
204 : ! **************************************************************************************************
205 336 : SUBROUTINE swap_acceptance_check(tree_elem, conf1, conf2, tmc_params, accept)
206 : TYPE(global_tree_type), POINTER :: tree_elem
207 : TYPE(tree_type), POINTER :: conf1, conf2
208 : TYPE(tmc_param_type), POINTER :: tmc_params
209 : LOGICAL :: accept
210 :
211 : CHARACTER(LEN=*), PARAMETER :: routineN = 'swap_acceptance_check'
212 :
213 : INTEGER :: handle
214 : REAL(KIND=dp) :: delta, elem1_ener, elem2_ener, kB, vol1, &
215 : vol2
216 :
217 168 : kB = boltzmann/joule
218 :
219 168 : CPASSERT(ASSOCIATED(tree_elem))
220 168 : CPASSERT(tree_elem%rnd_nr >= 0.0_dp)
221 168 : CPASSERT(ASSOCIATED(conf1))
222 168 : CPASSERT(ASSOCIATED(conf2))
223 168 : CPASSERT(ASSOCIATED(tmc_params))
224 :
225 : ! start the timing
226 168 : CALL timeset(routineN, handle)
227 :
228 168 : IF (tmc_params%pressure > 0.0_dp) THEN
229 : ! pt-NVT
230 : elem1_ener = conf1%potential &
231 18 : + conf1%ekin
232 : elem2_ener = conf2%potential &
233 18 : + conf2%ekin
234 : ! the swap is done with prob: exp((\beta_i-\beta_j)(U_i-U_j)),
235 : ! BUT because they are already swaped we exchange the energies.
236 18 : IF (tree_elem%rnd_nr < EXP((1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf)) - &
237 : 1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf + 1))) &
238 : *(elem2_ener - elem1_ener))) THEN
239 18 : accept = .TRUE.
240 : ELSE
241 0 : accept = .FALSE.
242 : END IF
243 : ELSE
244 : ! pt-NpT (parallel Tempering with constant pressure, temperature and num. particle)
245 : CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf1%box_scale, &
246 150 : vol=vol1)
247 : CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf2%box_scale, &
248 150 : vol=vol2)
249 : ! delta= (beta_m-beta_n)(U_n-U_m) + (beta_m*P_m-beta_n*P_n)(V_n-V_m)
250 : delta = (1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf)) &
251 : - 1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf + 1)))* &
252 : ((conf2%potential + conf2%ekin) - (conf1%potential + conf1%ekin)) + &
253 : (1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf))*tmc_params%pressure &
254 : - 1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf + 1))*tmc_params%pressure)* &
255 150 : (vol2 - vol1)
256 150 : IF (tree_elem%rnd_nr < EXP(delta)) THEN
257 112 : accept = .TRUE.
258 : ELSE
259 38 : accept = .FALSE.
260 : END IF
261 : END IF
262 : ! end the timing
263 168 : CALL timestop(handle)
264 168 : END SUBROUTINE swap_acceptance_check
265 :
266 : ! **************************************************************************************************
267 : !> \brief volume move acceptance check
268 : !> sampling NPT, we need to do volume moves. Their acceptance are
269 : !> checked regarding the old and new energy, the volume difference
270 : !> and the constant pressure
271 : !> \param parent_elem old tree element with old box size
272 : !> \param new_elem actual tree element with new box size
273 : !> \param temperature actual temperature
274 : !> \param rnd_nr random number to check with
275 : !> \param tmc_params TMC environment parameters
276 : !> \param accept Monte Carlo move acceptance (result)
277 : !> \author Mandes 12.2012
278 : ! **************************************************************************************************
279 340 : SUBROUTINE volume_acceptance_check(parent_elem, new_elem, temperature, &
280 : rnd_nr, tmc_params, accept)
281 : TYPE(tree_type), POINTER :: parent_elem, new_elem
282 : REAL(KIND=dp) :: temperature, rnd_nr
283 : TYPE(tmc_param_type), POINTER :: tmc_params
284 : LOGICAL :: accept
285 :
286 : CHARACTER(LEN=*), PARAMETER :: routineN = 'volume_acceptance_check'
287 :
288 : INTEGER :: handle
289 : REAL(KIND=dp) :: d_enthalpy, kB, n_vol, p_vol
290 :
291 170 : kB = boltzmann/joule
292 :
293 170 : CPASSERT(ASSOCIATED(parent_elem))
294 170 : CPASSERT(ASSOCIATED(new_elem))
295 170 : CPASSERT(temperature > 0.0_dp)
296 170 : CPASSERT(rnd_nr > 0.0_dp)
297 170 : CPASSERT(ASSOCIATED(tmc_params))
298 170 : CPASSERT(tmc_params%pressure >= 0.0_dp)
299 :
300 : ! start the timing
301 170 : CALL timeset(routineN, handle)
302 :
303 : CALL get_scaled_cell(cell=tmc_params%cell, box_scale=parent_elem%box_scale, &
304 170 : vol=p_vol)
305 : CALL get_scaled_cell(cell=tmc_params%cell, box_scale=new_elem%box_scale, &
306 170 : vol=n_vol)
307 : ! H=enthalpy, U=energy, P=pressure, V=volume, kB=Boltzmann constant, T=temperature, N=amount particle
308 : ! delta_H = delta_U + P*delta_V - kB*T*N*ln(V_n/V_p)
309 : IF (.FALSE.) THEN
310 : ! the volume move in volume space (dV)
311 : d_enthalpy = (new_elem%potential - parent_elem%potential) + &
312 : tmc_params%pressure*(n_vol - p_vol) - &
313 : kB*temperature*(SIZE(new_elem%pos)/ &
314 : tmc_params%dim_per_elem)* &
315 : LOG(n_vol/p_vol)
316 : ELSE
317 170 : IF (tmc_params%v_isotropic) THEN
318 : d_enthalpy = (new_elem%potential - parent_elem%potential) + &
319 : tmc_params%pressure*(n_vol - p_vol) - &
320 : kB*temperature*((SIZE(new_elem%pos)/ &
321 : tmc_params%dim_per_elem) + 2/REAL(3, KIND=dp))* &
322 170 : LOG(n_vol/p_vol)
323 : ELSE
324 : d_enthalpy = (new_elem%potential - parent_elem%potential) + &
325 : tmc_params%pressure*(n_vol - p_vol) - &
326 : kB*temperature*(SIZE(new_elem%pos)/ &
327 : tmc_params%dim_per_elem)* &
328 0 : LOG(n_vol/p_vol)
329 : END IF
330 : END IF
331 : ! acc(o->n) = min(1, exp(-beta*delta_H))
332 170 : IF (d_enthalpy <= 0.0_dp) THEN
333 53 : accept = .TRUE.
334 : ELSE
335 117 : IF (rnd_nr < EXP(-1.0_dp/(kB*temperature)*d_enthalpy)) THEN
336 9 : accept = .TRUE.
337 : ELSE
338 108 : accept = .FALSE.
339 : END IF
340 : END IF
341 : ! end the timing
342 170 : CALL timestop(handle)
343 170 : END SUBROUTINE volume_acceptance_check
344 :
345 : !============================================================================
346 : ! tree node acceptance
347 : !============================================================================
348 : ! **************************************************************************************************
349 : !> \brief check acceptance of energy calculated element and related childs,
350 : !> when ready
351 : !> \param tree_elem actual tree element with calculated energy
352 : !> \param tmc_env TMC environment parameters
353 : !> \author Mandes 12.2012
354 : ! **************************************************************************************************
355 8890 : SUBROUTINE check_acceptance_of_depending_subtree_nodes(tree_elem, tmc_env)
356 : TYPE(tree_type), POINTER :: tree_elem
357 : TYPE(tmc_env_type), POINTER :: tmc_env
358 :
359 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_acceptance_of_depending_subtree_nodes'
360 :
361 : INTEGER :: handle
362 : LOGICAL :: parent_ready
363 : TYPE(tree_type), POINTER :: act_elem, parent_elem
364 :
365 4445 : NULLIFY (parent_elem, act_elem)
366 :
367 4445 : CPASSERT(ASSOCIATED(tmc_env))
368 4445 : CPASSERT(ASSOCIATED(tree_elem))
369 4445 : CPASSERT(tree_elem%stat == status_calculated)
370 4445 : CPASSERT(ASSOCIATED(tree_elem%parent))
371 :
372 : ! start the timing
373 4445 : CALL timeset(routineN, handle)
374 :
375 4445 : act_elem => tree_elem
376 :
377 : ! ------------------------------------------------------
378 : ! check this element
379 4445 : parent_elem => search_parent_element(act_elem)
380 4445 : CPASSERT(.NOT. ASSOCIATED(act_elem, parent_elem))
381 :
382 : ! check status of parent element
383 4445 : SELECT CASE (parent_elem%stat)
384 : CASE (status_created, status_calculate_energy, status_calculate_MD, &
385 : status_calculate_NMC_steps, status_canceled_ener, &
386 : status_canceled_nmc, status_cancel_nmc, status_cancel_ener)
387 : parent_ready = .FALSE.
388 : CASE (status_accepted_result, status_rejected_result, &
389 : status_accepted, status_rejected, status_calculated)
390 0 : parent_ready = .TRUE.
391 : CASE DEFAULT
392 4445 : CPABORT("unknown parent element status"//cp_to_string(parent_elem%stat))
393 : END SELECT
394 :
395 : ! ready ?
396 0 : IF (parent_ready) THEN
397 : ! acceptance check
398 : CALL check_and_change_status_of_subtree_elem(act_elem=act_elem, &
399 4445 : parent_elem=parent_elem, tmc_env=tmc_env)
400 : END IF
401 : !------------------------------------------------------
402 : ! check acc child
403 4445 : parent_elem => tree_elem
404 4445 : IF (ASSOCIATED(tree_elem%acc)) THEN
405 0 : act_elem => tree_elem%acc
406 0 : IF (act_elem%stat == status_calculated) THEN
407 : CALL check_and_change_status_of_subtree_elem(act_elem=act_elem, &
408 0 : parent_elem=parent_elem, tmc_env=tmc_env)
409 : END IF
410 : !------------------------------------------------------
411 : ! check all nacc childs
412 : nacc_loop: DO
413 0 : IF (.NOT. ASSOCIATED(act_elem%nacc)) EXIT nacc_loop
414 0 : act_elem => act_elem%nacc
415 0 : IF (act_elem%stat == status_calculated) THEN
416 : CALL check_and_change_status_of_subtree_elem(act_elem=act_elem, &
417 0 : parent_elem=parent_elem, tmc_env=tmc_env)
418 : END IF
419 : END DO nacc_loop
420 : END IF
421 :
422 : ! end the timing
423 4445 : CALL timestop(handle)
424 4445 : END SUBROUTINE check_acceptance_of_depending_subtree_nodes
425 :
426 : ! **************************************************************************************************
427 : !> \brief checking the elements and change the status
428 : !> \param act_elem actual tree element (new configuration)
429 : !> \param parent_elem parent tree element (old configuration)
430 : !> \param tmc_env TMC environment parameters
431 : !> \author Mandes 12.2012
432 : ! **************************************************************************************************
433 4445 : SUBROUTINE check_and_change_status_of_subtree_elem(act_elem, parent_elem, tmc_env)
434 : TYPE(tree_type), POINTER :: act_elem, parent_elem
435 : TYPE(tmc_env_type), POINTER :: tmc_env
436 :
437 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_and_change_status_of_subtree_elem'
438 :
439 : INTEGER :: handle
440 : LOGICAL :: flag, result_acc
441 : TYPE(global_tree_type), POINTER :: tmp_gt_ptr
442 : TYPE(gt_elem_list_type), POINTER :: tmp_gt_list_ptr
443 :
444 4445 : NULLIFY (tmp_gt_list_ptr, tmp_gt_ptr)
445 :
446 : ! start the timing
447 4445 : CALL timeset(routineN, handle)
448 :
449 4445 : CPASSERT(ASSOCIATED(tmc_env))
450 4445 : CPASSERT(ASSOCIATED(tmc_env%params))
451 4445 : CPASSERT(ASSOCIATED(act_elem))
452 4445 : CPASSERT(ASSOCIATED(parent_elem))
453 4445 : CPASSERT(act_elem%stat == status_calculated)
454 : MARK_USED(parent_elem)
455 :
456 4445 : flag = .FALSE.
457 :
458 4445 : tmp_gt_list_ptr => act_elem%gt_nodes_references
459 : ! check for all global tree elements referring to this subtree element
460 : ! same subtree element can be accepted at a certain temperature and rejected at another one
461 8890 : DO WHILE (ASSOCIATED(tmp_gt_list_ptr))
462 4445 : CPASSERT(tmp_gt_list_ptr%gt_elem%stat /= status_accepted_result)
463 4445 : CPASSERT(tmp_gt_list_ptr%gt_elem%stat /= status_rejected_result)
464 :
465 : CALL check_elements(gt_act_elem=tmp_gt_list_ptr%gt_elem, tmc_env=tmc_env, &
466 4445 : check_done=flag, result_acc=result_acc)
467 4445 : IF (flag) THEN
468 : ! probability update
469 : CALL prob_update(move_types=tmc_env%params%move_types, elem=act_elem, &
470 4445 : acc=result_acc, prob_opt=tmc_env%params%esimate_acc_prob)
471 :
472 : ! change status
473 : ! accepted case: delete (and cancel) not accepted branch
474 4445 : IF (result_acc .EQV. .TRUE.) THEN
475 684 : tmp_gt_list_ptr%gt_elem%stat = status_accepted
476 684 : IF (ASSOCIATED(tmp_gt_list_ptr%gt_elem%nacc)) THEN
477 0 : tmp_gt_ptr => tmp_gt_list_ptr%gt_elem%nacc
478 : CALL remove_unused_g_tree(begin_ptr=tmp_gt_ptr, &
479 : end_ptr=tmp_gt_list_ptr%gt_elem, removed=flag, &
480 0 : tmc_env=tmc_env)
481 : END IF
482 : ELSE
483 : ! not accepted case: delete (and cancel) accepted branch
484 3761 : tmp_gt_list_ptr%gt_elem%stat = status_rejected
485 3761 : IF (ASSOCIATED(tmp_gt_list_ptr%gt_elem%acc)) THEN
486 0 : tmp_gt_ptr => tmp_gt_list_ptr%gt_elem%acc
487 : CALL remove_unused_g_tree(begin_ptr=tmp_gt_ptr, &
488 : end_ptr=tmp_gt_list_ptr%gt_elem, removed=flag, &
489 0 : tmc_env=tmc_env)
490 : END IF
491 : END IF
492 :
493 4445 : IF (tmc_env%params%DRAW_TREE) &
494 : CALL create_global_tree_dot_color(gt_tree_element=tmp_gt_list_ptr%gt_elem, &
495 36 : tmc_params=tmc_env%params)
496 : END IF
497 4445 : tmp_gt_list_ptr => tmp_gt_list_ptr%next
498 : END DO
499 : ! end the timing
500 4445 : CALL timestop(handle)
501 4445 : END SUBROUTINE check_and_change_status_of_subtree_elem
502 :
503 : ! **************************************************************************************************
504 : !> \brief change status flag of all subtree elements
505 : !> most important to draw the right subtrees
506 : !> \param gt_ptr global tree pointer, the related/changed sub tree element
507 : !> should change status
508 : !> \param stat the status of the global tree element
509 : !> \param tmc_params TMC emvironment parameters for drawing
510 : !> \author Mandes 12.2012
511 : ! **************************************************************************************************
512 8890 : SUBROUTINE subtree_configuration_stat_change(gt_ptr, stat, tmc_params)
513 : TYPE(global_tree_type), POINTER :: gt_ptr
514 : INTEGER :: stat
515 : TYPE(tmc_param_type), POINTER :: tmc_params
516 :
517 : CHARACTER(LEN=*), PARAMETER :: routineN = 'subtree_configuration_stat_change'
518 :
519 : INTEGER :: handle
520 : TYPE(tree_type), POINTER :: ptr
521 :
522 4445 : NULLIFY (ptr)
523 :
524 4445 : CPASSERT(ASSOCIATED(gt_ptr))
525 4445 : CPASSERT(ASSOCIATED(tmc_params))
526 :
527 : ! start the timing
528 4445 : CALL timeset(routineN, handle)
529 :
530 : ! don't respect swaped nodes (subtree nodes are already in the list,
531 : ! and we dont want them in marked in subtrees)
532 4445 : IF (.NOT. gt_ptr%swaped) THEN
533 4445 : ptr => gt_ptr%conf(gt_ptr%mv_conf)%elem
534 : ! dependent on the status (don't map each status 1:1,
535 : ! e.g. not the result ones)
536 5129 : SELECT CASE (stat)
537 : CASE (status_accepted_result)
538 684 : ptr%stat = status_accepted
539 : CASE (status_rejected_result)
540 3761 : ptr%stat = status_rejected
541 : CASE (status_accepted, status_rejected)
542 0 : ptr%stat = stat
543 : CASE DEFAULT
544 : CALL cp_abort(__LOCATION__, &
545 : "unknown global tree status"// &
546 4445 : cp_to_string(stat))
547 : END SELECT
548 :
549 4445 : IF (tmc_params%DRAW_TREE) &
550 36 : CALL create_dot_color(tree_element=ptr, tmc_params=tmc_params)
551 : END IF
552 :
553 : ! end the timing
554 4445 : CALL timestop(handle)
555 4445 : END SUBROUTINE subtree_configuration_stat_change
556 :
557 : ! **************************************************************************************************
558 : !> \brief checks if the element is ready for an acceptance check
559 : !> \param elem sub tree element
560 : !> \param ready return value
561 : !> \author Mandes 12.2012
562 : ! **************************************************************************************************
563 1113548 : SUBROUTINE elem_ready_to_check(elem, ready)
564 : TYPE(tree_type), POINTER :: elem
565 : LOGICAL :: ready
566 :
567 : CHARACTER(LEN=*), PARAMETER :: routineN = 'elem_ready_to_check'
568 :
569 : INTEGER :: handle
570 :
571 556774 : CPASSERT(ASSOCIATED(elem))
572 :
573 : ! start the timing
574 556774 : CALL timeset(routineN, handle)
575 :
576 556774 : ready = .FALSE.
577 :
578 556774 : SELECT CASE (elem%stat)
579 : CASE (status_created, status_calculate_energy, &
580 : status_calculate_MD, status_calculate_NMC_steps, status_calc_approx_ener, &
581 : status_cancel_nmc, status_cancel_ener, status_canceled_nmc, status_canceled_ener)
582 266766 : ready = .FALSE.
583 : CASE (status_calculated, status_accepted_result, status_accepted, status_rejected, status_rejected_result)
584 266766 : ready = .TRUE.
585 : CASE DEFAULT
586 : CALL cp_abort(__LOCATION__, &
587 : "status of actual tree node is unknown" &
588 556774 : //cp_to_string(elem%stat))
589 : END SELECT
590 : ! end the timing
591 556774 : CALL timestop(handle)
592 556774 : END SUBROUTINE elem_ready_to_check
593 :
594 : ! **************************************************************************************************
595 : !> \brief checking the elements (find element and do acceptance check)
596 : !> \param gt_act_elem actual global tree element
597 : !> \param tmc_env TMC environment
598 : !> \param check_done successful checked? (result)
599 : !> \param result_acc checked configuration accepted? (result)
600 : !> \author Mandes 12.2012
601 : ! **************************************************************************************************
602 565664 : SUBROUTINE check_elements(gt_act_elem, tmc_env, check_done, result_acc)
603 : TYPE(global_tree_type), POINTER :: gt_act_elem
604 : TYPE(tmc_env_type), POINTER :: tmc_env
605 : LOGICAL :: check_done, result_acc
606 :
607 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_elements'
608 :
609 : INTEGER :: handle
610 : LOGICAL :: act_ready, tmp_ready
611 : TYPE(tree_type), POINTER :: act_element, tmp_element
612 :
613 282832 : NULLIFY (act_element, tmp_element)
614 282832 : check_done = .FALSE.
615 :
616 282832 : CPASSERT(ASSOCIATED(gt_act_elem))
617 282832 : CPASSERT(ASSOCIATED(tmc_env))
618 282832 : CPASSERT(ASSOCIATED(tmc_env%m_env))
619 :
620 : ! start the timing
621 282832 : CALL timeset(routineN, handle)
622 :
623 : !====================================================================================
624 : ! check if global tree element is already checked
625 287277 : SELECT CASE (gt_act_elem%stat)
626 : CASE (status_accepted, status_rejected)
627 4445 : check_done = .TRUE.
628 282832 : IF (gt_act_elem%stat == status_accepted) THEN
629 684 : result_acc = .TRUE.
630 3761 : ELSE IF (gt_act_elem%stat == status_rejected) THEN
631 3761 : result_acc = .FALSE.
632 : ELSE
633 : CALL cp_abort(__LOCATION__, &
634 : "undefinite status of already checked elements:" &
635 0 : //cp_to_string(gt_act_elem%stat))
636 : END IF
637 : CASE DEFAULT
638 : END SELECT
639 :
640 282832 : IF (.NOT. check_done) THEN
641 : !====================================================================================
642 : ! get elements related to global tree element to check
643 : CALL get_subtree_elements_to_check(gt_act_elem=gt_act_elem, elem1=act_element, &
644 278387 : elem2=tmp_element)
645 278387 : CPASSERT(ASSOCIATED(act_element))
646 278387 : CPASSERT(ASSOCIATED(tmp_element))
647 :
648 : ! check status of both elements (if they are already calculated and hence ready to check)
649 278387 : CALL elem_ready_to_check(elem=act_element, ready=act_ready)
650 278387 : CALL elem_ready_to_check(elem=tmp_element, ready=tmp_ready)
651 :
652 : ! both ready ? check
653 278387 : IF (act_ready .AND. tmp_ready) THEN
654 : ! 2 temperature (swap) check
655 4613 : IF (gt_act_elem%swaped) THEN
656 : CALL swap_acceptance_check(tree_elem=gt_act_elem, conf1=act_element, &
657 : conf2=tmp_element, tmc_params=tmc_env%params, &
658 168 : accept=result_acc)
659 : ! volume move check
660 4445 : ELSE IF (act_element%move_type == mv_type_volume_move) THEN
661 : CALL volume_acceptance_check(parent_elem=tmp_element, new_elem=act_element, &
662 : temperature=tmc_env%params%Temp(gt_act_elem%mv_conf), &
663 : rnd_nr=gt_act_elem%rnd_nr, &
664 170 : tmc_params=tmc_env%params, accept=result_acc)
665 : ELSE
666 4275 : IF (tmc_env%m_env%temp_decrease /= 1.0_dp) THEN
667 : CALL acceptance_check(tree_element=act_element, parent_element=tmp_element, &
668 : tmc_params=tmc_env%params, temperature=gt_act_elem%Temp, &
669 : diff_pot_check=.TRUE., accept=result_acc, &
670 0 : rnd_nr=gt_act_elem%rnd_nr)
671 : ELSE
672 : CALL acceptance_check(tree_element=act_element, parent_element=tmp_element, &
673 : tmc_params=tmc_env%params, &
674 : temperature=tmc_env%params%Temp(gt_act_elem%mv_conf), &
675 : diff_pot_check=.TRUE., accept=result_acc, &
676 4275 : rnd_nr=gt_act_elem%rnd_nr)
677 : END IF
678 : END IF
679 4613 : check_done = .TRUE.
680 : ELSE
681 : ! if element is not ready, update status of global tree element
682 : !TODO maybe update this part (how much status we need in global tree from sub tree??)
683 547548 : SELECT CASE (gt_act_elem%stat)
684 : !-- check if at least the MD move is already done
685 : !-- hence new configurations can be created on basis of this configuration
686 : CASE (status_calculate_MD, status_calculate_NMC_steps, status_calculate_energy, &
687 : status_created, status_calc_approx_ener)
688 273774 : IF (gt_act_elem%conf(gt_act_elem%mv_conf)%elem%stat /= gt_act_elem%stat) &
689 4477 : gt_act_elem%stat = gt_act_elem%conf(gt_act_elem%mv_conf)%elem%stat
690 : CASE (status_calculated)
691 : CASE DEFAULT
692 : CALL cp_abort(__LOCATION__, &
693 : "status of actual checked node is unknown" &
694 273774 : //cp_to_string(gt_act_elem%stat))
695 : END SELECT
696 : END IF
697 : END IF
698 : ! end the timing
699 282832 : CALL timestop(handle)
700 282832 : END SUBROUTINE check_elements
701 :
702 : ! **************************************************************************************************
703 : !> \brief searching tree nodes to check for Markov Chain, elements are marked
704 : !> and stored in lists ... (main entry point)
705 : !> \param tmc_env TMC environment
706 : !> \param result_acc checked configuration accepted? (result)
707 : !> \param something_updated ...
708 : !> \param
709 : !> \param
710 : !> \author Mandes 12.2012
711 : ! **************************************************************************************************
712 556802 : SUBROUTINE tree_update(tmc_env, result_acc, something_updated)
713 : TYPE(tmc_env_type), POINTER :: tmc_env
714 : LOGICAL :: result_acc, something_updated
715 :
716 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tree_update'
717 :
718 : INTEGER :: handle, itmp
719 : LOGICAL :: found
720 : TYPE(global_tree_type), POINTER :: gt_act_elem
721 : TYPE(tree_type), POINTER :: act_element, tmp_element
722 :
723 278401 : NULLIFY (gt_act_elem, tmp_element, act_element)
724 :
725 278401 : CPASSERT(ASSOCIATED(tmc_env))
726 :
727 : ! start the timing
728 278401 : CALL timeset(routineN, handle)
729 :
730 278401 : result_acc = .FALSE.
731 278401 : something_updated = .FALSE.
732 :
733 278401 : gt_act_elem => tmc_env%m_env%gt_act
734 : search_calculated_element_loop: DO
735 283014 : NULLIFY (act_element, tmp_element)
736 : ! search for element to check
737 283014 : CALL search_next_gt_element_to_check(ptr=gt_act_elem, found=found)
738 283014 : IF (.NOT. found) EXIT search_calculated_element_loop
739 :
740 : ! check the elements status end, if possible, do acceptance check
741 : CALL check_elements(gt_act_elem=gt_act_elem, tmc_env=tmc_env, &
742 278387 : check_done=found, result_acc=result_acc)
743 : ! if check was not possible, update the status of the global tree element and return
744 278387 : IF (.NOT. found) EXIT search_calculated_element_loop
745 :
746 : ! get elements related to global tree element, which were checked
747 : CALL get_subtree_elements_to_check(gt_act_elem=gt_act_elem, &
748 4613 : elem1=act_element, elem2=tmp_element)
749 :
750 : !========================================================================
751 : !-- set result counters
752 : ! counter for certain temperature
753 : tmc_env%m_env%result_count(gt_act_elem%mv_conf) = &
754 4613 : tmc_env%m_env%result_count(gt_act_elem%mv_conf) + 1
755 : ! in case of swapped also count the result for
756 : ! the other swapped temperature
757 4613 : IF (gt_act_elem%swaped) &
758 : tmc_env%m_env%result_count(gt_act_elem%mv_conf + 1) = &
759 168 : tmc_env%m_env%result_count(gt_act_elem%mv_conf + 1) + 1
760 : ! count also for global tree Markov Chain
761 4613 : tmc_env%m_env%result_count(0) = tmc_env%m_env%result_count(0) + 1
762 :
763 : ! flag for doing tree cleaning with canceling and deallocation
764 4613 : something_updated = .TRUE.
765 :
766 : !IF((.NOT.gt_act_elem%swaped) .AND. act_element%temp_created/=0 .AND. &
767 : ! act_element%temp_created/=gt_act_elem%mv_conf)&
768 : ! count_wrong_temp_move(gt_act_elem%mv_conf) = &
769 : ! count_wrong_temp_move(gt_act_elem%mv_conf)+1
770 :
771 : !========================================================================
772 : !-- sort element in lists
773 : !-- case NOT ACCEPTED
774 4613 : IF (.NOT. result_acc) THEN !result NOT accepted
775 3799 : IF (gt_act_elem%prob_acc > 0.0_dp) THEN
776 3771 : IF (gt_act_elem%prob_acc >= 0.5_dp) THEN
777 : ! wrong estimate (estimated accepted)
778 197 : tmc_env%m_env%estim_corr_wrong(4) = tmc_env%m_env%estim_corr_wrong(4) + 1
779 : IF (DEBUG > 0) WRITE (tmc_env%m_env%io_unit, *) &
780 : "Wrong guess for NACC (elem/estim acc prob)", gt_act_elem%nr, gt_act_elem%prob_acc
781 : ELSE
782 3574 : tmc_env%m_env%estim_corr_wrong(3) = tmc_env%m_env%estim_corr_wrong(3) + 1
783 : END IF
784 : END IF
785 3799 : gt_act_elem%stat = status_rejected_result
786 3799 : gt_act_elem%prob_acc = 0.0_dp
787 3799 : itmp = gt_act_elem%conf(gt_act_elem%mv_conf)%elem%sub_tree_nr
788 :
789 : !-- case ACCEPTED
790 : ELSE
791 814 : IF (gt_act_elem%prob_acc > 0.0_dp) THEN
792 799 : IF (gt_act_elem%prob_acc <= 0.5_dp) THEN
793 : ! wrong estimate (estimated NOT accepted)
794 347 : tmc_env%m_env%estim_corr_wrong(2) = tmc_env%m_env%estim_corr_wrong(2) + 1
795 : IF (DEBUG > 0) WRITE (tmc_env%m_env%io_unit, *) &
796 : "wrong guess for ACC (elem/estim acc prob)", gt_act_elem%nr, gt_act_elem%prob_acc
797 : ELSE
798 452 : tmc_env%m_env%estim_corr_wrong(1) = tmc_env%m_env%estim_corr_wrong(1) + 1
799 : END IF
800 : END IF
801 814 : gt_act_elem%stat = status_accepted_result
802 814 : gt_act_elem%prob_acc = 1.0_dp
803 :
804 : ! save result in the result list
805 814 : IF (.NOT. gt_act_elem%swaped) THEN
806 : !-- saving moved/changed configuration of result node
807 : tmc_env%m_env%result_list(gt_act_elem%mv_conf)%elem => &
808 684 : gt_act_elem%conf(gt_act_elem%mv_conf)%elem
809 : ELSE
810 : ! ATTENTION: act_element != gt_act_elem%conf(gt_act_elem%mv_conf),
811 : ! because we take the last accepted conf
812 130 : tmc_env%m_env%result_list(gt_act_elem%mv_conf)%elem => act_element
813 130 : tmc_env%m_env%result_list(gt_act_elem%mv_conf + 1)%elem => tmp_element
814 : END IF
815 814 : tmc_env%m_env%gt_act => gt_act_elem
816 : END IF ! result acceptance check
817 : !======================================================================
818 : !-- set status accepted or rejected of certain subtree elements
819 4613 : IF (.NOT. gt_act_elem%swaped) &
820 : CALL subtree_configuration_stat_change(gt_ptr=gt_act_elem, &
821 : stat=gt_act_elem%stat, &
822 4445 : tmc_params=tmc_env%params)
823 :
824 4613 : IF (tmc_env%params%DRAW_TREE) &
825 : CALL create_global_tree_dot_color(gt_tree_element=gt_act_elem, &
826 74 : tmc_params=tmc_env%params)
827 :
828 : ! probability update
829 : CALL prob_update(move_types=tmc_env%params%move_types, &
830 : pt_el=gt_act_elem, &
831 4613 : prob_opt=tmc_env%params%esimate_acc_prob)
832 :
833 : !writes only different configurations with repetition file if possible
834 : CALL write_result_list_element(result_list=tmc_env%m_env%result_list, &
835 : result_count=tmc_env%m_env%result_count, &
836 : conf_updated=gt_act_elem%mv_conf, accepted=result_acc, &
837 4613 : tmc_params=tmc_env%params)
838 4613 : IF (gt_act_elem%swaped) &
839 : CALL write_result_list_element(result_list=tmc_env%m_env%result_list, &
840 : result_count=tmc_env%m_env%result_count, &
841 : conf_updated=gt_act_elem%mv_conf + 1, accepted=result_acc, &
842 168 : tmc_params=tmc_env%params)
843 :
844 : ! save for analysis
845 283014 : IF (tmc_env%tmc_comp_set%para_env_m_ana%num_pe > 1 .AND. result_acc) THEN
846 : CALL add_to_list(elem=tmc_env%m_env%result_list(gt_act_elem%mv_conf)%elem, &
847 : list=tmc_env%m_env%analysis_list, &
848 : temp_ind=gt_act_elem%mv_conf, &
849 0 : nr=tmc_env%m_env%result_count(gt_act_elem%mv_conf))
850 0 : IF (gt_act_elem%swaped) THEN
851 : CALL add_to_list(elem=tmc_env%m_env%result_list(gt_act_elem%mv_conf + 1)%elem, &
852 : list=tmc_env%m_env%analysis_list, &
853 : temp_ind=gt_act_elem%mv_conf + 1, &
854 0 : nr=tmc_env%m_env%result_count(gt_act_elem%mv_conf + 1))
855 : END IF
856 : END IF
857 : END DO search_calculated_element_loop
858 :
859 : ! end the timing
860 278401 : CALL timestop(handle)
861 :
862 278401 : END SUBROUTINE tree_update
863 :
864 : !============================================================================
865 : ! tree node acceptance probability
866 : !============================================================================
867 : ! **************************************************************************************************
868 : !> \brief checks if element is ready for accaptance probability update
869 : !> checks status and the amount of already received intermediate energies
870 : !> \param elem sub tree element to update
871 : !> \return return value
872 : !> \author Mandes 12.2012
873 : ! **************************************************************************************************
874 0 : FUNCTION ready_for_update_acc_prob(elem) RESULT(ready)
875 : TYPE(tree_type), POINTER :: elem
876 : LOGICAL :: ready
877 :
878 0 : CPASSERT(ASSOCIATED(elem))
879 0 : ready = .FALSE.
880 : IF ((elem%scf_energies_count >= 4) &
881 : .AND. (elem%stat /= status_deleted) .AND. (elem%stat /= status_deleted_result) &
882 0 : .AND. (elem%stat /= status_canceled_ener)) &
883 0 : ready = .TRUE.
884 0 : END FUNCTION ready_for_update_acc_prob
885 :
886 : ! **************************************************************************************************
887 : !> \brief updates the subtree acceptance probability
888 : !> the swap probabilities are handled within the certain checks of the
889 : !> global tree elements (pt references)
890 : !> \param tree_elem sub tree element to update
891 : !> \param tmc_env TMC environment
892 : !> \author Mandes 12.2012
893 : ! **************************************************************************************************
894 0 : SUBROUTINE check_elements_for_acc_prob_update(tree_elem, tmc_env)
895 : TYPE(tree_type), POINTER :: tree_elem
896 : TYPE(tmc_env_type), POINTER :: tmc_env
897 :
898 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_elements_for_acc_prob_update'
899 :
900 : INTEGER :: handle
901 : TYPE(tree_type), POINTER :: act_elem, parent_elem
902 :
903 0 : NULLIFY (parent_elem, act_elem)
904 :
905 0 : CPASSERT(ASSOCIATED(tree_elem))
906 0 : CPASSERT(ASSOCIATED(tmc_env))
907 :
908 : ! start the timing
909 0 : CALL timeset(routineN, handle)
910 :
911 0 : act_elem => tree_elem
912 : !-- nothing to do if option is disabled or element not ready
913 0 : IF (tmc_env%params%esimate_acc_prob .AND. &
914 : ready_for_update_acc_prob(act_elem)) THEN
915 : ! ------------------------------------------------------
916 : ! check this element
917 : ! for usual moves and swapping
918 : ! get parent subtree elment for case of not swapped configurations
919 0 : parent_elem => search_parent_element(act_elem)
920 : ! update the prob of accaptance
921 : CALL update_prob_gt_node_list(reference_list=act_elem%gt_nodes_references, &
922 : act_elem=act_elem, parent_elem=parent_elem, act_parent=.TRUE., &
923 0 : tmc_env=tmc_env)
924 :
925 : !------------------------------------------------------
926 : ! check the childs of this element
927 0 : parent_elem => tree_elem
928 :
929 : ! check ACC child (parent element is the entered/updated element)
930 0 : IF (ASSOCIATED(tree_elem%acc)) THEN
931 0 : act_elem => tree_elem%acc
932 : ! update the prob of accaptance
933 : CALL update_prob_gt_node_list(reference_list=act_elem%gt_nodes_references, &
934 : act_elem=act_elem, parent_elem=parent_elem, act_parent=.FALSE., &
935 0 : tmc_env=tmc_env)
936 : END IF
937 :
938 : ! check all NACC childs of next accepted one
939 0 : nacc_loop: DO
940 0 : IF (.NOT. ASSOCIATED(act_elem%nacc)) EXIT nacc_loop
941 0 : act_elem => act_elem%nacc
942 : CALL update_prob_gt_node_list(reference_list=act_elem%gt_nodes_references, &
943 : act_elem=act_elem, parent_elem=parent_elem, act_parent=.FALSE., &
944 0 : tmc_env=tmc_env)
945 : END DO nacc_loop
946 : END IF
947 : ! end the timing
948 0 : CALL timestop(handle)
949 0 : END SUBROUTINE check_elements_for_acc_prob_update
950 :
951 : ! **************************************************************************************************
952 : !> \brief search all pt nodes in reference list and update the estimated
953 : !> acceptance probabilities
954 : !> \param reference_list list of global tree element referring to the updated
955 : !> subtree element
956 : !> \param act_elem actual sub tree element updated or child of the updated one
957 : !> \param parent_elem parent or the actual updated subtree element
958 : !> \param act_parent flag if updated element is the actual (TRUE) or
959 : !> the parent (FALSE) element
960 : !> \param tmc_env TMC environment
961 : !> \author Mandes 12.2012
962 : ! **************************************************************************************************
963 0 : SUBROUTINE update_prob_gt_node_list(reference_list, act_elem, parent_elem, act_parent, tmc_env)
964 : TYPE(gt_elem_list_type), POINTER :: reference_list
965 : TYPE(tree_type), POINTER :: act_elem, parent_elem
966 : LOGICAL :: act_parent
967 : TYPE(tmc_env_type), POINTER :: tmc_env
968 :
969 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_prob_gt_node_list'
970 :
971 : INTEGER :: handle, integration_steps
972 : REAL(KIND=dp) :: kB, tmp_prob
973 : TYPE(gt_elem_list_type), POINTER :: tmp_pt_ptr
974 : TYPE(tree_type), POINTER :: tmp_elem, tmp_parent_elem
975 :
976 0 : kB = boltzmann/joule
977 :
978 0 : NULLIFY (tmp_pt_ptr, tmp_elem, tmp_parent_elem)
979 :
980 0 : IF (.NOT. ASSOCIATED(reference_list)) RETURN ! could be canceled already
981 :
982 0 : CPASSERT(ASSOCIATED(reference_list%gt_elem))
983 0 : CPASSERT(ASSOCIATED(act_elem))
984 0 : CPASSERT(ASSOCIATED(parent_elem))
985 :
986 : ! start the timing
987 0 : CALL timeset(routineN, handle)
988 :
989 : ! check if element is ready
990 0 : IF (ready_for_update_acc_prob(act_elem)) THEN
991 : ! set the number of integration steps used for 3 point approximation
992 : ! of the exact energy, using the sub step energies
993 : ! still fixed value
994 0 : integration_steps = 100
995 :
996 0 : tmp_pt_ptr => reference_list
997 : ! do for all references using the updated subtree node
998 0 : reference_loop: DO WHILE (ASSOCIATED(tmp_pt_ptr))
999 0 : tmp_prob = -1.0_dp
1000 0 : IF (tmp_pt_ptr%gt_elem%swaped) THEN
1001 0 : NULLIFY (tmp_elem, tmp_parent_elem)
1002 : ! in case of swap use the other swaped configuration as related one
1003 : CALL get_subtree_elements_to_check(gt_act_elem=tmp_pt_ptr%gt_elem, &
1004 0 : elem1=tmp_elem, elem2=tmp_parent_elem)
1005 : ! NOT if parent is the updated one, and check for correct elements ready
1006 0 : IF (act_parent .AND. ASSOCIATED(act_elem, tmp_elem) .AND. &
1007 : ready_for_update_acc_prob(elem=tmp_parent_elem)) THEN
1008 : ! using ln(rnd)/-(beta_i-beta_j) < U_j-U_i)
1009 : tmp_prob = compute_estimated_prob(elem_old=tmp_parent_elem, &
1010 : elem_new=act_elem, &
1011 : E_classical_diff=0.0_dp, &
1012 : rnd_nr=tmp_pt_ptr%gt_elem%rnd_nr, &
1013 : beta=1.0_dp/(kB*(tmc_env%params%Temp(tmp_pt_ptr%gt_elem%mv_conf) - &
1014 : tmc_env%params%Temp(tmp_pt_ptr%gt_elem%mv_conf + 1))), &
1015 0 : tmc_params=tmc_env%params)
1016 : ELSE
1017 0 : tmp_pt_ptr => tmp_pt_ptr%next
1018 0 : CYCLE reference_loop
1019 : END IF
1020 : ELSE
1021 : ! if no swap, use the parent configuration as related one
1022 0 : IF (ready_for_update_acc_prob(parent_elem)) THEN
1023 : tmp_prob = compute_estimated_prob(elem_old=parent_elem, &
1024 : elem_new=act_elem, &
1025 : E_classical_diff=act_elem%e_pot_approx - parent_elem%e_pot_approx, &
1026 : rnd_nr=tmp_pt_ptr%gt_elem%rnd_nr, &
1027 : beta=1.0_dp/(kB*tmc_env%params%Temp(tmp_pt_ptr%gt_elem%mv_conf)), &
1028 0 : tmc_params=tmc_env%params)
1029 : END IF
1030 : END IF
1031 : !successful probability update
1032 0 : IF (tmp_prob >= 0.0_dp) THEN
1033 0 : tmp_pt_ptr%gt_elem%prob_acc = tmp_prob
1034 : !-- speculative canceling for the related direction
1035 0 : IF (tmc_env%params%SPECULATIVE_CANCELING) &
1036 : CALL search_canceling_elements(pt_elem_in=tmp_pt_ptr%gt_elem, &
1037 0 : prob=tmp_pt_ptr%gt_elem%prob_acc, tmc_env=tmc_env)
1038 : END IF
1039 :
1040 : ! get next related global tree pointer
1041 0 : tmp_pt_ptr => tmp_pt_ptr%next
1042 : END DO reference_loop ! global tree references of actual subtree element
1043 : END IF
1044 : ! end the timing
1045 0 : CALL timestop(handle)
1046 : END SUBROUTINE update_prob_gt_node_list
1047 :
1048 : END MODULE tmc_tree_acceptance
|
