Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief acceptance ratio handling of the different Monte Carlo Moves types
10 : !> For each move type and each temperature average acceptance is
11 : !> determined.
12 : !> For each move is a weight (mv_weight) defined, which defines the
13 : !> probability to perform the move.
14 : !> We distinguish between moves performed on the exact potential
15 : !> (move on the master, energy on the energy worker) and
16 : !> NMC moves, which are performed on the worker using the approximate
17 : !> potential. The energies are calculated as usual on the energy worker
18 : !> with the exact potential.
19 : !> The move probabilities to perform a NMC is stored in the NMC move.
20 : !> The probilities of the single move types (performed with the
21 : !> approximate potential) are only compared within the NMC move
22 : !> \par History
23 : !> 11.2012 created [Mandes Schoenherr]
24 : !> \author Mandes
25 : ! **************************************************************************************************
26 :
27 : MODULE tmc_move_handle
28 : USE cp_log_handling, ONLY: cp_to_string
29 : USE input_section_types, ONLY: section_vals_get,&
30 : section_vals_get_subs_vals,&
31 : section_vals_type,&
32 : section_vals_val_get
33 : USE kinds, ONLY: default_string_length,&
34 : dp
35 : USE mathconstants, ONLY: pi
36 : USE physcon, ONLY: au2a => angstrom
37 : USE string_utilities, ONLY: uppercase
38 : USE tmc_move_types, ONLY: &
39 : move_types_create, move_types_release, mv_type_MD, mv_type_NMC_moves, mv_type_atom_swap, &
40 : mv_type_atom_trans, mv_type_gausian_adapt, mv_type_mol_rot, mv_type_mol_trans, &
41 : mv_type_proton_reorder, mv_type_swap_conf, mv_type_volume_move, tmc_move_type
42 : USE tmc_stati, ONLY: task_type_MC,&
43 : task_type_gaussian_adaptation,&
44 : task_type_ideal_gas
45 : USE tmc_tree_types, ONLY: global_tree_type,&
46 : status_accepted_result,&
47 : status_rejected_result,&
48 : tree_type
49 : USE tmc_types, ONLY: tmc_param_type
50 : #include "../base/base_uses.f90"
51 :
52 : IMPLICIT NONE
53 :
54 : PRIVATE
55 :
56 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_move_handle'
57 :
58 : PUBLIC :: finalize_mv_types, print_move_types, read_init_move_types
59 : PUBLIC :: check_moves
60 : PUBLIC :: select_random_move_type
61 : PUBLIC :: prob_update, add_mv_prob
62 : PUBLIC :: clear_move_probs
63 :
64 : CONTAINS
65 :
66 : ! **************************************************************************************************
67 : !> \brief initialization of the different moves, with sizes and probabilities
68 : !> \param tmc_params ...
69 : !> \param tmc_section ...
70 : !> \author Mandes 10.2013
71 : ! **************************************************************************************************
72 84 : SUBROUTINE read_init_move_types(tmc_params, tmc_section)
73 : TYPE(tmc_param_type), POINTER :: tmc_params
74 : TYPE(section_vals_type), POINTER :: tmc_section
75 :
76 : CHARACTER(LEN=default_string_length) :: inp_kind_name
77 : INTEGER :: i, i_rep, i_tmp, ind, n_items, &
78 : n_NMC_items, n_rep_val, nmc_steps
79 : LOGICAL :: explicit, flag
80 : REAL(KIND=dp) :: delta_x, init_acc_prob, mv_prob, &
81 : mv_prob_sum, nmc_init_acc_prob, &
82 : nmc_prob, nmc_prob_sum, prob_ex
83 : TYPE(section_vals_type), POINTER :: move_type_section, nmc_section
84 : TYPE(tmc_move_type), POINTER :: move_types
85 :
86 28 : NULLIFY (move_types, move_type_section, nmc_section)
87 :
88 28 : n_items = 0
89 28 : n_NMC_items = 0
90 : delta_x = 0.0_dp
91 28 : nmc_prob = 0.0_dp
92 : mv_prob = 0.0_dp
93 : nmc_prob = 0.0_dp
94 28 : mv_prob_sum = 0.0_dp
95 28 : nmc_prob_sum = 0.0_dp
96 28 : prob_ex = 0.0_dp
97 : init_acc_prob = 0.0_dp
98 :
99 : ! the move types on exact potential
100 28 : move_type_section => section_vals_get_subs_vals(tmc_section, "MOVE_TYPE")
101 28 : CALL section_vals_get(move_type_section, explicit=explicit)
102 28 : IF (explicit) THEN
103 22 : CALL section_vals_get(move_type_section, n_repetition=n_items)
104 22 : mv_prob_sum = 0.0_dp
105 82 : DO i_rep = 1, n_items
106 : CALL section_vals_val_get(move_type_section, "PROB", i_rep_section=i_rep, &
107 60 : r_val=mv_prob)
108 82 : mv_prob_sum = mv_prob_sum + mv_prob
109 : END DO
110 : END IF
111 :
112 : ! get the NMC prameters
113 28 : nmc_section => section_vals_get_subs_vals(tmc_section, "NMC_MOVES")
114 28 : CALL section_vals_get(nmc_section, explicit=explicit)
115 28 : IF (explicit) THEN
116 : ! check the approx potential file, already read
117 10 : IF (tmc_params%NMC_inp_file == "") THEN
118 0 : CPABORT("Please specify a valid approximate potential.")
119 : END IF
120 :
121 : CALL section_vals_val_get(nmc_section, "NR_NMC_STEPS", &
122 10 : i_val=nmc_steps)
123 10 : IF (nmc_steps <= 0) THEN
124 0 : CPABORT("Please specify a valid amount of NMC steps (NR_NMC_STEPS {INTEGER}).")
125 : END IF
126 :
127 10 : CALL section_vals_val_get(nmc_section, "PROB", r_val=nmc_prob)
128 :
129 : CALL section_vals_val_get(move_type_section, "INIT_ACC_PROB", &
130 10 : r_val=nmc_init_acc_prob)
131 10 : IF (nmc_init_acc_prob <= 0.0_dp) THEN
132 : CALL cp_abort(__LOCATION__, &
133 : "Please select a valid initial acceptance probability (>0.0) "// &
134 0 : "for INIT_ACC_PROB")
135 : END IF
136 :
137 10 : move_type_section => section_vals_get_subs_vals(nmc_section, "MOVE_TYPE")
138 10 : CALL section_vals_get(move_type_section, n_repetition=n_NMC_items)
139 :
140 : ! get the NMC move probability sum
141 10 : nmc_prob_sum = 0.0_dp
142 36 : DO i_rep = 1, n_NMC_items
143 : CALL section_vals_val_get(move_type_section, "PROB", i_rep_section=i_rep, &
144 26 : r_val=mv_prob)
145 36 : nmc_prob_sum = nmc_prob_sum + mv_prob
146 : END DO
147 : END IF
148 :
149 : ! get the total weight/amount of move probabilities
150 28 : mv_prob_sum = mv_prob_sum + nmc_prob
151 :
152 28 : IF (n_items + n_NMC_items > 0) THEN
153 : ! initilaize the move array with related sizes, probs, etc.
154 28 : CALL move_types_create(tmc_params%move_types, tmc_params%nr_temp)
155 :
156 28 : IF (mv_prob_sum <= 0.0) THEN
157 : CALL cp_abort(__LOCATION__, &
158 : "The probabilities to perform the moves are "// &
159 0 : "in total less equal 0")
160 : END IF
161 :
162 : ! get the sizes, probs, etc. for each move type and convert units
163 114 : DO i_tmp = 1, n_items + n_NMC_items
164 : ! select the correct section
165 86 : IF (i_tmp > n_items) THEN
166 26 : i_rep = i_tmp - n_items
167 26 : IF (i_rep == 1) THEN
168 : ! set the NMC stuff (approx potential)
169 : tmc_params%move_types%mv_weight(mv_type_NMC_moves) = &
170 10 : nmc_prob/REAL(mv_prob_sum, KIND=dp)
171 24 : tmc_params%move_types%mv_size(mv_type_NMC_moves, :) = nmc_steps
172 24 : tmc_params%move_types%acc_prob(mv_type_NMC_moves, :) = nmc_init_acc_prob
173 :
174 10 : move_type_section => section_vals_get_subs_vals(tmc_section, "NMC_MOVES%MOVE_TYPE")
175 10 : mv_prob_sum = nmc_prob_sum
176 : ! allocate the NMC move types
177 10 : CALL move_types_create(tmc_params%nmc_move_types, tmc_params%nr_temp)
178 10 : move_types => tmc_params%nmc_move_types
179 : END IF
180 : ELSE
181 : ! the moves on exact potential
182 60 : move_type_section => section_vals_get_subs_vals(tmc_section, "MOVE_TYPE")
183 60 : i_rep = i_tmp
184 60 : move_types => tmc_params%move_types
185 : END IF
186 :
187 : CALL section_vals_val_get(move_type_section, "_SECTION_PARAMETERS_", &
188 86 : c_val=inp_kind_name, i_rep_section=i_rep)
189 86 : CALL uppercase(inp_kind_name)
190 : CALL section_vals_val_get(move_type_section, "SIZE", i_rep_section=i_rep, &
191 86 : r_val=delta_x)
192 : ! move sizes are checked afterwards, because not all moves require a valid move size
193 : CALL section_vals_val_get(move_type_section, "PROB", i_rep_section=i_rep, &
194 86 : r_val=mv_prob)
195 86 : IF (mv_prob < 0.0_dp) THEN
196 : CALL cp_abort(__LOCATION__, &
197 : "Please select a valid move probability (>0.0) "// &
198 0 : "for the move type "//inp_kind_name)
199 : END IF
200 : CALL section_vals_val_get(move_type_section, "INIT_ACC_PROB", i_rep_section=i_rep, &
201 86 : r_val=init_acc_prob)
202 86 : IF (init_acc_prob < 0.0_dp) THEN
203 : CALL cp_abort(__LOCATION__, &
204 : "Please select a valid initial acceptance probability (>0.0) "// &
205 0 : "for the move type "//inp_kind_name)
206 : END IF
207 : ! set the related index and perform unit conversion of move sizes
208 44 : SELECT CASE (inp_kind_name)
209 : ! atom / molecule translation
210 : CASE ("ATOM_TRANS", "MOL_TRANS")
211 : SELECT CASE (inp_kind_name)
212 : CASE ("ATOM_TRANS")
213 16 : ind = mv_type_atom_trans
214 : CASE ("MOL_TRANS")
215 16 : ind = mv_type_mol_trans
216 : CASE DEFAULT
217 44 : CPABORT("move type is not defined in the translation types")
218 : END SELECT
219 : ! convert units
220 104 : SELECT CASE (tmc_params%task_type)
221 : CASE (task_type_MC, task_type_ideal_gas)
222 44 : delta_x = delta_x/au2a
223 : CASE (task_type_gaussian_adaptation)
224 : !nothing to do (no unit conversion)
225 : CASE DEFAULT
226 44 : CPABORT("move type atom / mol trans is not defined for this TMC run type")
227 : END SELECT
228 : ! molecule rotation
229 : CASE ("MOL_ROT")
230 16 : ind = mv_type_mol_rot
231 : ! convert units
232 34 : SELECT CASE (tmc_params%task_type)
233 : CASE (task_type_MC, task_type_ideal_gas)
234 16 : delta_x = delta_x*PI/180.0_dp
235 : CASE DEFAULT
236 16 : CPABORT("move type MOL_ROT is not defined for this TMC run type")
237 : END SELECT
238 : ! proton reordering
239 : CASE ("PROT_REORDER")
240 2 : ind = mv_type_proton_reorder
241 : ! the move size is not necessary
242 2 : delta_x = 0.0_dp
243 : ! Hybrid MC (MD)
244 : CASE ("HYBRID_MC")
245 0 : ind = mv_type_MD
246 0 : delta_x = delta_x*Pi/180.0_dp !input in degree, calculating in rad
247 0 : tmc_params%print_forces = .TRUE.
248 : ! parallel tempering swap move
249 : CASE ("PT_SWAP")
250 8 : ind = mv_type_swap_conf
251 : ! the move size is not necessary
252 8 : delta_x = 0.0_dp
253 8 : IF (tmc_params%nr_temp <= 1) THEN
254 : ! no configurational swapping if only one temperature
255 0 : mv_prob = 0.0_dp
256 : CALL cp_warn(__LOCATION__, &
257 : "Configurational swap disabled, because "// &
258 0 : "Parallel Tempering requires more than one temperature.")
259 : END IF
260 : ! volume moves
261 : CASE ("VOL_MOVE")
262 12 : ind = mv_type_volume_move
263 : ! check the selected pressure
264 12 : IF (tmc_params%pressure >= 0.0_dp) THEN
265 12 : delta_x = delta_x/au2a
266 12 : tmc_params%print_cell = .TRUE. ! print the cell sizes by default
267 : ELSE
268 : CALL cp_warn(__LOCATION__, &
269 : "no valid pressure defined, but volume move defined. "// &
270 0 : "Consequently, the volume move is disabled.")
271 0 : mv_prob = 0.0_dp
272 : END IF
273 : ! parallel tempering swap move
274 : CASE ("ATOM_SWAP")
275 4 : ind = mv_type_atom_swap
276 : ! the move size is not necessary
277 4 : delta_x = 0.0_dp
278 : ! select the types of atoms swapped
279 : CALL section_vals_val_get(move_type_section, "ATOMS", i_rep_section=i_rep, &
280 4 : n_rep_val=n_rep_val)
281 4 : IF (n_rep_val > 0) THEN
282 10 : ALLOCATE (move_types%atom_lists(n_rep_val))
283 6 : DO i = 1, n_rep_val
284 : CALL section_vals_val_get(move_type_section, "ATOMS", &
285 : i_rep_section=i_rep, i_rep_val=i, &
286 4 : c_vals=move_types%atom_lists(i)%atoms)
287 6 : IF (SIZE(move_types%atom_lists(i)%atoms) <= 1) THEN
288 0 : CPABORT("ATOM_SWAP requires minimum two atom kinds selected. ")
289 : END IF
290 : END DO
291 : END IF
292 : ! gaussian adaptation
293 : CASE ("GAUSS_ADAPT")
294 0 : ind = mv_type_gausian_adapt
295 0 : init_acc_prob = 0.5_dp
296 : CASE DEFAULT
297 86 : CPABORT("A unknown move type is selected: "//inp_kind_name)
298 : END SELECT
299 : ! check for valid move sizes
300 86 : IF (delta_x < 0.0_dp) THEN
301 : CALL cp_abort(__LOCATION__, &
302 : "Please select a valid move size (>0.0) "// &
303 0 : "for the move type "//inp_kind_name)
304 : END IF
305 : ! check if not already set
306 86 : IF (move_types%mv_weight(ind) > 0.0) THEN
307 0 : CPABORT("TMC: Each move type can be set only once. ")
308 : END IF
309 :
310 : ! set the move size
311 264 : move_types%mv_size(ind, :) = delta_x
312 : ! set the probability to perform move
313 86 : move_types%mv_weight(ind) = mv_prob/mv_prob_sum
314 : ! set the initial acceptance probability
315 464 : move_types%acc_prob(ind, :) = init_acc_prob
316 : END DO
317 : ELSE
318 0 : CPABORT("No move type selected, please select at least one.")
319 : END IF
320 308 : mv_prob_sum = SUM(tmc_params%move_types%mv_weight(:))
321 28 : flag = .TRUE.
322 28 : CPASSERT(ABS(mv_prob_sum - 1.0_dp) < 0.01_dp)
323 28 : IF (ASSOCIATED(tmc_params%nmc_move_types)) THEN
324 110 : mv_prob_sum = SUM(tmc_params%nmc_move_types%mv_weight(:))
325 10 : CPASSERT(ABS(mv_prob_sum - 1.0_dp) < 10*EPSILON(1.0_dp))
326 : END IF
327 28 : END SUBROUTINE read_init_move_types
328 :
329 : ! **************************************************************************************************
330 : !> \brief checks if the moves are possible
331 : !> \param tmc_params ...
332 : !> \param move_types ...
333 : !> \param mol_array ...
334 : !> \author Mandes 10.2013
335 : ! **************************************************************************************************
336 19 : SUBROUTINE check_moves(tmc_params, move_types, mol_array)
337 : TYPE(tmc_param_type), POINTER :: tmc_params
338 : TYPE(tmc_move_type), POINTER :: move_types
339 : INTEGER, DIMENSION(:), POINTER :: mol_array
340 :
341 : INTEGER :: atom_j, list_i, ref_k
342 : LOGICAL :: found
343 :
344 19 : CPASSERT(ASSOCIATED(tmc_params))
345 19 : CPASSERT(ASSOCIATED(move_types))
346 :
347 : ! molecule moves need molecule info
348 19 : IF (move_types%mv_weight(mv_type_mol_trans) > 0.0_dp .OR. &
349 : move_types%mv_weight(mv_type_mol_rot) > 0.0_dp) THEN
350 : ! if there is no molecule information available,
351 : ! molecules moves can not be performed
352 8 : IF (mol_array(SIZE(mol_array)) == SIZE(mol_array)) THEN
353 : CALL cp_abort(__LOCATION__, &
354 : "molecule move: there is no molecule "// &
355 : "information available. Please specify molecules when "// &
356 0 : "using molecule moves.")
357 : END IF
358 : END IF
359 :
360 : ! for the atom swap move
361 19 : IF (move_types%mv_weight(mv_type_atom_swap) > 0.0_dp) THEN
362 : ! check if the selected atom swaps are possible
363 2 : IF (ASSOCIATED(move_types%atom_lists)) THEN
364 3 : DO list_i = 1, SIZE(move_types%atom_lists(:))
365 7 : DO atom_j = 1, SIZE(move_types%atom_lists(list_i)%atoms(:))
366 : ! check if atoms exists
367 4 : found = .FALSE.
368 11 : ref_loop: DO ref_k = 1, SIZE(tmc_params%atoms(:))
369 11 : IF (move_types%atom_lists(list_i)%atoms(atom_j) == &
370 0 : tmc_params%atoms(ref_k)%name) THEN
371 : found = .TRUE.
372 : EXIT ref_loop
373 : END IF
374 : END DO ref_loop
375 4 : IF (.NOT. found) THEN
376 : CALL cp_abort(__LOCATION__, &
377 : "ATOM_SWAP: The selected atom type ("// &
378 : TRIM(move_types%atom_lists(list_i)%atoms(atom_j))// &
379 0 : ") is not contained in the system. ")
380 : END IF
381 : ! check if not be swapped with the same atom type
382 6 : IF (ANY(move_types%atom_lists(list_i)%atoms(atom_j) == &
383 2 : move_types%atom_lists(list_i)%atoms(atom_j + 1:))) THEN
384 : CALL cp_abort(__LOCATION__, &
385 : "ATOM_SWAP can not swap two atoms of same kind ("// &
386 : TRIM(move_types%atom_lists(list_i)%atoms(atom_j))// &
387 0 : ")")
388 : END IF
389 : END DO
390 : END DO
391 : ELSE
392 : ! check if there exisit different atoms
393 1 : found = .FALSE.
394 1 : IF (SIZE(tmc_params%atoms(:)) > 1) THEN
395 1 : ref_lop: DO ref_k = 2, SIZE(tmc_params%atoms(:))
396 1 : IF (tmc_params%atoms(1)%name /= tmc_params%atoms(ref_k)%name) THEN
397 : found = .TRUE.
398 : EXIT ref_lop
399 : END IF
400 : END DO ref_lop
401 : END IF
402 1 : IF (.NOT. found) THEN
403 : CALL cp_abort(__LOCATION__, &
404 : "The system contains only a single atom type,"// &
405 0 : " atom_swap is not possible.")
406 : END IF
407 : END IF
408 : END IF
409 19 : END SUBROUTINE check_moves
410 :
411 : ! **************************************************************************************************
412 : !> \brief deallocating the module variables
413 : !> \param tmc_params ...
414 : !> \author Mandes 11.2012
415 : !> \note deallocating the module variables
416 : ! **************************************************************************************************
417 28 : SUBROUTINE finalize_mv_types(tmc_params)
418 : TYPE(tmc_param_type), POINTER :: tmc_params
419 :
420 28 : CPASSERT(ASSOCIATED(tmc_params))
421 28 : CALL move_types_release(tmc_params%move_types)
422 28 : IF (ASSOCIATED(tmc_params%nmc_move_types)) THEN
423 10 : CALL move_types_release(tmc_params%nmc_move_types)
424 : END IF
425 28 : END SUBROUTINE finalize_mv_types
426 :
427 : ! **************************************************************************************************
428 : !> \brief routine pronts out the probabilities and sized for each type and
429 : !> temperature the output is divided into two parts the init,
430 : !> which is printed out at the beginning of the programm and
431 : !> .NOT.init which are the probabilites and counter printed out every
432 : !> print cycle
433 : !> \param init ...
434 : !> \param file_io ...
435 : !> \param tmc_params ...
436 : !> \author Mandes 11.2012
437 : ! **************************************************************************************************
438 209 : SUBROUTINE print_move_types(init, file_io, tmc_params)
439 : LOGICAL :: init
440 : INTEGER :: file_io
441 : TYPE(tmc_param_type), POINTER :: tmc_params
442 :
443 : CHARACTER(LEN=10) :: c_t
444 : CHARACTER(LEN=50) :: FMT_c, FMT_i, FMT_r
445 : CHARACTER(LEN=500) :: c_a, c_b, c_c, c_d, c_e, c_tit, c_tmp
446 : INTEGER :: column_size, move, nr_nmc_moves, temper, &
447 : typ
448 : LOGICAL :: subbox_out, type_title
449 : TYPE(tmc_move_type), POINTER :: move_types
450 :
451 209 : NULLIFY (move_types)
452 :
453 209 : c_a = ""; c_b = ""; c_c = ""
454 209 : c_d = ""; c_e = ""; c_tit = ""
455 209 : column_size = 10
456 209 : subbox_out = .FALSE.
457 209 : type_title = .FALSE.
458 209 : CPASSERT(file_io > 0)
459 209 : CPASSERT(ASSOCIATED(tmc_params%move_types))
460 :
461 209 : FLUSH (file_io)
462 :
463 : IF (.NOT. init .AND. &
464 785 : tmc_params%move_types%mv_weight(mv_type_NMC_moves) > 0 .AND. &
465 16 : ANY(tmc_params%sub_box_size > 0.0_dp)) subbox_out = .TRUE.
466 :
467 : ! set the format for each typ to add one column
468 209 : WRITE (FMT_c, '("(A,1X,A", I0, ")")') column_size
469 209 : WRITE (FMT_i, '("(A,1X,I", I0, ")")') column_size
470 209 : WRITE (FMT_r, '("(A,1X,F", I0, ".3)")') column_size
471 : !IF(init) &
472 209 : type_title = .TRUE.
473 :
474 209 : nr_nmc_moves = 0
475 209 : IF (ASSOCIATED(tmc_params%nmc_move_types)) THEN
476 38 : nr_nmc_moves = SIZE(tmc_params%nmc_move_types%mv_weight(1:))
477 : END IF
478 :
479 666 : temp_loop: DO temper = 1, tmc_params%nr_temp
480 463 : c_tit = ""; c_a = ""; c_b = ""; c_c = ""
481 463 : IF (init .AND. temper > 1) EXIT temp_loop
482 457 : WRITE (c_t, "(F10.2)") tmc_params%Temp(temper)
483 5904 : typ_loop: DO move = 0, SIZE(tmc_params%move_types%mv_weight) + nr_nmc_moves
484 : ! the NMC moves
485 5447 : IF (move <= SIZE(tmc_params%move_types%mv_weight)) THEN
486 5027 : typ = move
487 5027 : move_types => tmc_params%move_types
488 : ELSE
489 420 : typ = move - SIZE(tmc_params%move_types%mv_weight)
490 420 : move_types => tmc_params%nmc_move_types
491 : END IF
492 : ! total average
493 5904 : IF (typ == 0) THEN
494 : ! line start
495 457 : IF (type_title) WRITE (c_tit, TRIM(FMT_c)) " type temperature |"
496 457 : IF (init) WRITE (c_b, TRIM(FMT_c)) " I I |"
497 457 : IF (init) WRITE (c_c, TRIM(FMT_c)) " V V |"
498 457 : IF (.NOT. init) WRITE (c_a, TRIM(FMT_c)) "probs T="//c_t//" |"
499 457 : IF (.NOT. init) WRITE (c_b, TRIM(FMT_c)) "counts T="//c_t//" |"
500 457 : IF (.NOT. init) WRITE (c_c, TRIM(FMT_c)) "nr_acc T="//c_t//" |"
501 457 : IF (subbox_out) THEN
502 16 : WRITE (c_d, TRIM(FMT_c)) "sb_acc T="//c_t//" |"
503 16 : WRITE (c_e, TRIM(FMT_c)) "sb_cou T="//c_t//" |"
504 : END IF
505 : ! overall column
506 : IF (type_title) THEN
507 457 : c_tmp = TRIM(c_tit)
508 457 : WRITE (c_tit, TRIM(FMT_c)) TRIM(c_tmp), " trajec"
509 : END IF
510 457 : IF (init) THEN
511 14 : c_tmp = TRIM(c_b)
512 14 : WRITE (c_b, TRIM(FMT_c)) TRIM(c_tmp), " weight->"
513 : END IF
514 457 : IF (init) THEN
515 14 : c_tmp = TRIM(c_c)
516 14 : WRITE (c_c, TRIM(FMT_c)) TRIM(c_tmp), " size ->"
517 : END IF
518 457 : IF (.NOT. init) THEN
519 443 : c_tmp = TRIM(c_a)
520 443 : WRITE (c_a, TRIM(FMT_r)) TRIM(c_tmp), &
521 886 : move_types%acc_prob(typ, temper)
522 : END IF
523 457 : IF (.NOT. init) THEN
524 443 : c_tmp = TRIM(c_b)
525 443 : WRITE (c_b, TRIM(FMT_i)) TRIM(c_tmp), &
526 886 : move_types%mv_count(typ, temper)
527 : END IF
528 457 : IF (.NOT. init) THEN
529 443 : c_tmp = TRIM(c_c)
530 443 : WRITE (c_c, TRIM(FMT_i)) TRIM(c_tmp), &
531 886 : move_types%acc_count(typ, temper)
532 : END IF
533 457 : IF (subbox_out) THEN
534 16 : c_tmp = TRIM(c_d)
535 16 : WRITE (c_d, TRIM(FMT_c)) TRIM(c_tmp), "."
536 16 : c_tmp = TRIM(c_e)
537 16 : WRITE (c_e, TRIM(FMT_c)) TRIM(c_tmp), "."
538 : END IF
539 : ELSE
540 : ! certain move types
541 4990 : IF (move_types%mv_weight(typ) > 0.0_dp) THEN
542 : ! INIT: the weights in the initialisation output
543 1717 : IF (init) THEN
544 48 : c_tmp = TRIM(c_b)
545 48 : WRITE (c_b, TRIM(FMT_r)) TRIM(c_tmp), move_types%mv_weight(typ)
546 : END IF
547 : ! acc probabilities
548 1717 : IF (typ == mv_type_swap_conf .AND. &
549 : temper == tmc_params%nr_temp) THEN
550 116 : IF (.NOT. init) THEN
551 116 : c_tmp = TRIM(c_a)
552 116 : WRITE (c_a, TRIM(FMT_c)) TRIM(c_tmp), "---"
553 : END IF
554 : ELSE
555 1601 : IF (.NOT. init) THEN
556 1553 : c_tmp = TRIM(c_a)
557 1553 : WRITE (c_a, TRIM(FMT_r)) TRIM(c_tmp), move_types%acc_prob(typ, temper)
558 : END IF
559 : END IF
560 1717 : IF (.NOT. init) THEN
561 1669 : c_tmp = TRIM(c_b)
562 1669 : WRITE (c_b, TRIM(FMT_i)) TRIM(c_tmp), move_types%mv_count(typ, temper)
563 : END IF
564 1717 : IF (.NOT. init) THEN
565 1669 : c_tmp = TRIM(c_c)
566 1669 : WRITE (c_c, TRIM(FMT_i)) TRIM(c_tmp), move_types%acc_count(typ, temper)
567 : END IF
568 : ! sub box
569 1717 : IF (subbox_out) THEN
570 64 : IF (move > SIZE(tmc_params%move_types%mv_weight)) THEN
571 48 : c_tmp = TRIM(c_d)
572 48 : WRITE (c_d, TRIM(FMT_r)) TRIM(c_tmp), &
573 : move_types%subbox_acc_count(typ, temper)/ &
574 96 : REAL(MAX(1, move_types%subbox_count(typ, temper)), KIND=dp)
575 48 : c_tmp = TRIM(c_e)
576 48 : WRITE (c_e, TRIM(FMT_i)) TRIM(c_tmp), &
577 96 : move_types%subbox_count(typ, temper)
578 : ELSE
579 16 : c_tmp = TRIM(c_d)
580 16 : WRITE (c_d, TRIM(FMT_c)) TRIM(c_tmp), "-"
581 16 : c_tmp = TRIM(c_e)
582 16 : WRITE (c_e, TRIM(FMT_c)) TRIM(c_tmp), "-"
583 : END IF
584 : END IF
585 :
586 : SELECT CASE (typ)
587 : CASE (mv_type_atom_trans)
588 : IF (type_title) THEN
589 457 : c_tmp = TRIM(c_tit)
590 457 : WRITE (c_tit, TRIM(FMT_c)) TRIM(c_tmp), "atom trans."
591 : END IF
592 457 : IF (init) THEN
593 14 : c_tmp = TRIM(c_c)
594 14 : WRITE (c_c, TRIM(FMT_r)) TRIM(c_tmp), &
595 28 : move_types%mv_size(typ, temper)*au2a
596 : END IF
597 : CASE (mv_type_mol_trans)
598 : IF (type_title) THEN
599 403 : c_tmp = TRIM(c_tit)
600 403 : WRITE (c_tit, TRIM(FMT_c)) TRIM(c_tmp), "mol trans"
601 : END IF
602 403 : IF (init) THEN
603 8 : c_tmp = TRIM(c_c)
604 8 : WRITE (c_c, TRIM(FMT_r)) TRIM(c_tmp), &
605 16 : move_types%mv_size(typ, temper)*au2a
606 : END IF
607 : CASE (mv_type_mol_rot)
608 : IF (type_title) THEN
609 403 : c_tmp = TRIM(c_tit)
610 403 : WRITE (c_tit, TRIM(FMT_c)) TRIM(c_tmp), "mol rot"
611 : END IF
612 403 : IF (init) THEN
613 8 : c_tmp = TRIM(c_c)
614 8 : WRITE (c_c, TRIM(FMT_r)) TRIM(c_tmp), &
615 16 : move_types%mv_size(typ, temper)/(PI/180.0_dp)
616 : END IF
617 : CASE (mv_type_MD)
618 0 : CPWARN("md_time_step and nr md_steps not implemented...")
619 : ! IF(type_title) WRITE(c_tit,TRIM(FMT_c)) TRIM(c_tit), "HybridMC"
620 : ! IF(init) WRITE(c_c,TRIM(FMT_c)) TRIM(c_c), "s.above"
621 : ! IF(init) THEN
622 : ! WRITE(file_io,*)" move type: molecular dynamics with file ",NMC_inp_file
623 : ! WRITE(file_io,*)" with time step [fs] ",md_time_step*au2fs
624 : ! WRITE(file_io,*)" with number of steps ",md_steps
625 : ! WRITE(file_io,*)" with velocity changes consists of old vel and ",&
626 : ! sin(move_types%mv_size(typ,1))*100.0_dp,"% random Gaussian with variance to temperature,"
627 : ! END IF
628 : CASE (mv_type_proton_reorder)
629 : IF (type_title) THEN
630 12 : c_tmp = TRIM(c_tit)
631 12 : WRITE (c_tit, TRIM(FMT_c)) TRIM(c_tmp), "H-Reorder"
632 : END IF
633 12 : IF (init) THEN
634 1 : c_tmp = TRIM(c_c)
635 1 : WRITE (c_c, TRIM(FMT_c)) TRIM(c_tmp), "XXX"
636 : END IF
637 : CASE (mv_type_swap_conf)
638 : IF (type_title) THEN
639 352 : c_tmp = TRIM(c_tit)
640 352 : WRITE (c_tit, TRIM(FMT_c)) TRIM(c_tmp), "PT(swap)"
641 : END IF
642 352 : IF (init) THEN
643 4 : c_tmp = TRIM(c_c)
644 4 : WRITE (c_c, TRIM(FMT_c)) TRIM(c_tmp), "XXX" !move_types%mv_size(mv_type_swap_conf,1)
645 : END IF
646 : CASE (mv_type_NMC_moves)
647 : IF (type_title) THEN
648 42 : c_tmp = TRIM(c_tit)
649 42 : WRITE (c_tit, TRIM(FMT_c)) TRIM(c_tmp), "NMC:"
650 : END IF
651 42 : IF (init) THEN
652 5 : c_tmp = TRIM(c_c)
653 5 : WRITE (c_c, TRIM(FMT_i)) TRIM(c_tmp), &
654 10 : INT(move_types%mv_size(typ, temper))
655 : END IF
656 : CASE (mv_type_volume_move)
657 : IF (type_title) THEN
658 42 : c_tmp = TRIM(c_tit)
659 42 : WRITE (c_tit, TRIM(FMT_c)) TRIM(c_tmp), "volume"
660 : END IF
661 42 : IF (init) THEN
662 6 : c_tmp = TRIM(c_c)
663 6 : WRITE (c_c, TRIM(FMT_r)) TRIM(c_tmp), &
664 12 : move_types%mv_size(typ, temper)*au2a
665 : END IF
666 : CASE (mv_type_atom_swap)
667 : IF (type_title) THEN
668 6 : c_tmp = TRIM(c_tit)
669 6 : WRITE (c_tit, TRIM(FMT_c)) TRIM(c_tmp), "atom swap"
670 : END IF
671 6 : IF (init) THEN
672 2 : c_tmp = TRIM(c_c)
673 2 : WRITE (c_c, TRIM(FMT_c)) TRIM(c_tmp), "XXX"
674 : END IF
675 : CASE (mv_type_gausian_adapt)
676 : IF (type_title) THEN
677 0 : c_tmp = TRIM(c_tit)
678 0 : WRITE (c_tit, TRIM(FMT_c)) TRIM(c_tmp), "gauss adap"
679 : END IF
680 0 : IF (init) THEN
681 0 : c_tmp = TRIM(c_c)
682 0 : WRITE (c_c, TRIM(FMT_r)) TRIM(c_tmp), &
683 0 : move_types%mv_size(typ, temper)
684 : END IF
685 : CASE DEFAULT
686 : CALL cp_warn(__LOCATION__, &
687 : "unknown move type "//cp_to_string(typ)//" with weight"// &
688 1717 : cp_to_string(move_types%mv_weight(typ)))
689 : END SELECT
690 : END IF
691 : END IF
692 : END DO typ_loop
693 457 : IF (init) WRITE (UNIT=file_io, FMT="(/,T2,A)") REPEAT("-", 79)
694 457 : IF (type_title .AND. temper <= 1) WRITE (file_io, *) TRIM(c_tit)
695 457 : IF (.NOT. init) WRITE (file_io, *) TRIM(c_a)
696 457 : WRITE (file_io, *) TRIM(c_b)
697 457 : WRITE (file_io, *) TRIM(c_c)
698 457 : IF (subbox_out) WRITE (file_io, *) TRIM(c_d)
699 16 : IF (subbox_out) WRITE (file_io, *) TRIM(c_e)
700 666 : IF (init) WRITE (UNIT=file_io, FMT="(/,T2,A)") REPEAT("-", 79)
701 : END DO temp_loop
702 209 : END SUBROUTINE print_move_types
703 :
704 : ! **************************************************************************************************
705 : !> \brief adaptation of acceptance probability of every kind of change/move
706 : !> and the overall acc prob,
707 : !> using the acceptance and rejectance information
708 : !> \param move_types structure for storing sizes and probabilities of moves
709 : !> \param pt_el global tree element
710 : !> \param elem sub tree element
711 : !> \param acc input if the element is accepted
712 : !> \param subbox logical if move was with respect to the sub box
713 : !> \param prob_opt if the average probability should be adapted
714 : !> \author Mandes 12.2012
715 : ! **************************************************************************************************
716 18128 : SUBROUTINE prob_update(move_types, pt_el, elem, acc, subbox, prob_opt)
717 : TYPE(tmc_move_type), POINTER :: move_types
718 : TYPE(global_tree_type), OPTIONAL, POINTER :: pt_el
719 : TYPE(tree_type), OPTIONAL, POINTER :: elem
720 : LOGICAL, INTENT(IN), OPTIONAL :: acc, subbox
721 : LOGICAL, INTENT(IN) :: prob_opt
722 :
723 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prob_update'
724 :
725 : INTEGER :: change_res, change_sb_type, change_type, &
726 : conf_moved, handle, mv_type
727 :
728 9064 : CPASSERT(ASSOCIATED(move_types))
729 9064 : CPASSERT(.NOT. (PRESENT(pt_el) .AND. PRESENT(subbox)))
730 :
731 : ! start the timing
732 9064 : CALL timeset(routineN, handle)
733 :
734 9064 : mv_type = -1
735 9064 : conf_moved = -1
736 :
737 9064 : change_type = 0
738 9064 : change_res = 0
739 9064 : change_sb_type = 0
740 : ! updating probability of the trajectory
741 9064 : IF (PRESENT(pt_el)) THEN
742 4479 : CPASSERT(ASSOCIATED(pt_el))
743 4479 : conf_moved = pt_el%mv_conf
744 5285 : SELECT CASE (pt_el%stat)
745 : CASE (status_accepted_result)
746 806 : change_res = 1
747 : !-- swaped move is not noted in subtree elements
748 806 : IF (pt_el%swaped) THEN
749 130 : mv_type = mv_type_swap_conf
750 130 : change_type = 1
751 : END IF
752 : CASE (status_rejected_result)
753 3673 : change_res = -1
754 : !-- swaped move is not noted in subtree elements
755 3673 : IF (pt_el%swaped) THEN
756 38 : mv_type = mv_type_swap_conf
757 38 : change_type = -1
758 : END IF
759 : CASE DEFAULT
760 : CALL cp_abort(__LOCATION__, &
761 : "global elem"//cp_to_string(pt_el%nr)// &
762 4479 : "has unknown status"//cp_to_string(pt_el%stat))
763 : END SELECT
764 : END IF
765 :
766 9064 : IF (PRESENT(elem)) THEN
767 4585 : CPASSERT(ASSOCIATED(elem))
768 : !conf_moved = elem%sub_tree_nr
769 4585 : conf_moved = elem%temp_created
770 4585 : mv_type = elem%move_type
771 : ! for NMC prob update the acceptance is needed
772 4585 : CPASSERT(PRESENT(acc))
773 4585 : IF (PRESENT(subbox)) THEN
774 : ! only update subbox acceptance
775 137 : IF (acc) THEN
776 126 : move_types%subbox_acc_count(mv_type, conf_moved) = move_types%subbox_acc_count(mv_type, conf_moved) + 1
777 : END IF
778 137 : move_types%subbox_count(mv_type, conf_moved) = move_types%subbox_count(mv_type, conf_moved) + 1
779 : ! No more to do
780 : change_type = 0
781 : change_res = 0
782 : conf_moved = 0
783 : ! RETURN
784 : ELSE
785 : ! update move type acceptance
786 4448 : IF (acc) THEN
787 : change_type = 1
788 : ELSE
789 : change_type = -1
790 : END IF
791 : END IF
792 : END IF
793 :
794 : !-- INcrease or DEcrease accaptance rate
795 : ! MOVE types
796 4479 : IF (change_type > 0) THEN
797 866 : move_types%acc_count(mv_type, conf_moved) = move_types%acc_count(mv_type, conf_moved) + 1
798 : END IF
799 :
800 : ! RESULTs
801 9064 : IF (change_res > 0) THEN
802 806 : move_types%acc_count(0, conf_moved) = move_types%acc_count(0, conf_moved) + 1
803 : END IF
804 :
805 9064 : IF (conf_moved > 0) move_types%mv_count(0, conf_moved) = move_types%mv_count(0, conf_moved) + ABS(change_res)
806 9064 : IF (mv_type >= 0 .AND. conf_moved > 0) THEN
807 4616 : move_types%mv_count(mv_type, conf_moved) = move_types%mv_count(mv_type, conf_moved) + ABS(change_type)
808 : END IF
809 :
810 9064 : IF (prob_opt) THEN
811 : WHERE (move_types%mv_count > 0) &
812 157528 : move_types%acc_prob(:, :) = move_types%acc_count(:, :)/REAL(move_types%mv_count(:, :), KIND=dp)
813 : END IF
814 : ! end the timing
815 9064 : CALL timestop(handle)
816 9064 : END SUBROUTINE prob_update
817 :
818 : ! **************************************************************************************************
819 : !> \brief add the actual moves to the average probabilities
820 : !> \param move_types structure with move counters and probabilities
821 : !> \param prob_opt ...
822 : !> \param mv_counter move counter for actual performed moves of certain types
823 : !> \param acc_counter counters of acceptance for these moves
824 : !> \param subbox_counter same for sub box moves
825 : !> \param subbox_acc_counter same for sub box moves
826 : !> \author Mandes 12.2012
827 : ! **************************************************************************************************
828 142 : SUBROUTINE add_mv_prob(move_types, prob_opt, mv_counter, acc_counter, &
829 142 : subbox_counter, subbox_acc_counter)
830 : TYPE(tmc_move_type), POINTER :: move_types
831 : LOGICAL :: prob_opt
832 : INTEGER, DIMENSION(:, :), OPTIONAL :: mv_counter, acc_counter, subbox_counter, &
833 : subbox_acc_counter
834 :
835 71 : CPASSERT(ASSOCIATED(move_types))
836 71 : CPASSERT(PRESENT(mv_counter) .OR. PRESENT(subbox_counter))
837 :
838 71 : IF (PRESENT(mv_counter)) THEN
839 57 : CPASSERT(PRESENT(acc_counter))
840 789 : move_types%mv_count(:, :) = move_types%mv_count(:, :) + mv_counter(:, :)
841 789 : move_types%acc_count(:, :) = move_types%acc_count(:, :) + acc_counter(:, :)
842 57 : IF (prob_opt) THEN
843 : WHERE (move_types%mv_count > 0) &
844 789 : move_types%acc_prob(:, :) = move_types%acc_count(:, :)/REAL(move_types%mv_count(:, :), KIND=dp)
845 : END IF
846 : END IF
847 :
848 71 : IF (PRESENT(subbox_counter)) THEN
849 14 : CPASSERT(PRESENT(subbox_acc_counter))
850 168 : move_types%subbox_count(:, :) = move_types%subbox_count(:, :) + subbox_counter(:, :)
851 168 : move_types%subbox_acc_count(:, :) = move_types%subbox_acc_count(:, :) + subbox_acc_counter(:, :)
852 : END IF
853 71 : END SUBROUTINE add_mv_prob
854 :
855 : ! **************************************************************************************************
856 : !> \brief clear the statistics of accepting/rejection moves
857 : !> because worker statistics will be add separately on masters counters
858 : !> \param move_types counters for acceptance/rejection
859 : !> \author Mandes 02.2013
860 : ! **************************************************************************************************
861 57 : SUBROUTINE clear_move_probs(move_types)
862 : TYPE(tmc_move_type), POINTER :: move_types
863 :
864 57 : CPASSERT(ASSOCIATED(move_types))
865 :
866 789 : move_types%acc_prob(:, :) = 0.5_dp
867 789 : move_types%acc_count(:, :) = 0
868 789 : move_types%mv_count(:, :) = 0
869 728 : move_types%subbox_acc_count(:, :) = 0
870 728 : move_types%subbox_count(:, :) = 0
871 57 : END SUBROUTINE clear_move_probs
872 :
873 : ! **************************************************************************************************
874 : !> \brief selects a move type related to the weighings and the entered rnd nr
875 : !> \param move_types structure for storing sizes and probabilities of moves
876 : !> \param rnd random number
877 : !> \return (result) move type
878 : !> \author Mandes 12.2012
879 : !> \note function returns a possible move type without the PT swap moves
880 : !> \note (are selected in global tree, this routine is for sub tree elements)
881 : ! **************************************************************************************************
882 8898 : FUNCTION select_random_move_type(move_types, rnd) RESULT(mv_type)
883 : TYPE(tmc_move_type), POINTER :: move_types
884 : REAL(KIND=dp) :: rnd
885 : INTEGER :: mv_type
886 :
887 : CHARACTER(LEN=*), PARAMETER :: routineN = 'select_random_move_type'
888 :
889 : INTEGER :: handle, i
890 : REAL(KIND=dp) :: rnd_mv, total_moves
891 :
892 4449 : CPASSERT(ASSOCIATED(move_types))
893 4449 : CPASSERT(rnd >= 0.0_dp .AND. rnd < 1.0_dp)
894 :
895 4449 : CALL timeset(routineN, handle)
896 :
897 44490 : total_moves = SUM(move_types%mv_weight(2:))
898 4449 : rnd_mv = total_moves*rnd
899 4449 : mv_type = 0
900 7754 : search_loop: DO i = 2, SIZE(move_types%mv_weight(:))
901 25176 : IF (SUM(move_types%mv_weight(2:i)) >= rnd_mv) THEN
902 : mv_type = i
903 : EXIT search_loop
904 : END IF
905 : END DO search_loop
906 :
907 4449 : CALL timestop(handle)
908 4449 : END FUNCTION select_random_move_type
909 :
910 : END MODULE tmc_move_handle
|
