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