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