Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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) .NE. 0) THEN
103 1636 : IF (connectivity(iatom, current_atom) .NE. 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 : ELSEIF (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 .GT. exp_max_val) THEN
340 : w = 10.0_dp
341 420 : ELSEIF (value .LT. 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 .GE. 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 .LT. 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 : ELSEIF (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 .GT. exp_max_val) THEN
583 : w = 10.0_dp
584 528 : ELSEIF (value .LT. 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 .GE. 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 .LT. 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 .EQ. 1) THEN
782 : lx = .TRUE.
783 334 : ELSEIF (dir .EQ. 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 : ELSEIF (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 .GT. exp_max_val) THEN
886 : w = 10.0_dp
887 424 : ELSEIF (value .LT. 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 .GE. 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 .LT. 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 .LT. 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) .EQ. 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) .EQ. 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) .LE. 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 .LE. 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 .EQ. use_fist_force) THEN
1198 13 : CALL get_mc_par(mc_par, rcut=rcut)
1199 13 : IF (new_cell_length(1) .LT. 2.0_dp*rcut) ltoo_small = .TRUE.
1200 13 : IF (new_cell_length(2) .LT. 2.0_dp*rcut) ltoo_small = .TRUE.
1201 13 : IF (new_cell_length(3) .LT. 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) &
1212 0 : CPABORT("Attempted a volume move where box size got too small.")
1213 :
1214 : ! now compute the energy
1215 34 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
1216 : CALL force_env_get(force_env, &
1217 34 : potential_energy=new_energy)
1218 :
1219 : ! accept or reject the move
1220 : ! to prevent overflows
1221 34 : energy_term = new_energy - old_energy
1222 : volume_term = -REAL(SUM(nchains(:, box_number)), dp)/BETA* &
1223 : LOG(new_cell_length(1)*new_cell_length(2)*new_cell_length(3)/ &
1224 94 : (old_cell_length(1)*old_cell_length(2)*old_cell_length(3)))
1225 34 : pressure_term = pressure*vol_dis
1226 :
1227 34 : value = -BETA*(energy_term + volume_term + pressure_term)
1228 34 : IF (value .GT. exp_max_val) THEN
1229 : w = 10.0_dp
1230 34 : ELSEIF (value .LT. exp_min_val) THEN
1231 : w = 0.0_dp
1232 : ELSE
1233 34 : w = EXP(value)
1234 : END IF
1235 :
1236 : !!!!!!!!!!!!!!!! for E_V curves
1237 : ! w=1.0E0_dp
1238 : ! w=0.0E0_dp
1239 :
1240 34 : IF (w .GE. 1.0E0_dp) THEN
1241 18 : w = 1.0E0_dp
1242 18 : rand = 0.0E0_dp
1243 : ELSE
1244 16 : IF (ionode) rand = rng_stream%next()
1245 16 : CALL group%bcast(rand, source)
1246 : END IF
1247 :
1248 34 : IF (rand .LT. w) THEN
1249 :
1250 : ! accept the move
1251 30 : moves%volume%successes = moves%volume%successes + 1
1252 30 : move_updates%volume%successes = move_updates%volume%successes + 1
1253 :
1254 : ! update energies
1255 30 : energy_check = energy_check + (new_energy - old_energy)
1256 30 : old_energy = new_energy
1257 :
1258 720 : DO iatom = 1, nunits_tot(box_number)
1259 2790 : r_old(1:3, iatom) = particles_old%els(iatom)%r(1:3)
1260 : END DO
1261 :
1262 : ! update discrete_array if we're doing a discrete volume move
1263 30 : IF (ldiscrete) THEN
1264 : CALL create_discrete_array(new_cell_length(:), &
1265 0 : discrete_array(:, :), discrete_step)
1266 : END IF
1267 :
1268 : ELSE
1269 :
1270 : ! reset the cell and particle positions
1271 4 : CALL cp_subsys_set(oldsys, cell=cell_old)
1272 40 : DO iatom = 1, nunits_tot(box_number)
1273 148 : particles_old%els(iatom)%r(1:3) = r_old(1:3, iatom)
1274 : END DO
1275 :
1276 : END IF
1277 :
1278 : ! deallocate some stuff
1279 34 : DEALLOCATE (r)
1280 34 : CALL cell_release(cell_test)
1281 34 : CALL cell_release(cell_old)
1282 :
1283 : ! end the timing
1284 34 : CALL timestop(handle)
1285 :
1286 68 : END SUBROUTINE mc_volume_move
1287 :
1288 : ! **************************************************************************************************
1289 : !> \brief alters the length of a random bond for the given molecule, using
1290 : !> a mass weighted scheme so the lightest atoms move the most
1291 : !> \param r_old the initial coordinates of all molecules in the system
1292 : !> \param r_new the new coordinates of all molecules in the system
1293 : !> \param mc_par the mc parameters for the force env
1294 : !> \param molecule_type the molecule type that we're moving
1295 : !> \param molecule_kind the structure containing the molecule information
1296 : !> \param dis_length the ratio of the new bond length to the old bond length,
1297 : !> used in the acceptance rule
1298 : !> \param particles the particle_list_type for all particles in the force_env..
1299 : !> used to grab the mass of each atom
1300 : !> \param rng_stream the random number stream that we draw from
1301 : !>
1302 : !> This subroutine is written to be parallel.
1303 : !> \author MJM
1304 : ! **************************************************************************************************
1305 312 : SUBROUTINE change_bond_length(r_old, r_new, mc_par, molecule_type, molecule_kind, &
1306 : dis_length, particles, rng_stream)
1307 :
1308 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_old
1309 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: r_new
1310 : TYPE(mc_simpar_type), POINTER :: mc_par
1311 : INTEGER, INTENT(IN) :: molecule_type
1312 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1313 : REAL(KIND=dp), INTENT(OUT) :: dis_length
1314 : TYPE(particle_list_type), POINTER :: particles
1315 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1316 :
1317 : CHARACTER(len=*), PARAMETER :: routineN = 'change_bond_length'
1318 :
1319 : INTEGER :: bond_number, handle, i, iatom, ibond, &
1320 : ipart, natom, nbond, source
1321 312 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_b, counter
1322 312 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity
1323 312 : INTEGER, DIMENSION(:), POINTER :: nunits
1324 : LOGICAL :: ionode
1325 312 : REAL(dp), DIMENSION(:), POINTER :: rmbond
1326 : REAL(KIND=dp) :: atom_mass, mass_a, mass_b, new_length_a, &
1327 : new_length_b, old_length, rand
1328 : REAL(KIND=dp), DIMENSION(1:3) :: bond_a, bond_b
1329 312 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
1330 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1331 : TYPE(mp_comm_type) :: group
1332 :
1333 : ! begin the timing of the subroutine
1334 :
1335 312 : CALL timeset(routineN, handle)
1336 :
1337 312 : NULLIFY (rmbond, mc_molecule_info)
1338 :
1339 : ! get some stuff from mc_par
1340 : CALL get_mc_par(mc_par, mc_molecule_info=mc_molecule_info, source=source, &
1341 312 : group=group, rmbond=rmbond, ionode=ionode)
1342 312 : CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
1343 :
1344 : ! copy the incoming coordinates so we can change them
1345 1248 : DO ipart = 1, nunits(molecule_type)
1346 4056 : r_new(1:3, ipart) = r_old(1:3, ipart)
1347 : END DO
1348 :
1349 : ! pick which bond in the molecule at random
1350 312 : IF (ionode) THEN
1351 156 : rand = rng_stream%next()
1352 : END IF
1353 312 : CALL group%bcast(rand, source)
1354 : CALL get_molecule_kind(molecule_kind, natom=natom, nbond=nbond, &
1355 312 : bond_list=bond_list)
1356 312 : bond_number = CEILING(rand*REAL(nbond, dp))
1357 :
1358 936 : ALLOCATE (connection(1:natom, 1:2))
1359 : ! assume at most six bonds per atom
1360 936 : ALLOCATE (connectivity(1:6, 1:natom))
1361 936 : ALLOCATE (counter(1:natom))
1362 624 : ALLOCATE (atom_a(1:natom))
1363 624 : ALLOCATE (atom_b(1:natom))
1364 2808 : connection(:, :) = 0
1365 6864 : connectivity(:, :) = 0
1366 1248 : counter(:) = 0
1367 1248 : atom_a(:) = 0
1368 1248 : atom_b(:) = 0
1369 :
1370 : ! now we need to find a list of atoms that each atom in this bond is connected
1371 : ! to
1372 1248 : DO iatom = 1, natom
1373 3120 : DO ibond = 1, nbond
1374 2808 : IF (bond_list(ibond)%a == iatom) THEN
1375 624 : counter(iatom) = counter(iatom) + 1
1376 624 : connectivity(counter(iatom), iatom) = bond_list(ibond)%b
1377 1248 : ELSEIF (bond_list(ibond)%b == iatom) THEN
1378 624 : counter(iatom) = counter(iatom) + 1
1379 624 : connectivity(counter(iatom), iatom) = bond_list(ibond)%a
1380 : END IF
1381 : END DO
1382 : END DO
1383 :
1384 : ! now I need to do a depth first search to figure out which atoms are on atom a's
1385 : ! side and which are on atom b's
1386 1248 : atom_a(:) = 0
1387 312 : atom_a(bond_list(bond_number)%a) = 1
1388 : CALL depth_first_search(bond_list(bond_number)%a, bond_list(bond_number)%b, &
1389 312 : connectivity(:, :), atom_a(:))
1390 1248 : atom_b(:) = 0
1391 312 : atom_b(bond_list(bond_number)%b) = 1
1392 : CALL depth_first_search(bond_list(bond_number)%b, bond_list(bond_number)%a, &
1393 312 : connectivity(:, :), atom_b(:))
1394 :
1395 : ! now figure out the masses of the various sides, so we can weight how far we move each
1396 : ! group of atoms
1397 312 : mass_a = 0.0_dp
1398 312 : mass_b = 0.0_dp
1399 1248 : DO iatom = 1, natom
1400 : CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
1401 936 : mass=atom_mass)
1402 1248 : IF (atom_a(iatom) == 1) THEN
1403 624 : mass_a = mass_a + atom_mass
1404 : ELSE
1405 312 : mass_b = mass_b + atom_mass
1406 : END IF
1407 : END DO
1408 :
1409 : ! choose a displacement
1410 312 : IF (ionode) rand = rng_stream%next()
1411 312 : CALL group%bcast(rand, source)
1412 :
1413 312 : dis_length = rmbond(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
1414 :
1415 : ! find the bond vector that atom a will be moving
1416 1248 : DO i = 1, 3
1417 : bond_a(i) = r_new(i, bond_list(bond_number)%a) - &
1418 936 : r_new(i, bond_list(bond_number)%b)
1419 1248 : bond_b(i) = -bond_a(i)
1420 : END DO
1421 :
1422 : ! notice we weight by the opposite masses...therefore lighter segments
1423 : ! will move further
1424 1248 : old_length = SQRT(DOT_PRODUCT(bond_a, bond_a))
1425 312 : new_length_a = dis_length*mass_b/(mass_a + mass_b)
1426 312 : new_length_b = dis_length*mass_a/(mass_a + mass_b)
1427 :
1428 1248 : DO i = 1, 3
1429 936 : bond_a(i) = bond_a(i)/old_length*new_length_a
1430 1248 : bond_b(i) = bond_b(i)/old_length*new_length_b
1431 : END DO
1432 :
1433 1248 : DO iatom = 1, natom
1434 1248 : IF (atom_a(iatom) == 1) THEN
1435 624 : r_new(1, iatom) = r_new(1, iatom) + bond_a(1)
1436 624 : r_new(2, iatom) = r_new(2, iatom) + bond_a(2)
1437 624 : r_new(3, iatom) = r_new(3, iatom) + bond_a(3)
1438 : ELSE
1439 312 : r_new(1, iatom) = r_new(1, iatom) + bond_b(1)
1440 312 : r_new(2, iatom) = r_new(2, iatom) + bond_b(2)
1441 312 : r_new(3, iatom) = r_new(3, iatom) + bond_b(3)
1442 : END IF
1443 : END DO
1444 :
1445 : ! correct the value of dis_length for the acceptance rule
1446 312 : dis_length = (old_length + dis_length)/old_length
1447 :
1448 312 : DEALLOCATE (connection)
1449 312 : DEALLOCATE (connectivity)
1450 312 : DEALLOCATE (counter)
1451 312 : DEALLOCATE (atom_a)
1452 312 : DEALLOCATE (atom_b)
1453 : ! end the timing
1454 312 : CALL timestop(handle)
1455 :
1456 624 : END SUBROUTINE change_bond_length
1457 :
1458 : ! **************************************************************************************************
1459 : !> \brief Alters the magnitude of a random angle in a molecule centered on atom C
1460 : !> (connected to atoms A and B). Atoms A and B are moved amounts related
1461 : !> to their masses (and masses of all connecting atoms), so that heavier
1462 : !> segments are moved less.
1463 : !> \param r_old the initial coordinates of all molecules in the system
1464 : !> \param r_new the new coordinates of all molecules in the system
1465 : !> \param mc_par the mc parameters for the force env
1466 : !> \param molecule_type the type of molecule we're playing with
1467 : !> \param molecule_kind the structure containing the molecule information
1468 : !> \param particles the particle_list_type for all particles in the force_env...
1469 : !> used to grab the mass of each atom
1470 : !> \param rng_stream the random number stream that we draw from
1471 : !> \author MJM
1472 : ! **************************************************************************************************
1473 194 : SUBROUTINE change_bond_angle(r_old, r_new, mc_par, molecule_type, molecule_kind, &
1474 : particles, rng_stream)
1475 :
1476 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_old
1477 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: r_new
1478 : TYPE(mc_simpar_type), POINTER :: mc_par
1479 : INTEGER, INTENT(IN) :: molecule_type
1480 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1481 : TYPE(particle_list_type), POINTER :: particles
1482 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1483 :
1484 : CHARACTER(len=*), PARAMETER :: routineN = 'change_bond_angle'
1485 :
1486 : INTEGER :: bend_number, handle, i, iatom, ibond, &
1487 : ipart, natom, nbend, nbond, source
1488 194 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_c, counter
1489 194 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity
1490 194 : INTEGER, DIMENSION(:), POINTER :: nunits
1491 : LOGICAL :: ionode
1492 194 : REAL(dp), DIMENSION(:), POINTER :: rmangle
1493 : REAL(KIND=dp) :: atom_mass, bis_length, dis_angle, dis_angle_a, dis_angle_c, mass_a, mass_c, &
1494 : new_angle_a, new_angle_c, old_angle, old_length_a, old_length_c, rand, temp_length
1495 : REAL(KIND=dp), DIMENSION(1:3) :: bisector, bond_a, bond_c, cross_prod, &
1496 : cross_prod_plane, temp
1497 194 : TYPE(bend_type), DIMENSION(:), POINTER :: bend_list
1498 194 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
1499 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1500 : TYPE(mp_comm_type) :: group
1501 :
1502 : ! begin the timing of the subroutine
1503 :
1504 194 : CALL timeset(routineN, handle)
1505 :
1506 194 : NULLIFY (bend_list, bond_list, rmangle, mc_molecule_info)
1507 :
1508 : ! get some stuff from mc_par
1509 : CALL get_mc_par(mc_par, rmangle=rmangle, source=source, &
1510 194 : group=group, ionode=ionode, mc_molecule_info=mc_molecule_info)
1511 194 : CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
1512 :
1513 : ! copy the incoming coordinates so we can change them
1514 776 : DO ipart = 1, nunits(molecule_type)
1515 2522 : r_new(1:3, ipart) = r_old(1:3, ipart)
1516 : END DO
1517 :
1518 : ! pick which bond in the molecule at random
1519 194 : IF (ionode) THEN
1520 97 : rand = rng_stream%next()
1521 : END IF
1522 194 : CALL group%bcast(rand, source)
1523 : CALL get_molecule_kind(molecule_kind, natom=natom, nbend=nbend, &
1524 194 : bend_list=bend_list, bond_list=bond_list, nbond=nbond)
1525 194 : bend_number = CEILING(rand*REAL(nbend, dp))
1526 :
1527 582 : ALLOCATE (connection(1:natom, 1:2))
1528 : ! assume at most six bonds per atom
1529 582 : ALLOCATE (connectivity(1:6, 1:natom))
1530 582 : ALLOCATE (counter(1:natom))
1531 388 : ALLOCATE (atom_a(1:natom))
1532 388 : ALLOCATE (atom_c(1:natom))
1533 1746 : connection(:, :) = 0
1534 4268 : connectivity(:, :) = 0
1535 776 : counter(:) = 0
1536 776 : atom_a(:) = 0
1537 776 : atom_c(:) = 0
1538 :
1539 : ! now we need to find a list of atoms that each atom in this bond is connected
1540 : ! to
1541 776 : DO iatom = 1, natom
1542 1940 : DO ibond = 1, nbond
1543 1746 : IF (bond_list(ibond)%a == iatom) THEN
1544 388 : counter(iatom) = counter(iatom) + 1
1545 388 : connectivity(counter(iatom), iatom) = bond_list(ibond)%b
1546 776 : ELSEIF (bond_list(ibond)%b == iatom) THEN
1547 388 : counter(iatom) = counter(iatom) + 1
1548 388 : connectivity(counter(iatom), iatom) = bond_list(ibond)%a
1549 : END IF
1550 : END DO
1551 : END DO
1552 :
1553 : ! now I need to do a depth first search to figure out which atoms are on atom a's
1554 : ! side and which are on atom c's
1555 776 : atom_a(:) = 0
1556 194 : atom_a(bend_list(bend_number)%a) = 1
1557 : CALL depth_first_search(bend_list(bend_number)%a, bend_list(bend_number)%b, &
1558 194 : connectivity(:, :), atom_a(:))
1559 776 : atom_c(:) = 0
1560 194 : atom_c(bend_list(bend_number)%c) = 1
1561 : CALL depth_first_search(bend_list(bend_number)%c, bend_list(bend_number)%b, &
1562 194 : connectivity(:, :), atom_c(:))
1563 :
1564 : ! now figure out the masses of the various sides, so we can weight how far we move each
1565 : ! group of atoms
1566 194 : mass_a = 0.0_dp
1567 194 : mass_c = 0.0_dp
1568 776 : DO iatom = 1, natom
1569 : CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
1570 582 : mass=atom_mass)
1571 582 : IF (atom_a(iatom) == 1) mass_a = mass_a + atom_mass
1572 1358 : IF (atom_c(iatom) == 1) mass_c = mass_c + atom_mass
1573 : END DO
1574 :
1575 : ! choose a displacement
1576 194 : IF (ionode) rand = rng_stream%next()
1577 194 : CALL group%bcast(rand, source)
1578 :
1579 194 : dis_angle = rmangle(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
1580 :
1581 : ! need to find the A-B-C bisector
1582 :
1583 : ! this going to be tough...we need to find the plane of the A-B-C bond and only shift
1584 : ! that component for all atoms connected to A and C...otherwise we change other
1585 : ! internal degrees of freedom
1586 :
1587 : ! find the bond vectors
1588 776 : DO i = 1, 3
1589 : bond_a(i) = r_new(i, bend_list(bend_number)%a) - &
1590 582 : r_new(i, bend_list(bend_number)%b)
1591 : bond_c(i) = r_new(i, bend_list(bend_number)%c) - &
1592 776 : r_new(i, bend_list(bend_number)%b)
1593 : END DO
1594 776 : old_length_a = SQRT(DOT_PRODUCT(bond_a, bond_a))
1595 776 : old_length_c = SQRT(DOT_PRODUCT(bond_c, bond_c))
1596 776 : old_angle = ACOS(DOT_PRODUCT(bond_a, bond_c)/(old_length_a*old_length_c))
1597 :
1598 776 : DO i = 1, 3
1599 : bisector(i) = bond_a(i)/old_length_a + & ! not yet normalized
1600 776 : bond_c(i)/old_length_c
1601 : END DO
1602 776 : bis_length = SQRT(DOT_PRODUCT(bisector, bisector))
1603 776 : bisector(1:3) = bisector(1:3)/bis_length
1604 :
1605 : ! now we need to find the cross product of the B-A and B-C vectors and normalize
1606 : ! it, so we have a vector that defines the bend plane
1607 194 : cross_prod(1) = bond_a(2)*bond_c(3) - bond_a(3)*bond_c(2)
1608 194 : cross_prod(2) = bond_a(3)*bond_c(1) - bond_a(1)*bond_c(3)
1609 194 : cross_prod(3) = bond_a(1)*bond_c(2) - bond_a(2)*bond_c(1)
1610 1358 : cross_prod(1:3) = cross_prod(1:3)/SQRT(DOT_PRODUCT(cross_prod, cross_prod))
1611 :
1612 : ! we have two axis of a coordinate system...let's get the third
1613 194 : cross_prod_plane(1) = cross_prod(2)*bisector(3) - cross_prod(3)*bisector(2)
1614 194 : cross_prod_plane(2) = cross_prod(3)*bisector(1) - cross_prod(1)*bisector(3)
1615 194 : cross_prod_plane(3) = cross_prod(1)*bisector(2) - cross_prod(2)*bisector(1)
1616 : cross_prod_plane(1:3) = cross_prod_plane(1:3)/ &
1617 1358 : SQRT(DOT_PRODUCT(cross_prod_plane, cross_prod_plane))
1618 :
1619 : ! now bisector is x, cross_prod_plane is the y vector (pointing towards c),
1620 : ! and cross_prod is z
1621 : ! shift the molecule so that atom b is at the origin
1622 776 : DO iatom = 1, natom
1623 : r_new(1:3, iatom) = r_new(1:3, iatom) - &
1624 2522 : r_old(1:3, bend_list(bend_number)%b)
1625 : END DO
1626 :
1627 : ! figure out how much we move each side, since we're mass-weighting, by the
1628 : ! opposite masses, so lighter moves farther..this angle is the angle between
1629 : ! the bond vector BA or BC and the bisector
1630 194 : dis_angle_a = dis_angle*mass_c/(mass_a + mass_c)
1631 194 : dis_angle_c = dis_angle*mass_a/(mass_a + mass_c)
1632 :
1633 : ! now loop through all the atoms, moving the ones that are connected to a or c
1634 776 : DO iatom = 1, natom
1635 : ! subtract out the z component (perpendicular to the angle plane)
1636 : temp(1:3) = r_new(1:3, iatom) - &
1637 : DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
1638 4074 : cross_prod(1:3)
1639 2328 : temp_length = SQRT(DOT_PRODUCT(temp, temp))
1640 :
1641 : ! we can now compute all three components of the new bond vector along the
1642 : ! axis defined above
1643 776 : IF (atom_a(iatom) == 1) THEN
1644 :
1645 : ! if the y-coordinate is less than zero, we need to switch the sign when we make the vector,
1646 : ! as the angle computed by the dot product can't distinguish between that
1647 776 : IF (DOT_PRODUCT(cross_prod_plane(1:3), r_new(1:3, iatom)) &
1648 : .LT. 0.0_dp) THEN
1649 :
1650 : ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is
1651 : new_angle_a = ACOS(DOT_PRODUCT(bisector, temp(1:3))/ &
1652 776 : (temp_length)) + dis_angle_a
1653 :
1654 : r_new(1:3, iatom) = COS(new_angle_a)*temp_length*bisector(1:3) - &
1655 : SIN(new_angle_a)*temp_length*cross_prod_plane(1:3) + &
1656 : DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
1657 1358 : cross_prod(1:3)
1658 : ELSE
1659 :
1660 : ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is
1661 : new_angle_a = ACOS(DOT_PRODUCT(bisector, temp(1:3))/ &
1662 0 : (temp_length)) - dis_angle_a
1663 :
1664 : r_new(1:3, iatom) = COS(new_angle_a)*temp_length*bisector(1:3) + &
1665 : SIN(new_angle_a)*temp_length*cross_prod_plane(1:3) + &
1666 : DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
1667 0 : cross_prod(1:3)
1668 : END IF
1669 :
1670 388 : ELSEIF (atom_c(iatom) == 1) THEN
1671 :
1672 : ! if the y-coordinate is less than zero, we need to switch the sign when we make the vector,
1673 : ! as the angle computed by the dot product can't distinguish between that
1674 776 : IF (DOT_PRODUCT(cross_prod_plane(1:3), r_new(1:3, iatom)) &
1675 : .LT. 0.0_dp) THEN
1676 : ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is
1677 : new_angle_c = ACOS(DOT_PRODUCT(bisector(1:3), temp(1:3))/ &
1678 0 : (temp_length)) - dis_angle_c
1679 :
1680 : r_new(1:3, iatom) = COS(new_angle_c)*temp_length*bisector(1:3) - &
1681 : SIN(new_angle_c)*temp_length*cross_prod_plane(1:3) + &
1682 : DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
1683 0 : cross_prod(1:3)
1684 : ELSE
1685 : new_angle_c = ACOS(DOT_PRODUCT(bisector(1:3), temp(1:3))/ &
1686 776 : (temp_length)) + dis_angle_c
1687 :
1688 : r_new(1:3, iatom) = COS(new_angle_c)*temp_length*bisector(1:3) + &
1689 : SIN(new_angle_c)*temp_length*cross_prod_plane(1:3) + &
1690 : DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
1691 1358 : cross_prod(1:3)
1692 : END IF
1693 : END IF
1694 :
1695 : END DO
1696 :
1697 776 : DO iatom = 1, natom
1698 : r_new(1:3, iatom) = r_new(1:3, iatom) + &
1699 2522 : r_old(1:3, bend_list(bend_number)%b)
1700 : END DO
1701 :
1702 : ! deallocate some stuff
1703 194 : DEALLOCATE (connection)
1704 194 : DEALLOCATE (connectivity)
1705 194 : DEALLOCATE (counter)
1706 194 : DEALLOCATE (atom_a)
1707 194 : DEALLOCATE (atom_c)
1708 :
1709 : ! end the timing
1710 194 : CALL timestop(handle)
1711 :
1712 388 : END SUBROUTINE change_bond_angle
1713 :
1714 : ! **************************************************************************************************
1715 : !> \brief Alters a dihedral (A-B-C-D) in the molecule so that all other internal
1716 : !> degrees of freedom remain the same. If other dihedrals are centered
1717 : !> on B-C, they rotate as well to keep the relationship between the
1718 : !> dihedrals the same. Atoms A and D are moved amounts related to their
1719 : !> masses (and masses of all connecting atoms), so that heavier segments
1720 : !> are moved less. All atoms except B and C are rotated around the
1721 : !> B-C bond vector (B and C are not moved).
1722 : !> \param r_old the initial coordinates of all molecules in the system
1723 : !> \param r_new the new coordinates of all molecules in the system
1724 : !> \param mc_par the mc parameters for the force env
1725 : !> \param molecule_type the type of molecule we're playing with
1726 : !> \param molecule_kind the structure containing the molecule information
1727 : !> \param particles the particle_list_type for all particles in the force_env..
1728 : !> used to grab the mass of each atom
1729 : !> \param rng_stream the random number stream that we draw from
1730 : !> \author MJM
1731 : ! **************************************************************************************************
1732 0 : SUBROUTINE change_dihedral(r_old, r_new, mc_par, molecule_type, molecule_kind, &
1733 : particles, rng_stream)
1734 :
1735 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_old
1736 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: r_new
1737 : TYPE(mc_simpar_type), POINTER :: mc_par
1738 : INTEGER, INTENT(IN) :: molecule_type
1739 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1740 : TYPE(particle_list_type), POINTER :: particles
1741 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1742 :
1743 : CHARACTER(len=*), PARAMETER :: routineN = 'change_dihedral'
1744 :
1745 : INTEGER :: handle, i, iatom, ibond, ipart, natom, &
1746 : nbond, ntorsion, source, torsion_number
1747 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_d, counter
1748 0 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity
1749 0 : INTEGER, DIMENSION(:), POINTER :: nunits
1750 : LOGICAL :: ionode
1751 0 : REAL(dp), DIMENSION(:), POINTER :: rmdihedral
1752 : REAL(KIND=dp) :: atom_mass, dis_angle, dis_angle_a, &
1753 : dis_angle_d, mass_a, mass_d, &
1754 : old_length_a, rand, u, v, w, x, y, z
1755 : REAL(KIND=dp), DIMENSION(1:3) :: bond_a, temp
1756 0 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
1757 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1758 : TYPE(mp_comm_type) :: group
1759 0 : TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list
1760 :
1761 : ! begin the timing of the subroutine
1762 :
1763 0 : CALL timeset(routineN, handle)
1764 :
1765 0 : NULLIFY (rmdihedral, torsion_list, bond_list, mc_molecule_info)
1766 :
1767 : ! get some stuff from mc_par
1768 : CALL get_mc_par(mc_par, rmdihedral=rmdihedral, &
1769 : source=source, group=group, ionode=ionode, &
1770 0 : mc_molecule_info=mc_molecule_info)
1771 0 : CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
1772 :
1773 : ! copy the incoming coordinates so we can change them
1774 0 : DO ipart = 1, nunits(molecule_type)
1775 0 : r_new(1:3, ipart) = r_old(1:3, ipart)
1776 : END DO
1777 :
1778 : ! pick which bond in the molecule at random
1779 0 : IF (ionode) THEN
1780 0 : rand = rng_stream%next()
1781 : ! CALL RANDOM_NUMBER(rand)
1782 : END IF
1783 0 : CALL group%bcast(rand, source)
1784 : CALL get_molecule_kind(molecule_kind, natom=natom, &
1785 : bond_list=bond_list, nbond=nbond, &
1786 0 : ntorsion=ntorsion, torsion_list=torsion_list)
1787 0 : torsion_number = CEILING(rand*REAL(ntorsion, dp))
1788 :
1789 0 : ALLOCATE (connection(1:natom, 1:2))
1790 : ! assume at most six bonds per atom
1791 0 : ALLOCATE (connectivity(1:6, 1:natom))
1792 0 : ALLOCATE (counter(1:natom))
1793 0 : ALLOCATE (atom_a(1:natom))
1794 0 : ALLOCATE (atom_d(1:natom))
1795 0 : connection(:, :) = 0
1796 0 : connectivity(:, :) = 0
1797 0 : counter(:) = 0
1798 0 : atom_a(:) = 0
1799 0 : atom_d(:) = 0
1800 :
1801 : ! now we need to find a list of atoms that each atom in this bond is connected
1802 : ! to
1803 0 : DO iatom = 1, natom
1804 0 : DO ibond = 1, nbond
1805 0 : IF (bond_list(ibond)%a == iatom) THEN
1806 0 : counter(iatom) = counter(iatom) + 1
1807 0 : connectivity(counter(iatom), iatom) = bond_list(ibond)%b
1808 0 : ELSEIF (bond_list(ibond)%b == iatom) THEN
1809 0 : counter(iatom) = counter(iatom) + 1
1810 0 : connectivity(counter(iatom), iatom) = bond_list(ibond)%a
1811 : END IF
1812 : END DO
1813 : END DO
1814 :
1815 : ! now I need to do a depth first search to figure out which atoms are on atom
1816 : ! a's side and which are on atom d's, but remember we're moving all atoms on a's
1817 : ! side of b, including atoms not in a's branch
1818 0 : atom_a(:) = 0
1819 0 : atom_a(torsion_list(torsion_number)%a) = 1
1820 : CALL depth_first_search(torsion_list(torsion_number)%b, &
1821 0 : torsion_list(torsion_number)%c, connectivity(:, :), atom_a(:))
1822 0 : atom_d(:) = 0
1823 0 : atom_d(torsion_list(torsion_number)%d) = 1
1824 : CALL depth_first_search(torsion_list(torsion_number)%c, &
1825 0 : torsion_list(torsion_number)%b, connectivity(:, :), atom_d(:))
1826 :
1827 : ! now figure out the masses of the various sides, so we can weight how far we
1828 : ! move each group of atoms
1829 0 : mass_a = 0.0_dp
1830 0 : mass_d = 0.0_dp
1831 0 : DO iatom = 1, natom
1832 : CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
1833 0 : mass=atom_mass)
1834 0 : IF (atom_a(iatom) == 1) mass_a = mass_a + atom_mass
1835 0 : IF (atom_d(iatom) == 1) mass_d = mass_d + atom_mass
1836 : END DO
1837 :
1838 : ! choose a displacement
1839 0 : IF (ionode) rand = rng_stream%next()
1840 0 : CALL group%bcast(rand, source)
1841 :
1842 0 : dis_angle = rmdihedral(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
1843 :
1844 : ! find the bond vectors, B-C, so we know what to rotate around
1845 0 : DO i = 1, 3
1846 : bond_a(i) = r_new(i, torsion_list(torsion_number)%c) - &
1847 0 : r_new(i, torsion_list(torsion_number)%b)
1848 : END DO
1849 0 : old_length_a = SQRT(DOT_PRODUCT(bond_a, bond_a))
1850 0 : bond_a(1:3) = bond_a(1:3)/old_length_a
1851 :
1852 : ! figure out how much we move each side, since we're mass-weighting, by the
1853 : ! opposite masses, so lighter moves farther...we take the opposite sign of d
1854 : ! so we're not rotating both angles in the same direction
1855 0 : dis_angle_a = dis_angle*mass_d/(mass_a + mass_d)
1856 0 : dis_angle_d = -dis_angle*mass_a/(mass_a + mass_d)
1857 :
1858 0 : DO iatom = 1, natom
1859 :
1860 0 : IF (atom_a(iatom) == 1) THEN
1861 : ! shift the coords so b is at the origin
1862 : r_new(1:3, iatom) = r_new(1:3, iatom) - &
1863 0 : r_new(1:3, torsion_list(torsion_number)%b)
1864 :
1865 : ! multiply by the rotation matrix
1866 0 : u = bond_a(1)
1867 0 : v = bond_a(2)
1868 0 : w = bond_a(3)
1869 0 : x = r_new(1, iatom)
1870 0 : y = r_new(2, iatom)
1871 0 : z = r_new(3, iatom)
1872 : temp(1) = (u*(u*x + v*y + w*z) + (x*(v**2 + w**2) - u*(v*y + w*z))*COS(dis_angle_a) + &
1873 0 : SQRT(u**2 + v**2 + w**2)*(v*z - w*y)*SIN(dis_angle_a))/(u**2 + v**2 + w**2)
1874 : temp(2) = (v*(u*x + v*y + w*z) + (y*(u**2 + w**2) - v*(u*x + w*z))*COS(dis_angle_a) + &
1875 0 : SQRT(u**2 + v**2 + w**2)*(w*x - u*z)*SIN(dis_angle_a))/(u**2 + v**2 + w**2)
1876 : temp(3) = (w*(u*x + v*y + w*z) + (z*(v**2 + u**2) - w*(u*x + v*y))*COS(dis_angle_a) + &
1877 0 : SQRT(u**2 + v**2 + w**2)*(u*y - v*x)*SIN(dis_angle_a))/(u**2 + v**2 + w**2)
1878 :
1879 : ! shift back to the original position
1880 0 : temp(1:3) = temp(1:3) + r_new(1:3, torsion_list(torsion_number)%b)
1881 0 : r_new(1:3, iatom) = temp(1:3)
1882 :
1883 0 : ELSEIF (atom_d(iatom) == 1) THEN
1884 :
1885 : ! shift the coords so c is at the origin
1886 : r_new(1:3, iatom) = r_new(1:3, iatom) - &
1887 0 : r_new(1:3, torsion_list(torsion_number)%c)
1888 :
1889 : ! multiply by the rotation matrix
1890 0 : u = bond_a(1)
1891 0 : v = bond_a(2)
1892 0 : w = bond_a(3)
1893 0 : x = r_new(1, iatom)
1894 0 : y = r_new(2, iatom)
1895 0 : z = r_new(3, iatom)
1896 : temp(1) = (u*(u*x + v*y + w*z) + (x*(v**2 + w**2) - u*(v*y + w*z))*COS(dis_angle_d) + &
1897 0 : SQRT(u**2 + v**2 + w**2)*(v*z - w*y)*SIN(dis_angle_d))/(u**2 + v**2 + w**2)
1898 : temp(2) = (v*(u*x + v*y + w*z) + (y*(u**2 + w**2) - v*(u*x + w*z))*COS(dis_angle_d) + &
1899 0 : SQRT(u**2 + v**2 + w**2)*(w*x - u*z)*SIN(dis_angle_d))/(u**2 + v**2 + w**2)
1900 : temp(3) = (w*(u*x + v*y + w*z) + (z*(v**2 + u**2) - w*(u*x + v*y))*COS(dis_angle_d) + &
1901 0 : SQRT(u**2 + v**2 + w**2)*(u*y - v*x)*SIN(dis_angle_d))/(u**2 + v**2 + w**2)
1902 :
1903 : ! shift back to the original position
1904 0 : temp(1:3) = temp(1:3) + r_new(1:3, torsion_list(torsion_number)%c)
1905 0 : r_new(1:3, iatom) = temp(1:3)
1906 : END IF
1907 : END DO
1908 :
1909 : ! deallocate some stuff
1910 0 : DEALLOCATE (connection)
1911 0 : DEALLOCATE (connectivity)
1912 0 : DEALLOCATE (counter)
1913 0 : DEALLOCATE (atom_a)
1914 0 : DEALLOCATE (atom_d)
1915 :
1916 : ! end the timing
1917 0 : CALL timestop(handle)
1918 :
1919 0 : END SUBROUTINE change_dihedral
1920 :
1921 : ! **************************************************************************************************
1922 : !> \brief performs either a bond or angle change move for a given molecule
1923 : !> \param mc_par the mc parameters for the force env
1924 : !> \param force_env the force environment used in the move
1925 : !> \param bias_env the force environment used to bias the move, if any (it may
1926 : !> be null if lbias=.false. in mc_par)
1927 : !> \param moves the structure that keeps track of how many moves have been
1928 : !> accepted/rejected
1929 : !> \param energy_check the running energy difference between now and the initial
1930 : !> energy
1931 : !> \param r_old the coordinates of force_env before the move
1932 : !> \param old_energy the energy of the force_env before the move
1933 : !> \param start_atom_swap the number of the swap molecule's first atom, assuming the rest of
1934 : !> the atoms follow sequentially
1935 : !> \param target_atom the number of the target atom for swapping
1936 : !> \param molecule_type the molecule type for the atom we're swapping
1937 : !> \param box_number the number of the box we're doing this move in
1938 : !> \param bias_energy_old the biased energy of the system before the move
1939 : !> \param last_bias_energy the last biased energy of the system
1940 : !> \param move_type dictates if we're moving to an "in" or "out" region
1941 : !> \param rng_stream the random number stream that we draw from
1942 : !> \author MJM
1943 : !> \note Designed for parallel.
1944 : ! **************************************************************************************************
1945 0 : SUBROUTINE mc_avbmc_move(mc_par, force_env, bias_env, moves, &
1946 0 : energy_check, r_old, old_energy, start_atom_swap, &
1947 : target_atom, &
1948 : molecule_type, box_number, bias_energy_old, last_bias_energy, &
1949 : move_type, rng_stream)
1950 :
1951 : TYPE(mc_simpar_type), POINTER :: mc_par
1952 : TYPE(force_env_type), POINTER :: force_env, bias_env
1953 : TYPE(mc_moves_type), POINTER :: moves
1954 : REAL(KIND=dp), INTENT(INOUT) :: energy_check
1955 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r_old
1956 : REAL(KIND=dp), INTENT(INOUT) :: old_energy
1957 : INTEGER, INTENT(IN) :: start_atom_swap, target_atom, &
1958 : molecule_type, box_number
1959 : REAL(KIND=dp), INTENT(INOUT) :: bias_energy_old, last_bias_energy
1960 : CHARACTER(LEN=*), INTENT(IN) :: move_type
1961 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1962 :
1963 : CHARACTER(len=*), PARAMETER :: routineN = 'mc_avbmc_move'
1964 :
1965 : INTEGER :: end_mol, handle, ipart, jbox, natom, &
1966 : nswapmoves, source, start_mol
1967 0 : INTEGER, DIMENSION(:), POINTER :: avbmc_atom, mol_type, nunits, nunits_tot
1968 0 : INTEGER, DIMENSION(:, :), POINTER :: nchains
1969 : LOGICAL :: ionode, lbias, ldum, lin, loverlap
1970 0 : REAL(dp), DIMENSION(:), POINTER :: avbmc_rmax, avbmc_rmin, pbias
1971 0 : REAL(dp), DIMENSION(:, :), POINTER :: mass
1972 : REAL(KIND=dp) :: BETA, bias_energy_new, del_quickstep_energy, distance, exp_max_val, &
1973 : exp_min_val, max_val, min_val, new_energy, prefactor, rand, rdum, volume_in, volume_out, &
1974 : w, weight_new, weight_old
1975 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r_new
1976 : REAL(KIND=dp), DIMENSION(1:3) :: abc, RIJ
1977 : TYPE(cell_type), POINTER :: cell
1978 : TYPE(cp_subsys_type), POINTER :: subsys, subsys_force
1979 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1980 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
1981 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1982 : TYPE(mp_comm_type) :: group
1983 : TYPE(particle_list_type), POINTER :: particles, particles_force
1984 :
1985 0 : rdum = 1.0_dp
1986 :
1987 : ! begin the timing of the subroutine
1988 0 : CALL timeset(routineN, handle)
1989 :
1990 : ! get a bunch of stuff from mc_par
1991 : CALL get_mc_par(mc_par, lbias=lbias, &
1992 : BETA=BETA, max_val=max_val, min_val=min_val, exp_max_val=exp_max_val, &
1993 : exp_min_val=exp_min_val, avbmc_atom=avbmc_atom, &
1994 : avbmc_rmin=avbmc_rmin, avbmc_rmax=avbmc_rmax, &
1995 : nswapmoves=nswapmoves, ionode=ionode, source=source, &
1996 0 : group=group, pbias=pbias, mc_molecule_info=mc_molecule_info)
1997 : CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, &
1998 0 : mass=mass, nunits=nunits, nunits_tot=nunits_tot, mol_type=mol_type)
1999 : ! figure out some bounds for mol_type
2000 0 : start_mol = 1
2001 0 : DO jbox = 1, box_number - 1
2002 0 : start_mol = start_mol + SUM(nchains(:, jbox))
2003 : END DO
2004 0 : end_mol = start_mol + SUM(nchains(:, box_number)) - 1
2005 :
2006 : ! nullify some pointers
2007 0 : NULLIFY (particles, subsys, molecule_kinds, molecule_kind, &
2008 0 : particles_force, subsys_force)
2009 :
2010 : ! do some allocation
2011 0 : ALLOCATE (r_new(1:3, 1:nunits_tot(box_number)))
2012 :
2013 : ! now we need to grab and save coordinates, in case we reject
2014 : ! are we biasing this move?
2015 0 : IF (lbias) THEN
2016 :
2017 : ! grab the coordinates
2018 0 : CALL force_env_get(bias_env, cell=cell, subsys=subsys)
2019 : CALL cp_subsys_get(subsys, &
2020 0 : particles=particles, molecule_kinds=molecule_kinds)
2021 0 : molecule_kind => molecule_kinds%els(1)
2022 0 : CALL get_molecule_kind(molecule_kind, natom=natom)
2023 0 : CALL get_cell(cell, abc=abc)
2024 :
2025 : ! save the energy
2026 : ! bias_energy_old=bias_energy
2027 :
2028 : ELSE
2029 :
2030 : ! grab the coordinates
2031 0 : CALL force_env_get(force_env, cell=cell, subsys=subsys)
2032 : CALL cp_subsys_get(subsys, &
2033 0 : particles=particles, molecule_kinds=molecule_kinds)
2034 0 : molecule_kind => molecule_kinds%els(1)
2035 0 : CALL get_molecule_kind(molecule_kind, natom=natom)
2036 0 : CALL get_cell(cell, abc=abc)
2037 :
2038 : END IF
2039 :
2040 : ! let's determine if the molecule to be moved is in the "in" region or the
2041 : ! "out" region of the target
2042 : RIJ(1) = particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(1) - &
2043 : particles%els(target_atom)%r(1) - abc(1)*ANINT( &
2044 : (particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(1) - &
2045 0 : particles%els(target_atom)%r(1))/abc(1))
2046 : RIJ(2) = particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(2) - &
2047 : particles%els(target_atom)%r(2) - abc(2)*ANINT( &
2048 : (particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(2) - &
2049 0 : particles%els(target_atom)%r(2))/abc(2))
2050 : RIJ(3) = particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(3) - &
2051 : particles%els(target_atom)%r(3) - abc(3)*ANINT( &
2052 : (particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(3) - &
2053 0 : particles%els(target_atom)%r(3))/abc(3))
2054 0 : distance = SQRT(RIJ(1)**2 + RIJ(2)**2 + RIJ(3)**2)
2055 0 : IF (distance .LE. avbmc_rmax(molecule_type) .AND. distance .GE. avbmc_rmin(molecule_type)) THEN
2056 : lin = .TRUE.
2057 : ELSE
2058 : lin = .FALSE.
2059 : END IF
2060 :
2061 : ! increment the counter of the particular move we've done
2062 : ! swapping into the "in" region of mol_target
2063 : IF (lin) THEN
2064 0 : IF (move_type == 'in') THEN
2065 : moves%avbmc_inin%attempts = &
2066 0 : moves%avbmc_inin%attempts + 1
2067 : ELSE
2068 : moves%avbmc_inout%attempts = &
2069 0 : moves%avbmc_inout%attempts + 1
2070 : END IF
2071 : ELSE
2072 0 : IF (move_type == 'in') THEN
2073 : moves%avbmc_outin%attempts = &
2074 0 : moves%avbmc_outin%attempts + 1
2075 : ELSE
2076 : moves%avbmc_outout%attempts = &
2077 0 : moves%avbmc_outout%attempts + 1
2078 : END IF
2079 : END IF
2080 :
2081 0 : IF (lbias) THEN
2082 :
2083 0 : IF (move_type == 'in') THEN
2084 :
2085 : ! do CBMC for the old config
2086 : CALL generate_cbmc_swap_config(bias_env, BETA, max_val, min_val, exp_max_val, &
2087 : exp_min_val, nswapmoves, &
2088 : weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2089 : mass(:, molecule_type), ldum, rdum, &
2090 : bias_energy_old, ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2091 : source, group, rng_stream, &
2092 : avbmc_atom=avbmc_atom(molecule_type), &
2093 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='out', &
2094 0 : target_atom=target_atom)
2095 :
2096 : ELSE
2097 :
2098 : ! do CBMC for the old config
2099 : CALL generate_cbmc_swap_config(bias_env, BETA, max_val, min_val, exp_max_val, &
2100 : exp_min_val, nswapmoves, &
2101 : weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2102 : mass(:, molecule_type), ldum, rdum, &
2103 : bias_energy_old, ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2104 : source, group, rng_stream, &
2105 : avbmc_atom=avbmc_atom(molecule_type), &
2106 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='in', &
2107 0 : target_atom=target_atom)
2108 :
2109 : END IF
2110 :
2111 : ! generate the new config
2112 : CALL generate_cbmc_swap_config(bias_env, BETA, max_val, min_val, exp_max_val, &
2113 : exp_min_val, nswapmoves, &
2114 : weight_new, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2115 : mass(:, molecule_type), loverlap, bias_energy_new, &
2116 : bias_energy_old, ionode, .FALSE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2117 : source, group, rng_stream, &
2118 : avbmc_atom=avbmc_atom(molecule_type), &
2119 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type=move_type, &
2120 0 : target_atom=target_atom)
2121 :
2122 : ! the energy that comes out of the above routine is the difference...we want
2123 : ! the real energy for the acceptance rule...we don't do this for the
2124 : ! lbias=.false. case because it doesn't appear in the acceptance rule, and
2125 : ! we compensate in case of acceptance
2126 0 : bias_energy_new = bias_energy_new + bias_energy_old
2127 :
2128 : ELSE
2129 :
2130 0 : IF (move_type == 'in') THEN
2131 :
2132 : ! find the weight of the old config
2133 : CALL generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, &
2134 : exp_min_val, nswapmoves, &
2135 : weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2136 : mass(:, molecule_type), ldum, rdum, old_energy, &
2137 : ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2138 : source, group, rng_stream, &
2139 : avbmc_atom=avbmc_atom(molecule_type), &
2140 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='out', &
2141 0 : target_atom=target_atom)
2142 :
2143 : ELSE
2144 :
2145 : ! find the weight of the old config
2146 : CALL generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, &
2147 : exp_min_val, nswapmoves, &
2148 : weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2149 : mass(:, molecule_type), ldum, rdum, old_energy, &
2150 : ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2151 : source, group, rng_stream, &
2152 : avbmc_atom=avbmc_atom(molecule_type), &
2153 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='in', &
2154 0 : target_atom=target_atom)
2155 :
2156 : END IF
2157 :
2158 : ! generate the new config...do this after, because it changes the force_env
2159 : CALL generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, &
2160 : exp_min_val, nswapmoves, &
2161 : weight_new, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
2162 : mass(:, molecule_type), loverlap, new_energy, old_energy, &
2163 : ionode, .FALSE., mol_type(start_mol:end_mol), nchains(:, box_number), &
2164 : source, group, rng_stream, &
2165 : avbmc_atom=avbmc_atom(molecule_type), &
2166 : rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type=move_type, &
2167 0 : target_atom=target_atom)
2168 :
2169 : END IF
2170 :
2171 0 : IF (loverlap) THEN
2172 0 : DEALLOCATE (r_new)
2173 :
2174 : ! need to reset the old coordinates
2175 0 : IF (lbias) THEN
2176 0 : CALL force_env_get(bias_env, subsys=subsys)
2177 0 : CALL cp_subsys_get(subsys, particles=particles)
2178 : ELSE
2179 0 : CALL force_env_get(force_env, subsys=subsys)
2180 0 : CALL cp_subsys_get(subsys, particles=particles)
2181 : END IF
2182 0 : DO ipart = 1, nunits_tot(box_number)
2183 0 : particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2184 : END DO
2185 :
2186 0 : CALL timestop(handle)
2187 :
2188 : RETURN
2189 : END IF
2190 :
2191 : ! if we're biasing, we need to compute the new energy with the full
2192 : ! potential
2193 0 : IF (lbias) THEN
2194 : ! need to give the force_env the coords from the bias_env
2195 0 : CALL force_env_get(force_env, subsys=subsys_force)
2196 0 : CALL cp_subsys_get(subsys_force, particles=particles_force)
2197 0 : CALL force_env_get(bias_env, subsys=subsys)
2198 0 : CALL cp_subsys_get(subsys, particles=particles)
2199 0 : DO ipart = 1, nunits_tot(box_number)
2200 0 : particles_force%els(ipart)%r(1:3) = particles%els(ipart)%r(1:3)
2201 : END DO
2202 :
2203 : CALL force_env_calc_energy_force(force_env, &
2204 0 : calc_force=.FALSE.)
2205 : CALL force_env_get(force_env, &
2206 0 : potential_energy=new_energy)
2207 :
2208 : END IF
2209 :
2210 0 : volume_in = 4.0_dp/3.0_dp*pi*(avbmc_rmax(molecule_type)**3 - avbmc_rmin(molecule_type)**3)
2211 0 : volume_out = abc(1)*abc(2)*abc(3) - volume_in
2212 :
2213 0 : IF (lin .AND. move_type == 'in' .OR. &
2214 : .NOT. lin .AND. move_type == 'out') THEN
2215 : ! standard Metropolis rule
2216 : prefactor = 1.0_dp
2217 0 : ELSEIF (.NOT. lin .AND. move_type == 'in') THEN
2218 0 : prefactor = (1.0_dp - pbias(molecule_type))*volume_in/(pbias(molecule_type)*volume_out)
2219 : ELSE
2220 0 : prefactor = pbias(molecule_type)*volume_out/((1.0_dp - pbias(molecule_type))*volume_in)
2221 : END IF
2222 :
2223 0 : IF (lbias) THEN
2224 : ! AVBMC with CBMC and a biasing potential...notice that if the biasing
2225 : ! potential equals the quickstep potential, this cancels out to the
2226 : ! acceptance below
2227 : del_quickstep_energy = (-BETA)*(new_energy - old_energy - &
2228 0 : (bias_energy_new - bias_energy_old))
2229 :
2230 0 : IF (del_quickstep_energy .GT. exp_max_val) THEN
2231 0 : del_quickstep_energy = max_val
2232 0 : ELSEIF (del_quickstep_energy .LT. exp_min_val) THEN
2233 : del_quickstep_energy = 0.0_dp
2234 : ELSE
2235 0 : del_quickstep_energy = EXP(del_quickstep_energy)
2236 : END IF
2237 :
2238 0 : w = prefactor*del_quickstep_energy*weight_new/weight_old
2239 :
2240 : ELSE
2241 :
2242 : ! AVBMC with CBMC
2243 0 : w = prefactor*weight_new/weight_old
2244 : END IF
2245 :
2246 : ! check if the move is accepted
2247 0 : IF (w .GE. 1.0E0_dp) THEN
2248 0 : rand = 0.0E0_dp
2249 : ELSE
2250 0 : IF (ionode) rand = rng_stream%next()
2251 0 : CALL group%bcast(rand, source)
2252 : END IF
2253 :
2254 0 : IF (rand .LT. w) THEN
2255 :
2256 : ! accept the move
2257 :
2258 0 : IF (lin) THEN
2259 0 : IF (move_type == 'in') THEN
2260 : moves%avbmc_inin%successes = &
2261 0 : moves%avbmc_inin%successes + 1
2262 : ELSE
2263 : moves%avbmc_inout%successes = &
2264 0 : moves%avbmc_inout%successes + 1
2265 : END IF
2266 : ELSE
2267 0 : IF (move_type == 'in') THEN
2268 : moves%avbmc_outin%successes = &
2269 0 : moves%avbmc_outin%successes + 1
2270 : ELSE
2271 : moves%avbmc_outout%successes = &
2272 0 : moves%avbmc_outout%successes + 1
2273 : END IF
2274 : END IF
2275 :
2276 : ! we need to compensate for the fact that we take the difference in
2277 : ! generate_cbmc_config to keep the exponetials small
2278 0 : IF (.NOT. lbias) THEN
2279 0 : new_energy = new_energy + old_energy
2280 : END IF
2281 :
2282 : ! update energies
2283 0 : energy_check = energy_check + (new_energy - old_energy)
2284 0 : old_energy = new_energy
2285 :
2286 : ! if we're biasing the update the biasing energy
2287 0 : IF (lbias) THEN
2288 : ! need to do this outside of the routine
2289 0 : last_bias_energy = bias_energy_new
2290 0 : bias_energy_old = bias_energy_new
2291 : END IF
2292 :
2293 : ! update coordinates
2294 0 : CALL force_env_get(force_env, subsys=subsys)
2295 0 : CALL cp_subsys_get(subsys, particles=particles)
2296 0 : DO ipart = 1, nunits_tot(box_number)
2297 0 : r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
2298 : END DO
2299 : ELSE
2300 : ! reject the move...need to restore the old coordinates
2301 0 : IF (lbias) THEN
2302 0 : CALL force_env_get(bias_env, subsys=subsys)
2303 0 : CALL cp_subsys_get(subsys, particles=particles)
2304 0 : DO ipart = 1, nunits_tot(box_number)
2305 0 : particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2306 : END DO
2307 0 : CALL cp_subsys_set(subsys, particles=particles)
2308 : END IF
2309 0 : CALL force_env_get(force_env, subsys=subsys)
2310 0 : CALL cp_subsys_get(subsys, particles=particles)
2311 0 : DO ipart = 1, nunits_tot(box_number)
2312 0 : particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2313 : END DO
2314 0 : CALL cp_subsys_set(subsys, particles=particles)
2315 :
2316 : END IF
2317 :
2318 : ! deallocate some stuff
2319 0 : DEALLOCATE (r_new)
2320 : ! end the timing
2321 0 : CALL timestop(handle)
2322 :
2323 0 : END SUBROUTINE mc_avbmc_move
2324 :
2325 : ! **************************************************************************************************
2326 : !> \brief performs a hybrid Monte Carlo move that runs a short MD sequence
2327 : !> \param mc_par the mc parameters for the force env
2328 : !> \param force_env the force environment whose cell we're changing
2329 : !> \param globenv ...
2330 : !> \param moves the structure that keeps track of how many moves have been
2331 : !> accepted/rejected
2332 : !> \param move_updates the structure that keeps track of how many moves have
2333 : !> been accepted/rejected since the last time the displacements
2334 : !> were updated
2335 : !> \param old_energy the energy of the last accepted move involving an
2336 : !> unbiased calculation
2337 : !> \param box_number the box we're changing the volume of
2338 : !> \param energy_check the running total of how much the energy has changed
2339 : !> since the initial configuration
2340 : !> \param r_old the coordinates of the last accepted move involving an
2341 : !> unbiased calculation
2342 : !> \param rng_stream the random number stream that we draw from
2343 : !> \author MJM
2344 : !> \note Designed for parallel use.
2345 : ! **************************************************************************************************
2346 20 : SUBROUTINE mc_hmc_move(mc_par, force_env, globenv, moves, move_updates, &
2347 : old_energy, box_number, &
2348 20 : energy_check, r_old, rng_stream)
2349 :
2350 : TYPE(mc_simpar_type), POINTER :: mc_par
2351 : TYPE(force_env_type), POINTER :: force_env
2352 : TYPE(global_environment_type), POINTER :: globenv
2353 : TYPE(mc_moves_type), POINTER :: moves, move_updates
2354 : REAL(KIND=dp), INTENT(INOUT) :: old_energy
2355 : INTEGER, INTENT(IN) :: box_number
2356 : REAL(KIND=dp), INTENT(INOUT) :: energy_check
2357 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r_old
2358 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
2359 :
2360 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_hmc_move'
2361 :
2362 : INTEGER :: handle, iatom, source
2363 20 : INTEGER, DIMENSION(:), POINTER :: nunits_tot
2364 : LOGICAL :: ionode
2365 : REAL(KIND=dp) :: BETA, energy_term, exp_max_val, &
2366 : exp_min_val, new_energy, rand, value, w
2367 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r
2368 : TYPE(cp_subsys_type), POINTER :: oldsys
2369 : TYPE(mc_ekin_type), POINTER :: hmc_ekin
2370 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
2371 : TYPE(mp_comm_type) :: group
2372 : TYPE(particle_list_type), POINTER :: particles_old
2373 :
2374 : ! begin the timing of the subroutine
2375 :
2376 20 : CALL timeset(routineN, handle)
2377 :
2378 : ! get a bunch of stuff from mc_par
2379 : CALL get_mc_par(mc_par, ionode=ionode, &
2380 : BETA=BETA, exp_max_val=exp_max_val, &
2381 : exp_min_val=exp_min_val, source=source, group=group, &
2382 20 : mc_molecule_info=mc_molecule_info)
2383 20 : CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot)
2384 :
2385 : ! nullify some pointers
2386 20 : NULLIFY (particles_old, oldsys, hmc_ekin)
2387 :
2388 : ! do some allocation
2389 60 : ALLOCATE (r(1:3, 1:nunits_tot(box_number)))
2390 20 : ALLOCATE (hmc_ekin)
2391 :
2392 : ! record the attempt
2393 20 : moves%hmc%attempts = moves%hmc%attempts + 1
2394 20 : move_updates%hmc%attempts = move_updates%hmc%attempts + 1
2395 :
2396 : ! now let's grab the particle positions
2397 20 : CALL force_env_get(force_env, subsys=oldsys)
2398 20 : CALL cp_subsys_get(oldsys, particles=particles_old)
2399 :
2400 : ! save the old coordinates
2401 21700 : DO iatom = 1, nunits_tot(box_number)
2402 86740 : r(1:3, iatom) = particles_old%els(iatom)%r(1:3)
2403 : END DO
2404 :
2405 : ! now run the MD simulation
2406 20 : CALL qs_mol_dyn(force_env, globenv, hmc_e_initial=hmc_ekin%initial_ekin, hmc_e_final=hmc_ekin%final_ekin)
2407 :
2408 : ! get the energy
2409 : CALL force_env_get(force_env, &
2410 20 : potential_energy=new_energy)
2411 :
2412 : ! accept or reject the move
2413 : ! to prevent overflows
2414 20 : energy_term = new_energy + hmc_ekin%final_ekin - old_energy - hmc_ekin%initial_ekin
2415 :
2416 20 : value = -BETA*(energy_term)
2417 20 : IF (value .GT. exp_max_val) THEN
2418 : w = 10.0_dp
2419 20 : ELSEIF (value .LT. exp_min_val) THEN
2420 : w = 0.0_dp
2421 : ELSE
2422 20 : w = EXP(value)
2423 : END IF
2424 :
2425 20 : IF (w .GE. 1.0E0_dp) THEN
2426 10 : w = 1.0E0_dp
2427 10 : rand = 0.0E0_dp
2428 : ELSE
2429 10 : IF (ionode) rand = rng_stream%next()
2430 10 : CALL group%bcast(rand, source)
2431 : END IF
2432 :
2433 20 : IF (rand .LT. w) THEN
2434 :
2435 : ! accept the move
2436 14 : moves%hmc%successes = moves%hmc%successes + 1
2437 14 : move_updates%hmc%successes = move_updates%hmc%successes + 1
2438 :
2439 : ! update energies
2440 14 : energy_check = energy_check + (new_energy - old_energy)
2441 14 : old_energy = new_energy
2442 :
2443 15190 : DO iatom = 1, nunits_tot(box_number)
2444 60718 : r_old(1:3, iatom) = particles_old%els(iatom)%r(1:3)
2445 : END DO
2446 :
2447 : ELSE
2448 :
2449 : ! reset the cell and particle positions
2450 6510 : DO iatom = 1, nunits_tot(box_number)
2451 26022 : particles_old%els(iatom)%r(1:3) = r_old(1:3, iatom)
2452 : END DO
2453 :
2454 : END IF
2455 :
2456 : ! deallocate some stuff
2457 20 : DEALLOCATE (r)
2458 20 : DEALLOCATE (hmc_ekin)
2459 :
2460 : ! end the timing
2461 20 : CALL timestop(handle)
2462 :
2463 40 : END SUBROUTINE mc_hmc_move
2464 :
2465 : ! *****************************************************************************
2466 : !> \brief translates the cluster randomly in either the x,y, or z
2467 : !>direction
2468 : !> \param mc_par the mc parameters for the force env
2469 : !> \param force_env the force environment used in the move
2470 : !> \param bias_env the force environment used to bias the move, if any (it may
2471 : !> be null if lbias=.false. in mc_par)
2472 : !> \param moves the structure that keeps track of how many moves have been
2473 : !> accepted/rejected
2474 : !> \param move_updates the structure that keeps track of how many moves have
2475 : !> been accepted/rejected since the last time the displacements
2476 : !> were updated
2477 : !> \param box_number ...
2478 : !> \param bias_energy the biased energy of the system before the move
2479 : !> \param lreject set to .true. if there is an overlap
2480 : !> \param rng_stream the random number stream that we draw from
2481 : !> \author Himanshu Goel
2482 : !> \note Designed for parallel use.
2483 : ! **************************************************************************************************
2484 :
2485 10 : SUBROUTINE mc_cluster_translation(mc_par, force_env, bias_env, moves, &
2486 : move_updates, box_number, bias_energy, lreject, rng_stream)
2487 :
2488 : TYPE(mc_simpar_type), POINTER :: mc_par
2489 : TYPE(force_env_type), POINTER :: force_env, bias_env
2490 : TYPE(mc_moves_type), POINTER :: moves, move_updates
2491 : INTEGER, INTENT(IN) :: box_number
2492 : REAL(KIND=dp), INTENT(INOUT) :: bias_energy
2493 : LOGICAL, INTENT(OUT) :: lreject
2494 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
2495 :
2496 : CHARACTER(len=*), PARAMETER :: routineN = 'mc_cluster_translation'
2497 :
2498 : INTEGER :: cstart, end_mol, handle, imol, ipart, iparticle, iunit, jbox, jpart, junit, &
2499 : move_direction, nend, nunit, source, start_mol, total_clus, total_clusafmo
2500 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: cluster
2501 10 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
2502 10 : INTEGER, DIMENSION(:, :), POINTER :: nchains
2503 : LOGICAL :: ionode, lbias, loverlap
2504 : REAL(KIND=dp) :: BETA, bias_energy_new, bias_energy_old, &
2505 : dis_mol, exp_max_val, exp_min_val, &
2506 : rand, rmcltrans, value, w
2507 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r_old
2508 : TYPE(cp_subsys_type), POINTER :: subsys
2509 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
2510 : TYPE(mp_comm_type) :: group
2511 : TYPE(particle_list_type), POINTER :: particles
2512 :
2513 : ! *** Local Counters ***
2514 : ! begin the timing of the subroutine
2515 :
2516 10 : CALL timeset(routineN, handle)
2517 :
2518 : ! nullify some pointers
2519 10 : NULLIFY (particles, subsys)
2520 :
2521 : ! get a bunch of stuff from mc_par
2522 : CALL get_mc_par(mc_par, lbias=lbias, &
2523 : BETA=BETA, exp_max_val=exp_max_val, &
2524 : exp_min_val=exp_min_val, rmcltrans=rmcltrans, ionode=ionode, source=source, &
2525 10 : group=group, mc_molecule_info=mc_molecule_info)
2526 : CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot, &
2527 10 : nchains=nchains, nunits=nunits, mol_type=mol_type)
2528 :
2529 : ! find out some bounds for mol_type
2530 10 : start_mol = 1
2531 10 : DO jbox = 1, box_number - 1
2532 10 : start_mol = start_mol + SUM(nchains(:, jbox))
2533 : END DO
2534 20 : end_mol = start_mol + SUM(nchains(:, box_number)) - 1
2535 :
2536 : ! do some allocation
2537 30 : ALLOCATE (r_old(1:3, 1:nunits_tot(box_number)))
2538 :
2539 : ! Allocating cluster matrix size
2540 20 : nend = SUM(nchains(:, box_number))
2541 40 : ALLOCATE (cluster(nend, nend))
2542 60 : DO ipart = 1, nend
2543 310 : DO jpart = 1, nend
2544 300 : cluster(ipart, jpart) = 0
2545 : END DO
2546 : END DO
2547 :
2548 : ! Get cluster information in cluster matrix from cluster_search subroutine
2549 10 : IF (lbias) THEN
2550 : CALL cluster_search(mc_par, bias_env, cluster, nchains(:, box_number), &
2551 0 : nunits, mol_type(start_mol:end_mol), total_clus)
2552 : ELSE
2553 : CALL cluster_search(mc_par, force_env, cluster, nchains(:, box_number), &
2554 10 : nunits, mol_type(start_mol:end_mol), total_clus)
2555 : END IF
2556 :
2557 10 : IF (lbias) THEN
2558 :
2559 : ! grab the coordinates
2560 0 : CALL force_env_get(bias_env, subsys=subsys)
2561 0 : CALL cp_subsys_get(subsys, particles=particles)
2562 :
2563 : ! save the coordinates
2564 0 : DO ipart = 1, nunits_tot(box_number)
2565 0 : r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
2566 : END DO
2567 :
2568 : ! save the energy
2569 0 : bias_energy_old = bias_energy
2570 : ELSE
2571 :
2572 : ! grab the coordinates
2573 10 : CALL force_env_get(force_env, subsys=subsys)
2574 10 : CALL cp_subsys_get(subsys, particles=particles)
2575 : END IF
2576 :
2577 : ! record the attempt
2578 10 : moves%cltrans%attempts = moves%cltrans%attempts + 1
2579 10 : move_updates%cltrans%attempts = move_updates%cltrans%attempts + 1
2580 10 : moves%bias_cltrans%attempts = moves%bias_cltrans%attempts + 1
2581 10 : move_updates%bias_cltrans%attempts = move_updates%bias_cltrans%attempts + 1
2582 10 : IF (.NOT. lbias) THEN
2583 10 : moves%cltrans%qsuccesses = moves%cltrans%qsuccesses + 1
2584 10 : move_updates%cltrans%qsuccesses = move_updates%cltrans%qsuccesses + 1
2585 10 : moves%bias_cltrans%qsuccesses = moves%bias_cltrans%qsuccesses + 1
2586 10 : move_updates%bias_cltrans%qsuccesses = move_updates%bias_cltrans%qsuccesses + 1
2587 : END IF
2588 :
2589 : ! call a random number to figure out which direction we're moving
2590 10 : IF (ionode) rand = rng_stream%next()
2591 10 : CALL group%bcast(rand, source)
2592 10 : move_direction = INT(3*rand) + 1
2593 :
2594 : ! call a random number to figure out how far we're moving
2595 10 : IF (ionode) rand = rng_stream%next()
2596 10 : CALL group%bcast(rand, source)
2597 10 : dis_mol = rmcltrans*(rand - 0.5E0_dp)*2.0E0_dp
2598 :
2599 : ! choosing cluster
2600 10 : IF (ionode) rand = rng_stream%next()
2601 10 : CALL group%bcast(rand, source)
2602 10 : jpart = INT(1 + rand*total_clus)
2603 :
2604 : ! do the cluster move
2605 60 : DO cstart = 1, nend
2606 50 : imol = 0
2607 60 : IF (cluster(jpart, cstart) .NE. 0) THEN
2608 34 : imol = cluster(jpart, cstart)
2609 : iunit = 1
2610 34 : DO ipart = 1, imol - 1
2611 24 : nunit = nunits(mol_type(ipart + start_mol - 1))
2612 34 : iunit = iunit + nunit
2613 : END DO
2614 10 : nunit = nunits(mol_type(imol + start_mol - 1))
2615 10 : junit = iunit + nunit - 1
2616 30 : DO iparticle = iunit, junit
2617 : particles%els(iparticle)%r(move_direction) = &
2618 30 : particles%els(iparticle)%r(move_direction) + dis_mol
2619 : END DO
2620 : END IF
2621 : END DO
2622 10 : CALL cp_subsys_set(subsys, particles=particles)
2623 :
2624 : !Make cluster matrix null
2625 60 : DO ipart = 1, nend
2626 310 : DO jpart = 1, nend
2627 300 : cluster(ipart, jpart) = 0
2628 : END DO
2629 : END DO
2630 :
2631 : ! checking the number of cluster are same or got changed after cluster translation move
2632 10 : IF (lbias) THEN
2633 : CALL cluster_search(mc_par, bias_env, cluster, nchains(:, box_number), &
2634 0 : nunits, mol_type(start_mol:end_mol), total_clusafmo)
2635 : ELSE
2636 : CALL cluster_search(mc_par, force_env, cluster, nchains(:, box_number), &
2637 10 : nunits, mol_type(start_mol:end_mol), total_clusafmo)
2638 : END IF
2639 :
2640 : ! figure out if there is any overlap...need the number of the molecule
2641 10 : lreject = .FALSE.
2642 10 : IF (lbias) THEN
2643 : CALL check_for_overlap(bias_env, nchains(:, box_number), &
2644 0 : nunits(:), loverlap, mol_type(start_mol:end_mol))
2645 : ELSE
2646 : CALL check_for_overlap(force_env, nchains(:, box_number), &
2647 10 : nunits(:), loverlap, mol_type(start_mol:end_mol))
2648 10 : IF (loverlap) lreject = .TRUE.
2649 : END IF
2650 :
2651 : ! check if cluster size changes then reject the move
2652 10 : IF (lbias) THEN
2653 0 : IF (total_clusafmo .NE. total_clus) THEN
2654 0 : loverlap = .TRUE.
2655 : END IF
2656 : ELSE
2657 10 : IF (total_clusafmo .NE. total_clus) THEN
2658 0 : loverlap = .TRUE.
2659 0 : lreject = .TRUE.
2660 : END IF
2661 : END IF
2662 :
2663 : ! if we're biasing with a cheaper potential, check for acceptance
2664 10 : IF (lbias) THEN
2665 :
2666 : ! here's where we bias the moves
2667 0 : IF (loverlap) THEN
2668 : w = 0.0E0_dp
2669 : ELSE
2670 0 : CALL force_env_calc_energy_force(bias_env, calc_force=.FALSE.)
2671 : CALL force_env_get(bias_env, &
2672 0 : potential_energy=bias_energy_new)
2673 : ! accept or reject the move based on the Metropolis rule
2674 0 : value = -BETA*(bias_energy_new - bias_energy_old)
2675 0 : IF (value .GT. exp_max_val) THEN
2676 : w = 10.0_dp
2677 0 : ELSEIF (value .LT. exp_min_val) THEN
2678 : w = 0.0_dp
2679 : ELSE
2680 0 : w = EXP(value)
2681 : END IF
2682 :
2683 : END IF
2684 :
2685 0 : IF (w .GE. 1.0E0_dp) THEN
2686 0 : w = 1.0E0_dp
2687 0 : rand = 0.0E0_dp
2688 : ELSE
2689 0 : IF (ionode) rand = rng_stream%next()
2690 0 : CALL group%bcast(rand, source)
2691 : END IF
2692 0 : IF (rand .LT. w) THEN
2693 :
2694 : ! accept the move
2695 0 : moves%bias_cltrans%successes = moves%bias_cltrans%successes + 1
2696 0 : move_updates%bias_cltrans%successes = move_updates%bias_cltrans%successes + 1
2697 0 : moves%cltrans%qsuccesses = moves%cltrans%qsuccesses + 1
2698 : move_updates%cltrans%successes = &
2699 0 : move_updates%cltrans%successes + 1
2700 0 : moves%qcltrans_dis = moves%qcltrans_dis + ABS(dis_mol)
2701 : bias_energy = bias_energy + bias_energy_new - &
2702 0 : bias_energy_old
2703 :
2704 : ELSE
2705 :
2706 : ! reject the move
2707 : ! restore the coordinates
2708 0 : CALL force_env_get(bias_env, subsys=subsys)
2709 0 : CALL cp_subsys_get(subsys, particles=particles)
2710 0 : DO ipart = 1, nunits_tot(box_number)
2711 0 : particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2712 : END DO
2713 0 : CALL cp_subsys_set(subsys, particles=particles)
2714 :
2715 : END IF
2716 :
2717 : END IF
2718 :
2719 : ! deallocate some stuff
2720 10 : DEALLOCATE (cluster)
2721 10 : DEALLOCATE (r_old)
2722 :
2723 : ! end the timing
2724 10 : CALL timestop(handle)
2725 :
2726 10 : END SUBROUTINE mc_cluster_translation
2727 :
2728 : END MODULE mc_moves
|
