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