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