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