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