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 Used to run the bulk of the MC simulation, doing things like
10 : !> choosing move types and writing data to files
11 : !> \author Matthew J. McGrath (09.26.2003)
12 : !>
13 : !> REVISIONS
14 : !> 09.10.05 MJM combined the two subroutines in this module into one
15 : ! **************************************************************************************************
16 : MODULE mc_ensembles
17 : USE cell_types, ONLY: cell_p_type,&
18 : get_cell
19 : USE cp_external_control, ONLY: external_control
20 : USE cp_files, ONLY: close_file,&
21 : open_file
22 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
23 : USE cp_subsys_types, ONLY: cp_subsys_get,&
24 : cp_subsys_p_type,&
25 : cp_subsys_type
26 : USE cp_units, ONLY: cp_unit_from_cp2k
27 : USE force_env_methods, ONLY: force_env_calc_energy_force
28 : USE force_env_types, ONLY: force_env_get,&
29 : force_env_p_type,&
30 : force_env_release
31 : USE global_types, ONLY: global_environment_type
32 : USE input_constants, ONLY: dump_xmol
33 : USE input_section_types, ONLY: section_type,&
34 : section_vals_type,&
35 : section_vals_val_get
36 : USE kinds, ONLY: default_string_length,&
37 : dp
38 : USE machine, ONLY: m_flush
39 : USE mathconstants, ONLY: pi
40 : USE mc_control, ONLY: mc_create_bias_force_env,&
41 : write_mc_restart
42 : USE mc_coordinates, ONLY: check_for_overlap,&
43 : create_discrete_array,&
44 : find_mc_test_molecule,&
45 : get_center_of_mass,&
46 : mc_coordinate_fold,&
47 : rotate_molecule
48 : USE mc_environment_types, ONLY: get_mc_env,&
49 : mc_environment_p_type,&
50 : set_mc_env
51 : USE mc_ge_moves, ONLY: mc_ge_swap_move,&
52 : mc_ge_volume_move,&
53 : mc_quickstep_move
54 : USE mc_misc, ONLY: final_mc_write,&
55 : mc_averages_create,&
56 : mc_averages_release
57 : USE mc_move_control, ONLY: init_mc_moves,&
58 : mc_move_update,&
59 : mc_moves_release,&
60 : write_move_stats
61 : USE mc_moves, ONLY: mc_avbmc_move,&
62 : mc_cluster_translation,&
63 : mc_conformation_change,&
64 : mc_hmc_move,&
65 : mc_molecule_rotation,&
66 : mc_molecule_translation,&
67 : mc_volume_move
68 : USE mc_types, ONLY: get_mc_molecule_info,&
69 : get_mc_par,&
70 : mc_averages_p_type,&
71 : mc_input_file_type,&
72 : mc_molecule_info_type,&
73 : mc_moves_p_type,&
74 : mc_simulation_parameters_p_type,&
75 : set_mc_par
76 : USE message_passing, ONLY: mp_comm_type,&
77 : mp_para_env_type
78 : USE parallel_rng_types, ONLY: rng_stream_type
79 : USE particle_list_types, ONLY: particle_list_p_type,&
80 : particle_list_type
81 : USE particle_methods, ONLY: write_particle_coordinates
82 : USE physcon, ONLY: angstrom,&
83 : boltzmann,&
84 : joule,&
85 : n_avogadro
86 : #include "../../base/base_uses.f90"
87 :
88 : IMPLICIT NONE
89 :
90 : PRIVATE
91 :
92 : ! *** Global parameters ***
93 :
94 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_ensembles'
95 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
96 :
97 : PUBLIC :: mc_run_ensemble, mc_compute_virial
98 :
99 : CONTAINS
100 :
101 : ! **************************************************************************************************
102 : !> \brief directs the program in running one or two box MC simulations
103 : !> \param mc_env a pointer that contains all mc_env for all the simulation
104 : !> boxes
105 : !> \param para_env ...
106 : !> \param globenv the global environment for the simulation
107 : !> \param input_declaration ...
108 : !> \param nboxes the number of simulation boxes
109 : !> \param rng_stream the stream we pull random numbers from
110 : !>
111 : !> Suitable for parallel.
112 : !> \author MJM
113 : ! **************************************************************************************************
114 18 : SUBROUTINE mc_run_ensemble(mc_env, para_env, globenv, input_declaration, nboxes, rng_stream)
115 :
116 : TYPE(mc_environment_p_type), DIMENSION(:), POINTER :: mc_env
117 : TYPE(mp_para_env_type), POINTER :: para_env
118 : TYPE(global_environment_type), POINTER :: globenv
119 : TYPE(section_type), POINTER :: input_declaration
120 : INTEGER, INTENT(IN) :: nboxes
121 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
122 :
123 : CHARACTER(len=*), PARAMETER :: routineN = 'mc_run_ensemble'
124 :
125 : CHARACTER(default_string_length), ALLOCATABLE, &
126 18 : DIMENSION(:) :: atom_names_box
127 : CHARACTER(default_string_length), &
128 18 : DIMENSION(:, :), POINTER :: atom_names
129 : CHARACTER(LEN=20) :: ensemble
130 : CHARACTER(LEN=40) :: cbox, cstep, fft_lib, move_type, &
131 : move_type_avbmc
132 18 : INTEGER, DIMENSION(:, :), POINTER :: nchains
133 18 : INTEGER, DIMENSION(:), POINTER :: avbmc_atom, mol_type, nchains_box, &
134 18 : nunits, nunits_tot
135 36 : INTEGER, DIMENSION(1:nboxes) :: box_flag, cl, data_unit, diff, istep, &
136 54 : move_unit, rm
137 : INTEGER, DIMENSION(1:3, 1:2) :: discrete_array
138 : INTEGER :: atom_number, box_number, cell_unit, com_crd, com_ene, com_mol, end_mol, handle, &
139 : ibox, idum, imol_type, imolecule, imove, iparticle, iprint, itype, iunit, iuptrans, &
140 : iupvolume, iw, jbox, jdum, molecule_type, molecule_type_swap, molecule_type_target, &
141 : nchain_total, nmol_types, nmoves, nnstep, nstart, nstep, source, start_atom, &
142 : start_atom_swap, start_atom_target, start_mol
143 : CHARACTER(LEN=default_string_length) :: unit_str
144 36 : CHARACTER(LEN=40), DIMENSION(1:nboxes) :: cell_file, coords_file, data_file, &
145 36 : displacement_file, energy_file, &
146 36 : molecules_file, moves_file
147 : LOGICAL :: ionode, lbias, ldiscrete, lhmc, &
148 : lnew_bias_env, loverlap, lreject, &
149 : lstop, print_kind, should_stop
150 18 : REAL(dp), DIMENSION(:), POINTER :: pbias, pmavbmc_mol, pmclus_box, &
151 18 : pmhmc_box, pmrot_mol, pmtraion_mol, &
152 18 : pmtrans_mol, pmvol_box
153 18 : REAL(dp), DIMENSION(:, :), POINTER :: conf_prob, mass
154 : REAL(KIND=dp) :: discrete_step, pmavbmc, pmcltrans, &
155 : pmhmc, pmswap, pmtraion, pmtrans, &
156 : pmvolume, rand, test_energy, unit_conv
157 18 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r_temp
158 18 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r_old
159 36 : REAL(KIND=dp), DIMENSION(1:3, 1:nboxes) :: abc
160 36 : REAL(KIND=dp), DIMENSION(1:nboxes) :: bias_energy, energy_check, final_energy, &
161 36 : initial_energy, last_bias_energy, &
162 36 : old_energy
163 18 : TYPE(cell_p_type), DIMENSION(:), POINTER :: cell
164 18 : TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: oldsys
165 : TYPE(cp_subsys_type), POINTER :: biassys
166 18 : TYPE(force_env_p_type), DIMENSION(:), POINTER :: bias_env, force_env
167 18 : TYPE(mc_averages_p_type), DIMENSION(:), POINTER :: averages
168 : TYPE(mc_input_file_type), POINTER :: mc_bias_file
169 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
170 18 : TYPE(mc_moves_p_type), DIMENSION(:), POINTER :: test_moves
171 18 : TYPE(mc_moves_p_type), DIMENSION(:, :), POINTER :: move_updates, moves
172 : TYPE(mc_simulation_parameters_p_type), &
173 18 : DIMENSION(:), POINTER :: mc_par
174 : TYPE(mp_comm_type) :: group
175 18 : TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles_old
176 : TYPE(particle_list_type), POINTER :: particles_bias
177 : TYPE(section_vals_type), POINTER :: root_section
178 :
179 18 : CALL timeset(routineN, handle)
180 :
181 : ! nullify some pointers
182 18 : NULLIFY (moves, move_updates, test_moves, root_section)
183 :
184 : ! allocate a whole bunch of stuff based on how many boxes we have
185 96 : ALLOCATE (force_env(1:nboxes))
186 60 : ALLOCATE (bias_env(1:nboxes))
187 60 : ALLOCATE (cell(1:nboxes))
188 60 : ALLOCATE (particles_old(1:nboxes))
189 60 : ALLOCATE (oldsys(1:nboxes))
190 60 : ALLOCATE (averages(1:nboxes))
191 60 : ALLOCATE (mc_par(1:nboxes))
192 36 : ALLOCATE (pmvol_box(1:nboxes))
193 36 : ALLOCATE (pmclus_box(1:nboxes))
194 36 : ALLOCATE (pmhmc_box(1:nboxes))
195 :
196 42 : DO ibox = 1, nboxes
197 : CALL get_mc_env(mc_env(ibox)%mc_env, &
198 : mc_par=mc_par(ibox)%mc_par, &
199 42 : force_env=force_env(ibox)%force_env)
200 : END DO
201 :
202 : ! Gather units of measure for output (if available)
203 18 : root_section => force_env(1)%force_env%root_section
204 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%UNIT", &
205 18 : c_val=unit_str)
206 18 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
207 : CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%PRINT_ATOM_KIND", &
208 18 : l_val=print_kind)
209 :
210 : ! get some data out of mc_par
211 : CALL get_mc_par(mc_par(1)%mc_par, &
212 : ionode=ionode, source=source, group=group, &
213 : data_file=data_file(1), moves_file=moves_file(1), &
214 : cell_file=cell_file(1), coords_file=coords_file(1), &
215 : energy_file=energy_file(1), displacement_file=displacement_file(1), &
216 : lstop=lstop, nstep=nstep, nstart=nstart, pmvolume=pmvolume, pmhmc=pmhmc, &
217 : molecules_file=molecules_file(1), pmswap=pmswap, nmoves=nmoves, &
218 : pmtraion=pmtraion, pmtrans=pmtrans, pmcltrans=pmcltrans, iuptrans=iuptrans, &
219 : iupvolume=iupvolume, ldiscrete=ldiscrete, pmtraion_mol=pmtraion_mol, &
220 : lbias=lbias, iprint=iprint, pmavbmc_mol=pmavbmc_mol, &
221 : discrete_step=discrete_step, fft_lib=fft_lib, avbmc_atom=avbmc_atom, &
222 : pmavbmc=pmavbmc, pbias=pbias, mc_molecule_info=mc_molecule_info, &
223 : pmrot_mol=pmrot_mol, pmtrans_mol=pmtrans_mol, pmvol_box=pmvol_box(1), &
224 18 : pmclus_box=pmclus_box(1), ensemble=ensemble, pmhmc_box=pmhmc_box(1), lhmc=lhmc)
225 :
226 : ! get some data from the molecule types
227 : CALL get_mc_molecule_info(mc_molecule_info, conf_prob=conf_prob, &
228 : nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
229 : mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
230 18 : atom_names=atom_names, mass=mass)
231 :
232 : ! allocate some stuff based on the number of molecule types we have
233 136 : ALLOCATE (moves(1:nmol_types, 1:nboxes))
234 118 : ALLOCATE (move_updates(1:nmol_types, 1:nboxes))
235 :
236 18 : IF (nboxes > 1) THEN
237 12 : DO ibox = 2, nboxes
238 : CALL get_mc_par(mc_par(ibox)%mc_par, &
239 : data_file=data_file(ibox), &
240 : moves_file=moves_file(ibox), &
241 : cell_file=cell_file(ibox), coords_file=coords_file(ibox), &
242 : energy_file=energy_file(ibox), &
243 : displacement_file=displacement_file(ibox), &
244 : molecules_file=molecules_file(ibox), pmvol_box=pmvol_box(ibox), &
245 12 : pmclus_box=pmclus_box(ibox), pmhmc_box=pmhmc_box(ibox))
246 : END DO
247 : END IF
248 :
249 : ! this is a check we can't do in the input checking
250 18 : IF (pmvol_box(nboxes) < 1.0E0_dp) THEN
251 0 : CPABORT('The last value of PMVOL_BOX needs to be 1.0')
252 : END IF
253 18 : IF (pmclus_box(nboxes) < 1.0E0_dp) THEN
254 0 : CPABORT('The last value of PMVOL_BOX needs to be 1.0')
255 : END IF
256 18 : IF (pmhmc_box(nboxes) < 1.0E0_dp) THEN
257 0 : CPABORT('The last value of PMHMC_BOX needs to be 1.0')
258 : END IF
259 :
260 : ! allocate the particle positions array for broadcasting
261 96 : ALLOCATE (r_old(3, SUM(nunits_tot), 1:nboxes))
262 :
263 : ! figure out what the default write unit is
264 18 : iw = cp_logger_get_default_io_unit()
265 :
266 18 : IF (iw > 0) THEN
267 9 : WRITE (iw, *)
268 9 : WRITE (iw, *)
269 9 : WRITE (iw, *) 'Beginning the Monte Carlo calculation.'
270 9 : WRITE (iw, *)
271 9 : WRITE (iw, *)
272 : END IF
273 :
274 : ! initialize running average variables
275 42 : energy_check(:) = 0.0E0_dp
276 42 : box_flag(:) = 0
277 42 : istep(:) = 0
278 :
279 42 : DO ibox = 1, nboxes
280 : ! initialize the moves array, the arrays for updating maximum move
281 : ! displacements, and the averages array
282 64 : DO itype = 1, nmol_types
283 40 : CALL init_mc_moves(moves(itype, ibox)%moves)
284 64 : CALL init_mc_moves(move_updates(itype, ibox)%moves)
285 : END DO
286 24 : CALL mc_averages_create(averages(ibox)%averages)
287 :
288 : ! find the energy of the initial configuration
289 64 : IF (SUM(nchains(:, ibox)) /= 0) THEN
290 : CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
291 24 : calc_force=.FALSE.)
292 : CALL force_env_get(force_env(ibox)%force_env, &
293 24 : potential_energy=old_energy(ibox))
294 : ELSE
295 0 : old_energy(ibox) = 0.0E0_dp
296 : END IF
297 24 : initial_energy(ibox) = old_energy(ibox)
298 :
299 : ! don't care about overlaps if we're only doing HMC
300 :
301 24 : IF (.NOT. lhmc) THEN
302 : ! check for overlaps
303 : start_mol = 1
304 28 : DO jbox = 1, ibox - 1
305 40 : start_mol = start_mol + SUM(nchains(:, jbox))
306 : END DO
307 60 : end_mol = start_mol + SUM(nchains(:, ibox)) - 1
308 : CALL check_for_overlap(force_env(ibox)%force_env, nchains(:, ibox), &
309 22 : nunits, loverlap, mol_type(start_mol:end_mol))
310 22 : IF (loverlap) CPABORT("overlap in an initial configuration")
311 : END IF
312 :
313 : ! get the subsystems and the cell information
314 : CALL force_env_get(force_env(ibox)%force_env, &
315 24 : subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
316 24 : CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
317 : CALL cp_subsys_get(oldsys(ibox)%subsys, &
318 24 : particles=particles_old(ibox)%list)
319 : ! record the old coordinates, in case a move is rejected
320 2656 : DO iparticle = 1, nunits_tot(ibox)
321 : r_old(1:3, iparticle, ibox) = &
322 10552 : particles_old(ibox)%list%els(iparticle)%r(1:3)
323 : END DO
324 :
325 : ! find the bias energy of the initial run
326 24 : IF (lbias) THEN
327 : ! determine the atom names of every particle
328 42 : ALLOCATE (atom_names_box(1:nunits_tot(ibox)))
329 :
330 14 : atom_number = 1
331 212 : DO imolecule = 1, SUM(nchains(:, ibox))
332 538 : DO iunit = 1, nunits(mol_type(imolecule + start_mol - 1))
333 : atom_names_box(atom_number) = &
334 354 : atom_names(iunit, mol_type(imolecule + start_mol - 1))
335 524 : atom_number = atom_number + 1
336 : END DO
337 : END DO
338 :
339 14 : CALL get_mc_par(mc_par(ibox)%mc_par, mc_bias_file=mc_bias_file)
340 14 : nchains_box => nchains(:, ibox)
341 : CALL mc_create_bias_force_env(bias_env(ibox)%force_env, &
342 : r_old(:, :, ibox), atom_names_box(:), nunits_tot(ibox), &
343 : para_env, abc(:, ibox), nchains_box, input_declaration, mc_bias_file, &
344 14 : ionode)
345 42 : IF (SUM(nchains(:, ibox)) /= 0) THEN
346 : CALL force_env_calc_energy_force(bias_env(ibox)%force_env, &
347 14 : calc_force=.FALSE.)
348 : CALL force_env_get(bias_env(ibox)%force_env, &
349 14 : potential_energy=last_bias_energy(ibox))
350 :
351 : ELSE
352 0 : last_bias_energy(ibox) = 0.0E0_dp
353 : END IF
354 14 : bias_energy(ibox) = last_bias_energy(ibox)
355 14 : DEALLOCATE (atom_names_box)
356 : END IF
357 42 : lnew_bias_env = .FALSE.
358 :
359 : END DO
360 :
361 : ! back to seriel for a bunch of I/O stuff
362 18 : IF (ionode) THEN
363 :
364 : ! record the combined energies,coordinates, and cell lengths
365 : CALL open_file(file_name='mc_cell_length', &
366 : unit_number=cell_unit, file_position='APPEND', &
367 9 : file_action='WRITE', file_status='UNKNOWN')
368 : CALL open_file(file_name='mc_energies', &
369 : unit_number=com_ene, file_position='APPEND', &
370 9 : file_action='WRITE', file_status='UNKNOWN')
371 : CALL open_file(file_name='mc_coordinates', &
372 : unit_number=com_crd, file_position='APPEND', &
373 9 : file_action='WRITE', file_status='UNKNOWN')
374 : CALL open_file(file_name='mc_molecules', &
375 : unit_number=com_mol, file_position='APPEND', &
376 9 : file_action='WRITE', file_status='UNKNOWN')
377 9 : WRITE (com_ene, *) 'Initial Energies: ', &
378 18 : old_energy(1:nboxes)
379 21 : DO ibox = 1, nboxes
380 12 : WRITE (com_mol, *) 'Initial Molecules: ', &
381 53 : nchains(:, ibox)
382 : END DO
383 21 : DO ibox = 1, nboxes
384 12 : WRITE (cell_unit, *) 'Initial: ', &
385 60 : abc(1:3, ibox)*angstrom
386 12 : WRITE (cbox, '(I4)') ibox
387 : CALL open_file(file_name='energy_differences_box'// &
388 : TRIM(ADJUSTL(cbox)), &
389 : unit_number=diff(ibox), file_position='APPEND', &
390 12 : file_action='WRITE', file_status='UNKNOWN')
391 32 : IF (SUM(nchains(:, ibox)) == 0) THEN
392 0 : WRITE (com_crd, *) ' 0'
393 0 : WRITE (com_crd, *) 'INITIAL BOX '//TRIM(ADJUSTL(cbox))
394 : ELSE
395 : CALL write_particle_coordinates(particles_old(ibox)%list%els, &
396 : com_crd, dump_xmol, 'POS', 'INITIAL BOX '//TRIM(ADJUSTL(cbox)), &
397 12 : unit_conv=unit_conv, print_kind=print_kind)
398 : END IF
399 : CALL open_file(file_name=data_file(ibox), &
400 : unit_number=data_unit(ibox), file_position='APPEND', &
401 12 : file_action='WRITE', file_status='UNKNOWN')
402 : CALL open_file(file_name=moves_file(ibox), &
403 : unit_number=move_unit(ibox), file_position='APPEND', &
404 12 : file_action='WRITE', file_status='UNKNOWN')
405 : CALL open_file(file_name=displacement_file(ibox), &
406 : unit_number=rm(ibox), file_position='APPEND', &
407 12 : file_action='WRITE', file_status='UNKNOWN')
408 : CALL open_file(file_name=cell_file(ibox), &
409 : unit_number=cl(ibox), file_position='APPEND', &
410 21 : file_action='WRITE', file_status='UNKNOWN')
411 :
412 : END DO
413 :
414 : ! back to parallel mode
415 : END IF
416 :
417 42 : DO ibox = 1, nboxes
418 24 : CALL group%bcast(cl(ibox), source)
419 24 : CALL group%bcast(rm(ibox), source)
420 24 : CALL group%bcast(diff(ibox), source)
421 : ! set all the units numbers that we just opened in the respective mc_par
422 : CALL set_mc_par(mc_par(ibox)%mc_par, cl=cl(ibox), rm=rm(ibox), &
423 42 : diff=diff(ibox))
424 : END DO
425 :
426 : ! if we're doing a discrete volume move, we need to set up the array
427 : ! that keeps track of which direction we can move in
428 18 : IF (ldiscrete) THEN
429 0 : IF (nboxes /= 1) THEN
430 0 : CPABORT('ldiscrete=.true. ONLY for systems with 1 box')
431 : END IF
432 : CALL create_discrete_array(abc(:, 1), discrete_array(:, :), &
433 0 : discrete_step)
434 : END IF
435 :
436 : ! find out how many steps we're doing...change the updates to be in cycles
437 : ! if the total number of steps is measured in cycles
438 18 : IF (.NOT. lstop) THEN
439 10 : nstep = nstep*nchain_total
440 10 : iuptrans = iuptrans*nchain_total
441 10 : iupvolume = iupvolume*nchain_total
442 : END IF
443 :
444 486 : DO nnstep = nstart + 1, nstart + nstep
445 :
446 468 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
447 15 : WRITE (iw, *)
448 15 : WRITE (iw, *) "------- On Monte Carlo Step ", nnstep
449 : END IF
450 :
451 468 : IF (ionode) rand = rng_stream%next()
452 : ! broadcast the random number, to make sure we're on the same move
453 468 : CALL group%bcast(rand, source)
454 :
455 468 : IF (rand < pmvolume) THEN
456 :
457 58 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
458 1 : WRITE (iw, *) "Attempting a volume move"
459 1 : WRITE (iw, *)
460 : END IF
461 :
462 8 : SELECT CASE (ensemble)
463 : CASE ("TRADITIONAL")
464 : CALL mc_volume_move(mc_par(1)%mc_par, &
465 : force_env(1)%force_env, &
466 : moves(1, 1)%moves, move_updates(1, 1)%moves, &
467 : old_energy(1), 1, &
468 : energy_check(1), r_old(:, :, 1), iw, discrete_array(:, :), &
469 8 : rng_stream)
470 : CASE ("GEMC_NVT")
471 : CALL mc_ge_volume_move(mc_par, force_env, moves, &
472 : move_updates, nnstep, old_energy, energy_check, &
473 24 : r_old, rng_stream)
474 : CASE ("GEMC_NPT")
475 : ! we need to select a box based on the probability given in the input file
476 26 : IF (ionode) rand = rng_stream%next()
477 26 : CALL group%bcast(rand, source)
478 :
479 38 : DO ibox = 1, nboxes
480 38 : IF (rand <= pmvol_box(ibox)) THEN
481 26 : box_number = ibox
482 26 : EXIT
483 : END IF
484 : END DO
485 :
486 : CALL mc_volume_move(mc_par(box_number)%mc_par, &
487 : force_env(box_number)%force_env, &
488 : moves(1, box_number)%moves, &
489 : move_updates(1, box_number)%moves, &
490 : old_energy(box_number), box_number, &
491 : energy_check(box_number), r_old(:, :, box_number), iw, &
492 : discrete_array(:, :), &
493 84 : rng_stream)
494 : END SELECT
495 :
496 : ! update all the pointers here, because otherwise we may pass wrong information when we're making a bias environment
497 166 : DO ibox = 1, nboxes
498 : CALL force_env_get(force_env(ibox)%force_env, &
499 108 : subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
500 108 : CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
501 : CALL cp_subsys_get(oldsys(ibox)%subsys, &
502 166 : particles=particles_old(ibox)%list)
503 : END DO
504 :
505 : ! we need a new biasing environment now, if we're into that sort of thing
506 58 : IF (lbias) THEN
507 150 : DO ibox = 1, nboxes
508 100 : CALL force_env_release(bias_env(ibox)%force_env)
509 : ! determine the atom names of every particle
510 300 : ALLOCATE (atom_names_box(1:nunits_tot(ibox)))
511 150 : start_mol = 1
512 150 : DO jbox = 1, ibox - 1
513 250 : start_mol = start_mol + SUM(nchains(:, jbox))
514 : END DO
515 300 : end_mol = start_mol + SUM(nchains(:, ibox)) - 1
516 300 : atom_number = 1
517 1500 : DO imolecule = 1, SUM(nchains(:, ibox))
518 3800 : DO iunit = 1, nunits(mol_type(imolecule + start_mol - 1))
519 : atom_names_box(atom_number) = &
520 2500 : atom_names(iunit, mol_type(imolecule + start_mol - 1))
521 3700 : atom_number = atom_number + 1
522 : END DO
523 : END DO
524 :
525 : ! need to find out what the cell lengths are
526 : CALL force_env_get(force_env(ibox)%force_env, &
527 100 : subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
528 100 : CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
529 :
530 : CALL get_mc_par(mc_par(ibox)%mc_par, &
531 100 : mc_bias_file=mc_bias_file)
532 100 : nchains_box => nchains(:, ibox)
533 :
534 : CALL mc_create_bias_force_env(bias_env(ibox)%force_env, &
535 : r_old(:, :, ibox), atom_names_box(:), nunits_tot(ibox), &
536 : para_env, abc(:, ibox), nchains_box, input_declaration, &
537 100 : mc_bias_file, ionode)
538 :
539 300 : IF (SUM(nchains(:, ibox)) /= 0) THEN
540 : CALL force_env_calc_energy_force( &
541 : bias_env(ibox)%force_env, &
542 100 : calc_force=.FALSE.)
543 : CALL force_env_get(bias_env(ibox)%force_env, &
544 100 : potential_energy=last_bias_energy(ibox))
545 : ELSE
546 0 : last_bias_energy(ibox) = 0.0E0_dp
547 : END IF
548 100 : bias_energy(ibox) = last_bias_energy(ibox)
549 150 : DEALLOCATE (atom_names_box)
550 : END DO
551 : END IF
552 :
553 410 : ELSE IF (rand < pmswap) THEN
554 :
555 : ! try a swap move
556 22 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
557 0 : WRITE (iw, *) "Attempting a swap move"
558 0 : WRITE (iw, *)
559 : END IF
560 :
561 : CALL mc_ge_swap_move(mc_par, force_env, bias_env, moves, &
562 : energy_check(:), r_old(:, :, :), old_energy(:), input_declaration, &
563 22 : para_env, bias_energy(:), last_bias_energy(:), rng_stream)
564 :
565 : ! the number of molecules may have changed, which deallocated the whole
566 : ! mc_molecule_info structure
567 22 : CALL get_mc_par(mc_par(1)%mc_par, mc_molecule_info=mc_molecule_info)
568 : CALL get_mc_molecule_info(mc_molecule_info, conf_prob=conf_prob, &
569 : nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
570 : mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
571 22 : atom_names=atom_names, mass=mass)
572 :
573 388 : ELSE IF (rand < pmhmc) THEN
574 : ! try hybrid Monte Carlo
575 20 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
576 2 : WRITE (iw, *) "Attempting a hybrid Monte Carlo move"
577 2 : WRITE (iw, *)
578 : END IF
579 :
580 : ! pick a box at random
581 20 : IF (ionode) rand = rng_stream%next()
582 20 : CALL group%bcast(rand, source)
583 :
584 20 : DO ibox = 1, nboxes
585 20 : IF (rand <= pmhmc_box(ibox)) THEN
586 20 : box_number = ibox
587 20 : EXIT
588 : END IF
589 : END DO
590 :
591 : CALL mc_hmc_move(mc_par(box_number)%mc_par, &
592 : force_env(box_number)%force_env, globenv, &
593 : moves(1, box_number)%moves, &
594 : move_updates(1, box_number)%moves, &
595 : old_energy(box_number), box_number, &
596 : energy_check(box_number), r_old(:, :, box_number), &
597 20 : rng_stream)
598 :
599 368 : ELSE IF (rand < pmavbmc) THEN
600 : ! try an AVBMC move
601 0 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
602 0 : WRITE (iw, *) "Attempting an AVBMC1 move"
603 0 : WRITE (iw, *)
604 : END IF
605 :
606 : ! first, pick a box to do it for
607 0 : IF (ionode) rand = rng_stream%next()
608 0 : CALL group%bcast(rand, source)
609 :
610 0 : IF (nboxes == 2) THEN
611 0 : IF (rand < 0.1E0_dp) THEN
612 0 : box_number = 1
613 : ELSE
614 0 : box_number = 2
615 : END IF
616 : ELSE
617 0 : box_number = 1
618 : END IF
619 :
620 : ! now pick a molecule type to do it for
621 0 : IF (ionode) rand = rng_stream%next()
622 0 : CALL group%bcast(rand, source)
623 0 : molecule_type_swap = 0
624 0 : DO imol_type = 1, nmol_types
625 0 : IF (rand < pmavbmc_mol(imol_type)) THEN
626 0 : molecule_type_swap = imol_type
627 0 : EXIT
628 : END IF
629 : END DO
630 0 : IF (molecule_type_swap == 0) THEN
631 0 : CPABORT('Did not choose a molecule type to swap...check AVBMC input')
632 : END IF
633 :
634 : ! now pick a molecule, automatically rejecting the move if the
635 : ! box is empty or only has one molecule
636 0 : IF (SUM(nchains(:, box_number)) <= 1) THEN
637 : ! indicate that we tried a move
638 : moves(molecule_type_swap, box_number)%moves%empty_avbmc = &
639 0 : moves(molecule_type_swap, box_number)%moves%empty_avbmc + 1
640 : ELSE
641 :
642 : ! pick a molecule to be swapped in the box
643 0 : IF (ionode) THEN
644 : CALL find_mc_test_molecule(mc_molecule_info, &
645 : start_atom_swap, idum, jdum, rng_stream, &
646 0 : box=box_number, molecule_type_old=molecule_type_swap)
647 :
648 : ! pick a molecule to act as the target in the box...we don't care what type
649 0 : DO
650 : CALL find_mc_test_molecule(mc_molecule_info, &
651 : start_atom_target, idum, molecule_type_target, &
652 0 : rng_stream, box=box_number)
653 0 : IF (start_atom_swap /= start_atom_target) THEN
654 : start_atom_target = start_atom_target + &
655 0 : avbmc_atom(molecule_type_target) - 1
656 : EXIT
657 : END IF
658 : END DO
659 :
660 : ! choose if we're swapping into the bonded region of mol_target, or
661 : ! into the nonbonded region
662 0 : rand = rng_stream%next()
663 :
664 : END IF
665 0 : CALL group%bcast(start_atom_swap, source)
666 0 : CALL group%bcast(box_number, source)
667 0 : CALL group%bcast(start_atom_target, source)
668 0 : CALL group%bcast(rand, source)
669 :
670 0 : IF (rand < pbias(molecule_type_swap)) THEN
671 0 : move_type_avbmc = 'in'
672 : ELSE
673 0 : move_type_avbmc = 'out'
674 : END IF
675 :
676 : CALL mc_avbmc_move(mc_par(box_number)%mc_par, &
677 : force_env(box_number)%force_env, &
678 : bias_env(box_number)%force_env, &
679 : moves(molecule_type_swap, box_number)%moves, &
680 : energy_check(box_number), &
681 : r_old(:, :, box_number), old_energy(box_number), &
682 : start_atom_swap, start_atom_target, molecule_type_swap, &
683 : box_number, bias_energy(box_number), &
684 : last_bias_energy(box_number), &
685 0 : move_type_avbmc, rng_stream)
686 :
687 : END IF
688 :
689 : ELSE
690 :
691 368 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
692 12 : WRITE (iw, *) "Attempting an inner move"
693 12 : WRITE (iw, *)
694 : END IF
695 :
696 2010 : DO imove = 1, nmoves
697 :
698 1642 : IF (ionode) rand = rng_stream%next()
699 1642 : CALL group%bcast(rand, source)
700 2010 : IF (rand < pmtraion) THEN
701 : ! change molecular conformation
702 : ! first, pick a box to do it for
703 506 : IF (ionode) rand = rng_stream%next()
704 506 : CALL group%bcast(rand, source)
705 506 : IF (nboxes == 2) THEN
706 0 : IF (rand < 0.75E0_dp) THEN
707 0 : box_number = 1
708 : ELSE
709 0 : box_number = 2
710 : END IF
711 : ELSE
712 506 : box_number = 1
713 : END IF
714 :
715 : ! figure out which molecule type we're looking for
716 506 : IF (ionode) rand = rng_stream%next()
717 506 : CALL group%bcast(rand, source)
718 506 : molecule_type = 0
719 506 : DO imol_type = 1, nmol_types
720 506 : IF (rand < pmtraion_mol(imol_type)) THEN
721 506 : molecule_type = imol_type
722 506 : EXIT
723 : END IF
724 : END DO
725 506 : IF (molecule_type == 0) CALL cp_abort( &
726 : __LOCATION__, &
727 0 : 'Did not choose a molecule type to conf change...PMTRAION_MOL should not be all 0.0')
728 :
729 : ! now pick a molecule, automatically rejecting the move if the
730 : ! box is empty
731 506 : IF (nchains(molecule_type, box_number) == 0) THEN
732 : ! indicate that we tried a move
733 : moves(molecule_type, box_number)%moves%empty_conf = &
734 0 : moves(molecule_type, box_number)%moves%empty_conf + 1
735 : ELSE
736 : ! pick a molecule in the box
737 506 : IF (ionode) THEN
738 : CALL find_mc_test_molecule(mc_molecule_info, &
739 : start_atom, idum, &
740 : jdum, rng_stream, &
741 253 : box=box_number, molecule_type_old=molecule_type)
742 :
743 : ! choose if we're changing a bond length or an angle
744 253 : rand = rng_stream%next()
745 : END IF
746 506 : CALL group%bcast(rand, source)
747 506 : CALL group%bcast(start_atom, source)
748 506 : CALL group%bcast(box_number, source)
749 506 : CALL group%bcast(molecule_type, source)
750 :
751 : ! figure out what kind of move we're doing
752 506 : IF (rand < conf_prob(1, molecule_type)) THEN
753 312 : move_type = 'bond'
754 194 : ELSE IF (rand < (conf_prob(1, molecule_type) + &
755 : conf_prob(2, molecule_type))) THEN
756 194 : move_type = 'angle'
757 : ELSE
758 0 : move_type = 'dihedral'
759 : END IF
760 506 : box_flag(box_number) = 1
761 : CALL mc_conformation_change(mc_par(box_number)%mc_par, &
762 : force_env(box_number)%force_env, &
763 : bias_env(box_number)%force_env, &
764 : moves(molecule_type, box_number)%moves, &
765 : move_updates(molecule_type, box_number)%moves, &
766 : start_atom, molecule_type, box_number, &
767 : bias_energy(box_number), &
768 506 : move_type, lreject, rng_stream)
769 506 : IF (lreject) EXIT
770 : END IF
771 1136 : ELSE IF (rand < pmtrans) THEN
772 : ! translate a whole molecule in the system
773 : ! pick a molecule type
774 624 : IF (ionode) rand = rng_stream%next()
775 624 : CALL group%bcast(rand, source)
776 624 : molecule_type = 0
777 924 : DO imol_type = 1, nmol_types
778 924 : IF (rand < pmtrans_mol(imol_type)) THEN
779 624 : molecule_type = imol_type
780 624 : EXIT
781 : END IF
782 : END DO
783 624 : IF (molecule_type == 0) CALL cp_abort( &
784 : __LOCATION__, &
785 0 : 'Did not choose a molecule type to translate...PMTRANS_MOL should not be all 0.0')
786 :
787 : ! now pick a molecule of that type
788 624 : IF (ionode) THEN
789 : CALL find_mc_test_molecule(mc_molecule_info, &
790 : start_atom, box_number, idum, rng_stream, &
791 312 : molecule_type_old=molecule_type)
792 : END IF
793 624 : CALL group%bcast(start_atom, source)
794 624 : CALL group%bcast(box_number, source)
795 624 : box_flag(box_number) = 1
796 : CALL mc_molecule_translation(mc_par(box_number)%mc_par, &
797 : force_env(box_number)%force_env, &
798 : bias_env(box_number)%force_env, &
799 : moves(molecule_type, box_number)%moves, &
800 : move_updates(molecule_type, box_number)%moves, &
801 : start_atom, box_number, bias_energy(box_number), &
802 624 : molecule_type, lreject, rng_stream)
803 624 : IF (lreject) EXIT
804 512 : ELSE IF (rand < pmcltrans) THEN
805 : ! translate a whole cluster in the system
806 : ! first, pick a box to do it for
807 10 : IF (ionode) rand = rng_stream%next()
808 10 : CALL group%bcast(rand, source)
809 :
810 10 : DO ibox = 1, nboxes
811 10 : IF (rand <= pmclus_box(ibox)) THEN
812 10 : box_number = ibox
813 10 : EXIT
814 : END IF
815 : END DO
816 10 : box_flag(box_number) = 1
817 : CALL mc_cluster_translation(mc_par(box_number)%mc_par, &
818 : force_env(box_number)%force_env, &
819 : bias_env(box_number)%force_env, &
820 : moves(1, box_number)%moves, &
821 : move_updates(1, box_number)%moves, &
822 : box_number, bias_energy(box_number), &
823 10 : lreject, rng_stream)
824 10 : IF (lreject) EXIT
825 : ELSE
826 : ! rotate a whole molecule in the system
827 : ! pick a molecule type
828 502 : IF (ionode) rand = rng_stream%next()
829 502 : CALL group%bcast(rand, source)
830 502 : molecule_type = 0
831 502 : DO imol_type = 1, nmol_types
832 502 : IF (rand < pmrot_mol(imol_type)) THEN
833 502 : molecule_type = imol_type
834 502 : EXIT
835 : END IF
836 : END DO
837 502 : IF (molecule_type == 0) CALL cp_abort( &
838 : __LOCATION__, &
839 0 : 'Did not choose a molecule type to rotate...PMROT_MOL should not be all 0.0')
840 :
841 502 : IF (ionode) THEN
842 : CALL find_mc_test_molecule(mc_molecule_info, &
843 : start_atom, box_number, idum, rng_stream, &
844 251 : molecule_type_old=molecule_type)
845 : END IF
846 502 : CALL group%bcast(start_atom, source)
847 502 : CALL group%bcast(box_number, source)
848 502 : box_flag(box_number) = 1
849 : CALL mc_molecule_rotation(mc_par(box_number)%mc_par, &
850 : force_env(box_number)%force_env, &
851 : bias_env(box_number)%force_env, &
852 : moves(molecule_type, box_number)%moves, &
853 : move_updates(molecule_type, box_number)%moves, &
854 : box_number, start_atom, &
855 : molecule_type, bias_energy(box_number), &
856 502 : lreject, rng_stream)
857 502 : IF (lreject) EXIT
858 : END IF
859 :
860 : END DO
861 :
862 : ! now do a Quickstep calculation to see if we accept the sequence
863 : CALL mc_Quickstep_move(mc_par, force_env, bias_env, &
864 : moves, lreject, move_updates, energy_check(:), r_old(:, :, :), &
865 : nnstep, old_energy(:), bias_energy(:), last_bias_energy(:), &
866 : nboxes, box_flag(:), oldsys, particles_old, &
867 368 : rng_stream, unit_conv)
868 :
869 : END IF
870 :
871 : ! make sure the pointers are pointing correctly since the subsys may
872 : ! have changed
873 1080 : DO ibox = 1, nboxes
874 : CALL force_env_get(force_env(ibox)%force_env, &
875 612 : subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
876 612 : CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
877 : CALL cp_subsys_get(oldsys(ibox)%subsys, &
878 1080 : particles=particles_old(ibox)%list)
879 : END DO
880 :
881 468 : IF (ionode) THEN
882 :
883 234 : IF (MOD(nnstep, iprint) == 0) THEN
884 15 : WRITE (com_ene, *) nnstep, old_energy(1:nboxes)
885 :
886 33 : DO ibox = 1, nboxes
887 :
888 : ! write the molecule information
889 47 : WRITE (com_mol, *) nnstep, nchains(:, ibox)
890 :
891 : ! write the move statistics to file
892 47 : DO itype = 1, nmol_types
893 : CALL write_move_stats(moves(itype, ibox)%moves, &
894 47 : nnstep, move_unit(ibox))
895 : END DO
896 :
897 : ! write a restart file
898 : CALL write_mc_restart(nnstep, mc_par(ibox)%mc_par, &
899 18 : nchains(:, ibox), force_env(ibox)%force_env)
900 :
901 : ! write cell lengths
902 72 : WRITE (cell_unit, *) nnstep, abc(1:3, ibox)*angstrom
903 :
904 : ! write particle coordinates
905 18 : WRITE (cbox, '(I4)') ibox
906 18 : WRITE (cstep, '(I8)') nnstep
907 62 : IF (SUM(nchains(:, ibox)) == 0) THEN
908 0 : WRITE (com_crd, *) ' 0'
909 : WRITE (com_crd, *) 'BOX '//TRIM(ADJUSTL(cbox))// &
910 0 : ', STEP '//TRIM(ADJUSTL(cstep))
911 : ELSE
912 : CALL write_particle_coordinates( &
913 : particles_old(ibox)%list%els, &
914 : com_crd, dump_xmol, 'POS', &
915 : 'BOX '//TRIM(ADJUSTL(cbox))// &
916 : ', STEP '//TRIM(ADJUSTL(cstep)), &
917 18 : unit_conv=unit_conv)
918 : END IF
919 : END DO
920 : END IF ! end the things we only do every iprint moves
921 :
922 540 : DO ibox = 1, nboxes
923 : ! compute some averages
924 : averages(ibox)%averages%ave_energy = &
925 : averages(ibox)%averages%ave_energy*REAL(nnstep - &
926 : nstart - 1, dp)/REAL(nnstep - nstart, dp) + &
927 306 : old_energy(ibox)/REAL(nnstep - nstart, dp)
928 : averages(ibox)%averages%molecules = &
929 : averages(ibox)%averages%molecules*REAL(nnstep - &
930 : nstart - 1, dp)/REAL(nnstep - nstart, dp) + &
931 899 : REAL(SUM(nchains(:, ibox)), dp)/REAL(nnstep - nstart, dp)
932 : averages(ibox)%averages%ave_volume = &
933 : averages(ibox)%averages%ave_volume* &
934 : REAL(nnstep - nstart - 1, dp)/REAL(nnstep - nstart, dp) + &
935 : abc(1, ibox)*abc(2, ibox)*abc(3, ibox)/ &
936 306 : REAL(nnstep - nstart, dp)
937 :
938 : ! flush the buffers to the files
939 306 : CALL m_flush(data_unit(ibox))
940 306 : CALL m_flush(diff(ibox))
941 306 : CALL m_flush(move_unit(ibox))
942 306 : CALL m_flush(cl(ibox))
943 540 : CALL m_flush(rm(ibox))
944 :
945 : END DO
946 :
947 : ! flush more buffers to the files
948 234 : CALL m_flush(cell_unit)
949 234 : CALL m_flush(com_ene)
950 234 : CALL m_flush(com_crd)
951 234 : CALL m_flush(com_mol)
952 :
953 : END IF
954 :
955 : ! reset the box flags
956 1080 : box_flag(:) = 0
957 :
958 : ! check to see if EXIT file exists...if so, end the calculation
959 468 : CALL external_control(should_stop, "MC", globenv=globenv)
960 468 : IF (should_stop) EXIT
961 :
962 : ! update the move displacements, if necessary
963 1080 : DO ibox = 1, nboxes
964 612 : IF (MOD(nnstep - nstart, iuptrans) == 0) THEN
965 0 : DO itype = 1, nmol_types
966 : CALL mc_move_update(mc_par(ibox)%mc_par, &
967 : move_updates(itype, ibox)%moves, itype, &
968 0 : "trans", nnstep, ionode)
969 : END DO
970 : END IF
971 :
972 1080 : IF (MOD(nnstep - nstart, iupvolume) == 0) THEN
973 : CALL mc_move_update(mc_par(ibox)%mc_par, &
974 : move_updates(1, ibox)%moves, 1337, &
975 0 : "volume", nnstep, ionode)
976 : END IF
977 : END DO
978 :
979 : ! check to see if there are any overlaps in the boxes, and fold coordinates
980 : ! don't care about overlaps if we're only doing HMC
981 468 : IF (.NOT. lhmc) THEN
982 1040 : DO ibox = 1, nboxes
983 2206 : IF (SUM(nchains(:, ibox)) /= 0) THEN
984 : start_mol = 1
985 736 : DO jbox = 1, ibox - 1
986 1024 : start_mol = start_mol + SUM(nchains(:, jbox))
987 : END DO
988 1758 : end_mol = start_mol + SUM(nchains(:, ibox)) - 1
989 : CALL check_for_overlap(force_env(ibox)%force_env, &
990 : nchains(:, ibox), nunits, loverlap, &
991 592 : mol_type(start_mol:end_mol))
992 592 : IF (loverlap) THEN
993 0 : IF (iw > 0) WRITE (iw, *) nnstep
994 0 : CPABORT('coordinate overlap at the end of the above step')
995 : ! now fold the coordinates...don't do this anywhere but here, because
996 : ! we can get screwed up with the mc_molecule_info stuff (like in swap move)...
997 : ! this is kind of ugly, with allocated and deallocating every time
998 0 : ALLOCATE (r_temp(1:3, 1:nunits_tot(ibox)))
999 :
1000 0 : DO iunit = 1, nunits_tot(ibox)
1001 : r_temp(1:3, iunit) = &
1002 0 : particles_old(ibox)%list%els(iunit)%r(1:3)
1003 : END DO
1004 :
1005 : CALL mc_coordinate_fold(r_temp(:, :), &
1006 : SUM(nchains(:, ibox)), mol_type(start_mol:end_mol), &
1007 0 : mass, nunits, abc(1:3, ibox))
1008 :
1009 : ! save the folded coordinates
1010 0 : DO iunit = 1, nunits_tot(ibox)
1011 0 : r_old(1:3, iunit, ibox) = r_temp(1:3, iunit)
1012 : particles_old(ibox)%list%els(iunit)%r(1:3) = &
1013 0 : r_temp(1:3, iunit)
1014 : END DO
1015 :
1016 : ! if we're biasing, we need to do the same
1017 0 : IF (lbias) THEN
1018 : CALL force_env_get(bias_env(ibox)%force_env, &
1019 0 : subsys=biassys)
1020 : CALL cp_subsys_get(biassys, &
1021 0 : particles=particles_bias)
1022 :
1023 0 : DO iunit = 1, nunits_tot(ibox)
1024 : particles_bias%els(iunit)%r(1:3) = &
1025 0 : r_temp(1:3, iunit)
1026 : END DO
1027 : END IF
1028 :
1029 0 : DEALLOCATE (r_temp)
1030 : END IF
1031 : END IF
1032 : END DO
1033 : END IF
1034 :
1035 : !debug code
1036 486 : IF (debug_this_module) THEN
1037 : DO ibox = 1, nboxes
1038 : IF (SUM(nchains(:, ibox)) /= 0) THEN
1039 : CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
1040 : calc_force=.FALSE.)
1041 : CALL force_env_get(force_env(ibox)%force_env, &
1042 : potential_energy=test_energy)
1043 : ELSE
1044 : test_energy = 0.0E0_dp
1045 : END IF
1046 :
1047 : IF (ABS(initial_energy(ibox) + energy_check(ibox) - &
1048 : test_energy) > 0.0000001E0_dp) THEN
1049 : IF (iw > 0) THEN
1050 : WRITE (iw, *) '!!!!!!! We have an energy problem. !!!!!!!!'
1051 : WRITE (iw, '(A,T64,F16.10)') 'Final Energy = ', test_energy
1052 : WRITE (iw, '(A,T64,F16.10)') 'Initial Energy+energy_check=', &
1053 : initial_energy(ibox) + energy_check(ibox)
1054 : WRITE (iw, *) 'Box ', ibox
1055 : WRITE (iw, *) 'nchains ', nchains(:, ibox)
1056 : END IF
1057 : CPABORT('!!!!!!! We have an energy problem. !!!!!!!!')
1058 : END IF
1059 : END DO
1060 : END IF
1061 : END DO
1062 :
1063 : ! write a restart file
1064 18 : IF (ionode) THEN
1065 21 : DO ibox = 1, nboxes
1066 : CALL write_mc_restart(nnstep, mc_par(ibox)%mc_par, &
1067 21 : nchains(:, ibox), force_env(ibox)%force_env)
1068 : END DO
1069 : END IF
1070 :
1071 : ! calculate the final energy
1072 42 : DO ibox = 1, nboxes
1073 64 : IF (SUM(nchains(:, ibox)) /= 0) THEN
1074 : CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
1075 24 : calc_force=.FALSE.)
1076 : CALL force_env_get(force_env(ibox)%force_env, &
1077 24 : potential_energy=final_energy(ibox))
1078 : ELSE
1079 0 : final_energy(ibox) = 0.0E0_dp
1080 : END IF
1081 42 : IF (lbias) THEN
1082 14 : CALL force_env_release(bias_env(ibox)%force_env)
1083 : END IF
1084 : END DO
1085 :
1086 : ! do some stuff in serial
1087 18 : IF (ionode .OR. (iw > 0)) THEN
1088 :
1089 9 : WRITE (com_ene, *) 'Final Energies: ', &
1090 18 : final_energy(1:nboxes)
1091 :
1092 21 : DO ibox = 1, nboxes
1093 12 : WRITE (cbox, '(I4)') ibox
1094 32 : IF (SUM(nchains(:, ibox)) == 0) THEN
1095 0 : WRITE (com_crd, *) ' 0'
1096 0 : WRITE (com_crd, *) 'BOX '//TRIM(ADJUSTL(cbox))
1097 : ELSE
1098 : CALL write_particle_coordinates( &
1099 : particles_old(ibox)%list%els, &
1100 : com_crd, dump_xmol, 'POS', &
1101 12 : 'FINAL BOX '//TRIM(ADJUSTL(cbox)), unit_conv=unit_conv)
1102 : END IF
1103 :
1104 : ! write a bunch of data to the screen
1105 : WRITE (iw, '(A)') &
1106 12 : '------------------------------------------------'
1107 : WRITE (iw, '(A,I1,A)') &
1108 12 : '| BOX ', ibox, &
1109 24 : ' |'
1110 : WRITE (iw, '(A)') &
1111 12 : '------------------------------------------------'
1112 12 : test_moves => moves(:, ibox)
1113 : CALL final_mc_write(mc_par(ibox)%mc_par, test_moves, &
1114 : iw, energy_check(ibox), &
1115 : initial_energy(ibox), final_energy(ibox), &
1116 12 : averages(ibox)%averages)
1117 :
1118 : ! close any open files
1119 12 : CALL close_file(unit_number=diff(ibox))
1120 12 : CALL close_file(unit_number=data_unit(ibox))
1121 12 : CALL close_file(unit_number=move_unit(ibox))
1122 12 : CALL close_file(unit_number=cl(ibox))
1123 21 : CALL close_file(unit_number=rm(ibox))
1124 : END DO
1125 :
1126 : ! close some more files
1127 9 : CALL close_file(unit_number=cell_unit)
1128 9 : CALL close_file(unit_number=com_ene)
1129 9 : CALL close_file(unit_number=com_crd)
1130 9 : CALL close_file(unit_number=com_mol)
1131 : END IF
1132 :
1133 42 : DO ibox = 1, nboxes
1134 : CALL set_mc_env(mc_env(ibox)%mc_env, &
1135 : mc_par=mc_par(ibox)%mc_par, &
1136 24 : force_env=force_env(ibox)%force_env)
1137 :
1138 : ! deallocate some stuff
1139 64 : DO itype = 1, nmol_types
1140 40 : CALL mc_moves_release(move_updates(itype, ibox)%moves)
1141 64 : CALL mc_moves_release(moves(itype, ibox)%moves)
1142 : END DO
1143 42 : CALL mc_averages_release(averages(ibox)%averages)
1144 : END DO
1145 :
1146 18 : DEALLOCATE (pmhmc_box)
1147 18 : DEALLOCATE (pmvol_box)
1148 18 : DEALLOCATE (pmclus_box)
1149 18 : DEALLOCATE (r_old)
1150 18 : DEALLOCATE (force_env)
1151 18 : DEALLOCATE (bias_env)
1152 18 : DEALLOCATE (cell)
1153 18 : DEALLOCATE (particles_old)
1154 18 : DEALLOCATE (oldsys)
1155 18 : DEALLOCATE (averages)
1156 18 : DEALLOCATE (moves)
1157 18 : DEALLOCATE (move_updates)
1158 18 : DEALLOCATE (mc_par)
1159 :
1160 : ! end the timing
1161 18 : CALL timestop(handle)
1162 :
1163 54 : END SUBROUTINE mc_run_ensemble
1164 :
1165 : ! **************************************************************************************************
1166 : !> \brief Computes the second virial coefficient of a molecule by using the integral form
1167 : !> of the second virial coefficient found in McQuarrie "Statistical Thermodynamics",
1168 : !> B2(T) = -2Pi Int 0toInf [ Exp[-beta*u(r)] -1] r^2 dr Eq. 15-25
1169 : !> I use trapazoidal integration with various step sizes
1170 : !> (the integral is broken up into three parts, currently, but that's easily
1171 : !> changed by the first variables found below). It generates nvirial configurations,
1172 : !> doing the integration for each one, and then averages all the B2(T) to produce
1173 : !> the final answer.
1174 : !> \param mc_env a pointer that contains all mc_env for all the simulation
1175 : !> boxes
1176 : !> \param rng_stream the stream we pull random numbers from
1177 : !>
1178 : !> Suitable for parallel.
1179 : !> \author MJM
1180 : ! **************************************************************************************************
1181 2 : SUBROUTINE mc_compute_virial(mc_env, rng_stream)
1182 :
1183 : TYPE(mc_environment_p_type), DIMENSION(:), POINTER :: mc_env
1184 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1185 :
1186 : INTEGER :: current_division, end_atom, ibin, idivision, iparticle, iprint, itemp, iunit, &
1187 : ivirial, iw, nbins, nchain_total, nintegral_divisions, nmol_types, nvirial, &
1188 : nvirial_temps, source, start_atom
1189 2 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
1190 2 : INTEGER, DIMENSION(:, :), POINTER :: nchains
1191 : LOGICAL :: ionode, loverlap
1192 2 : REAL(dp), DIMENSION(:), POINTER :: BETA, virial_cutoffs, virial_stepsize, &
1193 2 : virial_temps
1194 2 : REAL(dp), DIMENSION(:, :), POINTER :: mass, mayer, r_old
1195 : REAL(KIND=dp) :: ave_virial, current_value, distance, exp_max_val, exp_min_val, exponent, &
1196 : integral, previous_value, square_value, trial_energy, triangle_value
1197 : REAL(KIND=dp), DIMENSION(1:3) :: abc, center_of_mass
1198 2 : TYPE(cell_p_type), DIMENSION(:), POINTER :: cell
1199 2 : TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsys
1200 2 : TYPE(force_env_p_type), DIMENSION(:), POINTER :: force_env
1201 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1202 : TYPE(mc_simulation_parameters_p_type), &
1203 2 : DIMENSION(:), POINTER :: mc_par
1204 : TYPE(mp_comm_type) :: group
1205 2 : TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1206 :
1207 : ! these are current magic numbers for how we compute the virial...
1208 : ! we break it up into three parts to integrate the function so provide
1209 : ! better statistics
1210 :
1211 2 : nintegral_divisions = 3
1212 0 : ALLOCATE (virial_cutoffs(1:nintegral_divisions))
1213 2 : ALLOCATE (virial_stepsize(1:nintegral_divisions))
1214 2 : virial_cutoffs(1) = 8.0 ! first distance, in bohr
1215 2 : virial_cutoffs(2) = 13.0 ! second distance, in bohr
1216 2 : virial_cutoffs(3) = 22.0 ! maximum distance, in bohr
1217 2 : virial_stepsize(1) = 0.04 ! stepsize from 0 to virial_cutoffs(1)
1218 2 : virial_stepsize(2) = 0.1
1219 2 : virial_stepsize(3) = 0.2
1220 :
1221 : nbins = CEILING(virial_cutoffs(1)/virial_stepsize(1) + (virial_cutoffs(2) - virial_cutoffs(1))/ &
1222 : virial_stepsize(2) + (virial_cutoffs(3) - virial_cutoffs(2))/virial_stepsize(3))
1223 :
1224 : ! figure out what the default write unit is
1225 2 : iw = cp_logger_get_default_io_unit()
1226 :
1227 : ! allocate a whole bunch of stuff based on how many boxes we have
1228 4 : ALLOCATE (force_env(1:1))
1229 4 : ALLOCATE (cell(1:1))
1230 4 : ALLOCATE (particles(1:1))
1231 4 : ALLOCATE (subsys(1:1))
1232 4 : ALLOCATE (mc_par(1:1))
1233 :
1234 : CALL get_mc_env(mc_env(1)%mc_env, &
1235 : mc_par=mc_par(1)%mc_par, &
1236 2 : force_env=force_env(1)%force_env)
1237 :
1238 : ! get some data out of mc_par
1239 : CALL get_mc_par(mc_par(1)%mc_par, &
1240 : exp_max_val=exp_max_val, &
1241 : exp_min_val=exp_min_val, nvirial=nvirial, &
1242 : ionode=ionode, source=source, group=group, &
1243 2 : mc_molecule_info=mc_molecule_info, virial_temps=virial_temps)
1244 :
1245 2 : IF (iw > 0) THEN
1246 1 : WRITE (iw, *)
1247 1 : WRITE (iw, *)
1248 1 : WRITE (iw, *) 'Beginning the calculation of the second virial coefficient'
1249 1 : WRITE (iw, *)
1250 1 : WRITE (iw, *)
1251 : END IF
1252 :
1253 : ! get some data from the molecule types
1254 : CALL get_mc_molecule_info(mc_molecule_info, &
1255 : nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
1256 : mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
1257 2 : mass=mass)
1258 :
1259 2 : nvirial_temps = SIZE(virial_temps)
1260 6 : ALLOCATE (BETA(1:nvirial_temps))
1261 :
1262 6 : DO itemp = 1, nvirial_temps
1263 6 : BETA(itemp) = 1/virial_temps(itemp)/boltzmann*joule
1264 : END DO
1265 :
1266 : ! get the subsystems and the cell information
1267 : CALL force_env_get(force_env(1)%force_env, &
1268 2 : subsys=subsys(1)%subsys, cell=cell(1)%cell)
1269 2 : CALL get_cell(cell(1)%cell, abc=abc(:))
1270 : CALL cp_subsys_get(subsys(1)%subsys, &
1271 2 : particles=particles(1)%list)
1272 :
1273 : ! check and make sure the box is big enough
1274 2 : IF (abc(1) /= abc(2) .OR. abc(2) /= abc(3)) THEN
1275 0 : CPABORT('The box needs to be cubic for a virial calculation (it is easiest).')
1276 : END IF
1277 2 : IF (virial_cutoffs(nintegral_divisions) > abc(1)/2.0E0_dp) THEN
1278 0 : IF (iw > 0) THEN
1279 0 : WRITE (iw, *) "Box length ", abc(1)*angstrom, " virial cutoff ", &
1280 0 : virial_cutoffs(nintegral_divisions)*angstrom
1281 : END IF
1282 0 : CPABORT('You need a bigger box to deal with this virial cutoff (see above).')
1283 : END IF
1284 :
1285 : ! store the coordinates of the molecules in an array so we can work with it
1286 6 : ALLOCATE (r_old(1:3, 1:nunits_tot(1)))
1287 :
1288 14 : DO iparticle = 1, nunits_tot(1)
1289 : r_old(1:3, iparticle) = &
1290 50 : particles(1)%list%els(iparticle)%r(1:3)
1291 : END DO
1292 :
1293 : ! move the center of mass of molecule 1 to the origin
1294 2 : start_atom = 1
1295 2 : end_atom = nunits(mol_type(1))
1296 : CALL get_center_of_mass(r_old(:, start_atom:end_atom), nunits(mol_type(1)), &
1297 2 : center_of_mass(:), mass(1:nunits(mol_type(1)), mol_type(1)))
1298 8 : DO iunit = start_atom, end_atom
1299 26 : r_old(:, iunit) = r_old(:, iunit) - center_of_mass(:)
1300 : END DO
1301 : ! set them in the force_env, so the first molecule is ready for the energy calc
1302 8 : DO iparticle = start_atom, end_atom
1303 44 : particles(1)%list%els(iparticle)%r(1:3) = r_old(1:3, iparticle)
1304 : END DO
1305 :
1306 : ! print out a notice every 1%
1307 2 : iprint = FLOOR(REAL(nvirial, KIND=dp)/100.0_dp)
1308 2 : IF (iprint == 0) iprint = 1
1309 :
1310 : ! we'll compute the average potential, and then integrate that, as opposed to
1311 : ! integrating every orientation and then averaging
1312 8 : ALLOCATE (mayer(1:nvirial_temps, 1:nbins))
1313 :
1314 1778 : mayer(:, :) = 0.0_dp
1315 :
1316 : ! loop over all nvirial random configurations
1317 22 : DO ivirial = 1, nvirial
1318 :
1319 : ! move molecule two back to the origin
1320 20 : start_atom = nunits(mol_type(1)) + 1
1321 20 : end_atom = nunits_tot(1)
1322 : CALL get_center_of_mass(r_old(:, start_atom:end_atom), nunits(mol_type(2)), &
1323 20 : center_of_mass(:), mass(1:nunits(mol_type(2)), mol_type(2)))
1324 80 : DO iunit = start_atom, end_atom
1325 260 : r_old(:, iunit) = r_old(:, iunit) - center_of_mass(:)
1326 : END DO
1327 :
1328 : ! now we need a random orientation for molecule 2...this routine is
1329 : ! only done in serial since it calls a random number
1330 20 : IF (ionode) THEN
1331 : CALL rotate_molecule(r_old(:, start_atom:end_atom), &
1332 : mass(1:nunits(mol_type(2)), mol_type(2)), &
1333 10 : nunits(mol_type(2)), rng_stream)
1334 : END IF
1335 980 : CALL group%bcast(r_old(:, :), source)
1336 :
1337 20 : distance = 0.0E0_dp
1338 20 : ibin = 1
1339 5900 : DO
1340 : ! find out what our stepsize is
1341 5920 : current_division = 0
1342 8780 : DO idivision = 1, nintegral_divisions
1343 8780 : IF (distance < virial_cutoffs(idivision) - virial_stepsize(idivision)/2.0E0_dp) THEN
1344 : current_division = idivision
1345 : EXIT
1346 : END IF
1347 : END DO
1348 5920 : IF (current_division == 0) EXIT
1349 5900 : distance = distance + virial_stepsize(current_division)
1350 :
1351 : ! move the second molecule only along the x direction
1352 23600 : DO iparticle = start_atom, end_atom
1353 17700 : particles(1)%list%els(iparticle)%r(1) = r_old(1, iparticle) + distance
1354 17700 : particles(1)%list%els(iparticle)%r(2) = r_old(2, iparticle)
1355 23600 : particles(1)%list%els(iparticle)%r(3) = r_old(3, iparticle)
1356 : END DO
1357 :
1358 : ! check for overlaps
1359 5900 : CALL check_for_overlap(force_env(1)%force_env, nchains(:, 1), nunits, loverlap, mol_type)
1360 :
1361 : ! compute the energy if there is no overlap
1362 : ! exponent is exp(-beta*energy)-1, also called the Mayer term
1363 5900 : IF (loverlap) THEN
1364 2202 : DO itemp = 1, nvirial_temps
1365 2202 : mayer(itemp, ibin) = mayer(itemp, ibin) - 1.0_dp
1366 : END DO
1367 : ELSE
1368 : CALL force_env_calc_energy_force(force_env(1)%force_env, &
1369 5166 : calc_force=.FALSE.)
1370 : CALL force_env_get(force_env(1)%force_env, &
1371 5166 : potential_energy=trial_energy)
1372 :
1373 15498 : DO itemp = 1, nvirial_temps
1374 :
1375 10332 : exponent = -BETA(itemp)*trial_energy
1376 :
1377 10332 : IF (exponent > exp_max_val) THEN
1378 : exponent = exp_max_val
1379 10332 : ELSE IF (exponent < exp_min_val) THEN
1380 : exponent = exp_min_val
1381 : END IF
1382 15498 : mayer(itemp, ibin) = mayer(itemp, ibin) + EXP(exponent) - 1.0_dp
1383 : END DO
1384 : END IF
1385 :
1386 5900 : ibin = ibin + 1
1387 : END DO
1388 : ! write out some info that keeps track of where we are
1389 22 : IF (iw > 0) THEN
1390 10 : IF (MOD(ivirial, iprint) == 0) THEN
1391 10 : WRITE (iw, '(A,I6,A,I6)') ' Done with config ', ivirial, ' out of ', nvirial
1392 : END IF
1393 : END IF
1394 : END DO
1395 :
1396 : ! now we integrate this average potential
1397 1778 : mayer(:, :) = mayer(:, :)/REAL(nvirial, dp)
1398 :
1399 6 : DO itemp = 1, nvirial_temps
1400 : integral = 0.0_dp
1401 : previous_value = 0.0_dp
1402 : distance = 0.0E0_dp
1403 : ibin = 1
1404 1180 : DO
1405 1184 : current_division = 0
1406 1756 : DO idivision = 1, nintegral_divisions
1407 1756 : IF (distance < virial_cutoffs(idivision) - virial_stepsize(idivision)/2.0E0_dp) THEN
1408 : current_division = idivision
1409 : EXIT
1410 : END IF
1411 : END DO
1412 1184 : IF (current_division == 0) EXIT
1413 1180 : distance = distance + virial_stepsize(current_division)
1414 :
1415 : ! now we need to integrate, using the trapazoidal method
1416 : ! first, find the value of the square
1417 1180 : current_value = mayer(itemp, ibin)*distance**2
1418 1180 : square_value = previous_value*virial_stepsize(current_division)
1419 : ! now the triangle that sits on top of it, which is half the size of this square...
1420 : ! notice this is negative if the current value is less than the previous value
1421 1180 : triangle_value = 0.5E0_dp*((current_value - previous_value)*virial_stepsize(current_division))
1422 :
1423 1180 : integral = integral + square_value + triangle_value
1424 1180 : previous_value = current_value
1425 1180 : ibin = ibin + 1
1426 : END DO
1427 :
1428 : ! now that the integration is done, compute the second virial that results
1429 4 : ave_virial = -2.0E0_dp*pi*integral
1430 :
1431 : ! convert from CP2K units to something else
1432 4 : ave_virial = ave_virial*n_avogadro*angstrom**3/1.0E8_dp**3
1433 :
1434 6 : IF (iw > 0) THEN
1435 2 : WRITE (iw, *)
1436 2 : WRITE (iw, *) '*********************************************************************'
1437 2 : WRITE (iw, '(A,F12.6,A)') ' *** Temperature = ', virial_temps(itemp), &
1438 4 : ' ***'
1439 2 : WRITE (iw, *) '*** ***'
1440 2 : WRITE (iw, '(A,E12.6,A)') ' *** B2(T) = ', ave_virial, &
1441 4 : ' cm**3/mol ***'
1442 2 : WRITE (iw, *) '*********************************************************************'
1443 2 : WRITE (iw, *)
1444 : END IF
1445 : END DO
1446 :
1447 : ! deallocate some stuff
1448 2 : DEALLOCATE (mc_par)
1449 2 : DEALLOCATE (subsys)
1450 2 : DEALLOCATE (force_env)
1451 2 : DEALLOCATE (particles)
1452 2 : DEALLOCATE (cell)
1453 2 : DEALLOCATE (virial_cutoffs)
1454 2 : DEALLOCATE (virial_stepsize)
1455 2 : DEALLOCATE (r_old)
1456 2 : DEALLOCATE (mayer)
1457 2 : DEALLOCATE (BETA)
1458 :
1459 6 : END SUBROUTINE mc_compute_virial
1460 :
1461 : END MODULE mc_ensembles
1462 :
|
