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