Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 .GT. 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) .LT. 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) .LT. 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) .LT. 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)) .NE. 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)) .NE. 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 33 : 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 .NE. 1) &
430 0 : CPABORT('ldiscrete=.true. ONLY for systems with 1 box')
431 : CALL create_discrete_array(abc(:, 1), discrete_array(:, :), &
432 0 : discrete_step)
433 : END IF
434 :
435 : ! find out how many steps we're doing...change the updates to be in cycles
436 : ! if the total number of steps is measured in cycles
437 18 : IF (.NOT. lstop) THEN
438 10 : nstep = nstep*nchain_total
439 10 : iuptrans = iuptrans*nchain_total
440 10 : iupvolume = iupvolume*nchain_total
441 : END IF
442 :
443 486 : DO nnstep = nstart + 1, nstart + nstep
444 :
445 468 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
446 15 : WRITE (iw, *)
447 15 : WRITE (iw, *) "------- On Monte Carlo Step ", nnstep
448 : END IF
449 :
450 468 : IF (ionode) rand = rng_stream%next()
451 : ! broadcast the random number, to make sure we're on the same move
452 468 : CALL group%bcast(rand, source)
453 :
454 468 : IF (rand .LT. pmvolume) THEN
455 :
456 58 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
457 1 : WRITE (iw, *) "Attempting a volume move"
458 1 : WRITE (iw, *)
459 : END IF
460 :
461 8 : SELECT CASE (ensemble)
462 : CASE ("TRADITIONAL")
463 : CALL mc_volume_move(mc_par(1)%mc_par, &
464 : force_env(1)%force_env, &
465 : moves(1, 1)%moves, move_updates(1, 1)%moves, &
466 : old_energy(1), 1, &
467 : energy_check(1), r_old(:, :, 1), iw, discrete_array(:, :), &
468 8 : rng_stream)
469 : CASE ("GEMC_NVT")
470 : CALL mc_ge_volume_move(mc_par, force_env, moves, &
471 : move_updates, nnstep, old_energy, energy_check, &
472 24 : r_old, rng_stream)
473 : CASE ("GEMC_NPT")
474 : ! we need to select a box based on the probability given in the input file
475 26 : IF (ionode) rand = rng_stream%next()
476 26 : CALL group%bcast(rand, source)
477 :
478 38 : DO ibox = 1, nboxes
479 38 : IF (rand .LE. pmvol_box(ibox)) THEN
480 26 : box_number = ibox
481 26 : EXIT
482 : END IF
483 : END DO
484 :
485 : CALL mc_volume_move(mc_par(box_number)%mc_par, &
486 : force_env(box_number)%force_env, &
487 : moves(1, box_number)%moves, &
488 : move_updates(1, box_number)%moves, &
489 : old_energy(box_number), box_number, &
490 : energy_check(box_number), r_old(:, :, box_number), iw, &
491 : discrete_array(:, :), &
492 84 : rng_stream)
493 : END SELECT
494 :
495 : ! update all the pointers here, because otherwise we may pass wrong information when we're making a bias environment
496 166 : DO ibox = 1, nboxes
497 : CALL force_env_get(force_env(ibox)%force_env, &
498 108 : subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
499 108 : CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
500 : CALL cp_subsys_get(oldsys(ibox)%subsys, &
501 166 : particles=particles_old(ibox)%list)
502 : END DO
503 :
504 : ! we need a new biasing environment now, if we're into that sort of thing
505 58 : IF (lbias) THEN
506 150 : DO ibox = 1, nboxes
507 100 : CALL force_env_release(bias_env(ibox)%force_env)
508 : ! determine the atom names of every particle
509 300 : ALLOCATE (atom_names_box(1:nunits_tot(ibox)))
510 150 : start_mol = 1
511 150 : DO jbox = 1, ibox - 1
512 250 : start_mol = start_mol + SUM(nchains(:, jbox))
513 : END DO
514 300 : end_mol = start_mol + SUM(nchains(:, ibox)) - 1
515 300 : atom_number = 1
516 1500 : DO imolecule = 1, SUM(nchains(:, ibox))
517 3800 : DO iunit = 1, nunits(mol_type(imolecule + start_mol - 1))
518 : atom_names_box(atom_number) = &
519 2500 : atom_names(iunit, mol_type(imolecule + start_mol - 1))
520 3700 : atom_number = atom_number + 1
521 : END DO
522 : END DO
523 :
524 : ! need to find out what the cell lengths are
525 : CALL force_env_get(force_env(ibox)%force_env, &
526 100 : subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
527 100 : CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
528 :
529 : CALL get_mc_par(mc_par(ibox)%mc_par, &
530 100 : mc_bias_file=mc_bias_file)
531 100 : nchains_box => nchains(:, ibox)
532 :
533 : CALL mc_create_bias_force_env(bias_env(ibox)%force_env, &
534 : r_old(:, :, ibox), atom_names_box(:), nunits_tot(ibox), &
535 : para_env, abc(:, ibox), nchains_box, input_declaration, &
536 100 : mc_bias_file, ionode)
537 :
538 300 : IF (SUM(nchains(:, ibox)) .NE. 0) THEN
539 : CALL force_env_calc_energy_force( &
540 : bias_env(ibox)%force_env, &
541 100 : calc_force=.FALSE.)
542 : CALL force_env_get(bias_env(ibox)%force_env, &
543 100 : potential_energy=last_bias_energy(ibox))
544 : ELSE
545 0 : last_bias_energy(ibox) = 0.0E0_dp
546 : END IF
547 100 : bias_energy(ibox) = last_bias_energy(ibox)
548 150 : DEALLOCATE (atom_names_box)
549 : END DO
550 : END IF
551 :
552 410 : ELSEIF (rand .LT. pmswap) THEN
553 :
554 : ! try a swap move
555 22 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
556 0 : WRITE (iw, *) "Attempting a swap move"
557 0 : WRITE (iw, *)
558 : END IF
559 :
560 : CALL mc_ge_swap_move(mc_par, force_env, bias_env, moves, &
561 : energy_check(:), r_old(:, :, :), old_energy(:), input_declaration, &
562 22 : para_env, bias_energy(:), last_bias_energy(:), rng_stream)
563 :
564 : ! the number of molecules may have changed, which deallocated the whole
565 : ! mc_molecule_info structure
566 22 : CALL get_mc_par(mc_par(1)%mc_par, mc_molecule_info=mc_molecule_info)
567 : CALL get_mc_molecule_info(mc_molecule_info, conf_prob=conf_prob, &
568 : nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
569 : mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
570 22 : atom_names=atom_names, mass=mass)
571 :
572 388 : ELSEIF (rand .LT. pmhmc) THEN
573 : ! try hybrid Monte Carlo
574 20 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
575 2 : WRITE (iw, *) "Attempting a hybrid Monte Carlo move"
576 2 : WRITE (iw, *)
577 : END IF
578 :
579 : ! pick a box at random
580 20 : IF (ionode) rand = rng_stream%next()
581 20 : CALL group%bcast(rand, source)
582 :
583 20 : DO ibox = 1, nboxes
584 20 : IF (rand .LE. pmhmc_box(ibox)) THEN
585 20 : box_number = ibox
586 20 : EXIT
587 : END IF
588 : END DO
589 :
590 : CALL mc_hmc_move(mc_par(box_number)%mc_par, &
591 : force_env(box_number)%force_env, globenv, &
592 : moves(1, box_number)%moves, &
593 : move_updates(1, box_number)%moves, &
594 : old_energy(box_number), box_number, &
595 : energy_check(box_number), r_old(:, :, box_number), &
596 20 : rng_stream)
597 :
598 368 : ELSEIF (rand .LT. pmavbmc) THEN
599 : ! try an AVBMC move
600 0 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
601 0 : WRITE (iw, *) "Attempting an AVBMC1 move"
602 0 : WRITE (iw, *)
603 : END IF
604 :
605 : ! first, pick a box to do it for
606 0 : IF (ionode) rand = rng_stream%next()
607 0 : CALL group%bcast(rand, source)
608 :
609 0 : IF (nboxes .EQ. 2) THEN
610 0 : IF (rand .LT. 0.1E0_dp) THEN
611 0 : box_number = 1
612 : ELSE
613 0 : box_number = 2
614 : END IF
615 : ELSE
616 0 : box_number = 1
617 : END IF
618 :
619 : ! now pick a molecule type to do it for
620 0 : IF (ionode) rand = rng_stream%next()
621 0 : CALL group%bcast(rand, source)
622 0 : molecule_type_swap = 0
623 0 : DO imol_type = 1, nmol_types
624 0 : IF (rand .LT. pmavbmc_mol(imol_type)) THEN
625 0 : molecule_type_swap = imol_type
626 0 : EXIT
627 : END IF
628 : END DO
629 0 : IF (molecule_type_swap == 0) &
630 0 : CPABORT('Did not choose a molecule type to swap...check AVBMC input')
631 :
632 : ! now pick a molecule, automatically rejecting the move if the
633 : ! box is empty or only has one molecule
634 0 : IF (SUM(nchains(:, box_number)) .LE. 1) THEN
635 : ! indicate that we tried a move
636 : moves(molecule_type_swap, box_number)%moves%empty_avbmc = &
637 0 : moves(molecule_type_swap, box_number)%moves%empty_avbmc + 1
638 : ELSE
639 :
640 : ! pick a molecule to be swapped in the box
641 0 : IF (ionode) THEN
642 : CALL find_mc_test_molecule(mc_molecule_info, &
643 : start_atom_swap, idum, jdum, rng_stream, &
644 0 : box=box_number, molecule_type_old=molecule_type_swap)
645 :
646 : ! pick a molecule to act as the target in the box...we don't care what type
647 0 : DO
648 : CALL find_mc_test_molecule(mc_molecule_info, &
649 : start_atom_target, idum, molecule_type_target, &
650 0 : rng_stream, box=box_number)
651 0 : IF (start_atom_swap .NE. start_atom_target) THEN
652 : start_atom_target = start_atom_target + &
653 0 : avbmc_atom(molecule_type_target) - 1
654 : EXIT
655 : END IF
656 : END DO
657 :
658 : ! choose if we're swapping into the bonded region of mol_target, or
659 : ! into the nonbonded region
660 0 : rand = rng_stream%next()
661 :
662 : END IF
663 0 : CALL group%bcast(start_atom_swap, source)
664 0 : CALL group%bcast(box_number, source)
665 0 : CALL group%bcast(start_atom_target, source)
666 0 : CALL group%bcast(rand, source)
667 :
668 0 : IF (rand .LT. pbias(molecule_type_swap)) THEN
669 0 : move_type_avbmc = 'in'
670 : ELSE
671 0 : move_type_avbmc = 'out'
672 : END IF
673 :
674 : CALL mc_avbmc_move(mc_par(box_number)%mc_par, &
675 : force_env(box_number)%force_env, &
676 : bias_env(box_number)%force_env, &
677 : moves(molecule_type_swap, box_number)%moves, &
678 : energy_check(box_number), &
679 : r_old(:, :, box_number), old_energy(box_number), &
680 : start_atom_swap, start_atom_target, molecule_type_swap, &
681 : box_number, bias_energy(box_number), &
682 : last_bias_energy(box_number), &
683 0 : move_type_avbmc, rng_stream)
684 :
685 : END IF
686 :
687 : ELSE
688 :
689 368 : IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
690 12 : WRITE (iw, *) "Attempting an inner move"
691 12 : WRITE (iw, *)
692 : END IF
693 :
694 2010 : DO imove = 1, nmoves
695 :
696 1642 : IF (ionode) rand = rng_stream%next()
697 1642 : CALL group%bcast(rand, source)
698 2010 : IF (rand .LT. pmtraion) THEN
699 : ! change molecular conformation
700 : ! first, pick a box to do it for
701 506 : IF (ionode) rand = rng_stream%next()
702 506 : CALL group%bcast(rand, source)
703 506 : IF (nboxes .EQ. 2) THEN
704 0 : IF (rand .LT. 0.75E0_dp) THEN
705 0 : box_number = 1
706 : ELSE
707 0 : box_number = 2
708 : END IF
709 : ELSE
710 506 : box_number = 1
711 : END IF
712 :
713 : ! figure out which molecule type we're looking for
714 506 : IF (ionode) rand = rng_stream%next()
715 506 : CALL group%bcast(rand, source)
716 506 : molecule_type = 0
717 506 : DO imol_type = 1, nmol_types
718 506 : IF (rand .LT. pmtraion_mol(imol_type)) THEN
719 506 : molecule_type = imol_type
720 506 : EXIT
721 : END IF
722 : END DO
723 506 : IF (molecule_type == 0) CALL cp_abort( &
724 : __LOCATION__, &
725 0 : 'Did not choose a molecule type to conf change...PMTRAION_MOL should not be all 0.0')
726 :
727 : ! now pick a molecule, automatically rejecting the move if the
728 : ! box is empty
729 506 : IF (nchains(molecule_type, box_number) == 0) THEN
730 : ! indicate that we tried a move
731 : moves(molecule_type, box_number)%moves%empty_conf = &
732 0 : moves(molecule_type, box_number)%moves%empty_conf + 1
733 : ELSE
734 : ! pick a molecule in the box
735 506 : IF (ionode) THEN
736 : CALL find_mc_test_molecule(mc_molecule_info, &
737 : start_atom, idum, &
738 : jdum, rng_stream, &
739 253 : box=box_number, molecule_type_old=molecule_type)
740 :
741 : ! choose if we're changing a bond length or an angle
742 253 : rand = rng_stream%next()
743 : END IF
744 506 : CALL group%bcast(rand, source)
745 506 : CALL group%bcast(start_atom, source)
746 506 : CALL group%bcast(box_number, source)
747 506 : CALL group%bcast(molecule_type, source)
748 :
749 : ! figure out what kind of move we're doing
750 506 : IF (rand .LT. conf_prob(1, molecule_type)) THEN
751 312 : move_type = 'bond'
752 194 : ELSEIF (rand .LT. (conf_prob(1, molecule_type) + &
753 : conf_prob(2, molecule_type))) THEN
754 194 : move_type = 'angle'
755 : ELSE
756 0 : move_type = 'dihedral'
757 : END IF
758 506 : box_flag(box_number) = 1
759 : CALL mc_conformation_change(mc_par(box_number)%mc_par, &
760 : force_env(box_number)%force_env, &
761 : bias_env(box_number)%force_env, &
762 : moves(molecule_type, box_number)%moves, &
763 : move_updates(molecule_type, box_number)%moves, &
764 : start_atom, molecule_type, box_number, &
765 : bias_energy(box_number), &
766 506 : move_type, lreject, rng_stream)
767 506 : IF (lreject) EXIT
768 : END IF
769 1136 : ELSEIF (rand .LT. pmtrans) THEN
770 : ! translate a whole molecule in the system
771 : ! pick a molecule type
772 624 : IF (ionode) rand = rng_stream%next()
773 624 : CALL group%bcast(rand, source)
774 624 : molecule_type = 0
775 924 : DO imol_type = 1, nmol_types
776 924 : IF (rand .LT. pmtrans_mol(imol_type)) THEN
777 624 : molecule_type = imol_type
778 624 : EXIT
779 : END IF
780 : END DO
781 624 : IF (molecule_type == 0) CALL cp_abort( &
782 : __LOCATION__, &
783 0 : 'Did not choose a molecule type to translate...PMTRANS_MOL should not be all 0.0')
784 :
785 : ! now pick a molecule of that type
786 624 : IF (ionode) &
787 : CALL find_mc_test_molecule(mc_molecule_info, &
788 : start_atom, box_number, idum, rng_stream, &
789 312 : molecule_type_old=molecule_type)
790 624 : CALL group%bcast(start_atom, source)
791 624 : CALL group%bcast(box_number, source)
792 624 : box_flag(box_number) = 1
793 : CALL mc_molecule_translation(mc_par(box_number)%mc_par, &
794 : force_env(box_number)%force_env, &
795 : bias_env(box_number)%force_env, &
796 : moves(molecule_type, box_number)%moves, &
797 : move_updates(molecule_type, box_number)%moves, &
798 : start_atom, box_number, bias_energy(box_number), &
799 624 : molecule_type, lreject, rng_stream)
800 624 : IF (lreject) EXIT
801 512 : ELSEIF (rand .LT. pmcltrans) THEN
802 : ! translate a whole cluster in the system
803 : ! first, pick a box to do it for
804 10 : IF (ionode) rand = rng_stream%next()
805 10 : CALL group%bcast(rand, source)
806 :
807 10 : DO ibox = 1, nboxes
808 10 : IF (rand .LE. pmclus_box(ibox)) THEN
809 10 : box_number = ibox
810 10 : EXIT
811 : END IF
812 : END DO
813 10 : box_flag(box_number) = 1
814 : CALL mc_cluster_translation(mc_par(box_number)%mc_par, &
815 : force_env(box_number)%force_env, &
816 : bias_env(box_number)%force_env, &
817 : moves(1, box_number)%moves, &
818 : move_updates(1, box_number)%moves, &
819 : box_number, bias_energy(box_number), &
820 10 : lreject, rng_stream)
821 10 : IF (lreject) EXIT
822 : ELSE
823 : ! rotate a whole molecule in the system
824 : ! pick a molecule type
825 502 : IF (ionode) rand = rng_stream%next()
826 502 : CALL group%bcast(rand, source)
827 502 : molecule_type = 0
828 502 : DO imol_type = 1, nmol_types
829 502 : IF (rand .LT. pmrot_mol(imol_type)) THEN
830 502 : molecule_type = imol_type
831 502 : EXIT
832 : END IF
833 : END DO
834 502 : IF (molecule_type == 0) CALL cp_abort( &
835 : __LOCATION__, &
836 0 : 'Did not choose a molecule type to rotate...PMROT_MOL should not be all 0.0')
837 :
838 502 : IF (ionode) &
839 : CALL find_mc_test_molecule(mc_molecule_info, &
840 : start_atom, box_number, idum, rng_stream, &
841 251 : molecule_type_old=molecule_type)
842 502 : CALL group%bcast(start_atom, source)
843 502 : CALL group%bcast(box_number, source)
844 502 : box_flag(box_number) = 1
845 : CALL mc_molecule_rotation(mc_par(box_number)%mc_par, &
846 : force_env(box_number)%force_env, &
847 : bias_env(box_number)%force_env, &
848 : moves(molecule_type, box_number)%moves, &
849 : move_updates(molecule_type, box_number)%moves, &
850 : box_number, start_atom, &
851 : molecule_type, bias_energy(box_number), &
852 502 : lreject, rng_stream)
853 502 : IF (lreject) EXIT
854 : END IF
855 :
856 : END DO
857 :
858 : ! now do a Quickstep calculation to see if we accept the sequence
859 : CALL mc_Quickstep_move(mc_par, force_env, bias_env, &
860 : moves, lreject, move_updates, energy_check(:), r_old(:, :, :), &
861 : nnstep, old_energy(:), bias_energy(:), last_bias_energy(:), &
862 : nboxes, box_flag(:), oldsys, particles_old, &
863 368 : rng_stream, unit_conv)
864 :
865 : END IF
866 :
867 : ! make sure the pointers are pointing correctly since the subsys may
868 : ! have changed
869 1080 : DO ibox = 1, nboxes
870 : CALL force_env_get(force_env(ibox)%force_env, &
871 612 : subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
872 612 : CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
873 : CALL cp_subsys_get(oldsys(ibox)%subsys, &
874 1080 : particles=particles_old(ibox)%list)
875 : END DO
876 :
877 468 : IF (ionode) THEN
878 :
879 234 : IF (MOD(nnstep, iprint) == 0) THEN
880 15 : WRITE (com_ene, *) nnstep, old_energy(1:nboxes)
881 :
882 33 : DO ibox = 1, nboxes
883 :
884 : ! write the molecule information
885 18 : WRITE (com_mol, *) nnstep, nchains(:, ibox)
886 :
887 : ! write the move statistics to file
888 47 : DO itype = 1, nmol_types
889 : CALL write_move_stats(moves(itype, ibox)%moves, &
890 47 : nnstep, move_unit(ibox))
891 : END DO
892 :
893 : ! write a restart file
894 : CALL write_mc_restart(nnstep, mc_par(ibox)%mc_par, &
895 18 : nchains(:, ibox), force_env(ibox)%force_env)
896 :
897 : ! write cell lengths
898 72 : WRITE (cell_unit, *) nnstep, abc(1:3, ibox)*angstrom
899 :
900 : ! write particle coordinates
901 18 : WRITE (cbox, '(I4)') ibox
902 18 : WRITE (cstep, '(I8)') nnstep
903 62 : IF (SUM(nchains(:, ibox)) == 0) THEN
904 0 : WRITE (com_crd, *) ' 0'
905 : WRITE (com_crd, *) 'BOX '//TRIM(ADJUSTL(cbox))// &
906 0 : ', STEP '//TRIM(ADJUSTL(cstep))
907 : ELSE
908 : CALL write_particle_coordinates( &
909 : particles_old(ibox)%list%els, &
910 : com_crd, dump_xmol, 'POS', &
911 : 'BOX '//TRIM(ADJUSTL(cbox))// &
912 : ', STEP '//TRIM(ADJUSTL(cstep)), &
913 18 : unit_conv=unit_conv)
914 : END IF
915 : END DO
916 : END IF ! end the things we only do every iprint moves
917 :
918 540 : DO ibox = 1, nboxes
919 : ! compute some averages
920 : averages(ibox)%averages%ave_energy = &
921 : averages(ibox)%averages%ave_energy*REAL(nnstep - &
922 : nstart - 1, dp)/REAL(nnstep - nstart, dp) + &
923 306 : old_energy(ibox)/REAL(nnstep - nstart, dp)
924 : averages(ibox)%averages%molecules = &
925 : averages(ibox)%averages%molecules*REAL(nnstep - &
926 : nstart - 1, dp)/REAL(nnstep - nstart, dp) + &
927 899 : REAL(SUM(nchains(:, ibox)), dp)/REAL(nnstep - nstart, dp)
928 : averages(ibox)%averages%ave_volume = &
929 : averages(ibox)%averages%ave_volume* &
930 : REAL(nnstep - nstart - 1, dp)/REAL(nnstep - nstart, dp) + &
931 : abc(1, ibox)*abc(2, ibox)*abc(3, ibox)/ &
932 306 : REAL(nnstep - nstart, dp)
933 :
934 : ! flush the buffers to the files
935 306 : CALL m_flush(data_unit(ibox))
936 306 : CALL m_flush(diff(ibox))
937 306 : CALL m_flush(move_unit(ibox))
938 306 : CALL m_flush(cl(ibox))
939 540 : CALL m_flush(rm(ibox))
940 :
941 : END DO
942 :
943 : ! flush more buffers to the files
944 234 : CALL m_flush(cell_unit)
945 234 : CALL m_flush(com_ene)
946 234 : CALL m_flush(com_crd)
947 234 : CALL m_flush(com_mol)
948 :
949 : END IF
950 :
951 : ! reset the box flags
952 1080 : box_flag(:) = 0
953 :
954 : ! check to see if EXIT file exists...if so, end the calculation
955 468 : CALL external_control(should_stop, "MC", globenv=globenv)
956 468 : IF (should_stop) EXIT
957 :
958 : ! update the move displacements, if necessary
959 1080 : DO ibox = 1, nboxes
960 612 : IF (MOD(nnstep - nstart, iuptrans) == 0) THEN
961 0 : DO itype = 1, nmol_types
962 : CALL mc_move_update(mc_par(ibox)%mc_par, &
963 : move_updates(itype, ibox)%moves, itype, &
964 0 : "trans", nnstep, ionode)
965 : END DO
966 : END IF
967 :
968 1080 : IF (MOD(nnstep - nstart, iupvolume) == 0) THEN
969 : CALL mc_move_update(mc_par(ibox)%mc_par, &
970 : move_updates(1, ibox)%moves, 1337, &
971 0 : "volume", nnstep, ionode)
972 : END IF
973 : END DO
974 :
975 : ! check to see if there are any overlaps in the boxes, and fold coordinates
976 : ! don't care about overlaps if we're only doing HMC
977 468 : IF (.NOT. lhmc) THEN
978 1040 : DO ibox = 1, nboxes
979 2206 : IF (SUM(nchains(:, ibox)) .NE. 0) THEN
980 : start_mol = 1
981 736 : DO jbox = 1, ibox - 1
982 1024 : start_mol = start_mol + SUM(nchains(:, jbox))
983 : END DO
984 1758 : end_mol = start_mol + SUM(nchains(:, ibox)) - 1
985 : CALL check_for_overlap(force_env(ibox)%force_env, &
986 : nchains(:, ibox), nunits, loverlap, &
987 592 : mol_type(start_mol:end_mol))
988 592 : IF (loverlap) THEN
989 0 : IF (iw > 0) WRITE (iw, *) nnstep
990 0 : CPABORT('coordinate overlap at the end of the above step')
991 : ! now fold the coordinates...don't do this anywhere but here, because
992 : ! we can get screwed up with the mc_molecule_info stuff (like in swap move)...
993 : ! this is kind of ugly, with allocated and deallocating every time
994 0 : ALLOCATE (r_temp(1:3, 1:nunits_tot(ibox)))
995 :
996 0 : DO iunit = 1, nunits_tot(ibox)
997 : r_temp(1:3, iunit) = &
998 0 : particles_old(ibox)%list%els(iunit)%r(1:3)
999 : END DO
1000 :
1001 : CALL mc_coordinate_fold(r_temp(:, :), &
1002 : SUM(nchains(:, ibox)), mol_type(start_mol:end_mol), &
1003 0 : mass, nunits, abc(1:3, ibox))
1004 :
1005 : ! save the folded coordinates
1006 0 : DO iunit = 1, nunits_tot(ibox)
1007 0 : r_old(1:3, iunit, ibox) = r_temp(1:3, iunit)
1008 : particles_old(ibox)%list%els(iunit)%r(1:3) = &
1009 0 : r_temp(1:3, iunit)
1010 : END DO
1011 :
1012 : ! if we're biasing, we need to do the same
1013 0 : IF (lbias) THEN
1014 : CALL force_env_get(bias_env(ibox)%force_env, &
1015 0 : subsys=biassys)
1016 : CALL cp_subsys_get(biassys, &
1017 0 : particles=particles_bias)
1018 :
1019 0 : DO iunit = 1, nunits_tot(ibox)
1020 : particles_bias%els(iunit)%r(1:3) = &
1021 0 : r_temp(1:3, iunit)
1022 : END DO
1023 : END IF
1024 :
1025 0 : DEALLOCATE (r_temp)
1026 : END IF
1027 : END IF
1028 : END DO
1029 : END IF
1030 :
1031 : !debug code
1032 486 : IF (debug_this_module) THEN
1033 : DO ibox = 1, nboxes
1034 : IF (SUM(nchains(:, ibox)) .NE. 0) THEN
1035 : CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
1036 : calc_force=.FALSE.)
1037 : CALL force_env_get(force_env(ibox)%force_env, &
1038 : potential_energy=test_energy)
1039 : ELSE
1040 : test_energy = 0.0E0_dp
1041 : END IF
1042 :
1043 : IF (ABS(initial_energy(ibox) + energy_check(ibox) - &
1044 : test_energy) .GT. 0.0000001E0_dp) THEN
1045 : IF (iw > 0) THEN
1046 : WRITE (iw, *) '!!!!!!! We have an energy problem. !!!!!!!!'
1047 : WRITE (iw, '(A,T64,F16.10)') 'Final Energy = ', test_energy
1048 : WRITE (iw, '(A,T64,F16.10)') 'Initial Energy+energy_check=', &
1049 : initial_energy(ibox) + energy_check(ibox)
1050 : WRITE (iw, *) 'Box ', ibox
1051 : WRITE (iw, *) 'nchains ', nchains(:, ibox)
1052 : END IF
1053 : CPABORT('!!!!!!! We have an energy problem. !!!!!!!!')
1054 : END IF
1055 : END DO
1056 : END IF
1057 : END DO
1058 :
1059 : ! write a restart file
1060 18 : IF (ionode) THEN
1061 21 : DO ibox = 1, nboxes
1062 : CALL write_mc_restart(nnstep, mc_par(ibox)%mc_par, &
1063 21 : nchains(:, ibox), force_env(ibox)%force_env)
1064 : END DO
1065 : END IF
1066 :
1067 : ! calculate the final energy
1068 42 : DO ibox = 1, nboxes
1069 64 : IF (SUM(nchains(:, ibox)) .NE. 0) THEN
1070 : CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
1071 24 : calc_force=.FALSE.)
1072 : CALL force_env_get(force_env(ibox)%force_env, &
1073 24 : potential_energy=final_energy(ibox))
1074 : ELSE
1075 0 : final_energy(ibox) = 0.0E0_dp
1076 : END IF
1077 42 : IF (lbias) THEN
1078 14 : CALL force_env_release(bias_env(ibox)%force_env)
1079 : END IF
1080 : END DO
1081 :
1082 : ! do some stuff in serial
1083 18 : IF (ionode .OR. (iw > 0)) THEN
1084 :
1085 9 : WRITE (com_ene, *) 'Final Energies: ', &
1086 18 : final_energy(1:nboxes)
1087 :
1088 21 : DO ibox = 1, nboxes
1089 12 : WRITE (cbox, '(I4)') ibox
1090 32 : IF (SUM(nchains(:, ibox)) == 0) THEN
1091 0 : WRITE (com_crd, *) ' 0'
1092 0 : WRITE (com_crd, *) 'BOX '//TRIM(ADJUSTL(cbox))
1093 : ELSE
1094 : CALL write_particle_coordinates( &
1095 : particles_old(ibox)%list%els, &
1096 : com_crd, dump_xmol, 'POS', &
1097 12 : 'FINAL BOX '//TRIM(ADJUSTL(cbox)), unit_conv=unit_conv)
1098 : END IF
1099 :
1100 : ! write a bunch of data to the screen
1101 : WRITE (iw, '(A)') &
1102 12 : '------------------------------------------------'
1103 : WRITE (iw, '(A,I1,A)') &
1104 12 : '| BOX ', ibox, &
1105 24 : ' |'
1106 : WRITE (iw, '(A)') &
1107 12 : '------------------------------------------------'
1108 12 : test_moves => moves(:, ibox)
1109 : CALL final_mc_write(mc_par(ibox)%mc_par, test_moves, &
1110 : iw, energy_check(ibox), &
1111 : initial_energy(ibox), final_energy(ibox), &
1112 12 : averages(ibox)%averages)
1113 :
1114 : ! close any open files
1115 12 : CALL close_file(unit_number=diff(ibox))
1116 12 : CALL close_file(unit_number=data_unit(ibox))
1117 12 : CALL close_file(unit_number=move_unit(ibox))
1118 12 : CALL close_file(unit_number=cl(ibox))
1119 21 : CALL close_file(unit_number=rm(ibox))
1120 : END DO
1121 :
1122 : ! close some more files
1123 9 : CALL close_file(unit_number=cell_unit)
1124 9 : CALL close_file(unit_number=com_ene)
1125 9 : CALL close_file(unit_number=com_crd)
1126 9 : CALL close_file(unit_number=com_mol)
1127 : END IF
1128 :
1129 42 : DO ibox = 1, nboxes
1130 : CALL set_mc_env(mc_env(ibox)%mc_env, &
1131 : mc_par=mc_par(ibox)%mc_par, &
1132 24 : force_env=force_env(ibox)%force_env)
1133 :
1134 : ! deallocate some stuff
1135 64 : DO itype = 1, nmol_types
1136 40 : CALL mc_moves_release(move_updates(itype, ibox)%moves)
1137 64 : CALL mc_moves_release(moves(itype, ibox)%moves)
1138 : END DO
1139 42 : CALL mc_averages_release(averages(ibox)%averages)
1140 : END DO
1141 :
1142 18 : DEALLOCATE (pmhmc_box)
1143 18 : DEALLOCATE (pmvol_box)
1144 18 : DEALLOCATE (pmclus_box)
1145 18 : DEALLOCATE (r_old)
1146 18 : DEALLOCATE (force_env)
1147 18 : DEALLOCATE (bias_env)
1148 18 : DEALLOCATE (cell)
1149 18 : DEALLOCATE (particles_old)
1150 18 : DEALLOCATE (oldsys)
1151 18 : DEALLOCATE (averages)
1152 18 : DEALLOCATE (moves)
1153 18 : DEALLOCATE (move_updates)
1154 18 : DEALLOCATE (mc_par)
1155 :
1156 : ! end the timing
1157 18 : CALL timestop(handle)
1158 :
1159 54 : END SUBROUTINE mc_run_ensemble
1160 :
1161 : ! **************************************************************************************************
1162 : !> \brief Computes the second virial coefficient of a molecule by using the integral form
1163 : !> of the second virial coefficient found in McQuarrie "Statistical Thermodynamics",
1164 : !> B2(T) = -2Pi Int 0toInf [ Exp[-beta*u(r)] -1] r^2 dr Eq. 15-25
1165 : !> I use trapazoidal integration with various step sizes
1166 : !> (the integral is broken up into three parts, currently, but that's easily
1167 : !> changed by the first variables found below). It generates nvirial configurations,
1168 : !> doing the integration for each one, and then averages all the B2(T) to produce
1169 : !> the final answer.
1170 : !> \param mc_env a pointer that contains all mc_env for all the simulation
1171 : !> boxes
1172 : !> \param rng_stream the stream we pull random numbers from
1173 : !>
1174 : !> Suitable for parallel.
1175 : !> \author MJM
1176 : ! **************************************************************************************************
1177 2 : SUBROUTINE mc_compute_virial(mc_env, rng_stream)
1178 :
1179 : TYPE(mc_environment_p_type), DIMENSION(:), POINTER :: mc_env
1180 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1181 :
1182 : INTEGER :: current_division, end_atom, ibin, idivision, iparticle, iprint, itemp, iunit, &
1183 : ivirial, iw, nbins, nchain_total, nintegral_divisions, nmol_types, nvirial, &
1184 : nvirial_temps, source, start_atom
1185 2 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot
1186 2 : INTEGER, DIMENSION(:, :), POINTER :: nchains
1187 : LOGICAL :: ionode, loverlap
1188 2 : REAL(dp), DIMENSION(:), POINTER :: BETA, virial_cutoffs, virial_stepsize, &
1189 2 : virial_temps
1190 2 : REAL(dp), DIMENSION(:, :), POINTER :: mass, mayer, r_old
1191 : REAL(KIND=dp) :: ave_virial, current_value, distance, exp_max_val, exp_min_val, exponent, &
1192 : integral, previous_value, square_value, trial_energy, triangle_value
1193 : REAL(KIND=dp), DIMENSION(1:3) :: abc, center_of_mass
1194 2 : TYPE(cell_p_type), DIMENSION(:), POINTER :: cell
1195 2 : TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsys
1196 2 : TYPE(force_env_p_type), DIMENSION(:), POINTER :: force_env
1197 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
1198 : TYPE(mc_simulation_parameters_p_type), &
1199 2 : DIMENSION(:), POINTER :: mc_par
1200 : TYPE(mp_comm_type) :: group
1201 2 : TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
1202 :
1203 : ! these are current magic numbers for how we compute the virial...
1204 : ! we break it up into three parts to integrate the function so provide
1205 : ! better statistics
1206 :
1207 2 : nintegral_divisions = 3
1208 0 : ALLOCATE (virial_cutoffs(1:nintegral_divisions))
1209 2 : ALLOCATE (virial_stepsize(1:nintegral_divisions))
1210 2 : virial_cutoffs(1) = 8.0 ! first distance, in bohr
1211 2 : virial_cutoffs(2) = 13.0 ! second distance, in bohr
1212 2 : virial_cutoffs(3) = 22.0 ! maximum distance, in bohr
1213 2 : virial_stepsize(1) = 0.04 ! stepsize from 0 to virial_cutoffs(1)
1214 2 : virial_stepsize(2) = 0.1
1215 2 : virial_stepsize(3) = 0.2
1216 :
1217 : nbins = CEILING(virial_cutoffs(1)/virial_stepsize(1) + (virial_cutoffs(2) - virial_cutoffs(1))/ &
1218 2 : virial_stepsize(2) + (virial_cutoffs(3) - virial_cutoffs(2))/virial_stepsize(3))
1219 :
1220 : ! figure out what the default write unit is
1221 2 : iw = cp_logger_get_default_io_unit()
1222 :
1223 : ! allocate a whole bunch of stuff based on how many boxes we have
1224 4 : ALLOCATE (force_env(1:1))
1225 4 : ALLOCATE (cell(1:1))
1226 4 : ALLOCATE (particles(1:1))
1227 4 : ALLOCATE (subsys(1:1))
1228 4 : ALLOCATE (mc_par(1:1))
1229 :
1230 : CALL get_mc_env(mc_env(1)%mc_env, &
1231 : mc_par=mc_par(1)%mc_par, &
1232 2 : force_env=force_env(1)%force_env)
1233 :
1234 : ! get some data out of mc_par
1235 : CALL get_mc_par(mc_par(1)%mc_par, &
1236 : exp_max_val=exp_max_val, &
1237 : exp_min_val=exp_min_val, nvirial=nvirial, &
1238 : ionode=ionode, source=source, group=group, &
1239 2 : mc_molecule_info=mc_molecule_info, virial_temps=virial_temps)
1240 :
1241 2 : IF (iw > 0) THEN
1242 1 : WRITE (iw, *)
1243 1 : WRITE (iw, *)
1244 1 : WRITE (iw, *) 'Beginning the calculation of the second virial coefficient'
1245 1 : WRITE (iw, *)
1246 1 : WRITE (iw, *)
1247 : END IF
1248 :
1249 : ! get some data from the molecule types
1250 : CALL get_mc_molecule_info(mc_molecule_info, &
1251 : nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
1252 : mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
1253 2 : mass=mass)
1254 :
1255 2 : nvirial_temps = SIZE(virial_temps)
1256 6 : ALLOCATE (BETA(1:nvirial_temps))
1257 :
1258 6 : DO itemp = 1, nvirial_temps
1259 6 : BETA(itemp) = 1/virial_temps(itemp)/boltzmann*joule
1260 : END DO
1261 :
1262 : ! get the subsystems and the cell information
1263 : CALL force_env_get(force_env(1)%force_env, &
1264 2 : subsys=subsys(1)%subsys, cell=cell(1)%cell)
1265 2 : CALL get_cell(cell(1)%cell, abc=abc(:))
1266 : CALL cp_subsys_get(subsys(1)%subsys, &
1267 2 : particles=particles(1)%list)
1268 :
1269 : ! check and make sure the box is big enough
1270 2 : IF (abc(1) .NE. abc(2) .OR. abc(2) .NE. abc(3)) THEN
1271 0 : CPABORT('The box needs to be cubic for a virial calculation (it is easiest).')
1272 : END IF
1273 2 : IF (virial_cutoffs(nintegral_divisions) .GT. abc(1)/2.0E0_dp) THEN
1274 0 : IF (iw > 0) &
1275 0 : WRITE (iw, *) "Box length ", abc(1)*angstrom, " virial cutoff ", &
1276 0 : virial_cutoffs(nintegral_divisions)*angstrom
1277 0 : CPABORT('You need a bigger box to deal with this virial cutoff (see above).')
1278 : END IF
1279 :
1280 : ! store the coordinates of the molecules in an array so we can work with it
1281 6 : ALLOCATE (r_old(1:3, 1:nunits_tot(1)))
1282 :
1283 14 : DO iparticle = 1, nunits_tot(1)
1284 : r_old(1:3, iparticle) = &
1285 50 : particles(1)%list%els(iparticle)%r(1:3)
1286 : END DO
1287 :
1288 : ! move the center of mass of molecule 1 to the origin
1289 2 : start_atom = 1
1290 2 : end_atom = nunits(mol_type(1))
1291 : CALL get_center_of_mass(r_old(:, start_atom:end_atom), nunits(mol_type(1)), &
1292 2 : center_of_mass(:), mass(1:nunits(mol_type(1)), mol_type(1)))
1293 8 : DO iunit = start_atom, end_atom
1294 26 : r_old(:, iunit) = r_old(:, iunit) - center_of_mass(:)
1295 : END DO
1296 : ! set them in the force_env, so the first molecule is ready for the energy calc
1297 8 : DO iparticle = start_atom, end_atom
1298 44 : particles(1)%list%els(iparticle)%r(1:3) = r_old(1:3, iparticle)
1299 : END DO
1300 :
1301 : ! print out a notice every 1%
1302 2 : iprint = FLOOR(REAL(nvirial, KIND=dp)/100.0_dp)
1303 2 : IF (iprint == 0) iprint = 1
1304 :
1305 : ! we'll compute the average potential, and then integrate that, as opposed to
1306 : ! integrating every orientation and then averaging
1307 8 : ALLOCATE (mayer(1:nvirial_temps, 1:nbins))
1308 :
1309 1778 : mayer(:, :) = 0.0_dp
1310 :
1311 : ! loop over all nvirial random configurations
1312 22 : DO ivirial = 1, nvirial
1313 :
1314 : ! move molecule two back to the origin
1315 20 : start_atom = nunits(mol_type(1)) + 1
1316 20 : end_atom = nunits_tot(1)
1317 : CALL get_center_of_mass(r_old(:, start_atom:end_atom), nunits(mol_type(2)), &
1318 20 : center_of_mass(:), mass(1:nunits(mol_type(2)), mol_type(2)))
1319 80 : DO iunit = start_atom, end_atom
1320 260 : r_old(:, iunit) = r_old(:, iunit) - center_of_mass(:)
1321 : END DO
1322 :
1323 : ! now we need a random orientation for molecule 2...this routine is
1324 : ! only done in serial since it calls a random number
1325 20 : IF (ionode) THEN
1326 : CALL rotate_molecule(r_old(:, start_atom:end_atom), &
1327 : mass(1:nunits(mol_type(2)), mol_type(2)), &
1328 10 : nunits(mol_type(2)), rng_stream)
1329 : END IF
1330 980 : CALL group%bcast(r_old(:, :), source)
1331 :
1332 20 : distance = 0.0E0_dp
1333 20 : ibin = 1
1334 5900 : DO
1335 : ! find out what our stepsize is
1336 5920 : current_division = 0
1337 8780 : DO idivision = 1, nintegral_divisions
1338 8780 : IF (distance .LT. virial_cutoffs(idivision) - virial_stepsize(idivision)/2.0E0_dp) THEN
1339 : current_division = idivision
1340 : EXIT
1341 : END IF
1342 : END DO
1343 5920 : IF (current_division == 0) EXIT
1344 5900 : distance = distance + virial_stepsize(current_division)
1345 :
1346 : ! move the second molecule only along the x direction
1347 23600 : DO iparticle = start_atom, end_atom
1348 17700 : particles(1)%list%els(iparticle)%r(1) = r_old(1, iparticle) + distance
1349 17700 : particles(1)%list%els(iparticle)%r(2) = r_old(2, iparticle)
1350 23600 : particles(1)%list%els(iparticle)%r(3) = r_old(3, iparticle)
1351 : END DO
1352 :
1353 : ! check for overlaps
1354 5900 : CALL check_for_overlap(force_env(1)%force_env, nchains(:, 1), nunits, loverlap, mol_type)
1355 :
1356 : ! compute the energy if there is no overlap
1357 : ! exponent is exp(-beta*energy)-1, also called the Mayer term
1358 5900 : IF (loverlap) THEN
1359 3834 : DO itemp = 1, nvirial_temps
1360 3834 : mayer(itemp, ibin) = mayer(itemp, ibin) - 1.0_dp
1361 : END DO
1362 : ELSE
1363 : CALL force_env_calc_energy_force(force_env(1)%force_env, &
1364 4622 : calc_force=.FALSE.)
1365 : CALL force_env_get(force_env(1)%force_env, &
1366 4622 : potential_energy=trial_energy)
1367 :
1368 13866 : DO itemp = 1, nvirial_temps
1369 :
1370 9244 : exponent = -BETA(itemp)*trial_energy
1371 :
1372 9244 : IF (exponent .GT. exp_max_val) THEN
1373 : exponent = exp_max_val
1374 9244 : ELSEIF (exponent .LT. exp_min_val) THEN
1375 1750 : exponent = exp_min_val
1376 : END IF
1377 13866 : mayer(itemp, ibin) = mayer(itemp, ibin) + EXP(exponent) - 1.0_dp
1378 : END DO
1379 : END IF
1380 :
1381 5900 : ibin = ibin + 1
1382 : END DO
1383 : ! write out some info that keeps track of where we are
1384 22 : IF (iw > 0) THEN
1385 10 : IF (MOD(ivirial, iprint) == 0) &
1386 10 : WRITE (iw, '(A,I6,A,I6)') ' Done with config ', ivirial, ' out of ', nvirial
1387 : END IF
1388 : END DO
1389 :
1390 : ! now we integrate this average potential
1391 1778 : mayer(:, :) = mayer(:, :)/REAL(nvirial, dp)
1392 :
1393 6 : DO itemp = 1, nvirial_temps
1394 : integral = 0.0_dp
1395 : previous_value = 0.0_dp
1396 : distance = 0.0E0_dp
1397 : ibin = 1
1398 1180 : DO
1399 1184 : current_division = 0
1400 1756 : DO idivision = 1, nintegral_divisions
1401 1756 : IF (distance .LT. virial_cutoffs(idivision) - virial_stepsize(idivision)/2.0E0_dp) THEN
1402 : current_division = idivision
1403 : EXIT
1404 : END IF
1405 : END DO
1406 1184 : IF (current_division == 0) EXIT
1407 1180 : distance = distance + virial_stepsize(current_division)
1408 :
1409 : ! now we need to integrate, using the trapazoidal method
1410 : ! first, find the value of the square
1411 1180 : current_value = mayer(itemp, ibin)*distance**2
1412 1180 : square_value = previous_value*virial_stepsize(current_division)
1413 : ! now the triangle that sits on top of it, which is half the size of this square...
1414 : ! notice this is negative if the current value is less than the previous value
1415 1180 : triangle_value = 0.5E0_dp*((current_value - previous_value)*virial_stepsize(current_division))
1416 :
1417 1180 : integral = integral + square_value + triangle_value
1418 1180 : previous_value = current_value
1419 1180 : ibin = ibin + 1
1420 : END DO
1421 :
1422 : ! now that the integration is done, compute the second virial that results
1423 4 : ave_virial = -2.0E0_dp*pi*integral
1424 :
1425 : ! convert from CP2K units to something else
1426 4 : ave_virial = ave_virial*n_avogadro*angstrom**3/1.0E8_dp**3
1427 :
1428 6 : IF (iw > 0) THEN
1429 2 : WRITE (iw, *)
1430 2 : WRITE (iw, *) '*********************************************************************'
1431 2 : WRITE (iw, '(A,F12.6,A)') ' *** Temperature = ', virial_temps(itemp), &
1432 4 : ' ***'
1433 2 : WRITE (iw, *) '*** ***'
1434 2 : WRITE (iw, '(A,E12.6,A)') ' *** B2(T) = ', ave_virial, &
1435 4 : ' cm**3/mol ***'
1436 2 : WRITE (iw, *) '*********************************************************************'
1437 2 : WRITE (iw, *)
1438 : END IF
1439 : END DO
1440 :
1441 : ! deallocate some stuff
1442 2 : DEALLOCATE (mc_par)
1443 2 : DEALLOCATE (subsys)
1444 2 : DEALLOCATE (force_env)
1445 2 : DEALLOCATE (particles)
1446 2 : DEALLOCATE (cell)
1447 2 : DEALLOCATE (virial_cutoffs)
1448 2 : DEALLOCATE (virial_stepsize)
1449 2 : DEALLOCATE (r_old)
1450 2 : DEALLOCATE (mayer)
1451 2 : DEALLOCATE (BETA)
1452 :
1453 6 : END SUBROUTINE mc_compute_virial
1454 :
1455 : END MODULE mc_ensembles
1456 :
|
