Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief the various moves in Monte Carlo (MC) simulations, including
10 : !> change of internal conformation, translation of a molecule, rotation
11 : !> of a molecule, and changing the size of the simulation box
12 : !> \par History
13 : !> none
14 : !> \author Matthew J. McGrath (10.16.2003)
15 : ! **************************************************************************************************
16 : MODULE mc_moves
17 : USE atomic_kind_types, ONLY: get_atomic_kind
18 : USE cell_methods, ONLY: cell_create
19 : USE cell_types, ONLY: cell_clone,&
20 : cell_release,&
21 : cell_type,&
22 : get_cell
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_get_default_io_unit,&
25 : cp_logger_type
26 : USE cp_subsys_types, ONLY: cp_subsys_get,&
27 : cp_subsys_set,&
28 : cp_subsys_type
29 : USE force_env_methods, ONLY: force_env_calc_energy_force
30 : USE force_env_types, ONLY: force_env_get,&
31 : force_env_type,&
32 : use_fist_force
33 : USE global_types, ONLY: global_environment_type
34 : USE kinds, ONLY: default_string_length,&
35 : dp
36 : USE mathconstants, ONLY: pi
37 : USE mc_coordinates, ONLY: check_for_overlap,&
38 : cluster_search,&
39 : create_discrete_array,&
40 : generate_cbmc_swap_config,&
41 : get_center_of_mass
42 : USE mc_types, ONLY: get_mc_molecule_info,&
43 : get_mc_par,&
44 : mc_ekin_type,&
45 : mc_molecule_info_type,&
46 : mc_moves_type,&
47 : mc_simpar_type
48 : USE md_run, ONLY: qs_mol_dyn
49 : USE message_passing, ONLY: mp_comm_type
50 : USE molecule_kind_list_types, ONLY: molecule_kind_list_type
51 : USE molecule_kind_types, ONLY: bend_type,&
52 : bond_type,&
53 : get_molecule_kind,&
54 : molecule_kind_type,&
55 : torsion_type
56 : USE parallel_rng_types, ONLY: rng_stream_type
57 : USE particle_list_types, ONLY: particle_list_type
58 : USE physcon, ONLY: angstrom
59 : #include "../../base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 :
63 : PRIVATE
64 :
65 : PRIVATE :: change_bond_angle, change_bond_length, depth_first_search, &
66 : change_dihedral
67 :
68 : ! *** Global parameters ***
69 :
70 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_moves'
71 :
72 : PUBLIC :: mc_conformation_change, mc_molecule_translation, &
73 : mc_molecule_rotation, mc_volume_move, mc_avbmc_move, &
74 : mc_hmc_move, mc_cluster_translation
75 :
76 : CONTAINS
77 :
78 : ! **************************************************************************************************
79 : !> \brief essentially performs a depth-first search of the molecule structure
80 : !> to find all atoms connected to a specific atom excluding one branch...
81 : !> for instance, if water is labelled 1-2-3 for O-H-H, calling this
82 : !> routine with current_atom=1,avoid_atom=2 returns the array
83 : !> atom=(0,0,1)
84 : !> \param current_atom the atom whose connections we're looking at
85 : !> \param avoid_atom the atom whose direction the search is not supposed to go
86 : !> \param connectivity an array telling us the neighbors of all atoms
87 : !> \param atom the array that tells us if one can get to a given atom by
88 : !> starting at current_atom and not going through avoid_atom...0 is no,
89 : !> 1 is yes
90 : !> \author MJM
91 : ! **************************************************************************************************
92 1324 : RECURSIVE SUBROUTINE depth_first_search(current_atom, avoid_atom, &
93 1324 : connectivity, atom)
94 :
95 : INTEGER, INTENT(IN) :: current_atom, avoid_atom
96 : INTEGER, DIMENSION(:, :), INTENT(IN) :: connectivity
97 : INTEGER, DIMENSION(:), INTENT(INOUT) :: atom
98 :
99 : INTEGER :: iatom
100 :
101 2960 : DO iatom = 1, 6
102 2960 : IF (connectivity(iatom, current_atom) /= 0) THEN
103 1636 : IF (connectivity(iatom, current_atom) /= avoid_atom) THEN
104 312 : atom(connectivity(iatom, current_atom)) = 1
105 : CALL depth_first_search(connectivity(iatom, current_atom), &
106 312 : current_atom, connectivity, atom)
107 : END IF
108 : ELSE
109 : RETURN
110 : END IF
111 : END DO
112 :
113 : END SUBROUTINE depth_first_search
114 :
115 : ! **************************************************************************************************
116 : !> \brief performs either a bond or angle change move for a given molecule
117 : !> \param mc_par the mc parameters for the force env
118 : !> \param force_env the force environment used in the move
119 : !> \param bias_env the force environment used to bias the move, if any (it may
120 : !> be null if lbias=.false. in mc_par)
121 : !> \param moves the structure that keeps track of how many moves have been
122 : !> accepted/rejected
123 : !> \param move_updates the structure that keeps track of how many moves have
124 : !> been accepted/rejected since the last time the displacements
125 : !> were updated
126 : !> \param start_atom the number of the molecule's first atom, assuming the rest
127 : !> of the atoms follow sequentially
128 : !> \param molecule_type the type of the molecule we're moving
129 : !> \param box_number the box the molecule is in
130 : !> \param bias_energy the biased energy of the system before the move
131 : !> \param move_type dictates what kind of conformational change we do
132 : !> \param lreject set to .true. if there is an overlap
133 : !> \param rng_stream the random number stream that we draw from
134 : !> \author MJM
135 : ! **************************************************************************************************
136 506 : SUBROUTINE mc_conformation_change(mc_par, force_env, bias_env, moves, &
137 : move_updates, start_atom, molecule_type, box_number, &
138 : bias_energy, move_type, lreject, &
139 : rng_stream)
140 :
141 : TYPE(mc_simpar_type), POINTER :: mc_par
142 : TYPE(force_env_type), POINTER :: force_env, bias_env
143 : TYPE(mc_moves_type), POINTER :: moves, move_updates
144 : INTEGER, INTENT(IN) :: start_atom, molecule_type, box_number
145 : REAL(KIND=dp), INTENT(INOUT) :: bias_energy
146 : CHARACTER(LEN=*), INTENT(IN) :: move_type
147 : LOGICAL, INTENT(OUT) :: lreject
148 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
149 :
150 : CHARACTER(len=*), PARAMETER :: routineN = 'mc_conformation_change'
151 :
152 : CHARACTER(default_string_length) :: name
153 : CHARACTER(default_string_length), DIMENSION(:), &
154 506 : POINTER :: names
155 : INTEGER :: atom_number, end_atom, end_mol, handle, imol_type, imolecule, ipart, jbox, &
156 : molecule_number, nunits_mol, source, start_mol
157 506 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits
158 506 : INTEGER, DIMENSION(:, :), POINTER :: nchains
159 : LOGICAL :: ionode, lbias, loverlap
160 : REAL(KIND=dp) :: BETA, bias_energy_new, bias_energy_old, &
161 : dis_length, exp_max_val, exp_min_val, &
162 : rand, value, w
163 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r_new, r_old
164 : TYPE(cp_subsys_type), POINTER :: subsys
165 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
166 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
167 : TYPE(molecule_kind_type), POINTER :: molecule_kind, molecule_kind_test
168 : TYPE(mp_comm_type) :: group
169 : TYPE(particle_list_type), POINTER :: particles
170 :
171 : ! begin the timing of the subroutine
172 :
173 506 : CALL timeset(routineN, handle)
174 :
175 : ! nullify some pointers
176 506 : NULLIFY (particles, subsys, molecule_kinds, molecule_kind, &
177 506 : molecule_kind_test)
178 :
179 : ! get a bunch of stuff from mc_par
180 : CALL get_mc_par(mc_par, lbias=lbias, mc_molecule_info=mc_molecule_info, &
181 : BETA=BETA, exp_max_val=exp_max_val, &
182 506 : exp_min_val=exp_min_val, group=group, source=source, ionode=ionode)
183 : CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, nunits=nunits, &
184 506 : mol_type=mol_type, names=names)
185 :
186 : ! do some allocation
187 506 : nunits_mol = nunits(molecule_type)
188 1518 : ALLOCATE (r_old(1:3, 1:nunits_mol))
189 1012 : ALLOCATE (r_new(1:3, 1:nunits_mol))
190 :
191 : ! find out some bounds for mol_type
192 506 : start_mol = 1
193 506 : DO jbox = 1, box_number - 1
194 506 : start_mol = start_mol + SUM(nchains(:, jbox))
195 : END DO
196 1518 : end_mol = start_mol + SUM(nchains(:, box_number)) - 1
197 :
198 : ! figure out which molecule number we are
199 506 : end_atom = start_atom + nunits_mol - 1
200 506 : molecule_number = 0
201 506 : atom_number = 1
202 2974 : DO imolecule = 1, SUM(nchains(:, box_number))
203 1962 : IF (atom_number == start_atom) THEN
204 506 : molecule_number = imolecule
205 506 : EXIT
206 : END IF
207 1456 : atom_number = atom_number + nunits(mol_type(imolecule + start_mol - 1))
208 : END DO
209 506 : IF (molecule_number == 0) CPABORT('Cannot find the molecule number')
210 :
211 : ! are we biasing this move?
212 506 : IF (lbias) THEN
213 :
214 : ! grab the coordinates
215 420 : CALL force_env_get(bias_env, subsys=subsys)
216 : ! save the energy
217 420 : bias_energy_old = bias_energy
218 :
219 : ELSE
220 :
221 : ! grab the coordinates
222 86 : CALL force_env_get(force_env, subsys=subsys)
223 : END IF
224 :
225 : ! now find the molecule type associated with this guy
226 : CALL cp_subsys_get(subsys, &
227 506 : particles=particles, molecule_kinds=molecule_kinds)
228 506 : DO imol_type = 1, SIZE(molecule_kinds%els(:))
229 506 : molecule_kind_test => molecule_kinds%els(imol_type)
230 506 : CALL get_molecule_kind(molecule_kind_test, name=name)
231 506 : IF (TRIM(ADJUSTL(name)) == TRIM(ADJUSTL(names(molecule_type)))) THEN
232 506 : molecule_kind => molecule_kinds%els(imol_type)
233 506 : EXIT
234 : END IF
235 : END DO
236 :
237 : ! save the coordinates
238 2024 : DO ipart = start_atom, end_atom
239 6578 : r_old(1:3, ipart - start_atom + 1) = particles%els(ipart)%r(1:3)
240 : END DO
241 :
242 506 : IF (.NOT. ASSOCIATED(molecule_kind)) CPABORT('Cannot find the molecule type')
243 : ! do the move
244 506 : IF (move_type == 'bond') THEN
245 :
246 : ! record the attempt
247 312 : moves%bond%attempts = moves%bond%attempts + 1
248 312 : move_updates%bond%attempts = move_updates%bond%attempts + 1
249 312 : moves%bias_bond%attempts = moves%bias_bond%attempts + 1
250 312 : move_updates%bias_bond%attempts = move_updates%bias_bond%attempts + 1
251 312 : IF (.NOT. lbias) THEN
252 48 : moves%bond%qsuccesses = moves%bond%qsuccesses + 1
253 : move_updates%bond%qsuccesses = &
254 48 : move_updates%bond%qsuccesses + 1
255 48 : moves%bias_bond%qsuccesses = moves%bias_bond%qsuccesses + 1
256 : move_updates%bias_bond%qsuccesses = &
257 48 : move_updates%bias_bond%qsuccesses + 1
258 : END IF
259 :
260 : ! do the move
261 : CALL change_bond_length(r_old, r_new, mc_par, molecule_type, &
262 312 : molecule_kind, dis_length, particles, rng_stream)
263 :
264 194 : ELSE IF (move_type == 'angle') THEN
265 :
266 : ! record the attempt
267 194 : moves%angle%attempts = moves%angle%attempts + 1
268 194 : move_updates%angle%attempts = move_updates%angle%attempts + 1
269 194 : moves%bias_angle%attempts = moves%bias_angle%attempts + 1
270 194 : move_updates%bias_angle%attempts = move_updates%bias_angle%attempts + 1
271 194 : IF (.NOT. lbias) THEN
272 38 : moves%angle%qsuccesses = moves%angle%qsuccesses + 1
273 : move_updates%angle%qsuccesses = &
274 38 : move_updates%angle%qsuccesses + 1
275 38 : moves%bias_angle%qsuccesses = moves%bias_angle%qsuccesses + 1
276 : move_updates%bias_angle%qsuccesses = &
277 38 : move_updates%bias_angle%qsuccesses + 1
278 : END IF
279 :
280 : ! do the move
281 : CALL change_bond_angle(r_old, r_new, mc_par, molecule_type, &
282 194 : molecule_kind, particles, rng_stream)
283 194 : dis_length = 1.0E0_dp
284 : ELSE
285 : ! record the attempt
286 0 : moves%dihedral%attempts = moves%dihedral%attempts + 1
287 0 : move_updates%dihedral%attempts = move_updates%dihedral%attempts + 1
288 0 : moves%bias_dihedral%attempts = moves%bias_dihedral%attempts + 1
289 0 : move_updates%bias_dihedral%attempts = move_updates%bias_dihedral%attempts + 1
290 0 : IF (.NOT. lbias) THEN
291 0 : moves%dihedral%qsuccesses = moves%dihedral%qsuccesses + 1
292 : move_updates%dihedral%qsuccesses = &
293 0 : move_updates%dihedral%qsuccesses + 1
294 0 : moves%bias_dihedral%qsuccesses = moves%bias_dihedral%qsuccesses + 1
295 : move_updates%bias_dihedral%qsuccesses = &
296 0 : move_updates%bias_dihedral%qsuccesses + 1
297 : END IF
298 :
299 : ! do the move
300 : CALL change_dihedral(r_old, r_new, mc_par, molecule_type, &
301 0 : molecule_kind, particles, rng_stream)
302 0 : dis_length = 1.0E0_dp
303 :
304 : END IF
305 :
306 : ! set the coordinates
307 2024 : DO ipart = start_atom, end_atom
308 6578 : particles%els(ipart)%r(1:3) = r_new(1:3, ipart - start_atom + 1)
309 : END DO
310 :
311 : ! check for overlap
312 506 : lreject = .FALSE.
313 506 : IF (lbias) THEN
314 : CALL check_for_overlap(bias_env, nchains(:, box_number), &
315 : nunits(:), loverlap, mol_type(start_mol:end_mol), &
316 420 : molecule_number=molecule_number)
317 : ELSE
318 : CALL check_for_overlap(force_env, nchains(:, box_number), &
319 : nunits(:), loverlap, mol_type(start_mol:end_mol), &
320 86 : molecule_number=molecule_number)
321 86 : IF (loverlap) lreject = .TRUE.
322 : END IF
323 :
324 : ! if we're biasing classical, check for acceptance
325 506 : IF (lbias) THEN
326 :
327 : ! here's where we bias the moves
328 :
329 420 : IF (loverlap) THEN
330 : w = 0.0E0_dp
331 : ELSE
332 420 : CALL force_env_calc_energy_force(bias_env, calc_force=.FALSE.)
333 : CALL force_env_get(bias_env, &
334 420 : potential_energy=bias_energy_new)
335 : ! accept or reject the move based on the Metropolis rule with a
336 : ! correction factor for the change in phase space...dis_length is
337 : ! made unitless in change_bond_length
338 420 : value = -BETA*(bias_energy_new - bias_energy_old)
339 420 : IF (value > exp_max_val) THEN
340 : w = 10.0_dp
341 420 : ELSE IF (value < exp_min_val) THEN
342 : w = 0.0_dp
343 : ELSE
344 420 : w = EXP(value)*dis_length**2
345 : END IF
346 :
347 : END IF
348 :
349 420 : IF (w >= 1.0E0_dp) THEN
350 194 : w = 1.0E0_dp
351 194 : rand = 0.0E0_dp
352 : ELSE
353 226 : IF (ionode) THEN
354 113 : rand = rng_stream%next()
355 : END IF
356 226 : CALL group%bcast(rand, source)
357 : END IF
358 :
359 420 : IF (rand < w) THEN
360 :
361 : ! accept the move
362 252 : IF (move_type == 'bond') THEN
363 140 : moves%bond%qsuccesses = moves%bond%qsuccesses + 1
364 : move_updates%bond%successes = &
365 140 : move_updates%bond%successes + 1
366 140 : moves%bias_bond%successes = moves%bias_bond%successes + 1
367 : move_updates%bias_bond%successes = &
368 140 : move_updates%bias_bond%successes + 1
369 112 : ELSE IF (move_type == 'angle') THEN
370 112 : moves%angle%qsuccesses = moves%angle%qsuccesses + 1
371 : move_updates%angle%successes = &
372 112 : move_updates%angle%successes + 1
373 112 : moves%bias_angle%successes = moves%bias_angle%successes + 1
374 : move_updates%bias_angle%successes = &
375 112 : move_updates%bias_angle%successes + 1
376 : ELSE
377 0 : moves%dihedral%qsuccesses = moves%dihedral%qsuccesses + 1
378 : move_updates%dihedral%successes = &
379 0 : move_updates%dihedral%successes + 1
380 0 : moves%bias_dihedral%successes = moves%bias_dihedral%successes + 1
381 : move_updates%bias_dihedral%successes = &
382 0 : move_updates%bias_dihedral%successes + 1
383 : END IF
384 :
385 : bias_energy = bias_energy + bias_energy_new - &
386 252 : bias_energy_old
387 :
388 : ELSE
389 :
390 : ! reject the move
391 : ! restore the coordinates
392 168 : CALL force_env_get(bias_env, subsys=subsys)
393 168 : CALL cp_subsys_get(subsys, particles=particles)
394 672 : DO ipart = start_atom, end_atom
395 2184 : particles%els(ipart)%r(1:3) = r_old(1:3, ipart - start_atom + 1)
396 : END DO
397 168 : CALL cp_subsys_set(subsys, particles=particles)
398 :
399 : END IF
400 :
401 : END IF
402 :
403 : ! deallocate some stuff
404 506 : DEALLOCATE (r_old)
405 506 : DEALLOCATE (r_new)
406 :
407 : ! end the timing
408 506 : CALL timestop(handle)
409 :
410 506 : END SUBROUTINE mc_conformation_change
411 :
412 : ! **************************************************************************************************
413 : !> \brief translates the given molecule randomly in either the x,y, or z direction
414 : !> \param mc_par the mc parameters for the force env
415 : !> \param force_env the force environment used in the move
416 : !> \param bias_env the force environment used to bias the move, if any (it may
417 : !> be null if lbias=.false. in mc_par)
418 : !> \param moves the structure that keeps track of how many moves have been
419 : !> accepted/rejected
420 : !> \param move_updates the structure that keeps track of how many moves have
421 : !> been accepted/rejected since the last time the displacements
422 : !> were updated
423 : !> \param start_atom the number of the molecule's first atom, assuming the rest of
424 : !> the atoms follow sequentially
425 : !> \param box_number the box the molecule is in
426 : !> \param bias_energy the biased energy of the system before the move
427 : !> \param molecule_type the type of molecule we're moving
428 : !> \param lreject set to .true. if there is an overlap
429 : !> \param rng_stream the random number stream that we draw from
430 : !> \author MJM
431 : ! **************************************************************************************************
432 624 : SUBROUTINE mc_molecule_translation(mc_par, force_env, bias_env, moves, &
433 : move_updates, start_atom, box_number, &
434 : bias_energy, molecule_type, &
435 : lreject, rng_stream)
436 :
437 : TYPE(mc_simpar_type), POINTER :: mc_par
438 : TYPE(force_env_type), POINTER :: force_env, bias_env
439 : TYPE(mc_moves_type), POINTER :: moves, move_updates
440 : INTEGER, INTENT(IN) :: start_atom, box_number
441 : REAL(KIND=dp), INTENT(INOUT) :: bias_energy
442 : INTEGER, INTENT(IN) :: molecule_type
443 : LOGICAL, INTENT(OUT) :: lreject
444 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
445 :
446 : CHARACTER(len=*), PARAMETER :: routineN = 'mc_molecule_translation'
447 :
448 : INTEGER :: atom_number, end_atom, end_mol, handle, imolecule, ipart, iparticle, jbox, &
449 : molecule_number, move_direction, nunits_mol, source, start_mol
450 624 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
451 624 : INTEGER, DIMENSION(:, :), POINTER :: nchains
452 : LOGICAL :: ionode, lbias, loverlap
453 624 : REAL(dp), DIMENSION(:), POINTER :: rmtrans
454 : REAL(KIND=dp) :: BETA, bias_energy_new, bias_energy_old, &
455 : dis_mol, exp_max_val, exp_min_val, &
456 : rand, value, w
457 624 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r_old
458 : TYPE(cp_subsys_type), POINTER :: subsys
459 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
460 : TYPE(mp_comm_type) :: group
461 : TYPE(particle_list_type), POINTER :: particles
462 :
463 : ! *** Local Counters ***
464 : ! begin the timing of the subroutine
465 :
466 624 : CALL timeset(routineN, handle)
467 :
468 : ! nullify some pointers
469 624 : NULLIFY (particles, subsys)
470 :
471 : ! get a bunch of stuff from mc_par
472 : CALL get_mc_par(mc_par, lbias=lbias, &
473 : BETA=BETA, exp_max_val=exp_max_val, &
474 : exp_min_val=exp_min_val, rmtrans=rmtrans, ionode=ionode, source=source, &
475 624 : group=group, mc_molecule_info=mc_molecule_info)
476 : CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot, &
477 624 : nchains=nchains, nunits=nunits, mol_type=mol_type)
478 :
479 : ! find out some bounds for mol_type
480 624 : start_mol = 1
481 640 : DO jbox = 1, box_number - 1
482 672 : start_mol = start_mol + SUM(nchains(:, jbox))
483 : END DO
484 1872 : end_mol = start_mol + SUM(nchains(:, box_number)) - 1
485 :
486 : ! do some allocation
487 1872 : ALLOCATE (r_old(1:3, 1:nunits_tot(box_number)))
488 :
489 : ! find the index of the last atom of this molecule, and the molecule number
490 624 : nunits_mol = nunits(molecule_type)
491 624 : end_atom = start_atom + nunits_mol - 1
492 624 : molecule_number = 0
493 624 : atom_number = 1
494 5770 : DO imolecule = 1, SUM(nchains(:, box_number))
495 4522 : IF (atom_number == start_atom) THEN
496 624 : molecule_number = imolecule
497 624 : EXIT
498 : END IF
499 3898 : atom_number = atom_number + nunits(mol_type(imolecule + start_mol - 1))
500 : END DO
501 624 : IF (molecule_number == 0) CPABORT('Cannot find the molecule number')
502 :
503 : ! are we biasing this move?
504 624 : IF (lbias) THEN
505 :
506 : ! grab the coordinates
507 528 : CALL force_env_get(bias_env, subsys=subsys)
508 528 : CALL cp_subsys_get(subsys, particles=particles)
509 :
510 : ! save the coordinates
511 14710 : DO ipart = 1, nunits_tot(box_number)
512 57256 : r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
513 : END DO
514 :
515 : ! save the energy
516 528 : bias_energy_old = bias_energy
517 :
518 : ELSE
519 :
520 : ! grab the coordinates
521 96 : CALL force_env_get(force_env, subsys=subsys)
522 96 : CALL cp_subsys_get(subsys, particles=particles)
523 : END IF
524 :
525 : ! record the attempt
526 624 : moves%trans%attempts = moves%trans%attempts + 1
527 624 : move_updates%trans%attempts = move_updates%trans%attempts + 1
528 624 : moves%bias_trans%attempts = moves%bias_trans%attempts + 1
529 624 : move_updates%bias_trans%attempts = move_updates%bias_trans%attempts + 1
530 624 : IF (.NOT. lbias) THEN
531 96 : moves%trans%qsuccesses = moves%trans%qsuccesses + 1
532 96 : move_updates%trans%qsuccesses = move_updates%trans%qsuccesses + 1
533 96 : moves%bias_trans%qsuccesses = moves%bias_trans%qsuccesses + 1
534 96 : move_updates%bias_trans%qsuccesses = move_updates%bias_trans%qsuccesses + 1
535 : END IF
536 :
537 : ! move one molecule in the system
538 :
539 : ! call a random number to figure out which direction we're moving
540 624 : IF (ionode) rand = rng_stream%next()
541 624 : CALL group%bcast(rand, source)
542 : ! 1,2,3 with equal prob
543 624 : move_direction = INT(3*rand) + 1
544 :
545 : ! call a random number to figure out how far we're moving
546 624 : IF (ionode) rand = rng_stream%next()
547 624 : CALL group%bcast(rand, source)
548 624 : dis_mol = rmtrans(molecule_type)*(rand - 0.5E0_dp)*2.0E0_dp
549 :
550 : ! do the move
551 1896 : DO iparticle = start_atom, end_atom
552 : particles%els(iparticle)%r(move_direction) = &
553 1896 : particles%els(iparticle)%r(move_direction) + dis_mol
554 : END DO
555 624 : CALL cp_subsys_set(subsys, particles=particles)
556 :
557 : ! figure out if there is any overlap...need the number of the molecule
558 624 : lreject = .FALSE.
559 624 : IF (lbias) THEN
560 : CALL check_for_overlap(bias_env, nchains(:, box_number), &
561 : nunits(:), loverlap, mol_type(start_mol:end_mol), &
562 528 : molecule_number=molecule_number)
563 : ELSE
564 : CALL check_for_overlap(force_env, nchains(:, box_number), &
565 : nunits(:), loverlap, mol_type(start_mol:end_mol), &
566 96 : molecule_number=molecule_number)
567 96 : IF (loverlap) lreject = .TRUE.
568 : END IF
569 :
570 : ! if we're biasing with a cheaper potential, check for acceptance
571 624 : IF (lbias) THEN
572 :
573 : ! here's where we bias the moves
574 528 : IF (loverlap) THEN
575 : w = 0.0E0_dp
576 : ELSE
577 528 : CALL force_env_calc_energy_force(bias_env, calc_force=.FALSE.)
578 : CALL force_env_get(bias_env, &
579 528 : potential_energy=bias_energy_new)
580 : ! accept or reject the move based on the Metropolis rule
581 528 : value = -BETA*(bias_energy_new - bias_energy_old)
582 528 : IF (value > exp_max_val) THEN
583 : w = 10.0_dp
584 528 : ELSE IF (value < exp_min_val) THEN
585 : w = 0.0_dp
586 : ELSE
587 528 : w = EXP(value)
588 : END IF
589 :
590 : END IF
591 :
592 528 : IF (w >= 1.0E0_dp) THEN
593 258 : w = 1.0E0_dp
594 258 : rand = 0.0E0_dp
595 : ELSE
596 270 : IF (ionode) rand = rng_stream%next()
597 270 : CALL group%bcast(rand, source)
598 : END IF
599 :
600 528 : IF (rand < w) THEN
601 :
602 : ! accept the move
603 454 : moves%bias_trans%successes = moves%bias_trans%successes + 1
604 454 : move_updates%bias_trans%successes = move_updates%bias_trans%successes + 1
605 454 : moves%trans%qsuccesses = moves%trans%qsuccesses + 1
606 : move_updates%trans%successes = &
607 454 : move_updates%trans%successes + 1
608 454 : moves%qtrans_dis = moves%qtrans_dis + ABS(dis_mol)
609 : bias_energy = bias_energy + bias_energy_new - &
610 454 : bias_energy_old
611 :
612 : ELSE
613 :
614 : ! reject the move
615 : ! restore the coordinates
616 74 : CALL force_env_get(bias_env, subsys=subsys)
617 74 : CALL cp_subsys_get(subsys, particles=particles)
618 2072 : DO ipart = 1, nunits_tot(box_number)
619 8066 : particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
620 : END DO
621 74 : CALL cp_subsys_set(subsys, particles=particles)
622 :
623 : END IF
624 :
625 : END IF
626 :
627 : ! deallocate some stuff
628 624 : DEALLOCATE (r_old)
629 :
630 : ! end the timing
631 624 : CALL timestop(handle)
632 :
633 624 : END SUBROUTINE mc_molecule_translation
634 :
635 : ! **************************************************************************************************
636 : !> \brief rotates the given molecule randomly around the x,y, or z axis...
637 : !> only works for water at the moment
638 : !> \param mc_par the mc parameters for the force env
639 : !> \param force_env the force environment used in the move
640 : !> \param bias_env the force environment used to bias the move, if any (it may
641 : !> be null if lbias=.false. in mc_par)
642 : !> \param moves the structure that keeps track of how many moves have been
643 : !> accepted/rejected
644 : !> \param move_updates the structure that keeps track of how many moves have
645 : !> been accepted/rejected since the last time the displacements
646 : !> were updated
647 : !> \param box_number the box the molecule is in
648 : !> \param start_atom the number of the molecule's first atom, assuming the rest of
649 : !> the atoms follow sequentially
650 : !> \param molecule_type the type of molecule we're moving
651 : !> \param bias_energy the biased energy of the system before the move
652 : !> \param lreject set to .true. if there is an overlap
653 : !> \param rng_stream the random number stream that we draw from
654 : !> \author MJM
655 : ! **************************************************************************************************
656 502 : SUBROUTINE mc_molecule_rotation(mc_par, force_env, bias_env, moves, &
657 : move_updates, box_number, &
658 : start_atom, molecule_type, bias_energy, lreject, &
659 : rng_stream)
660 :
661 : TYPE(mc_simpar_type), POINTER :: mc_par
662 : TYPE(force_env_type), POINTER :: force_env, bias_env
663 : TYPE(mc_moves_type), POINTER :: moves, move_updates
664 : INTEGER, INTENT(IN) :: box_number, start_atom, molecule_type
665 : REAL(KIND=dp), INTENT(INOUT) :: bias_energy
666 : LOGICAL, INTENT(OUT) :: lreject
667 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
668 :
669 : CHARACTER(len=*), PARAMETER :: routineN = 'mc_molecule_rotation'
670 :
671 : INTEGER :: atom_number, dir, end_atom, end_mol, handle, ii, imolecule, ipart, iunit, jbox, &
672 : molecule_number, nunits_mol, source, start_mol
673 502 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
674 502 : INTEGER, DIMENSION(:, :), POINTER :: nchains
675 : LOGICAL :: ionode, lbias, loverlap, lx, ly
676 502 : REAL(dp), DIMENSION(:), POINTER :: rmrot
677 502 : REAL(dp), DIMENSION(:, :), POINTER :: mass
678 : REAL(KIND=dp) :: BETA, bias_energy_new, bias_energy_old, cosdg, dgamma, exp_max_val, &
679 : exp_min_val, masstot, nxcm, nycm, nzcm, rand, rx, rxnew, ry, rynew, rz, rznew, sindg, &
680 : value, w
681 502 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r_old
682 : TYPE(cp_subsys_type), POINTER :: subsys
683 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
684 : TYPE(mp_comm_type) :: group
685 : TYPE(particle_list_type), POINTER :: particles
686 :
687 : ! begin the timing of the subroutine
688 :
689 502 : CALL timeset(routineN, handle)
690 :
691 502 : NULLIFY (rmrot, subsys, particles)
692 :
693 : ! get a bunch of stuff from mc_par
694 : CALL get_mc_par(mc_par, lbias=lbias, &
695 : BETA=BETA, exp_max_val=exp_max_val, &
696 : exp_min_val=exp_min_val, rmrot=rmrot, mc_molecule_info=mc_molecule_info, &
697 502 : ionode=ionode, group=group, source=source)
698 : CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits, &
699 : nunits_tot=nunits_tot, nchains=nchains, mass=mass, &
700 502 : mol_type=mol_type)
701 :
702 : ! figure out some bounds for mol_type
703 502 : start_mol = 1
704 516 : DO jbox = 1, box_number - 1
705 544 : start_mol = start_mol + SUM(nchains(:, jbox))
706 : END DO
707 1506 : end_mol = start_mol + SUM(nchains(:, box_number)) - 1
708 :
709 502 : nunits_mol = nunits(molecule_type)
710 :
711 : ! nullify some pointers
712 502 : NULLIFY (particles, subsys)
713 :
714 : ! do some allocation
715 1506 : ALLOCATE (r_old(1:3, 1:nunits_tot(box_number)))
716 :
717 : ! initialize some stuff
718 502 : lx = .FALSE.
719 502 : ly = .FALSE.
720 :
721 : ! determine what the final atom in the molecule is numbered, and which
722 : ! molecule number this is
723 502 : end_atom = start_atom + nunits_mol - 1
724 502 : molecule_number = 0
725 502 : atom_number = 1
726 2956 : DO imolecule = 1, SUM(nchains(:, box_number))
727 1952 : IF (atom_number == start_atom) THEN
728 502 : molecule_number = imolecule
729 502 : EXIT
730 : END IF
731 1450 : atom_number = atom_number + nunits(mol_type(imolecule + start_mol - 1))
732 : END DO
733 502 : IF (molecule_number == 0) CPABORT('Cannot find the molecule number')
734 :
735 : ! are we biasing this move?
736 502 : IF (lbias) THEN
737 :
738 : ! grab the coordinates
739 424 : CALL force_env_get(bias_env, subsys=subsys)
740 424 : CALL cp_subsys_get(subsys, particles=particles)
741 :
742 : ! save the coordinates
743 11808 : DO ipart = 1, nunits_tot(box_number)
744 45960 : r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
745 : END DO
746 :
747 : ! save the energy
748 424 : bias_energy_old = bias_energy
749 :
750 : ELSE
751 :
752 : ! grab the coordinates
753 78 : CALL force_env_get(force_env, subsys=subsys)
754 78 : CALL cp_subsys_get(subsys, particles=particles)
755 : END IF
756 :
757 : ! grab the masses
758 2008 : masstot = SUM(mass(1:nunits(molecule_type), molecule_type))
759 :
760 : ! record the attempt
761 502 : moves%bias_rot%attempts = moves%bias_rot%attempts + 1
762 502 : move_updates%bias_rot%attempts = move_updates%bias_rot%attempts + 1
763 502 : moves%rot%attempts = moves%rot%attempts + 1
764 502 : move_updates%rot%attempts = move_updates%rot%attempts + 1
765 502 : IF (.NOT. lbias) THEN
766 78 : moves%rot%qsuccesses = moves%rot%qsuccesses + 1
767 78 : move_updates%rot%qsuccesses = move_updates%rot%qsuccesses + 1
768 78 : moves%bias_rot%qsuccesses = moves%bias_rot%qsuccesses + 1
769 78 : move_updates%bias_rot%qsuccesses = move_updates%bias_rot%qsuccesses + 1
770 : END IF
771 :
772 : ! rotate one molecule in the system
773 :
774 : ! call a random number to figure out which direction we're moving
775 502 : IF (ionode) rand = rng_stream%next()
776 : ! CALL RANDOM_NUMBER(rand)
777 502 : CALL group%bcast(rand, source)
778 : ! 1,2,3 with equal prob
779 502 : dir = INT(3*rand) + 1
780 :
781 502 : IF (dir == 1) THEN
782 : lx = .TRUE.
783 334 : ELSE IF (dir == 2) THEN
784 176 : ly = .TRUE.
785 : END IF
786 :
787 : ! Determine new center of mass for chain i by finding the sum
788 : ! of m*r for each unit, then dividing by the total mass of the chain
789 502 : nxcm = 0.0E0_dp
790 502 : nycm = 0.0E0_dp
791 502 : nzcm = 0.0E0_dp
792 2008 : DO ii = 1, nunits_mol
793 1506 : nxcm = nxcm + particles%els(start_atom - 1 + ii)%r(1)*mass(ii, molecule_type)
794 1506 : nycm = nycm + particles%els(start_atom - 1 + ii)%r(2)*mass(ii, molecule_type)
795 2008 : nzcm = nzcm + particles%els(start_atom - 1 + ii)%r(3)*mass(ii, molecule_type)
796 : END DO
797 502 : nxcm = nxcm/masstot
798 502 : nycm = nycm/masstot
799 502 : nzcm = nzcm/masstot
800 :
801 : ! call a random number to figure out how far we're moving
802 502 : IF (ionode) rand = rng_stream%next()
803 502 : CALL group%bcast(rand, source)
804 502 : dgamma = rmrot(molecule_type)*(rand - 0.5E0_dp)*2.0E0_dp
805 :
806 : ! *** set up the rotation matrix ***
807 :
808 502 : cosdg = COS(dgamma)
809 502 : sindg = SIN(dgamma)
810 :
811 502 : IF (lx) THEN
812 :
813 : ! *** ROTATE UNITS OF I AROUND X-AXIS ***
814 :
815 672 : DO iunit = start_atom, end_atom
816 504 : ry = particles%els(iunit)%r(2) - nycm
817 504 : rz = particles%els(iunit)%r(3) - nzcm
818 504 : rynew = cosdg*ry - sindg*rz
819 504 : rznew = cosdg*rz + sindg*ry
820 :
821 504 : particles%els(iunit)%r(2) = rynew + nycm
822 672 : particles%els(iunit)%r(3) = rznew + nzcm
823 :
824 : END DO
825 334 : ELSE IF (ly) THEN
826 :
827 : ! *** ROTATE UNITS OF I AROUND y-AXIS ***
828 :
829 704 : DO iunit = start_atom, end_atom
830 528 : rx = particles%els(iunit)%r(1) - nxcm
831 528 : rz = particles%els(iunit)%r(3) - nzcm
832 528 : rxnew = cosdg*rx + sindg*rz
833 528 : rznew = cosdg*rz - sindg*rx
834 :
835 528 : particles%els(iunit)%r(1) = rxnew + nxcm
836 704 : particles%els(iunit)%r(3) = rznew + nzcm
837 :
838 : END DO
839 :
840 : ELSE
841 :
842 : ! *** ROTATE UNITS OF I AROUND z-AXIS ***
843 :
844 632 : DO iunit = start_atom, end_atom
845 474 : rx = particles%els(iunit)%r(1) - nxcm
846 474 : ry = particles%els(iunit)%r(2) - nycm
847 :
848 474 : rxnew = cosdg*rx - sindg*ry
849 474 : rynew = cosdg*ry + sindg*rx
850 :
851 474 : particles%els(iunit)%r(1) = rxnew + nxcm
852 632 : particles%els(iunit)%r(2) = rynew + nycm
853 :
854 : END DO
855 :
856 : END IF
857 502 : CALL cp_subsys_set(subsys, particles=particles)
858 :
859 : ! check for overlap
860 502 : lreject = .FALSE.
861 502 : IF (lbias) THEN
862 : CALL check_for_overlap(bias_env, nchains(:, box_number), &
863 : nunits(:), loverlap, mol_type(start_mol:end_mol), &
864 424 : molecule_number=molecule_number)
865 : ELSE
866 : CALL check_for_overlap(force_env, nchains(:, box_number), &
867 : nunits(:), loverlap, mol_type(start_mol:end_mol), &
868 78 : molecule_number=molecule_number)
869 78 : IF (loverlap) lreject = .TRUE.
870 : END IF
871 :
872 : ! if we're biasing classical, check for acceptance
873 502 : IF (lbias) THEN
874 :
875 : ! here's where we bias the moves
876 :
877 424 : IF (loverlap) THEN
878 : w = 0.0E0_dp
879 : ELSE
880 424 : CALL force_env_calc_energy_force(bias_env, calc_force=.FALSE.)
881 : CALL force_env_get(bias_env, &
882 424 : potential_energy=bias_energy_new)
883 : ! accept or reject the move based on the Metropolis rule
884 424 : value = -BETA*(bias_energy_new - bias_energy_old)
885 424 : IF (value > exp_max_val) THEN
886 : w = 10.0_dp
887 424 : ELSE IF (value < exp_min_val) THEN
888 : w = 0.0_dp
889 : ELSE
890 424 : w = EXP(value)
891 : END IF
892 :
893 : END IF
894 :
895 424 : IF (w >= 1.0E0_dp) THEN
896 180 : w = 1.0E0_dp
897 180 : rand = 0.0E0_dp
898 : ELSE
899 244 : IF (ionode) rand = rng_stream%next()
900 244 : CALL group%bcast(rand, source)
901 : END IF
902 :
903 424 : IF (rand < w) THEN
904 :
905 : ! accept the move
906 340 : moves%bias_rot%successes = moves%bias_rot%successes + 1
907 340 : move_updates%bias_rot%successes = move_updates%bias_rot%successes + 1
908 340 : moves%rot%qsuccesses = moves%rot%qsuccesses + 1
909 340 : move_updates%rot%successes = move_updates%rot%successes + 1
910 : bias_energy = bias_energy + bias_energy_new - &
911 340 : bias_energy_old
912 :
913 : ELSE
914 :
915 : ! reject the move
916 : ! restore the coordinates
917 84 : CALL force_env_get(bias_env, subsys=subsys)
918 84 : CALL cp_subsys_get(subsys, particles=particles)
919 2316 : DO ipart = 1, nunits_tot(box_number)
920 9012 : particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
921 : END DO
922 84 : CALL cp_subsys_set(subsys, particles=particles)
923 :
924 : END IF
925 :
926 : END IF
927 :
928 : ! deallocate some stuff
929 502 : DEALLOCATE (r_old)
930 :
931 : ! end the timing
932 502 : CALL timestop(handle)
933 :
934 502 : END SUBROUTINE mc_molecule_rotation
935 :
936 : ! **************************************************************************************************
937 : !> \brief performs a Monte Carlo move that alters the volume of the simulation box
938 : !> \param mc_par the mc parameters for the force env
939 : !> \param force_env the force environment whose cell we're changing
940 : !> \param moves the structure that keeps track of how many moves have been
941 : !> accepted/rejected
942 : !> \param move_updates the structure that keeps track of how many moves have
943 : !> been accepted/rejected since the last time the displacements
944 : !> were updated
945 : !> \param old_energy the energy of the last accepted move involving an
946 : !> unbiased calculation
947 : !> \param box_number the box we're changing the volume of
948 : !> \param energy_check the running total of how much the energy has changed
949 : !> since the initial configuration
950 : !> \param r_old the coordinates of the last accepted move involving an
951 : !> unbiased calculation
952 : !> \param iw the unit number that writes to the screen
953 : !> \param discrete_array tells use which volumes we can do for the discrete
954 : !> case
955 : !> \param rng_stream the random number stream that we draw from
956 : !> \author MJM
957 : !> \note Designed for parallel use.
958 : ! **************************************************************************************************
959 34 : SUBROUTINE mc_volume_move(mc_par, force_env, moves, move_updates, &
960 : old_energy, box_number, &
961 34 : energy_check, r_old, iw, discrete_array, rng_stream)
962 :
963 : TYPE(mc_simpar_type), POINTER :: mc_par
964 : TYPE(force_env_type), POINTER :: force_env
965 : TYPE(mc_moves_type), POINTER :: moves, move_updates
966 : REAL(KIND=dp), INTENT(INOUT) :: old_energy
967 : INTEGER, INTENT(IN) :: box_number
968 : REAL(KIND=dp), INTENT(INOUT) :: energy_check
969 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r_old
970 : INTEGER, INTENT(IN) :: iw
971 : INTEGER, DIMENSION(1:3, 1:2), INTENT(INOUT) :: discrete_array
972 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
973 :
974 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_volume_move'
975 :
976 : CHARACTER(LEN=200) :: fft_lib
977 : CHARACTER(LEN=40) :: dat_file
978 : INTEGER :: cl, end_atom, end_mol, handle, iatom, idim, imolecule, iside, iside_change, &
979 : iunit, jbox, nunits_mol, output_unit, print_level, source, start_atom, start_mol
980 34 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
981 34 : INTEGER, DIMENSION(:, :), POINTER :: nchains
982 : LOGICAL :: ionode, ldiscrete, lincrease, loverlap, &
983 : ltoo_small
984 34 : REAL(dp), DIMENSION(:, :), POINTER :: mass
985 : REAL(KIND=dp) :: BETA, discrete_step, energy_term, exp_max_val, exp_min_val, new_energy, &
986 : pressure, pressure_term, rand, rcut, rmvolume, temp_var, value, vol_dis, volume_term, w
987 34 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r
988 : REAL(KIND=dp), DIMENSION(1:3) :: abc, center_of_mass, center_of_mass_new, &
989 : diff, new_cell_length, old_cell_length
990 : REAL(KIND=dp), DIMENSION(1:3, 1:3) :: hmat_test
991 : TYPE(cell_type), POINTER :: cell, cell_old, cell_test
992 : TYPE(cp_logger_type), POINTER :: logger
993 : TYPE(cp_subsys_type), POINTER :: oldsys
994 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
995 : TYPE(mp_comm_type) :: group
996 : TYPE(particle_list_type), POINTER :: particles_old
997 :
998 : ! begin the timing of the subroutine
999 :
1000 34 : CALL timeset(routineN, handle)
1001 :
1002 : ! get a bunch of stuff from mc_par
1003 : CALL get_mc_par(mc_par, ionode=ionode, &
1004 : BETA=BETA, exp_max_val=exp_max_val, &
1005 : exp_min_val=exp_min_val, source=source, group=group, &
1006 : dat_file=dat_file, rmvolume=rmvolume, pressure=pressure, cl=cl, &
1007 : fft_lib=fft_lib, discrete_step=discrete_step, &
1008 34 : ldiscrete=ldiscrete, mc_molecule_info=mc_molecule_info)
1009 : CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, &
1010 : nunits=nunits, nunits_tot=nunits_tot, mol_type=mol_type, &
1011 34 : mass=mass)
1012 : ! figure out some bounds for mol_type
1013 34 : start_mol = 1
1014 46 : DO jbox = 1, box_number - 1
1015 70 : start_mol = start_mol + SUM(nchains(:, jbox))
1016 : END DO
1017 94 : end_mol = start_mol + SUM(nchains(:, box_number)) - 1
1018 :
1019 34 : print_level = 1 ! hack, printlevel is for print_keys
1020 :
1021 : ! nullify some pointers
1022 34 : NULLIFY (particles_old, cell_old, oldsys, cell_test, cell)
1023 :
1024 : ! do some allocation
1025 102 : ALLOCATE (r(1:3, 1:nunits_tot(box_number)))
1026 :
1027 : ! record the attempt
1028 34 : moves%volume%attempts = moves%volume%attempts + 1
1029 34 : move_updates%volume%attempts = move_updates%volume%attempts + 1
1030 :
1031 : ! now let's grab the cell length and particle positions
1032 34 : CALL force_env_get(force_env, subsys=oldsys, cell=cell)
1033 34 : CALL get_cell(cell, abc=abc)
1034 34 : CALL cell_create(cell_old)
1035 34 : CALL cell_clone(cell, cell_old, tag="CELL_OLD")
1036 34 : CALL cp_subsys_get(oldsys, particles=particles_old)
1037 :
1038 : ! find the old cell length
1039 34 : old_cell_length(1) = abc(1)
1040 34 : old_cell_length(2) = abc(2)
1041 34 : old_cell_length(3) = abc(3)
1042 :
1043 : ! save the old coordinates
1044 760 : DO iatom = 1, nunits_tot(box_number)
1045 2938 : r(1:3, iatom) = particles_old%els(iatom)%r(1:3)
1046 : END DO
1047 :
1048 : ! now do the move
1049 :
1050 : ! call a random number to figure out how far we're moving
1051 34 : IF (ionode) rand = rng_stream%next()
1052 34 : CALL group%bcast(rand, source)
1053 :
1054 : ! find the test cell lengths for the discrete volume move
1055 34 : IF (ldiscrete) THEN
1056 0 : IF (rand < 0.5_dp) THEN
1057 : lincrease = .TRUE.
1058 : ELSE
1059 0 : lincrease = .FALSE.
1060 : END IF
1061 :
1062 0 : new_cell_length(1:3) = old_cell_length(1:3)
1063 :
1064 : ! if we're increasing the volume, we need to find a side we can increase
1065 0 : IF (lincrease) THEN
1066 : DO
1067 0 : IF (ionode) rand = rng_stream%next()
1068 0 : CALL group%bcast(rand, source)
1069 0 : iside_change = CEILING(3.0_dp*rand)
1070 0 : IF (discrete_array(iside_change, 1) == 1) THEN
1071 : new_cell_length(iside_change) = &
1072 0 : new_cell_length(iside_change) + discrete_step
1073 : EXIT
1074 : END IF
1075 : END DO
1076 : ELSE
1077 : DO
1078 0 : IF (ionode) rand = rng_stream%next()
1079 0 : CALL group%bcast(rand, source)
1080 0 : iside_change = CEILING(3.0_dp*rand)
1081 0 : IF (discrete_array(iside_change, 2) == 1) THEN
1082 : new_cell_length(iside_change) = &
1083 0 : new_cell_length(iside_change) - discrete_step
1084 0 : EXIT
1085 : END IF
1086 : END DO
1087 : END IF
1088 : vol_dis = (new_cell_length(1)*new_cell_length(2)*new_cell_length(3)) &
1089 0 : - old_cell_length(1)*old_cell_length(2)*old_cell_length(3)
1090 : ELSE
1091 : ! now for the not discrete volume move
1092 : !!!!!!!!!!!!!!!! for E_V curves
1093 34 : vol_dis = rmvolume*(rand - 0.5E0_dp)*2.0E0_dp
1094 : ! WRITE(output_unit,*) '************************ be sure to change back!',&
1095 : ! old_cell_length(1),14.64_dp/angstrom
1096 : ! vol_dis=-56.423592_dp/angstrom**3
1097 : ! IF(old_cell_length(1) <= 14.64_dp/angstrom) THEN
1098 : ! vol_dis=0.0_dp
1099 : ! WRITE(output_unit,*) 'Found the correct box length!'
1100 : ! ENDIF
1101 :
1102 : temp_var = vol_dis + &
1103 : old_cell_length(1)*old_cell_length(2)* &
1104 34 : old_cell_length(3)
1105 :
1106 34 : IF (temp_var <= 0.0E0_dp) THEN
1107 0 : loverlap = .TRUE. ! cannot have a negative volume
1108 : ELSE
1109 34 : new_cell_length(1) = (temp_var)**(1.0E0_dp/3.0E0_dp)
1110 34 : new_cell_length(2) = new_cell_length(1)
1111 34 : new_cell_length(3) = new_cell_length(1)
1112 34 : loverlap = .FALSE.
1113 : END IF
1114 : END IF
1115 34 : CALL group%bcast(loverlap, source)
1116 :
1117 34 : IF (loverlap) THEN
1118 : ! deallocate some stuff
1119 0 : DEALLOCATE (r)
1120 0 : logger => cp_get_default_logger()
1121 0 : output_unit = cp_logger_get_default_io_unit(logger)
1122 0 : IF (output_unit > 0) WRITE (output_unit, *) &
1123 0 : "Volume move rejected because we tried to make too small of box.", vol_dis
1124 : ! end the timing
1125 0 : CALL timestop(handle)
1126 0 : RETURN
1127 : END IF
1128 :
1129 : ! now we need to make the new cell
1130 34 : hmat_test(:, :) = 0.0e0_dp
1131 34 : hmat_test(1, 1) = new_cell_length(1)
1132 34 : hmat_test(2, 2) = new_cell_length(2)
1133 34 : hmat_test(3, 3) = new_cell_length(3)
1134 34 : CALL cell_create(cell_test, hmat=hmat_test(:, :), periodic=cell%perd)
1135 34 : CALL cp_subsys_set(oldsys, cell=cell_test)
1136 :
1137 : ! now we need to scale the coordinates of all the molecules by the
1138 : ! center of mass, using the minimum image (not all molecules are in
1139 : ! the central box)
1140 :
1141 : ! now we need to scale the coordinates of all the molecules by the
1142 : ! center of mass
1143 34 : end_atom = 0
1144 432 : DO imolecule = 1, SUM(nchains(:, box_number))
1145 338 : nunits_mol = nunits(mol_type(imolecule + start_mol - 1))
1146 338 : start_atom = end_atom + 1
1147 338 : end_atom = start_atom + nunits_mol - 1
1148 : ! now find the center of mass
1149 : CALL get_center_of_mass(r(:, start_atom:end_atom), nunits_mol, &
1150 338 : center_of_mass(:), mass(:, mol_type(imolecule + start_mol - 1)))
1151 :
1152 : ! scale the center of mass and determine the vector that points from the
1153 : ! old COM to the new one
1154 1352 : DO iside = 1, 3
1155 : center_of_mass_new(iside) = center_of_mass(iside)* &
1156 1352 : new_cell_length(iside)/old_cell_length(iside)
1157 : END DO
1158 :
1159 1386 : DO idim = 1, 3
1160 1014 : diff(idim) = center_of_mass_new(idim) - center_of_mass(idim)
1161 : ! now change the particle positions
1162 3530 : DO iunit = start_atom, end_atom
1163 : particles_old%els(iunit)%r(idim) = &
1164 3192 : particles_old%els(iunit)%r(idim) + diff(idim)
1165 : END DO
1166 : END DO
1167 : END DO
1168 :
1169 : ! check for overlap
1170 : CALL check_for_overlap(force_env, nchains(:, box_number), &
1171 : nunits(:), loverlap, mol_type(start_mol:end_mol), &
1172 34 : cell_length=new_cell_length)
1173 :
1174 : ! figure out if we have overlap problems
1175 34 : CALL group%bcast(loverlap, source)
1176 34 : IF (loverlap) THEN
1177 : ! deallocate some stuff
1178 0 : DEALLOCATE (r)
1179 :
1180 0 : logger => cp_get_default_logger()
1181 0 : output_unit = cp_logger_get_default_io_unit(logger)
1182 0 : IF (output_unit > 0) WRITE (output_unit, *) &
1183 0 : "Volume move rejected due to overlap.", vol_dis
1184 : ! end the timing
1185 0 : CALL timestop(handle)
1186 : ! reset the cell and particle positions
1187 0 : CALL cp_subsys_set(oldsys, cell=cell_old)
1188 0 : DO iatom = 1, nunits_tot(box_number)
1189 0 : particles_old%els(iatom)%r(1:3) = r_old(1:3, iatom)
1190 : END DO
1191 : RETURN
1192 : END IF
1193 :
1194 : ! stop if we're trying to change a box to a boxlength smaller than rcut
1195 34 : IF (ionode) THEN
1196 17 : ltoo_small = .FALSE.
1197 17 : IF (force_env%in_use == use_fist_force) THEN
1198 13 : CALL get_mc_par(mc_par, rcut=rcut)
1199 13 : IF (new_cell_length(1) < 2.0_dp*rcut) ltoo_small = .TRUE.
1200 13 : IF (new_cell_length(2) < 2.0_dp*rcut) ltoo_small = .TRUE.
1201 13 : IF (new_cell_length(3) < 2.0_dp*rcut) ltoo_small = .TRUE.
1202 :
1203 13 : IF (ltoo_small) THEN
1204 0 : WRITE (iw, *) 'new_cell_lengths ', &
1205 0 : new_cell_length(1:3)/angstrom
1206 0 : WRITE (iw, *) 'rcut ', rcut/angstrom
1207 : END IF
1208 : END IF
1209 : END IF
1210 34 : CALL group%bcast(ltoo_small, source)
1211 34 : IF (ltoo_small) THEN
1212 0 : CPABORT("Attempted a volume move where box size got too small.")
1213 : END IF
1214 :
1215 : ! now compute the energy
1216 34 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
1217 : CALL force_env_get(force_env, &
1218 34 : potential_energy=new_energy)
1219 :
1220 : ! accept or reject the move
1221 : ! to prevent overflows
1222 34 : energy_term = new_energy - old_energy
1223 : volume_term = -REAL(SUM(nchains(:, box_number)), dp)/BETA* &
1224 : LOG(new_cell_length(1)*new_cell_length(2)*new_cell_length(3)/ &
1225 94 : (old_cell_length(1)*old_cell_length(2)*old_cell_length(3)))
1226 34 : pressure_term = pressure*vol_dis
1227 :
1228 34 : value = -BETA*(energy_term + volume_term + pressure_term)
1229 34 : IF (value > exp_max_val) THEN
1230 : w = 10.0_dp
1231 34 : ELSE IF (value < exp_min_val) THEN
1232 : w = 0.0_dp
1233 : ELSE
1234 34 : w = EXP(value)
1235 : END IF
1236 :
1237 : !!!!!!!!!!!!!!!! for E_V curves
1238 : ! w=1.0E0_dp
1239 : ! w=0.0E0_dp
1240 :
1241 34 : IF (w >= 1.0E0_dp) THEN
1242 18 : w = 1.0E0_dp
1243 18 : rand = 0.0E0_dp
1244 : ELSE
1245 16 : IF (ionode) rand = rng_stream%next()
1246 16 : CALL group%bcast(rand, source)
1247 : END IF
1248 :
1249 34 : IF (rand < w) THEN
1250 :
1251 : ! accept the move
1252 30 : moves%volume%successes = moves%volume%successes + 1
1253 30 : move_updates%volume%successes = move_updates%volume%successes + 1
1254 :
1255 : ! update energies
1256 30 : energy_check = energy_check + (new_energy - old_energy)
1257 30 : old_energy = new_energy
1258 :
1259 720 : DO iatom = 1, nunits_tot(box_number)
1260 2790 : r_old(1:3, iatom) = particles_old%els(iatom)%r(1:3)
1261 : END DO
1262 :
1263 : ! update discrete_array if we're doing a discrete volume move
1264 30 : IF (ldiscrete) THEN
1265 : CALL create_discrete_array(new_cell_length(:), &
1266 0 : discrete_array(:, :), discrete_step)
1267 : END IF
1268 :
1269 : ELSE
1270 :
1271 : ! reset the cell and particle positions
1272 4 : CALL cp_subsys_set(oldsys, cell=cell_old)
1273 40 : DO iatom = 1, nunits_tot(box_number)
1274 148 : particles_old%els(iatom)%r(1:3) = r_old(1:3, iatom)
1275 : END DO
1276 :
1277 : END IF
1278 :
1279 : ! deallocate some stuff
1280 34 : DEALLOCATE (r)
1281 34 : CALL cell_release(cell_test)
1282 34 : CALL cell_release(cell_old)
1283 :
1284 : ! end the timing
1285 34 : CALL timestop(handle)
1286 :
1287 68 : END SUBROUTINE mc_volume_move
1288 :
1289 : ! **************************************************************************************************
1290 : !> \brief alters the length of a random bond for the given molecule, using
1291 : !> a mass weighted scheme so the lightest atoms move the most
1292 : !> \param r_old the initial coordinates of all molecules in the system
1293 : !> \param r_new the new coordinates of all molecules in the system
1294 : !> \param mc_par the mc parameters for the force env
1295 : !> \param molecule_type the molecule type that we're moving
1296 : !> \param molecule_kind the structure containing the molecule information
1297 : !> \param dis_length the ratio of the new bond length to the old bond length,
1298 : !> used in the acceptance rule
1299 : !> \param particles the particle_list_type for all particles in the force_env..
1300 : !> used to grab the mass of each atom
1301 : !> \param rng_stream the random number stream that we draw from
1302 : !>
1303 : !> This subroutine is written to be parallel.
1304 : !> \author MJM
1305 : ! **************************************************************************************************
1306 312 : SUBROUTINE change_bond_length(r_old, r_new, mc_par, molecule_type, molecule_kind, &
1307 : dis_length, particles, rng_stream)
1308 :
1309 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_old
1310 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: r_new
1311 : TYPE(mc_simpar_type), POINTER :: mc_par
1312 : INTEGER, INTENT(IN) :: molecule_type
1313 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1314 : REAL(KIND=dp), INTENT(OUT) :: dis_length
1315 : TYPE(particle_list_type), POINTER :: particles
1316 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1317 :
1318 : CHARACTER(len=*), PARAMETER :: routineN = 'change_bond_length'
1319 :
1320 : INTEGER :: bond_number, handle, i, iatom, ibond, &
1321 : ipart, natom, nbond, source
1322 312 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_b, counter
1323 312 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity
1324 312 : INTEGER, DIMENSION(:), POINTER :: nunits
1325 : LOGICAL :: ionode
1326 312 : REAL(dp), DIMENSION(:), POINTER :: rmbond
1327 : REAL(KIND=dp) :: atom_mass, mass_a, mass_b, new_length_a, &
1328 : new_length_b, old_length, rand
1329 : REAL(KIND=dp), DIMENSION(1:3) :: bond_a, bond_b
1330 312 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
1331 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1332 : TYPE(mp_comm_type) :: group
1333 :
1334 : ! begin the timing of the subroutine
1335 :
1336 312 : CALL timeset(routineN, handle)
1337 :
1338 312 : NULLIFY (rmbond, mc_molecule_info)
1339 :
1340 : ! get some stuff from mc_par
1341 : CALL get_mc_par(mc_par, mc_molecule_info=mc_molecule_info, source=source, &
1342 312 : group=group, rmbond=rmbond, ionode=ionode)
1343 312 : CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
1344 :
1345 : ! copy the incoming coordinates so we can change them
1346 1248 : DO ipart = 1, nunits(molecule_type)
1347 4056 : r_new(1:3, ipart) = r_old(1:3, ipart)
1348 : END DO
1349 :
1350 : ! pick which bond in the molecule at random
1351 312 : IF (ionode) THEN
1352 156 : rand = rng_stream%next()
1353 : END IF
1354 312 : CALL group%bcast(rand, source)
1355 : CALL get_molecule_kind(molecule_kind, natom=natom, nbond=nbond, &
1356 312 : bond_list=bond_list)
1357 312 : bond_number = CEILING(rand*REAL(nbond, dp))
1358 :
1359 936 : ALLOCATE (connection(1:natom, 1:2))
1360 : ! assume at most six bonds per atom
1361 936 : ALLOCATE (connectivity(1:6, 1:natom))
1362 936 : ALLOCATE (counter(1:natom))
1363 624 : ALLOCATE (atom_a(1:natom))
1364 624 : ALLOCATE (atom_b(1:natom))
1365 312 : connection(:, :) = 0
1366 312 : connectivity(:, :) = 0
1367 312 : counter(:) = 0
1368 312 : atom_a(:) = 0
1369 312 : atom_b(:) = 0
1370 :
1371 : ! now we need to find a list of atoms that each atom in this bond is connected
1372 : ! to
1373 1248 : DO iatom = 1, natom
1374 3120 : DO ibond = 1, nbond
1375 2808 : IF (bond_list(ibond)%a == iatom) THEN
1376 624 : counter(iatom) = counter(iatom) + 1
1377 624 : connectivity(counter(iatom), iatom) = bond_list(ibond)%b
1378 1248 : ELSE IF (bond_list(ibond)%b == iatom) THEN
1379 624 : counter(iatom) = counter(iatom) + 1
1380 624 : connectivity(counter(iatom), iatom) = bond_list(ibond)%a
1381 : END IF
1382 : END DO
1383 : END DO
1384 :
1385 : ! now I need to do a depth first search to figure out which atoms are on atom a's
1386 : ! side and which are on atom b's
1387 312 : atom_a(:) = 0
1388 312 : atom_a(bond_list(bond_number)%a) = 1
1389 : CALL depth_first_search(bond_list(bond_number)%a, bond_list(bond_number)%b, &
1390 312 : connectivity(:, :), atom_a(:))
1391 312 : atom_b(:) = 0
1392 312 : atom_b(bond_list(bond_number)%b) = 1
1393 : CALL depth_first_search(bond_list(bond_number)%b, bond_list(bond_number)%a, &
1394 312 : connectivity(:, :), atom_b(:))
1395 :
1396 : ! now figure out the masses of the various sides, so we can weight how far we move each
1397 : ! group of atoms
1398 312 : mass_a = 0.0_dp
1399 312 : mass_b = 0.0_dp
1400 1248 : DO iatom = 1, natom
1401 : CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
1402 936 : mass=atom_mass)
1403 1248 : IF (atom_a(iatom) == 1) THEN
1404 624 : mass_a = mass_a + atom_mass
1405 : ELSE
1406 312 : mass_b = mass_b + atom_mass
1407 : END IF
1408 : END DO
1409 :
1410 : ! choose a displacement
1411 312 : IF (ionode) rand = rng_stream%next()
1412 312 : CALL group%bcast(rand, source)
1413 :
1414 312 : dis_length = rmbond(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
1415 :
1416 : ! find the bond vector that atom a will be moving
1417 1248 : DO i = 1, 3
1418 : bond_a(i) = r_new(i, bond_list(bond_number)%a) - &
1419 936 : r_new(i, bond_list(bond_number)%b)
1420 1248 : bond_b(i) = -bond_a(i)
1421 : END DO
1422 :
1423 : ! notice we weight by the opposite masses...therefore lighter segments
1424 : ! will move further
1425 1248 : old_length = NORM2(bond_a)
1426 312 : new_length_a = dis_length*mass_b/(mass_a + mass_b)
1427 312 : new_length_b = dis_length*mass_a/(mass_a + mass_b)
1428 :
1429 1248 : DO i = 1, 3
1430 936 : bond_a(i) = bond_a(i)/old_length*new_length_a
1431 1248 : bond_b(i) = bond_b(i)/old_length*new_length_b
1432 : END DO
1433 :
1434 1248 : DO iatom = 1, natom
1435 1248 : IF (atom_a(iatom) == 1) THEN
1436 624 : r_new(1, iatom) = r_new(1, iatom) + bond_a(1)
1437 624 : r_new(2, iatom) = r_new(2, iatom) + bond_a(2)
1438 624 : r_new(3, iatom) = r_new(3, iatom) + bond_a(3)
1439 : ELSE
1440 312 : r_new(1, iatom) = r_new(1, iatom) + bond_b(1)
1441 312 : r_new(2, iatom) = r_new(2, iatom) + bond_b(2)
1442 312 : r_new(3, iatom) = r_new(3, iatom) + bond_b(3)
1443 : END IF
1444 : END DO
1445 :
1446 : ! correct the value of dis_length for the acceptance rule
1447 312 : dis_length = (old_length + dis_length)/old_length
1448 :
1449 312 : DEALLOCATE (connection)
1450 312 : DEALLOCATE (connectivity)
1451 312 : DEALLOCATE (counter)
1452 312 : DEALLOCATE (atom_a)
1453 312 : DEALLOCATE (atom_b)
1454 : ! end the timing
1455 312 : CALL timestop(handle)
1456 :
1457 624 : END SUBROUTINE change_bond_length
1458 :
1459 : ! **************************************************************************************************
1460 : !> \brief Alters the magnitude of a random angle in a molecule centered on atom C
1461 : !> (connected to atoms A and B). Atoms A and B are moved amounts related
1462 : !> to their masses (and masses of all connecting atoms), so that heavier
1463 : !> segments are moved less.
1464 : !> \param r_old the initial coordinates of all molecules in the system
1465 : !> \param r_new the new coordinates of all molecules in the system
1466 : !> \param mc_par the mc parameters for the force env
1467 : !> \param molecule_type the type of molecule we're playing with
1468 : !> \param molecule_kind the structure containing the molecule information
1469 : !> \param particles the particle_list_type for all particles in the force_env...
1470 : !> used to grab the mass of each atom
1471 : !> \param rng_stream the random number stream that we draw from
1472 : !> \author MJM
1473 : ! **************************************************************************************************
1474 194 : SUBROUTINE change_bond_angle(r_old, r_new, mc_par, molecule_type, molecule_kind, &
1475 : particles, rng_stream)
1476 :
1477 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_old
1478 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: r_new
1479 : TYPE(mc_simpar_type), POINTER :: mc_par
1480 : INTEGER, INTENT(IN) :: molecule_type
1481 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1482 : TYPE(particle_list_type), POINTER :: particles
1483 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1484 :
1485 : CHARACTER(len=*), PARAMETER :: routineN = 'change_bond_angle'
1486 :
1487 : INTEGER :: bend_number, handle, i, iatom, ibond, &
1488 : ipart, natom, nbend, nbond, source
1489 194 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_c, counter
1490 194 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity
1491 194 : INTEGER, DIMENSION(:), POINTER :: nunits
1492 : LOGICAL :: ionode
1493 194 : REAL(dp), DIMENSION(:), POINTER :: rmangle
1494 : REAL(KIND=dp) :: atom_mass, bis_length, dis_angle, dis_angle_a, dis_angle_c, mass_a, mass_c, &
1495 : new_angle_a, new_angle_c, old_angle, old_length_a, old_length_c, rand, temp_length
1496 : REAL(KIND=dp), DIMENSION(1:3) :: bisector, bond_a, bond_c, cross_prod, &
1497 : cross_prod_plane, temp
1498 194 : TYPE(bend_type), DIMENSION(:), POINTER :: bend_list
1499 194 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
1500 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1501 : TYPE(mp_comm_type) :: group
1502 :
1503 : ! begin the timing of the subroutine
1504 :
1505 194 : CALL timeset(routineN, handle)
1506 :
1507 194 : NULLIFY (bend_list, bond_list, rmangle, mc_molecule_info)
1508 :
1509 : ! get some stuff from mc_par
1510 : CALL get_mc_par(mc_par, rmangle=rmangle, source=source, &
1511 194 : group=group, ionode=ionode, mc_molecule_info=mc_molecule_info)
1512 194 : CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
1513 :
1514 : ! copy the incoming coordinates so we can change them
1515 776 : DO ipart = 1, nunits(molecule_type)
1516 2522 : r_new(1:3, ipart) = r_old(1:3, ipart)
1517 : END DO
1518 :
1519 : ! pick which bond in the molecule at random
1520 194 : IF (ionode) THEN
1521 97 : rand = rng_stream%next()
1522 : END IF
1523 194 : CALL group%bcast(rand, source)
1524 : CALL get_molecule_kind(molecule_kind, natom=natom, nbend=nbend, &
1525 194 : bend_list=bend_list, bond_list=bond_list, nbond=nbond)
1526 194 : bend_number = CEILING(rand*REAL(nbend, dp))
1527 :
1528 582 : ALLOCATE (connection(1:natom, 1:2))
1529 : ! assume at most six bonds per atom
1530 582 : ALLOCATE (connectivity(1:6, 1:natom))
1531 582 : ALLOCATE (counter(1:natom))
1532 388 : ALLOCATE (atom_a(1:natom))
1533 388 : ALLOCATE (atom_c(1:natom))
1534 194 : connection(:, :) = 0
1535 194 : connectivity(:, :) = 0
1536 194 : counter(:) = 0
1537 194 : atom_a(:) = 0
1538 194 : atom_c(:) = 0
1539 :
1540 : ! now we need to find a list of atoms that each atom in this bond is connected
1541 : ! to
1542 776 : DO iatom = 1, natom
1543 1940 : DO ibond = 1, nbond
1544 1746 : IF (bond_list(ibond)%a == iatom) THEN
1545 388 : counter(iatom) = counter(iatom) + 1
1546 388 : connectivity(counter(iatom), iatom) = bond_list(ibond)%b
1547 776 : ELSE IF (bond_list(ibond)%b == iatom) THEN
1548 388 : counter(iatom) = counter(iatom) + 1
1549 388 : connectivity(counter(iatom), iatom) = bond_list(ibond)%a
1550 : END IF
1551 : END DO
1552 : END DO
1553 :
1554 : ! now I need to do a depth first search to figure out which atoms are on atom a's
1555 : ! side and which are on atom c's
1556 194 : atom_a(:) = 0
1557 194 : atom_a(bend_list(bend_number)%a) = 1
1558 : CALL depth_first_search(bend_list(bend_number)%a, bend_list(bend_number)%b, &
1559 194 : connectivity(:, :), atom_a(:))
1560 194 : atom_c(:) = 0
1561 194 : atom_c(bend_list(bend_number)%c) = 1
1562 : CALL depth_first_search(bend_list(bend_number)%c, bend_list(bend_number)%b, &
1563 194 : connectivity(:, :), atom_c(:))
1564 :
1565 : ! now figure out the masses of the various sides, so we can weight how far we move each
1566 : ! group of atoms
1567 194 : mass_a = 0.0_dp
1568 194 : mass_c = 0.0_dp
1569 776 : DO iatom = 1, natom
1570 : CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
1571 582 : mass=atom_mass)
1572 582 : IF (atom_a(iatom) == 1) mass_a = mass_a + atom_mass
1573 1358 : IF (atom_c(iatom) == 1) mass_c = mass_c + atom_mass
1574 : END DO
1575 :
1576 : ! choose a displacement
1577 194 : IF (ionode) rand = rng_stream%next()
1578 194 : CALL group%bcast(rand, source)
1579 :
1580 194 : dis_angle = rmangle(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
1581 :
1582 : ! need to find the A-B-C bisector
1583 :
1584 : ! this going to be tough...we need to find the plane of the A-B-C bond and only shift
1585 : ! that component for all atoms connected to A and C...otherwise we change other
1586 : ! internal degrees of freedom
1587 :
1588 : ! find the bond vectors
1589 776 : DO i = 1, 3
1590 : bond_a(i) = r_new(i, bend_list(bend_number)%a) - &
1591 582 : r_new(i, bend_list(bend_number)%b)
1592 : bond_c(i) = r_new(i, bend_list(bend_number)%c) - &
1593 776 : r_new(i, bend_list(bend_number)%b)
1594 : END DO
1595 776 : old_length_a = NORM2(bond_a)
1596 776 : old_length_c = NORM2(bond_c)
1597 776 : old_angle = ACOS(DOT_PRODUCT(bond_a, bond_c)/(old_length_a*old_length_c))
1598 :
1599 776 : DO i = 1, 3
1600 : bisector(i) = bond_a(i)/old_length_a + & ! not yet normalized
1601 776 : bond_c(i)/old_length_c
1602 : END DO
1603 776 : bis_length = NORM2(bisector)
1604 776 : bisector(1:3) = bisector(1:3)/bis_length
1605 :
1606 : ! now we need to find the cross product of the B-A and B-C vectors and normalize
1607 : ! it, so we have a vector that defines the bend plane
1608 194 : cross_prod(1) = bond_a(2)*bond_c(3) - bond_a(3)*bond_c(2)
1609 194 : cross_prod(2) = bond_a(3)*bond_c(1) - bond_a(1)*bond_c(3)
1610 194 : cross_prod(3) = bond_a(1)*bond_c(2) - bond_a(2)*bond_c(1)
1611 1358 : cross_prod(1:3) = cross_prod(1:3)/NORM2(cross_prod)
1612 :
1613 : ! we have two axis of a coordinate system...let's get the third
1614 194 : cross_prod_plane(1) = cross_prod(2)*bisector(3) - cross_prod(3)*bisector(2)
1615 194 : cross_prod_plane(2) = cross_prod(3)*bisector(1) - cross_prod(1)*bisector(3)
1616 194 : cross_prod_plane(3) = cross_prod(1)*bisector(2) - cross_prod(2)*bisector(1)
1617 : cross_prod_plane(1:3) = cross_prod_plane(1:3)/ &
1618 1358 : NORM2(cross_prod_plane)
1619 :
1620 : ! now bisector is x, cross_prod_plane is the y vector (pointing towards c),
1621 : ! and cross_prod is z
1622 : ! shift the molecule so that atom b is at the origin
1623 776 : DO iatom = 1, natom
1624 : r_new(1:3, iatom) = r_new(1:3, iatom) - &
1625 2522 : r_old(1:3, bend_list(bend_number)%b)
1626 : END DO
1627 :
1628 : ! figure out how much we move each side, since we're mass-weighting, by the
1629 : ! opposite masses, so lighter moves farther..this angle is the angle between
1630 : ! the bond vector BA or BC and the bisector
1631 194 : dis_angle_a = dis_angle*mass_c/(mass_a + mass_c)
1632 194 : dis_angle_c = dis_angle*mass_a/(mass_a + mass_c)
1633 :
1634 : ! now loop through all the atoms, moving the ones that are connected to a or c
1635 776 : DO iatom = 1, natom
1636 : ! subtract out the z component (perpendicular to the angle plane)
1637 : temp(1:3) = r_new(1:3, iatom) - &
1638 : DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
1639 4074 : cross_prod(1:3)
1640 2328 : temp_length = NORM2(temp)
1641 :
1642 : ! we can now compute all three components of the new bond vector along the
1643 : ! axis defined above
1644 776 : IF (atom_a(iatom) == 1) THEN
1645 :
1646 : ! if the y-coordinate is less than zero, we need to switch the sign when we make the vector,
1647 : ! as the angle computed by the dot product can't distinguish between that
1648 776 : IF (DOT_PRODUCT(cross_prod_plane(1:3), r_new(1:3, iatom)) &
1649 : < 0.0_dp) THEN
1650 :
1651 : ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is
1652 : new_angle_a = ACOS(DOT_PRODUCT(bisector, temp(1:3))/ &
1653 776 : (temp_length)) + dis_angle_a
1654 :
1655 : r_new(1:3, iatom) = COS(new_angle_a)*temp_length*bisector(1:3) - &
1656 : SIN(new_angle_a)*temp_length*cross_prod_plane(1:3) + &
1657 : DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
1658 1358 : cross_prod(1:3)
1659 : ELSE
1660 :
1661 : ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is
1662 : new_angle_a = ACOS(DOT_PRODUCT(bisector, temp(1:3))/ &
1663 0 : (temp_length)) - dis_angle_a
1664 :
1665 : r_new(1:3, iatom) = COS(new_angle_a)*temp_length*bisector(1:3) + &
1666 : SIN(new_angle_a)*temp_length*cross_prod_plane(1:3) + &
1667 : DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
1668 0 : cross_prod(1:3)
1669 : END IF
1670 :
1671 388 : ELSE IF (atom_c(iatom) == 1) THEN
1672 :
1673 : ! if the y-coordinate is less than zero, we need to switch the sign when we make the vector,
1674 : ! as the angle computed by the dot product can't distinguish between that
1675 776 : IF (DOT_PRODUCT(cross_prod_plane(1:3), r_new(1:3, iatom)) &
1676 : < 0.0_dp) THEN
1677 : ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is
1678 : new_angle_c = ACOS(DOT_PRODUCT(bisector(1:3), temp(1:3))/ &
1679 0 : (temp_length)) - dis_angle_c
1680 :
1681 : r_new(1:3, iatom) = COS(new_angle_c)*temp_length*bisector(1:3) - &
1682 : SIN(new_angle_c)*temp_length*cross_prod_plane(1:3) + &
1683 : DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
1684 0 : cross_prod(1:3)
1685 : ELSE
1686 : new_angle_c = ACOS(DOT_PRODUCT(bisector(1:3), temp(1:3))/ &
1687 776 : (temp_length)) + dis_angle_c
1688 :
1689 : r_new(1:3, iatom) = COS(new_angle_c)*temp_length*bisector(1:3) + &
1690 : SIN(new_angle_c)*temp_length*cross_prod_plane(1:3) + &
1691 : DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
1692 1358 : cross_prod(1:3)
1693 : END IF
1694 : END IF
1695 :
1696 : END DO
1697 :
1698 776 : DO iatom = 1, natom
1699 : r_new(1:3, iatom) = r_new(1:3, iatom) + &
1700 2522 : r_old(1:3, bend_list(bend_number)%b)
1701 : END DO
1702 :
1703 : ! deallocate some stuff
1704 194 : DEALLOCATE (connection)
1705 194 : DEALLOCATE (connectivity)
1706 194 : DEALLOCATE (counter)
1707 194 : DEALLOCATE (atom_a)
1708 194 : DEALLOCATE (atom_c)
1709 :
1710 : ! end the timing
1711 194 : CALL timestop(handle)
1712 :
1713 388 : END SUBROUTINE change_bond_angle
1714 :
1715 : ! **************************************************************************************************
1716 : !> \brief Alters a dihedral (A-B-C-D) in the molecule so that all other internal
1717 : !> degrees of freedom remain the same. If other dihedrals are centered
1718 : !> on B-C, they rotate as well to keep the relationship between the
1719 : !> dihedrals the same. Atoms A and D are moved amounts related to their
1720 : !> masses (and masses of all connecting atoms), so that heavier segments
1721 : !> are moved less. All atoms except B and C are rotated around the
1722 : !> B-C bond vector (B and C are not moved).
1723 : !> \param r_old the initial coordinates of all molecules in the system
1724 : !> \param r_new the new coordinates of all molecules in the system
1725 : !> \param mc_par the mc parameters for the force env
1726 : !> \param molecule_type the type of molecule we're playing with
1727 : !> \param molecule_kind the structure containing the molecule information
1728 : !> \param particles the particle_list_type for all particles in the force_env..
1729 : !> used to grab the mass of each atom
1730 : !> \param rng_stream the random number stream that we draw from
1731 : !> \author MJM
1732 : ! **************************************************************************************************
1733 0 : SUBROUTINE change_dihedral(r_old, r_new, mc_par, molecule_type, molecule_kind, &
1734 : particles, rng_stream)
1735 :
1736 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_old
1737 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: r_new
1738 : TYPE(mc_simpar_type), POINTER :: mc_par
1739 : INTEGER, INTENT(IN) :: molecule_type
1740 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1741 : TYPE(particle_list_type), POINTER :: particles
1742 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1743 :
1744 : CHARACTER(len=*), PARAMETER :: routineN = 'change_dihedral'
1745 :
1746 : INTEGER :: handle, i, iatom, ibond, ipart, natom, &
1747 : nbond, ntorsion, source, torsion_number
1748 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_d, counter
1749 0 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity
1750 0 : INTEGER, DIMENSION(:), POINTER :: nunits
1751 : LOGICAL :: ionode
1752 0 : REAL(dp), DIMENSION(:), POINTER :: rmdihedral
1753 : REAL(KIND=dp) :: atom_mass, dis_angle, dis_angle_a, &
1754 : dis_angle_d, mass_a, mass_d, &
1755 : old_length_a, rand, u, v, w, x, y, z
1756 : REAL(KIND=dp), DIMENSION(1:3) :: bond_a, temp
1757 0 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
1758 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1759 : TYPE(mp_comm_type) :: group
1760 0 : TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list
1761 :
1762 : ! begin the timing of the subroutine
1763 :
1764 0 : CALL timeset(routineN, handle)
1765 :
1766 0 : NULLIFY (rmdihedral, torsion_list, bond_list, mc_molecule_info)
1767 :
1768 : ! get some stuff from mc_par
1769 : CALL get_mc_par(mc_par, rmdihedral=rmdihedral, &
1770 : source=source, group=group, ionode=ionode, &
1771 0 : mc_molecule_info=mc_molecule_info)
1772 0 : CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
1773 :
1774 : ! copy the incoming coordinates so we can change them
1775 0 : DO ipart = 1, nunits(molecule_type)
1776 0 : r_new(1:3, ipart) = r_old(1:3, ipart)
1777 : END DO
1778 :
1779 : ! pick which bond in the molecule at random
1780 0 : IF (ionode) THEN
1781 0 : rand = rng_stream%next()
1782 : ! CALL RANDOM_NUMBER(rand)
1783 : END IF
1784 0 : CALL group%bcast(rand, source)
1785 : CALL get_molecule_kind(molecule_kind, natom=natom, &
1786 : bond_list=bond_list, nbond=nbond, &
1787 0 : ntorsion=ntorsion, torsion_list=torsion_list)
1788 0 : torsion_number = CEILING(rand*REAL(ntorsion, dp))
1789 :
1790 0 : ALLOCATE (connection(1:natom, 1:2))
1791 : ! assume at most six bonds per atom
1792 0 : ALLOCATE (connectivity(1:6, 1:natom))
1793 0 : ALLOCATE (counter(1:natom))
1794 0 : ALLOCATE (atom_a(1:natom))
1795 0 : ALLOCATE (atom_d(1:natom))
1796 0 : connection(:, :) = 0
1797 0 : connectivity(:, :) = 0
1798 0 : counter(:) = 0
1799 0 : atom_a(:) = 0
1800 0 : atom_d(:) = 0
1801 :
1802 : ! now we need to find a list of atoms that each atom in this bond is connected
1803 : ! to
1804 0 : DO iatom = 1, natom
1805 0 : DO ibond = 1, nbond
1806 0 : IF (bond_list(ibond)%a == iatom) THEN
1807 0 : counter(iatom) = counter(iatom) + 1
1808 0 : connectivity(counter(iatom), iatom) = bond_list(ibond)%b
1809 0 : ELSE IF (bond_list(ibond)%b == iatom) THEN
1810 0 : counter(iatom) = counter(iatom) + 1
1811 0 : connectivity(counter(iatom), iatom) = bond_list(ibond)%a
1812 : END IF
1813 : END DO
1814 : END DO
1815 :
1816 : ! now I need to do a depth first search to figure out which atoms are on atom
1817 : ! a's side and which are on atom d's, but remember we're moving all atoms on a's
1818 : ! side of b, including atoms not in a's branch
1819 0 : atom_a(:) = 0
1820 0 : atom_a(torsion_list(torsion_number)%a) = 1
1821 : CALL depth_first_search(torsion_list(torsion_number)%b, &
1822 0 : torsion_list(torsion_number)%c, connectivity(:, :), atom_a(:))
1823 0 : atom_d(:) = 0
1824 0 : atom_d(torsion_list(torsion_number)%d) = 1
1825 : CALL depth_first_search(torsion_list(torsion_number)%c, &
1826 0 : torsion_list(torsion_number)%b, connectivity(:, :), atom_d(:))
1827 :
1828 : ! now figure out the masses of the various sides, so we can weight how far we
1829 : ! move each group of atoms
1830 0 : mass_a = 0.0_dp
1831 0 : mass_d = 0.0_dp
1832 0 : DO iatom = 1, natom
1833 : CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
1834 0 : mass=atom_mass)
1835 0 : IF (atom_a(iatom) == 1) mass_a = mass_a + atom_mass
1836 0 : IF (atom_d(iatom) == 1) mass_d = mass_d + atom_mass
1837 : END DO
1838 :
1839 : ! choose a displacement
1840 0 : IF (ionode) rand = rng_stream%next()
1841 0 : CALL group%bcast(rand, source)
1842 :
1843 0 : dis_angle = rmdihedral(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
1844 :
1845 : ! find the bond vectors, B-C, so we know what to rotate around
1846 0 : DO i = 1, 3
1847 : bond_a(i) = r_new(i, torsion_list(torsion_number)%c) - &
1848 0 : r_new(i, torsion_list(torsion_number)%b)
1849 : END DO
1850 0 : old_length_a = NORM2(bond_a)
1851 0 : bond_a(1:3) = bond_a(1:3)/old_length_a
1852 :
1853 : ! figure out how much we move each side, since we're mass-weighting, by the
1854 : ! opposite masses, so lighter moves farther...we take the opposite sign of d
1855 : ! so we're not rotating both angles in the same direction
1856 0 : dis_angle_a = dis_angle*mass_d/(mass_a + mass_d)
1857 0 : dis_angle_d = -dis_angle*mass_a/(mass_a + mass_d)
1858 :
1859 0 : DO iatom = 1, natom
1860 :
1861 0 : IF (atom_a(iatom) == 1) THEN
1862 : ! shift the coords so b is at the origin
1863 : r_new(1:3, iatom) = r_new(1:3, iatom) - &
1864 0 : r_new(1:3, torsion_list(torsion_number)%b)
1865 :
1866 : ! multiply by the rotation matrix
1867 0 : u = bond_a(1)
1868 0 : v = bond_a(2)
1869 0 : w = bond_a(3)
1870 0 : x = r_new(1, iatom)
1871 0 : y = r_new(2, iatom)
1872 0 : z = r_new(3, iatom)
1873 : temp(1) = (u*(u*x + v*y + w*z) + (x*(v**2 + w**2) - u*(v*y + w*z))*COS(dis_angle_a) + &
1874 0 : SQRT(u**2 + v**2 + w**2)*(v*z - w*y)*SIN(dis_angle_a))/(u**2 + v**2 + w**2)
1875 : temp(2) = (v*(u*x + v*y + w*z) + (y*(u**2 + w**2) - v*(u*x + w*z))*COS(dis_angle_a) + &
1876 0 : SQRT(u**2 + v**2 + w**2)*(w*x - u*z)*SIN(dis_angle_a))/(u**2 + v**2 + w**2)
1877 : temp(3) = (w*(u*x + v*y + w*z) + (z*(v**2 + u**2) - w*(u*x + v*y))*COS(dis_angle_a) + &
1878 0 : SQRT(u**2 + v**2 + w**2)*(u*y - v*x)*SIN(dis_angle_a))/(u**2 + v**2 + w**2)
1879 :
1880 : ! shift back to the original position
1881 0 : temp(1:3) = temp(1:3) + r_new(1:3, torsion_list(torsion_number)%b)
1882 0 : r_new(1:3, iatom) = temp(1:3)
1883 :
1884 0 : ELSE IF (atom_d(iatom) == 1) THEN
1885 :
1886 : ! shift the coords so c is at the origin
1887 : r_new(1:3, iatom) = r_new(1:3, iatom) - &
1888 0 : r_new(1:3, torsion_list(torsion_number)%c)
1889 :
1890 : ! multiply by the rotation matrix
1891 0 : u = bond_a(1)
1892 0 : v = bond_a(2)
1893 0 : w = bond_a(3)
1894 0 : x = r_new(1, iatom)
1895 0 : y = r_new(2, iatom)
1896 0 : z = r_new(3, iatom)
1897 : temp(1) = (u*(u*x + v*y + w*z) + (x*(v**2 + w**2) - u*(v*y + w*z))*COS(dis_angle_d) + &
1898 0 : SQRT(u**2 + v**2 + w**2)*(v*z - w*y)*SIN(dis_angle_d))/(u**2 + v**2 + w**2)
1899 : temp(2) = (v*(u*x + v*y + w*z) + (y*(u**2 + w**2) - v*(u*x + w*z))*COS(dis_angle_d) + &
1900 0 : SQRT(u**2 + v**2 + w**2)*(w*x - u*z)*SIN(dis_angle_d))/(u**2 + v**2 + w**2)
1901 : temp(3) = (w*(u*x + v*y + w*z) + (z*(v**2 + u**2) - w*(u*x + v*y))*COS(dis_angle_d) + &
1902 0 : SQRT(u**2 + v**2 + w**2)*(u*y - v*x)*SIN(dis_angle_d))/(u**2 + v**2 + w**2)
1903 :
1904 : ! shift back to the original position
1905 0 : temp(1:3) = temp(1:3) + r_new(1:3, torsion_list(torsion_number)%c)
1906 0 : r_new(1:3, iatom) = temp(1:3)
1907 : END IF
1908 : END DO
1909 :
1910 : ! deallocate some stuff
1911 0 : DEALLOCATE (connection)
1912 0 : DEALLOCATE (connectivity)
1913 0 : DEALLOCATE (counter)
1914 0 : DEALLOCATE (atom_a)
1915 0 : DEALLOCATE (atom_d)
1916 :
1917 : ! end the timing
1918 0 : CALL timestop(handle)
1919 :
1920 0 : END SUBROUTINE change_dihedral
1921 :
1922 : ! **************************************************************************************************
1923 : !> \brief performs either a bond or angle change move for a given molecule
1924 : !> \param mc_par the mc parameters for the force env
1925 : !> \param force_env the force environment used in the move
1926 : !> \param bias_env the force environment used to bias the move, if any (it may
1927 : !> be null if lbias=.false. in mc_par)
1928 : !> \param moves the structure that keeps track of how many moves have been
1929 : !> accepted/rejected
1930 : !> \param energy_check the running energy difference between now and the initial
1931 : !> energy
1932 : !> \param r_old the coordinates of force_env before the move
1933 : !> \param old_energy the energy of the force_env before the move
1934 : !> \param start_atom_swap the number of the swap molecule's first atom, assuming the rest of
1935 : !> the atoms follow sequentially
1936 : !> \param target_atom the number of the target atom for swapping
1937 : !> \param molecule_type the molecule type for the atom we're swapping
1938 : !> \param box_number the number of the box we're doing this move in
1939 : !> \param bias_energy_old the biased energy of the system before the move
1940 : !> \param last_bias_energy the last biased energy of the system
1941 : !> \param move_type dictates if we're moving to an "in" or "out" region
1942 : !> \param rng_stream the random number stream that we draw from
1943 : !> \author MJM
1944 : !> \note Designed for parallel.
1945 : ! **************************************************************************************************
1946 0 : SUBROUTINE mc_avbmc_move(mc_par, force_env, bias_env, moves, &
1947 0 : energy_check, r_old, old_energy, start_atom_swap, &
1948 : target_atom, &
1949 : molecule_type, box_number, bias_energy_old, last_bias_energy, &
1950 : move_type, rng_stream)
1951 :
1952 : TYPE(mc_simpar_type), POINTER :: mc_par
1953 : TYPE(force_env_type), POINTER :: force_env, bias_env
1954 : TYPE(mc_moves_type), POINTER :: moves
1955 : REAL(KIND=dp), INTENT(INOUT) :: energy_check
1956 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r_old
1957 : REAL(KIND=dp), INTENT(INOUT) :: old_energy
1958 : INTEGER, INTENT(IN) :: start_atom_swap, target_atom, &
1959 : molecule_type, box_number
1960 : REAL(KIND=dp), INTENT(INOUT) :: bias_energy_old, last_bias_energy
1961 : CHARACTER(LEN=*), INTENT(IN) :: move_type
1962 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1963 :
1964 : CHARACTER(len=*), PARAMETER :: routineN = 'mc_avbmc_move'
1965 :
1966 : INTEGER :: end_mol, handle, ipart, jbox, natom, &
1967 : nswapmoves, source, start_mol
1968 0 : INTEGER, DIMENSION(:), POINTER :: avbmc_atom, mol_type, nunits, nunits_tot
1969 0 : INTEGER, DIMENSION(:, :), POINTER :: nchains
1970 : LOGICAL :: ionode, lbias, ldum, lin, loverlap
1971 0 : REAL(dp), DIMENSION(:), POINTER :: avbmc_rmax, avbmc_rmin, pbias
1972 0 : REAL(dp), DIMENSION(:, :), POINTER :: mass
1973 : REAL(KIND=dp) :: BETA, bias_energy_new, del_quickstep_energy, distance, exp_max_val, &
1974 : exp_min_val, max_val, min_val, new_energy, prefactor, rand, rdum, volume_in, volume_out, &
1975 : w, weight_new, weight_old
1976 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r_new
1977 : REAL(KIND=dp), DIMENSION(1:3) :: abc, RIJ
1978 : TYPE(cell_type), POINTER :: cell
1979 : TYPE(cp_subsys_type), POINTER :: subsys, subsys_force
1980 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1981 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
1982 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1983 : TYPE(mp_comm_type) :: group
1984 : TYPE(particle_list_type), POINTER :: particles, particles_force
1985 :
1986 0 : rdum = 1.0_dp
1987 :
1988 : ! begin the timing of the subroutine
1989 0 : CALL timeset(routineN, handle)
1990 :
1991 : ! get a bunch of stuff from mc_par
1992 : CALL get_mc_par(mc_par, lbias=lbias, &
1993 : BETA=BETA, max_val=max_val, min_val=min_val, exp_max_val=exp_max_val, &
1994 : exp_min_val=exp_min_val, avbmc_atom=avbmc_atom, &
1995 : avbmc_rmin=avbmc_rmin, avbmc_rmax=avbmc_rmax, &
1996 : nswapmoves=nswapmoves, ionode=ionode, source=source, &
1997 0 : group=group, pbias=pbias, mc_molecule_info=mc_molecule_info)
1998 : CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, &
1999 0 : mass=mass, nunits=nunits, nunits_tot=nunits_tot, mol_type=mol_type)
2000 : ! figure out some bounds for mol_type
2001 0 : start_mol = 1
2002 0 : DO jbox = 1, box_number - 1
2003 0 : start_mol = start_mol + SUM(nchains(:, jbox))
2004 : END DO
2005 0 : end_mol = start_mol + SUM(nchains(:, box_number)) - 1
2006 :
2007 : ! nullify some pointers
2008 0 : NULLIFY (particles, subsys, molecule_kinds, molecule_kind, &
2009 0 : particles_force, subsys_force)
2010 :
2011 : ! do some allocation
2012 0 : ALLOCATE (r_new(1:3, 1:nunits_tot(box_number)))
2013 :
2014 : ! now we need to grab and save coordinates, in case we reject
2015 : ! are we biasing this move?
2016 0 : IF (lbias) THEN
2017 :
2018 : ! grab the coordinates
2019 0 : CALL force_env_get(bias_env, cell=cell, subsys=subsys)
2020 : CALL cp_subsys_get(subsys, &
2021 0 : particles=particles, molecule_kinds=molecule_kinds)
2022 0 : molecule_kind => molecule_kinds%els(1)
2023 0 : CALL get_molecule_kind(molecule_kind, natom=natom)
2024 0 : CALL get_cell(cell, abc=abc)
2025 :
2026 : ! save the energy
2027 : ! bias_energy_old=bias_energy
2028 :
2029 : ELSE
2030 :
2031 : ! grab the coordinates
2032 0 : CALL force_env_get(force_env, cell=cell, subsys=subsys)
2033 : CALL cp_subsys_get(subsys, &
2034 0 : particles=particles, molecule_kinds=molecule_kinds)
2035 0 : molecule_kind => molecule_kinds%els(1)
2036 0 : CALL get_molecule_kind(molecule_kind, natom=natom)
2037 0 : CALL get_cell(cell, abc=abc)
2038 :
2039 : END IF
2040 :
2041 : ! let's determine if the molecule to be moved is in the "in" region or the
2042 : ! "out" region of the target
2043 : RIJ(1) = particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(1) - &
2044 : particles%els(target_atom)%r(1) - abc(1)*ANINT( &
2045 : (particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(1) - &
2046 0 : particles%els(target_atom)%r(1))/abc(1))
2047 : RIJ(2) = particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(2) - &
2048 : particles%els(target_atom)%r(2) - abc(2)*ANINT( &
2049 : (particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(2) - &
2050 0 : particles%els(target_atom)%r(2))/abc(2))
2051 : RIJ(3) = particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(3) - &
2052 : particles%els(target_atom)%r(3) - abc(3)*ANINT( &
2053 : (particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(3) - &
2054 0 : particles%els(target_atom)%r(3))/abc(3))
2055 0 : distance = SQRT(RIJ(1)**2 + RIJ(2)**2 + RIJ(3)**2)
2056 0 : IF (distance <= avbmc_rmax(molecule_type) .AND. distance >= avbmc_rmin(molecule_type)) THEN
2057 : lin = .TRUE.
2058 : ELSE
2059 : lin = .FALSE.
2060 : END IF
2061 :
2062 : ! increment the counter of the particular move we've done
2063 : ! swapping into the "in" region of mol_target
2064 : IF (lin) THEN
2065 0 : IF (move_type == 'in') THEN
2066 : moves%avbmc_inin%attempts = &
2067 0 : moves%avbmc_inin%attempts + 1
2068 : ELSE
2069 : moves%avbmc_inout%attempts = &
2070 0 : moves%avbmc_inout%attempts + 1
2071 : END IF
2072 : ELSE
2073 0 : IF (move_type == 'in') THEN
2074 : moves%avbmc_outin%attempts = &
2075 0 : moves%avbmc_outin%attempts + 1
2076 : ELSE
2077 : moves%avbmc_outout%attempts = &
2078 0 : moves%avbmc_outout%attempts + 1
2079 : END IF
2080 : END IF
2081 :
2082 0 : IF (lbias) THEN
2083 :
2084 0 : IF (move_type == 'in') THEN
2085 :
2086 : ! do CBMC for the old config
2087 : CALL generate_cbmc_swap_config(bias_env, BETA, max_val, min_val, exp_max_val, &
2088 : exp_min_val, nswapmoves, &
2089 : weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2090 : mass(:, molecule_type), ldum, rdum, &
2091 : bias_energy_old, ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2092 : source, group, rng_stream, &
2093 : avbmc_atom=avbmc_atom(molecule_type), &
2094 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='out', &
2095 0 : target_atom=target_atom)
2096 :
2097 : ELSE
2098 :
2099 : ! do CBMC for the old config
2100 : CALL generate_cbmc_swap_config(bias_env, BETA, max_val, min_val, exp_max_val, &
2101 : exp_min_val, nswapmoves, &
2102 : weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2103 : mass(:, molecule_type), ldum, rdum, &
2104 : bias_energy_old, ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2105 : source, group, rng_stream, &
2106 : avbmc_atom=avbmc_atom(molecule_type), &
2107 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='in', &
2108 0 : target_atom=target_atom)
2109 :
2110 : END IF
2111 :
2112 : ! generate the new config
2113 : CALL generate_cbmc_swap_config(bias_env, BETA, max_val, min_val, exp_max_val, &
2114 : exp_min_val, nswapmoves, &
2115 : weight_new, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2116 : mass(:, molecule_type), loverlap, bias_energy_new, &
2117 : bias_energy_old, ionode, .FALSE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2118 : source, group, rng_stream, &
2119 : avbmc_atom=avbmc_atom(molecule_type), &
2120 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type=move_type, &
2121 0 : target_atom=target_atom)
2122 :
2123 : ! the energy that comes out of the above routine is the difference...we want
2124 : ! the real energy for the acceptance rule...we don't do this for the
2125 : ! lbias=.false. case because it doesn't appear in the acceptance rule, and
2126 : ! we compensate in case of acceptance
2127 0 : bias_energy_new = bias_energy_new + bias_energy_old
2128 :
2129 : ELSE
2130 :
2131 0 : IF (move_type == 'in') THEN
2132 :
2133 : ! find the weight of the old config
2134 : CALL generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, &
2135 : exp_min_val, nswapmoves, &
2136 : weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2137 : mass(:, molecule_type), ldum, rdum, old_energy, &
2138 : ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2139 : source, group, rng_stream, &
2140 : avbmc_atom=avbmc_atom(molecule_type), &
2141 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='out', &
2142 0 : target_atom=target_atom)
2143 :
2144 : ELSE
2145 :
2146 : ! find the weight of the old config
2147 : CALL generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, &
2148 : exp_min_val, nswapmoves, &
2149 : weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2150 : mass(:, molecule_type), ldum, rdum, old_energy, &
2151 : ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2152 : source, group, rng_stream, &
2153 : avbmc_atom=avbmc_atom(molecule_type), &
2154 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='in', &
2155 0 : target_atom=target_atom)
2156 :
2157 : END IF
2158 :
2159 : ! generate the new config...do this after, because it changes the force_env
2160 : CALL generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, &
2161 : exp_min_val, nswapmoves, &
2162 : weight_new, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2163 : mass(:, molecule_type), loverlap, new_energy, old_energy, &
2164 : ionode, .FALSE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2165 : source, group, rng_stream, &
2166 : avbmc_atom=avbmc_atom(molecule_type), &
2167 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type=move_type, &
2168 0 : target_atom=target_atom)
2169 :
2170 : END IF
2171 :
2172 0 : IF (loverlap) THEN
2173 0 : DEALLOCATE (r_new)
2174 :
2175 : ! need to reset the old coordinates
2176 0 : IF (lbias) THEN
2177 0 : CALL force_env_get(bias_env, subsys=subsys)
2178 0 : CALL cp_subsys_get(subsys, particles=particles)
2179 : ELSE
2180 0 : CALL force_env_get(force_env, subsys=subsys)
2181 0 : CALL cp_subsys_get(subsys, particles=particles)
2182 : END IF
2183 0 : DO ipart = 1, nunits_tot(box_number)
2184 0 : particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2185 : END DO
2186 :
2187 0 : CALL timestop(handle)
2188 :
2189 : RETURN
2190 : END IF
2191 :
2192 : ! if we're biasing, we need to compute the new energy with the full
2193 : ! potential
2194 0 : IF (lbias) THEN
2195 : ! need to give the force_env the coords from the bias_env
2196 0 : CALL force_env_get(force_env, subsys=subsys_force)
2197 0 : CALL cp_subsys_get(subsys_force, particles=particles_force)
2198 0 : CALL force_env_get(bias_env, subsys=subsys)
2199 0 : CALL cp_subsys_get(subsys, particles=particles)
2200 0 : DO ipart = 1, nunits_tot(box_number)
2201 0 : particles_force%els(ipart)%r(1:3) = particles%els(ipart)%r(1:3)
2202 : END DO
2203 :
2204 : CALL force_env_calc_energy_force(force_env, &
2205 0 : calc_force=.FALSE.)
2206 : CALL force_env_get(force_env, &
2207 0 : potential_energy=new_energy)
2208 :
2209 : END IF
2210 :
2211 0 : volume_in = 4.0_dp/3.0_dp*pi*(avbmc_rmax(molecule_type)**3 - avbmc_rmin(molecule_type)**3)
2212 0 : volume_out = abc(1)*abc(2)*abc(3) - volume_in
2213 :
2214 0 : IF (lin .AND. move_type == 'in' .OR. &
2215 : .NOT. lin .AND. move_type == 'out') THEN
2216 : ! standard Metropolis rule
2217 : prefactor = 1.0_dp
2218 0 : ELSE IF (.NOT. lin .AND. move_type == 'in') THEN
2219 0 : prefactor = (1.0_dp - pbias(molecule_type))*volume_in/(pbias(molecule_type)*volume_out)
2220 : ELSE
2221 0 : prefactor = pbias(molecule_type)*volume_out/((1.0_dp - pbias(molecule_type))*volume_in)
2222 : END IF
2223 :
2224 0 : IF (lbias) THEN
2225 : ! AVBMC with CBMC and a biasing potential...notice that if the biasing
2226 : ! potential equals the quickstep potential, this cancels out to the
2227 : ! acceptance below
2228 : del_quickstep_energy = (-BETA)*(new_energy - old_energy - &
2229 0 : (bias_energy_new - bias_energy_old))
2230 :
2231 0 : IF (del_quickstep_energy > exp_max_val) THEN
2232 0 : del_quickstep_energy = max_val
2233 0 : ELSE IF (del_quickstep_energy < exp_min_val) THEN
2234 : del_quickstep_energy = 0.0_dp
2235 : ELSE
2236 0 : del_quickstep_energy = EXP(del_quickstep_energy)
2237 : END IF
2238 :
2239 0 : w = prefactor*del_quickstep_energy*weight_new/weight_old
2240 :
2241 : ELSE
2242 :
2243 : ! AVBMC with CBMC
2244 0 : w = prefactor*weight_new/weight_old
2245 : END IF
2246 :
2247 : ! check if the move is accepted
2248 0 : IF (w >= 1.0E0_dp) THEN
2249 0 : rand = 0.0E0_dp
2250 : ELSE
2251 0 : IF (ionode) rand = rng_stream%next()
2252 0 : CALL group%bcast(rand, source)
2253 : END IF
2254 :
2255 0 : IF (rand < w) THEN
2256 :
2257 : ! accept the move
2258 :
2259 0 : IF (lin) THEN
2260 0 : IF (move_type == 'in') THEN
2261 : moves%avbmc_inin%successes = &
2262 0 : moves%avbmc_inin%successes + 1
2263 : ELSE
2264 : moves%avbmc_inout%successes = &
2265 0 : moves%avbmc_inout%successes + 1
2266 : END IF
2267 : ELSE
2268 0 : IF (move_type == 'in') THEN
2269 : moves%avbmc_outin%successes = &
2270 0 : moves%avbmc_outin%successes + 1
2271 : ELSE
2272 : moves%avbmc_outout%successes = &
2273 0 : moves%avbmc_outout%successes + 1
2274 : END IF
2275 : END IF
2276 :
2277 : ! we need to compensate for the fact that we take the difference in
2278 : ! generate_cbmc_config to keep the exponetials small
2279 0 : IF (.NOT. lbias) THEN
2280 0 : new_energy = new_energy + old_energy
2281 : END IF
2282 :
2283 : ! update energies
2284 0 : energy_check = energy_check + (new_energy - old_energy)
2285 0 : old_energy = new_energy
2286 :
2287 : ! if we're biasing the update the biasing energy
2288 0 : IF (lbias) THEN
2289 : ! need to do this outside of the routine
2290 0 : last_bias_energy = bias_energy_new
2291 0 : bias_energy_old = bias_energy_new
2292 : END IF
2293 :
2294 : ! update coordinates
2295 0 : CALL force_env_get(force_env, subsys=subsys)
2296 0 : CALL cp_subsys_get(subsys, particles=particles)
2297 0 : DO ipart = 1, nunits_tot(box_number)
2298 0 : r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
2299 : END DO
2300 : ELSE
2301 : ! reject the move...need to restore the old coordinates
2302 0 : IF (lbias) THEN
2303 0 : CALL force_env_get(bias_env, subsys=subsys)
2304 0 : CALL cp_subsys_get(subsys, particles=particles)
2305 0 : DO ipart = 1, nunits_tot(box_number)
2306 0 : particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2307 : END DO
2308 0 : CALL cp_subsys_set(subsys, particles=particles)
2309 : END IF
2310 0 : CALL force_env_get(force_env, subsys=subsys)
2311 0 : CALL cp_subsys_get(subsys, particles=particles)
2312 0 : DO ipart = 1, nunits_tot(box_number)
2313 0 : particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2314 : END DO
2315 0 : CALL cp_subsys_set(subsys, particles=particles)
2316 :
2317 : END IF
2318 :
2319 : ! deallocate some stuff
2320 0 : DEALLOCATE (r_new)
2321 : ! end the timing
2322 0 : CALL timestop(handle)
2323 :
2324 0 : END SUBROUTINE mc_avbmc_move
2325 :
2326 : ! **************************************************************************************************
2327 : !> \brief performs a hybrid Monte Carlo move that runs a short MD sequence
2328 : !> \param mc_par the mc parameters for the force env
2329 : !> \param force_env the force environment whose cell we're changing
2330 : !> \param globenv ...
2331 : !> \param moves the structure that keeps track of how many moves have been
2332 : !> accepted/rejected
2333 : !> \param move_updates the structure that keeps track of how many moves have
2334 : !> been accepted/rejected since the last time the displacements
2335 : !> were updated
2336 : !> \param old_energy the energy of the last accepted move involving an
2337 : !> unbiased calculation
2338 : !> \param box_number the box we're changing the volume of
2339 : !> \param energy_check the running total of how much the energy has changed
2340 : !> since the initial configuration
2341 : !> \param r_old the coordinates of the last accepted move involving an
2342 : !> unbiased calculation
2343 : !> \param rng_stream the random number stream that we draw from
2344 : !> \author MJM
2345 : !> \note Designed for parallel use.
2346 : ! **************************************************************************************************
2347 20 : SUBROUTINE mc_hmc_move(mc_par, force_env, globenv, moves, move_updates, &
2348 : old_energy, box_number, &
2349 20 : energy_check, r_old, rng_stream)
2350 :
2351 : TYPE(mc_simpar_type), POINTER :: mc_par
2352 : TYPE(force_env_type), POINTER :: force_env
2353 : TYPE(global_environment_type), POINTER :: globenv
2354 : TYPE(mc_moves_type), POINTER :: moves, move_updates
2355 : REAL(KIND=dp), INTENT(INOUT) :: old_energy
2356 : INTEGER, INTENT(IN) :: box_number
2357 : REAL(KIND=dp), INTENT(INOUT) :: energy_check
2358 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r_old
2359 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
2360 :
2361 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_hmc_move'
2362 :
2363 : INTEGER :: handle, iatom, source
2364 20 : INTEGER, DIMENSION(:), POINTER :: nunits_tot
2365 : LOGICAL :: ionode
2366 : REAL(KIND=dp) :: BETA, energy_term, exp_max_val, &
2367 : exp_min_val, new_energy, rand, value, w
2368 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r
2369 : TYPE(cp_subsys_type), POINTER :: oldsys
2370 : TYPE(mc_ekin_type), POINTER :: hmc_ekin
2371 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
2372 : TYPE(mp_comm_type) :: group
2373 : TYPE(particle_list_type), POINTER :: particles_old
2374 :
2375 : ! begin the timing of the subroutine
2376 :
2377 20 : CALL timeset(routineN, handle)
2378 :
2379 : ! get a bunch of stuff from mc_par
2380 : CALL get_mc_par(mc_par, ionode=ionode, &
2381 : BETA=BETA, exp_max_val=exp_max_val, &
2382 : exp_min_val=exp_min_val, source=source, group=group, &
2383 20 : mc_molecule_info=mc_molecule_info)
2384 20 : CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot)
2385 :
2386 : ! nullify some pointers
2387 20 : NULLIFY (particles_old, oldsys, hmc_ekin)
2388 :
2389 : ! do some allocation
2390 60 : ALLOCATE (r(1:3, 1:nunits_tot(box_number)))
2391 20 : ALLOCATE (hmc_ekin)
2392 :
2393 : ! record the attempt
2394 20 : moves%hmc%attempts = moves%hmc%attempts + 1
2395 20 : move_updates%hmc%attempts = move_updates%hmc%attempts + 1
2396 :
2397 : ! now let's grab the particle positions
2398 20 : CALL force_env_get(force_env, subsys=oldsys)
2399 20 : CALL cp_subsys_get(oldsys, particles=particles_old)
2400 :
2401 : ! save the old coordinates
2402 21700 : DO iatom = 1, nunits_tot(box_number)
2403 86740 : r(1:3, iatom) = particles_old%els(iatom)%r(1:3)
2404 : END DO
2405 :
2406 : ! now run the MD simulation
2407 20 : CALL qs_mol_dyn(force_env, globenv, hmc_e_initial=hmc_ekin%initial_ekin, hmc_e_final=hmc_ekin%final_ekin)
2408 :
2409 : ! get the energy
2410 : CALL force_env_get(force_env, &
2411 20 : potential_energy=new_energy)
2412 :
2413 : ! accept or reject the move
2414 : ! to prevent overflows
2415 20 : energy_term = new_energy + hmc_ekin%final_ekin - old_energy - hmc_ekin%initial_ekin
2416 :
2417 20 : value = -BETA*(energy_term)
2418 20 : IF (value > exp_max_val) THEN
2419 : w = 10.0_dp
2420 20 : ELSE IF (value < exp_min_val) THEN
2421 : w = 0.0_dp
2422 : ELSE
2423 20 : w = EXP(value)
2424 : END IF
2425 :
2426 20 : IF (w >= 1.0E0_dp) THEN
2427 10 : w = 1.0E0_dp
2428 10 : rand = 0.0E0_dp
2429 : ELSE
2430 10 : IF (ionode) rand = rng_stream%next()
2431 10 : CALL group%bcast(rand, source)
2432 : END IF
2433 :
2434 20 : IF (rand < w) THEN
2435 :
2436 : ! accept the move
2437 14 : moves%hmc%successes = moves%hmc%successes + 1
2438 14 : move_updates%hmc%successes = move_updates%hmc%successes + 1
2439 :
2440 : ! update energies
2441 14 : energy_check = energy_check + (new_energy - old_energy)
2442 14 : old_energy = new_energy
2443 :
2444 15190 : DO iatom = 1, nunits_tot(box_number)
2445 60718 : r_old(1:3, iatom) = particles_old%els(iatom)%r(1:3)
2446 : END DO
2447 :
2448 : ELSE
2449 :
2450 : ! reset the cell and particle positions
2451 6510 : DO iatom = 1, nunits_tot(box_number)
2452 26022 : particles_old%els(iatom)%r(1:3) = r_old(1:3, iatom)
2453 : END DO
2454 :
2455 : END IF
2456 :
2457 : ! deallocate some stuff
2458 20 : DEALLOCATE (r)
2459 20 : DEALLOCATE (hmc_ekin)
2460 :
2461 : ! end the timing
2462 20 : CALL timestop(handle)
2463 :
2464 40 : END SUBROUTINE mc_hmc_move
2465 :
2466 : ! *****************************************************************************
2467 : !> \brief translates the cluster randomly in either the x,y, or z
2468 : !>direction
2469 : !> \param mc_par the mc parameters for the force env
2470 : !> \param force_env the force environment used in the move
2471 : !> \param bias_env the force environment used to bias the move, if any (it may
2472 : !> be null if lbias=.false. in mc_par)
2473 : !> \param moves the structure that keeps track of how many moves have been
2474 : !> accepted/rejected
2475 : !> \param move_updates the structure that keeps track of how many moves have
2476 : !> been accepted/rejected since the last time the displacements
2477 : !> were updated
2478 : !> \param box_number ...
2479 : !> \param bias_energy the biased energy of the system before the move
2480 : !> \param lreject set to .true. if there is an overlap
2481 : !> \param rng_stream the random number stream that we draw from
2482 : !> \author Himanshu Goel
2483 : !> \note Designed for parallel use.
2484 : ! **************************************************************************************************
2485 :
2486 10 : SUBROUTINE mc_cluster_translation(mc_par, force_env, bias_env, moves, &
2487 : move_updates, box_number, bias_energy, lreject, rng_stream)
2488 :
2489 : TYPE(mc_simpar_type), POINTER :: mc_par
2490 : TYPE(force_env_type), POINTER :: force_env, bias_env
2491 : TYPE(mc_moves_type), POINTER :: moves, move_updates
2492 : INTEGER, INTENT(IN) :: box_number
2493 : REAL(KIND=dp), INTENT(INOUT) :: bias_energy
2494 : LOGICAL, INTENT(OUT) :: lreject
2495 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
2496 :
2497 : CHARACTER(len=*), PARAMETER :: routineN = 'mc_cluster_translation'
2498 :
2499 : INTEGER :: cstart, end_mol, handle, imol, ipart, iparticle, iunit, jbox, jpart, junit, &
2500 : move_direction, nend, nunit, source, start_mol, total_clus, total_clusafmo
2501 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: cluster
2502 10 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
2503 10 : INTEGER, DIMENSION(:, :), POINTER :: nchains
2504 : LOGICAL :: ionode, lbias, loverlap
2505 : REAL(KIND=dp) :: BETA, bias_energy_new, bias_energy_old, &
2506 : dis_mol, exp_max_val, exp_min_val, &
2507 : rand, rmcltrans, value, w
2508 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r_old
2509 : TYPE(cp_subsys_type), POINTER :: subsys
2510 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
2511 : TYPE(mp_comm_type) :: group
2512 : TYPE(particle_list_type), POINTER :: particles
2513 :
2514 : ! *** Local Counters ***
2515 : ! begin the timing of the subroutine
2516 :
2517 10 : CALL timeset(routineN, handle)
2518 :
2519 : ! nullify some pointers
2520 10 : NULLIFY (particles, subsys)
2521 :
2522 : ! get a bunch of stuff from mc_par
2523 : CALL get_mc_par(mc_par, lbias=lbias, &
2524 : BETA=BETA, exp_max_val=exp_max_val, &
2525 : exp_min_val=exp_min_val, rmcltrans=rmcltrans, ionode=ionode, source=source, &
2526 10 : group=group, mc_molecule_info=mc_molecule_info)
2527 : CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot, &
2528 10 : nchains=nchains, nunits=nunits, mol_type=mol_type)
2529 :
2530 : ! find out some bounds for mol_type
2531 10 : start_mol = 1
2532 10 : DO jbox = 1, box_number - 1
2533 10 : start_mol = start_mol + SUM(nchains(:, jbox))
2534 : END DO
2535 20 : end_mol = start_mol + SUM(nchains(:, box_number)) - 1
2536 :
2537 : ! do some allocation
2538 30 : ALLOCATE (r_old(1:3, 1:nunits_tot(box_number)))
2539 :
2540 : ! Allocating cluster matrix size
2541 20 : nend = SUM(nchains(:, box_number))
2542 40 : ALLOCATE (cluster(nend, nend))
2543 60 : DO ipart = 1, nend
2544 310 : DO jpart = 1, nend
2545 300 : cluster(ipart, jpart) = 0
2546 : END DO
2547 : END DO
2548 :
2549 : ! Get cluster information in cluster matrix from cluster_search subroutine
2550 10 : IF (lbias) THEN
2551 : CALL cluster_search(mc_par, bias_env, cluster, nchains(:, box_number), &
2552 0 : nunits, mol_type(start_mol:end_mol), total_clus)
2553 : ELSE
2554 : CALL cluster_search(mc_par, force_env, cluster, nchains(:, box_number), &
2555 10 : nunits, mol_type(start_mol:end_mol), total_clus)
2556 : END IF
2557 :
2558 10 : IF (lbias) THEN
2559 :
2560 : ! grab the coordinates
2561 0 : CALL force_env_get(bias_env, subsys=subsys)
2562 0 : CALL cp_subsys_get(subsys, particles=particles)
2563 :
2564 : ! save the coordinates
2565 0 : DO ipart = 1, nunits_tot(box_number)
2566 0 : r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
2567 : END DO
2568 :
2569 : ! save the energy
2570 0 : bias_energy_old = bias_energy
2571 : ELSE
2572 :
2573 : ! grab the coordinates
2574 10 : CALL force_env_get(force_env, subsys=subsys)
2575 10 : CALL cp_subsys_get(subsys, particles=particles)
2576 : END IF
2577 :
2578 : ! record the attempt
2579 10 : moves%cltrans%attempts = moves%cltrans%attempts + 1
2580 10 : move_updates%cltrans%attempts = move_updates%cltrans%attempts + 1
2581 10 : moves%bias_cltrans%attempts = moves%bias_cltrans%attempts + 1
2582 10 : move_updates%bias_cltrans%attempts = move_updates%bias_cltrans%attempts + 1
2583 10 : IF (.NOT. lbias) THEN
2584 10 : moves%cltrans%qsuccesses = moves%cltrans%qsuccesses + 1
2585 10 : move_updates%cltrans%qsuccesses = move_updates%cltrans%qsuccesses + 1
2586 10 : moves%bias_cltrans%qsuccesses = moves%bias_cltrans%qsuccesses + 1
2587 10 : move_updates%bias_cltrans%qsuccesses = move_updates%bias_cltrans%qsuccesses + 1
2588 : END IF
2589 :
2590 : ! call a random number to figure out which direction we're moving
2591 10 : IF (ionode) rand = rng_stream%next()
2592 10 : CALL group%bcast(rand, source)
2593 10 : move_direction = INT(3*rand) + 1
2594 :
2595 : ! call a random number to figure out how far we're moving
2596 10 : IF (ionode) rand = rng_stream%next()
2597 10 : CALL group%bcast(rand, source)
2598 10 : dis_mol = rmcltrans*(rand - 0.5E0_dp)*2.0E0_dp
2599 :
2600 : ! choosing cluster
2601 10 : IF (ionode) rand = rng_stream%next()
2602 10 : CALL group%bcast(rand, source)
2603 10 : jpart = INT(1 + rand*total_clus)
2604 :
2605 : ! do the cluster move
2606 60 : DO cstart = 1, nend
2607 50 : imol = 0
2608 60 : IF (cluster(jpart, cstart) /= 0) THEN
2609 34 : imol = cluster(jpart, cstart)
2610 : iunit = 1
2611 34 : DO ipart = 1, imol - 1
2612 24 : nunit = nunits(mol_type(ipart + start_mol - 1))
2613 34 : iunit = iunit + nunit
2614 : END DO
2615 10 : nunit = nunits(mol_type(imol + start_mol - 1))
2616 10 : junit = iunit + nunit - 1
2617 30 : DO iparticle = iunit, junit
2618 : particles%els(iparticle)%r(move_direction) = &
2619 30 : particles%els(iparticle)%r(move_direction) + dis_mol
2620 : END DO
2621 : END IF
2622 : END DO
2623 10 : CALL cp_subsys_set(subsys, particles=particles)
2624 :
2625 : !Make cluster matrix null
2626 60 : DO ipart = 1, nend
2627 310 : DO jpart = 1, nend
2628 300 : cluster(ipart, jpart) = 0
2629 : END DO
2630 : END DO
2631 :
2632 : ! checking the number of cluster are same or got changed after cluster translation move
2633 10 : IF (lbias) THEN
2634 : CALL cluster_search(mc_par, bias_env, cluster, nchains(:, box_number), &
2635 0 : nunits, mol_type(start_mol:end_mol), total_clusafmo)
2636 : ELSE
2637 : CALL cluster_search(mc_par, force_env, cluster, nchains(:, box_number), &
2638 10 : nunits, mol_type(start_mol:end_mol), total_clusafmo)
2639 : END IF
2640 :
2641 : ! figure out if there is any overlap...need the number of the molecule
2642 10 : lreject = .FALSE.
2643 10 : IF (lbias) THEN
2644 : CALL check_for_overlap(bias_env, nchains(:, box_number), &
2645 0 : nunits(:), loverlap, mol_type(start_mol:end_mol))
2646 : ELSE
2647 : CALL check_for_overlap(force_env, nchains(:, box_number), &
2648 10 : nunits(:), loverlap, mol_type(start_mol:end_mol))
2649 10 : IF (loverlap) lreject = .TRUE.
2650 : END IF
2651 :
2652 : ! check if cluster size changes then reject the move
2653 10 : IF (lbias) THEN
2654 0 : IF (total_clusafmo /= total_clus) THEN
2655 0 : loverlap = .TRUE.
2656 : END IF
2657 : ELSE
2658 10 : IF (total_clusafmo /= total_clus) THEN
2659 0 : loverlap = .TRUE.
2660 0 : lreject = .TRUE.
2661 : END IF
2662 : END IF
2663 :
2664 : ! if we're biasing with a cheaper potential, check for acceptance
2665 10 : IF (lbias) THEN
2666 :
2667 : ! here's where we bias the moves
2668 0 : IF (loverlap) THEN
2669 : w = 0.0E0_dp
2670 : ELSE
2671 0 : CALL force_env_calc_energy_force(bias_env, calc_force=.FALSE.)
2672 : CALL force_env_get(bias_env, &
2673 0 : potential_energy=bias_energy_new)
2674 : ! accept or reject the move based on the Metropolis rule
2675 0 : value = -BETA*(bias_energy_new - bias_energy_old)
2676 0 : IF (value > exp_max_val) THEN
2677 : w = 10.0_dp
2678 0 : ELSE IF (value < exp_min_val) THEN
2679 : w = 0.0_dp
2680 : ELSE
2681 0 : w = EXP(value)
2682 : END IF
2683 :
2684 : END IF
2685 :
2686 0 : IF (w >= 1.0E0_dp) THEN
2687 0 : w = 1.0E0_dp
2688 0 : rand = 0.0E0_dp
2689 : ELSE
2690 0 : IF (ionode) rand = rng_stream%next()
2691 0 : CALL group%bcast(rand, source)
2692 : END IF
2693 0 : IF (rand < w) THEN
2694 :
2695 : ! accept the move
2696 0 : moves%bias_cltrans%successes = moves%bias_cltrans%successes + 1
2697 0 : move_updates%bias_cltrans%successes = move_updates%bias_cltrans%successes + 1
2698 0 : moves%cltrans%qsuccesses = moves%cltrans%qsuccesses + 1
2699 : move_updates%cltrans%successes = &
2700 0 : move_updates%cltrans%successes + 1
2701 0 : moves%qcltrans_dis = moves%qcltrans_dis + ABS(dis_mol)
2702 : bias_energy = bias_energy + bias_energy_new - &
2703 0 : bias_energy_old
2704 :
2705 : ELSE
2706 :
2707 : ! reject the move
2708 : ! restore the coordinates
2709 0 : CALL force_env_get(bias_env, subsys=subsys)
2710 0 : CALL cp_subsys_get(subsys, particles=particles)
2711 0 : DO ipart = 1, nunits_tot(box_number)
2712 0 : particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2713 : END DO
2714 0 : CALL cp_subsys_set(subsys, particles=particles)
2715 :
2716 : END IF
2717 :
2718 : END IF
2719 :
2720 : ! deallocate some stuff
2721 10 : DEALLOCATE (cluster)
2722 10 : DEALLOCATE (r_old)
2723 :
2724 : ! end the timing
2725 10 : CALL timestop(handle)
2726 :
2727 10 : END SUBROUTINE mc_cluster_translation
2728 :
2729 : END MODULE mc_moves
|
