Line data Source code
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 : ! **************************************************************************************************
18 : MODULE 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 :
68 : CONTAINS
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 368 : SUBROUTINE mc_Quickstep_move(mc_par, force_env, bias_env, moves, &
101 368 : lreject, move_updates, energy_check, r_old, &
102 368 : nnstep, old_energy, bias_energy_new, last_bias_energy, &
103 368 : 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 368 : INTEGER, DIMENSION(:, :), POINTER :: nchains
130 368 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
131 736 : INTEGER, DIMENSION(1:nboxes) :: diff
132 : LOGICAL :: ionode, lbias, loverlap
133 : REAL(KIND=dp) :: BETA, energies, rand, w
134 736 : REAL(KIND=dp), DIMENSION(1:nboxes) :: bias_energy_old, new_energy
135 368 : 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 368 : TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles_bias
139 :
140 : ! begin the timing of the subroutine
141 :
142 368 : CALL timeset(routineN, handle)
143 :
144 368 : 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 368 : mc_molecule_info=mc_molecule_info)
151 : CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types, &
152 368 : nchains=nchains, nunits_tot=nunits_tot, nunits=nunits, mol_type=mol_type)
153 :
154 368 : IF (nboxes .GT. 1) THEN
155 144 : DO ibox = 2, nboxes
156 144 : CALL get_mc_par(mc_par(ibox)%mc_par, diff=diff(ibox))
157 : END DO
158 : END IF
159 :
160 : ! allocate some stuff
161 1912 : ALLOCATE (subsys_bias(1:nboxes))
162 1544 : 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 368 : moves(1, 1)%moves%Quickstep%attempts + 1
168 :
169 : ! grab the coordinates for the force_env
170 808 : DO ibox = 1, nboxes
171 : CALL force_env_get(force_env(ibox)%force_env, &
172 440 : subsys=subsys(ibox)%subsys)
173 : CALL cp_subsys_get(subsys(ibox)%subsys, &
174 808 : 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 808 : DO ibox = 1, nboxes
180 808 : IF (box_flag(ibox) == 1) THEN
181 368 : 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 98 : subsys=subsys_bias(ibox)%subsys)
185 : CALL cp_subsys_get(subsys_bias(ibox)%subsys, &
186 98 : particles=particles_bias(ibox)%list)
187 :
188 2606 : DO iparticle = 1, nunits_tot(ibox)
189 : particles(ibox)%list%els(iparticle)%r(1:3) = &
190 17654 : 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 98 : calc_force=.FALSE.)
195 : CALL force_env_get(force_env(ibox)%force_env, &
196 98 : potential_energy=new_energy(ibox))
197 : ELSE
198 270 : IF (.NOT. lreject) THEN
199 : CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
200 270 : calc_force=.FALSE.)
201 : CALL force_env_get(force_env(ibox)%force_env, &
202 270 : potential_energy=new_energy(ibox))
203 : END IF
204 : END IF
205 : ELSE
206 72 : 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 368 : IF (ionode) THEN
213 :
214 : ! write them out in case something bad happens
215 184 : IF (MOD(nnstep, iprint) == 0) THEN
216 26 : DO ibox = 1, nboxes
217 49 : IF (SUM(nchains(:, ibox)) == 0) THEN
218 0 : WRITE (diff(ibox), *) nnstep
219 0 : WRITE (diff(ibox), *) nchains(:, ibox)
220 : ELSE
221 14 : WRITE (diff(ibox), *) nnstep
222 : CALL write_particle_coordinates( &
223 : particles(ibox)%list%els, &
224 : diff(ibox), dump_xmol, 'POS', 'TRIAL', &
225 14 : unit_conv=unit_conv)
226 : END IF
227 : END DO
228 : END IF
229 : END IF
230 :
231 368 : IF (.NOT. lreject) THEN
232 368 : IF (lbias) THEN
233 :
234 268 : DO ibox = 1, nboxes
235 : ! look for overlap
236 510 : IF (SUM(nchains(:, ibox)) .NE. 0) THEN
237 : ! find the molecule bounds
238 : start_mol = 1
239 242 : DO jbox = 1, ibox - 1
240 386 : start_mol = start_mol + SUM(nchains(:, jbox))
241 : END DO
242 510 : end_mol = start_mol + SUM(nchains(:, ibox)) - 1
243 : CALL check_for_overlap(bias_env(ibox)%force_env, &
244 170 : nchains(:, ibox), nunits(:), loverlap, mol_type(start_mol:end_mol))
245 170 : IF (loverlap) &
246 0 : CPABORT('Quickstep move found an overlap in the old config')
247 : END IF
248 268 : bias_energy_old(ibox) = last_bias_energy(ibox)
249 : END DO
250 :
251 : energies = -BETA*((SUM(new_energy(:)) - SUM(bias_energy_new(:))) &
252 778 : - (SUM(old_energy(:)) - SUM(bias_energy_old(:))))
253 :
254 : ! used to prevent over and underflows
255 98 : IF (energies .GE. -1.0E-8) THEN
256 : w = 1.0_dp
257 34 : ELSEIF (energies .LE. -500.0_dp) THEN
258 : w = 0.0_dp
259 : ELSE
260 34 : w = EXP(energies)
261 : END IF
262 :
263 98 : IF (ionode) THEN
264 134 : DO ibox = 1, nboxes
265 85 : WRITE (diff(ibox), *) nnstep, new_energy(ibox) - &
266 85 : old_energy(ibox), &
267 219 : bias_energy_new(ibox) - bias_energy_old(ibox)
268 : END DO
269 : END IF
270 : ELSE
271 810 : energies = -BETA*(SUM(new_energy(:)) - SUM(old_energy(:)))
272 : ! used to prevent over and underflows
273 270 : IF (energies .GE. 0.0_dp) THEN
274 : w = 1.0_dp
275 156 : ELSEIF (energies .LE. -500.0_dp) THEN
276 : w = 0.0_dp
277 : ELSE
278 156 : w = EXP(energies)
279 : END IF
280 : END IF
281 : ELSE
282 : w = 0.0E0_dp
283 : END IF
284 368 : IF (w .GE. 1.0E0_dp) THEN
285 178 : w = 1.0E0_dp
286 178 : rand = 0.0E0_dp
287 : ELSE
288 190 : IF (ionode) rand = rng_stream%next()
289 190 : CALL group%bcast(rand, source)
290 : END IF
291 :
292 368 : IF (rand .LT. w) THEN
293 :
294 : ! accept the move
295 : moves(1, 1)%moves%Quickstep%successes = &
296 286 : moves(1, 1)%moves%Quickstep%successes + 1
297 :
298 644 : DO ibox = 1, nboxes
299 : ! remember what kind of move we did for lbias=.false.
300 358 : IF (.NOT. lbias) THEN
301 554 : DO itype = 1, nmol_types
302 366 : CALL q_move_accept(moves(itype, ibox)%moves, .TRUE.)
303 366 : CALL q_move_accept(move_updates(itype, ibox)%moves, .TRUE.)
304 :
305 : ! reset the counters
306 366 : CALL move_q_reinit(moves(itype, ibox)%moves, .TRUE.)
307 554 : CALL move_q_reinit(move_updates(itype, ibox)%moves, .TRUE.)
308 : END DO
309 : END IF
310 :
311 1064 : DO itype = 1, nmol_types
312 : ! we need to record all accepted moves since last Quickstep calculation
313 706 : CALL q_move_accept(moves(itype, ibox)%moves, .FALSE.)
314 706 : CALL q_move_accept(move_updates(itype, ibox)%moves, .FALSE.)
315 :
316 : ! reset the counters
317 706 : CALL move_q_reinit(moves(itype, ibox)%moves, .FALSE.)
318 1064 : 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 358 : (new_energy(ibox) - old_energy(ibox))
324 644 : old_energy(ibox) = new_energy(ibox)
325 :
326 : END DO
327 :
328 286 : IF (lbias) THEN
329 268 : DO ibox = 1, nboxes
330 268 : last_bias_energy(ibox) = bias_energy_new(ibox)
331 : END DO
332 : END IF
333 :
334 : ! update coordinates
335 644 : DO ibox = 1, nboxes
336 644 : IF (nunits_tot(ibox) .NE. 0) THEN
337 9566 : DO iparticle = 1, nunits_tot(ibox)
338 : r_old(1:3, iparticle, ibox) = &
339 37190 : 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 164 : DO ibox = 1, nboxes
348 328 : DO itype = 1, nmol_types
349 164 : CALL move_q_reinit(moves(itype, ibox)%moves, .FALSE.)
350 164 : CALL move_q_reinit(move_updates(itype, ibox)%moves, .FALSE.)
351 246 : IF (.NOT. lbias) THEN
352 : ! reset the counters
353 164 : CALL move_q_reinit(moves(itype, ibox)%moves, .TRUE.)
354 164 : CALL move_q_reinit(move_updates(itype, ibox)%moves, .TRUE.)
355 : END IF
356 : END DO
357 :
358 : END DO
359 :
360 4551 : 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 17958 : CALL group%bcast(r_old, source)
365 :
366 164 : DO ibox = 1, nboxes
367 2378 : DO iparticle = 1, nunits_tot(ibox)
368 : particles(ibox)%list%els(iparticle)%r(1:3) = &
369 8856 : r_old(1:3, iparticle, ibox)
370 2214 : IF (lbias .AND. box_flag(ibox) == 1) &
371 : particles_bias(ibox)%list%els(iparticle)%r(1:3) = &
372 82 : r_old(1:3, iparticle, ibox)
373 : END DO
374 : END DO
375 :
376 : ! need to reset the energies of the biasing potential
377 82 : IF (lbias) THEN
378 0 : DO ibox = 1, nboxes
379 0 : 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 808 : DO ibox = 1, nboxes
387 : CALL cp_subsys_set(subsys(ibox)%subsys, &
388 440 : particles=particles(ibox)%list)
389 440 : IF (lbias .AND. box_flag(ibox) == 1) &
390 : CALL cp_subsys_set(subsys_bias(ibox)%subsys, &
391 466 : particles=particles_bias(ibox)%list)
392 : END DO
393 :
394 : ! deallocate some stuff
395 368 : DEALLOCATE (subsys_bias)
396 368 : DEALLOCATE (particles_bias)
397 :
398 : ! end the timing
399 368 : CALL timestop(handle)
400 :
401 368 : 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 22 : SUBROUTINE mc_ge_swap_move(mc_par, force_env, bias_env, moves, &
425 22 : 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 22 : DIMENSION(:) :: atom_names_insert, atom_names_remove
445 : CHARACTER(default_string_length), &
446 22 : 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 22 : INTEGER, DIMENSION(:), POINTER :: mol_type, mol_type_test, nunits, &
453 22 : nunits_tot
454 22 : INTEGER, DIMENSION(:, :), POINTER :: nchains, nchains_test
455 : LOGICAL :: ionode, lbias, loverlap, loverlap_ins, &
456 : loverlap_rem
457 22 : REAL(dp), DIMENSION(:), POINTER :: eta_insert, eta_remove, pmswap_mol
458 22 : 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 22 : 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 22 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mass
466 : TYPE(cell_type), POINTER :: cell_insert, cell_remove
467 22 : TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: oldsys
468 : TYPE(cp_subsys_type), POINTER :: insert_sys, remove_sys
469 22 : 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 22 : 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 22 : CALL timeset(routineN, handle)
479 :
480 : ! reset the overlap flag
481 22 : loverlap = .FALSE.
482 :
483 : ! nullify some pointers
484 22 : NULLIFY (particles_old, mol_type, mol_type_test, mc_input_file, mc_bias_file)
485 22 : NULLIFY (oldsys, atom_names, pmswap_mol, insert_coords, remove_coords)
486 22 : 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 22 : 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 22 : atom_names=atom_names, mass=mass, mol_type=mol_type)
497 :
498 22 : print_level = 1
499 :
500 22 : CALL get_mc_par(mc_par(2)%mc_par, dat_file=dat_file(2))
501 :
502 : ! allocate some stuff
503 66 : ALLOCATE (oldsys(1:2))
504 66 : ALLOCATE (particles_old(1:2))
505 :
506 : ! get the old coordinates
507 66 : DO ibox = 1, 2
508 : CALL force_env_get(force_env(ibox)%force_env, &
509 44 : subsys=oldsys(ibox)%subsys)
510 : CALL cp_subsys_get(oldsys(ibox)%subsys, &
511 66 : particles=particles_old(ibox)%list)
512 : END DO
513 :
514 : ! choose a direction to swap
515 22 : IF (ionode) rand = rng_stream%next()
516 22 : CALL group%bcast(rand, source)
517 :
518 22 : IF (rand .LE. 0.50E0_dp) THEN
519 12 : remove_box = 1
520 12 : insert_box = 2
521 : ELSE
522 10 : remove_box = 2
523 10 : insert_box = 1
524 : END IF
525 :
526 : ! now assign the eta values for the insert and remove boxes
527 22 : CALL get_mc_par(mc_par(remove_box)%mc_par, eta=eta_remove)
528 22 : 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 22 : IF (ionode) rand = rng_stream%next()
536 22 : CALL group%bcast(rand, source)
537 34 : DO itype = 1, nmol_types
538 34 : 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 22 : 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 22 : IF (nchains(molecule_type, remove_box) == 0) THEN
551 0 : loverlap = .TRUE.
552 : moves(molecule_type, insert_box)%moves%empty = &
553 0 : moves(molecule_type, insert_box)%moves%empty + 1
554 : ELSE
555 :
556 22 : IF (ionode) rand = rng_stream%next()
557 22 : CALL group%bcast(rand, source)
558 22 : imolecule = CEILING(rand*nchains(molecule_type, remove_box))
559 : ! figure out the atom number this molecule starts on
560 22 : start_atom_rem = 1
561 34 : DO itype = 1, nmol_types
562 34 : IF (itype == molecule_type) THEN
563 22 : start_atom_rem = start_atom_rem + (imolecule - 1)*nunits(itype)
564 22 : EXIT
565 : ELSE
566 12 : start_atom_rem = start_atom_rem + nchains(itype, remove_box)*nunits(itype)
567 : END IF
568 : END DO
569 :
570 : ! check for overlap
571 22 : start_mol = 1
572 32 : DO jbox = 1, remove_box - 1
573 52 : start_mol = start_mol + SUM(nchains(:, jbox))
574 : END DO
575 66 : end_mol = start_mol + SUM(nchains(:, remove_box)) - 1
576 : CALL check_for_overlap(force_env(remove_box)%force_env, &
577 22 : nchains(:, remove_box), nunits, loverlap, mol_type(start_mol:end_mol))
578 22 : IF (loverlap) CALL cp_abort(__LOCATION__, &
579 0 : 'CBMC swap move found an overlap in the old remove config')
580 22 : start_mol = 1
581 34 : DO jbox = 1, insert_box - 1
582 58 : start_mol = start_mol + SUM(nchains(:, jbox))
583 : END DO
584 66 : end_mol = start_mol + SUM(nchains(:, insert_box)) - 1
585 : CALL check_for_overlap(force_env(insert_box)%force_env, &
586 22 : nchains(:, insert_box), nunits, loverlap, mol_type(start_mol:end_mol))
587 22 : IF (loverlap) CALL cp_abort(__LOCATION__, &
588 0 : 'CBMC swap move found an overlap in the old insert config')
589 : END IF
590 :
591 22 : IF (loverlap) THEN
592 0 : DEALLOCATE (oldsys)
593 0 : DEALLOCATE (particles_old)
594 0 : CALL timestop(handle)
595 0 : RETURN
596 : END IF
597 :
598 : ! figure out how many atoms will be in each box after the move
599 22 : ins_atoms = nunits_tot(insert_box) + nunits(molecule_type)
600 22 : 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 22 : IF (rem_atoms == 0) THEN
604 0 : ALLOCATE (remove_coords(1:3, 1:nunits(1)))
605 0 : ALLOCATE (atom_names_remove(1:nunits(1)))
606 : ELSE
607 66 : ALLOCATE (remove_coords(1:3, 1:rem_atoms))
608 66 : ALLOCATE (atom_names_remove(1:rem_atoms))
609 : END IF
610 66 : ALLOCATE (insert_coords(1:3, 1:ins_atoms))
611 66 : ALLOCATE (atom_names_insert(1:ins_atoms))
612 :
613 : ! grab the cells for later...acceptance and insertion
614 22 : IF (lbias) THEN
615 : CALL force_env_get(bias_env(insert_box)%force_env, &
616 22 : cell=cell_insert)
617 : CALL force_env_get(bias_env(remove_box)%force_env, &
618 22 : cell=cell_remove)
619 : ELSE
620 : CALL force_env_get(force_env(insert_box)%force_env, &
621 0 : cell=cell_insert)
622 : CALL force_env_get(force_env(remove_box)%force_env, &
623 0 : cell=cell_remove)
624 : END IF
625 22 : CALL get_cell(cell_remove, abc=abc_remove, deth=vol_remove)
626 22 : CALL get_cell(cell_insert, abc=abc_insert, deth=vol_insert)
627 :
628 22 : IF (ionode) THEN
629 : ! choose an insertion point
630 44 : DO idim = 1, 3
631 33 : rand = rng_stream%next()
632 44 : pos_insert(idim) = rand*abc_insert(idim)
633 : END DO
634 : END IF
635 22 : CALL group%bcast(pos_insert, source)
636 :
637 : ! allocate some arrays we'll be using
638 66 : ALLOCATE (r_insert_mol(1:3, 1:nunits(molecule_type)))
639 :
640 22 : iiatom = 1
641 64 : DO iatom = start_atom_rem, start_atom_rem + nunits(molecule_type) - 1
642 : r_insert_mol(1:3, iiatom) = &
643 168 : particles_old(remove_box)%list%els(iatom)%r(1:3)
644 64 : 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 22 : center_of_mass(:), mass(:, molecule_type))
650 :
651 : ! move the center of mass to the insertion point
652 88 : displace_molecule(1:3) = pos_insert(1:3) - center_of_mass(1:3)
653 64 : DO iatom = 1, nunits(molecule_type)
654 : r_insert_mol(1:3, iatom) = r_insert_mol(1:3, iatom) + &
655 190 : 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 66 : IF (SUM(nchains(:, insert_box)) == 0) THEN
662 0 : DO iatom = 1, nunits(molecule_type)
663 0 : insert_coords(1:3, iatom) = r_insert_mol(1:3, iatom)
664 : atom_names_insert(iatom) = &
665 0 : particles_old(remove_box)%list%els(start_atom_rem + iatom - 1)%atomic_kind%name
666 : END DO
667 0 : 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 22 : start_atom_ins = 1
675 34 : DO itype = 1, nmol_types
676 : start_atom_ins = start_atom_ins + &
677 34 : nchains(itype, insert_box)*nunits(itype)
678 34 : IF (itype == molecule_type) EXIT
679 : END DO
680 :
681 504 : DO iatom = 1, start_atom_ins - 1
682 : insert_coords(1:3, iatom) = &
683 1928 : particles_old(insert_box)%list%els(iatom)%r(1:3)
684 : atom_names_insert(iatom) = &
685 504 : particles_old(insert_box)%list%els(iatom)%atomic_kind%name
686 : END DO
687 22 : iiatom = 1
688 64 : DO iatom = start_atom_ins, start_atom_ins + nunits(molecule_type) - 1
689 168 : insert_coords(1:3, iatom) = r_insert_mol(1:3, iiatom)
690 42 : atom_names_insert(iatom) = atom_names(iiatom, molecule_type)
691 64 : iiatom = iiatom + 1
692 : END DO
693 76 : DO iatom = start_atom_ins + nunits(molecule_type), ins_atoms
694 : insert_coords(1:3, iatom) = &
695 216 : particles_old(insert_box)%list%els(iatom - nunits(molecule_type))%r(1:3)
696 : atom_names_insert(iatom) = &
697 76 : 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 22 : start_mol = 1
703 34 : DO jbox = 1, insert_box - 1
704 58 : start_mol = start_mol + SUM(nchains(:, jbox))
705 : END DO
706 66 : end_mol = start_mol + SUM(nchains(:, insert_box)) - 1
707 :
708 : ! make the .dat file
709 22 : IF (ionode) THEN
710 :
711 11 : nchains(molecule_type, insert_box) = nchains(molecule_type, insert_box) + 1
712 11 : IF (lbias) THEN
713 11 : 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 11 : mc_bias_file)
717 : ELSE
718 0 : 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 0 : mc_input_file)
722 : END IF
723 11 : 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 22 : IF (rem_atoms == 0) THEN
729 0 : DO iatom = 1, nunits(molecule_type)
730 0 : remove_coords(1:3, iatom) = r_insert_mol(1:3, iatom)
731 0 : 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 0 : nchains(molecule_type, remove_box) = nchains(molecule_type, remove_box) - 1
738 0 : IF (ionode) THEN
739 0 : IF (lbias) THEN
740 0 : 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 0 : mc_bias_file)
744 : ELSE
745 0 : 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 0 : mc_input_file)
749 : END IF
750 :
751 : END IF
752 0 : nchains(molecule_type, remove_box) = nchains(molecule_type, remove_box) + 1
753 :
754 : ELSE
755 368 : DO iatom = 1, start_atom_rem - 1
756 : remove_coords(1:3, iatom) = &
757 1384 : particles_old(remove_box)%list%els(iatom)%r(1:3)
758 : atom_names_remove(iatom) = &
759 368 : particles_old(remove_box)%list%els(iatom)%atomic_kind%name
760 : END DO
761 198 : DO iatom = start_atom_rem + nunits(molecule_type), nunits_tot(remove_box)
762 : remove_coords(1:3, iatom - nunits(molecule_type)) = &
763 704 : particles_old(remove_box)%list%els(iatom)%r(1:3)
764 : atom_names_remove(iatom - nunits(molecule_type)) = &
765 198 : particles_old(remove_box)%list%els(iatom)%atomic_kind%name
766 : END DO
767 :
768 : ! make the .dat file
769 22 : IF (ionode) THEN
770 11 : nchains(molecule_type, remove_box) = nchains(molecule_type, remove_box) - 1
771 11 : IF (lbias) THEN
772 11 : 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 11 : mc_bias_file)
776 : ELSE
777 0 : 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 0 : mc_input_file)
781 : END IF
782 11 : 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 22 : DEALLOCATE (r_insert_mol)
789 :
790 : ! now let's create the two new environments with the different number
791 : ! of molecules
792 66 : ALLOCATE (test_env(1:2))
793 : CALL mc_create_force_env(test_env(insert_box)%force_env, input_declaration, &
794 22 : para_env, dat_file(insert_box))
795 : CALL mc_create_force_env(test_env(remove_box)%force_env, input_declaration, &
796 22 : para_env, dat_file(remove_box))
797 :
798 : ! allocate an array we'll need
799 66 : ALLOCATE (r_cbmc(1:3, 1:ins_atoms))
800 66 : ALLOCATE (cbmc_energies(1:nswapmoves, 1:2))
801 :
802 22 : loverlap_ins = .FALSE.
803 22 : loverlap_rem = .FALSE.
804 :
805 : ! compute the new molecule information...we need this for the CBMC part
806 22 : IF (rem_atoms == 0) THEN
807 : CALL mc_determine_molecule_info(test_env, mc_molecule_info_test, &
808 0 : box_number=remove_box)
809 : ELSE
810 22 : 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 22 : mol_type=mol_type_test)
814 :
815 : ! figure out the position of the molecule we're inserting, and the
816 : ! Rosenbluth weight
817 22 : start_mol = 1
818 34 : DO jbox = 1, insert_box - 1
819 58 : start_mol = start_mol + SUM(nchains_test(:, jbox))
820 : END DO
821 66 : end_mol = start_mol + SUM(nchains_test(:, insert_box)) - 1
822 :
823 22 : 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 22 : 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 22 : 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 0 : nchains_test(:, insert_box), source, group, rng_stream)
844 : END IF
845 :
846 : CALL force_env_get(test_env(insert_box)%force_env, &
847 22 : subsys=insert_sys)
848 : CALL cp_subsys_get(insert_sys, &
849 22 : particles=particles_insert)
850 :
851 600 : DO iatom = 1, ins_atoms
852 2334 : 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 22 : IF (loverlap_ins .OR. loverlap_rem) THEN
858 : ! deallocate some stuff
859 0 : CALL mc_molecule_info_destroy(mc_molecule_info_test)
860 0 : CALL force_env_release(test_env(insert_box)%force_env)
861 0 : CALL force_env_release(test_env(remove_box)%force_env)
862 0 : DEALLOCATE (insert_coords)
863 0 : DEALLOCATE (remove_coords)
864 0 : DEALLOCATE (r_cbmc)
865 0 : DEALLOCATE (cbmc_energies)
866 0 : DEALLOCATE (oldsys)
867 0 : DEALLOCATE (particles_old)
868 0 : DEALLOCATE (test_env)
869 0 : CALL timestop(handle)
870 0 : 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 22 : subsys=insert_sys)
877 : CALL cp_subsys_get(insert_sys, &
878 22 : particles=particles_insert)
879 :
880 600 : DO iatom = 1, ins_atoms
881 : particles_insert%els(iatom)%r(1:3) = &
882 2334 : 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 22 : 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 22 : IF (lbias) THEN
892 :
893 66 : ALLOCATE (test_env_bias(1:2))
894 :
895 : ! first, the environment to which we added a molecule
896 22 : CALL get_mc_par(mc_par(insert_box)%mc_par, mc_input_file=mc_input_file)
897 22 : 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 11 : mc_input_file)
900 22 : test_env_bias(insert_box)%force_env => test_env(insert_box)%force_env
901 22 : NULLIFY (test_env(insert_box)%force_env)
902 : CALL mc_create_force_env(test_env(insert_box)%force_env, input_declaration, &
903 22 : para_env, dat_file(insert_box))
904 :
905 : CALL force_env_calc_energy_force(test_env(insert_box)%force_env, &
906 22 : calc_force=.FALSE.)
907 : CALL force_env_get(test_env(insert_box)%force_env, &
908 22 : potential_energy=new_energy(insert_box))
909 :
910 : ! now the environment that has one less molecule
911 66 : IF (SUM(nchains_test(:, remove_box)) == 0) THEN
912 0 : CALL get_mc_par(mc_par(remove_box)%mc_par, mc_input_file=mc_input_file)
913 0 : 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 0 : mc_input_file)
916 0 : test_env_bias(remove_box)%force_env => test_env(remove_box)%force_env
917 0 : NULLIFY (test_env(remove_box)%force_env)
918 : CALL mc_create_force_env(test_env(remove_box)%force_env, input_declaration, &
919 0 : para_env, dat_file(remove_box))
920 0 : new_energy(remove_box) = 0.0E0_dp
921 0 : bias_energy_new(remove_box) = 0.0E0_dp
922 : ELSE
923 22 : CALL get_mc_par(mc_par(remove_box)%mc_par, mc_input_file=mc_input_file)
924 22 : 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 11 : mc_input_file)
927 22 : test_env_bias(remove_box)%force_env => test_env(remove_box)%force_env
928 22 : NULLIFY (test_env(remove_box)%force_env)
929 : CALL mc_create_force_env(test_env(remove_box)%force_env, input_declaration, &
930 22 : para_env, dat_file(remove_box))
931 : CALL force_env_calc_energy_force(test_env(remove_box)%force_env, &
932 22 : calc_force=.FALSE.)
933 : CALL force_env_get(test_env(remove_box)%force_env, &
934 22 : potential_energy=new_energy(remove_box))
935 : CALL force_env_calc_energy_force(test_env_bias(remove_box)%force_env, &
936 22 : calc_force=.FALSE.)
937 : CALL force_env_get(test_env_bias(remove_box)%force_env, &
938 22 : potential_energy=bias_energy_new(remove_box))
939 : END IF
940 : ELSE
941 0 : IF (SUM(nchains_test(:, remove_box)) == 0) THEN
942 0 : new_energy(remove_box) = 0.0E0_dp
943 : ELSE
944 : CALL force_env_calc_energy_force(test_env(remove_box)%force_env, &
945 0 : calc_force=.FALSE.)
946 : CALL force_env_get(test_env(remove_box)%force_env, &
947 0 : 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 22 : start_mol = 1
957 32 : DO jbox = 1, remove_box - 1
958 52 : start_mol = start_mol + SUM(nchains(:, jbox))
959 : END DO
960 66 : end_mol = start_mol + SUM(nchains(:, remove_box)) - 1
961 22 : 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 22 : 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 0 : 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 22 : vol_insert/vol_remove
983 :
984 22 : 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 22 : - bias_energy_old(remove_box)))
991 :
992 22 : IF (del_quickstep_energy .GT. exp_max_val) THEN
993 0 : del_quickstep_energy = max_val
994 22 : ELSEIF (del_quickstep_energy .LT. exp_min_val) THEN
995 0 : del_quickstep_energy = min_val
996 : ELSE
997 22 : del_quickstep_energy = EXP(del_quickstep_energy)
998 : END IF
999 : w = prefactor*del_quickstep_energy*weight_new/weight_old &
1000 22 : *EXP(BETA*(eta_remove(molecule_type) - eta_insert(molecule_type)))
1001 :
1002 : ELSE
1003 : w = prefactor*weight_new/weight_old &
1004 0 : *EXP(BETA*(eta_remove(molecule_type) - eta_insert(molecule_type)))
1005 :
1006 : END IF
1007 :
1008 : ! check if the move is accepted
1009 22 : IF (w .GE. 1.0E0_dp) THEN
1010 8 : rand = 0.0E0_dp
1011 : ELSE
1012 14 : IF (ionode) rand = rng_stream%next()
1013 14 : CALL group%bcast(rand, source)
1014 : END IF
1015 :
1016 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1017 22 : 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 14 : 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 14 : IF (.NOT. lbias) THEN
1028 : new_energy(insert_box) = new_energy(insert_box) + &
1029 0 : old_energy(insert_box)
1030 : END IF
1031 :
1032 42 : DO ibox = 1, 2
1033 : ! update energies
1034 : energy_check(ibox) = energy_check(ibox) + (new_energy(ibox) - &
1035 28 : old_energy(ibox))
1036 28 : old_energy(ibox) = new_energy(ibox)
1037 : ! if we're biasing the update the biasing energy
1038 42 : IF (lbias) THEN
1039 28 : last_bias_energy(ibox) = bias_energy_new(ibox)
1040 28 : 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 14 : CALL mc_molecule_info_destroy(mc_molecule_info)
1049 14 : CALL set_mc_par(mc_par(insert_box)%mc_par, mc_molecule_info=mc_molecule_info_test)
1050 14 : 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 14 : subsys=insert_sys)
1055 : CALL cp_subsys_get(insert_sys, &
1056 14 : particles=particles_insert)
1057 376 : DO ipart = 1, ins_atoms
1058 1462 : 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 14 : subsys=remove_sys)
1062 : CALL cp_subsys_get(remove_sys, &
1063 14 : particles=particles_remove)
1064 352 : DO ipart = 1, rem_atoms
1065 1366 : r_old(1:3, ipart, remove_box) = particles_remove%els(ipart)%r(1:3)
1066 : END DO
1067 :
1068 : ! insertion box
1069 14 : CALL force_env_release(force_env(insert_box)%force_env)
1070 14 : force_env(insert_box)%force_env => test_env(insert_box)%force_env
1071 :
1072 : ! removal box
1073 14 : CALL force_env_release(force_env(remove_box)%force_env)
1074 14 : force_env(remove_box)%force_env => test_env(remove_box)%force_env
1075 :
1076 : ! if we're biasing, update the bias_env
1077 14 : IF (lbias) THEN
1078 14 : CALL force_env_release(bias_env(insert_box)%force_env)
1079 14 : bias_env(insert_box)%force_env => test_env_bias(insert_box)%force_env
1080 14 : CALL force_env_release(bias_env(remove_box)%force_env)
1081 14 : bias_env(remove_box)%force_env => test_env_bias(remove_box)%force_env
1082 14 : DEALLOCATE (test_env_bias)
1083 : END IF
1084 :
1085 : ELSE
1086 :
1087 : ! reject the move
1088 8 : CALL mc_molecule_info_destroy(mc_molecule_info_test)
1089 8 : CALL force_env_release(test_env(insert_box)%force_env)
1090 8 : CALL force_env_release(test_env(remove_box)%force_env)
1091 8 : IF (lbias) THEN
1092 8 : CALL force_env_release(test_env_bias(insert_box)%force_env)
1093 8 : CALL force_env_release(test_env_bias(remove_box)%force_env)
1094 8 : DEALLOCATE (test_env_bias)
1095 : END IF
1096 : END IF
1097 :
1098 : ! deallocate some stuff
1099 22 : DEALLOCATE (insert_coords)
1100 22 : DEALLOCATE (remove_coords)
1101 22 : DEALLOCATE (test_env)
1102 22 : DEALLOCATE (cbmc_energies)
1103 22 : DEALLOCATE (r_cbmc)
1104 22 : DEALLOCATE (oldsys)
1105 22 : DEALLOCATE (particles_old)
1106 :
1107 : ! end the timing
1108 22 : CALL timestop(handle)
1109 :
1110 66 : 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 24 : SUBROUTINE mc_ge_volume_move(mc_par, force_env, moves, move_updates, &
1133 24 : 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 24 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
1151 24 : INTEGER, DIMENSION(:, :), POINTER :: nchains
1152 : LOGICAL :: ionode
1153 24 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: loverlap
1154 : LOGICAL, DIMENSION(1:2) :: lempty
1155 24 : REAL(dp), DIMENSION(:, :), POINTER :: mass
1156 : REAL(KIND=dp) :: BETA, prefactor, rand, rmvolume, &
1157 : vol_dis, w
1158 24 : 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 24 : TYPE(cell_p_type), DIMENSION(:), POINTER :: cell, cell_old, cell_test
1164 24 : 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 24 : TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles_old
1169 :
1170 : ! begin the timing of the subroutine
1171 :
1172 24 : CALL timeset(routineN, handle)
1173 :
1174 : ! nullify some pointers
1175 24 : 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 24 : mc_molecule_info=mc_molecule_info)
1182 : CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot, &
1183 24 : mass=mass, nchains=nchains, nunits=nunits, mol_type=mol_type)
1184 :
1185 24 : print_level = 1
1186 24 : CALL get_mc_par(mc_par(2)%mc_par, dat_file=dat_file(2))
1187 :
1188 : ! allocate some stuff
1189 24 : max_atoms = MAX(nunits_tot(1), nunits_tot(2))
1190 96 : ALLOCATE (r(1:3, max_atoms, 1:2))
1191 72 : ALLOCATE (oldsys(1:2))
1192 72 : ALLOCATE (particles_old(1:2))
1193 72 : ALLOCATE (cell(1:2))
1194 72 : ALLOCATE (cell_test(1:2))
1195 72 : ALLOCATE (cell_old(1:2))
1196 24 : 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 72 : DO ibox = 1, 2
1201 48 : lempty(ibox) = .FALSE.
1202 168 : IF (SUM(nchains(:, ibox)) == 0) THEN
1203 0 : lempty(ibox) = .TRUE.
1204 : END IF
1205 : END DO
1206 :
1207 : ! record the attempt
1208 72 : DO ibox = 1, 2
1209 : moves(1, ibox)%moves%volume%attempts = &
1210 48 : moves(1, ibox)%moves%volume%attempts + 1
1211 : move_updates(1, ibox)%moves%volume%attempts = &
1212 72 : move_updates(1, ibox)%moves%volume%attempts + 1
1213 : END DO
1214 :
1215 : ! now let's grab the cell length and particle positions
1216 72 : DO ibox = 1, 2
1217 : CALL force_env_get(force_env(ibox)%force_env, &
1218 48 : subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
1219 48 : CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
1220 48 : NULLIFY (cell_old(ibox)%cell)
1221 48 : CALL cell_create(cell_old(ibox)%cell)
1222 48 : CALL cell_clone(cell(ibox)%cell, cell_old(ibox)%cell, tag="CELL_OLD")
1223 : CALL cp_subsys_get(oldsys(ibox)%subsys, &
1224 48 : particles=particles_old(ibox)%list)
1225 :
1226 : ! find the old cell length
1227 216 : old_cell_length(1:3, ibox) = abc(1:3, ibox)
1228 :
1229 : END DO
1230 :
1231 72 : DO ibox = 1, 2
1232 :
1233 : ! save the old coordinates
1234 1272 : DO iatom = 1, nunits_tot(ibox)
1235 4848 : 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 24 : IF (ionode) rand = rng_stream%next()
1242 24 : CALL group%bcast(rand, source)
1243 :
1244 24 : vol_dis = rmvolume*(rand - 0.5E0_dp)*2.0E0_dp
1245 :
1246 : ! add to one box, subtract from the other
1247 24 : 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 0 : CPABORT('GE_volume moves are trying to make box 1 smaller than 3')
1250 24 : 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 0 : CPABORT('GE_volume moves are trying to make box 2 smaller than 3')
1253 :
1254 96 : DO iside = 1, 3
1255 : new_cell_length(iside, 1) = (old_cell_length(1, 1)**3 + &
1256 72 : vol_dis)**(1.0E0_dp/3.0E0_dp)
1257 : new_cell_length(iside, 2) = (old_cell_length(1, 2)**3 - &
1258 96 : vol_dis)**(1.0E0_dp/3.0E0_dp)
1259 : END DO
1260 :
1261 : ! now we need to make the new cells
1262 72 : DO ibox = 1, 2
1263 624 : hmat_test(:, :, ibox) = 0.0e0_dp
1264 48 : hmat_test(1, 1, ibox) = new_cell_length(1, ibox)
1265 48 : hmat_test(2, 2, ibox) = new_cell_length(2, ibox)
1266 48 : hmat_test(3, 3, ibox) = new_cell_length(3, ibox)
1267 48 : NULLIFY (cell_test(ibox)%cell)
1268 : CALL cell_create(cell_test(ibox)%cell, hmat=hmat_test(:, :, ibox), &
1269 48 : periodic=cell(ibox)%cell%perd)
1270 48 : CALL force_env_get(force_env(ibox)%force_env, subsys=subsys)
1271 72 : CALL cp_subsys_set(subsys, cell=cell_test(ibox)%cell)
1272 : END DO
1273 :
1274 72 : DO ibox = 1, 2
1275 :
1276 : ! save the coords
1277 1248 : DO iatom = 1, nunits_tot(ibox)
1278 4848 : 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 72 : start_atom = 1
1284 : molecule_index = 1
1285 72 : DO jbox = 1, ibox - 1
1286 : IF (jbox == ibox) EXIT
1287 120 : molecule_index = molecule_index + SUM(nchains(:, jbox))
1288 : END DO
1289 720 : DO imolecule = 1, SUM(nchains(:, ibox))
1290 576 : molecule_type = mol_type(imolecule + molecule_index - 1)
1291 576 : IF (imolecule .NE. 1) THEN
1292 528 : start_atom = start_atom + nunits(mol_type(imolecule + molecule_index - 2))
1293 : END IF
1294 576 : 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 576 : 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 2304 : new_cell_length(1:3, ibox)/old_cell_length(1:3, ibox)
1304 2352 : DO j = 1, 3
1305 1728 : diff(j) = center_of_mass_new(j) - center_of_mass(j)
1306 : ! now change the particle positions
1307 5904 : DO jatom = start_atom, end_atom
1308 : particles_old(ibox)%list%els(jatom)%r(j) = &
1309 5328 : 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 72 : DO jbox = 1, ibox - 1
1318 120 : start_mol = start_mol + SUM(nchains(:, jbox))
1319 : END DO
1320 144 : 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 72 : cell_length=new_cell_length(:, ibox))
1324 :
1325 : END DO
1326 :
1327 : ! determine the overall energy difference
1328 :
1329 72 : DO ibox = 1, 2
1330 48 : IF (loverlap(ibox)) CYCLE
1331 : ! remake the force environment and calculate the energy
1332 72 : IF (lempty(ibox)) THEN
1333 0 : new_energy(ibox) = 0.0E0_dp
1334 : ELSE
1335 :
1336 : CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
1337 48 : calc_force=.FALSE.)
1338 : CALL force_env_get(force_env(ibox)%force_env, &
1339 48 : potential_energy=new_energy(ibox))
1340 :
1341 : END IF
1342 : END DO
1343 :
1344 : ! accept or reject the move
1345 72 : DO ibox = 1, 2
1346 : volume_new(ibox) = new_cell_length(1, ibox)* &
1347 48 : new_cell_length(2, ibox)*new_cell_length(3, ibox)
1348 : volume_old(ibox) = old_cell_length(1, ibox)* &
1349 72 : 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 120 : (volume_new(2)/volume_old(2))**(SUM(nchains(:, 2)))
1353 :
1354 24 : IF (loverlap(1) .OR. loverlap(2)) THEN
1355 0 : w = 0.0E0_dp
1356 : ELSE
1357 : w = prefactor*EXP(-BETA* &
1358 : (new_energy(1) + new_energy(2) - &
1359 24 : old_energy(1) - old_energy(2)))
1360 :
1361 : END IF
1362 :
1363 24 : IF (w .GE. 1.0E0_dp) THEN
1364 6 : w = 1.0E0_dp
1365 6 : rand = 0.0E0_dp
1366 : ELSE
1367 18 : IF (ionode) rand = rng_stream%next()
1368 18 : CALL group%bcast(rand, source)
1369 : END IF
1370 :
1371 24 : IF (rand .LT. w) THEN
1372 :
1373 : ! write cell length, volume, density, and trial displacement to a file
1374 18 : IF (ionode) THEN
1375 :
1376 9 : WRITE (cl, *) nnstep, new_cell_length(1, 1)* &
1377 9 : angstrom, vol_dis*(angstrom)**3, new_cell_length(1, 2)* &
1378 18 : angstrom
1379 9 : WRITE (cl, *) nnstep, new_energy(1), &
1380 18 : old_energy(1), new_energy(2), old_energy(2)
1381 9 : WRITE (cl, *) prefactor, w
1382 : END IF
1383 :
1384 54 : DO ibox = 1, 2
1385 : ! accept the move
1386 : moves(1, ibox)%moves%volume%successes = &
1387 36 : moves(1, ibox)%moves%volume%successes + 1
1388 : move_updates(1, ibox)%moves%volume%successes = &
1389 36 : move_updates(1, ibox)%moves%volume%successes + 1
1390 :
1391 : ! update energies
1392 : energy_check(ibox) = energy_check(ibox) + (new_energy(ibox) - &
1393 36 : old_energy(ibox))
1394 36 : old_energy(ibox) = new_energy(ibox)
1395 :
1396 : ! and the new "old" coordinates
1397 954 : DO iatom = 1, nunits_tot(ibox)
1398 : r_old(1:3, iatom, ibox) = &
1399 3636 : 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 6 : IF (ionode) THEN
1408 :
1409 3 : WRITE (cl, *) nnstep, new_cell_length(1, 1)* &
1410 3 : angstrom, vol_dis*(angstrom)**3, new_cell_length(1, 2)* &
1411 6 : angstrom
1412 3 : WRITE (cl, *) nnstep, new_energy(1), &
1413 6 : old_energy(1), new_energy(2), old_energy(2)
1414 3 : WRITE (cl, *) prefactor, w
1415 :
1416 : END IF
1417 :
1418 : ! reset the cell and particle positions
1419 18 : DO ibox = 1, 2
1420 12 : CALL force_env_get(force_env(ibox)%force_env, subsys=subsys)
1421 12 : CALL cp_subsys_set(subsys, cell=cell_old(ibox)%cell)
1422 318 : DO iatom = 1, nunits_tot(ibox)
1423 1212 : 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 72 : DO ibox = 1, 2
1431 48 : CALL cell_release(cell_test(ibox)%cell)
1432 72 : CALL cell_release(cell_old(ibox)%cell)
1433 : END DO
1434 24 : DEALLOCATE (r)
1435 24 : DEALLOCATE (oldsys)
1436 24 : DEALLOCATE (particles_old)
1437 24 : DEALLOCATE (cell)
1438 24 : DEALLOCATE (cell_old)
1439 24 : DEALLOCATE (cell_test)
1440 24 : DEALLOCATE (loverlap)
1441 :
1442 : ! end the timing
1443 24 : CALL timestop(handle)
1444 :
1445 24 : END SUBROUTINE mc_ge_volume_move
1446 :
1447 : END MODULE mc_ge_moves
1448 :
|
