1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief different move types are applied 8!> \par History 9!> 11.2012 created [Mandes Schoenherr] 10!> \author Mandes 11/2012 11! ************************************************************************************************** 12 13MODULE tmc_moves 14 USE cell_types, ONLY: cell_type,& 15 get_cell,& 16 pbc 17 USE cp_log_handling, ONLY: cp_to_string 18 USE kinds, ONLY: dp 19 USE mathconstants, ONLY: pi 20 USE mathlib, ONLY: dihedral_angle,& 21 rotate_vector 22 USE parallel_rng_types, ONLY: rng_stream_type 23 USE physcon, ONLY: boltzmann,& 24 joule 25 USE tmc_calculations, ONLY: center_of_mass,& 26 geometrical_center,& 27 get_scaled_cell,& 28 nearest_distance 29 USE tmc_move_types, ONLY: & 30 mv_type_MD, mv_type_atom_swap, mv_type_atom_trans, mv_type_gausian_adapt, mv_type_mol_rot, & 31 mv_type_mol_trans, mv_type_proton_reorder, mv_type_volume_move, tmc_move_type 32 USE tmc_tree_types, ONLY: status_frozen,& 33 status_ok,& 34 tree_type 35 USE tmc_types, ONLY: tmc_atom_type,& 36 tmc_param_type 37#include "../base/base_uses.f90" 38 39 IMPLICIT NONE 40 41 PRIVATE 42 43 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_moves' 44 45 PUBLIC :: change_pos 46 PUBLIC :: elements_in_new_subbox 47 48 INTEGER, PARAMETER :: not_selected = 0 49 INTEGER, PARAMETER :: proton_donor = -1 50 INTEGER, PARAMETER :: proton_acceptor = 1 51 52CONTAINS 53! ************************************************************************************************** 54!> \brief applying the preselected move type 55!> \param tmc_params TMC parameters with dimensions ... 56!> \param move_types ... 57!> \param rng_stream random number stream 58!> \param elem configuration to change 59!> \param mv_conf temperature index for determinig the move size 60!> \param new_subbox flag if new sub box should be crated 61!> \param move_rejected return flag if during configurational change 62!> configuration should still be accepted (not if e.g. atom/molecule 63!> leave the sub box 64!> \author Mandes 12.2012 65! ************************************************************************************************** 66 SUBROUTINE change_pos(tmc_params, move_types, rng_stream, elem, mv_conf, & 67 new_subbox, move_rejected) 68 TYPE(tmc_param_type), POINTER :: tmc_params 69 TYPE(tmc_move_type), POINTER :: move_types 70 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 71 TYPE(tree_type), POINTER :: elem 72 INTEGER :: mv_conf 73 LOGICAL :: new_subbox, move_rejected 74 75 CHARACTER(LEN=*), PARAMETER :: routineN = 'change_pos', routineP = moduleN//':'//routineN 76 77 INTEGER :: act_nr_elem_mv, counter, d, i, ind, & 78 ind_e, m, nr_molec, nr_sub_box_elem 79 INTEGER, DIMENSION(:), POINTER :: mol_in_sb 80 REAL(KIND=dp) :: rnd 81 REAL(KIND=dp), DIMENSION(:), POINTER :: direction, elem_center 82 83 NULLIFY (direction, elem_center, mol_in_sb) 84 85 CPASSERT(ASSOCIATED(tmc_params)) 86 CPASSERT(ASSOCIATED(move_types)) 87 CPASSERT(ASSOCIATED(elem)) 88 89 move_rejected = .FALSE. 90 91 CALL rng_stream%set(bg=elem%rng_seed(:, :, 1), & 92 cg=elem%rng_seed(:, :, 2), ig=elem%rng_seed(:, :, 3)) 93 94 IF (new_subbox) THEN 95 IF (ALL(tmc_params%sub_box_size .GT. 0.0_dp)) THEN 96 CALL elements_in_new_subbox(tmc_params=tmc_params, & 97 rng_stream=rng_stream, elem=elem, & 98 nr_of_sub_box_elements=nr_sub_box_elem) 99 ELSE 100 elem%elem_stat(:) = status_ok 101 END IF 102 END IF 103 104 ! at least one atom should be in the sub box 105 CPASSERT(ANY(elem%elem_stat(:) .EQ. status_ok)) 106 IF (tmc_params%nr_elem_mv .EQ. 0) THEN 107 ! move all elements (could be all atoms or all molecules) 108 act_nr_elem_mv = 0 109 ELSE 110 act_nr_elem_mv = tmc_params%nr_elem_mv 111 END IF 112 !-- select the type of move (looked up in list, using the move type index) 113 !-- for each move type exist single moves of certain number of elements 114 !-- or move of all elements 115 !-- one element is a position or velocity of an atom. 116 !-- Always all dimension are changed. 117 SELECT CASE (elem%move_type) 118 CASE (mv_type_gausian_adapt) 119 ! just for Gaussian Adaptation 120 CPABORT("gaussian adaptation is not imlemented yet.") 121!TODO CALL new_pos_gauss_adapt(acc=ASSOCIATED(elem%parent%acc, elem), & 122! pos=elem%pos, covari=elem%frc, pot=elem%potential, & 123! step_size=elem%ekin, pos_aver=elem%vel, temp=elem%ekin_before_md, & 124! rng_seed=elem%rng_seed, rng_seed_last_acc=last_acc_elem%rng_seed) 125 !-- atom translation 126 CASE (mv_type_atom_trans) 127 IF (act_nr_elem_mv .EQ. 0) & 128 act_nr_elem_mv = SIZE(elem%pos)/tmc_params%dim_per_elem 129 ALLOCATE (elem_center(tmc_params%dim_per_elem)) 130 i = 1 131 move_elements_loop: DO 132 ! select atom 133 IF (tmc_params%nr_elem_mv .EQ. 0) THEN 134 ind = (i - 1)*(tmc_params%dim_per_elem) + 1 135 ELSE 136 rnd = rng_stream%next() 137 ind = tmc_params%dim_per_elem* & 138 INT(rnd*(SIZE(elem%pos)/tmc_params%dim_per_elem)) + 1 139 END IF 140 ! apply move 141 IF (elem%elem_stat(ind) .EQ. status_ok) THEN 142 ! displace atom 143 DO d = 0, tmc_params%dim_per_elem - 1 144 rnd = rng_stream%next() 145 elem%pos(ind + d) = elem%pos(ind + d) + (rnd - 0.5)*2.0* & 146 move_types%mv_size(mv_type_atom_trans, mv_conf) 147 END DO 148 ! check if new position is in subbox 149 elem_center = elem%pos(ind:ind + tmc_params%dim_per_elem - 1) 150 IF (.NOT. check_pos_in_subbox(pos=elem_center, & 151 subbox_center=elem%subbox_center, & 152 box_scale=elem%box_scale, tmc_params=tmc_params) & 153 ) THEN 154 move_rejected = .TRUE. 155 EXIT move_elements_loop 156 END IF 157 ELSE 158 ! element was not in sub box, search new one instead 159 IF (tmc_params%nr_elem_mv .GT. 0) i = i - 1 160 END IF 161 i = i + 1 162 IF (i .GT. act_nr_elem_mv) EXIT move_elements_loop 163 END DO move_elements_loop 164 DEALLOCATE (elem_center) 165 166 !-- molecule translation 167 CASE (mv_type_mol_trans) 168 nr_molec = MAXVAL(elem%mol(:)) 169 ! if all particles should be displaced, set the amount of molecules 170 IF (act_nr_elem_mv .EQ. 0) & 171 act_nr_elem_mv = nr_molec 172 ALLOCATE (mol_in_sb(nr_molec)) 173 ALLOCATE (elem_center(tmc_params%dim_per_elem)) 174 mol_in_sb(:) = status_frozen 175 ! check if any molecule is in sub_box 176 DO m = 1, nr_molec 177 CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, & 178 start_ind=ind, end_ind=ind_e) 179 CALL geometrical_center(pos=elem%pos(ind:ind_e + tmc_params%dim_per_elem - 1), & 180 center=elem_center) 181 IF (check_pos_in_subbox(pos=elem_center, & 182 subbox_center=elem%subbox_center, & 183 box_scale=elem%box_scale, tmc_params=tmc_params) & 184 ) & 185 mol_in_sb(m) = status_ok 186 END DO 187 ! displace the selected amount of molecules 188 IF (ANY(mol_in_sb(:) .EQ. status_ok)) THEN 189 ALLOCATE (direction(tmc_params%dim_per_elem)) 190 counter = 1 191 move_molecule_loop: DO 192 ! select molecule 193 IF (tmc_params%nr_elem_mv .EQ. 0) THEN 194 m = counter 195 ELSE 196 rnd = rng_stream%next() 197 m = INT(rnd*nr_molec) + 1 198 END IF 199 CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, & 200 start_ind=ind, end_ind=ind_e) 201 ! when "molecule" is single atom, search a new one 202 IF (ind .EQ. ind_e) CYCLE move_molecule_loop 203 204 ! calculate displacement 205 ! move only molecules, with geom. center in subbox 206 IF (mol_in_sb(m) .EQ. status_ok) THEN 207 ! calculate displacement 208 DO d = 1, tmc_params%dim_per_elem 209 rnd = rng_stream%next() 210 direction(d) = (rnd - 0.5)*2.0_dp*move_types%mv_size( & 211 mv_type_mol_trans, mv_conf) 212 END DO 213 ! check if displaced position is still in subbox 214 elem_center(:) = elem_center(:) + direction(:) 215 IF (check_pos_in_subbox(pos=elem_center, & 216 subbox_center=elem%subbox_center, & 217 box_scale=elem%box_scale, tmc_params=tmc_params) & 218 ) THEN 219 ! apply move 220 atom_in_mol_loop: DO i = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem 221 dim_loop: DO d = 0, tmc_params%dim_per_elem - 1 222 elem%pos(i + d) = elem%pos(i + d) + direction(d + 1) 223 END DO dim_loop 224 END DO atom_in_mol_loop 225 ELSE 226 ! the whole move is rejected, because one element is outside the subbox 227 move_rejected = .TRUE. 228 EXIT move_molecule_loop 229 END IF 230 ELSE 231 ! element was not in sub box, search new one instead 232 IF (tmc_params%nr_elem_mv .GT. 0) counter = counter - 1 233 END IF 234 counter = counter + 1 235 IF (counter .GT. act_nr_elem_mv) EXIT move_molecule_loop 236 END DO move_molecule_loop 237 DEALLOCATE (direction) 238 END IF 239 DEALLOCATE (elem_center) 240 DEALLOCATE (mol_in_sb) 241 242 !-- molecule rotation 243 CASE (mv_type_mol_rot) 244 nr_molec = MAXVAL(elem%mol(:)) 245 IF (act_nr_elem_mv .EQ. 0) & 246 act_nr_elem_mv = nr_molec 247 ALLOCATE (mol_in_sb(nr_molec)) 248 ALLOCATE (elem_center(tmc_params%dim_per_elem)) 249 mol_in_sb(:) = status_frozen 250 ! check if any molecule is in sub_box 251 DO m = 1, nr_molec 252 CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, & 253 start_ind=ind, end_ind=ind_e) 254 CALL geometrical_center(pos=elem%pos(ind:ind_e + tmc_params%dim_per_elem - 1), & 255 center=elem_center) 256 IF (check_pos_in_subbox(pos=elem_center, & 257 subbox_center=elem%subbox_center, & 258 box_scale=elem%box_scale, tmc_params=tmc_params) & 259 ) & 260 mol_in_sb(m) = status_ok 261 END DO 262 ! rotate the selected amount of molecules 263 IF (ANY(mol_in_sb(:) .EQ. status_ok)) THEN 264 counter = 1 265 rot_molecule_loop: DO 266 ! select molecule 267 IF (tmc_params%nr_elem_mv .EQ. 0) THEN 268 m = counter 269 ELSE 270 rnd = rng_stream%next() 271 m = INT(rnd*nr_molec) + 1 272 END IF 273 CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, & 274 start_ind=ind, end_ind=ind_e) 275 ! when "molecule" is single atom, search a new one 276 IF (ind .EQ. ind_e) CYCLE rot_molecule_loop 277 278 ! apply move 279 IF (mol_in_sb(m) .EQ. status_ok) THEN 280 CALL do_mol_rot(pos=elem%pos, ind_start=ind, ind_end=ind_e, & 281 max_angle=move_types%mv_size( & 282 mv_type_mol_rot, mv_conf), & 283 move_types=move_types, rng_stream=rng_stream, & 284 dim_per_elem=tmc_params%dim_per_elem) 285 ! update sub box status of single atom 286 DO i = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem 287 elem_center = elem%pos(i:i + tmc_params%dim_per_elem - 1) 288 IF (check_pos_in_subbox(pos=elem_center, & 289 subbox_center=elem%subbox_center, & 290 box_scale=elem%box_scale, tmc_params=tmc_params) & 291 ) THEN 292 elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_ok 293 ELSE 294 elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_frozen 295 END IF 296 END DO 297 ELSE 298 ! element was not in sub box, search new one instead 299 IF (tmc_params%nr_elem_mv .GT. 0) counter = counter - 1 300 END IF 301 counter = counter + 1 302 IF (counter .GT. act_nr_elem_mv) EXIT rot_molecule_loop 303 END DO rot_molecule_loop 304 END IF 305 DEALLOCATE (elem_center) 306 DEALLOCATE (mol_in_sb) 307 308 !-- velocity changes for MD 309 !-- here all velocities are changed 310 CASE (mv_type_MD) 311 CPASSERT(ASSOCIATED(tmc_params%atoms)) 312 change_all_velocities_loop: DO i = 1, SIZE(elem%pos) 313 !-- attention, move type size is in atomic units of velocity 314 IF (elem%elem_stat(i) .NE. status_frozen) THEN 315 CALL vel_change(vel=elem%vel(i), & 316 atom_kind=tmc_params%atoms(INT(i/REAL(tmc_params%dim_per_elem, KIND=dp)) + 1), & 317 phi=move_types%mv_size(mv_type_MD, 1), & ! TODO parallel tempering move sizes for vel_change 318 temp=tmc_params%Temp(mv_conf), & 319 rnd_sign_change=.TRUE., & ! MD_vel_invert, & 320 rng_stream=rng_stream) 321 END IF 322 END DO change_all_velocities_loop 323 324 !-- proton order and disorder 325 ! a loop of molecules is build an in this loop proton acceptors become proton donators 326 ! Therefor the molecules are rotated along the not involved O-H bond 327 CASE (mv_type_proton_reorder) 328 CALL search_and_do_proton_displace_loop(elem=elem, & 329 short_loop=move_rejected, rng_stream=rng_stream, & 330 tmc_params=tmc_params) 331 332 !-- volume move 333 ! the box is increased or decreased and with it the coordinates 334 CASE (mv_type_volume_move) 335 CALL change_volume(conf=elem, T_ind=mv_conf, move_types=move_types, & 336 rng_stream=rng_stream, tmc_params=tmc_params, & 337 mv_cen_of_mass=tmc_params%mv_cen_of_mass) 338 339 !-- atom swap 340 ! two atoms of different types are swapped 341 CASE (mv_type_atom_swap) 342 CALL swap_atoms(conf=elem, move_types=move_types, rng_stream=rng_stream, & 343 tmc_params=tmc_params) 344 345 CASE DEFAULT 346 CALL cp_abort(__LOCATION__, & 347 "unknown move type "// & 348 cp_to_string(elem%move_type)) 349 END SELECT 350 351 CALL rng_stream%get(bg=elem%rng_seed(:, :, 1), & 352 cg=elem%rng_seed(:, :, 2), ig=elem%rng_seed(:, :, 3)) 353 354 END SUBROUTINE change_pos 355 356! ************************************************************************************************** 357!> \brief gets the index of the first molecule element position and the size 358!> \param tmc_params TMC parameters with dim_per_elem 359!> \param mol_arr array with molecule information (which atom attend which mol) 360!> \param mol the selected molecule number 361!> \param start_ind start index of the first atom in molecule 362!> \param end_ind index of the last atom in molecule 363!> \author Mandes 10.2013 364! ************************************************************************************************** 365 SUBROUTINE get_mol_indeces(tmc_params, mol_arr, mol, start_ind, end_ind) 366 TYPE(tmc_param_type), POINTER :: tmc_params 367 INTEGER, DIMENSION(:), INTENT(IN), POINTER :: mol_arr 368 INTEGER, INTENT(IN) :: mol 369 INTEGER, INTENT(OUT) :: start_ind, end_ind 370 371 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mol_indeces', & 372 routineP = moduleN//':'//routineN 373 374 INTEGER :: i 375 376 start_ind = -1 377 end_ind = -1 378 379 CPASSERT(ASSOCIATED(mol_arr)) 380 CPASSERT(mol .LE. MAXVAL(mol_arr(:))) 381 ! get start index 382 loop_start: DO i = 1, SIZE(mol_arr) 383 IF (mol_arr(i) .EQ. mol) THEN 384 start_ind = i 385 EXIT loop_start 386 END IF 387 END DO loop_start 388 ! get end index 389 loop_end: DO i = SIZE(mol_arr), i, -1 390 IF (mol_arr(i) .EQ. mol) THEN 391 end_ind = i 392 EXIT loop_end 393 END IF 394 END DO loop_end 395 ! check if all atoms inbetween attend to molecule 396 CPASSERT(ALL(mol_arr(start_ind:end_ind) .EQ. mol)) 397 CPASSERT(start_ind .GT. 0) 398 CPASSERT(end_ind .GT. 0) 399 ! convert to indeces mapped for the position array (multiple dim per atom) 400 start_ind = (start_ind - 1)*tmc_params%dim_per_elem + 1 401 end_ind = (end_ind - 1)*tmc_params%dim_per_elem + 1 402 END SUBROUTINE 403 404! ************************************************************************************************** 405!> \brief checks if a position is within the sub box 406!> returns true if position is inside 407!> \param pos array with positions 408!> \param subbox_center actual center of sub box 409!> \param box_scale scaling factors for the cell 410!> \param tmc_params TMC parameters with sub box size and cell 411!> \return ... 412!> \author Mandes 11.2012 413! ************************************************************************************************** 414 FUNCTION check_pos_in_subbox(pos, subbox_center, box_scale, tmc_params) & 415 RESULT(inside) 416 REAL(KIND=dp), DIMENSION(:), POINTER :: pos, subbox_center, box_scale 417 TYPE(tmc_param_type), POINTER :: tmc_params 418 LOGICAL :: inside 419 420 CHARACTER(LEN=*), PARAMETER :: routineN = 'check_pos_in_subbox', & 421 routineP = moduleN//':'//routineN 422 423 INTEGER :: handle 424 LOGICAL :: flag 425 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pos_tmp 426 427 CPASSERT(ASSOCIATED(pos)) 428 CPASSERT(ASSOCIATED(subbox_center)) 429 CPASSERT(ASSOCIATED(box_scale)) 430 ! if pressure is defined, no scale should be 0 431 flag = .NOT. ((tmc_params%pressure .GT. 0.0_dp) .AND. (ANY(box_scale .EQ. 0.0_dp))) 432 CPASSERT(flag) 433 CPASSERT(SIZE(pos) .EQ. 3) 434 CPASSERT(SIZE(pos) .EQ. SIZE(subbox_center)) 435 436 ! start the timing 437 CALL timeset(routineN, handle) 438 439 ALLOCATE (pos_tmp(SIZE(pos))) 440 441 inside = .TRUE. 442 ! return if no subbox is defined 443 IF (.NOT. ANY(tmc_params%sub_box_size(:) .LE. 0.1_dp)) THEN 444 pos_tmp(:) = pos(:) - subbox_center(:) 445 CALL get_scaled_cell(cell=tmc_params%cell, box_scale=box_scale, & 446 vec=pos_tmp) 447 ! check 448 IF (ANY(pos_tmp(:) .GE. tmc_params%sub_box_size(:)/2.0) .OR. & 449 ANY(pos_tmp(:) .LE. -tmc_params%sub_box_size(:)/2.0)) THEN 450 inside = .FALSE. 451 END IF 452 END IF 453 DEALLOCATE (pos_tmp) 454 ! end the timing 455 CALL timestop(handle) 456 END FUNCTION check_pos_in_subbox 457 458! ************************************************************************************************** 459!> \brief set a new random sub box center and counte the number of atoms in it 460!> \param tmc_params ... 461!> \param rng_stream ... 462!> \param elem ... 463!> \param nr_of_sub_box_elements ... 464!> \param 465!> \param 466!> \author Mandes 11.2012 467! ************************************************************************************************** 468 SUBROUTINE elements_in_new_subbox(tmc_params, rng_stream, elem, & 469 nr_of_sub_box_elements) 470 TYPE(tmc_param_type), POINTER :: tmc_params 471 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 472 TYPE(tree_type), POINTER :: elem 473 INTEGER, INTENT(OUT) :: nr_of_sub_box_elements 474 475 CHARACTER(LEN=*), PARAMETER :: routineN = 'elements_in_new_subbox', & 476 routineP = moduleN//':'//routineN 477 478 INTEGER :: handle, i 479 REAL(KIND=dp) :: rnd 480 REAL(KIND=dp), DIMENSION(3) :: box_size 481 REAL(KIND=dp), DIMENSION(:), POINTER :: atom_tmp, center_of_sub_box 482 483 NULLIFY (center_of_sub_box, atom_tmp) 484 485 CPASSERT(ASSOCIATED(tmc_params)) 486 CPASSERT(ASSOCIATED(elem)) 487 488 ! start the timing 489 CALL timeset(routineN, handle) 490 491 IF (ANY(tmc_params%sub_box_size(:) .LE. 0.1_dp)) THEN 492 !CPWARN("try to count elements in sub box without sub box.") 493 elem%elem_stat = status_ok 494 nr_of_sub_box_elements = SIZE(elem%elem_stat) 495 ELSE 496 ALLOCATE (center_of_sub_box(tmc_params%dim_per_elem)) 497 ALLOCATE (atom_tmp(tmc_params%dim_per_elem)) 498 nr_of_sub_box_elements = 0 499 ! -- define the center of the sub box 500 CALL rng_stream%set(bg=elem%rng_seed(:, :, 1), cg=elem%rng_seed(:, :, 2), & 501 ig=elem%rng_seed(:, :, 3)) 502 503 CALL get_cell(cell=tmc_params%cell, abc=box_size) 504 DO i = 1, SIZE(tmc_params%sub_box_size) 505 rnd = rng_stream%next() 506 center_of_sub_box(i) = rnd*box_size(i) 507 END DO 508 elem%subbox_center(:) = center_of_sub_box(:) 509 510 CALL rng_stream%get(bg=elem%rng_seed(:, :, 1), cg=elem%rng_seed(:, :, 2), & 511 ig=elem%rng_seed(:, :, 3)) 512 513 ! check all elements if they are in subbox 514 DO i = 1, SIZE(elem%pos), tmc_params%dim_per_elem 515 atom_tmp(:) = elem%pos(i:i + tmc_params%dim_per_elem - 1) 516 IF (check_pos_in_subbox(pos=atom_tmp, & 517 subbox_center=center_of_sub_box, box_scale=elem%box_scale, & 518 tmc_params=tmc_params)) THEN 519 elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_ok 520 nr_of_sub_box_elements = nr_of_sub_box_elements + 1 521 ELSE 522 elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_frozen 523 END IF 524 END DO 525 DEALLOCATE (atom_tmp) 526 DEALLOCATE (center_of_sub_box) 527 END IF 528 ! end the timing 529 CALL timestop(handle) 530 END SUBROUTINE elements_in_new_subbox 531 532! ************************************************************************************************** 533!> \brief molecule rotation using quaternions 534!> \param pos atom positions 535!> \param ind_start starting index in the array 536!> \param ind_end index of last atom in the array 537!> \param max_angle maximal angle in each direction 538!> \param move_types ... 539!> \param rng_stream ramdon stream 540!> \param dim_per_elem dimension per atom 541!> \author Mandes 11.2012 542! ************************************************************************************************** 543 SUBROUTINE do_mol_rot(pos, ind_start, ind_end, max_angle, move_types, & 544 rng_stream, dim_per_elem) 545 REAL(KIND=dp), DIMENSION(:), POINTER :: pos 546 INTEGER :: ind_start, ind_end 547 REAL(KIND=dp) :: max_angle 548 TYPE(tmc_move_type), POINTER :: move_types 549 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 550 INTEGER :: dim_per_elem 551 552 CHARACTER(LEN=*), PARAMETER :: routineN = 'do_mol_rot', routineP = moduleN//':'//routineN 553 554 INTEGER :: i 555 REAL(KIND=dp) :: a1, a2, a3, q0, q1, q2, q3, rnd 556 REAL(KIND=dp), DIMENSION(3, 3) :: rot 557 REAL(KIND=dp), DIMENSION(:), POINTER :: elem_center 558 559 NULLIFY (elem_center) 560 561 CPASSERT(ASSOCIATED(pos)) 562 CPASSERT(dim_per_elem .EQ. 3) 563 CPASSERT(ind_start .GT. 0 .AND. ind_start .LT. SIZE(pos)) 564 CPASSERT(ind_end .GT. 0 .AND. ind_end .LT. SIZE(pos)) 565 CPASSERT(ASSOCIATED(move_types)) 566 567 ! calculate rotation matrix (using quanternions) 568 rnd = rng_stream%next() 569 a1 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf) 570 rnd = rng_stream%next() 571 a2 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf) 572 rnd = rng_stream%next() 573 a3 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf) 574 q0 = COS(a2/2)*COS((a1 + a3)/2.0_dp) 575 q1 = SIN(a2/2)*COS((a1 - a3)/2.0_dp) 576 q2 = SIN(a2/2)*SIN((a1 - a3)/2.0_dp) 577 q3 = COS(a2/2)*SIN((a1 + a3)/2.0_dp) 578 rot = RESHAPE((/q0*q0 + q1*q1 - q2*q2 - q3*q3, 2*(q1*q2 - q0*q3), 2*(q1*q3 + q0*q2), & 579 2*(q1*q2 + q0*q3), q0*q0 - q1*q1 + q2*q2 - q3*q3, 2*(q2*q3 - q0*q1), & 580 2*(q1*q3 - q0*q2), 2*(q2*q3 + q0*q1), q0*q0 - q1*q1 - q2*q2 + q3*q3/), (/3, 3/)) 581 582 ALLOCATE (elem_center(dim_per_elem)) 583 ! calculate geometrical center 584 CALL geometrical_center(pos=pos(ind_start:ind_end + dim_per_elem - 1), & 585 center=elem_center) 586 587 ! proceed rotation 588 atom_loop: DO i = ind_start, ind_end + dim_per_elem - 1, dim_per_elem 589 pos(i:i + 2) = MATMUL(pos(i:i + 2) - elem_center(:), rot) + elem_center(:) 590 END DO atom_loop 591 DEALLOCATE (elem_center) 592 END SUBROUTINE do_mol_rot 593 594! ************************************************************************************************** 595!> \brief velocity change should be gaussian distributed 596!> around the old velocity with respect to kB*T/m 597!> \param vel velocity of atom (one direction) 598!> \param atom_kind ... 599!> \param phi angle for mixing old with random gaussian distributed velocity 600!> phi =90 degree -> only gaussian velocity around 0 601!> phi = 0 degree -> only old velocity (with sign change) 602!> \param temp temperature for gaussian distributed velocity 603!> \param rnd_sign_change if sign of old velocity should change randomly 604!> \param rng_stream random number stream 605!> \author Mandes 11.2012 606! ************************************************************************************************** 607 SUBROUTINE vel_change(vel, atom_kind, phi, temp, rnd_sign_change, rng_stream) 608 REAL(KIND=dp), INTENT(INOUT) :: vel 609 TYPE(tmc_atom_type) :: atom_kind 610 REAL(KIND=dp), INTENT(IN) :: phi, temp 611 LOGICAL :: rnd_sign_change 612 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 613 614 CHARACTER(LEN=*), PARAMETER :: routineN = 'vel_change', routineP = moduleN//':'//routineN 615 616 INTEGER :: d 617 REAL(KIND=dp) :: delta_vel, kB, rnd1, rnd2, rnd3, rnd_g 618 619 kB = boltzmann/joule 620 621 !phi = move_types%mv_size(mv_type_MD,1) ! TODO parallel tempering move sizes for vel_change 622 ! hence first producing a gaussian random number 623 rnd1 = rng_stream%next() 624 rnd2 = rng_stream%next() 625 626 rnd_g = SQRT(-2.0_dp*LOG(rnd1))*COS(2.0_dp*PI*rnd2) 627 !we can also produce a second one in the same step: 628 !rnd_g2 = SQRT(-2.0_dp*LOG(rnd1))*SIN(2.0_dp*PI*rnd2) 629 630 ! adapting the variance with respect to kB*T/m 631 delta_vel = SQRT(kB*temp/atom_kind%mass)*rnd_g 632 ! check if TODO random velocity sign change 633 ! using detailed balance, velocity sign changes are necessary, 634 ! which are done randomly and 635 ! can be switched of using MD_vel_invert 636 ! without still the balance condition should be fulfilled 637 638 rnd3 = rng_stream%next() 639 IF (rnd3 .GE. 0.5 .AND. rnd_sign_change) THEN 640 d = -1 641 ELSE 642 d = 1 643 END IF 644 vel = SIN(phi)*delta_vel + COS(phi)*vel*d*1.0_dp 645 END SUBROUTINE vel_change 646 647! ************************************************************************************************** 648!> \brief proton order and disorder (specialized move for ice Ih) 649!> a loop of molecules is build an 650!> in this loop proton acceptors become proton donators 651!> Therefor the molecules are rotated along the not involved O-H bond 652!> \param elem sub tree element with actual positions 653!> \param short_loop return if the a loop shorter than 6 molecules is found 654!> (should not be in ice structure) 655!> \param rng_stream random number stream 656!> \param tmc_params TMC parameters with numbers of dimensions per element 657!> number of atoms per molecule 658!> \author Mandes 11.2012 659! ************************************************************************************************** 660 SUBROUTINE search_and_do_proton_displace_loop(elem, short_loop, rng_stream, & 661 tmc_params) 662 TYPE(tree_type), POINTER :: elem 663 LOGICAL :: short_loop 664 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 665 TYPE(tmc_param_type), POINTER :: tmc_params 666 667 CHARACTER(LEN=*), PARAMETER :: routineN = 'search_and_do_proton_displace_loop', & 668 routineP = moduleN//':'//routineN 669 670 CHARACTER(LEN=1000) :: tmp_chr 671 INTEGER :: counter, donor_acceptor, handle, k, mol, & 672 nr_mol 673 INTEGER, DIMENSION(:), POINTER :: mol_arr 674 REAL(KIND=dp) :: rnd 675 676 NULLIFY (mol_arr) 677 678 CPASSERT(ASSOCIATED(elem)) 679 CPASSERT(ASSOCIATED(tmc_params)) 680 681 ! start the timing 682 CALL timeset(routineN, handle) 683 684 short_loop = .FALSE. 685 counter = 0 686 nr_mol = MAXVAL(elem%mol(:)) 687 ! ind_arr: one array element for each molecule 688 ALLOCATE (mol_arr(nr_mol)) 689 mol_arr(:) = -1 690 donor_acceptor = not_selected 691 ! select randomly if neighboring molecule is donor / acceptor 692 IF (rng_stream%next() .LT. 0.5_dp) THEN 693 donor_acceptor = proton_acceptor 694 ELSE 695 donor_acceptor = proton_donor 696 END IF 697 698 ! first step build loop 699 ! select randomly one atom 700 rnd = rng_stream%next() 701 ! the randomly selected first atom 702 mol = INT(rnd*nr_mol) + 1 703 counter = counter + 1 704 mol_arr(counter) = mol 705 706 ! do until the loop is closed 707 ! (until path connects back to any spot of the path) 708 chain_completition_loop: DO 709 counter = counter + 1 710 ! find nearest neighbor 711 ! (with same state, in the chain, proton donator or proton accptor) 712 CALL find_nearest_proton_acceptor_donator(elem=elem, mol=mol, & 713 donor_acceptor=donor_acceptor, tmc_params=tmc_params, & 714 rng_stream=rng_stream) 715 IF (ANY(mol_arr(:) .EQ. mol)) & 716 EXIT chain_completition_loop 717 mol_arr(counter) = mol 718 END DO chain_completition_loop 719 counter = counter - 1 ! last searched element is equal to one other in list 720 721 ! just take the loop of molecules out of the chain 722 DO k = 1, counter 723 IF (mol_arr(k) .EQ. mol) & 724 EXIT 725 END DO 726 mol_arr(1:counter - k + 1) = mol_arr(k:counter) 727 counter = counter - k + 1 728 729 ! check if loop is minimum size of 6 molecules 730 IF (counter .LT. 6) THEN 731 CALL cp_warn(__LOCATION__, & 732 "short proton loop with"//cp_to_string(counter)// & 733 "molecules.") 734 tmp_chr = "" 735 WRITE (tmp_chr, *) mol_arr(1:counter) 736 CPWARN("selected molecules:"//TRIM(tmp_chr)) 737 short_loop = .TRUE. 738 ENDIF 739 740 ! rotate the molecule along the not involved O-H bond 741 ! (about the angle in of the neighboring chain elements) 742 CALL rotate_molecules_in_chain(tmc_params=tmc_params, elem=elem, & 743 mol_arr_in=mol_arr(1:counter), donor_acceptor=donor_acceptor) 744 DEALLOCATE (mol_arr) 745 746 ! end the timing 747 CALL timestop(handle) 748 END SUBROUTINE search_and_do_proton_displace_loop 749 750! ************************************************************************************************** 751!> \brief searches the next (first atom of) neighboring molecule 752!> which is proton donor / acceptor 753!> \param elem sub tree element with actual positions 754!> \param mol (in_out) actual regarded molecule, which neighbor is searched for 755!> \param donor_acceptor type of searched neighbor 756!> (proton donor or proton acceptor) 757!> \param tmc_params TMC parameters with numbers of dimensions per element 758!> number of atoms per molecule 759!> \param rng_stream random number stream 760!> \author Mandes 12.2012 761! ************************************************************************************************** 762 SUBROUTINE find_nearest_proton_acceptor_donator(elem, mol, donor_acceptor, & 763 tmc_params, rng_stream) 764 TYPE(tree_type), POINTER :: elem 765 INTEGER :: mol, donor_acceptor 766 TYPE(tmc_param_type), POINTER :: tmc_params 767 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 768 769 CHARACTER(LEN=*), PARAMETER :: routineN = 'find_nearest_proton_acceptor_donator', & 770 routineP = moduleN//':'//routineN 771 772 INTEGER :: handle, ind, ind_e, ind_n, mol_tmp, & 773 nr_mol 774 INTEGER, DIMENSION(2) :: neighbor_mol 775 REAL(KIND=dp) :: dist_tmp, rnd 776 REAL(KIND=dp), DIMENSION(:), POINTER :: distH1, distH2, distO 777 778 NULLIFY (distO, distH1, distH2) 779 CPASSERT(ASSOCIATED(elem)) 780 CPASSERT(ASSOCIATED(tmc_params)) 781 782 ! start the timing 783 CALL timeset(routineN, handle) 784 785 nr_mol = MAXVAL(elem%mol) 786 ALLOCATE (distO(nr_mol)) 787 ALLOCATE (distH1(nr_mol)) 788 ALLOCATE (distH2(nr_mol)) 789 !-- initialize the distances to huge values 790 ! distance of nearest proton of certain molecule to preselected O 791 distO(:) = HUGE(distO(1)) 792 ! distance of (first) proton of preselected molecule to certain molecule 793 distH1(:) = HUGE(distH1(1)) 794 ! distance of (second) proton of preselected molecule to certain molecule 795 distH2(:) = HUGE(distH2(1)) 796 797 ! get the indices of the old O atom (assuming the first atom of the molecule the first atom) 798 CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=mol, & 799 start_ind=ind, end_ind=ind_e) 800 801 ! calculate distances to all molecules 802 list_distances: DO mol_tmp = 1, nr_mol 803 IF (mol_tmp .EQ. mol) CYCLE list_distances 804 ! index of the molecule (the O atom) 805 ! assume the first atom of the molecule the first atom 806 CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, & 807 mol=mol_tmp, start_ind=ind_n, end_ind=ind_e) 808 ! check if selected molecule is water respectively consists of 3 atoms 809 IF (MOD(ind_e - ind_n, 3) .GT. 0) THEN 810 CALL cp_warn(__LOCATION__, & 811 "selected a molecule with more than 3 atoms, "// & 812 "the proton reordering does not support, skip molecule") 813 CYCLE list_distances 814 END IF 815 IF (donor_acceptor .EQ. proton_acceptor) THEN 816 IF (check_donor_acceptor(elem=elem, i_orig=ind, i_neighbor=ind_n, & 817 tmc_params=tmc_params) .EQ. proton_acceptor) THEN 818 !distance of fist proton to certain O 819 distH1(mol_tmp) = nearest_distance( & 820 x1=elem%pos(ind + tmc_params%dim_per_elem: & 821 ind + 2*tmc_params%dim_per_elem - 1), & 822 x2=elem%pos(ind_n:ind_n + tmc_params%dim_per_elem - 1), & 823 cell=tmc_params%cell, box_scale=elem%box_scale) 824 !distance of second proton to certain O 825 distH2(mol_tmp) = nearest_distance( & 826 x1=elem%pos(ind + 2*tmc_params%dim_per_elem: & 827 ind + 3*tmc_params%dim_per_elem - 1), & 828 x2=elem%pos(ind_n:ind_n + tmc_params%dim_per_elem - 1), & 829 cell=tmc_params%cell, box_scale=elem%box_scale) 830 END IF 831 END IF 832 !check for neighboring proton donors 833 IF (donor_acceptor .EQ. proton_donor) THEN 834 IF (check_donor_acceptor(elem=elem, i_orig=ind, i_neighbor=ind_n, & 835 tmc_params=tmc_params) .EQ. proton_donor) THEN 836 !distance of selected O to all first protons of other melecules 837 distO(mol_tmp) = nearest_distance( & 838 x1=elem%pos(ind:ind + tmc_params%dim_per_elem - 1), & 839 x2=elem%pos(ind_n + tmc_params%dim_per_elem: & 840 ind_n + 2*tmc_params%dim_per_elem - 1), & 841 cell=tmc_params%cell, box_scale=elem%box_scale) 842 dist_tmp = nearest_distance( & 843 x1=elem%pos(ind:ind + tmc_params%dim_per_elem - 1), & 844 x2=elem%pos(ind_n + 2*tmc_params%dim_per_elem: & 845 ind_n + 3*tmc_params%dim_per_elem - 1), & 846 cell=tmc_params%cell, box_scale=elem%box_scale) 847 IF (dist_tmp .LT. distO(mol_tmp)) distO(mol_tmp) = dist_tmp 848 END IF 849 END IF 850 END DO list_distances 851 852 mol_tmp = 1 853 ! select the nearest neighbors 854 !check for neighboring proton acceptors 855 IF (donor_acceptor .EQ. proton_acceptor) THEN 856 neighbor_mol(mol_tmp) = MINLOC(distH1(:), 1) 857 neighbor_mol(mol_tmp + 1) = MINLOC(distH2(:), 1) 858 ! if both smallest distances points to the shortest molecule search also the second next shortest distance 859 IF (neighbor_mol(mol_tmp) .EQ. neighbor_mol(mol_tmp + 1)) THEN 860 distH1(neighbor_mol(mol_tmp)) = HUGE(distH1(1)) 861 distH2(neighbor_mol(mol_tmp + 1)) = HUGE(distH2(1)) 862 IF (MINVAL(distH1(:), 1) .LT. MINVAL(distH2(:), 1)) THEN 863 neighbor_mol(mol_tmp) = MINLOC(distH1(:), 1) 864 ELSE 865 neighbor_mol(mol_tmp + 1) = MINLOC(distH2(:), 1) 866 END IF 867 END IF 868 mol_tmp = mol_tmp + 2 869 END IF 870 871 !check for neighboring proton donors 872 IF (donor_acceptor .EQ. proton_donor) THEN 873 neighbor_mol(mol_tmp) = MINLOC(distO(:), 1) 874 distO(neighbor_mol(mol_tmp)) = HUGE(distO(1)) 875 neighbor_mol(mol_tmp + 1) = MINLOC(distO(:), 1) 876 END IF 877 878 ! select randomly the next neighboring molecule 879 rnd = rng_stream%next() 880 ! the randomly selected atom: return value! 881 mol_tmp = neighbor_mol(INT(rnd*SIZE(neighbor_mol(:))) + 1) 882 mol = mol_tmp 883 884 DEALLOCATE (distO) 885 DEALLOCATE (distH1) 886 DEALLOCATE (distH2) 887 888 ! end the timing 889 CALL timestop(handle) 890 END SUBROUTINE find_nearest_proton_acceptor_donator 891 892! ************************************************************************************************** 893!> \brief checks if neighbor of the selected/orig element 894!> is a proron donator or acceptor 895!> \param elem ... 896!> \param i_orig ... 897!> \param i_neighbor ... 898!> \param tmc_params ... 899!> \return ... 900!> \author Mandes 11.2012 901! ************************************************************************************************** 902 FUNCTION check_donor_acceptor(elem, i_orig, i_neighbor, tmc_params) & 903 RESULT(donor_acceptor) 904 TYPE(tree_type), POINTER :: elem 905 INTEGER :: i_orig, i_neighbor 906 TYPE(tmc_param_type), POINTER :: tmc_params 907 INTEGER :: donor_acceptor 908 909 CHARACTER(LEN=*), PARAMETER :: routineN = 'check_donor_acceptor', & 910 routineP = moduleN//':'//routineN 911 912 REAL(KIND=dp), DIMENSION(4) :: distances 913 914 CPASSERT(ASSOCIATED(elem)) 915 CPASSERT(i_orig .GE. 1 .AND. i_orig .LE. SIZE(elem%pos)) 916 CPASSERT(i_neighbor .GE. 1 .AND. i_neighbor .LE. SIZE(elem%pos)) 917 CPASSERT(ASSOCIATED(tmc_params)) 918 919 ! 1. proton of orig with neighbor O 920 distances(1) = nearest_distance( & 921 x1=elem%pos(i_neighbor:i_neighbor + tmc_params%dim_per_elem - 1), & 922 x2=elem%pos(i_orig + tmc_params%dim_per_elem: & 923 i_orig + 2*tmc_params%dim_per_elem - 1), & 924 cell=tmc_params%cell, box_scale=elem%box_scale) 925 ! 2. proton of orig with neighbor O 926 distances(2) = nearest_distance( & 927 x1=elem%pos(i_neighbor:i_neighbor + tmc_params%dim_per_elem - 1), & 928 x2=elem%pos(i_orig + 2*tmc_params%dim_per_elem: & 929 i_orig + 3*tmc_params%dim_per_elem - 1), & 930 cell=tmc_params%cell, box_scale=elem%box_scale) 931 ! 1. proton of neighbor with orig O 932 distances(3) = nearest_distance( & 933 x1=elem%pos(i_orig:i_orig + tmc_params%dim_per_elem - 1), & 934 x2=elem%pos(i_neighbor + tmc_params%dim_per_elem: & 935 i_neighbor + 2*tmc_params%dim_per_elem - 1), & 936 cell=tmc_params%cell, box_scale=elem%box_scale) 937 ! 2. proton of neigbor with orig O 938 distances(4) = nearest_distance( & 939 x1=elem%pos(i_orig:i_orig + tmc_params%dim_per_elem - 1), & 940 x2=elem%pos(i_neighbor + 2*tmc_params%dim_per_elem: & 941 i_neighbor + 3*tmc_params%dim_per_elem - 1), & 942 cell=tmc_params%cell, box_scale=elem%box_scale) 943 944 IF (MINLOC(distances(:), 1) .LE. 2) THEN 945 donor_acceptor = proton_acceptor 946 ELSE 947 donor_acceptor = proton_donor 948 END IF 949 END FUNCTION check_donor_acceptor 950 951! ************************************************************************************************** 952!> \brief rotates all the molecules in the chain 953!> the protons were flipped from the donor to the acceptor 954!> \param tmc_params TMC environment parameters 955!> \param elem sub tree element the pos of the molecules in chain should be 956!> changed by rotating 957!> \param mol_arr_in array of indeces of molecules, should be rotated 958!> \param donor_acceptor gives the direction of rotation 959!> \author Mandes 11.2012 960! ************************************************************************************************** 961 SUBROUTINE rotate_molecules_in_chain(tmc_params, elem, mol_arr_in, & 962 donor_acceptor) 963 TYPE(tmc_param_type), POINTER :: tmc_params 964 TYPE(tree_type), POINTER :: elem 965 INTEGER, DIMENSION(:) :: mol_arr_in 966 INTEGER :: donor_acceptor 967 968 CHARACTER(LEN=*), PARAMETER :: routineN = 'rotate_molecules_in_chain', & 969 routineP = moduleN//':'//routineN 970 971 INTEGER :: H_offset, handle, i, ind 972 INTEGER, DIMENSION(:), POINTER :: ind_arr 973 REAL(KIND=dp) :: dihe_angle, dist_near, tmp 974 REAL(KIND=dp), DIMENSION(3) :: rot_axis, tmp_1, tmp_2, vec_1O, & 975 vec_2H_f, vec_2H_m, vec_2O, vec_3O, & 976 vec_4O, vec_rotated 977 TYPE(cell_type), POINTER :: tmp_cell 978 979 NULLIFY (ind_arr, tmp_cell) 980 981 CPASSERT(ASSOCIATED(tmc_params)) 982 CPASSERT(ASSOCIATED(elem)) 983 984 ! start the timing 985 CALL timeset(routineN, handle) 986 987 ALLOCATE (ind_arr(0:SIZE(mol_arr_in) + 1)) 988 DO i = 1, SIZE(mol_arr_in) 989 CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, & 990 mol=mol_arr_in(i), & 991 start_ind=ind_arr(i), end_ind=ind) 992 END DO 993 ind_arr(0) = ind_arr(SIZE(ind_arr) - 2) 994 ind_arr(SIZE(ind_arr) - 1) = ind_arr(1) 995 996 ! get the scaled cell 997 ALLOCATE (tmp_cell) 998 CALL get_scaled_cell(cell=tmc_params%cell, box_scale=elem%box_scale, & 999 scaled_cell=tmp_cell) 1000 1001 ! rotate single molecules 1002 DO i = 1, SIZE(ind_arr) - 2 1003 ! the 3 O atoms 1004 vec_1O(:) = elem%pos(ind_arr(i - 1):ind_arr(i - 1) + tmc_params%dim_per_elem - 1) 1005 vec_2O(:) = elem%pos(ind_arr(i):ind_arr(i) + tmc_params%dim_per_elem - 1) 1006 vec_3O(:) = elem%pos(ind_arr(i + 1):ind_arr(i + 1) + tmc_params%dim_per_elem - 1) 1007 ! the H atoms 1008 ! distinguished between the one fixed (rotation axis with 2 O) 1009 ! and the moved one 1010 ! if true the first H atom is between the O atoms 1011 IF (nearest_distance( & 1012 x1=elem%pos(ind_arr(i + donor_acceptor): & 1013 ind_arr(i + donor_acceptor) + tmc_params%dim_per_elem - 1), & 1014 x2=elem%pos(ind_arr(i) + tmc_params%dim_per_elem: & 1015 ind_arr(i) + 2*tmc_params%dim_per_elem - 1), & 1016 cell=tmc_params%cell, box_scale=elem%box_scale) & 1017 .LT. & 1018 nearest_distance( & 1019 x1=elem%pos(ind_arr(i + donor_acceptor): & 1020 ind_arr(i + donor_acceptor) + tmc_params%dim_per_elem - 1), & 1021 x2=elem%pos(ind_arr(i) + 2*tmc_params%dim_per_elem: & 1022 ind_arr(i) + 3*tmc_params%dim_per_elem - 1), & 1023 cell=tmc_params%cell, box_scale=elem%box_scale) & 1024 ) THEN 1025 vec_2H_m = elem%pos(ind_arr(i) + tmc_params%dim_per_elem: & 1026 ind_arr(i) + 2*tmc_params%dim_per_elem - 1) 1027 vec_2H_f = elem%pos(ind_arr(i) + 2*tmc_params%dim_per_elem: & 1028 ind_arr(i) + 3*tmc_params%dim_per_elem - 1) 1029 H_offset = 1 1030 ELSE 1031 vec_2H_f = elem%pos(ind_arr(i) + tmc_params%dim_per_elem: & 1032 ind_arr(i) + 2*tmc_params%dim_per_elem - 1) 1033 vec_2H_m = elem%pos(ind_arr(i) + 2*tmc_params%dim_per_elem: & 1034 ind_arr(i) + 3*tmc_params%dim_per_elem - 1) 1035 H_offset = 2 1036 END IF 1037 1038 IF (.TRUE.) THEN !TODO find a better switch for the pauling model 1039 1040 ! do rotation (NOT pauling model) 1041 tmp_1 = pbc(vec_2O - vec_1O, tmp_cell) 1042 tmp_2 = pbc(vec_3O - vec_2H_f, tmp_cell) 1043 1044 dihe_angle = donor_acceptor*dihedral_angle(tmp_1, vec_2H_f - vec_2O, tmp_2) 1045 DO ind = ind_arr(i), ind_arr(i) + tmc_params%dim_per_elem*3 - 1, tmc_params%dim_per_elem 1046 ! set rotation vector 1047 !vec_rotated = rotate_vector(vec_2H_m-vec_2O, dihe_angle, vec_2H_f-vec_2O) 1048 vec_rotated = rotate_vector(elem%pos(ind: & 1049 ind + tmc_params%dim_per_elem - 1) - vec_2O, & 1050 dihe_angle, vec_2H_f - vec_2O) 1051 1052 ! set new position 1053 !elem%pos(ind_arr(i)+H_offset*dim_per_elem:ind_arr(i)+(H_offset+1)*dim_per_elem-1) = vec_2O+vec_rotated 1054 elem%pos(ind:ind + tmc_params%dim_per_elem - 1) = vec_2O + vec_rotated 1055 END DO 1056 ELSE 1057 ! using the pauling model 1058 ! (see Aragones and Vega: Dielectric constant of ices...) 1059 ! the rotation axis is defined using the 4th not involved O 1060 ! (next to the not involved H) 1061 ! O atom next to not involved proton for axis calculation 1062 dist_near = HUGE(dist_near) 1063 search_O_loop: DO ind = 1, SIZE(elem%pos), & 1064 tmc_params%dim_per_elem*3 1065 IF (ind .EQ. ind_arr(i)) CYCLE search_O_loop 1066 tmp = nearest_distance(x1=vec_2H_f, & 1067 x2=elem%pos(ind:ind + tmc_params%dim_per_elem - 1), & 1068 cell=tmc_params%cell, box_scale=elem%box_scale) 1069 IF (dist_near .GT. tmp) THEN 1070 dist_near = tmp 1071 vec_4O = elem%pos(ind:ind + tmc_params%dim_per_elem - 1) 1072 END IF 1073 END DO search_O_loop 1074 rot_axis = pbc(-vec_2O(:) + vec_4O(:), tmp_cell) 1075 tmp_1 = pbc(vec_2O - vec_1O, tmp_cell) 1076 tmp_2 = pbc(vec_3O - vec_4O, tmp_cell) 1077 dihe_angle = donor_acceptor*dihedral_angle(tmp_1, rot_axis, tmp_2) 1078 vec_rotated = rotate_vector(vec_2H_m - vec_2O, dihe_angle, rot_axis) 1079 ! set new position 1080 elem%pos(ind_arr(i) + H_offset*tmc_params%dim_per_elem: & 1081 ind_arr(i) + (H_offset + 1)*tmc_params%dim_per_elem - 1) & 1082 = vec_2O + vec_rotated 1083 vec_rotated = rotate_vector(vec_2H_f - vec_2O, dihe_angle, rot_axis) 1084 IF (H_offset .EQ. 1) THEN 1085 H_offset = 2 1086 ELSE 1087 H_offset = 1 1088 END IF 1089 elem%pos(ind_arr(i) + H_offset*tmc_params%dim_per_elem: & 1090 ind_arr(i) + (H_offset + 1)*tmc_params%dim_per_elem - 1) & 1091 = vec_2O + vec_rotated 1092 END IF 1093 END DO 1094 DEALLOCATE (tmp_cell) 1095 DEALLOCATE (ind_arr) 1096 ! end the timing 1097 CALL timestop(handle) 1098 END SUBROUTINE rotate_molecules_in_chain 1099 1100! ************************************************************************************************** 1101!> \brief volume move, the box size is increased or decreased, 1102!> using the mv_size a the factor. 1103!> the coordinated are scaled moleculewise 1104!> (the is moved like the center of mass is moves) 1105!> \param conf configuration to change with positions 1106!> \param T_ind temperature index, to select the correct temperature 1107!> for move size 1108!> \param move_types ... 1109!> \param rng_stream random number generator stream 1110!> \param tmc_params TMC parameters with e.g. dimensions of atoms and molecules 1111!> \param mv_cen_of_mass ... 1112!> \author Mandes 11.2012 1113! ************************************************************************************************** 1114 SUBROUTINE change_volume(conf, T_ind, move_types, rng_stream, tmc_params, & 1115 mv_cen_of_mass) 1116 TYPE(tree_type), POINTER :: conf 1117 INTEGER :: T_ind 1118 TYPE(tmc_move_type), POINTER :: move_types 1119 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 1120 TYPE(tmc_param_type), POINTER :: tmc_params 1121 LOGICAL :: mv_cen_of_mass 1122 1123 CHARACTER(LEN=*), PARAMETER :: routineN = 'change_volume', routineP = moduleN//':'//routineN 1124 1125 INTEGER :: atom, dir, handle, ind, ind_e, mol 1126 REAL(KIND=dp) :: rnd, vol 1127 REAL(KIND=dp), DIMENSION(3) :: box_length_new, box_length_orig, & 1128 box_scale_old 1129 REAL(KIND=dp), DIMENSION(:), POINTER :: disp, scaling 1130 1131 NULLIFY (scaling, disp) 1132 1133 CPASSERT(ASSOCIATED(conf)) 1134 CPASSERT(ASSOCIATED(move_types)) 1135 CPASSERT(ASSOCIATED(tmc_params)) 1136 CPASSERT(T_ind .GT. 0 .AND. T_ind .LE. tmc_params%nr_temp) 1137 CPASSERT(tmc_params%dim_per_elem .EQ. 3) 1138 CPASSERT(tmc_params%cell%orthorhombic) 1139 1140 ! start the timing 1141 CALL timeset(routineN, handle) 1142 1143 ALLOCATE (scaling(tmc_params%dim_per_elem)) 1144 ALLOCATE (disp(tmc_params%dim_per_elem)) 1145 1146 box_scale_old(:) = conf%box_scale 1147 ! get the cell vector length of the configuration (before move) 1148 CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf%box_scale, & 1149 abc=box_length_new) 1150 1151 IF (.FALSE.) THEN 1152 ! the volume move in volume space (dV) 1153 IF (tmc_params%v_isotropic) THEN 1154 CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf%box_scale, & 1155 abc=box_length_new, vol=vol) 1156 rnd = rng_stream%next() 1157 vol = vol + (rnd - 0.5_dp)*2.0_dp*move_types%mv_size(mv_type_volume_move, T_ind) 1158 box_length_new(:) = vol**(1/REAL(3, KIND=dp)) 1159 ELSE 1160 CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf%box_scale, & 1161 abc=box_length_new, vol=vol) 1162 rnd = rng_stream%next() 1163 vol = vol + (rnd - 0.5_dp)*2.0_dp*move_types%mv_size(mv_type_volume_move, T_ind) 1164 rnd = rng_stream%next() 1165 dir = 1 + INT(rnd*3) 1166 box_length_new(dir) = 1.0_dp 1167 box_length_new(dir) = vol/PRODUCT(box_length_new(:)) 1168 END IF 1169 ELSE 1170 ! the volume move in box length space (dL) 1171 ! increase / decrease box length in this direction 1172 ! l_n = l_o +- rnd * mv_size 1173 IF (tmc_params%v_isotropic) THEN 1174 rnd = rng_stream%next() 1175 box_length_new(:) = box_length_new(:) + & 1176 (rnd - 0.5_dp)*2.0_dp* & 1177 move_types%mv_size(mv_type_volume_move, T_ind) 1178 ELSE 1179 ! select a random direction 1180 rnd = rng_stream%next() 1181 dir = 1 + INT(rnd*3) 1182 rnd = rng_stream%next() 1183 box_length_new(dir) = box_length_new(dir) + & 1184 (rnd - 0.5_dp)*2.0_dp* & 1185 move_types%mv_size(mv_type_volume_move, T_ind) 1186 END IF 1187 END IF 1188 1189 ! get the original box length 1190 scaling(:) = 1.0_dp 1191 CALL get_scaled_cell(cell=tmc_params%cell, & 1192 box_scale=scaling, & 1193 abc=box_length_orig) 1194 ! get the new box scale 1195 conf%box_scale(:) = box_length_new(:)/box_length_orig(:) 1196 ! molecule scaling 1197 scaling(:) = conf%box_scale(:)/box_scale_old(:) 1198 1199 IF (mv_cen_of_mass .EQV. .FALSE.) THEN 1200 ! homogene scaling of atomic coordinates 1201 DO atom = 1, SIZE(conf%pos), tmc_params%dim_per_elem 1202 conf%pos(atom:atom + tmc_params%dim_per_elem - 1) = & 1203 conf%pos(atom:atom + tmc_params%dim_per_elem - 1)*scaling(:) 1204 END DO 1205 ELSE 1206 DO mol = 1, MAXVAL(conf%mol(:)) 1207 ! move the molecule related to the molecule center of mass 1208 ! get center of mass 1209 CPASSERT(ASSOCIATED(tmc_params%atoms)) 1210 1211 CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=conf%mol, mol=mol, & 1212 start_ind=ind, end_ind=ind_e) 1213 CALL center_of_mass( & 1214 pos=conf%pos(ind:ind_e + tmc_params%dim_per_elem - 1), & 1215 atoms=tmc_params%atoms(INT(ind/REAL(tmc_params%dim_per_elem, KIND=dp)) + 1: & 1216 INT(ind_e/REAL(tmc_params%dim_per_elem, KIND=dp)) + 1), & 1217 center=disp) 1218 ! calculate the center of mass DISPLACEMENT 1219 disp(:) = disp(:)*(scaling(:) - 1.0_dp) 1220 ! displace all atoms of the molecule 1221 DO atom = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem 1222 conf%pos(atom:atom + tmc_params%dim_per_elem - 1) = & 1223 conf%pos(atom:atom + tmc_params%dim_per_elem - 1) + disp(:) 1224 END DO 1225 END DO 1226 END IF 1227 1228 DEALLOCATE (scaling) 1229 DEALLOCATE (disp) 1230 1231 ! end the timing 1232 CALL timestop(handle) 1233 END SUBROUTINE change_volume 1234 1235! ************************************************************************************************** 1236!> \brief volume move, two atoms of different types are swapped, both selected 1237!> randomly 1238!> \param conf configuration to change with positions 1239!> \param move_types ... 1240!> \param rng_stream random number generator stream 1241!> \param tmc_params TMC parameters with e.g. dimensions of atoms and molecules 1242!> \author Mandes 11.2012 1243! ************************************************************************************************** 1244 SUBROUTINE swap_atoms(conf, move_types, rng_stream, tmc_params) 1245 TYPE(tree_type), POINTER :: conf 1246 TYPE(tmc_move_type), POINTER :: move_types 1247 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 1248 TYPE(tmc_param_type), POINTER :: tmc_params 1249 1250 CHARACTER(LEN=*), PARAMETER :: routineN = 'swap_atoms', routineP = moduleN//':'//routineN 1251 1252 INTEGER :: a_1, a_2, ind_1, ind_2 1253 LOGICAL :: found 1254 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pos_tmp 1255 1256 CPASSERT(ASSOCIATED(conf)) 1257 CPASSERT(ASSOCIATED(move_types)) 1258 CPASSERT(ASSOCIATED(tmc_params)) 1259 CPASSERT(ASSOCIATED(tmc_params%atoms)) 1260 1261 ! loop until two different atoms are found 1262 atom_search_loop: DO 1263 ! select one atom randomly 1264 a_1 = INT(SIZE(conf%pos)/REAL(tmc_params%dim_per_elem, KIND=dp)* & 1265 rng_stream%next()) + 1 1266 ! select the second atom randomly 1267 a_2 = INT(SIZE(conf%pos)/REAL(tmc_params%dim_per_elem, KIND=dp)* & 1268 rng_stream%next()) + 1 1269 ! check if they have different kinds 1270 IF (tmc_params%atoms(a_1)%name .NE. tmc_params%atoms(a_2)%name) THEN 1271 ! if present, check if atoms have different type related to the specified table 1272 IF (ASSOCIATED(move_types%atom_lists)) THEN 1273 DO ind_1 = 1, SIZE(move_types%atom_lists) 1274 IF (ANY(move_types%atom_lists(ind_1)%atoms(:) .EQ. & 1275 tmc_params%atoms(a_1)%name) .AND. & 1276 ANY(move_types%atom_lists(ind_1)%atoms(:) .EQ. & 1277 tmc_params%atoms(a_2)%name)) THEN 1278 found = .TRUE. 1279 EXIT atom_search_loop 1280 END IF 1281 END DO 1282 ELSE 1283 found = .TRUE. 1284 EXIT atom_search_loop 1285 END IF 1286 END IF 1287 END DO atom_search_loop 1288 IF (found) THEN 1289 ! perform coordinate exchange 1290 ALLOCATE (pos_tmp(tmc_params%dim_per_elem)) 1291 ind_1 = (a_1 - 1)*tmc_params%dim_per_elem + 1 1292 pos_tmp(:) = conf%pos(ind_1:ind_1 + tmc_params%dim_per_elem - 1) 1293 ind_2 = (a_2 - 1)*tmc_params%dim_per_elem + 1 1294 conf%pos(ind_1:ind_1 + tmc_params%dim_per_elem - 1) = & 1295 conf%pos(ind_2:ind_2 + tmc_params%dim_per_elem - 1) 1296 conf%pos(ind_2:ind_2 + tmc_params%dim_per_elem - 1) = pos_tmp(:) 1297 DEALLOCATE (pos_tmp) 1298 END IF 1299 END SUBROUTINE swap_atoms 1300 1301END MODULE tmc_moves 1302