(git:79ce675)
Loading...
Searching...
No Matches
mc_ge_moves.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief contains the Monte Carlo moves that can handle more than one
10!> box, including the Quickstep move, a volume swap between boxes,
11!> and a particle swap between boxes
12!> \par History
13!> MJM (07.28.2005): make the Quickstep move general, and changed
14!> the swap and volume moves to work with the
15!> CP2K classical routines
16!> \author Matthew J. McGrath (01.25.2004)
17! **************************************************************************************************
18MODULE mc_ge_moves
19 USE cell_methods, ONLY: cell_create
20 USE cell_types, ONLY: cell_clone,&
21 cell_p_type,&
22 cell_release,&
23 cell_type,&
24 get_cell
25 USE cp_subsys_types, ONLY: cp_subsys_get,&
26 cp_subsys_p_type,&
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_p_type,&
32 force_env_release
33 USE input_constants, ONLY: dump_xmol
34 USE input_section_types, ONLY: section_type
35 USE kinds, ONLY: default_string_length,&
36 dp
37 USE mc_control, ONLY: mc_create_force_env
38 USE mc_coordinates, ONLY: check_for_overlap,&
39 generate_cbmc_swap_config,&
40 get_center_of_mass
41 USE mc_misc, ONLY: mc_make_dat_file_new
42 USE mc_move_control, ONLY: move_q_reinit,&
43 q_move_accept
44 USE mc_types, ONLY: &
45 get_mc_molecule_info, get_mc_par, mc_determine_molecule_info, mc_input_file_type, &
46 mc_molecule_info_destroy, mc_molecule_info_type, mc_moves_p_type, &
47 mc_simulation_parameters_p_type, set_mc_par
48 USE message_passing, ONLY: mp_comm_type,&
49 mp_para_env_type
50 USE parallel_rng_types, ONLY: rng_stream_type
51 USE particle_list_types, ONLY: particle_list_p_type,&
52 particle_list_type
53 USE particle_methods, ONLY: write_particle_coordinates
54 USE physcon, ONLY: angstrom
55#include "../../base/base_uses.f90"
56
57 IMPLICIT NONE
58
59 PRIVATE
60
61! *** Global parameters ***
62
63 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_ge_moves'
64
65 PUBLIC :: mc_ge_volume_move, mc_ge_swap_move, &
66 mc_quickstep_move
67
68CONTAINS
69
70! **************************************************************************************************
71!> \brief computes the acceptance of a series of biased or unbiased moves
72!> (translation, rotation, conformational changes)
73!> \param mc_par the mc parameters for the force envs of the boxes
74!> \param force_env the force environments for the boxes
75!> \param bias_env the force environments with the biasing potential for the boxes
76!> \param moves the structure that keeps track of how many moves have been
77!> accepted/rejected for both boxes
78!> \param lreject automatically rejects the move (used when an overlap occurs in
79!> the sequence of moves)
80!> \param move_updates the structure that keeps track of how many moves have
81!> been accepted/rejected since the last time the displacements
82!> were updated for both boxes
83!> \param energy_check the running total of how much the energy has changed
84!> since the initial configuration
85!> \param r_old the coordinates of the last accepted move before the sequence
86!> whose acceptance is determined by this call
87!> \param nnstep the Monte Carlo step we're on
88!> \param old_energy the energy of the last accepted move involving the full potential
89!> \param bias_energy_new the energy of the current configuration involving the bias potential
90!> \param last_bias_energy ...
91!> \param nboxes the number of boxes (force environments) in the system
92!> \param box_flag indicates if a move has been tried in a given box..if not, we don't
93!> recompute the energy
94!> \param subsys the pointers for the particle subsystems of both boxes
95!> \param particles the pointers for the particle sets
96!> \param rng_stream the stream we pull random numbers from
97!> \param unit_conv ...
98!> \author MJM
99! **************************************************************************************************
100 SUBROUTINE mc_quickstep_move(mc_par, force_env, bias_env, moves, &
101 lreject, move_updates, energy_check, r_old, &
102 nnstep, old_energy, bias_energy_new, last_bias_energy, &
103 nboxes, box_flag, subsys, particles, rng_stream, &
104 unit_conv)
105
106 TYPE(mc_simulation_parameters_p_type), &
107 DIMENSION(:), POINTER :: mc_par
108 TYPE(force_env_p_type), DIMENSION(:), POINTER :: force_env, bias_env
109 TYPE(mc_moves_p_type), DIMENSION(:, :), POINTER :: moves
110 LOGICAL, INTENT(IN) :: lreject
111 TYPE(mc_moves_p_type), DIMENSION(:, :), POINTER :: move_updates
112 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: energy_check
113 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: r_old
114 INTEGER, INTENT(IN) :: nnstep
115 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: old_energy, bias_energy_new, &
116 last_bias_energy
117 INTEGER, INTENT(IN) :: nboxes
118 INTEGER, DIMENSION(:), INTENT(IN) :: box_flag
119 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsys
120 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
121 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
122 REAL(kind=dp), INTENT(IN) :: unit_conv
123
124 CHARACTER(len=*), PARAMETER :: routinen = 'mc_Quickstep_move'
125
126 INTEGER :: end_mol, handle, ibox, iparticle, &
127 iprint, itype, jbox, nmol_types, &
128 source, start_mol
129 INTEGER, DIMENSION(:, :), POINTER :: nchains
130 INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
131 INTEGER, DIMENSION(1:nboxes) :: diff
132 LOGICAL :: ionode, lbias, loverlap
133 REAL(kind=dp) :: beta, energies, rand, w
134 REAL(kind=dp), DIMENSION(1:nboxes) :: bias_energy_old, new_energy
135 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsys_bias
136 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
137 TYPE(mp_comm_type) :: group
138 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles_bias
139
140! begin the timing of the subroutine
141
142 CALL timeset(routinen, handle)
143
144 NULLIFY (subsys_bias, particles_bias)
145
146! get a bunch of data from mc_par
147 CALL get_mc_par(mc_par(1)%mc_par, ionode=ionode, lbias=lbias, &
148 beta=beta, diff=diff(1), source=source, group=group, &
149 iprint=iprint, &
150 mc_molecule_info=mc_molecule_info)
151 CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types, &
152 nchains=nchains, nunits_tot=nunits_tot, nunits=nunits, mol_type=mol_type)
153
154 IF (nboxes > 1) THEN
155 DO ibox = 2, nboxes
156 CALL get_mc_par(mc_par(ibox)%mc_par, diff=diff(ibox))
157 END DO
158 END IF
159
160! allocate some stuff
161 ALLOCATE (subsys_bias(1:nboxes))
162 ALLOCATE (particles_bias(1:nboxes))
163
164! record the attempt...we really only care about molecule type 1 and box
165! type 1, since the acceptance will be identical for all boxes and molecules
166 moves(1, 1)%moves%Quickstep%attempts = &
167 moves(1, 1)%moves%Quickstep%attempts + 1
168
169! grab the coordinates for the force_env
170 DO ibox = 1, nboxes
171 CALL force_env_get(force_env(ibox)%force_env, &
172 subsys=subsys(ibox)%subsys)
173 CALL cp_subsys_get(subsys(ibox)%subsys, &
174 particles=particles(ibox)%list)
175 END DO
176
177! calculate the new energy of the system...if we're biasing,
178! force_env hasn't changed but bias_env has
179 DO ibox = 1, nboxes
180 IF (box_flag(ibox) == 1) THEN
181 IF (lbias) THEN
182! grab the coords from bias_env and put them into force_env
183 CALL force_env_get(bias_env(ibox)%force_env, &
184 subsys=subsys_bias(ibox)%subsys)
185 CALL cp_subsys_get(subsys_bias(ibox)%subsys, &
186 particles=particles_bias(ibox)%list)
187
188 DO iparticle = 1, nunits_tot(ibox)
189 particles(ibox)%list%els(iparticle)%r(1:3) = &
190 particles_bias(ibox)%list%els(iparticle)%r(1:3)
191 END DO
192
193 CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
194 calc_force=.false.)
195 CALL force_env_get(force_env(ibox)%force_env, &
196 potential_energy=new_energy(ibox))
197 ELSE
198 IF (.NOT. lreject) THEN
199 CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
200 calc_force=.false.)
201 CALL force_env_get(force_env(ibox)%force_env, &
202 potential_energy=new_energy(ibox))
203 END IF
204 END IF
205 ELSE
206 new_energy(ibox) = old_energy(ibox)
207 END IF
208
209 END DO
210
211! accept or reject the move based on Metropolis or the Iftimie rule
212 IF (ionode) THEN
213
214! write them out in case something bad happens
215 IF (mod(nnstep, iprint) == 0) THEN
216 DO ibox = 1, nboxes
217 IF (sum(nchains(:, ibox)) == 0) THEN
218 WRITE (diff(ibox), *) nnstep
219 WRITE (diff(ibox), *) nchains(:, ibox)
220 ELSE
221 WRITE (diff(ibox), *) nnstep
222 CALL write_particle_coordinates( &
223 particles(ibox)%list%els, &
224 diff(ibox), dump_xmol, 'POS', 'TRIAL', &
225 unit_conv=unit_conv)
226 END IF
227 END DO
228 END IF
229 END IF
230
231 IF (.NOT. lreject) THEN
232 IF (lbias) THEN
233
234 DO ibox = 1, nboxes
235! look for overlap
236 IF (sum(nchains(:, ibox)) /= 0) THEN
237! find the molecule bounds
238 start_mol = 1
239 DO jbox = 1, ibox - 1
240 start_mol = start_mol + sum(nchains(:, jbox))
241 END DO
242 end_mol = start_mol + sum(nchains(:, ibox)) - 1
243 CALL check_for_overlap(bias_env(ibox)%force_env, &
244 nchains(:, ibox), nunits(:), loverlap, mol_type(start_mol:end_mol))
245 IF (loverlap) THEN
246 cpabort('Quickstep move found an overlap in the old config')
247 END IF
248 END IF
249 bias_energy_old(ibox) = last_bias_energy(ibox)
250 END DO
251
252 energies = -beta*((sum(new_energy(:)) - sum(bias_energy_new(:))) &
253 - (sum(old_energy(:)) - sum(bias_energy_old(:))))
254
255! used to prevent over and underflows
256 IF (energies >= -1.0e-8) THEN
257 w = 1.0_dp
258 ELSE IF (energies <= -500.0_dp) THEN
259 w = 0.0_dp
260 ELSE
261 w = exp(energies)
262 END IF
263
264 IF (ionode) THEN
265 DO ibox = 1, nboxes
266 WRITE (diff(ibox), *) nnstep, new_energy(ibox) - &
267 old_energy(ibox), &
268 bias_energy_new(ibox) - bias_energy_old(ibox)
269 END DO
270 END IF
271 ELSE
272 energies = -beta*(sum(new_energy(:)) - sum(old_energy(:)))
273! used to prevent over and underflows
274 IF (energies >= 0.0_dp) THEN
275 w = 1.0_dp
276 ELSE IF (energies <= -500.0_dp) THEN
277 w = 0.0_dp
278 ELSE
279 w = exp(energies)
280 END IF
281 END IF
282 ELSE
283 w = 0.0e0_dp
284 END IF
285 IF (w >= 1.0e0_dp) THEN
286 w = 1.0e0_dp
287 rand = 0.0e0_dp
288 ELSE
289 IF (ionode) rand = rng_stream%next()
290 CALL group%bcast(rand, source)
291 END IF
292
293 IF (rand < w) THEN
294
295! accept the move
296 moves(1, 1)%moves%Quickstep%successes = &
297 moves(1, 1)%moves%Quickstep%successes + 1
298
299 DO ibox = 1, nboxes
300! remember what kind of move we did for lbias=.false.
301 IF (.NOT. lbias) THEN
302 DO itype = 1, nmol_types
303 CALL q_move_accept(moves(itype, ibox)%moves, .true.)
304 CALL q_move_accept(move_updates(itype, ibox)%moves, .true.)
305
306! reset the counters
307 CALL move_q_reinit(moves(itype, ibox)%moves, .true.)
308 CALL move_q_reinit(move_updates(itype, ibox)%moves, .true.)
309 END DO
310 END IF
311
312 DO itype = 1, nmol_types
313! we need to record all accepted moves since last Quickstep calculation
314 CALL q_move_accept(moves(itype, ibox)%moves, .false.)
315 CALL q_move_accept(move_updates(itype, ibox)%moves, .false.)
316
317! reset the counters
318 CALL move_q_reinit(moves(itype, ibox)%moves, .false.)
319 CALL move_q_reinit(move_updates(itype, ibox)%moves, .false.)
320 END DO
321
322! update energies
323 energy_check(ibox) = energy_check(ibox) + &
324 (new_energy(ibox) - old_energy(ibox))
325 old_energy(ibox) = new_energy(ibox)
326
327 END DO
328
329 IF (lbias) THEN
330 DO ibox = 1, nboxes
331 last_bias_energy(ibox) = bias_energy_new(ibox)
332 END DO
333 END IF
334
335! update coordinates
336 DO ibox = 1, nboxes
337 IF (nunits_tot(ibox) /= 0) THEN
338 DO iparticle = 1, nunits_tot(ibox)
339 r_old(1:3, iparticle, ibox) = &
340 particles(ibox)%list%els(iparticle)%r(1:3)
341 END DO
342 END IF
343 END DO
344
345 ELSE
346
347 ! reject the move
348 DO ibox = 1, nboxes
349 DO itype = 1, nmol_types
350 CALL move_q_reinit(moves(itype, ibox)%moves, .false.)
351 CALL move_q_reinit(move_updates(itype, ibox)%moves, .false.)
352 IF (.NOT. lbias) THEN
353! reset the counters
354 CALL move_q_reinit(moves(itype, ibox)%moves, .true.)
355 CALL move_q_reinit(move_updates(itype, ibox)%moves, .true.)
356 END IF
357 END DO
358
359 END DO
360
361 IF (.NOT. ionode) r_old(:, :, :) = 0.0e0_dp
362
363! coodinates changed, so we need to broadcast those, even for the lbias
364! case since bias_env needs to have the same coords as force_env
365 CALL group%bcast(r_old, source)
366
367 DO ibox = 1, nboxes
368 DO iparticle = 1, nunits_tot(ibox)
369 particles(ibox)%list%els(iparticle)%r(1:3) = &
370 r_old(1:3, iparticle, ibox)
371 IF (lbias .AND. box_flag(ibox) == 1) THEN
372 particles_bias(ibox)%list%els(iparticle)%r(1:3) = &
373 r_old(1:3, iparticle, ibox)
374 END IF
375 END DO
376 END DO
377
378! need to reset the energies of the biasing potential
379 IF (lbias) THEN
380 DO ibox = 1, nboxes
381 bias_energy_new(ibox) = last_bias_energy(ibox)
382 END DO
383 END IF
384
385 END IF
386
387! make sure the coordinates are transferred
388 DO ibox = 1, nboxes
389 CALL cp_subsys_set(subsys(ibox)%subsys, &
390 particles=particles(ibox)%list)
391 IF (lbias .AND. box_flag(ibox) == 1) THEN
392 CALL cp_subsys_set(subsys_bias(ibox)%subsys, &
393 particles=particles_bias(ibox)%list)
394 END IF
395 END DO
396
397 ! deallocate some stuff
398 DEALLOCATE (subsys_bias)
399 DEALLOCATE (particles_bias)
400
401! end the timing
402 CALL timestop(handle)
403
404 END SUBROUTINE mc_quickstep_move
405
406! **************************************************************************************************
407!> \brief attempts a swap move between two simulation boxes
408!> \param mc_par the mc parameters for the force envs of the boxes
409!> \param force_env the force environments for the boxes
410!> \param bias_env the force environments used to bias moves for the boxes
411!> \param moves the structure that keeps track of how many moves have been
412!> accepted/rejected for both boxes
413!> \param energy_check the running total of how much the energy has changed
414!> since the initial configuration
415!> \param r_old the coordinates of the last accepted move involving a
416!> full potential calculation for both boxes
417!> \param old_energy the energy of the last accepted move involving a
418!> a full potential calculation
419!> \param input_declaration ...
420!> \param para_env the parallel environment for this simulation
421!> \param bias_energy_old the energies of both boxes computed using the biasing
422!> potential
423!> \param last_bias_energy the energy for the biased simulations
424!> \param rng_stream the stream we pull random numbers from
425!> \author MJM
426! **************************************************************************************************
427 SUBROUTINE mc_ge_swap_move(mc_par, force_env, bias_env, moves, &
428 energy_check, r_old, old_energy, input_declaration, &
429 para_env, bias_energy_old, last_bias_energy, &
430 rng_stream)
431
432 TYPE(mc_simulation_parameters_p_type), &
433 DIMENSION(:), POINTER :: mc_par
434 TYPE(force_env_p_type), DIMENSION(:), POINTER :: force_env, bias_env
435 TYPE(mc_moves_p_type), DIMENSION(:, :), POINTER :: moves
436 REAL(kind=dp), DIMENSION(1:2), INTENT(INOUT) :: energy_check
437 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: r_old
438 REAL(kind=dp), DIMENSION(1:2), INTENT(INOUT) :: old_energy
439 TYPE(section_type), POINTER :: input_declaration
440 TYPE(mp_para_env_type), POINTER :: para_env
441 REAL(kind=dp), DIMENSION(1:2), INTENT(INOUT) :: bias_energy_old, last_bias_energy
442 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
443
444 CHARACTER(len=*), PARAMETER :: routinen = 'mc_ge_swap_move'
445
446 CHARACTER(default_string_length), ALLOCATABLE, &
447 DIMENSION(:) :: atom_names_insert, atom_names_remove
448 CHARACTER(default_string_length), &
449 DIMENSION(:, :), POINTER :: atom_names
450 CHARACTER(LEN=200) :: fft_lib
451 CHARACTER(LEN=40), DIMENSION(1:2) :: dat_file
452 INTEGER :: end_mol, handle, iatom, ibox, idim, iiatom, imolecule, ins_atoms, insert_box, &
453 ipart, itype, jbox, molecule_type, nmol_types, nswapmoves, print_level, rem_atoms, &
454 remove_box, source, start_atom_ins, start_atom_rem, start_mol
455 INTEGER, DIMENSION(:), POINTER :: mol_type, mol_type_test, nunits, &
456 nunits_tot
457 INTEGER, DIMENSION(:, :), POINTER :: nchains, nchains_test
458 LOGICAL :: ionode, lbias, loverlap, loverlap_ins, &
459 loverlap_rem
460 REAL(dp), DIMENSION(:), POINTER :: eta_insert, eta_remove, pmswap_mol
461 REAL(dp), DIMENSION(:, :), POINTER :: insert_coords, remove_coords
462 REAL(kind=dp) :: beta, del_quickstep_energy, exp_max_val, exp_min_val, max_val, min_val, &
463 prefactor, rand, rdum, vol_insert, vol_remove, w, weight_new, weight_old
464 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cbmc_energies, r_cbmc, r_insert_mol
465 REAL(kind=dp), DIMENSION(1:2) :: bias_energy_new, new_energy
466 REAL(kind=dp), DIMENSION(1:3) :: abc_insert, abc_remove, center_of_mass, &
467 displace_molecule, pos_insert
468 REAL(kind=dp), DIMENSION(:, :), POINTER :: mass
469 TYPE(cell_type), POINTER :: cell_insert, cell_remove
470 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: oldsys
471 TYPE(cp_subsys_type), POINTER :: insert_sys, remove_sys
472 TYPE(force_env_p_type), DIMENSION(:), POINTER :: test_env, test_env_bias
473 TYPE(mc_input_file_type), POINTER :: mc_bias_file, mc_input_file
474 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info, mc_molecule_info_test
475 TYPE(mp_comm_type) :: group
476 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles_old
477 TYPE(particle_list_type), POINTER :: particles_insert, particles_remove
478
479! begin the timing of the subroutine
480
481 CALL timeset(routinen, handle)
482
483! reset the overlap flag
484 loverlap = .false.
485
486! nullify some pointers
487 NULLIFY (particles_old, mol_type, mol_type_test, mc_input_file, mc_bias_file)
488 NULLIFY (oldsys, atom_names, pmswap_mol, insert_coords, remove_coords)
489 NULLIFY (eta_insert, eta_remove)
490
491! grab some stuff from mc_par
492 CALL get_mc_par(mc_par(1)%mc_par, ionode=ionode, beta=beta, &
493 max_val=max_val, min_val=min_val, exp_max_val=exp_max_val, &
494 exp_min_val=exp_min_val, nswapmoves=nswapmoves, group=group, source=source, &
495 lbias=lbias, dat_file=dat_file(1), fft_lib=fft_lib, &
496 mc_molecule_info=mc_molecule_info, pmswap_mol=pmswap_mol)
497 CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, &
498 nunits=nunits, nunits_tot=nunits_tot, nmol_types=nmol_types, &
499 atom_names=atom_names, mass=mass, mol_type=mol_type)
500
501 print_level = 1
502
503 CALL get_mc_par(mc_par(2)%mc_par, dat_file=dat_file(2))
504
505! allocate some stuff
506 ALLOCATE (oldsys(1:2))
507 ALLOCATE (particles_old(1:2))
508
509! get the old coordinates
510 DO ibox = 1, 2
511 CALL force_env_get(force_env(ibox)%force_env, &
512 subsys=oldsys(ibox)%subsys)
513 CALL cp_subsys_get(oldsys(ibox)%subsys, &
514 particles=particles_old(ibox)%list)
515 END DO
516
517! choose a direction to swap
518 IF (ionode) rand = rng_stream%next()
519 CALL group%bcast(rand, source)
520
521 IF (rand <= 0.50e0_dp) THEN
522 remove_box = 1
523 insert_box = 2
524 ELSE
525 remove_box = 2
526 insert_box = 1
527 END IF
528
529! now assign the eta values for the insert and remove boxes
530 CALL get_mc_par(mc_par(remove_box)%mc_par, eta=eta_remove)
531 CALL get_mc_par(mc_par(insert_box)%mc_par, eta=eta_insert)
532
533!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! testing
534! remove_box=2
535! insert_box=1
536
537! now choose a molecule type at random
538 IF (ionode) rand = rng_stream%next()
539 CALL group%bcast(rand, source)
540 DO itype = 1, nmol_types
541 IF (rand < pmswap_mol(itype)) THEN
542 molecule_type = itype
543 EXIT
544 END IF
545 END DO
546
547! record the attempt for the box the particle is to be inserted into
548 moves(molecule_type, insert_box)%moves%swap%attempts = &
549 moves(molecule_type, insert_box)%moves%swap%attempts + 1
550
551! now choose a random molecule to remove from the removal box, checking
552! to make sure the box isn't empty
553 IF (nchains(molecule_type, remove_box) == 0) THEN
554 loverlap = .true.
555 moves(molecule_type, insert_box)%moves%empty = &
556 moves(molecule_type, insert_box)%moves%empty + 1
557 ELSE
558
559 IF (ionode) rand = rng_stream%next()
560 CALL group%bcast(rand, source)
561 imolecule = ceiling(rand*nchains(molecule_type, remove_box))
562! figure out the atom number this molecule starts on
563 start_atom_rem = 1
564 DO itype = 1, nmol_types
565 IF (itype == molecule_type) THEN
566 start_atom_rem = start_atom_rem + (imolecule - 1)*nunits(itype)
567 EXIT
568 ELSE
569 start_atom_rem = start_atom_rem + nchains(itype, remove_box)*nunits(itype)
570 END IF
571 END DO
572
573! check for overlap
574 start_mol = 1
575 DO jbox = 1, remove_box - 1
576 start_mol = start_mol + sum(nchains(:, jbox))
577 END DO
578 end_mol = start_mol + sum(nchains(:, remove_box)) - 1
579 CALL check_for_overlap(force_env(remove_box)%force_env, &
580 nchains(:, remove_box), nunits, loverlap, mol_type(start_mol:end_mol))
581 IF (loverlap) CALL cp_abort(__location__, &
582 'CBMC swap move found an overlap in the old remove config')
583 start_mol = 1
584 DO jbox = 1, insert_box - 1
585 start_mol = start_mol + sum(nchains(:, jbox))
586 END DO
587 end_mol = start_mol + sum(nchains(:, insert_box)) - 1
588 CALL check_for_overlap(force_env(insert_box)%force_env, &
589 nchains(:, insert_box), nunits, loverlap, mol_type(start_mol:end_mol))
590 IF (loverlap) CALL cp_abort(__location__, &
591 'CBMC swap move found an overlap in the old insert config')
592 END IF
593
594 IF (loverlap) THEN
595 DEALLOCATE (oldsys)
596 DEALLOCATE (particles_old)
597 CALL timestop(handle)
598 RETURN
599 END IF
600
601! figure out how many atoms will be in each box after the move
602 ins_atoms = nunits_tot(insert_box) + nunits(molecule_type)
603 rem_atoms = nunits_tot(remove_box) - nunits(molecule_type)
604! now allocate the arrays that will hold the coordinates and the
605! atom name, for writing to the dat file
606 IF (rem_atoms == 0) THEN
607 ALLOCATE (remove_coords(1:3, 1:nunits(1)))
608 ALLOCATE (atom_names_remove(1:nunits(1)))
609 ELSE
610 ALLOCATE (remove_coords(1:3, 1:rem_atoms))
611 ALLOCATE (atom_names_remove(1:rem_atoms))
612 END IF
613 ALLOCATE (insert_coords(1:3, 1:ins_atoms))
614 ALLOCATE (atom_names_insert(1:ins_atoms))
615
616! grab the cells for later...acceptance and insertion
617 IF (lbias) THEN
618 CALL force_env_get(bias_env(insert_box)%force_env, &
619 cell=cell_insert)
620 CALL force_env_get(bias_env(remove_box)%force_env, &
621 cell=cell_remove)
622 ELSE
623 CALL force_env_get(force_env(insert_box)%force_env, &
624 cell=cell_insert)
625 CALL force_env_get(force_env(remove_box)%force_env, &
626 cell=cell_remove)
627 END IF
628 CALL get_cell(cell_remove, abc=abc_remove, deth=vol_remove)
629 CALL get_cell(cell_insert, abc=abc_insert, deth=vol_insert)
630
631 IF (ionode) THEN
632! choose an insertion point
633 DO idim = 1, 3
634 rand = rng_stream%next()
635 pos_insert(idim) = rand*abc_insert(idim)
636 END DO
637 END IF
638 CALL group%bcast(pos_insert, source)
639
640! allocate some arrays we'll be using
641 ALLOCATE (r_insert_mol(1:3, 1:nunits(molecule_type)))
642
643 iiatom = 1
644 DO iatom = start_atom_rem, start_atom_rem + nunits(molecule_type) - 1
645 r_insert_mol(1:3, iiatom) = &
646 particles_old(remove_box)%list%els(iatom)%r(1:3)
647 iiatom = iiatom + 1
648 END DO
649
650! find the center of mass of the molecule
651 CALL get_center_of_mass(r_insert_mol(:, :), nunits(molecule_type), &
652 center_of_mass(:), mass(:, molecule_type))
653
654! move the center of mass to the insertion point
655 displace_molecule(1:3) = pos_insert(1:3) - center_of_mass(1:3)
656 DO iatom = 1, nunits(molecule_type)
657 r_insert_mol(1:3, iatom) = r_insert_mol(1:3, iatom) + &
658 displace_molecule(1:3)
659 END DO
660
661! prepare the insertion coordinates to be written to the .dat file so
662! we can create a new force environment...remember there is still a particle
663! in the box even if nchain=0
664 IF (sum(nchains(:, insert_box)) == 0) THEN
665 DO iatom = 1, nunits(molecule_type)
666 insert_coords(1:3, iatom) = r_insert_mol(1:3, iatom)
667 atom_names_insert(iatom) = &
668 particles_old(remove_box)%list%els(start_atom_rem + iatom - 1)%atomic_kind%name
669 END DO
670 start_atom_ins = 1
671 ELSE
672! the problem is I can't just tack the new molecule on to the end,
673! because of reading in the dat_file...the topology stuff will crash
674! if the molecules aren't all grouped together, so I have to insert it
675! at the end of the section of molecules with the same type, then
676! remember the start number for the CBMC stuff
677 start_atom_ins = 1
678 DO itype = 1, nmol_types
679 start_atom_ins = start_atom_ins + &
680 nchains(itype, insert_box)*nunits(itype)
681 IF (itype == molecule_type) EXIT
682 END DO
683
684 DO iatom = 1, start_atom_ins - 1
685 insert_coords(1:3, iatom) = &
686 particles_old(insert_box)%list%els(iatom)%r(1:3)
687 atom_names_insert(iatom) = &
688 particles_old(insert_box)%list%els(iatom)%atomic_kind%name
689 END DO
690 iiatom = 1
691 DO iatom = start_atom_ins, start_atom_ins + nunits(molecule_type) - 1
692 insert_coords(1:3, iatom) = r_insert_mol(1:3, iiatom)
693 atom_names_insert(iatom) = atom_names(iiatom, molecule_type)
694 iiatom = iiatom + 1
695 END DO
696 DO iatom = start_atom_ins + nunits(molecule_type), ins_atoms
697 insert_coords(1:3, iatom) = &
698 particles_old(insert_box)%list%els(iatom - nunits(molecule_type))%r(1:3)
699 atom_names_insert(iatom) = &
700 particles_old(insert_box)%list%els(iatom - nunits(molecule_type))%atomic_kind%name
701 END DO
702 END IF
703
704! fold the coordinates into the box and check for overlaps
705 start_mol = 1
706 DO jbox = 1, insert_box - 1
707 start_mol = start_mol + sum(nchains(:, jbox))
708 END DO
709 end_mol = start_mol + sum(nchains(:, insert_box)) - 1
710
711! make the .dat file
712 IF (ionode) THEN
713
714 nchains(molecule_type, insert_box) = nchains(molecule_type, insert_box) + 1
715 IF (lbias) THEN
716 CALL get_mc_par(mc_par(insert_box)%mc_par, mc_bias_file=mc_bias_file)
717 CALL mc_make_dat_file_new(insert_coords(:, :), atom_names_insert(:), ins_atoms, &
718 abc_insert(:), dat_file(insert_box), nchains(:, insert_box), &
719 mc_bias_file)
720 ELSE
721 CALL get_mc_par(mc_par(insert_box)%mc_par, mc_input_file=mc_input_file)
722 CALL mc_make_dat_file_new(insert_coords(:, :), atom_names_insert(:), ins_atoms, &
723 abc_insert(:), dat_file(insert_box), nchains(:, insert_box), &
724 mc_input_file)
725 END IF
726 nchains(molecule_type, insert_box) = nchains(molecule_type, insert_box) - 1
727
728 END IF
729
730! now do the same for the removal box...be careful not to make an empty box
731 IF (rem_atoms == 0) THEN
732 DO iatom = 1, nunits(molecule_type)
733 remove_coords(1:3, iatom) = r_insert_mol(1:3, iatom)
734 atom_names_remove(iatom) = atom_names(iatom, molecule_type)
735 END DO
736
737! need to adjust nchains, because otherwise if we are removing a molecule type
738! that is not the first molecule, the dat file will have two molecules in it but
739! only the coordinates for one
740 nchains(molecule_type, remove_box) = nchains(molecule_type, remove_box) - 1
741 IF (ionode) THEN
742 IF (lbias) THEN
743 CALL get_mc_par(mc_par(remove_box)%mc_par, mc_bias_file=mc_bias_file)
744 CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
745 abc_remove(:), dat_file(remove_box), nchains(:, remove_box), &
746 mc_bias_file)
747 ELSE
748 CALL get_mc_par(mc_par(remove_box)%mc_par, mc_input_file=mc_input_file)
749 CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
750 abc_remove(:), dat_file(remove_box), nchains(:, remove_box), &
751 mc_input_file)
752 END IF
753
754 END IF
755 nchains(molecule_type, remove_box) = nchains(molecule_type, remove_box) + 1
756
757 ELSE
758 DO iatom = 1, start_atom_rem - 1
759 remove_coords(1:3, iatom) = &
760 particles_old(remove_box)%list%els(iatom)%r(1:3)
761 atom_names_remove(iatom) = &
762 particles_old(remove_box)%list%els(iatom)%atomic_kind%name
763 END DO
764 DO iatom = start_atom_rem + nunits(molecule_type), nunits_tot(remove_box)
765 remove_coords(1:3, iatom - nunits(molecule_type)) = &
766 particles_old(remove_box)%list%els(iatom)%r(1:3)
767 atom_names_remove(iatom - nunits(molecule_type)) = &
768 particles_old(remove_box)%list%els(iatom)%atomic_kind%name
769 END DO
770
771! make the .dat file
772 IF (ionode) THEN
773 nchains(molecule_type, remove_box) = nchains(molecule_type, remove_box) - 1
774 IF (lbias) THEN
775 CALL get_mc_par(mc_par(remove_box)%mc_par, mc_bias_file=mc_bias_file)
776 CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
777 abc_remove(:), dat_file(remove_box), nchains(:, remove_box), &
778 mc_bias_file)
779 ELSE
780 CALL get_mc_par(mc_par(remove_box)%mc_par, mc_input_file=mc_input_file)
781 CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
782 abc_remove(:), dat_file(remove_box), nchains(:, remove_box), &
783 mc_input_file)
784 END IF
785 nchains(molecule_type, remove_box) = nchains(molecule_type, remove_box) + 1
786
787 END IF
788 END IF
789
790! deallocate r_insert_mol
791 DEALLOCATE (r_insert_mol)
792
793! now let's create the two new environments with the different number
794! of molecules
795 ALLOCATE (test_env(1:2))
796 CALL mc_create_force_env(test_env(insert_box)%force_env, input_declaration, &
797 para_env, dat_file(insert_box))
798 CALL mc_create_force_env(test_env(remove_box)%force_env, input_declaration, &
799 para_env, dat_file(remove_box))
800
801! allocate an array we'll need
802 ALLOCATE (r_cbmc(1:3, 1:ins_atoms))
803 ALLOCATE (cbmc_energies(1:nswapmoves, 1:2))
804
805 loverlap_ins = .false.
806 loverlap_rem = .false.
807
808! compute the new molecule information...we need this for the CBMC part
809 IF (rem_atoms == 0) THEN
810 CALL mc_determine_molecule_info(test_env, mc_molecule_info_test, &
811 box_number=remove_box)
812 ELSE
813 CALL mc_determine_molecule_info(test_env, mc_molecule_info_test)
814 END IF
815 CALL get_mc_molecule_info(mc_molecule_info_test, nchains=nchains_test, &
816 mol_type=mol_type_test)
817
818! figure out the position of the molecule we're inserting, and the
819! Rosenbluth weight
820 start_mol = 1
821 DO jbox = 1, insert_box - 1
822 start_mol = start_mol + sum(nchains_test(:, jbox))
823 END DO
824 end_mol = start_mol + sum(nchains_test(:, insert_box)) - 1
825
826 IF (lbias) THEN
827 CALL generate_cbmc_swap_config(test_env(insert_box)%force_env, &
828 beta, max_val, min_val, exp_max_val, &
829 exp_min_val, nswapmoves, weight_new, start_atom_ins, ins_atoms, nunits(:), &
830 nunits(molecule_type), mass(:, molecule_type), loverlap_ins, bias_energy_new(insert_box), &
831 bias_energy_old(insert_box), ionode, .false., mol_type_test(start_mol:end_mol), &
832 nchains_test(:, insert_box), source, group, rng_stream)
833
834! the energy that comes out of the above routine is the difference...we want
835! the real energy for the acceptance rule...we don't do this for the
836! lbias=.false. case because it doesn't appear in the acceptance rule, and
837! we compensate in case of acceptance
838 bias_energy_new(insert_box) = bias_energy_new(insert_box) + &
839 bias_energy_old(insert_box)
840 ELSE
841 CALL generate_cbmc_swap_config(test_env(insert_box)%force_env, &
842 beta, max_val, min_val, exp_max_val, &
843 exp_min_val, nswapmoves, weight_new, start_atom_ins, ins_atoms, nunits(:), &
844 nunits(molecule_type), mass(:, molecule_type), loverlap_ins, new_energy(insert_box), &
845 old_energy(insert_box), ionode, .false., mol_type_test(start_mol:end_mol), &
846 nchains_test(:, insert_box), source, group, rng_stream)
847 END IF
848
849 CALL force_env_get(test_env(insert_box)%force_env, &
850 subsys=insert_sys)
851 CALL cp_subsys_get(insert_sys, &
852 particles=particles_insert)
853
854 DO iatom = 1, ins_atoms
855 r_cbmc(1:3, iatom) = particles_insert%els(iatom)%r(1:3)
856 END DO
857
858! make sure there is no overlap
859
860 IF (loverlap_ins .OR. loverlap_rem) THEN
861! deallocate some stuff
862 CALL mc_molecule_info_destroy(mc_molecule_info_test)
863 CALL force_env_release(test_env(insert_box)%force_env)
864 CALL force_env_release(test_env(remove_box)%force_env)
865 DEALLOCATE (insert_coords)
866 DEALLOCATE (remove_coords)
867 DEALLOCATE (r_cbmc)
868 DEALLOCATE (cbmc_energies)
869 DEALLOCATE (oldsys)
870 DEALLOCATE (particles_old)
871 DEALLOCATE (test_env)
872 CALL timestop(handle)
873 RETURN
874 END IF
875
876! broadcast the chosen coordinates to all processors
877
878 CALL force_env_get(test_env(insert_box)%force_env, &
879 subsys=insert_sys)
880 CALL cp_subsys_get(insert_sys, &
881 particles=particles_insert)
882
883 DO iatom = 1, ins_atoms
884 particles_insert%els(iatom)%r(1:3) = &
885 r_cbmc(1:3, iatom)
886 END DO
887
888! if we made it this far, we have no overlaps
889 moves(molecule_type, insert_box)%moves%grown = &
890 moves(molecule_type, insert_box)%moves%grown + 1
891
892! if we're biasing, we need to make environments with the non-biasing
893! potentials, and calculate the energies
894 IF (lbias) THEN
895
896 ALLOCATE (test_env_bias(1:2))
897
898! first, the environment to which we added a molecule
899 CALL get_mc_par(mc_par(insert_box)%mc_par, mc_input_file=mc_input_file)
900 IF (ionode) CALL mc_make_dat_file_new(r_cbmc(:, :), atom_names_insert(:), ins_atoms, &
901 abc_insert(:), dat_file(insert_box), nchains_test(:, insert_box), &
902 mc_input_file)
903 test_env_bias(insert_box)%force_env => test_env(insert_box)%force_env
904 NULLIFY (test_env(insert_box)%force_env)
905 CALL mc_create_force_env(test_env(insert_box)%force_env, input_declaration, &
906 para_env, dat_file(insert_box))
907
908 CALL force_env_calc_energy_force(test_env(insert_box)%force_env, &
909 calc_force=.false.)
910 CALL force_env_get(test_env(insert_box)%force_env, &
911 potential_energy=new_energy(insert_box))
912
913! now the environment that has one less molecule
914 IF (sum(nchains_test(:, remove_box)) == 0) THEN
915 CALL get_mc_par(mc_par(remove_box)%mc_par, mc_input_file=mc_input_file)
916 IF (ionode) CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
917 abc_remove(:), dat_file(remove_box), nchains_test(:, remove_box), &
918 mc_input_file)
919 test_env_bias(remove_box)%force_env => test_env(remove_box)%force_env
920 NULLIFY (test_env(remove_box)%force_env)
921 CALL mc_create_force_env(test_env(remove_box)%force_env, input_declaration, &
922 para_env, dat_file(remove_box))
923 new_energy(remove_box) = 0.0e0_dp
924 bias_energy_new(remove_box) = 0.0e0_dp
925 ELSE
926 CALL get_mc_par(mc_par(remove_box)%mc_par, mc_input_file=mc_input_file)
927 IF (ionode) CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
928 abc_remove(:), dat_file(remove_box), nchains_test(:, remove_box), &
929 mc_input_file)
930 test_env_bias(remove_box)%force_env => test_env(remove_box)%force_env
931 NULLIFY (test_env(remove_box)%force_env)
932 CALL mc_create_force_env(test_env(remove_box)%force_env, input_declaration, &
933 para_env, dat_file(remove_box))
934 CALL force_env_calc_energy_force(test_env(remove_box)%force_env, &
935 calc_force=.false.)
936 CALL force_env_get(test_env(remove_box)%force_env, &
937 potential_energy=new_energy(remove_box))
938 CALL force_env_calc_energy_force(test_env_bias(remove_box)%force_env, &
939 calc_force=.false.)
940 CALL force_env_get(test_env_bias(remove_box)%force_env, &
941 potential_energy=bias_energy_new(remove_box))
942 END IF
943 ELSE
944 IF (sum(nchains_test(:, remove_box)) == 0) THEN
945 new_energy(remove_box) = 0.0e0_dp
946 ELSE
947 CALL force_env_calc_energy_force(test_env(remove_box)%force_env, &
948 calc_force=.false.)
949 CALL force_env_get(test_env(remove_box)%force_env, &
950 potential_energy=new_energy(remove_box))
951 END IF
952 END IF
953
954! now we need to figure out the rosenbluth weight for the old configuration...
955! we wait until now to do that because we need the energy of the box that
956! has had a molecule removed...notice we use the environment that has not
957! had a molecule removed for the CBMC configurations, and therefore nchains
958! and mol_type instead of nchains_test and mol_type_test
959 start_mol = 1
960 DO jbox = 1, remove_box - 1
961 start_mol = start_mol + sum(nchains(:, jbox))
962 END DO
963 end_mol = start_mol + sum(nchains(:, remove_box)) - 1
964 IF (lbias) THEN
965 CALL generate_cbmc_swap_config(bias_env(remove_box)%force_env, &
966 beta, max_val, min_val, exp_max_val, &
967 exp_min_val, nswapmoves, weight_old, start_atom_rem, nunits_tot(remove_box), &
968 nunits(:), nunits(molecule_type), mass(:, molecule_type), loverlap_rem, rdum, &
969 bias_energy_new(remove_box), ionode, .true., mol_type(start_mol:end_mol), &
970 nchains(:, remove_box), source, group, rng_stream)
971 ELSE
972 CALL generate_cbmc_swap_config(force_env(remove_box)%force_env, &
973 beta, max_val, min_val, exp_max_val, &
974 exp_min_val, nswapmoves, weight_old, start_atom_rem, nunits_tot(remove_box), &
975 nunits(:), nunits(molecule_type), mass(:, molecule_type), loverlap_rem, rdum, &
976 new_energy(remove_box), ionode, .true., mol_type(start_mol:end_mol), &
977 nchains(:, remove_box), source, group, rng_stream)
978 END IF
979
980! figure out the prefactor to the boltzmann weight in the acceptance
981! rule, based on numbers of particles and volumes
982
983 prefactor = real(nchains(molecule_type, remove_box), dp)/ &
984 REAL(nchains(molecule_type, insert_box) + 1, dp)* &
985 vol_insert/vol_remove
986
987 IF (lbias) THEN
988
989 del_quickstep_energy = (-beta)*(new_energy(insert_box) - &
990 old_energy(insert_box) + new_energy(remove_box) - &
991 old_energy(remove_box) - (bias_energy_new(insert_box) + &
992 bias_energy_new(remove_box) - bias_energy_old(insert_box) &
993 - bias_energy_old(remove_box)))
994
995 IF (del_quickstep_energy > exp_max_val) THEN
996 del_quickstep_energy = max_val
997 ELSE IF (del_quickstep_energy < exp_min_val) THEN
998 del_quickstep_energy = min_val
999 ELSE
1000 del_quickstep_energy = exp(del_quickstep_energy)
1001 END IF
1002 w = prefactor*del_quickstep_energy*weight_new/weight_old &
1003 *exp(beta*(eta_remove(molecule_type) - eta_insert(molecule_type)))
1004
1005 ELSE
1006 w = prefactor*weight_new/weight_old &
1007 *exp(beta*(eta_remove(molecule_type) - eta_insert(molecule_type)))
1008
1009 END IF
1010
1011! check if the move is accepted
1012 IF (w >= 1.0e0_dp) THEN
1013 rand = 0.0e0_dp
1014 ELSE
1015 IF (ionode) rand = rng_stream%next()
1016 CALL group%bcast(rand, source)
1017 END IF
1018
1019!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1020 IF (rand < w) THEN
1021
1022! accept the move
1023
1024! accept the move
1025 moves(molecule_type, insert_box)%moves%swap%successes = &
1026 moves(molecule_type, insert_box)%moves%swap%successes + 1
1027
1028! we need to compensate for the fact that we take the difference in
1029! generate_cbmc_config to keep the exponetials small
1030 IF (.NOT. lbias) THEN
1031 new_energy(insert_box) = new_energy(insert_box) + &
1032 old_energy(insert_box)
1033 END IF
1034
1035 DO ibox = 1, 2
1036! update energies
1037 energy_check(ibox) = energy_check(ibox) + (new_energy(ibox) - &
1038 old_energy(ibox))
1039 old_energy(ibox) = new_energy(ibox)
1040! if we're biasing the update the biasing energy
1041 IF (lbias) THEN
1042 last_bias_energy(ibox) = bias_energy_new(ibox)
1043 bias_energy_old(ibox) = bias_energy_new(ibox)
1044 END IF
1045
1046 END DO
1047
1048! change particle numbers...basically destroy the old mc_molecule_info and attach
1049! the new stuff to the mc_pars
1050! figure out the molecule information for the new environments
1051 CALL mc_molecule_info_destroy(mc_molecule_info)
1052 CALL set_mc_par(mc_par(insert_box)%mc_par, mc_molecule_info=mc_molecule_info_test)
1053 CALL set_mc_par(mc_par(remove_box)%mc_par, mc_molecule_info=mc_molecule_info_test)
1054
1055! update coordinates
1056 CALL force_env_get(test_env(insert_box)%force_env, &
1057 subsys=insert_sys)
1058 CALL cp_subsys_get(insert_sys, &
1059 particles=particles_insert)
1060 DO ipart = 1, ins_atoms
1061 r_old(1:3, ipart, insert_box) = particles_insert%els(ipart)%r(1:3)
1062 END DO
1063 CALL force_env_get(test_env(remove_box)%force_env, &
1064 subsys=remove_sys)
1065 CALL cp_subsys_get(remove_sys, &
1066 particles=particles_remove)
1067 DO ipart = 1, rem_atoms
1068 r_old(1:3, ipart, remove_box) = particles_remove%els(ipart)%r(1:3)
1069 END DO
1070
1071 ! insertion box
1072 CALL force_env_release(force_env(insert_box)%force_env)
1073 force_env(insert_box)%force_env => test_env(insert_box)%force_env
1074
1075 ! removal box
1076 CALL force_env_release(force_env(remove_box)%force_env)
1077 force_env(remove_box)%force_env => test_env(remove_box)%force_env
1078
1079! if we're biasing, update the bias_env
1080 IF (lbias) THEN
1081 CALL force_env_release(bias_env(insert_box)%force_env)
1082 bias_env(insert_box)%force_env => test_env_bias(insert_box)%force_env
1083 CALL force_env_release(bias_env(remove_box)%force_env)
1084 bias_env(remove_box)%force_env => test_env_bias(remove_box)%force_env
1085 DEALLOCATE (test_env_bias)
1086 END IF
1087
1088 ELSE
1089
1090! reject the move
1091 CALL mc_molecule_info_destroy(mc_molecule_info_test)
1092 CALL force_env_release(test_env(insert_box)%force_env)
1093 CALL force_env_release(test_env(remove_box)%force_env)
1094 IF (lbias) THEN
1095 CALL force_env_release(test_env_bias(insert_box)%force_env)
1096 CALL force_env_release(test_env_bias(remove_box)%force_env)
1097 DEALLOCATE (test_env_bias)
1098 END IF
1099 END IF
1100
1101! deallocate some stuff
1102 DEALLOCATE (insert_coords)
1103 DEALLOCATE (remove_coords)
1104 DEALLOCATE (test_env)
1105 DEALLOCATE (cbmc_energies)
1106 DEALLOCATE (r_cbmc)
1107 DEALLOCATE (oldsys)
1108 DEALLOCATE (particles_old)
1109
1110! end the timing
1111 CALL timestop(handle)
1112
1113 END SUBROUTINE mc_ge_swap_move
1114
1115! **************************************************************************************************
1116!> \brief performs a Monte Carlo move that alters the volume of the simulation boxes,
1117!> keeping the total volume of the two boxes the same
1118!> \param mc_par the mc parameters for the force env
1119!> \param force_env the force environments used in the move
1120!> \param moves the structure that keeps track of how many moves have been
1121!> accepted/rejected
1122!> \param move_updates the structure that keeps track of how many moves have
1123!> been accepted/rejected since the last time the displacements
1124!> were updated
1125!> \param nnstep the total number of Monte Carlo moves already performed
1126!> \param old_energy the energy of the last accepted move involving an
1127!> unbiased potential calculation
1128!> \param energy_check the running total of how much the energy has changed
1129!> since the initial configuration
1130!> \param r_old the coordinates of the last accepted move involving a
1131!> Quickstep calculation
1132!> \param rng_stream the stream we pull random numbers from
1133!> \author MJM
1134! **************************************************************************************************
1135 SUBROUTINE mc_ge_volume_move(mc_par, force_env, moves, move_updates, &
1136 nnstep, old_energy, energy_check, r_old, rng_stream)
1137
1138 TYPE(mc_simulation_parameters_p_type), &
1139 DIMENSION(:), POINTER :: mc_par
1140 TYPE(force_env_p_type), DIMENSION(:), POINTER :: force_env
1141 TYPE(mc_moves_p_type), DIMENSION(:, :), POINTER :: moves, move_updates
1142 INTEGER, INTENT(IN) :: nnstep
1143 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: old_energy, energy_check
1144 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: r_old
1145 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1146
1147 CHARACTER(len=*), PARAMETER :: routinen = 'mc_ge_volume_move'
1148
1149 CHARACTER(LEN=200) :: fft_lib
1150 CHARACTER(LEN=40), DIMENSION(1:2) :: dat_file
1151 INTEGER :: cl, end_atom, end_mol, handle, iatom, ibox, imolecule, iside, j, jatom, jbox, &
1152 max_atoms, molecule_index, molecule_type, print_level, source, start_atom, start_mol
1153 INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
1154 INTEGER, DIMENSION(:, :), POINTER :: nchains
1155 LOGICAL :: ionode
1156 LOGICAL, ALLOCATABLE, DIMENSION(:) :: loverlap
1157 LOGICAL, DIMENSION(1:2) :: lempty
1158 REAL(dp), DIMENSION(:, :), POINTER :: mass
1159 REAL(kind=dp) :: beta, prefactor, rand, rmvolume, &
1160 vol_dis, w
1161 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r
1162 REAL(kind=dp), DIMENSION(1:2) :: new_energy, volume_new, volume_old
1163 REAL(kind=dp), DIMENSION(1:3) :: center_of_mass, center_of_mass_new, diff
1164 REAL(kind=dp), DIMENSION(1:3, 1:2) :: abc, new_cell_length, old_cell_length
1165 REAL(kind=dp), DIMENSION(1:3, 1:3, 1:2) :: hmat_test
1166 TYPE(cell_p_type), DIMENSION(:), POINTER :: cell, cell_old, cell_test
1167 TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: oldsys
1168 TYPE(cp_subsys_type), POINTER :: subsys
1169 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1170 TYPE(mp_comm_type) :: group
1171 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles_old
1172
1173! begin the timing of the subroutine
1174
1175 CALL timeset(routinen, handle)
1176
1177! nullify some pointers
1178 NULLIFY (particles_old, cell, oldsys, cell_old, cell_test, subsys)
1179
1180! get some data from mc_par
1181 CALL get_mc_par(mc_par(1)%mc_par, ionode=ionode, source=source, &
1182 group=group, dat_file=dat_file(1), rmvolume=rmvolume, &
1183 beta=beta, cl=cl, fft_lib=fft_lib, &
1184 mc_molecule_info=mc_molecule_info)
1185 CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot, &
1186 mass=mass, nchains=nchains, nunits=nunits, mol_type=mol_type)
1187
1188 print_level = 1
1189 CALL get_mc_par(mc_par(2)%mc_par, dat_file=dat_file(2))
1190
1191! allocate some stuff
1192 max_atoms = max(nunits_tot(1), nunits_tot(2))
1193 ALLOCATE (r(1:3, max_atoms, 1:2))
1194 ALLOCATE (oldsys(1:2))
1195 ALLOCATE (particles_old(1:2))
1196 ALLOCATE (cell(1:2))
1197 ALLOCATE (cell_test(1:2))
1198 ALLOCATE (cell_old(1:2))
1199 ALLOCATE (loverlap(1:2))
1200
1201! check for empty boxes...need to be careful because we can't build
1202! a force_env with no particles
1203 DO ibox = 1, 2
1204 lempty(ibox) = .false.
1205 IF (sum(nchains(:, ibox)) == 0) THEN
1206 lempty(ibox) = .true.
1207 END IF
1208 END DO
1209
1210! record the attempt
1211 DO ibox = 1, 2
1212 moves(1, ibox)%moves%volume%attempts = &
1213 moves(1, ibox)%moves%volume%attempts + 1
1214 move_updates(1, ibox)%moves%volume%attempts = &
1215 move_updates(1, ibox)%moves%volume%attempts + 1
1216 END DO
1217
1218! now let's grab the cell length and particle positions
1219 DO ibox = 1, 2
1220 CALL force_env_get(force_env(ibox)%force_env, &
1221 subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
1222 CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
1223 NULLIFY (cell_old(ibox)%cell)
1224 CALL cell_create(cell_old(ibox)%cell)
1225 CALL cell_clone(cell(ibox)%cell, cell_old(ibox)%cell, tag="CELL_OLD")
1226 CALL cp_subsys_get(oldsys(ibox)%subsys, &
1227 particles=particles_old(ibox)%list)
1228
1229! find the old cell length
1230 old_cell_length(1:3, ibox) = abc(1:3, ibox)
1231
1232 END DO
1233
1234 DO ibox = 1, 2
1235
1236! save the old coordinates
1237 DO iatom = 1, nunits_tot(ibox)
1238 r(1:3, iatom, ibox) = particles_old(ibox)%list%els(iatom)%r(1:3)
1239 END DO
1240
1241 END DO
1242
1243! call a random number to figure out how far we're moving
1244 IF (ionode) rand = rng_stream%next()
1245 CALL group%bcast(rand, source)
1246
1247 vol_dis = rmvolume*(rand - 0.5e0_dp)*2.0e0_dp
1248
1249! add to one box, subtract from the other
1250 IF (old_cell_length(1, 1)*old_cell_length(2, 1)* &
1251 old_cell_length(3, 1) + vol_dis <= (3.0e0_dp/angstrom)**3) THEN
1252 cpabort('GE_volume moves are trying to make box 1 smaller than 3')
1253 END IF
1254 IF (old_cell_length(1, 2)*old_cell_length(2, 2)* &
1255 old_cell_length(3, 2) + vol_dis <= (3.0e0_dp/angstrom)**3) THEN
1256 cpabort('GE_volume moves are trying to make box 2 smaller than 3')
1257 END IF
1258
1259 DO iside = 1, 3
1260 new_cell_length(iside, 1) = (old_cell_length(1, 1)**3 + &
1261 vol_dis)**(1.0e0_dp/3.0e0_dp)
1262 new_cell_length(iside, 2) = (old_cell_length(1, 2)**3 - &
1263 vol_dis)**(1.0e0_dp/3.0e0_dp)
1264 END DO
1265
1266! now we need to make the new cells
1267 DO ibox = 1, 2
1268 hmat_test(:, :, ibox) = 0.0e0_dp
1269 hmat_test(1, 1, ibox) = new_cell_length(1, ibox)
1270 hmat_test(2, 2, ibox) = new_cell_length(2, ibox)
1271 hmat_test(3, 3, ibox) = new_cell_length(3, ibox)
1272 NULLIFY (cell_test(ibox)%cell)
1273 CALL cell_create(cell_test(ibox)%cell, hmat=hmat_test(:, :, ibox), &
1274 periodic=cell(ibox)%cell%perd)
1275 CALL force_env_get(force_env(ibox)%force_env, subsys=subsys)
1276 CALL cp_subsys_set(subsys, cell=cell_test(ibox)%cell)
1277 END DO
1278
1279 DO ibox = 1, 2
1280
1281! save the coords
1282 DO iatom = 1, nunits_tot(ibox)
1283 r(1:3, iatom, ibox) = particles_old(ibox)%list%els(iatom)%r(1:3)
1284 END DO
1285
1286! now we need to scale the coordinates of all the molecules by the
1287! center of mass
1288 start_atom = 1
1289 molecule_index = 1
1290 DO jbox = 1, ibox - 1
1291 IF (jbox == ibox) EXIT
1292 molecule_index = molecule_index + sum(nchains(:, jbox))
1293 END DO
1294 DO imolecule = 1, sum(nchains(:, ibox))
1295 molecule_type = mol_type(imolecule + molecule_index - 1)
1296 IF (imolecule /= 1) THEN
1297 start_atom = start_atom + nunits(mol_type(imolecule + molecule_index - 2))
1298 END IF
1299 end_atom = start_atom + nunits(molecule_type) - 1
1300
1301! now find the center of mass
1302 CALL get_center_of_mass(r(:, start_atom:end_atom, ibox), &
1303 nunits(molecule_type), center_of_mass(:), mass(:, molecule_type))
1304
1305! scale the center of mass and determine the vector that points from the
1306! old COM to the new one
1307 center_of_mass_new(1:3) = center_of_mass(1:3)* &
1308 new_cell_length(1:3, ibox)/old_cell_length(1:3, ibox)
1309 DO j = 1, 3
1310 diff(j) = center_of_mass_new(j) - center_of_mass(j)
1311! now change the particle positions
1312 DO jatom = start_atom, end_atom
1313 particles_old(ibox)%list%els(jatom)%r(j) = &
1314 particles_old(ibox)%list%els(jatom)%r(j) + diff(j)
1315 END DO
1316
1317 END DO
1318 END DO
1319
1320! check for any overlaps we might have
1321 start_mol = 1
1322 DO jbox = 1, ibox - 1
1323 start_mol = start_mol + sum(nchains(:, jbox))
1324 END DO
1325 end_mol = start_mol + sum(nchains(:, ibox)) - 1
1326 CALL check_for_overlap(force_env(ibox)%force_env, &
1327 nchains(:, ibox), nunits, loverlap(ibox), mol_type(start_mol:end_mol), &
1328 cell_length=new_cell_length(:, ibox))
1329
1330 END DO
1331
1332! determine the overall energy difference
1333
1334 DO ibox = 1, 2
1335 IF (loverlap(ibox)) cycle
1336! remake the force environment and calculate the energy
1337 IF (lempty(ibox)) THEN
1338 new_energy(ibox) = 0.0e0_dp
1339 ELSE
1340
1341 CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
1342 calc_force=.false.)
1343 CALL force_env_get(force_env(ibox)%force_env, &
1344 potential_energy=new_energy(ibox))
1345
1346 END IF
1347 END DO
1348
1349! accept or reject the move
1350 DO ibox = 1, 2
1351 volume_new(ibox) = new_cell_length(1, ibox)* &
1352 new_cell_length(2, ibox)*new_cell_length(3, ibox)
1353 volume_old(ibox) = old_cell_length(1, ibox)* &
1354 old_cell_length(2, ibox)*old_cell_length(3, ibox)
1355 END DO
1356 prefactor = (volume_new(1)/volume_old(1))**(sum(nchains(:, 1)))* &
1357 (volume_new(2)/volume_old(2))**(sum(nchains(:, 2)))
1358
1359 IF (loverlap(1) .OR. loverlap(2)) THEN
1360 w = 0.0e0_dp
1361 ELSE
1362 w = prefactor*exp(-beta* &
1363 (new_energy(1) + new_energy(2) - &
1364 old_energy(1) - old_energy(2)))
1365
1366 END IF
1367
1368 IF (w >= 1.0e0_dp) THEN
1369 w = 1.0e0_dp
1370 rand = 0.0e0_dp
1371 ELSE
1372 IF (ionode) rand = rng_stream%next()
1373 CALL group%bcast(rand, source)
1374 END IF
1375
1376 IF (rand < w) THEN
1377
1378! write cell length, volume, density, and trial displacement to a file
1379 IF (ionode) THEN
1380
1381 WRITE (cl, *) nnstep, new_cell_length(1, 1)* &
1382 angstrom, vol_dis*(angstrom)**3, new_cell_length(1, 2)* &
1383 angstrom
1384 WRITE (cl, *) nnstep, new_energy(1), &
1385 old_energy(1), new_energy(2), old_energy(2)
1386 WRITE (cl, *) prefactor, w
1387 END IF
1388
1389 DO ibox = 1, 2
1390! accept the move
1391 moves(1, ibox)%moves%volume%successes = &
1392 moves(1, ibox)%moves%volume%successes + 1
1393 move_updates(1, ibox)%moves%volume%successes = &
1394 move_updates(1, ibox)%moves%volume%successes + 1
1395
1396! update energies
1397 energy_check(ibox) = energy_check(ibox) + (new_energy(ibox) - &
1398 old_energy(ibox))
1399 old_energy(ibox) = new_energy(ibox)
1400
1401! and the new "old" coordinates
1402 DO iatom = 1, nunits_tot(ibox)
1403 r_old(1:3, iatom, ibox) = &
1404 particles_old(ibox)%list%els(iatom)%r(1:3)
1405 END DO
1406
1407 END DO
1408 ELSE
1409
1410! reject the move
1411! write cell length, volume, density, and trial displacement to a file
1412 IF (ionode) THEN
1413
1414 WRITE (cl, *) nnstep, new_cell_length(1, 1)* &
1415 angstrom, vol_dis*(angstrom)**3, new_cell_length(1, 2)* &
1416 angstrom
1417 WRITE (cl, *) nnstep, new_energy(1), &
1418 old_energy(1), new_energy(2), old_energy(2)
1419 WRITE (cl, *) prefactor, w
1420
1421 END IF
1422
1423! reset the cell and particle positions
1424 DO ibox = 1, 2
1425 CALL force_env_get(force_env(ibox)%force_env, subsys=subsys)
1426 CALL cp_subsys_set(subsys, cell=cell_old(ibox)%cell)
1427 DO iatom = 1, nunits_tot(ibox)
1428 particles_old(ibox)%list%els(iatom)%r(1:3) = r_old(1:3, iatom, ibox)
1429 END DO
1430 END DO
1431
1432 END IF
1433
1434! free up some memory
1435 DO ibox = 1, 2
1436 CALL cell_release(cell_test(ibox)%cell)
1437 CALL cell_release(cell_old(ibox)%cell)
1438 END DO
1439 DEALLOCATE (r)
1440 DEALLOCATE (oldsys)
1441 DEALLOCATE (particles_old)
1442 DEALLOCATE (cell)
1443 DEALLOCATE (cell_old)
1444 DEALLOCATE (cell_test)
1445 DEALLOCATE (loverlap)
1446
1447! end the timing
1448 CALL timestop(handle)
1449
1450 END SUBROUTINE mc_ge_volume_move
1451
1452END MODULE mc_ge_moves
1453
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:94
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:608
cell_types::cell_clone
subroutine, public cell_clone(cell_in, cell_out, tag)
Clone cell variable.
Definition cell_types.F:118
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:210
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_set
subroutine, public cp_subsys_set(subsys, atomic_kinds, particles, local_particles, molecules, molecule_kinds, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, results, cell, cell_ref, use_ref_cell)
sets various propreties of the subsys
Definition cp_subsys_types.F:210
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
fft_lib
Definition fft_lib.F:7
force_env_methods
Interface for the force calculations.
Definition force_env_methods.F:17
force_env_methods::force_env_calc_energy_force
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
Definition force_env_methods.F:218
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:344
force_env_types::force_env_release
recursive subroutine, public force_env_release(force_env)
releases the given force env
Definition force_env_types.F:207
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::dump_xmol
integer, parameter, public dump_xmol
Definition input_constants.F:386
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
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
mc_control
contains some general routines for dealing with the restart files and creating force_env for MC use
Definition mc_control.F:15
mc_control::mc_create_force_env
subroutine, public mc_create_force_env(force_env, input_declaration, para_env, input_file_name, globenv_new)
creates a force environment for any of the different kinds of MC simulations we can do (FIST,...
Definition mc_control.F:351
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_ge_moves
contains the Monte Carlo moves that can handle more than one box, including the Quickstep move,...
Definition mc_ge_moves.F:18
mc_ge_moves::mc_ge_volume_move
subroutine, public mc_ge_volume_move(mc_par, force_env, moves, move_updates, nnstep, old_energy, energy_check, r_old, rng_stream)
performs a Monte Carlo move that alters the volume of the simulation boxes, keeping the total volume ...
Definition mc_ge_moves.F:1137
mc_ge_moves::mc_quickstep_move
subroutine, public mc_quickstep_move(mc_par, force_env, bias_env, moves, lreject, move_updates, energy_check, r_old, nnstep, old_energy, bias_energy_new, last_bias_energy, nboxes, box_flag, subsys, particles, rng_stream, unit_conv)
computes the acceptance of a series of biased or unbiased moves (translation, rotation,...
Definition mc_ge_moves.F:105
mc_ge_moves::mc_ge_swap_move
subroutine, public mc_ge_swap_move(mc_par, force_env, bias_env, moves, energy_check, r_old, old_energy, input_declaration, para_env, bias_energy_old, last_bias_energy, rng_stream)
attempts a swap move between two simulation boxes
Definition mc_ge_moves.F:431
mc_misc
contains miscellaneous subroutines used in the Monte Carlo runs, mostly I/O stuff
Definition mc_misc.F:13
mc_misc::mc_make_dat_file_new
subroutine, public mc_make_dat_file_new(coordinates, atom_names, nunits_tot, box_length, filename, nchains, mc_input_file)
writes a new input file that CP2K can read in for when we want to change a force env (change molecule...
Definition mc_misc.F:461
mc_move_control
control the handling of the move data in Monte Carlo (MC) simulations
Definition mc_move_control.F:14
mc_move_control::q_move_accept
subroutine, public q_move_accept(moves, lbias)
updates accepted moves in the given structure...assumes you've been recording all successful moves in...
Definition mc_move_control.F:200
mc_move_control::move_q_reinit
subroutine, public move_q_reinit(moves, lbias)
sets all qsuccess counters back to zero
Definition mc_move_control.F:150
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
mc_types::mc_determine_molecule_info
subroutine, public mc_determine_molecule_info(force_env, mc_molecule_info, box_number, coordinates_empty)
figures out the number of total molecules, the number of atoms in each molecule, an array with the mo...
Definition mc_types.F:1788
mc_types::mc_molecule_info_destroy
subroutine, public mc_molecule_info_destroy(mc_molecule_info)
deallocates all the arrays in the mc_molecule_info_type
Definition mc_types.F:2003
mc_types::set_mc_par
subroutine, public set_mc_par(mc_par, rm, cl, diff, nstart, rmvolume, rmcltrans, rmbond, rmangle, rmdihedral, rmrot, rmtrans, program, nmoves, nswapmoves, lstop, temperature, pressure, rclus, iuptrans, iupcltrans, iupvolume, pmswap, pmvolume, pmtraion, pmtrans, pmcltrans, beta, rcut, iprint, lbias, nstep, lrestart, ldiscrete, discrete_step, pmavbmc, mc_molecule_info, pmavbmc_mol, pmtrans_mol, pmrot_mol, pmtraion_mol, pmswap_mol, avbmc_rmin, avbmc_rmax, avbmc_atom, pbias, ensemble, pmvol_box, pmclus_box, eta, mc_input_file, mc_bias_file, exp_max_val, exp_min_val, min_val, max_val, pmhmc, pmhmc_box, lhmc, ionode, source, group, rand2skip)
changes the private elements of the mc_parameters_type
Definition mc_types.F:667
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
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
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::write_particle_coordinates
subroutine, public write_particle_coordinates(particle_set, iunit, output_format, content, title, cell, array, unit_conv, charge_occup, charge_beta, charge_extended, print_kind)
Should be able to write a few formats e.g. xmol, and some binary format (dcd) some format can be used...
Definition particle_methods.F:226
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
cell_types::cell_p_type
Definition cell_types.F:75
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_subsys_types::cp_subsys_p_type
represent a pointer to a subsys, to be able to create arrays of pointers
Definition cp_subsys_types.F:123
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_p_type
allows for the creation of an array of force_env
Definition force_env_types.F:174
input_section_types::section_type
represent a section of the input file
Definition input_section_types.F:102
mc_types::mc_input_file_type
Definition mc_types.F:155
mc_types::mc_molecule_info_type
Definition mc_types.F:169
mc_types::mc_moves_p_type
Definition mc_types.F:239
mc_types::mc_simulation_parameters_p_type
Definition mc_types.F:184
message_passing::mp_comm_type
Definition message_passing.F:158
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:746
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157
particle_list_types::particle_list_p_type
represents a pointer to a list
Definition particle_list_types.F:59
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46