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