(git:58e3e09)
mc_moves.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief 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
18  USE cell_methods, ONLY: cell_create
19  USE cell_types, ONLY: cell_clone,&
20  cell_release,&
21  cell_type,&
22  get_cell
25  cp_logger_type
26  USE cp_subsys_types, ONLY: cp_subsys_get,&
28  cp_subsys_type
30  USE force_env_types, ONLY: force_env_get,&
31  force_env_type,&
33  USE global_types, ONLY: global_environment_type
34  USE kinds, ONLY: default_string_length,&
35  dp
36  USE mathconstants, ONLY: pi
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,&
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 
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  RECURSIVE SUBROUTINE depth_first_search(current_atom, avoid_atom, &
93  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  DO iatom = 1, 6
102  IF (connectivity(iatom, current_atom) .NE. 0) THEN
103  IF (connectivity(iatom, current_atom) .NE. avoid_atom) THEN
104  atom(connectivity(iatom, current_atom)) = 1
105  CALL depth_first_search(connectivity(iatom, current_atom), &
106  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  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  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  INTEGER, DIMENSION(:), POINTER :: mol_type, nunits
158  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  CALL timeset(routinen, handle)
174 
175 ! nullify some pointers
176  NULLIFY (particles, subsys, molecule_kinds, molecule_kind, &
177  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  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  mol_type=mol_type, names=names)
185 
186 ! do some allocation
187  nunits_mol = nunits(molecule_type)
188  ALLOCATE (r_old(1:3, 1:nunits_mol))
189  ALLOCATE (r_new(1:3, 1:nunits_mol))
190 
191 ! find out some bounds for mol_type
192  start_mol = 1
193  DO jbox = 1, box_number - 1
194  start_mol = start_mol + sum(nchains(:, jbox))
195  END DO
196  end_mol = start_mol + sum(nchains(:, box_number)) - 1
197 
198 ! figure out which molecule number we are
199  end_atom = start_atom + nunits_mol - 1
200  molecule_number = 0
201  atom_number = 1
202  DO imolecule = 1, sum(nchains(:, box_number))
203  IF (atom_number == start_atom) THEN
204  molecule_number = imolecule
205  EXIT
206  END IF
207  atom_number = atom_number + nunits(mol_type(imolecule + start_mol - 1))
208  END DO
209  IF (molecule_number == 0) cpabort('Cannot find the molecule number')
210 
211 ! are we biasing this move?
212  IF (lbias) THEN
213 
214 ! grab the coordinates
215  CALL force_env_get(bias_env, subsys=subsys)
216 ! save the energy
217  bias_energy_old = bias_energy
218 
219  ELSE
220 
221 ! grab the coordinates
222  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  particles=particles, molecule_kinds=molecule_kinds)
228  DO imol_type = 1, SIZE(molecule_kinds%els(:))
229  molecule_kind_test => molecule_kinds%els(imol_type)
230  CALL get_molecule_kind(molecule_kind_test, name=name)
231  IF (trim(adjustl(name)) == trim(adjustl(names(molecule_type)))) THEN
232  molecule_kind => molecule_kinds%els(imol_type)
233  EXIT
234  END IF
235  END DO
236 
237 ! save the coordinates
238  DO ipart = start_atom, end_atom
239  r_old(1:3, ipart - start_atom + 1) = particles%els(ipart)%r(1:3)
240  END DO
241 
242  IF (.NOT. ASSOCIATED(molecule_kind)) cpabort('Cannot find the molecule type')
243 ! do the move
244  IF (move_type == 'bond') THEN
245 
246 ! record the attempt
247  moves%bond%attempts = moves%bond%attempts + 1
248  move_updates%bond%attempts = move_updates%bond%attempts + 1
249  moves%bias_bond%attempts = moves%bias_bond%attempts + 1
250  move_updates%bias_bond%attempts = move_updates%bias_bond%attempts + 1
251  IF (.NOT. lbias) THEN
252  moves%bond%qsuccesses = moves%bond%qsuccesses + 1
253  move_updates%bond%qsuccesses = &
254  move_updates%bond%qsuccesses + 1
255  moves%bias_bond%qsuccesses = moves%bias_bond%qsuccesses + 1
256  move_updates%bias_bond%qsuccesses = &
257  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  molecule_kind, dis_length, particles, rng_stream)
263 
264  ELSEIF (move_type == 'angle') THEN
265 
266 ! record the attempt
267  moves%angle%attempts = moves%angle%attempts + 1
268  move_updates%angle%attempts = move_updates%angle%attempts + 1
269  moves%bias_angle%attempts = moves%bias_angle%attempts + 1
270  move_updates%bias_angle%attempts = move_updates%bias_angle%attempts + 1
271  IF (.NOT. lbias) THEN
272  moves%angle%qsuccesses = moves%angle%qsuccesses + 1
273  move_updates%angle%qsuccesses = &
274  move_updates%angle%qsuccesses + 1
275  moves%bias_angle%qsuccesses = moves%bias_angle%qsuccesses + 1
276  move_updates%bias_angle%qsuccesses = &
277  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  molecule_kind, particles, rng_stream)
283  dis_length = 1.0e0_dp
284  ELSE
285 ! record the attempt
286  moves%dihedral%attempts = moves%dihedral%attempts + 1
287  move_updates%dihedral%attempts = move_updates%dihedral%attempts + 1
288  moves%bias_dihedral%attempts = moves%bias_dihedral%attempts + 1
289  move_updates%bias_dihedral%attempts = move_updates%bias_dihedral%attempts + 1
290  IF (.NOT. lbias) THEN
291  moves%dihedral%qsuccesses = moves%dihedral%qsuccesses + 1
292  move_updates%dihedral%qsuccesses = &
293  move_updates%dihedral%qsuccesses + 1
294  moves%bias_dihedral%qsuccesses = moves%bias_dihedral%qsuccesses + 1
295  move_updates%bias_dihedral%qsuccesses = &
296  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  molecule_kind, particles, rng_stream)
302  dis_length = 1.0e0_dp
303 
304  END IF
305 
306 ! set the coordinates
307  DO ipart = start_atom, end_atom
308  particles%els(ipart)%r(1:3) = r_new(1:3, ipart - start_atom + 1)
309  END DO
310 
311 ! check for overlap
312  lreject = .false.
313  IF (lbias) THEN
314  CALL check_for_overlap(bias_env, nchains(:, box_number), &
315  nunits(:), loverlap, mol_type(start_mol:end_mol), &
316  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  molecule_number=molecule_number)
321  IF (loverlap) lreject = .true.
322  END IF
323 
324 ! if we're biasing classical, check for acceptance
325  IF (lbias) THEN
326 
327 ! here's where we bias the moves
328 
329  IF (loverlap) THEN
330  w = 0.0e0_dp
331  ELSE
332  CALL force_env_calc_energy_force(bias_env, calc_force=.false.)
333  CALL force_env_get(bias_env, &
334  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  value = -beta*(bias_energy_new - bias_energy_old)
339  IF (value .GT. exp_max_val) THEN
340  w = 10.0_dp
341  ELSEIF (value .LT. exp_min_val) THEN
342  w = 0.0_dp
343  ELSE
344  w = exp(value)*dis_length**2
345  END IF
346 
347  END IF
348 
349  IF (w .GE. 1.0e0_dp) THEN
350  w = 1.0e0_dp
351  rand = 0.0e0_dp
352  ELSE
353  IF (ionode) THEN
354  rand = rng_stream%next()
355  END IF
356  CALL group%bcast(rand, source)
357  END IF
358 
359  IF (rand .LT. w) THEN
360 
361 ! accept the move
362  IF (move_type == 'bond') THEN
363  moves%bond%qsuccesses = moves%bond%qsuccesses + 1
364  move_updates%bond%successes = &
365  move_updates%bond%successes + 1
366  moves%bias_bond%successes = moves%bias_bond%successes + 1
367  move_updates%bias_bond%successes = &
368  move_updates%bias_bond%successes + 1
369  ELSEIF (move_type == 'angle') THEN
370  moves%angle%qsuccesses = moves%angle%qsuccesses + 1
371  move_updates%angle%successes = &
372  move_updates%angle%successes + 1
373  moves%bias_angle%successes = moves%bias_angle%successes + 1
374  move_updates%bias_angle%successes = &
375  move_updates%bias_angle%successes + 1
376  ELSE
377  moves%dihedral%qsuccesses = moves%dihedral%qsuccesses + 1
378  move_updates%dihedral%successes = &
379  move_updates%dihedral%successes + 1
380  moves%bias_dihedral%successes = moves%bias_dihedral%successes + 1
381  move_updates%bias_dihedral%successes = &
382  move_updates%bias_dihedral%successes + 1
383  END IF
384 
385  bias_energy = bias_energy + bias_energy_new - &
386  bias_energy_old
387 
388  ELSE
389 
390 ! reject the move
391 ! restore the coordinates
392  CALL force_env_get(bias_env, subsys=subsys)
393  CALL cp_subsys_get(subsys, particles=particles)
394  DO ipart = start_atom, end_atom
395  particles%els(ipart)%r(1:3) = r_old(1:3, ipart - start_atom + 1)
396  END DO
397  CALL cp_subsys_set(subsys, particles=particles)
398 
399  END IF
400 
401  END IF
402 
403 ! deallocate some stuff
404  DEALLOCATE (r_old)
405  DEALLOCATE (r_new)
406 
407 ! end the timing
408  CALL timestop(handle)
409 
410  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  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  INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
451  INTEGER, DIMENSION(:, :), POINTER :: nchains
452  LOGICAL :: ionode, lbias, loverlap
453  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  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  CALL timeset(routinen, handle)
467 
468 ! nullify some pointers
469  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  group=group, mc_molecule_info=mc_molecule_info)
476  CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot, &
477  nchains=nchains, nunits=nunits, mol_type=mol_type)
478 
479 ! find out some bounds for mol_type
480  start_mol = 1
481  DO jbox = 1, box_number - 1
482  start_mol = start_mol + sum(nchains(:, jbox))
483  END DO
484  end_mol = start_mol + sum(nchains(:, box_number)) - 1
485 
486 ! do some allocation
487  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  nunits_mol = nunits(molecule_type)
491  end_atom = start_atom + nunits_mol - 1
492  molecule_number = 0
493  atom_number = 1
494  DO imolecule = 1, sum(nchains(:, box_number))
495  IF (atom_number == start_atom) THEN
496  molecule_number = imolecule
497  EXIT
498  END IF
499  atom_number = atom_number + nunits(mol_type(imolecule + start_mol - 1))
500  END DO
501  IF (molecule_number == 0) cpabort('Cannot find the molecule number')
502 
503 ! are we biasing this move?
504  IF (lbias) THEN
505 
506 ! grab the coordinates
507  CALL force_env_get(bias_env, subsys=subsys)
508  CALL cp_subsys_get(subsys, particles=particles)
509 
510 ! save the coordinates
511  DO ipart = 1, nunits_tot(box_number)
512  r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
513  END DO
514 
515 ! save the energy
516  bias_energy_old = bias_energy
517 
518  ELSE
519 
520 ! grab the coordinates
521  CALL force_env_get(force_env, subsys=subsys)
522  CALL cp_subsys_get(subsys, particles=particles)
523  END IF
524 
525 ! record the attempt
526  moves%trans%attempts = moves%trans%attempts + 1
527  move_updates%trans%attempts = move_updates%trans%attempts + 1
528  moves%bias_trans%attempts = moves%bias_trans%attempts + 1
529  move_updates%bias_trans%attempts = move_updates%bias_trans%attempts + 1
530  IF (.NOT. lbias) THEN
531  moves%trans%qsuccesses = moves%trans%qsuccesses + 1
532  move_updates%trans%qsuccesses = move_updates%trans%qsuccesses + 1
533  moves%bias_trans%qsuccesses = moves%bias_trans%qsuccesses + 1
534  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  IF (ionode) rand = rng_stream%next()
541  CALL group%bcast(rand, source)
542  ! 1,2,3 with equal prob
543  move_direction = int(3*rand) + 1
544 
545 ! call a random number to figure out how far we're moving
546  IF (ionode) rand = rng_stream%next()
547  CALL group%bcast(rand, source)
548  dis_mol = rmtrans(molecule_type)*(rand - 0.5e0_dp)*2.0e0_dp
549 
550 ! do the move
551  DO iparticle = start_atom, end_atom
552  particles%els(iparticle)%r(move_direction) = &
553  particles%els(iparticle)%r(move_direction) + dis_mol
554  END DO
555  CALL cp_subsys_set(subsys, particles=particles)
556 
557 ! figure out if there is any overlap...need the number of the molecule
558  lreject = .false.
559  IF (lbias) THEN
560  CALL check_for_overlap(bias_env, nchains(:, box_number), &
561  nunits(:), loverlap, mol_type(start_mol:end_mol), &
562  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  molecule_number=molecule_number)
567  IF (loverlap) lreject = .true.
568  END IF
569 
570 ! if we're biasing with a cheaper potential, check for acceptance
571  IF (lbias) THEN
572 
573 ! here's where we bias the moves
574  IF (loverlap) THEN
575  w = 0.0e0_dp
576  ELSE
577  CALL force_env_calc_energy_force(bias_env, calc_force=.false.)
578  CALL force_env_get(bias_env, &
579  potential_energy=bias_energy_new)
580 ! accept or reject the move based on the Metropolis rule
581  value = -beta*(bias_energy_new - bias_energy_old)
582  IF (value .GT. exp_max_val) THEN
583  w = 10.0_dp
584  ELSEIF (value .LT. exp_min_val) THEN
585  w = 0.0_dp
586  ELSE
587  w = exp(value)
588  END IF
589 
590  END IF
591 
592  IF (w .GE. 1.0e0_dp) THEN
593  w = 1.0e0_dp
594  rand = 0.0e0_dp
595  ELSE
596  IF (ionode) rand = rng_stream%next()
597  CALL group%bcast(rand, source)
598  END IF
599 
600  IF (rand .LT. w) THEN
601 
602 ! accept the move
603  moves%bias_trans%successes = moves%bias_trans%successes + 1
604  move_updates%bias_trans%successes = move_updates%bias_trans%successes + 1
605  moves%trans%qsuccesses = moves%trans%qsuccesses + 1
606  move_updates%trans%successes = &
607  move_updates%trans%successes + 1
608  moves%qtrans_dis = moves%qtrans_dis + abs(dis_mol)
609  bias_energy = bias_energy + bias_energy_new - &
610  bias_energy_old
611 
612  ELSE
613 
614 ! reject the move
615 ! restore the coordinates
616  CALL force_env_get(bias_env, subsys=subsys)
617  CALL cp_subsys_get(subsys, particles=particles)
618  DO ipart = 1, nunits_tot(box_number)
619  particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
620  END DO
621  CALL cp_subsys_set(subsys, particles=particles)
622 
623  END IF
624 
625  END IF
626 
627 ! deallocate some stuff
628  DEALLOCATE (r_old)
629 
630 ! end the timing
631  CALL timestop(handle)
632 
633  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  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  INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
674  INTEGER, DIMENSION(:, :), POINTER :: nchains
675  LOGICAL :: ionode, lbias, loverlap, lx, ly
676  REAL(dp), DIMENSION(:), POINTER :: rmrot
677  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  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  CALL timeset(routinen, handle)
690 
691  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  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  mol_type=mol_type)
701 
702 ! figure out some bounds for mol_type
703  start_mol = 1
704  DO jbox = 1, box_number - 1
705  start_mol = start_mol + sum(nchains(:, jbox))
706  END DO
707  end_mol = start_mol + sum(nchains(:, box_number)) - 1
708 
709  nunits_mol = nunits(molecule_type)
710 
711 ! nullify some pointers
712  NULLIFY (particles, subsys)
713 
714 ! do some allocation
715  ALLOCATE (r_old(1:3, 1:nunits_tot(box_number)))
716 
717 ! initialize some stuff
718  lx = .false.
719  ly = .false.
720 
721 ! determine what the final atom in the molecule is numbered, and which
722 ! molecule number this is
723  end_atom = start_atom + nunits_mol - 1
724  molecule_number = 0
725  atom_number = 1
726  DO imolecule = 1, sum(nchains(:, box_number))
727  IF (atom_number == start_atom) THEN
728  molecule_number = imolecule
729  EXIT
730  END IF
731  atom_number = atom_number + nunits(mol_type(imolecule + start_mol - 1))
732  END DO
733  IF (molecule_number == 0) cpabort('Cannot find the molecule number')
734 
735 ! are we biasing this move?
736  IF (lbias) THEN
737 
738 ! grab the coordinates
739  CALL force_env_get(bias_env, subsys=subsys)
740  CALL cp_subsys_get(subsys, particles=particles)
741 
742 ! save the coordinates
743  DO ipart = 1, nunits_tot(box_number)
744  r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
745  END DO
746 
747 ! save the energy
748  bias_energy_old = bias_energy
749 
750  ELSE
751 
752 ! grab the coordinates
753  CALL force_env_get(force_env, subsys=subsys)
754  CALL cp_subsys_get(subsys, particles=particles)
755  END IF
756 
757 ! grab the masses
758  masstot = sum(mass(1:nunits(molecule_type), molecule_type))
759 
760 ! record the attempt
761  moves%bias_rot%attempts = moves%bias_rot%attempts + 1
762  move_updates%bias_rot%attempts = move_updates%bias_rot%attempts + 1
763  moves%rot%attempts = moves%rot%attempts + 1
764  move_updates%rot%attempts = move_updates%rot%attempts + 1
765  IF (.NOT. lbias) THEN
766  moves%rot%qsuccesses = moves%rot%qsuccesses + 1
767  move_updates%rot%qsuccesses = move_updates%rot%qsuccesses + 1
768  moves%bias_rot%qsuccesses = moves%bias_rot%qsuccesses + 1
769  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  IF (ionode) rand = rng_stream%next()
776 ! CALL RANDOM_NUMBER(rand)
777  CALL group%bcast(rand, source)
778  ! 1,2,3 with equal prob
779  dir = int(3*rand) + 1
780 
781  IF (dir .EQ. 1) THEN
782  lx = .true.
783  ELSEIF (dir .EQ. 2) THEN
784  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  nxcm = 0.0e0_dp
790  nycm = 0.0e0_dp
791  nzcm = 0.0e0_dp
792  DO ii = 1, nunits_mol
793  nxcm = nxcm + particles%els(start_atom - 1 + ii)%r(1)*mass(ii, molecule_type)
794  nycm = nycm + particles%els(start_atom - 1 + ii)%r(2)*mass(ii, molecule_type)
795  nzcm = nzcm + particles%els(start_atom - 1 + ii)%r(3)*mass(ii, molecule_type)
796  END DO
797  nxcm = nxcm/masstot
798  nycm = nycm/masstot
799  nzcm = nzcm/masstot
800 
801 ! call a random number to figure out how far we're moving
802  IF (ionode) rand = rng_stream%next()
803  CALL group%bcast(rand, source)
804  dgamma = rmrot(molecule_type)*(rand - 0.5e0_dp)*2.0e0_dp
805 
806 ! *** set up the rotation matrix ***
807 
808  cosdg = cos(dgamma)
809  sindg = sin(dgamma)
810 
811  IF (lx) THEN
812 
813 ! *** ROTATE UNITS OF I AROUND X-AXIS ***
814 
815  DO iunit = start_atom, end_atom
816  ry = particles%els(iunit)%r(2) - nycm
817  rz = particles%els(iunit)%r(3) - nzcm
818  rynew = cosdg*ry - sindg*rz
819  rznew = cosdg*rz + sindg*ry
820 
821  particles%els(iunit)%r(2) = rynew + nycm
822  particles%els(iunit)%r(3) = rznew + nzcm
823 
824  END DO
825  ELSEIF (ly) THEN
826 
827 ! *** ROTATE UNITS OF I AROUND y-AXIS ***
828 
829  DO iunit = start_atom, end_atom
830  rx = particles%els(iunit)%r(1) - nxcm
831  rz = particles%els(iunit)%r(3) - nzcm
832  rxnew = cosdg*rx + sindg*rz
833  rznew = cosdg*rz - sindg*rx
834 
835  particles%els(iunit)%r(1) = rxnew + nxcm
836  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  DO iunit = start_atom, end_atom
845  rx = particles%els(iunit)%r(1) - nxcm
846  ry = particles%els(iunit)%r(2) - nycm
847 
848  rxnew = cosdg*rx - sindg*ry
849  rynew = cosdg*ry + sindg*rx
850 
851  particles%els(iunit)%r(1) = rxnew + nxcm
852  particles%els(iunit)%r(2) = rynew + nycm
853 
854  END DO
855 
856  END IF
857  CALL cp_subsys_set(subsys, particles=particles)
858 
859 ! check for overlap
860  lreject = .false.
861  IF (lbias) THEN
862  CALL check_for_overlap(bias_env, nchains(:, box_number), &
863  nunits(:), loverlap, mol_type(start_mol:end_mol), &
864  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  molecule_number=molecule_number)
869  IF (loverlap) lreject = .true.
870  END IF
871 
872 ! if we're biasing classical, check for acceptance
873  IF (lbias) THEN
874 
875 ! here's where we bias the moves
876 
877  IF (loverlap) THEN
878  w = 0.0e0_dp
879  ELSE
880  CALL force_env_calc_energy_force(bias_env, calc_force=.false.)
881  CALL force_env_get(bias_env, &
882  potential_energy=bias_energy_new)
883 ! accept or reject the move based on the Metropolis rule
884  value = -beta*(bias_energy_new - bias_energy_old)
885  IF (value .GT. exp_max_val) THEN
886  w = 10.0_dp
887  ELSEIF (value .LT. exp_min_val) THEN
888  w = 0.0_dp
889  ELSE
890  w = exp(value)
891  END IF
892 
893  END IF
894 
895  IF (w .GE. 1.0e0_dp) THEN
896  w = 1.0e0_dp
897  rand = 0.0e0_dp
898  ELSE
899  IF (ionode) rand = rng_stream%next()
900  CALL group%bcast(rand, source)
901  END IF
902 
903  IF (rand .LT. w) THEN
904 
905 ! accept the move
906  moves%bias_rot%successes = moves%bias_rot%successes + 1
907  move_updates%bias_rot%successes = move_updates%bias_rot%successes + 1
908  moves%rot%qsuccesses = moves%rot%qsuccesses + 1
909  move_updates%rot%successes = move_updates%rot%successes + 1
910  bias_energy = bias_energy + bias_energy_new - &
911  bias_energy_old
912 
913  ELSE
914 
915 ! reject the move
916 ! restore the coordinates
917  CALL force_env_get(bias_env, subsys=subsys)
918  CALL cp_subsys_get(subsys, particles=particles)
919  DO ipart = 1, nunits_tot(box_number)
920  particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
921  END DO
922  CALL cp_subsys_set(subsys, particles=particles)
923 
924  END IF
925 
926  END IF
927 
928 ! deallocate some stuff
929  DEALLOCATE (r_old)
930 
931 ! end the timing
932  CALL timestop(handle)
933 
934  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  SUBROUTINE mc_volume_move(mc_par, force_env, moves, move_updates, &
960  old_energy, box_number, &
961  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  INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
981  INTEGER, DIMENSION(:, :), POINTER :: nchains
982  LOGICAL :: ionode, ldiscrete, lincrease, loverlap, &
983  ltoo_small
984  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  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  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  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  mass=mass)
1012 ! figure out some bounds for mol_type
1013  start_mol = 1
1014  DO jbox = 1, box_number - 1
1015  start_mol = start_mol + sum(nchains(:, jbox))
1016  END DO
1017  end_mol = start_mol + sum(nchains(:, box_number)) - 1
1018 
1019  print_level = 1 ! hack, printlevel is for print_keys
1020 
1021 ! nullify some pointers
1022  NULLIFY (particles_old, cell_old, oldsys, cell_test, cell)
1023 
1024 ! do some allocation
1025  ALLOCATE (r(1:3, 1:nunits_tot(box_number)))
1026 
1027 ! record the attempt
1028  moves%volume%attempts = moves%volume%attempts + 1
1029  move_updates%volume%attempts = move_updates%volume%attempts + 1
1030 
1031 ! now let's grab the cell length and particle positions
1032  CALL force_env_get(force_env, subsys=oldsys, cell=cell)
1033  CALL get_cell(cell, abc=abc)
1034  CALL cell_create(cell_old)
1035  CALL cell_clone(cell, cell_old, tag="CELL_OLD")
1036  CALL cp_subsys_get(oldsys, particles=particles_old)
1037 
1038 ! find the old cell length
1039  old_cell_length(1) = abc(1)
1040  old_cell_length(2) = abc(2)
1041  old_cell_length(3) = abc(3)
1042 
1043 ! save the old coordinates
1044  DO iatom = 1, nunits_tot(box_number)
1045  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  IF (ionode) rand = rng_stream%next()
1052  CALL group%bcast(rand, source)
1053 
1054 ! find the test cell lengths for the discrete volume move
1055  IF (ldiscrete) THEN
1056  IF (rand .LT. 0.5_dp) THEN
1057  lincrease = .true.
1058  ELSE
1059  lincrease = .false.
1060  END IF
1061 
1062  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  IF (lincrease) THEN
1066  DO
1067  IF (ionode) rand = rng_stream%next()
1068  CALL group%bcast(rand, source)
1069  iside_change = ceiling(3.0_dp*rand)
1070  IF (discrete_array(iside_change, 1) .EQ. 1) THEN
1071  new_cell_length(iside_change) = &
1072  new_cell_length(iside_change) + discrete_step
1073  EXIT
1074  END IF
1075  END DO
1076  ELSE
1077  DO
1078  IF (ionode) rand = rng_stream%next()
1079  CALL group%bcast(rand, source)
1080  iside_change = ceiling(3.0_dp*rand)
1081  IF (discrete_array(iside_change, 2) .EQ. 1) THEN
1082  new_cell_length(iside_change) = &
1083  new_cell_length(iside_change) - discrete_step
1084  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  - 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  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  old_cell_length(3)
1105 
1106  IF (temp_var .LE. 0.0e0_dp) THEN
1107  loverlap = .true. ! cannot have a negative volume
1108  ELSE
1109  new_cell_length(1) = (temp_var)**(1.0e0_dp/3.0e0_dp)
1110  new_cell_length(2) = new_cell_length(1)
1111  new_cell_length(3) = new_cell_length(1)
1112  loverlap = .false.
1113  END IF
1114  END IF
1115  CALL group%bcast(loverlap, source)
1116 
1117  IF (loverlap) THEN
1118 ! deallocate some stuff
1119  DEALLOCATE (r)
1120  logger => cp_get_default_logger()
1121  output_unit = cp_logger_get_default_io_unit(logger)
1122  IF (output_unit > 0) WRITE (output_unit, *) &
1123  "Volume move rejected because we tried to make too small of box.", vol_dis
1124 ! end the timing
1125  CALL timestop(handle)
1126  RETURN
1127  END IF
1128 
1129 ! now we need to make the new cell
1130  hmat_test(:, :) = 0.0e0_dp
1131  hmat_test(1, 1) = new_cell_length(1)
1132  hmat_test(2, 2) = new_cell_length(2)
1133  hmat_test(3, 3) = new_cell_length(3)
1134  CALL cell_create(cell_test, hmat=hmat_test(:, :), periodic=cell%perd)
1135  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  end_atom = 0
1144  DO imolecule = 1, sum(nchains(:, box_number))
1145  nunits_mol = nunits(mol_type(imolecule + start_mol - 1))
1146  start_atom = end_atom + 1
1147  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  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  DO iside = 1, 3
1155  center_of_mass_new(iside) = center_of_mass(iside)* &
1156  new_cell_length(iside)/old_cell_length(iside)
1157  END DO
1158 
1159  DO idim = 1, 3
1160  diff(idim) = center_of_mass_new(idim) - center_of_mass(idim)
1161 ! now change the particle positions
1162  DO iunit = start_atom, end_atom
1163  particles_old%els(iunit)%r(idim) = &
1164  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  cell_length=new_cell_length)
1173 
1174 ! figure out if we have overlap problems
1175  CALL group%bcast(loverlap, source)
1176  IF (loverlap) THEN
1177 ! deallocate some stuff
1178  DEALLOCATE (r)
1179 
1180  logger => cp_get_default_logger()
1181  output_unit = cp_logger_get_default_io_unit(logger)
1182  IF (output_unit > 0) WRITE (output_unit, *) &
1183  "Volume move rejected due to overlap.", vol_dis
1184 ! end the timing
1185  CALL timestop(handle)
1186 ! reset the cell and particle positions
1187  CALL cp_subsys_set(oldsys, cell=cell_old)
1188  DO iatom = 1, nunits_tot(box_number)
1189  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  IF (ionode) THEN
1196  ltoo_small = .false.
1197  IF (force_env%in_use .EQ. use_fist_force) THEN
1198  CALL get_mc_par(mc_par, rcut=rcut)
1199  IF (new_cell_length(1) .LT. 2.0_dp*rcut) ltoo_small = .true.
1200  IF (new_cell_length(2) .LT. 2.0_dp*rcut) ltoo_small = .true.
1201  IF (new_cell_length(3) .LT. 2.0_dp*rcut) ltoo_small = .true.
1202 
1203  IF (ltoo_small) THEN
1204  WRITE (iw, *) 'new_cell_lengths ', &
1205  new_cell_length(1:3)/angstrom
1206  WRITE (iw, *) 'rcut ', rcut/angstrom
1207  END IF
1208  END IF
1209  END IF
1210  CALL group%bcast(ltoo_small, source)
1211  IF (ltoo_small) &
1212  cpabort("Attempted a volume move where box size got too small.")
1213 
1214 ! now compute the energy
1215  CALL force_env_calc_energy_force(force_env, calc_force=.false.)
1216  CALL force_env_get(force_env, &
1217  potential_energy=new_energy)
1218 
1219 ! accept or reject the move
1220 ! to prevent overflows
1221  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  (old_cell_length(1)*old_cell_length(2)*old_cell_length(3)))
1225  pressure_term = pressure*vol_dis
1226 
1227  value = -beta*(energy_term + volume_term + pressure_term)
1228  IF (value .GT. exp_max_val) THEN
1229  w = 10.0_dp
1230  ELSEIF (value .LT. exp_min_val) THEN
1231  w = 0.0_dp
1232  ELSE
1233  w = exp(value)
1234  END IF
1235 
1236 !!!!!!!!!!!!!!!! for E_V curves
1237 ! w=1.0E0_dp
1238 ! w=0.0E0_dp
1239 
1240  IF (w .GE. 1.0e0_dp) THEN
1241  w = 1.0e0_dp
1242  rand = 0.0e0_dp
1243  ELSE
1244  IF (ionode) rand = rng_stream%next()
1245  CALL group%bcast(rand, source)
1246  END IF
1247 
1248  IF (rand .LT. w) THEN
1249 
1250 ! accept the move
1251  moves%volume%successes = moves%volume%successes + 1
1252  move_updates%volume%successes = move_updates%volume%successes + 1
1253 
1254 ! update energies
1255  energy_check = energy_check + (new_energy - old_energy)
1256  old_energy = new_energy
1257 
1258  DO iatom = 1, nunits_tot(box_number)
1259  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  IF (ldiscrete) THEN
1264  CALL create_discrete_array(new_cell_length(:), &
1265  discrete_array(:, :), discrete_step)
1266  END IF
1267 
1268  ELSE
1269 
1270 ! reset the cell and particle positions
1271  CALL cp_subsys_set(oldsys, cell=cell_old)
1272  DO iatom = 1, nunits_tot(box_number)
1273  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  DEALLOCATE (r)
1280  CALL cell_release(cell_test)
1281  CALL cell_release(cell_old)
1282 
1283 ! end the timing
1284  CALL timestop(handle)
1285 
1286  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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_b, counter
1322  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity
1323  INTEGER, DIMENSION(:), POINTER :: nunits
1324  LOGICAL :: ionode
1325  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  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  CALL timeset(routinen, handle)
1336 
1337  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  group=group, rmbond=rmbond, ionode=ionode)
1342  CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
1343 
1344 ! copy the incoming coordinates so we can change them
1345  DO ipart = 1, nunits(molecule_type)
1346  r_new(1:3, ipart) = r_old(1:3, ipart)
1347  END DO
1348 
1349 ! pick which bond in the molecule at random
1350  IF (ionode) THEN
1351  rand = rng_stream%next()
1352  END IF
1353  CALL group%bcast(rand, source)
1354  CALL get_molecule_kind(molecule_kind, natom=natom, nbond=nbond, &
1355  bond_list=bond_list)
1356  bond_number = ceiling(rand*real(nbond, dp))
1357 
1358  ALLOCATE (connection(1:natom, 1:2))
1359 ! assume at most six bonds per atom
1360  ALLOCATE (connectivity(1:6, 1:natom))
1361  ALLOCATE (counter(1:natom))
1362  ALLOCATE (atom_a(1:natom))
1363  ALLOCATE (atom_b(1:natom))
1364  connection(:, :) = 0
1365  connectivity(:, :) = 0
1366  counter(:) = 0
1367  atom_a(:) = 0
1368  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  DO iatom = 1, natom
1373  DO ibond = 1, nbond
1374  IF (bond_list(ibond)%a == iatom) THEN
1375  counter(iatom) = counter(iatom) + 1
1376  connectivity(counter(iatom), iatom) = bond_list(ibond)%b
1377  ELSEIF (bond_list(ibond)%b == iatom) THEN
1378  counter(iatom) = counter(iatom) + 1
1379  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  atom_a(:) = 0
1387  atom_a(bond_list(bond_number)%a) = 1
1388  CALL depth_first_search(bond_list(bond_number)%a, bond_list(bond_number)%b, &
1389  connectivity(:, :), atom_a(:))
1390  atom_b(:) = 0
1391  atom_b(bond_list(bond_number)%b) = 1
1392  CALL depth_first_search(bond_list(bond_number)%b, bond_list(bond_number)%a, &
1393  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  mass_a = 0.0_dp
1398  mass_b = 0.0_dp
1399  DO iatom = 1, natom
1400  CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
1401  mass=atom_mass)
1402  IF (atom_a(iatom) == 1) THEN
1403  mass_a = mass_a + atom_mass
1404  ELSE
1405  mass_b = mass_b + atom_mass
1406  END IF
1407  END DO
1408 
1409 ! choose a displacement
1410  IF (ionode) rand = rng_stream%next()
1411  CALL group%bcast(rand, source)
1412 
1413  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  DO i = 1, 3
1417  bond_a(i) = r_new(i, bond_list(bond_number)%a) - &
1418  r_new(i, bond_list(bond_number)%b)
1419  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  old_length = sqrt(dot_product(bond_a, bond_a))
1425  new_length_a = dis_length*mass_b/(mass_a + mass_b)
1426  new_length_b = dis_length*mass_a/(mass_a + mass_b)
1427 
1428  DO i = 1, 3
1429  bond_a(i) = bond_a(i)/old_length*new_length_a
1430  bond_b(i) = bond_b(i)/old_length*new_length_b
1431  END DO
1432 
1433  DO iatom = 1, natom
1434  IF (atom_a(iatom) == 1) THEN
1435  r_new(1, iatom) = r_new(1, iatom) + bond_a(1)
1436  r_new(2, iatom) = r_new(2, iatom) + bond_a(2)
1437  r_new(3, iatom) = r_new(3, iatom) + bond_a(3)
1438  ELSE
1439  r_new(1, iatom) = r_new(1, iatom) + bond_b(1)
1440  r_new(2, iatom) = r_new(2, iatom) + bond_b(2)
1441  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  dis_length = (old_length + dis_length)/old_length
1447 
1448  DEALLOCATE (connection)
1449  DEALLOCATE (connectivity)
1450  DEALLOCATE (counter)
1451  DEALLOCATE (atom_a)
1452  DEALLOCATE (atom_b)
1453 ! end the timing
1454  CALL timestop(handle)
1455 
1456  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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_c, counter
1489  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity
1490  INTEGER, DIMENSION(:), POINTER :: nunits
1491  LOGICAL :: ionode
1492  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  TYPE(bend_type), DIMENSION(:), POINTER :: bend_list
1498  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  CALL timeset(routinen, handle)
1505 
1506  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  group=group, ionode=ionode, mc_molecule_info=mc_molecule_info)
1511  CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
1512 
1513 ! copy the incoming coordinates so we can change them
1514  DO ipart = 1, nunits(molecule_type)
1515  r_new(1:3, ipart) = r_old(1:3, ipart)
1516  END DO
1517 
1518 ! pick which bond in the molecule at random
1519  IF (ionode) THEN
1520  rand = rng_stream%next()
1521  END IF
1522  CALL group%bcast(rand, source)
1523  CALL get_molecule_kind(molecule_kind, natom=natom, nbend=nbend, &
1524  bend_list=bend_list, bond_list=bond_list, nbond=nbond)
1525  bend_number = ceiling(rand*real(nbend, dp))
1526 
1527  ALLOCATE (connection(1:natom, 1:2))
1528 ! assume at most six bonds per atom
1529  ALLOCATE (connectivity(1:6, 1:natom))
1530  ALLOCATE (counter(1:natom))
1531  ALLOCATE (atom_a(1:natom))
1532  ALLOCATE (atom_c(1:natom))
1533  connection(:, :) = 0
1534  connectivity(:, :) = 0
1535  counter(:) = 0
1536  atom_a(:) = 0
1537  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  DO iatom = 1, natom
1542  DO ibond = 1, nbond
1543  IF (bond_list(ibond)%a == iatom) THEN
1544  counter(iatom) = counter(iatom) + 1
1545  connectivity(counter(iatom), iatom) = bond_list(ibond)%b
1546  ELSEIF (bond_list(ibond)%b == iatom) THEN
1547  counter(iatom) = counter(iatom) + 1
1548  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  atom_a(:) = 0
1556  atom_a(bend_list(bend_number)%a) = 1
1557  CALL depth_first_search(bend_list(bend_number)%a, bend_list(bend_number)%b, &
1558  connectivity(:, :), atom_a(:))
1559  atom_c(:) = 0
1560  atom_c(bend_list(bend_number)%c) = 1
1561  CALL depth_first_search(bend_list(bend_number)%c, bend_list(bend_number)%b, &
1562  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  mass_a = 0.0_dp
1567  mass_c = 0.0_dp
1568  DO iatom = 1, natom
1569  CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
1570  mass=atom_mass)
1571  IF (atom_a(iatom) == 1) mass_a = mass_a + atom_mass
1572  IF (atom_c(iatom) == 1) mass_c = mass_c + atom_mass
1573  END DO
1574 
1575 ! choose a displacement
1576  IF (ionode) rand = rng_stream%next()
1577  CALL group%bcast(rand, source)
1578 
1579  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  DO i = 1, 3
1589  bond_a(i) = r_new(i, bend_list(bend_number)%a) - &
1590  r_new(i, bend_list(bend_number)%b)
1591  bond_c(i) = r_new(i, bend_list(bend_number)%c) - &
1592  r_new(i, bend_list(bend_number)%b)
1593  END DO
1594  old_length_a = sqrt(dot_product(bond_a, bond_a))
1595  old_length_c = sqrt(dot_product(bond_c, bond_c))
1596  old_angle = acos(dot_product(bond_a, bond_c)/(old_length_a*old_length_c))
1597 
1598  DO i = 1, 3
1599  bisector(i) = bond_a(i)/old_length_a + & ! not yet normalized
1600  bond_c(i)/old_length_c
1601  END DO
1602  bis_length = sqrt(dot_product(bisector, bisector))
1603  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  cross_prod(1) = bond_a(2)*bond_c(3) - bond_a(3)*bond_c(2)
1608  cross_prod(2) = bond_a(3)*bond_c(1) - bond_a(1)*bond_c(3)
1609  cross_prod(3) = bond_a(1)*bond_c(2) - bond_a(2)*bond_c(1)
1610  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  cross_prod_plane(1) = cross_prod(2)*bisector(3) - cross_prod(3)*bisector(2)
1614  cross_prod_plane(2) = cross_prod(3)*bisector(1) - cross_prod(1)*bisector(3)
1615  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  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  DO iatom = 1, natom
1623  r_new(1:3, iatom) = r_new(1:3, iatom) - &
1624  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  dis_angle_a = dis_angle*mass_c/(mass_a + mass_c)
1631  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  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  cross_prod(1:3)
1639  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  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  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  (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  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  (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  cross_prod(1:3)
1668  END IF
1669 
1670  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  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  (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  cross_prod(1:3)
1684  ELSE
1685  new_angle_c = acos(dot_product(bisector(1:3), temp(1:3))/ &
1686  (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  cross_prod(1:3)
1692  END IF
1693  END IF
1694 
1695  END DO
1696 
1697  DO iatom = 1, natom
1698  r_new(1:3, iatom) = r_new(1:3, iatom) + &
1699  r_old(1:3, bend_list(bend_number)%b)
1700  END DO
1701 
1702 ! deallocate some stuff
1703  DEALLOCATE (connection)
1704  DEALLOCATE (connectivity)
1705  DEALLOCATE (counter)
1706  DEALLOCATE (atom_a)
1707  DEALLOCATE (atom_c)
1708 
1709 ! end the timing
1710  CALL timestop(handle)
1711 
1712  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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_d, counter
1748  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity
1749  INTEGER, DIMENSION(:), POINTER :: nunits
1750  LOGICAL :: ionode
1751  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  TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
1757  TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1758  TYPE(mp_comm_type) :: group
1759  TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list
1760 
1761 ! begin the timing of the subroutine
1762 
1763  CALL timeset(routinen, handle)
1764 
1765  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  mc_molecule_info=mc_molecule_info)
1771  CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
1772 
1773 ! copy the incoming coordinates so we can change them
1774  DO ipart = 1, nunits(molecule_type)
1775  r_new(1:3, ipart) = r_old(1:3, ipart)
1776  END DO
1777 
1778 ! pick which bond in the molecule at random
1779  IF (ionode) THEN
1780  rand = rng_stream%next()
1781 ! CALL RANDOM_NUMBER(rand)
1782  END IF
1783  CALL group%bcast(rand, source)
1784  CALL get_molecule_kind(molecule_kind, natom=natom, &
1785  bond_list=bond_list, nbond=nbond, &
1786  ntorsion=ntorsion, torsion_list=torsion_list)
1787  torsion_number = ceiling(rand*real(ntorsion, dp))
1788 
1789  ALLOCATE (connection(1:natom, 1:2))
1790 ! assume at most six bonds per atom
1791  ALLOCATE (connectivity(1:6, 1:natom))
1792  ALLOCATE (counter(1:natom))
1793  ALLOCATE (atom_a(1:natom))
1794  ALLOCATE (atom_d(1:natom))
1795  connection(:, :) = 0
1796  connectivity(:, :) = 0
1797  counter(:) = 0
1798  atom_a(:) = 0
1799  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  DO iatom = 1, natom
1804  DO ibond = 1, nbond
1805  IF (bond_list(ibond)%a == iatom) THEN
1806  counter(iatom) = counter(iatom) + 1
1807  connectivity(counter(iatom), iatom) = bond_list(ibond)%b
1808  ELSEIF (bond_list(ibond)%b == iatom) THEN
1809  counter(iatom) = counter(iatom) + 1
1810  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  atom_a(:) = 0
1819  atom_a(torsion_list(torsion_number)%a) = 1
1820  CALL depth_first_search(torsion_list(torsion_number)%b, &
1821  torsion_list(torsion_number)%c, connectivity(:, :), atom_a(:))
1822  atom_d(:) = 0
1823  atom_d(torsion_list(torsion_number)%d) = 1
1824  CALL depth_first_search(torsion_list(torsion_number)%c, &
1825  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  mass_a = 0.0_dp
1830  mass_d = 0.0_dp
1831  DO iatom = 1, natom
1832  CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
1833  mass=atom_mass)
1834  IF (atom_a(iatom) == 1) mass_a = mass_a + atom_mass
1835  IF (atom_d(iatom) == 1) mass_d = mass_d + atom_mass
1836  END DO
1837 
1838 ! choose a displacement
1839  IF (ionode) rand = rng_stream%next()
1840  CALL group%bcast(rand, source)
1841 
1842  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  DO i = 1, 3
1846  bond_a(i) = r_new(i, torsion_list(torsion_number)%c) - &
1847  r_new(i, torsion_list(torsion_number)%b)
1848  END DO
1849  old_length_a = sqrt(dot_product(bond_a, bond_a))
1850  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  dis_angle_a = dis_angle*mass_d/(mass_a + mass_d)
1856  dis_angle_d = -dis_angle*mass_a/(mass_a + mass_d)
1857 
1858  DO iatom = 1, natom
1859 
1860  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  r_new(1:3, torsion_list(torsion_number)%b)
1864 
1865 ! multiply by the rotation matrix
1866  u = bond_a(1)
1867  v = bond_a(2)
1868  w = bond_a(3)
1869  x = r_new(1, iatom)
1870  y = r_new(2, iatom)
1871  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  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  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  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  temp(1:3) = temp(1:3) + r_new(1:3, torsion_list(torsion_number)%b)
1881  r_new(1:3, iatom) = temp(1:3)
1882 
1883  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  r_new(1:3, torsion_list(torsion_number)%c)
1888 
1889 ! multiply by the rotation matrix
1890  u = bond_a(1)
1891  v = bond_a(2)
1892  w = bond_a(3)
1893  x = r_new(1, iatom)
1894  y = r_new(2, iatom)
1895  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  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  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  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  temp(1:3) = temp(1:3) + r_new(1:3, torsion_list(torsion_number)%c)
1905  r_new(1:3, iatom) = temp(1:3)
1906  END IF
1907  END DO
1908 
1909 ! deallocate some stuff
1910  DEALLOCATE (connection)
1911  DEALLOCATE (connectivity)
1912  DEALLOCATE (counter)
1913  DEALLOCATE (atom_a)
1914  DEALLOCATE (atom_d)
1915 
1916 ! end the timing
1917  CALL timestop(handle)
1918 
1919  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  SUBROUTINE mc_avbmc_move(mc_par, force_env, bias_env, moves, &
1946  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  INTEGER, DIMENSION(:), POINTER :: avbmc_atom, mol_type, nunits, nunits_tot
1968  INTEGER, DIMENSION(:, :), POINTER :: nchains
1969  LOGICAL :: ionode, lbias, ldum, lin, loverlap
1970  REAL(dp), DIMENSION(:), POINTER :: avbmc_rmax, avbmc_rmin, pbias
1971  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  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  rdum = 1.0_dp
1986 
1987 ! begin the timing of the subroutine
1988  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  group=group, pbias=pbias, mc_molecule_info=mc_molecule_info)
1997  CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, &
1998  mass=mass, nunits=nunits, nunits_tot=nunits_tot, mol_type=mol_type)
1999 ! figure out some bounds for mol_type
2000  start_mol = 1
2001  DO jbox = 1, box_number - 1
2002  start_mol = start_mol + sum(nchains(:, jbox))
2003  END DO
2004  end_mol = start_mol + sum(nchains(:, box_number)) - 1
2005 
2006 ! nullify some pointers
2007  NULLIFY (particles, subsys, molecule_kinds, molecule_kind, &
2008  particles_force, subsys_force)
2009 
2010 ! do some allocation
2011  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  IF (lbias) THEN
2016 
2017 ! grab the coordinates
2018  CALL force_env_get(bias_env, cell=cell, subsys=subsys)
2019  CALL cp_subsys_get(subsys, &
2020  particles=particles, molecule_kinds=molecule_kinds)
2021  molecule_kind => molecule_kinds%els(1)
2022  CALL get_molecule_kind(molecule_kind, natom=natom)
2023  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  CALL force_env_get(force_env, cell=cell, subsys=subsys)
2032  CALL cp_subsys_get(subsys, &
2033  particles=particles, molecule_kinds=molecule_kinds)
2034  molecule_kind => molecule_kinds%els(1)
2035  CALL get_molecule_kind(molecule_kind, natom=natom)
2036  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  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  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  particles%els(target_atom)%r(3))/abc(3))
2054  distance = sqrt(rij(1)**2 + rij(2)**2 + rij(3)**2)
2055  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  IF (move_type == 'in') THEN
2065  moves%avbmc_inin%attempts = &
2066  moves%avbmc_inin%attempts + 1
2067  ELSE
2068  moves%avbmc_inout%attempts = &
2069  moves%avbmc_inout%attempts + 1
2070  END IF
2071  ELSE
2072  IF (move_type == 'in') THEN
2073  moves%avbmc_outin%attempts = &
2074  moves%avbmc_outin%attempts + 1
2075  ELSE
2076  moves%avbmc_outout%attempts = &
2077  moves%avbmc_outout%attempts + 1
2078  END IF
2079  END IF
2080 
2081  IF (lbias) THEN
2082 
2083  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  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  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  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  bias_energy_new = bias_energy_new + bias_energy_old
2127 
2128  ELSE
2129 
2130  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  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  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  target_atom=target_atom)
2168 
2169  END IF
2170 
2171  IF (loverlap) THEN
2172  DEALLOCATE (r_new)
2173 
2174 ! need to reset the old coordinates
2175  IF (lbias) THEN
2176  CALL force_env_get(bias_env, subsys=subsys)
2177  CALL cp_subsys_get(subsys, particles=particles)
2178  ELSE
2179  CALL force_env_get(force_env, subsys=subsys)
2180  CALL cp_subsys_get(subsys, particles=particles)
2181  END IF
2182  DO ipart = 1, nunits_tot(box_number)
2183  particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2184  END DO
2185 
2186  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  IF (lbias) THEN
2194 ! need to give the force_env the coords from the bias_env
2195  CALL force_env_get(force_env, subsys=subsys_force)
2196  CALL cp_subsys_get(subsys_force, particles=particles_force)
2197  CALL force_env_get(bias_env, subsys=subsys)
2198  CALL cp_subsys_get(subsys, particles=particles)
2199  DO ipart = 1, nunits_tot(box_number)
2200  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  calc_force=.false.)
2205  CALL force_env_get(force_env, &
2206  potential_energy=new_energy)
2207 
2208  END IF
2209 
2210  volume_in = 4.0_dp/3.0_dp*pi*(avbmc_rmax(molecule_type)**3 - avbmc_rmin(molecule_type)**3)
2211  volume_out = abc(1)*abc(2)*abc(3) - volume_in
2212 
2213  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  ELSEIF (.NOT. lin .AND. move_type == 'in') THEN
2218  prefactor = (1.0_dp - pbias(molecule_type))*volume_in/(pbias(molecule_type)*volume_out)
2219  ELSE
2220  prefactor = pbias(molecule_type)*volume_out/((1.0_dp - pbias(molecule_type))*volume_in)
2221  END IF
2222 
2223  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  (bias_energy_new - bias_energy_old))
2229 
2230  IF (del_quickstep_energy .GT. exp_max_val) THEN
2231  del_quickstep_energy = max_val
2232  ELSEIF (del_quickstep_energy .LT. exp_min_val) THEN
2233  del_quickstep_energy = 0.0_dp
2234  ELSE
2235  del_quickstep_energy = exp(del_quickstep_energy)
2236  END IF
2237 
2238  w = prefactor*del_quickstep_energy*weight_new/weight_old
2239 
2240  ELSE
2241 
2242 ! AVBMC with CBMC
2243  w = prefactor*weight_new/weight_old
2244  END IF
2245 
2246 ! check if the move is accepted
2247  IF (w .GE. 1.0e0_dp) THEN
2248  rand = 0.0e0_dp
2249  ELSE
2250  IF (ionode) rand = rng_stream%next()
2251  CALL group%bcast(rand, source)
2252  END IF
2253 
2254  IF (rand .LT. w) THEN
2255 
2256 ! accept the move
2257 
2258  IF (lin) THEN
2259  IF (move_type == 'in') THEN
2260  moves%avbmc_inin%successes = &
2261  moves%avbmc_inin%successes + 1
2262  ELSE
2263  moves%avbmc_inout%successes = &
2264  moves%avbmc_inout%successes + 1
2265  END IF
2266  ELSE
2267  IF (move_type == 'in') THEN
2268  moves%avbmc_outin%successes = &
2269  moves%avbmc_outin%successes + 1
2270  ELSE
2271  moves%avbmc_outout%successes = &
2272  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  IF (.NOT. lbias) THEN
2279  new_energy = new_energy + old_energy
2280  END IF
2281 
2282 ! update energies
2283  energy_check = energy_check + (new_energy - old_energy)
2284  old_energy = new_energy
2285 
2286 ! if we're biasing the update the biasing energy
2287  IF (lbias) THEN
2288 ! need to do this outside of the routine
2289  last_bias_energy = bias_energy_new
2290  bias_energy_old = bias_energy_new
2291  END IF
2292 
2293 ! update coordinates
2294  CALL force_env_get(force_env, subsys=subsys)
2295  CALL cp_subsys_get(subsys, particles=particles)
2296  DO ipart = 1, nunits_tot(box_number)
2297  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  IF (lbias) THEN
2302  CALL force_env_get(bias_env, subsys=subsys)
2303  CALL cp_subsys_get(subsys, particles=particles)
2304  DO ipart = 1, nunits_tot(box_number)
2305  particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2306  END DO
2307  CALL cp_subsys_set(subsys, particles=particles)
2308  END IF
2309  CALL force_env_get(force_env, subsys=subsys)
2310  CALL cp_subsys_get(subsys, particles=particles)
2311  DO ipart = 1, nunits_tot(box_number)
2312  particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2313  END DO
2314  CALL cp_subsys_set(subsys, particles=particles)
2315 
2316  END IF
2317 
2318 ! deallocate some stuff
2319  DEALLOCATE (r_new)
2320 ! end the timing
2321  CALL timestop(handle)
2322 
2323  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  SUBROUTINE mc_hmc_move(mc_par, force_env, globenv, moves, move_updates, &
2347  old_energy, box_number, &
2348  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  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  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  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  mc_molecule_info=mc_molecule_info)
2383  CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot)
2384 
2385 ! nullify some pointers
2386  NULLIFY (particles_old, oldsys, hmc_ekin)
2387 
2388 ! do some allocation
2389  ALLOCATE (r(1:3, 1:nunits_tot(box_number)))
2390  ALLOCATE (hmc_ekin)
2391 
2392 ! record the attempt
2393  moves%hmc%attempts = moves%hmc%attempts + 1
2394  move_updates%hmc%attempts = move_updates%hmc%attempts + 1
2395 
2396 ! now let's grab the particle positions
2397  CALL force_env_get(force_env, subsys=oldsys)
2398  CALL cp_subsys_get(oldsys, particles=particles_old)
2399 
2400 ! save the old coordinates
2401  DO iatom = 1, nunits_tot(box_number)
2402  r(1:3, iatom) = particles_old%els(iatom)%r(1:3)
2403  END DO
2404 
2405 ! now run the MD simulation
2406  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  potential_energy=new_energy)
2411 
2412 ! accept or reject the move
2413 ! to prevent overflows
2414  energy_term = new_energy + hmc_ekin%final_ekin - old_energy - hmc_ekin%initial_ekin
2415 
2416  value = -beta*(energy_term)
2417  IF (value .GT. exp_max_val) THEN
2418  w = 10.0_dp
2419  ELSEIF (value .LT. exp_min_val) THEN
2420  w = 0.0_dp
2421  ELSE
2422  w = exp(value)
2423  END IF
2424 
2425  IF (w .GE. 1.0e0_dp) THEN
2426  w = 1.0e0_dp
2427  rand = 0.0e0_dp
2428  ELSE
2429  IF (ionode) rand = rng_stream%next()
2430  CALL group%bcast(rand, source)
2431  END IF
2432 
2433  IF (rand .LT. w) THEN
2434 
2435 ! accept the move
2436  moves%hmc%successes = moves%hmc%successes + 1
2437  move_updates%hmc%successes = move_updates%hmc%successes + 1
2438 
2439 ! update energies
2440  energy_check = energy_check + (new_energy - old_energy)
2441  old_energy = new_energy
2442 
2443  DO iatom = 1, nunits_tot(box_number)
2444  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  DO iatom = 1, nunits_tot(box_number)
2451  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  DEALLOCATE (r)
2458  DEALLOCATE (hmc_ekin)
2459 
2460 ! end the timing
2461  CALL timestop(handle)
2462 
2463  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  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  INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
2502  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  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  CALL timeset(routinen, handle)
2517 
2518 ! nullify some pointers
2519  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  group=group, mc_molecule_info=mc_molecule_info)
2526  CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot, &
2527  nchains=nchains, nunits=nunits, mol_type=mol_type)
2528 
2529 ! find out some bounds for mol_type
2530  start_mol = 1
2531  DO jbox = 1, box_number - 1
2532  start_mol = start_mol + sum(nchains(:, jbox))
2533  END DO
2534  end_mol = start_mol + sum(nchains(:, box_number)) - 1
2535 
2536 ! do some allocation
2537  ALLOCATE (r_old(1:3, 1:nunits_tot(box_number)))
2538 
2539 ! Allocating cluster matrix size
2540  nend = sum(nchains(:, box_number))
2541  ALLOCATE (cluster(nend, nend))
2542  DO ipart = 1, nend
2543  DO jpart = 1, nend
2544  cluster(ipart, jpart) = 0
2545  END DO
2546  END DO
2547 
2548 ! Get cluster information in cluster matrix from cluster_search subroutine
2549  IF (lbias) THEN
2550  CALL cluster_search(mc_par, bias_env, cluster, nchains(:, box_number), &
2551  nunits, mol_type(start_mol:end_mol), total_clus)
2552  ELSE
2553  CALL cluster_search(mc_par, force_env, cluster, nchains(:, box_number), &
2554  nunits, mol_type(start_mol:end_mol), total_clus)
2555  END IF
2556 
2557  IF (lbias) THEN
2558 
2559 ! grab the coordinates
2560  CALL force_env_get(bias_env, subsys=subsys)
2561  CALL cp_subsys_get(subsys, particles=particles)
2562 
2563 ! save the coordinates
2564  DO ipart = 1, nunits_tot(box_number)
2565  r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
2566  END DO
2567 
2568 ! save the energy
2569  bias_energy_old = bias_energy
2570  ELSE
2571 
2572 ! grab the coordinates
2573  CALL force_env_get(force_env, subsys=subsys)
2574  CALL cp_subsys_get(subsys, particles=particles)
2575  END IF
2576 
2577 ! record the attempt
2578  moves%cltrans%attempts = moves%cltrans%attempts + 1
2579  move_updates%cltrans%attempts = move_updates%cltrans%attempts + 1
2580  moves%bias_cltrans%attempts = moves%bias_cltrans%attempts + 1
2581  move_updates%bias_cltrans%attempts = move_updates%bias_cltrans%attempts + 1
2582  IF (.NOT. lbias) THEN
2583  moves%cltrans%qsuccesses = moves%cltrans%qsuccesses + 1
2584  move_updates%cltrans%qsuccesses = move_updates%cltrans%qsuccesses + 1
2585  moves%bias_cltrans%qsuccesses = moves%bias_cltrans%qsuccesses + 1
2586  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  IF (ionode) rand = rng_stream%next()
2591  CALL group%bcast(rand, source)
2592  move_direction = int(3*rand) + 1
2593 
2594 ! call a random number to figure out how far we're moving
2595  IF (ionode) rand = rng_stream%next()
2596  CALL group%bcast(rand, source)
2597  dis_mol = rmcltrans*(rand - 0.5e0_dp)*2.0e0_dp
2598 
2599 ! choosing cluster
2600  IF (ionode) rand = rng_stream%next()
2601  CALL group%bcast(rand, source)
2602  jpart = int(1 + rand*total_clus)
2603 
2604 ! do the cluster move
2605  DO cstart = 1, nend
2606  imol = 0
2607  IF (cluster(jpart, cstart) .NE. 0) THEN
2608  imol = cluster(jpart, cstart)
2609  iunit = 1
2610  DO ipart = 1, imol - 1
2611  nunit = nunits(mol_type(ipart + start_mol - 1))
2612  iunit = iunit + nunit
2613  END DO
2614  nunit = nunits(mol_type(imol + start_mol - 1))
2615  junit = iunit + nunit - 1
2616  DO iparticle = iunit, junit
2617  particles%els(iparticle)%r(move_direction) = &
2618  particles%els(iparticle)%r(move_direction) + dis_mol
2619  END DO
2620  END IF
2621  END DO
2622  CALL cp_subsys_set(subsys, particles=particles)
2623 
2624 !Make cluster matrix null
2625  DO ipart = 1, nend
2626  DO jpart = 1, nend
2627  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  IF (lbias) THEN
2633  CALL cluster_search(mc_par, bias_env, cluster, nchains(:, box_number), &
2634  nunits, mol_type(start_mol:end_mol), total_clusafmo)
2635  ELSE
2636  CALL cluster_search(mc_par, force_env, cluster, nchains(:, box_number), &
2637  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  lreject = .false.
2642  IF (lbias) THEN
2643  CALL check_for_overlap(bias_env, nchains(:, box_number), &
2644  nunits(:), loverlap, mol_type(start_mol:end_mol))
2645  ELSE
2646  CALL check_for_overlap(force_env, nchains(:, box_number), &
2647  nunits(:), loverlap, mol_type(start_mol:end_mol))
2648  IF (loverlap) lreject = .true.
2649  END IF
2650 
2651 ! check if cluster size changes then reject the move
2652  IF (lbias) THEN
2653  IF (total_clusafmo .NE. total_clus) THEN
2654  loverlap = .true.
2655  END IF
2656  ELSE
2657  IF (total_clusafmo .NE. total_clus) THEN
2658  loverlap = .true.
2659  lreject = .true.
2660  END IF
2661  END IF
2662 
2663 ! if we're biasing with a cheaper potential, check for acceptance
2664  IF (lbias) THEN
2665 
2666 ! here's where we bias the moves
2667  IF (loverlap) THEN
2668  w = 0.0e0_dp
2669  ELSE
2670  CALL force_env_calc_energy_force(bias_env, calc_force=.false.)
2671  CALL force_env_get(bias_env, &
2672  potential_energy=bias_energy_new)
2673 ! accept or reject the move based on the Metropolis rule
2674  value = -beta*(bias_energy_new - bias_energy_old)
2675  IF (value .GT. exp_max_val) THEN
2676  w = 10.0_dp
2677  ELSEIF (value .LT. exp_min_val) THEN
2678  w = 0.0_dp
2679  ELSE
2680  w = exp(value)
2681  END IF
2682 
2683  END IF
2684 
2685  IF (w .GE. 1.0e0_dp) THEN
2686  w = 1.0e0_dp
2687  rand = 0.0e0_dp
2688  ELSE
2689  IF (ionode) rand = rng_stream%next()
2690  CALL group%bcast(rand, source)
2691  END IF
2692  IF (rand .LT. w) THEN
2693 
2694 ! accept the move
2695  moves%bias_cltrans%successes = moves%bias_cltrans%successes + 1
2696  move_updates%bias_cltrans%successes = move_updates%bias_cltrans%successes + 1
2697  moves%cltrans%qsuccesses = moves%cltrans%qsuccesses + 1
2698  move_updates%cltrans%successes = &
2699  move_updates%cltrans%successes + 1
2700  moves%qcltrans_dis = moves%qcltrans_dis + abs(dis_mol)
2701  bias_energy = bias_energy + bias_energy_new - &
2702  bias_energy_old
2703 
2704  ELSE
2705 
2706 ! reject the move
2707 ! restore the coordinates
2708  CALL force_env_get(bias_env, subsys=subsys)
2709  CALL cp_subsys_get(subsys, particles=particles)
2710  DO ipart = 1, nunits_tot(box_number)
2711  particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
2712  END DO
2713  CALL cp_subsys_set(subsys, particles=particles)
2714 
2715  END IF
2716 
2717  END IF
2718 
2719 ! deallocate some stuff
2720  DEALLOCATE (cluster)
2721  DEALLOCATE (r_old)
2722 
2723 ! end the timing
2724  CALL timestop(handle)
2725 
2726  END SUBROUTINE mc_cluster_translation
2727 
2728 END MODULE mc_moves
Definition: atom.F:9
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Handles all functions related to the CELL.
Definition: cell_methods.F:15
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Definition: cell_methods.F:85
Handles all functions related to the CELL.
Definition: cell_types.F:15
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition: cell_types.F:559
subroutine, public cell_clone(cell_in, cell_out, tag)
Clone cell variable.
Definition: cell_types.F:107
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition: cell_types.F:195
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
types that represent a subsys, i.e. a part of the system
subroutine, public cp_subsys_set(subsys, atomic_kinds, particles, local_particles, molecules, molecule_kinds, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, results, cell)
sets various propreties of the subsys
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition: fft_lib.F:7
Interface for the force calculations.
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
Interface for the force calculations.
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
integer, parameter, public use_fist_force
Define type storing the global information of a run. Keep the amount of stored data small....
Definition: global_types.F:21
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
contains miscellaneous subroutines used in the Monte Carlo runs,mostly geared towards changes in coor...
subroutine, public get_center_of_mass(coordinates, natom, center_of_mass, mass)
calculates the center of mass of a given molecule
subroutine, public generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, exp_min_val, nswapmoves, rosenbluth_weight, start_atom, natoms_tot, nunits, nunits_mol, mass, loverlap, choosen_energy, old_energy, ionode, lremove, mol_type, nchains, source, group, rng_stream, avbmc_atom, rmin, rmax, move_type, target_atom)
takes the last molecule in a force environment and moves it around to different center of mass positi...
subroutine, public check_for_overlap(force_env, nchains, nunits, loverlap, mol_type, cell_length, molecule_number)
looks for overlaps (intermolecular distances less than rmin)
subroutine, public create_discrete_array(cell, discrete_array, step_size)
generates an array that tells us which sides of the simulation cell we can increase or decrease using...
subroutine, public cluster_search(mc_par, force_env, cluster, nchains, nunits, mol_type, total_clus)
determine the number of cluster present in the given configuration based on the rclus value
the various moves in Monte Carlo (MC) simulations, including change of internal conformation,...
Definition: mc_moves.F:16
subroutine, public mc_molecule_rotation(mc_par, force_env, bias_env, moves, move_updates, box_number, start_atom, molecule_type, bias_energy, lreject, rng_stream)
rotates the given molecule randomly around the x,y, or z axis... only works for water at the moment
Definition: mc_moves.F:660
subroutine, public mc_avbmc_move(mc_par, force_env, bias_env, moves, energy_check, r_old, old_energy, start_atom_swap, target_atom, molecule_type, box_number, bias_energy_old, last_bias_energy, move_type, rng_stream)
performs either a bond or angle change move for a given molecule
Definition: mc_moves.F:1950
subroutine, public mc_volume_move(mc_par, force_env, moves, move_updates, old_energy, box_number, energy_check, r_old, iw, discrete_array, rng_stream)
performs a Monte Carlo move that alters the volume of the simulation box
Definition: mc_moves.F:962
subroutine, public mc_molecule_translation(mc_par, force_env, bias_env, moves, move_updates, start_atom, box_number, bias_energy, molecule_type, lreject, rng_stream)
translates the given molecule randomly in either the x,y, or z direction
Definition: mc_moves.F:436
subroutine, public mc_cluster_translation(mc_par, force_env, bias_env, moves, move_updates, box_number, bias_energy, lreject, rng_stream)
translates the cluster randomly in either the x,y, or z direction
Definition: mc_moves.F:2487
subroutine, public mc_hmc_move(mc_par, force_env, globenv, moves, move_updates, old_energy, box_number, energy_check, r_old, rng_stream)
performs a hybrid Monte Carlo move that runs a short MD sequence
Definition: mc_moves.F:2349
subroutine, public mc_conformation_change(mc_par, force_env, bias_env, moves, move_updates, start_atom, molecule_type, box_number, bias_energy, move_type, lreject, rng_stream)
performs either a bond or angle change move for a given molecule
Definition: mc_moves.F:140
holds all the structure types needed for Monte Carlo, except the mc_environment_type
Definition: mc_types.F:15
subroutine, public get_mc_molecule_info(mc_molecule_info, nmol_types, nchain_total, nboxes, names, conf_prob, nchains, nunits, mol_type, nunits_tot, in_box, atom_names, mass)
...
Definition: mc_types.F:554
subroutine, public get_mc_par(mc_par, nstep, nvirial, iuptrans, iupcltrans, iupvolume, nmoves, nswapmoves, rm, cl, diff, nstart, source, group, lbias, ionode, lrestart, lstop, rmvolume, rmcltrans, rmbond, rmangle, rmrot, rmtrans, temperature, pressure, rclus, BETA, pmswap, pmvolume, pmtraion, pmtrans, pmcltrans, ensemble, PROGRAM, restart_file_name, molecules_file, moves_file, coords_file, energy_file, displacement_file, cell_file, dat_file, data_file, box2_file, fft_lib, iprint, rcut, ldiscrete, discrete_step, pmavbmc, pbias, avbmc_atom, avbmc_rmin, avbmc_rmax, rmdihedral, input_file, mc_molecule_info, pmswap_mol, pmavbmc_mol, pmtrans_mol, pmrot_mol, pmtraion_mol, mc_input_file, mc_bias_file, pmvol_box, pmclus_box, virial_temps, exp_min_val, exp_max_val, min_val, max_val, eta, pmhmc, pmhmc_box, lhmc, rand2skip)
...
Definition: mc_types.F:405
Perform a molecular dynamics (MD) run using QUICKSTEP.
Definition: md_run.F:14
subroutine, public qs_mol_dyn(force_env, globenv, averages, rm_restart_info, hmc_e_initial, hmc_e_final, mdctrl)
Main driver module for Molecular Dynamics.
Definition: md_run.F:121
Interface to the message passing library MPI.
represent a simple array based list of the given type
Define the molecule kind structure types and the corresponding functionality.
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
represent a simple array based list of the given type
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144