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