(git:b2ae3e3)
Loading...
Searching...
No Matches
tmc_moves.F
Go to the documentation of this file.
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
15MODULE 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
54CONTAINS
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 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 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 > 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(:) == status_ok))
106 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 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 == 0) THEN
128 act_nr_elem_mv = SIZE(elem%pos)/tmc_params%dim_per_elem
129 END IF
130 ALLOCATE (elem_center(tmc_params%dim_per_elem))
131 i = 1
132 move_elements_loop: DO
133 ! select atom
134 IF (tmc_params%nr_elem_mv == 0) THEN
135 ind = (i - 1)*(tmc_params%dim_per_elem) + 1
136 ELSE
137 rnd = rng_stream%next()
138 ind = tmc_params%dim_per_elem* &
139 int(rnd*(SIZE(elem%pos)/tmc_params%dim_per_elem)) + 1
140 END IF
141 ! apply move
142 IF (elem%elem_stat(ind) == status_ok) THEN
143 ! displace atom
144 DO d = 0, tmc_params%dim_per_elem - 1
145 rnd = rng_stream%next()
146 elem%pos(ind + d) = elem%pos(ind + d) + (rnd - 0.5)*2.0* &
147 move_types%mv_size(mv_type_atom_trans, mv_conf)
148 END DO
149 ! check if new position is in subbox
150 elem_center = elem%pos(ind:ind + tmc_params%dim_per_elem - 1)
151 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 move_rejected = .true.
156 EXIT move_elements_loop
157 END IF
158 ELSE
159 ! element was not in sub box, search new one instead
160 IF (tmc_params%nr_elem_mv > 0) i = i - 1
161 END IF
162 i = i + 1
163 IF (i > act_nr_elem_mv) EXIT move_elements_loop
164 END DO move_elements_loop
165 DEALLOCATE (elem_center)
166
167 !-- molecule translation
168 CASE (mv_type_mol_trans)
169 nr_molec = maxval(elem%mol(:))
170 ! if all particles should be displaced, set the amount of molecules
171 IF (act_nr_elem_mv == 0) THEN
172 act_nr_elem_mv = nr_molec
173 END IF
174 ALLOCATE (mol_in_sb(nr_molec))
175 ALLOCATE (elem_center(tmc_params%dim_per_elem))
176 mol_in_sb(:) = status_frozen
177 ! check if any molecule is in sub_box
178 DO m = 1, nr_molec
179 CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
180 start_ind=ind, end_ind=ind_e)
181 CALL geometrical_center(pos=elem%pos(ind:ind_e + tmc_params%dim_per_elem - 1), &
182 center=elem_center)
183 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 ) THEN
187 mol_in_sb(m) = status_ok
188 END IF
189 END DO
190 ! displace the selected amount of molecules
191 IF (any(mol_in_sb(:) == status_ok)) THEN
192 ALLOCATE (direction(tmc_params%dim_per_elem))
193 counter = 1
194 move_molecule_loop: DO
195 ! select molecule
196 IF (tmc_params%nr_elem_mv == 0) THEN
197 m = counter
198 ELSE
199 rnd = rng_stream%next()
200 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 start_ind=ind, end_ind=ind_e)
204 ! when "molecule" is single atom, search a new one
205 IF (ind == ind_e) cycle move_molecule_loop
206
207 ! calculate displacement
208 ! move only molecules, with geom. center in subbox
209 IF (mol_in_sb(m) == status_ok) THEN
210 ! calculate displacement
211 DO d = 1, tmc_params%dim_per_elem
212 rnd = rng_stream%next()
213 direction(d) = (rnd - 0.5)*2.0_dp*move_types%mv_size( &
214 mv_type_mol_trans, mv_conf)
215 END DO
216 ! check if displaced position is still in subbox
217 elem_center(:) = elem_center(:) + direction(:)
218 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 atom_in_mol_loop: DO i = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem
224 dim_loop: DO d = 0, tmc_params%dim_per_elem - 1
225 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 move_rejected = .true.
231 EXIT move_molecule_loop
232 END IF
233 ELSE
234 ! element was not in sub box, search new one instead
235 IF (tmc_params%nr_elem_mv > 0) counter = counter - 1
236 END IF
237 counter = counter + 1
238 IF (counter > act_nr_elem_mv) EXIT move_molecule_loop
239 END DO move_molecule_loop
240 DEALLOCATE (direction)
241 END IF
242 DEALLOCATE (elem_center)
243 DEALLOCATE (mol_in_sb)
244
245 !-- molecule rotation
246 CASE (mv_type_mol_rot)
247 nr_molec = maxval(elem%mol(:))
248 IF (act_nr_elem_mv == 0) THEN
249 act_nr_elem_mv = nr_molec
250 END IF
251 ALLOCATE (mol_in_sb(nr_molec))
252 ALLOCATE (elem_center(tmc_params%dim_per_elem))
253 mol_in_sb(:) = status_frozen
254 ! check if any molecule is in sub_box
255 DO m = 1, nr_molec
256 CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
257 start_ind=ind, end_ind=ind_e)
258 CALL geometrical_center(pos=elem%pos(ind:ind_e + tmc_params%dim_per_elem - 1), &
259 center=elem_center)
260 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 ) THEN
264 mol_in_sb(m) = status_ok
265 END IF
266 END DO
267 ! rotate the selected amount of molecules
268 IF (any(mol_in_sb(:) == status_ok)) THEN
269 counter = 1
270 rot_molecule_loop: DO
271 ! select molecule
272 IF (tmc_params%nr_elem_mv == 0) THEN
273 m = counter
274 ELSE
275 rnd = rng_stream%next()
276 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 start_ind=ind, end_ind=ind_e)
280 ! when "molecule" is single atom, search a new one
281 IF (ind == ind_e) cycle rot_molecule_loop
282
283 ! apply move
284 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 dim_per_elem=tmc_params%dim_per_elem)
290 ! update sub box status of single atom
291 DO i = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem
292 elem_center = elem%pos(i:i + tmc_params%dim_per_elem - 1)
293 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 ) THEN
297 elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_ok
298 ELSE
299 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 IF (tmc_params%nr_elem_mv > 0) counter = counter - 1
305 END IF
306 counter = counter + 1
307 IF (counter > act_nr_elem_mv) EXIT rot_molecule_loop
308 END DO rot_molecule_loop
309 END IF
310 DEALLOCATE (elem_center)
311 DEALLOCATE (mol_in_sb)
312
313 !-- velocity changes for MD
314 !-- here all velocities are changed
315 CASE (mv_type_md)
316 cpassert(ASSOCIATED(tmc_params%atoms))
317 change_all_velocities_loop: DO i = 1, SIZE(elem%pos)
318 !-- attention, move type size is in atomic units of velocity
319 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 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 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 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 tmc_params=tmc_params)
349
350 CASE DEFAULT
351 CALL cp_abort(__location__, &
352 "unknown move type "// &
353 cp_to_string(elem%move_type))
354 END SELECT
355
356 CALL rng_stream%get(bg=elem%rng_seed(:, :, 1), &
357 cg=elem%rng_seed(:, :, 2), ig=elem%rng_seed(:, :, 3))
358
359 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 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 start_ind = -1
379 end_ind = -1
380
381 cpassert(ASSOCIATED(mol_arr))
382 cpassert(mol <= maxval(mol_arr(:)))
383 ! get start index
384 loop_start: DO i = 1, SIZE(mol_arr)
385 IF (mol_arr(i) == mol) THEN
386 start_ind = i
387 EXIT loop_start
388 END IF
389 END DO loop_start
390 ! get end index
391 loop_end: DO i = SIZE(mol_arr), i, -1
392 IF (mol_arr(i) == mol) THEN
393 end_ind = i
394 EXIT loop_end
395 END IF
396 END DO loop_end
397 ! check if all atoms inbetween attend to molecule
398 cpassert(all(mol_arr(start_ind:end_ind) == mol))
399 cpassert(start_ind > 0)
400 cpassert(end_ind > 0)
401 ! convert to indeces mapped for the position array (multiple dim per atom)
402 start_ind = (start_ind - 1)*tmc_params%dim_per_elem + 1
403 end_ind = (end_ind - 1)*tmc_params%dim_per_elem + 1
404 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: pos_tmp
427
428 cpassert(ASSOCIATED(pos))
429 cpassert(ASSOCIATED(subbox_center))
430 cpassert(ASSOCIATED(box_scale))
431 ! if pressure is defined, no scale should be 0
432 flag = .NOT. ((tmc_params%pressure > 0.0_dp) .AND. (any(box_scale == 0.0_dp)))
433 cpassert(flag)
434 cpassert(SIZE(pos) == 3)
435 cpassert(SIZE(pos) == SIZE(subbox_center))
436
437 ! start the timing
438 CALL timeset(routinen, handle)
439
440 ALLOCATE (pos_tmp(SIZE(pos)))
441
442 inside = .true.
443 ! return if no subbox is defined
444 IF (.NOT. any(tmc_params%sub_box_size(:) <= 0.1_dp)) THEN
445 pos_tmp(:) = pos(:) - subbox_center(:)
446 CALL get_scaled_cell(cell=tmc_params%cell, box_scale=box_scale, &
447 vec=pos_tmp)
448 ! check
449 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 inside = .false.
452 END IF
453 END IF
454 DEALLOCATE (pos_tmp)
455 ! end the timing
456 CALL timestop(handle)
457 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 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 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(:) <= 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 INTEGER :: i
553 REAL(kind=dp) :: a1, a2, a3, q0, q1, q2, q3, rnd
554 REAL(kind=dp), DIMENSION(3, 3) :: rot
555 REAL(kind=dp), DIMENSION(:), POINTER :: elem_center
556
557 NULLIFY (elem_center)
558
559 cpassert(ASSOCIATED(pos))
560 cpassert(dim_per_elem == 3)
561 cpassert(ind_start > 0 .AND. ind_start < SIZE(pos))
562 cpassert(ind_end > 0 .AND. ind_end < SIZE(pos))
563 cpassert(ASSOCIATED(move_types))
564 mark_used(move_types)
565
566 ! calculate rotation matrix (using quanternions)
567 rnd = rng_stream%next()
568 a1 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
569 rnd = rng_stream%next()
570 a2 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
571 rnd = rng_stream%next()
572 a3 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
573 q0 = cos(a2/2)*cos((a1 + a3)/2.0_dp)
574 q1 = sin(a2/2)*cos((a1 - a3)/2.0_dp)
575 q2 = sin(a2/2)*sin((a1 - a3)/2.0_dp)
576 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 2*(q1*q3 - q0*q2), 2*(q2*q3 + q0*q1), q0*q0 - q1*q1 - q2*q2 + q3*q3], [3, 3])
580
581 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 center=elem_center)
585
586 ! proceed rotation
587 atom_loop: DO i = ind_start, ind_end + dim_per_elem - 1, dim_per_elem
588 pos(i:i + 2) = matmul(pos(i:i + 2) - elem_center(:), rot) + elem_center(:)
589 END DO atom_loop
590 DEALLOCATE (elem_center)
591 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 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 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 rnd1 = rng_stream%next()
621 rnd2 = rng_stream%next()
622
623 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 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 rnd3 = rng_stream%next()
636 IF (rnd3 >= 0.5 .AND. rnd_sign_change) THEN
637 d = -1
638 ELSE
639 d = 1
640 END IF
641 vel = sin(phi)*delta_vel + cos(phi)*vel*d*1.0_dp
642 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 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 INTEGER, DIMENSION(:), POINTER :: mol_arr
670 REAL(kind=dp) :: rnd
671
672 NULLIFY (mol_arr)
673
674 cpassert(ASSOCIATED(elem))
675 cpassert(ASSOCIATED(tmc_params))
676
677 ! start the timing
678 CALL timeset(routinen, handle)
679
680 short_loop = .false.
681 counter = 0
682 nr_mol = maxval(elem%mol(:))
683 ! ind_arr: one array element for each molecule
684 ALLOCATE (mol_arr(nr_mol))
685 mol_arr(:) = -1
686 donor_acceptor = not_selected
687 ! select randomly if neighboring molecule is donor / acceptor
688 IF (rng_stream%next() < 0.5_dp) THEN
689 donor_acceptor = proton_acceptor
690 ELSE
691 donor_acceptor = proton_donor
692 END IF
693
694 ! first step build loop
695 ! select randomly one atom
696 rnd = rng_stream%next()
697 ! the randomly selected first atom
698 mol = int(rnd*nr_mol) + 1
699 counter = counter + 1
700 mol_arr(counter) = mol
701
702 ! do until the loop is closed
703 ! (until path connects back to any spot of the path)
704 chain_completition_loop: DO
705 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 rng_stream=rng_stream)
711 IF (any(mol_arr(:) == mol)) THEN
712 EXIT chain_completition_loop
713 END IF
714 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 DO k = 1, counter
720 IF (mol_arr(k) == mol) THEN
721 EXIT
722 END IF
723 END DO
724 mol_arr(1:counter - k + 1) = mol_arr(k:counter)
725 counter = counter - k + 1
726
727 ! check if loop is minimum size of 6 molecules
728 IF (counter < 6) THEN
729 CALL cp_warn(__location__, &
730 "short proton loop with"//cp_to_string(counter)// &
731 "molecules.")
732 tmp_chr = ""
733 WRITE (tmp_chr, *) mol_arr(1:counter)
734 cpwarn("selected molecules:"//trim(tmp_chr))
735 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 mol_arr_in=mol_arr(1:counter), donor_acceptor=donor_acceptor)
742 DEALLOCATE (mol_arr)
743
744 ! end the timing
745 CALL timestop(handle)
746 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: disth1, disth2, disto
774
775 NULLIFY (disto, disth1, disth2)
776 cpassert(ASSOCIATED(elem))
777 cpassert(ASSOCIATED(tmc_params))
778
779 ! start the timing
780 CALL timeset(routinen, handle)
781
782 nr_mol = maxval(elem%mol)
783 ALLOCATE (disto(nr_mol))
784 ALLOCATE (disth1(nr_mol))
785 ALLOCATE (disth2(nr_mol))
786 !-- initialize the distances to huge values
787 ! distance of nearest proton of certain molecule to preselected O
788 disto(:) = huge(disto(1))
789 ! distance of (first) proton of preselected molecule to certain molecule
790 disth1(:) = huge(disth1(1))
791 ! distance of (second) proton of preselected molecule to certain molecule
792 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 start_ind=ind, end_ind=ind_e)
797
798 ! calculate distances to all molecules
799 list_distances: DO mol_tmp = 1, nr_mol
800 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 mol=mol_tmp, start_ind=ind_n, end_ind=ind_e)
805 ! check if selected molecule is water respectively consists of 3 atoms
806 IF (mod(ind_e - ind_n, 3) > 0) THEN
807 CALL cp_warn(__location__, &
808 "selected a molecule with more than 3 atoms, "// &
809 "the proton reordering does not support, skip molecule")
810 cycle list_distances
811 END IF
812 IF (donor_acceptor == proton_acceptor) THEN
813 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 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 cell=tmc_params%cell, box_scale=elem%box_scale)
827 END IF
828 END IF
829 !check for neighboring proton donors
830 IF (donor_acceptor == proton_donor) THEN
831 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 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 cell=tmc_params%cell, box_scale=elem%box_scale)
844 IF (dist_tmp < disto(mol_tmp)) disto(mol_tmp) = dist_tmp
845 END IF
846 END IF
847 END DO list_distances
848
849 mol_tmp = 1
850 ! select the nearest neighbors
851 !check for neighboring proton acceptors
852 IF (donor_acceptor == proton_acceptor) THEN
853 neighbor_mol(mol_tmp) = minloc(disth1(:), 1)
854 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 IF (neighbor_mol(mol_tmp) == neighbor_mol(mol_tmp + 1)) THEN
857 disth1(neighbor_mol(mol_tmp)) = huge(disth1(1))
858 disth2(neighbor_mol(mol_tmp + 1)) = huge(disth2(1))
859 IF (minval(disth1(:), 1) < minval(disth2(:), 1)) THEN
860 neighbor_mol(mol_tmp) = minloc(disth1(:), 1)
861 ELSE
862 neighbor_mol(mol_tmp + 1) = minloc(disth2(:), 1)
863 END IF
864 END IF
865 mol_tmp = mol_tmp + 2
866 END IF
867
868 !check for neighboring proton donors
869 IF (donor_acceptor == proton_donor) THEN
870 neighbor_mol(mol_tmp) = minloc(disto(:), 1)
871 disto(neighbor_mol(mol_tmp)) = huge(disto(1))
872 neighbor_mol(mol_tmp + 1) = minloc(disto(:), 1)
873 END IF
874
875 ! select randomly the next neighboring molecule
876 rnd = rng_stream%next()
877 ! the randomly selected atom: return value!
878 mol_tmp = neighbor_mol(int(rnd*SIZE(neighbor_mol(:))) + 1)
879 mol = mol_tmp
880
881 DEALLOCATE (disto)
882 DEALLOCATE (disth1)
883 DEALLOCATE (disth2)
884
885 ! end the timing
886 CALL timestop(handle)
887 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 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 cpassert(ASSOCIATED(elem))
909 cpassert(i_orig >= 1 .AND. i_orig <= SIZE(elem%pos))
910 cpassert(i_neighbor >= 1 .AND. i_neighbor <= SIZE(elem%pos))
911 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 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 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 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 cell=tmc_params%cell, box_scale=elem%box_scale)
937
938 IF (minloc(distances(:), 1) <= 2) THEN
939 donor_acceptor = proton_acceptor
940 ELSE
941 donor_acceptor = proton_donor
942 END IF
943 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 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 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 NULLIFY (ind_arr, tmp_cell)
973
974 cpassert(ASSOCIATED(tmc_params))
975 cpassert(ASSOCIATED(elem))
976
977 ! start the timing
978 CALL timeset(routinen, handle)
979
980 ALLOCATE (ind_arr(0:SIZE(mol_arr_in) + 1))
981 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 start_ind=ind_arr(i), end_ind=ind)
985 END DO
986 ind_arr(0) = ind_arr(SIZE(ind_arr) - 2)
987 ind_arr(SIZE(ind_arr) - 1) = ind_arr(1)
988
989 ! get the scaled cell
990 ALLOCATE (tmp_cell)
991 CALL get_scaled_cell(cell=tmc_params%cell, box_scale=elem%box_scale, &
992 scaled_cell=tmp_cell)
993
994 ! rotate single molecules
995 DO i = 1, SIZE(ind_arr) - 2
996 ! the 3 O atoms
997 vec_1o(:) = elem%pos(ind_arr(i - 1):ind_arr(i - 1) + tmc_params%dim_per_elem - 1)
998 vec_2o(:) = elem%pos(ind_arr(i):ind_arr(i) + tmc_params%dim_per_elem - 1)
999 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 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 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 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 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 ind_arr(i) + 3*tmc_params%dim_per_elem - 1)
1028 h_offset = 2
1029 END IF
1030
1031 IF (.true.) THEN !TODO find a better switch for the pauling model
1032
1033 ! do rotation (NOT pauling model)
1034 tmp_1 = pbc(vec_2o - vec_1o, tmp_cell)
1035 tmp_2 = pbc(vec_3o - vec_2h_f, tmp_cell)
1036
1037 dihe_angle = donor_acceptor*dihedral_angle(tmp_1, vec_2h_f - vec_2o, tmp_2)
1038 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 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 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 DEALLOCATE (tmp_cell)
1088 DEALLOCATE (ind_arr)
1089 ! end the timing
1090 CALL timestop(handle)
1091 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: disp, scaling
1123
1124 NULLIFY (scaling, disp)
1125
1126 cpassert(ASSOCIATED(conf))
1127 cpassert(ASSOCIATED(move_types))
1128 cpassert(ASSOCIATED(tmc_params))
1129 cpassert(t_ind > 0 .AND. t_ind <= tmc_params%nr_temp)
1130 cpassert(tmc_params%dim_per_elem == 3)
1131 cpassert(tmc_params%cell%orthorhombic)
1132
1133 ! start the timing
1134 CALL timeset(routinen, handle)
1135
1136 ALLOCATE (scaling(tmc_params%dim_per_elem))
1137 ALLOCATE (disp(tmc_params%dim_per_elem))
1138
1139 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 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 IF (tmc_params%v_isotropic) THEN
1167 rnd = rng_stream%next()
1168 box_length_new(:) = box_length_new(:) + &
1169 (rnd - 0.5_dp)*2.0_dp* &
1170 move_types%mv_size(mv_type_volume_move, t_ind)
1171 ELSE
1172 ! select a random direction
1173 rnd = rng_stream%next()
1174 dir = 1 + int(rnd*3)
1175 rnd = rng_stream%next()
1176 box_length_new(dir) = box_length_new(dir) + &
1177 (rnd - 0.5_dp)*2.0_dp* &
1178 move_types%mv_size(mv_type_volume_move, t_ind)
1179 END IF
1180 END IF
1181
1182 ! get the original box length
1183 scaling(:) = 1.0_dp
1184 CALL get_scaled_cell(cell=tmc_params%cell, &
1185 box_scale=scaling, &
1186 abc=box_length_orig)
1187 ! get the new box scale
1188 conf%box_scale(:) = box_length_new(:)/box_length_orig(:)
1189 ! molecule scaling
1190 scaling(:) = conf%box_scale(:)/box_scale_old(:)
1191
1192 IF (mv_cen_of_mass .EQV. .false.) THEN
1193 ! homogene scaling of atomic coordinates
1194 DO atom = 1, SIZE(conf%pos), tmc_params%dim_per_elem
1195 conf%pos(atom:atom + tmc_params%dim_per_elem - 1) = &
1196 conf%pos(atom:atom + tmc_params%dim_per_elem - 1)*scaling(:)
1197 END DO
1198 ELSE
1199 DO mol = 1, maxval(conf%mol(:))
1200 ! move the molecule related to the molecule center of mass
1201 ! get center of mass
1202 cpassert(ASSOCIATED(tmc_params%atoms))
1203
1204 CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=conf%mol, mol=mol, &
1205 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 center=disp)
1211 ! calculate the center of mass DISPLACEMENT
1212 disp(:) = disp(:)*(scaling(:) - 1.0_dp)
1213 ! displace all atoms of the molecule
1214 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 conf%pos(atom:atom + tmc_params%dim_per_elem - 1) + disp(:)
1217 END DO
1218 END DO
1219 END IF
1220
1221 DEALLOCATE (scaling)
1222 DEALLOCATE (disp)
1223
1224 ! end the timing
1225 CALL timestop(handle)
1226 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: pos_tmp
1246
1247 cpassert(ASSOCIATED(conf))
1248 cpassert(ASSOCIATED(move_types))
1249 cpassert(ASSOCIATED(tmc_params))
1250 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 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 rng_stream%next()) + 1
1260 ! check if they have different kinds
1261 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 IF (ASSOCIATED(move_types%atom_lists)) THEN
1264 DO ind_1 = 1, SIZE(move_types%atom_lists)
1265 IF (any(move_types%atom_lists(ind_1)%atoms(:) == &
1266 tmc_params%atoms(a_1)%name) .AND. &
1267 any(move_types%atom_lists(ind_1)%atoms(:) == &
1268 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 ALLOCATE (pos_tmp(tmc_params%dim_per_elem))
1282 ind_1 = (a_1 - 1)*tmc_params%dim_per_elem + 1
1283 pos_tmp(:) = conf%pos(ind_1:ind_1 + tmc_params%dim_per_elem - 1)
1284 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 conf%pos(ind_2:ind_2 + tmc_params%dim_per_elem - 1)
1287 conf%pos(ind_2:ind_2 + tmc_params%dim_per_elem - 1) = pos_tmp(:)
1288 DEALLOCATE (pos_tmp)
1289 END IF
1290 END SUBROUTINE swap_atoms
1291
1292END MODULE tmc_moves
cell_types::pbc
Definition cell_types.F:103
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
atom
Definition atom.F:9
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:210
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::rotate_vector
pure real(kind=dp) function, dimension(3), public rotate_vector(a, phi, b)
Rotation of the vector a about an rotation axis defined by the vector b. The rotation angle is phi (r...
Definition mathlib.F:1133
mathlib::dihedral_angle
pure real(kind=dp) function, public dihedral_angle(ab, bc, cd)
Returns the dihedral angle, i.e. the angle between the planes defined by the vectors (-ab,...
Definition mathlib.F:476
parallel_rng_types
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition parallel_rng_types.F:51
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::boltzmann
real(kind=dp), parameter, public boltzmann
Definition physcon.F:129
physcon::joule
real(kind=dp), parameter, public joule
Definition physcon.F:159
tmc_calculations
calculation section for TreeMonteCarlo
Definition tmc_calculations.F:15
tmc_calculations::geometrical_center
subroutine, public geometrical_center(pos, center)
calculate the geometrical center of an amount of atoms array size should be multiple of dim_per_elem
Definition tmc_calculations.F:269
tmc_calculations::get_scaled_cell
subroutine, public get_scaled_cell(cell, box_scale, scaled_hmat, scaled_cell, vol, abc, vec)
handles properties and calculations of a scaled cell
Definition tmc_calculations.F:153
tmc_calculations::center_of_mass
subroutine, public center_of_mass(pos, atoms, center)
calculate the center of mass of an amount of atoms array size should be multiple of dim_per_elem
Definition tmc_calculations.F:302
tmc_calculations::nearest_distance
real(kind=dp) function, public nearest_distance(x1, x2, cell, box_scale)
neares distance of atoms within the periodic boundary condition
Definition tmc_calculations.F:234
tmc_move_types
tree nodes creation, searching, deallocation, references etc.
Definition tmc_move_types.F:15
tmc_move_types::mv_type_mol_rot
integer, parameter, public mv_type_mol_rot
Definition tmc_move_types.F:31
tmc_move_types::mv_type_volume_move
integer, parameter, public mv_type_volume_move
Definition tmc_move_types.F:35
tmc_move_types::mv_type_proton_reorder
integer, parameter, public mv_type_proton_reorder
Definition tmc_move_types.F:32
tmc_move_types::mv_type_md
integer, parameter, public mv_type_md
Definition tmc_move_types.F:34
tmc_move_types::mv_type_mol_trans
integer, parameter, public mv_type_mol_trans
Definition tmc_move_types.F:30
tmc_move_types::mv_type_atom_swap
integer, parameter, public mv_type_atom_swap
Definition tmc_move_types.F:33
tmc_move_types::mv_type_gausian_adapt
integer, parameter, public mv_type_gausian_adapt
Definition tmc_move_types.F:36
tmc_move_types::mv_type_atom_trans
integer, parameter, public mv_type_atom_trans
Definition tmc_move_types.F:29
tmc_moves
different move types are applied
Definition tmc_moves.F:15
tmc_moves::elements_in_new_subbox
subroutine, public elements_in_new_subbox(tmc_params, rng_stream, elem, nr_of_sub_box_elements)
set a new random sub box center and counte the number of atoms in it
Definition tmc_moves.F:471
tmc_moves::change_pos
subroutine, public change_pos(tmc_params, move_types, rng_stream, elem, mv_conf, new_subbox, move_rejected)
applying the preselected move type
Definition tmc_moves.F:70
tmc_tree_types
module handles definition of the tree nodes for the global and the subtrees binary tree parent elemen...
Definition tmc_tree_types.F:32
tmc_tree_types::status_ok
integer, parameter, public status_ok
Definition tmc_tree_types.F:70
tmc_tree_types::status_frozen
integer, parameter, public status_frozen
Definition tmc_tree_types.F:71
tmc_types
module handles definition of the tree nodes for the global and the subtrees binary tree parent elemen...
Definition tmc_types.F:32
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157
tmc_move_types::tmc_move_type
Definition tmc_move_types.F:42
tmc_tree_types::tree_type
Definition tmc_tree_types.F:75
tmc_types::tmc_atom_type
Definition tmc_types.F:124
tmc_types::tmc_param_type
Definition tmc_types.F:90