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