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