Line data Source code
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 :
15 : MODULE tmc_moves
16 : USE cell_types, ONLY: cell_type,&
17 : get_cell,&
18 : pbc
19 : USE cp_log_handling, ONLY: cp_to_string
20 : USE kinds, ONLY: dp
21 : USE mathconstants, ONLY: pi
22 : USE mathlib, ONLY: dihedral_angle,&
23 : rotate_vector
24 : USE parallel_rng_types, ONLY: rng_stream_type
25 : USE physcon, ONLY: boltzmann,&
26 : joule
27 : USE tmc_calculations, ONLY: center_of_mass,&
28 : geometrical_center,&
29 : get_scaled_cell,&
30 : nearest_distance
31 : USE tmc_move_types, ONLY: &
32 : mv_type_MD, mv_type_atom_swap, mv_type_atom_trans, mv_type_gausian_adapt, mv_type_mol_rot, &
33 : mv_type_mol_trans, mv_type_proton_reorder, mv_type_volume_move, tmc_move_type
34 : USE tmc_tree_types, ONLY: status_frozen,&
35 : status_ok,&
36 : tree_type
37 : USE tmc_types, ONLY: tmc_atom_type,&
38 : tmc_param_type
39 : #include "../base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_moves'
46 :
47 : PUBLIC :: change_pos
48 : PUBLIC :: elements_in_new_subbox
49 :
50 : INTEGER, PARAMETER :: not_selected = 0
51 : INTEGER, PARAMETER :: proton_donor = -1
52 : INTEGER, PARAMETER :: proton_acceptor = 1
53 :
54 : CONTAINS
55 : ! **************************************************************************************************
56 : !> \brief applying the preselected move type
57 : !> \param tmc_params TMC parameters with dimensions ...
58 : !> \param move_types ...
59 : !> \param rng_stream random number stream
60 : !> \param elem configuration to change
61 : !> \param mv_conf temperature index for determinig the move size
62 : !> \param new_subbox flag if new sub box should be crated
63 : !> \param move_rejected return flag if during configurational change
64 : !> configuration should still be accepted (not if e.g. atom/molecule
65 : !> leave the sub box
66 : !> \author Mandes 12.2012
67 : ! **************************************************************************************************
68 4400 : 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 4400 : INTEGER, DIMENSION(:), POINTER :: mol_in_sb
80 : REAL(KIND=dp) :: rnd
81 4400 : REAL(KIND=dp), DIMENSION(:), POINTER :: direction, elem_center
82 :
83 4400 : NULLIFY (direction, elem_center, mol_in_sb)
84 :
85 0 : CPASSERT(ASSOCIATED(tmc_params))
86 4400 : CPASSERT(ASSOCIATED(move_types))
87 4400 : CPASSERT(ASSOCIATED(elem))
88 :
89 4400 : move_rejected = .FALSE.
90 :
91 : CALL rng_stream%set(bg=elem%rng_seed(:, :, 1), &
92 4400 : cg=elem%rng_seed(:, :, 2), ig=elem%rng_seed(:, :, 3))
93 :
94 4400 : IF (new_subbox) THEN
95 4093 : 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 0 : nr_of_sub_box_elements=nr_sub_box_elem)
99 : ELSE
100 266446 : elem%elem_stat(:) = status_ok
101 : END IF
102 : END IF
103 :
104 : ! at least one atom should be in the sub box
105 17738 : CPASSERT(ANY(elem%elem_stat(:) .EQ. status_ok))
106 4400 : 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 4400 : SELECT CASE (elem%move_type)
118 : CASE (mv_type_gausian_adapt)
119 : ! just for Gaussian Adaptation
120 0 : CPABORT("gaussian adaptation is not imlemented yet.")
121 : !TODO CALL new_pos_gauss_adapt(acc=ASSOCIATED(elem%parent%acc, elem), &
122 : ! pos=elem%pos, covari=elem%frc, pot=elem%potential, &
123 : ! step_size=elem%ekin, pos_aver=elem%vel, temp=elem%ekin_before_md, &
124 : ! rng_seed=elem%rng_seed, rng_seed_last_acc=last_acc_elem%rng_seed)
125 : !-- atom translation
126 : CASE (mv_type_atom_trans)
127 3511 : IF (act_nr_elem_mv .EQ. 0) &
128 264 : act_nr_elem_mv = SIZE(elem%pos)/tmc_params%dim_per_elem
129 10533 : ALLOCATE (elem_center(tmc_params%dim_per_elem))
130 3511 : i = 1
131 : move_elements_loop: DO
132 : ! select atom
133 12396 : IF (tmc_params%nr_elem_mv .EQ. 0) THEN
134 9149 : ind = (i - 1)*(tmc_params%dim_per_elem) + 1
135 : ELSE
136 3247 : rnd = rng_stream%next()
137 : ind = tmc_params%dim_per_elem* &
138 3247 : INT(rnd*(SIZE(elem%pos)/tmc_params%dim_per_elem)) + 1
139 : END IF
140 : ! apply move
141 12396 : IF (elem%elem_stat(ind) .EQ. status_ok) THEN
142 : ! displace atom
143 33088 : DO d = 0, tmc_params%dim_per_elem - 1
144 24816 : rnd = rng_stream%next()
145 : elem%pos(ind + d) = elem%pos(ind + d) + (rnd - 0.5)*2.0* &
146 33088 : move_types%mv_size(mv_type_atom_trans, mv_conf)
147 : END DO
148 : ! check if new position is in subbox
149 66176 : elem_center = elem%pos(ind:ind + tmc_params%dim_per_elem - 1)
150 8272 : 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 6 : move_rejected = .TRUE.
155 6 : EXIT move_elements_loop
156 : END IF
157 : ELSE
158 : ! element was not in sub box, search new one instead
159 4124 : IF (tmc_params%nr_elem_mv .GT. 0) i = i - 1
160 : END IF
161 12390 : i = i + 1
162 12390 : IF (i .GT. act_nr_elem_mv) EXIT move_elements_loop
163 : END DO move_elements_loop
164 3511 : DEALLOCATE (elem_center)
165 :
166 : !-- molecule translation
167 : CASE (mv_type_mol_trans)
168 8052 : nr_molec = MAXVAL(elem%mol(:))
169 : ! if all particles should be displaced, set the amount of molecules
170 207 : IF (act_nr_elem_mv .EQ. 0) &
171 201 : act_nr_elem_mv = nr_molec
172 621 : ALLOCATE (mol_in_sb(nr_molec))
173 621 : ALLOCATE (elem_center(tmc_params%dim_per_elem))
174 2812 : mol_in_sb(:) = status_frozen
175 : ! check if any molecule is in sub_box
176 2812 : DO m = 1, nr_molec
177 : CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
178 2605 : start_ind=ind, end_ind=ind_e)
179 : CALL geometrical_center(pos=elem%pos(ind:ind_e + tmc_params%dim_per_elem - 1), &
180 2605 : center=elem_center)
181 2605 : 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 4680 : mol_in_sb(m) = status_ok
186 : END DO
187 : ! displace the selected amount of molecules
188 558 : IF (ANY(mol_in_sb(:) .EQ. status_ok)) THEN
189 621 : ALLOCATE (direction(tmc_params%dim_per_elem))
190 1638 : counter = 1
191 1638 : move_molecule_loop: DO
192 : ! select molecule
193 1638 : IF (tmc_params%nr_elem_mv .EQ. 0) THEN
194 1632 : m = counter
195 : ELSE
196 6 : rnd = rng_stream%next()
197 6 : 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 1638 : start_ind=ind, end_ind=ind_e)
201 : ! when "molecule" is single atom, search a new one
202 1638 : IF (ind .EQ. ind_e) CYCLE move_molecule_loop
203 :
204 : ! calculate displacement
205 : ! move only molecules, with geom. center in subbox
206 1638 : IF (mol_in_sb(m) .EQ. status_ok) THEN
207 : ! calculate displacement
208 5124 : DO d = 1, tmc_params%dim_per_elem
209 3843 : rnd = rng_stream%next()
210 : direction(d) = (rnd - 0.5)*2.0_dp*move_types%mv_size( &
211 5124 : mv_type_mol_trans, mv_conf)
212 : END DO
213 : ! check if displaced position is still in subbox
214 10248 : elem_center(:) = elem_center(:) + direction(:)
215 1281 : 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 5134 : atom_in_mol_loop: DO i = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem
221 16708 : dim_loop: DO d = 0, tmc_params%dim_per_elem - 1
222 15432 : 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 5 : move_rejected = .TRUE.
228 5 : EXIT move_molecule_loop
229 : END IF
230 : ELSE
231 : ! element was not in sub box, search new one instead
232 357 : IF (tmc_params%nr_elem_mv .GT. 0) counter = counter - 1
233 : END IF
234 1633 : counter = counter + 1
235 1633 : IF (counter .GT. act_nr_elem_mv) EXIT move_molecule_loop
236 : END DO move_molecule_loop
237 207 : DEALLOCATE (direction)
238 : END IF
239 207 : DEALLOCATE (elem_center)
240 207 : DEALLOCATE (mol_in_sb)
241 :
242 : !-- molecule rotation
243 : CASE (mv_type_mol_rot)
244 10810 : nr_molec = MAXVAL(elem%mol(:))
245 261 : IF (act_nr_elem_mv .EQ. 0) &
246 257 : act_nr_elem_mv = nr_molec
247 783 : ALLOCATE (mol_in_sb(nr_molec))
248 783 : ALLOCATE (elem_center(tmc_params%dim_per_elem))
249 3766 : mol_in_sb(:) = status_frozen
250 : ! check if any molecule is in sub_box
251 3766 : DO m = 1, nr_molec
252 : CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
253 3505 : start_ind=ind, end_ind=ind_e)
254 : CALL geometrical_center(pos=elem%pos(ind:ind_e + tmc_params%dim_per_elem - 1), &
255 3505 : center=elem_center)
256 3505 : 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 5796 : mol_in_sb(m) = status_ok
261 : END DO
262 : ! rotate the selected amount of molecules
263 961 : IF (ANY(mol_in_sb(:) .EQ. status_ok)) THEN
264 : counter = 1
265 3125 : rot_molecule_loop: DO
266 : ! select molecule
267 3125 : IF (tmc_params%nr_elem_mv .EQ. 0) THEN
268 3121 : m = counter
269 : ELSE
270 4 : rnd = rng_stream%next()
271 4 : 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 3125 : start_ind=ind, end_ind=ind_e)
275 : ! when "molecule" is single atom, search a new one
276 3125 : IF (ind .EQ. ind_e) CYCLE rot_molecule_loop
277 :
278 : ! apply move
279 3125 : 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 1650 : dim_per_elem=tmc_params%dim_per_elem)
285 : ! update sub box status of single atom
286 6634 : DO i = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem
287 39872 : elem_center = elem%pos(i:i + tmc_params%dim_per_elem - 1)
288 4984 : 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 1650 : ) THEN
292 19772 : elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_ok
293 : ELSE
294 164 : 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 1475 : IF (tmc_params%nr_elem_mv .GT. 0) counter = counter - 1
300 : END IF
301 3125 : counter = counter + 1
302 3125 : IF (counter .GT. act_nr_elem_mv) EXIT rot_molecule_loop
303 : END DO rot_molecule_loop
304 : END IF
305 261 : DEALLOCATE (elem_center)
306 261 : DEALLOCATE (mol_in_sb)
307 :
308 : !-- velocity changes for MD
309 : !-- here all velocities are changed
310 : CASE (mv_type_MD)
311 0 : CPASSERT(ASSOCIATED(tmc_params%atoms))
312 0 : change_all_velocities_loop: DO i = 1, SIZE(elem%pos)
313 : !-- attention, move type size is in atomic units of velocity
314 0 : 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 0 : 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 12 : 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 224 : 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 185 : tmc_params=tmc_params)
344 :
345 : CASE DEFAULT
346 : CALL cp_abort(__LOCATION__, &
347 : "unknown move type "// &
348 4400 : cp_to_string(elem%move_type))
349 : END SELECT
350 :
351 : CALL rng_stream%get(bg=elem%rng_seed(:, :, 1), &
352 4400 : cg=elem%rng_seed(:, :, 2), ig=elem%rng_seed(:, :, 3))
353 :
354 4400 : 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 27693 : 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 27693 : start_ind = -1
374 27693 : end_ind = -1
375 :
376 27693 : CPASSERT(ASSOCIATED(mol_arr))
377 6479417 : CPASSERT(mol .LE. MAXVAL(mol_arr(:)))
378 : ! get start index
379 3201629 : loop_start: DO i = 1, SIZE(mol_arr)
380 3201629 : IF (mol_arr(i) .EQ. mol) THEN
381 27693 : start_ind = i
382 27693 : EXIT loop_start
383 : END IF
384 : END DO loop_start
385 : ! get end index
386 3222274 : loop_end: DO i = SIZE(mol_arr), i, -1
387 3222274 : IF (mol_arr(i) .EQ. mol) THEN
388 27693 : end_ind = i
389 27693 : EXIT loop_end
390 : END IF
391 : END DO loop_end
392 : ! check if all atoms inbetween attend to molecule
393 110900 : CPASSERT(ALL(mol_arr(start_ind:end_ind) .EQ. mol))
394 27693 : CPASSERT(start_ind .GT. 0)
395 27693 : CPASSERT(end_ind .GT. 0)
396 : ! convert to indeces mapped for the position array (multiple dim per atom)
397 27693 : start_ind = (start_ind - 1)*tmc_params%dim_per_elem + 1
398 27693 : end_ind = (end_ind - 1)*tmc_params%dim_per_elem + 1
399 27693 : 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 24679 : 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 24679 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pos_tmp
422 :
423 24679 : CPASSERT(ASSOCIATED(pos))
424 24679 : CPASSERT(ASSOCIATED(subbox_center))
425 24679 : CPASSERT(ASSOCIATED(box_scale))
426 : ! if pressure is defined, no scale should be 0
427 98716 : flag = .NOT. ((tmc_params%pressure .GT. 0.0_dp) .AND. (ANY(box_scale .EQ. 0.0_dp)))
428 0 : CPASSERT(flag)
429 24679 : CPASSERT(SIZE(pos) .EQ. 3)
430 24679 : CPASSERT(SIZE(pos) .EQ. SIZE(subbox_center))
431 :
432 : ! start the timing
433 24679 : CALL timeset(routineN, handle)
434 :
435 74037 : ALLOCATE (pos_tmp(SIZE(pos)))
436 :
437 24679 : inside = .TRUE.
438 : ! return if no subbox is defined
439 44695 : IF (.NOT. ANY(tmc_params%sub_box_size(:) .LE. 0.1_dp)) THEN
440 26688 : pos_tmp(:) = pos(:) - subbox_center(:)
441 : CALL get_scaled_cell(cell=tmc_params%cell, box_scale=box_scale, &
442 6672 : vec=pos_tmp)
443 : ! check
444 33520 : 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 6142 : inside = .FALSE.
447 : END IF
448 : END IF
449 24679 : DEALLOCATE (pos_tmp)
450 : ! end the timing
451 24679 : CALL timestop(handle)
452 24679 : 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 114 : 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 57 : REAL(KIND=dp), DIMENSION(:), POINTER :: atom_tmp, center_of_sub_box
477 :
478 57 : NULLIFY (center_of_sub_box, atom_tmp)
479 :
480 0 : CPASSERT(ASSOCIATED(tmc_params))
481 57 : CPASSERT(ASSOCIATED(elem))
482 :
483 : ! start the timing
484 57 : CALL timeset(routineN, handle)
485 :
486 99 : 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 37195 : elem%elem_stat = status_ok
489 43 : nr_of_sub_box_elements = SIZE(elem%elem_stat)
490 : ELSE
491 42 : ALLOCATE (center_of_sub_box(tmc_params%dim_per_elem))
492 42 : ALLOCATE (atom_tmp(tmc_params%dim_per_elem))
493 14 : 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 14 : ig=elem%rng_seed(:, :, 3))
497 :
498 14 : CALL get_cell(cell=tmc_params%cell, abc=box_size)
499 56 : DO i = 1, SIZE(tmc_params%sub_box_size)
500 42 : rnd = rng_stream%next()
501 56 : center_of_sub_box(i) = rnd*box_size(i)
502 : END DO
503 112 : elem%subbox_center(:) = center_of_sub_box(:)
504 :
505 : CALL rng_stream%get(bg=elem%rng_seed(:, :, 1), cg=elem%rng_seed(:, :, 2), &
506 14 : ig=elem%rng_seed(:, :, 3))
507 :
508 : ! check all elements if they are in subbox
509 4046 : DO i = 1, SIZE(elem%pos), tmc_params%dim_per_elem
510 32256 : atom_tmp(:) = elem%pos(i:i + tmc_params%dim_per_elem - 1)
511 4032 : IF (check_pos_in_subbox(pos=atom_tmp, &
512 : subbox_center=center_of_sub_box, box_scale=elem%box_scale, &
513 14 : tmc_params=tmc_params)) THEN
514 616 : elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_ok
515 154 : nr_of_sub_box_elements = nr_of_sub_box_elements + 1
516 : ELSE
517 15512 : elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_frozen
518 : END IF
519 : END DO
520 14 : DEALLOCATE (atom_tmp)
521 14 : DEALLOCATE (center_of_sub_box)
522 : END IF
523 : ! end the timing
524 57 : CALL timestop(handle)
525 57 : 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 1650 : 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 1650 : REAL(KIND=dp), DIMENSION(:), POINTER :: elem_center
551 :
552 1650 : NULLIFY (elem_center)
553 :
554 0 : CPASSERT(ASSOCIATED(pos))
555 1650 : CPASSERT(dim_per_elem .EQ. 3)
556 1650 : CPASSERT(ind_start .GT. 0 .AND. ind_start .LT. SIZE(pos))
557 1650 : CPASSERT(ind_end .GT. 0 .AND. ind_end .LT. SIZE(pos))
558 1650 : CPASSERT(ASSOCIATED(move_types))
559 : MARK_USED(move_types)
560 :
561 : ! calculate rotation matrix (using quanternions)
562 1650 : rnd = rng_stream%next()
563 1650 : a1 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
564 1650 : rnd = rng_stream%next()
565 1650 : a2 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
566 1650 : rnd = rng_stream%next()
567 1650 : a3 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
568 1650 : q0 = COS(a2/2)*COS((a1 + a3)/2.0_dp)
569 1650 : q1 = SIN(a2/2)*COS((a1 - a3)/2.0_dp)
570 1650 : q2 = SIN(a2/2)*SIN((a1 - a3)/2.0_dp)
571 1650 : 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 16500 : 2*(q1*q3 - q0*q2), 2*(q2*q3 + q0*q1), q0*q0 - q1*q1 - q2*q2 + q3*q3/), (/3, 3/))
575 :
576 4950 : 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 1650 : center=elem_center)
580 :
581 : ! proceed rotation
582 1650 : atom_loop: DO i = ind_start, ind_end + dim_per_elem - 1, dim_per_elem
583 109648 : pos(i:i + 2) = MATMUL(pos(i:i + 2) - elem_center(:), rot) + elem_center(:)
584 : END DO atom_loop
585 1650 : DEALLOCATE (elem_center)
586 1650 : 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 0 : 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 0 : 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 0 : rnd1 = rng_stream%next()
616 0 : rnd2 = rng_stream%next()
617 :
618 0 : 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 0 : 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 0 : rnd3 = rng_stream%next()
631 0 : IF (rnd3 .GE. 0.5 .AND. rnd_sign_change) THEN
632 : d = -1
633 : ELSE
634 0 : d = 1
635 : END IF
636 0 : vel = SIN(phi)*delta_vel + COS(phi)*vel*d*1.0_dp
637 0 : 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 12 : 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 12 : INTEGER, DIMENSION(:), POINTER :: mol_arr
665 : REAL(KIND=dp) :: rnd
666 :
667 12 : NULLIFY (mol_arr)
668 :
669 0 : CPASSERT(ASSOCIATED(elem))
670 12 : CPASSERT(ASSOCIATED(tmc_params))
671 :
672 : ! start the timing
673 12 : CALL timeset(routineN, handle)
674 :
675 12 : short_loop = .FALSE.
676 12 : counter = 0
677 3468 : nr_mol = MAXVAL(elem%mol(:))
678 : ! ind_arr: one array element for each molecule
679 36 : ALLOCATE (mol_arr(nr_mol))
680 1164 : mol_arr(:) = -1
681 : donor_acceptor = not_selected
682 : ! select randomly if neighboring molecule is donor / acceptor
683 12 : IF (rng_stream%next() .LT. 0.5_dp) THEN
684 7 : donor_acceptor = proton_acceptor
685 : ELSE
686 5 : donor_acceptor = proton_donor
687 : END IF
688 :
689 : ! first step build loop
690 : ! select randomly one atom
691 12 : rnd = rng_stream%next()
692 : ! the randomly selected first atom
693 12 : mol = INT(rnd*nr_mol) + 1
694 12 : counter = counter + 1
695 12 : mol_arr(counter) = mol
696 :
697 : ! do until the loop is closed
698 : ! (until path connects back to any spot of the path)
699 162 : chain_completition_loop: DO
700 174 : 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 174 : rng_stream=rng_stream)
706 15784 : IF (ANY(mol_arr(:) .EQ. mol)) &
707 : EXIT chain_completition_loop
708 174 : 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 70 : DO k = 1, counter
714 70 : IF (mol_arr(k) .EQ. mol) &
715 12 : EXIT
716 : END DO
717 256 : mol_arr(1:counter - k + 1) = mol_arr(k:counter)
718 12 : counter = counter - k + 1
719 :
720 : ! check if loop is minimum size of 6 molecules
721 12 : IF (counter .LT. 6) THEN
722 : CALL cp_warn(__LOCATION__, &
723 : "short proton loop with"//cp_to_string(counter)// &
724 0 : "molecules.")
725 0 : tmp_chr = ""
726 0 : WRITE (tmp_chr, *) mol_arr(1:counter)
727 0 : CPWARN("selected molecules:"//TRIM(tmp_chr))
728 0 : 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 12 : mol_arr_in=mol_arr(1:counter), donor_acceptor=donor_acceptor)
735 12 : DEALLOCATE (mol_arr)
736 :
737 : ! end the timing
738 12 : CALL timestop(handle)
739 24 : 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 174 : 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 174 : REAL(KIND=dp), DIMENSION(:), POINTER :: distH1, distH2, distO
767 :
768 174 : NULLIFY (distO, distH1, distH2)
769 0 : CPASSERT(ASSOCIATED(elem))
770 174 : CPASSERT(ASSOCIATED(tmc_params))
771 :
772 : ! start the timing
773 174 : CALL timeset(routineN, handle)
774 :
775 50286 : nr_mol = MAXVAL(elem%mol)
776 522 : ALLOCATE (distO(nr_mol))
777 348 : ALLOCATE (distH1(nr_mol))
778 348 : ALLOCATE (distH2(nr_mol))
779 : !-- initialize the distances to huge values
780 : ! distance of nearest proton of certain molecule to preselected O
781 16878 : distO(:) = HUGE(distO(1))
782 : ! distance of (first) proton of preselected molecule to certain molecule
783 16878 : distH1(:) = HUGE(distH1(1))
784 : ! distance of (second) proton of preselected molecule to certain molecule
785 16878 : 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 174 : start_ind=ind, end_ind=ind_e)
790 :
791 : ! calculate distances to all molecules
792 16878 : list_distances: DO mol_tmp = 1, nr_mol
793 16704 : 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 16530 : mol=mol_tmp, start_ind=ind_n, end_ind=ind_e)
798 : ! check if selected molecule is water respectively consists of 3 atoms
799 16530 : 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 0 : "the proton reordering does not support, skip molecule")
803 0 : CYCLE list_distances
804 : END IF
805 16530 : IF (donor_acceptor .EQ. proton_acceptor) THEN
806 9785 : 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 4979 : 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 4979 : cell=tmc_params%cell, box_scale=elem%box_scale)
820 : END IF
821 : END IF
822 : !check for neighboring proton donors
823 33234 : IF (donor_acceptor .EQ. proton_donor) THEN
824 6745 : 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 3315 : 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 3315 : cell=tmc_params%cell, box_scale=elem%box_scale)
837 3315 : 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 174 : mol_tmp = 1
843 : ! select the nearest neighbors
844 : !check for neighboring proton acceptors
845 174 : IF (donor_acceptor .EQ. proton_acceptor) THEN
846 10094 : neighbor_mol(mol_tmp) = MINLOC(distH1(:), 1)
847 10094 : 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 103 : IF (neighbor_mol(mol_tmp) .EQ. neighbor_mol(mol_tmp + 1)) THEN
850 0 : distH1(neighbor_mol(mol_tmp)) = HUGE(distH1(1))
851 0 : distH2(neighbor_mol(mol_tmp + 1)) = HUGE(distH2(1))
852 0 : IF (MINVAL(distH1(:), 1) .LT. MINVAL(distH2(:), 1)) THEN
853 0 : neighbor_mol(mol_tmp) = MINLOC(distH1(:), 1)
854 : ELSE
855 0 : neighbor_mol(mol_tmp + 1) = MINLOC(distH2(:), 1)
856 : END IF
857 : END IF
858 103 : mol_tmp = mol_tmp + 2
859 : END IF
860 :
861 : !check for neighboring proton donors
862 174 : IF (donor_acceptor .EQ. proton_donor) THEN
863 6958 : neighbor_mol(mol_tmp) = MINLOC(distO(:), 1)
864 71 : distO(neighbor_mol(mol_tmp)) = HUGE(distO(1))
865 6958 : neighbor_mol(mol_tmp + 1) = MINLOC(distO(:), 1)
866 : END IF
867 :
868 : ! select randomly the next neighboring molecule
869 174 : rnd = rng_stream%next()
870 : ! the randomly selected atom: return value!
871 174 : mol_tmp = neighbor_mol(INT(rnd*SIZE(neighbor_mol(:))) + 1)
872 174 : mol = mol_tmp
873 :
874 174 : DEALLOCATE (distO)
875 174 : DEALLOCATE (distH1)
876 174 : DEALLOCATE (distH2)
877 :
878 : ! end the timing
879 174 : CALL timestop(handle)
880 522 : 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 16530 : 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 16530 : CPASSERT(ASSOCIATED(elem))
902 16530 : CPASSERT(i_orig .GE. 1 .AND. i_orig .LE. SIZE(elem%pos))
903 16530 : CPASSERT(i_neighbor .GE. 1 .AND. i_neighbor .LE. SIZE(elem%pos))
904 16530 : 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 16530 : 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 16530 : 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 16530 : 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 16530 : cell=tmc_params%cell, box_scale=elem%box_scale)
930 :
931 99180 : IF (MINLOC(distances(:), 1) .LE. 2) THEN
932 : donor_acceptor = proton_acceptor
933 : ELSE
934 8121 : donor_acceptor = proton_donor
935 : END IF
936 16530 : 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 12 : 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 12 : 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 12 : NULLIFY (ind_arr, tmp_cell)
966 :
967 0 : CPASSERT(ASSOCIATED(tmc_params))
968 12 : CPASSERT(ASSOCIATED(elem))
969 :
970 : ! start the timing
971 12 : CALL timeset(routineN, handle)
972 :
973 36 : ALLOCATE (ind_arr(0:SIZE(mol_arr_in) + 1))
974 128 : 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 128 : start_ind=ind_arr(i), end_ind=ind)
978 : END DO
979 12 : ind_arr(0) = ind_arr(SIZE(ind_arr) - 2)
980 12 : ind_arr(SIZE(ind_arr) - 1) = ind_arr(1)
981 :
982 : ! get the scaled cell
983 336 : ALLOCATE (tmp_cell)
984 : CALL get_scaled_cell(cell=tmc_params%cell, box_scale=elem%box_scale, &
985 12 : scaled_cell=tmp_cell)
986 :
987 : ! rotate single molecules
988 128 : DO i = 1, SIZE(ind_arr) - 2
989 : ! the 3 O atoms
990 464 : vec_1O(:) = elem%pos(ind_arr(i - 1):ind_arr(i - 1) + tmc_params%dim_per_elem - 1)
991 464 : vec_2O(:) = elem%pos(ind_arr(i):ind_arr(i) + tmc_params%dim_per_elem - 1)
992 464 : 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 116 : 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 276 : 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 276 : 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 188 : 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 188 : ind_arr(i) + 3*tmc_params%dim_per_elem - 1)
1021 : H_offset = 2
1022 : END IF
1023 :
1024 12 : IF (.TRUE.) THEN !TODO find a better switch for the pauling model
1025 :
1026 : ! do rotation (NOT pauling model)
1027 464 : tmp_1 = pbc(vec_2O - vec_1O, tmp_cell)
1028 464 : tmp_2 = pbc(vec_3O - vec_2H_f, tmp_cell)
1029 :
1030 464 : dihe_angle = donor_acceptor*dihedral_angle(tmp_1, vec_2H_f - vec_2O, tmp_2)
1031 464 : 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 2436 : 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 1508 : 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 12 : DEALLOCATE (tmp_cell)
1081 12 : DEALLOCATE (ind_arr)
1082 : ! end the timing
1083 12 : CALL timestop(handle)
1084 24 : 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 224 : 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 224 : REAL(KIND=dp), DIMENSION(:), POINTER :: disp, scaling
1116 :
1117 224 : NULLIFY (scaling, disp)
1118 :
1119 0 : CPASSERT(ASSOCIATED(conf))
1120 224 : CPASSERT(ASSOCIATED(move_types))
1121 224 : CPASSERT(ASSOCIATED(tmc_params))
1122 224 : CPASSERT(T_ind .GT. 0 .AND. T_ind .LE. tmc_params%nr_temp)
1123 224 : CPASSERT(tmc_params%dim_per_elem .EQ. 3)
1124 224 : CPASSERT(tmc_params%cell%orthorhombic)
1125 :
1126 : ! start the timing
1127 224 : CALL timeset(routineN, handle)
1128 :
1129 672 : ALLOCATE (scaling(tmc_params%dim_per_elem))
1130 672 : ALLOCATE (disp(tmc_params%dim_per_elem))
1131 :
1132 896 : 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 224 : 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 224 : IF (tmc_params%v_isotropic) THEN
1160 224 : rnd = rng_stream%next()
1161 : box_length_new(:) = box_length_new(:) + &
1162 : (rnd - 0.5_dp)*2.0_dp* &
1163 896 : move_types%mv_size(mv_type_volume_move, T_ind)
1164 : ELSE
1165 : ! select a random direction
1166 0 : rnd = rng_stream%next()
1167 0 : dir = 1 + INT(rnd*3)
1168 0 : rnd = rng_stream%next()
1169 : box_length_new(dir) = box_length_new(dir) + &
1170 : (rnd - 0.5_dp)*2.0_dp* &
1171 0 : move_types%mv_size(mv_type_volume_move, T_ind)
1172 : END IF
1173 : END IF
1174 :
1175 : ! get the original box length
1176 896 : scaling(:) = 1.0_dp
1177 : CALL get_scaled_cell(cell=tmc_params%cell, &
1178 : box_scale=scaling, &
1179 224 : abc=box_length_orig)
1180 : ! get the new box scale
1181 896 : conf%box_scale(:) = box_length_new(:)/box_length_orig(:)
1182 : ! molecule scaling
1183 1792 : scaling(:) = conf%box_scale(:)/box_scale_old(:)
1184 :
1185 224 : IF (mv_cen_of_mass .EQV. .FALSE.) THEN
1186 : ! homogene scaling of atomic coordinates
1187 224 : DO atom = 1, SIZE(conf%pos), tmc_params%dim_per_elem
1188 : conf%pos(atom:atom + tmc_params%dim_per_elem - 1) = &
1189 182880 : conf%pos(atom:atom + tmc_params%dim_per_elem - 1)*scaling(:)
1190 : END DO
1191 : ELSE
1192 0 : DO mol = 1, MAXVAL(conf%mol(:))
1193 : ! move the molecule related to the molecule center of mass
1194 : ! get center of mass
1195 0 : CPASSERT(ASSOCIATED(tmc_params%atoms))
1196 :
1197 : CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=conf%mol, mol=mol, &
1198 0 : 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 0 : center=disp)
1204 : ! calculate the center of mass DISPLACEMENT
1205 0 : disp(:) = disp(:)*(scaling(:) - 1.0_dp)
1206 : ! displace all atoms of the molecule
1207 0 : 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 0 : conf%pos(atom:atom + tmc_params%dim_per_elem - 1) + disp(:)
1210 : END DO
1211 : END DO
1212 : END IF
1213 :
1214 224 : DEALLOCATE (scaling)
1215 224 : DEALLOCATE (disp)
1216 :
1217 : ! end the timing
1218 224 : CALL timestop(handle)
1219 448 : 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 185 : 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 185 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pos_tmp
1239 :
1240 185 : CPASSERT(ASSOCIATED(conf))
1241 185 : CPASSERT(ASSOCIATED(move_types))
1242 185 : CPASSERT(ASSOCIATED(tmc_params))
1243 185 : 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 532 : 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 532 : rng_stream%next()) + 1
1253 : ! check if they have different kinds
1254 532 : 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 234 : IF (ASSOCIATED(move_types%atom_lists)) THEN
1257 338 : DO ind_1 = 1, SIZE(move_types%atom_lists)
1258 : IF (ANY(move_types%atom_lists(ind_1)%atoms(:) .EQ. &
1259 1073 : tmc_params%atoms(a_1)%name) .AND. &
1260 : ANY(move_types%atom_lists(ind_1)%atoms(:) .EQ. &
1261 49 : 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 555 : ALLOCATE (pos_tmp(tmc_params%dim_per_elem))
1275 185 : ind_1 = (a_1 - 1)*tmc_params%dim_per_elem + 1
1276 740 : pos_tmp(:) = conf%pos(ind_1:ind_1 + tmc_params%dim_per_elem - 1)
1277 185 : 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 1295 : conf%pos(ind_2:ind_2 + tmc_params%dim_per_elem - 1)
1280 740 : conf%pos(ind_2:ind_2 + tmc_params%dim_per_elem - 1) = pos_tmp(:)
1281 185 : DEALLOCATE (pos_tmp)
1282 : END IF
1283 185 : END SUBROUTINE swap_atoms
1284 :
1285 : END MODULE tmc_moves
|
