1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief holds all the structure types needed for Monte Carlo, except 8!> the mc_environment_type 9!> \par History 10!> none 11!> \author MJM 12! ************************************************************************************************** 13MODULE mc_types 14 USE atomic_kind_types, ONLY: atomic_kind_type,& 15 get_atomic_kind 16 USE cell_types, ONLY: cell_type,& 17 get_cell 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_type 24 USE cp_units, ONLY: cp_unit_to_cp2k 25 USE fist_environment_types, ONLY: fist_env_get,& 26 fist_environment_type 27 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get,& 28 fist_nonbond_env_type 29 USE force_env_types, ONLY: force_env_get,& 30 force_env_p_type,& 31 force_env_type,& 32 use_fist_force,& 33 use_qs_force 34 USE global_types, ONLY: global_environment_type 35 USE input_constants, ONLY: do_fist,& 36 do_qs 37 USE input_section_types, ONLY: section_vals_get_subs_vals,& 38 section_vals_type,& 39 section_vals_val_get,& 40 section_vals_val_set 41 USE kinds, ONLY: default_path_length,& 42 default_string_length,& 43 dp 44 USE mathconstants, ONLY: pi 45 USE molecule_kind_list_types, ONLY: molecule_kind_list_p_type 46 USE molecule_kind_types, ONLY: atom_type,& 47 get_molecule_kind,& 48 molecule_kind_type 49 USE pair_potential_types, ONLY: pair_potential_pp_type 50 USE particle_list_types, ONLY: particle_list_p_type 51 USE physcon, ONLY: boltzmann,& 52 joule 53 USE string_utilities, ONLY: uppercase,& 54 xstring 55#include "../../base/base_uses.f90" 56 57 IMPLICIT NONE 58 59 PRIVATE 60 61! *** Global parameters *** 62 63 PUBLIC :: mc_simpar_type, & 64 mc_simulation_parameters_p_type, & 65 mc_averages_type, mc_averages_p_type, & 66 mc_moves_type, mc_moves_p_type, accattempt, & 67 get_mc_par, set_mc_par, read_mc_section, & 68 find_mc_rcut, mc_molecule_info_type, & 69 mc_determine_molecule_info, mc_molecule_info_destroy, & 70 mc_sim_par_create, mc_sim_par_destroy, & 71 get_mc_molecule_info, & 72 mc_input_file_type, mc_input_file_create, & 73 get_mc_input_file, mc_input_file_destroy, & 74 mc_input_parameters_check, & 75 mc_ekin_type 76 77! ************************************************************************************************** 78 TYPE mc_simpar_type 79! contains all the parameters needed for running Monte Carlo simulations 80 PRIVATE 81 INTEGER, DIMENSION(:), POINTER :: avbmc_atom 82 INTEGER :: nstep 83 INTEGER :: iupvolume 84 INTEGER :: iuptrans 85 INTEGER :: iupcltrans 86 INTEGER :: nvirial 87 INTEGER :: nbox 88 INTEGER :: nmoves 89 INTEGER :: nswapmoves 90 INTEGER :: rm 91 INTEGER :: cl 92 INTEGER :: diff 93 INTEGER :: nstart 94 INTEGER :: source 95 INTEGER :: group 96 INTEGER :: iprint 97 LOGICAL :: ldiscrete 98 LOGICAL :: lbias 99 LOGICAL :: ionode 100 LOGICAL :: lrestart 101 LOGICAL :: lstop 102 LOGICAL :: lhmc 103 CHARACTER(LEN=20) :: ensemble 104 CHARACTER(LEN=default_path_length) :: restart_file_name 105 CHARACTER(LEN=default_path_length) :: molecules_file 106 CHARACTER(LEN=default_path_length) :: moves_file 107 CHARACTER(LEN=default_path_length) :: coords_file 108 CHARACTER(LEN=default_path_length) :: energy_file 109 CHARACTER(LEN=default_path_length) :: displacement_file 110 CHARACTER(LEN=default_path_length) :: cell_file 111 CHARACTER(LEN=default_path_length) :: dat_file 112 CHARACTER(LEN=default_path_length) :: data_file 113 CHARACTER(LEN=default_path_length) :: box2_file 114 CHARACTER(LEN=200) :: fft_lib 115 CHARACTER(LEN=50) :: PROGRAM 116 REAL(dp), DIMENSION(:), POINTER :: pmtraion_mol, pmtrans_mol, pmrot_mol, pmswap_mol, & 117 pbias, pmavbmc_mol 118 REAL(dp) :: discrete_step 119 REAL(dp) :: rmvolume, rmcltrans 120 REAL(dp), DIMENSION(:), POINTER :: rmbond, rmangle, rmdihedral, rmrot, rmtrans 121 REAL(dp), DIMENSION(:), POINTER :: eta 122 REAL(dp) :: temperature 123 REAL(dp) :: pressure 124 REAL(dp) :: rclus 125 REAL(dp) :: pmavbmc 126 REAL(dp) :: pmswap 127 REAL(dp) :: pmvolume, pmvol_box, pmclus_box 128 REAL(dp) :: pmhmc, pmhmc_box 129 REAL(dp) :: pmtraion 130 REAL(dp) :: pmtrans 131 REAL(dp) :: pmcltrans 132 REAL(dp) :: BETA 133 REAL(dp) :: rcut 134 REAL(dp), DIMENSION(:), POINTER :: avbmc_rmin, avbmc_rmax 135 REAL(dp), DIMENSION(:), POINTER :: virial_temps 136 TYPE(mc_input_file_type), POINTER :: & 137 mc_input_file, mc_bias_file 138 TYPE(section_vals_type), POINTER :: input_file 139 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 140 REAL(dp) :: exp_min_val, exp_max_val, min_val, max_val 141 INTEGER :: rand2skip 142 END TYPE mc_simpar_type 143 144! ************************************************************************************************** 145 TYPE mc_ekin_type 146! contains the kinetic energy of an MD sequence from hybrid MC 147 REAL(dp) :: initial_ekin, final_ekin 148 END TYPE mc_ekin_type 149! ************************************************************************************************** 150 TYPE mc_input_file_type 151! contains all the text of the input file 152 PRIVATE 153 INTEGER :: run_type_row, run_type_column, coord_row_start, coord_row_end, & 154 cell_row, cell_column, global_row_end, in_use, nunits_empty, & 155 motion_row_end, motion_row_start 156 INTEGER, DIMENSION(:), POINTER :: mol_set_nmol_row, mol_set_nmol_column 157 CHARACTER(LEN=default_string_length), DIMENSION(:), POINTER :: text 158 CHARACTER(LEN=default_string_length), DIMENSION(:), POINTER :: atom_names_empty 159 REAL(dp), DIMENSION(:, :), POINTER :: coordinates_empty 160 END TYPE mc_input_file_type 161 162! ************************************************************************************************** 163 TYPE mc_molecule_info_type 164! contains information on molecules...I had to do this because I use multiple force 165! environments, and they know nothing about each other 166 PRIVATE 167 INTEGER :: nboxes 168 CHARACTER(LEN=default_string_length), DIMENSION(:), POINTER :: names 169 CHARACTER(LEN=default_string_length), & 170 DIMENSION(:, :), POINTER :: atom_names 171 REAL(dp), DIMENSION(:, :), POINTER :: conf_prob, mass 172 INTEGER, DIMENSION(:, :), POINTER :: nchains 173 INTEGER :: nmol_types, nchain_total 174 INTEGER, DIMENSION(:), POINTER :: nunits, mol_type, nunits_tot, in_box 175 END TYPE mc_molecule_info_type 176 177! ************************************************************************************************** 178 TYPE mc_simulation_parameters_p_type 179! a pointer for multiple boxes 180 TYPE(mc_simpar_type), POINTER :: mc_par 181 END TYPE mc_simulation_parameters_p_type 182 183! ************************************************************************************************** 184 TYPE mc_averages_type 185! contains some data that is averaged throughout the course of a run 186 REAL(KIND=dp) :: ave_energy 187 REAL(KIND=dp) :: ave_energy_squared 188 REAL(KIND=dp) :: ave_volume 189 REAL(KIND=dp) :: molecules 190 END TYPE mc_averages_type 191 192! ************************************************************************************************** 193 TYPE mc_averages_p_type 194 TYPE(mc_averages_type), POINTER :: averages 195 END TYPE mc_averages_p_type 196 197! ************************************************************************************************** 198 TYPE mc_moves_type 199! contains data for how many moves of a give type we've accepted/rejected 200 TYPE(accattempt), POINTER :: bias_bond 201 TYPE(accattempt), POINTER :: bias_angle 202 TYPE(accattempt), POINTER :: bias_dihedral 203 TYPE(accattempt), POINTER :: bias_trans 204 TYPE(accattempt), POINTER :: bias_cltrans 205 TYPE(accattempt), POINTER :: bias_rot 206 TYPE(accattempt), POINTER :: bond 207 TYPE(accattempt), POINTER :: angle 208 TYPE(accattempt), POINTER :: dihedral 209 TYPE(accattempt), POINTER :: trans 210 TYPE(accattempt), POINTER :: cltrans 211 TYPE(accattempt), POINTER :: rot 212 TYPE(accattempt), POINTER :: swap 213 TYPE(accattempt), POINTER :: volume 214 TYPE(accattempt), POINTER :: hmc 215 TYPE(accattempt), POINTER :: avbmc_inin 216 TYPE(accattempt), POINTER :: avbmc_outin 217 TYPE(accattempt), POINTER :: avbmc_inout 218 TYPE(accattempt), POINTER :: avbmc_outout 219 TYPE(accattempt), POINTER :: Quickstep 220 REAL(KIND=dp) :: trans_dis, qtrans_dis 221 REAL(KIND=dp) :: cltrans_dis, qcltrans_dis 222 INTEGER :: empty, grown, empty_conf, empty_avbmc 223 END TYPE mc_moves_type 224 225! ************************************************************************************************** 226 TYPE accattempt 227 INTEGER :: successes 228 INTEGER :: qsuccesses 229 INTEGER :: attempts 230 END TYPE accattempt 231 232! ************************************************************************************************** 233 TYPE mc_moves_p_type 234 TYPE(mc_moves_type), POINTER :: moves 235 END TYPE mc_moves_p_type 236 237! *** Global parameters *** 238 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_types' 239 240CONTAINS 241 242! ************************************************************************************************** 243!> \brief accesses the private elements of the mc_input_file_type 244!> \param mc_input_file the input file you want data for 245!> 246!> Suitable for parallel. 247!> \param run_type_row ... 248!> \param run_type_column ... 249!> \param coord_row_start ... 250!> \param coord_row_end ... 251!> \param cell_row ... 252!> \param cell_column ... 253!> \param global_row_end ... 254!> \param mol_set_nmol_row ... 255!> \param mol_set_nmol_column ... 256!> \param in_use ... 257!> \param text ... 258!> \param atom_names_empty ... 259!> \param nunits_empty ... 260!> \param coordinates_empty ... 261!> \param motion_row_end ... 262!> \param motion_row_start ... 263!> \author MJM 264! ************************************************************************************************** 265 SUBROUTINE get_mc_input_file(mc_input_file, run_type_row, run_type_column, & 266 coord_row_start, coord_row_end, cell_row, cell_column, global_row_end, & 267 mol_set_nmol_row, mol_set_nmol_column, in_use, text, atom_names_empty, & 268 nunits_empty, coordinates_empty, motion_row_end, motion_row_start) 269 270 TYPE(mc_input_file_type), POINTER :: mc_input_file 271 INTEGER, INTENT(OUT), OPTIONAL :: run_type_row, run_type_column, & 272 coord_row_start, coord_row_end, & 273 cell_row, cell_column, global_row_end 274 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: mol_set_nmol_row, mol_set_nmol_column 275 INTEGER, INTENT(OUT), OPTIONAL :: in_use 276 CHARACTER(LEN=default_string_length), & 277 DIMENSION(:), OPTIONAL, POINTER :: text, atom_names_empty 278 INTEGER, INTENT(OUT), OPTIONAL :: nunits_empty 279 REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: coordinates_empty 280 INTEGER, INTENT(OUT), OPTIONAL :: motion_row_end, motion_row_start 281 282 IF (PRESENT(in_use)) in_use = mc_input_file%in_use 283 IF (PRESENT(run_type_row)) run_type_row = mc_input_file%run_type_row 284 IF (PRESENT(run_type_column)) run_type_column = mc_input_file%run_type_column 285 IF (PRESENT(coord_row_start)) coord_row_start = mc_input_file%coord_row_start 286 IF (PRESENT(coord_row_end)) coord_row_end = mc_input_file%coord_row_end 287 IF (PRESENT(cell_row)) cell_row = mc_input_file%cell_row 288 IF (PRESENT(cell_column)) cell_column = mc_input_file%cell_column 289 IF (PRESENT(global_row_end)) global_row_end = mc_input_file%global_row_end 290 IF (PRESENT(nunits_empty)) nunits_empty = mc_input_file%nunits_empty 291 IF (PRESENT(motion_row_start)) motion_row_start = mc_input_file%motion_row_start 292 IF (PRESENT(motion_row_end)) motion_row_end = mc_input_file%motion_row_end 293 294 IF (PRESENT(mol_set_nmol_row)) mol_set_nmol_row => mc_input_file%mol_set_nmol_row 295 IF (PRESENT(mol_set_nmol_column)) mol_set_nmol_column => mc_input_file%mol_set_nmol_column 296 IF (PRESENT(text)) text => mc_input_file%text 297 IF (PRESENT(atom_names_empty)) atom_names_empty => mc_input_file%atom_names_empty 298 IF (PRESENT(coordinates_empty)) coordinates_empty => mc_input_file%coordinates_empty 299 300 END SUBROUTINE GET_MC_INPUT_FILE 301! ************************************************************************************************** 302!> \brief ... 303!> \param mc_par ... 304!> \param nstep ... 305!> \param nvirial ... 306!> \param iuptrans ... 307!> \param iupcltrans ... 308!> \param iupvolume ... 309!> \param nmoves ... 310!> \param nswapmoves ... 311!> \param rm ... 312!> \param cl ... 313!> \param diff ... 314!> \param nstart ... 315!> \param source ... 316!> \param group ... 317!> \param lbias ... 318!> \param ionode ... 319!> \param lrestart ... 320!> \param lstop ... 321!> \param rmvolume ... 322!> \param rmcltrans ... 323!> \param rmbond ... 324!> \param rmangle ... 325!> \param rmrot ... 326!> \param rmtrans ... 327!> \param temperature ... 328!> \param pressure ... 329!> \param rclus ... 330!> \param BETA ... 331!> \param pmswap ... 332!> \param pmvolume ... 333!> \param pmtraion ... 334!> \param pmtrans ... 335!> \param pmcltrans ... 336!> \param ensemble ... 337!> \param PROGRAM ... 338!> \param restart_file_name ... 339!> \param molecules_file ... 340!> \param moves_file ... 341!> \param coords_file ... 342!> \param energy_file ... 343!> \param displacement_file ... 344!> \param cell_file ... 345!> \param dat_file ... 346!> \param data_file ... 347!> \param box2_file ... 348!> \param fft_lib ... 349!> \param iprint ... 350!> \param rcut ... 351!> \param ldiscrete ... 352!> \param discrete_step ... 353!> \param pmavbmc ... 354!> \param pbias ... 355!> \param avbmc_atom ... 356!> \param avbmc_rmin ... 357!> \param avbmc_rmax ... 358!> \param rmdihedral ... 359!> \param input_file ... 360!> \param mc_molecule_info ... 361!> \param pmswap_mol ... 362!> \param pmavbmc_mol ... 363!> \param pmtrans_mol ... 364!> \param pmrot_mol ... 365!> \param pmtraion_mol ... 366!> \param mc_input_file ... 367!> \param mc_bias_file ... 368!> \param pmvol_box ... 369!> \param pmclus_box ... 370!> \param virial_temps ... 371!> \param exp_min_val ... 372!> \param exp_max_val ... 373!> \param min_val ... 374!> \param max_val ... 375!> \param eta ... 376!> \param pmhmc ... 377!> \param pmhmc_box ... 378!> \param lhmc ... 379!> \param rand2skip ... 380! ************************************************************************************************** 381 SUBROUTINE get_mc_par(mc_par, nstep, nvirial, iuptrans, iupcltrans, iupvolume, & 382 nmoves, nswapmoves, rm, cl, diff, nstart, & 383 source, group, lbias, ionode, lrestart, lstop, rmvolume, rmcltrans, rmbond, rmangle, & 384 rmrot, rmtrans, temperature, pressure, rclus, BETA, pmswap, pmvolume, pmtraion, pmtrans, & 385 pmcltrans, ensemble, PROGRAM, restart_file_name, molecules_file, moves_file, coords_file, & 386 energy_file, displacement_file, cell_file, dat_file, data_file, box2_file, & 387 fft_lib, iprint, rcut, ldiscrete, discrete_step, & 388 pmavbmc, pbias, avbmc_atom, avbmc_rmin, avbmc_rmax, rmdihedral, & 389 input_file, mc_molecule_info, pmswap_mol, pmavbmc_mol, pmtrans_mol, pmrot_mol, & 390 pmtraion_mol, mc_input_file, mc_bias_file, pmvol_box, pmclus_box, virial_temps, & 391 exp_min_val, exp_max_val, min_val, max_val, eta, pmhmc, pmhmc_box, lhmc, rand2skip) 392 393! ************************************************************************************************** 394!> \brief accesses the private elements of the mc_parameters_type 395!> \param mc_par the structure mc parameters you want 396!> 397!> Suitable for parallel. 398!> \author MJM 399 TYPE(mc_simpar_type), POINTER :: mc_par 400 INTEGER, INTENT(OUT), OPTIONAL :: nstep, nvirial, iuptrans, iupcltrans, & 401 iupvolume, nmoves, nswapmoves, rm, cl, & 402 diff, nstart, source, group 403 LOGICAL, INTENT(OUT), OPTIONAL :: lbias, ionode, lrestart, lstop 404 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: rmvolume, rmcltrans 405 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: rmbond, rmangle, rmrot, rmtrans 406 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: temperature, pressure, rclus, BETA, & 407 pmswap, pmvolume, pmtraion, pmtrans, & 408 pmcltrans 409 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: ensemble, PROGRAM, restart_file_name, & 410 molecules_file, moves_file, coords_file, energy_file, displacement_file, cell_file, & 411 dat_file, data_file, box2_file, fft_lib 412 INTEGER, INTENT(OUT), OPTIONAL :: iprint 413 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: rcut 414 LOGICAL, INTENT(OUT), OPTIONAL :: ldiscrete 415 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: discrete_step, pmavbmc 416 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: pbias 417 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: avbmc_atom 418 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: avbmc_rmin, avbmc_rmax, rmdihedral 419 TYPE(section_vals_type), OPTIONAL, POINTER :: input_file 420 TYPE(mc_molecule_info_type), OPTIONAL, POINTER :: mc_molecule_info 421 REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: pmswap_mol 422 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: pmavbmc_mol, pmtrans_mol, pmrot_mol, & 423 pmtraion_mol 424 TYPE(mc_input_file_type), OPTIONAL, POINTER :: mc_input_file, mc_bias_file 425 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: pmvol_box, pmclus_box 426 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: virial_temps 427 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: exp_min_val, exp_max_val, min_val, & 428 max_val 429 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eta 430 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: pmhmc, pmhmc_box 431 LOGICAL, INTENT(OUT), OPTIONAL :: lhmc 432 INTEGER, INTENT(OUT), OPTIONAL :: rand2skip 433 434 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mc_par', routineP = moduleN//':'//routineN 435 436 IF (PRESENT(nstep)) nstep = mc_par%nstep 437 IF (PRESENT(nvirial)) nvirial = mc_par%nvirial 438 IF (PRESENT(iuptrans)) iuptrans = mc_par%iuptrans 439 IF (PRESENT(iupcltrans)) iupcltrans = mc_par%iupcltrans 440 IF (PRESENT(iupvolume)) iupvolume = mc_par%iupvolume 441 IF (PRESENT(nmoves)) nmoves = mc_par%nmoves 442 IF (PRESENT(nswapmoves)) nswapmoves = mc_par%nswapmoves 443 IF (PRESENT(rm)) rm = mc_par%rm 444 IF (PRESENT(cl)) cl = mc_par%cl 445 IF (PRESENT(diff)) diff = mc_par%diff 446 IF (PRESENT(nstart)) nstart = mc_par%nstart 447 IF (PRESENT(source)) source = mc_par%source 448 IF (PRESENT(group)) group = mc_par%group 449 IF (PRESENT(iprint)) iprint = mc_par%iprint 450 451 IF (PRESENT(lbias)) lbias = mc_par%lbias 452 IF (PRESENT(ionode)) ionode = mc_par%ionode 453 IF (PRESENT(lrestart)) lrestart = mc_par%lrestart 454 IF (PRESENT(lstop)) lstop = mc_par%lstop 455 IF (PRESENT(ldiscrete)) ldiscrete = mc_par%ldiscrete 456 457 IF (PRESENT(virial_temps)) virial_temps => mc_par%virial_temps 458 IF (PRESENT(avbmc_atom)) avbmc_atom => mc_par%avbmc_atom 459 IF (PRESENT(avbmc_rmin)) avbmc_rmin => mc_par%avbmc_rmin 460 IF (PRESENT(avbmc_rmax)) avbmc_rmax => mc_par%avbmc_rmax 461 IF (PRESENT(rcut)) rcut = mc_par%rcut 462 IF (PRESENT(discrete_step)) discrete_step = mc_par%discrete_step 463 IF (PRESENT(rmvolume)) rmvolume = mc_par%rmvolume 464 IF (PRESENT(rmcltrans)) rmcltrans = mc_par%rmcltrans 465 IF (PRESENT(temperature)) temperature = mc_par%temperature 466 IF (PRESENT(pressure)) pressure = mc_par%pressure 467 IF (PRESENT(rclus)) rclus = mc_par%rclus 468 IF (PRESENT(BETA)) BETA = mc_par%BETA 469 IF (PRESENT(exp_max_val)) exp_max_val = mc_par%exp_max_val 470 IF (PRESENT(exp_min_val)) exp_min_val = mc_par%exp_min_val 471 IF (PRESENT(max_val)) max_val = mc_par%max_val 472 IF (PRESENT(min_val)) min_val = mc_par%min_val 473 IF (PRESENT(pmswap)) pmswap = mc_par%pmswap 474 IF (PRESENT(pmvolume)) pmvolume = mc_par%pmvolume 475 IF (PRESENT(lhmc)) lhmc = mc_par%lhmc 476 IF (PRESENT(pmhmc)) pmhmc = mc_par%pmhmc 477 IF (PRESENT(pmhmc_box)) pmhmc_box = mc_par%pmhmc_box 478 IF (PRESENT(pmvol_box)) pmvol_box = mc_par%pmvol_box 479 IF (PRESENT(pmclus_box)) pmclus_box = mc_par%pmclus_box 480 IF (PRESENT(pmtraion)) pmtraion = mc_par%pmtraion 481 IF (PRESENT(pmtrans)) pmtrans = mc_par%pmtrans 482 IF (PRESENT(pmcltrans)) pmcltrans = mc_par%pmcltrans 483 IF (PRESENT(pmavbmc)) pmavbmc = mc_par%pmavbmc 484 IF (PRESENT(pmrot_mol)) pmrot_mol => mc_par%pmrot_mol 485 IF (PRESENT(pmtrans_mol)) pmtrans_mol => mc_par%pmtrans_mol 486 IF (PRESENT(pmtraion_mol)) pmtraion_mol => mc_par%pmtraion_mol 487 IF (PRESENT(pmavbmc_mol)) pmavbmc_mol => mc_par%pmavbmc_mol 488 IF (PRESENT(pbias)) pbias => mc_par%pbias 489 490 IF (PRESENT(rmbond)) rmbond => mc_par%rmbond 491 IF (PRESENT(rmangle)) rmangle => mc_par%rmangle 492 IF (PRESENT(rmdihedral)) rmdihedral => mc_par%rmdihedral 493 IF (PRESENT(rmrot)) rmrot => mc_par%rmrot 494 IF (PRESENT(rmtrans)) rmtrans => mc_par%rmtrans 495 496 IF (PRESENT(eta)) eta => mc_par%eta 497 498 IF (PRESENT(ensemble)) ensemble = mc_par%ensemble 499 IF (PRESENT(PROGRAM)) PROGRAM = mc_par%program 500 IF (PRESENT(restart_file_name)) restart_file_name = & 501 mc_par%restart_file_name 502 IF (PRESENT(moves_file)) moves_file = mc_par%moves_file 503 IF (PRESENT(coords_file)) coords_file = mc_par%coords_file 504 IF (PRESENT(molecules_file)) molecules_file = mc_par%molecules_file 505 IF (PRESENT(energy_file)) energy_file = mc_par%energy_file 506 IF (PRESENT(displacement_file)) displacement_file = & 507 mc_par%displacement_file 508 IF (PRESENT(cell_file)) cell_file = mc_par%cell_file 509 IF (PRESENT(dat_file)) dat_file = mc_par%dat_file 510 IF (PRESENT(data_file)) data_file = mc_par%data_file 511 IF (PRESENT(box2_file)) box2_file = mc_par%box2_file 512 IF (PRESENT(fft_lib)) fft_lib = mc_par%fft_lib 513 514 IF (PRESENT(input_file)) input_file => mc_par%input_file 515 IF (PRESENT(mc_molecule_info)) mc_molecule_info => mc_par%mc_molecule_info 516 IF (PRESENT(mc_input_file)) mc_input_file => mc_par%mc_input_file 517 IF (PRESENT(mc_bias_file)) mc_bias_file => mc_par%mc_bias_file 518 519 IF (PRESENT(pmswap_mol)) pmswap_mol => mc_par%pmswap_mol 520 IF (PRESENT(rand2skip)) rand2skip = mc_par%rand2skip 521 END SUBROUTINE get_mc_par 522 523! ************************************************************************************************** 524!> \brief ... 525!> \param mc_molecule_info ... 526!> \param nmol_types ... 527!> \param nchain_total ... 528!> \param nboxes ... 529!> \param names ... 530!> \param conf_prob ... 531!> \param nchains ... 532!> \param nunits ... 533!> \param mol_type ... 534!> \param nunits_tot ... 535!> \param in_box ... 536!> \param atom_names ... 537!> \param mass ... 538! ************************************************************************************************** 539 SUBROUTINE get_mc_molecule_info(mc_molecule_info, nmol_types, nchain_total, nboxes, & 540 names, conf_prob, nchains, nunits, mol_type, nunits_tot, in_box, atom_names, & 541 mass) 542 543! ************************************************************************************************** 544!> \brief accesses the private elements of the mc_molecule_info_type 545!> \param mc_molecule_info the structure you want the parameters for 546!> 547!> Suitable for parallel. 548!> \author MJM 549 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 550 INTEGER, INTENT(OUT), OPTIONAL :: nmol_types, nchain_total, nboxes 551 CHARACTER(LEN=default_string_length), & 552 DIMENSION(:), OPTIONAL, POINTER :: names 553 REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: conf_prob 554 INTEGER, DIMENSION(:, :), OPTIONAL, POINTER :: nchains 555 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nunits, mol_type, nunits_tot, in_box 556 CHARACTER(LEN=default_string_length), & 557 DIMENSION(:, :), OPTIONAL, POINTER :: atom_names 558 REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: mass 559 560 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mc_molecule_info', & 561 routineP = moduleN//':'//routineN 562 563 IF (PRESENT(nboxes)) nboxes = mc_molecule_info%nboxes 564 IF (PRESENT(nchain_total)) nchain_total = mc_molecule_info%nchain_total 565 IF (PRESENT(nmol_types)) nmol_types = mc_molecule_info%nmol_types 566 567 IF (PRESENT(names)) names => mc_molecule_info%names 568 IF (PRESENT(atom_names)) atom_names => mc_molecule_info%atom_names 569 570 IF (PRESENT(conf_prob)) conf_prob => mc_molecule_info%conf_prob 571 IF (PRESENT(mass)) mass => mc_molecule_info%mass 572 573 IF (PRESENT(nchains)) nchains => mc_molecule_info%nchains 574 IF (PRESENT(nunits)) nunits => mc_molecule_info%nunits 575 IF (PRESENT(mol_type)) mol_type => mc_molecule_info%mol_type 576 IF (PRESENT(nunits_tot)) nunits_tot => mc_molecule_info%nunits_tot 577 IF (PRESENT(in_box)) in_box => mc_molecule_info%in_box 578 579 END SUBROUTINE get_mc_molecule_info 580 581! ************************************************************************************************** 582!> \brief ... 583!> \param mc_par ... 584!> \param rm ... 585!> \param cl ... 586!> \param diff ... 587!> \param nstart ... 588!> \param rmvolume ... 589!> \param rmcltrans ... 590!> \param rmbond ... 591!> \param rmangle ... 592!> \param rmdihedral ... 593!> \param rmrot ... 594!> \param rmtrans ... 595!> \param PROGRAM ... 596!> \param nmoves ... 597!> \param nswapmoves ... 598!> \param lstop ... 599!> \param temperature ... 600!> \param pressure ... 601!> \param rclus ... 602!> \param iuptrans ... 603!> \param iupcltrans ... 604!> \param iupvolume ... 605!> \param pmswap ... 606!> \param pmvolume ... 607!> \param pmtraion ... 608!> \param pmtrans ... 609!> \param pmcltrans ... 610!> \param BETA ... 611!> \param rcut ... 612!> \param iprint ... 613!> \param lbias ... 614!> \param nstep ... 615!> \param lrestart ... 616!> \param ldiscrete ... 617!> \param discrete_step ... 618!> \param pmavbmc ... 619!> \param mc_molecule_info ... 620!> \param pmavbmc_mol ... 621!> \param pmtrans_mol ... 622!> \param pmrot_mol ... 623!> \param pmtraion_mol ... 624!> \param pmswap_mol ... 625!> \param avbmc_rmin ... 626!> \param avbmc_rmax ... 627!> \param avbmc_atom ... 628!> \param pbias ... 629!> \param ensemble ... 630!> \param pmvol_box ... 631!> \param pmclus_box ... 632!> \param eta ... 633!> \param mc_input_file ... 634!> \param mc_bias_file ... 635!> \param exp_max_val ... 636!> \param exp_min_val ... 637!> \param min_val ... 638!> \param max_val ... 639!> \param pmhmc ... 640!> \param pmhmc_box ... 641!> \param lhmc ... 642!> \param ionode ... 643!> \param source ... 644!> \param group ... 645!> \param rand2skip ... 646! ************************************************************************************************** 647 SUBROUTINE set_mc_par(mc_par, rm, cl, & 648 diff, nstart, rmvolume, rmcltrans, rmbond, rmangle, rmdihedral, rmrot, rmtrans, PROGRAM, & 649 nmoves, nswapmoves, lstop, temperature, pressure, rclus, iuptrans, iupcltrans, iupvolume, & 650 pmswap, pmvolume, pmtraion, pmtrans, pmcltrans, BETA, rcut, iprint, lbias, nstep, & 651 lrestart, ldiscrete, discrete_step, pmavbmc, mc_molecule_info, & 652 pmavbmc_mol, pmtrans_mol, pmrot_mol, pmtraion_mol, pmswap_mol, & 653 avbmc_rmin, avbmc_rmax, avbmc_atom, pbias, ensemble, pmvol_box, pmclus_box, eta, & 654 mc_input_file, mc_bias_file, exp_max_val, exp_min_val, min_val, max_val, & 655 pmhmc, pmhmc_box, lhmc, ionode, source, group, rand2skip) 656 657! ************************************************************************************************** 658!> \brief changes the private elements of the mc_parameters_type 659!> \param mc_par the structure mc parameters you want 660!> 661!> Suitable for parallel. 662!> \author MJM 663 TYPE(mc_simpar_type), POINTER :: mc_par 664 INTEGER, INTENT(IN), OPTIONAL :: rm, cl, diff, nstart 665 REAL(KIND=dp), INTENT(IN), OPTIONAL :: rmvolume, rmcltrans 666 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: rmbond, rmangle, rmdihedral, rmrot, & 667 rmtrans 668 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: PROGRAM 669 INTEGER, INTENT(IN), OPTIONAL :: nmoves, nswapmoves 670 LOGICAL, INTENT(IN), OPTIONAL :: lstop 671 REAL(KIND=dp), INTENT(IN), OPTIONAL :: temperature, pressure, rclus 672 INTEGER, INTENT(IN), OPTIONAL :: iuptrans, iupcltrans, iupvolume 673 REAL(KIND=dp), INTENT(IN), OPTIONAL :: pmswap, pmvolume, pmtraion, pmtrans, & 674 pmcltrans, BETA, rcut 675 INTEGER, INTENT(IN), OPTIONAL :: iprint 676 LOGICAL, INTENT(IN), OPTIONAL :: lbias 677 INTEGER, INTENT(IN), OPTIONAL :: nstep 678 LOGICAL, INTENT(IN), OPTIONAL :: lrestart, ldiscrete 679 REAL(KIND=dp), INTENT(IN), OPTIONAL :: discrete_step, pmavbmc 680 TYPE(mc_molecule_info_type), OPTIONAL, POINTER :: mc_molecule_info 681 REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: pmavbmc_mol, pmtrans_mol, pmrot_mol, & 682 pmtraion_mol, pmswap_mol, avbmc_rmin, & 683 avbmc_rmax 684 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: avbmc_atom 685 REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: pbias 686 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ensemble 687 REAL(KIND=dp), INTENT(IN), OPTIONAL :: pmvol_box, pmclus_box 688 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eta 689 TYPE(mc_input_file_type), OPTIONAL, POINTER :: mc_input_file, mc_bias_file 690 REAL(KIND=dp), INTENT(IN), OPTIONAL :: exp_max_val, exp_min_val, min_val, & 691 max_val, pmhmc, pmhmc_box 692 LOGICAL, INTENT(IN), OPTIONAL :: lhmc, ionode 693 INTEGER, INTENT(IN), OPTIONAL :: source, group, rand2skip 694 695 CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mc_par', routineP = moduleN//':'//routineN 696 697 INTEGER :: iparm 698 699! These are the only values that change during the course of the simulation 700! or are computed outside of this module 701 702 IF (PRESENT(nstep)) mc_par%nstep = nstep 703 IF (PRESENT(rm)) mc_par%rm = rm 704 IF (PRESENT(cl)) mc_par%cl = cl 705 IF (PRESENT(diff)) mc_par%diff = diff 706 IF (PRESENT(nstart)) mc_par%nstart = nstart 707 IF (PRESENT(nmoves)) mc_par%nmoves = nmoves 708 IF (PRESENT(nswapmoves)) mc_par%nswapmoves = nswapmoves 709 IF (PRESENT(iprint)) mc_par%iprint = iprint 710 IF (PRESENT(iuptrans)) mc_par%iuptrans = iuptrans 711 IF (PRESENT(iupcltrans)) mc_par%iupcltrans = iupcltrans 712 IF (PRESENT(iupvolume)) mc_par%iupvolume = iupvolume 713 714 IF (PRESENT(ldiscrete)) mc_par%ldiscrete = ldiscrete 715 IF (PRESENT(lstop)) mc_par%lstop = lstop 716 IF (PRESENT(lbias)) mc_par%lbias = lbias 717 IF (PRESENT(lrestart)) mc_par%lrestart = lrestart 718 719 IF (PRESENT(exp_max_val)) mc_par%exp_max_val = exp_max_val 720 IF (PRESENT(exp_min_val)) mc_par%exp_min_val = exp_min_val 721 IF (PRESENT(max_val)) mc_par%max_val = max_val 722 IF (PRESENT(min_val)) mc_par%min_val = min_val 723 IF (PRESENT(BETA)) mc_par%BETA = BETA 724 IF (PRESENT(temperature)) mc_par%temperature = temperature 725 IF (PRESENT(rcut)) mc_par%rcut = rcut 726 IF (PRESENT(pressure)) mc_par%pressure = pressure 727 IF (PRESENT(rclus)) mc_par%rclus = rclus 728 IF (PRESENT(pmvolume)) mc_par%pmvolume = pmvolume 729 IF (PRESENT(lhmc)) mc_par%lhmc = lhmc 730 IF (PRESENT(pmhmc)) mc_par%pmhmc = pmhmc 731 IF (PRESENT(pmhmc_box)) mc_par%pmhmc_box = pmhmc_box 732 IF (PRESENT(pmvol_box)) mc_par%pmvol_box = pmvol_box 733 IF (PRESENT(pmclus_box)) mc_par%pmclus_box = pmclus_box 734 IF (PRESENT(pmswap)) mc_par%pmswap = pmswap 735 IF (PRESENT(pmtrans)) mc_par%pmtrans = pmtrans 736 IF (PRESENT(pmcltrans)) mc_par%pmcltrans = pmcltrans 737 IF (PRESENT(pmtraion)) mc_par%pmtraion = pmtraion 738 IF (PRESENT(pmavbmc)) mc_par%pmavbmc = pmavbmc 739 740 IF (PRESENT(discrete_step)) mc_par%discrete_step = discrete_step 741 IF (PRESENT(rmvolume)) mc_par%rmvolume = rmvolume 742 IF (PRESENT(rmcltrans)) mc_par%rmcltrans = rmcltrans 743 744 IF (PRESENT(mc_input_file)) mc_par%mc_input_file => mc_input_file 745 IF (PRESENT(mc_bias_file)) mc_par%mc_bias_file => mc_bias_file 746 747! cannot just change the pointers, because then we have memory leaks 748! and the program crashes if we try to release it (in certain cases) 749 IF (PRESENT(eta)) THEN 750 DO iparm = 1, SIZE(eta) 751 mc_par%eta(iparm) = eta(iparm) 752 ENDDO 753 ENDIF 754 IF (PRESENT(rmbond)) THEN 755 DO iparm = 1, SIZE(rmbond) 756 mc_par%rmbond(iparm) = rmbond(iparm) 757 ENDDO 758 ENDIF 759 IF (PRESENT(rmangle)) THEN 760 DO iparm = 1, SIZE(rmangle) 761 mc_par%rmangle(iparm) = rmangle(iparm) 762 ENDDO 763 ENDIF 764 IF (PRESENT(rmdihedral)) THEN 765 DO iparm = 1, SIZE(rmdihedral) 766 mc_par%rmdihedral(iparm) = rmdihedral(iparm) 767 ENDDO 768 ENDIF 769 IF (PRESENT(rmrot)) THEN 770 DO iparm = 1, SIZE(rmrot) 771 mc_par%rmrot(iparm) = rmrot(iparm) 772 ENDDO 773 ENDIF 774 IF (PRESENT(rmtrans)) THEN 775 DO iparm = 1, SIZE(rmtrans) 776 mc_par%rmtrans(iparm) = rmtrans(iparm) 777 ENDDO 778 ENDIF 779 780 IF (PRESENT(pmavbmc_mol)) THEN 781 DO iparm = 1, SIZE(pmavbmc_mol) 782 mc_par%pmavbmc_mol(iparm) = pmavbmc_mol(iparm) 783 ENDDO 784 ENDIF 785 IF (PRESENT(pmtrans_mol)) THEN 786 DO iparm = 1, SIZE(pmtrans_mol) 787 mc_par%pmtrans_mol(iparm) = pmtrans_mol(iparm) 788 ENDDO 789 ENDIF 790 IF (PRESENT(pmrot_mol)) THEN 791 DO iparm = 1, SIZE(pmrot_mol) 792 mc_par%pmrot_mol(iparm) = pmrot_mol(iparm) 793 ENDDO 794 ENDIF 795 IF (PRESENT(pmtraion_mol)) THEN 796 DO iparm = 1, SIZE(pmtraion_mol) 797 mc_par%pmtraion_mol(iparm) = pmtraion_mol(iparm) 798 ENDDO 799 ENDIF 800 IF (PRESENT(pmswap_mol)) THEN 801 DO iparm = 1, SIZE(pmswap_mol) 802 mc_par%pmswap_mol(iparm) = pmswap_mol(iparm) 803 ENDDO 804 ENDIF 805 806 IF (PRESENT(avbmc_atom)) THEN 807 DO iparm = 1, SIZE(avbmc_atom) 808 mc_par%avbmc_atom(iparm) = avbmc_atom(iparm) 809 ENDDO 810 ENDIF 811 IF (PRESENT(avbmc_rmin)) THEN 812 DO iparm = 1, SIZE(avbmc_rmin) 813 mc_par%avbmc_rmin(iparm) = avbmc_rmin(iparm) 814 ENDDO 815 ENDIF 816 IF (PRESENT(avbmc_rmax)) THEN 817 DO iparm = 1, SIZE(avbmc_rmax) 818 mc_par%avbmc_rmax(iparm) = avbmc_rmax(iparm) 819 ENDDO 820 ENDIF 821 IF (PRESENT(pbias)) THEN 822 DO iparm = 1, SIZE(pbias) 823 mc_par%pbias(iparm) = pbias(iparm) 824 ENDDO 825 ENDIF 826 827 IF (PRESENT(PROGRAM)) mc_par%program = PROGRAM 828 IF (PRESENT(ensemble)) mc_par%ensemble = ensemble 829 830 IF (PRESENT(mc_molecule_info)) mc_par%mc_molecule_info => mc_molecule_info 831 IF (PRESENT(source)) mc_par%source = source 832 IF (PRESENT(group)) mc_par%group = group 833 IF (PRESENT(ionode)) mc_par%ionode = ionode 834 835 IF (PRESENT(rand2skip)) mc_par%rand2skip = rand2skip 836 END SUBROUTINE set_mc_par 837 838! ************************************************************************************************** 839!> \brief creates (allocates) the mc_simulation_parameters type 840!> \param mc_par the structure that will be created (allocated) 841!> \param nmol_types the number of molecule types in the system 842!> \author MJM 843! ************************************************************************************************** 844 SUBROUTINE mc_sim_par_create(mc_par, nmol_types) 845 846 TYPE(mc_simpar_type), POINTER :: mc_par 847 INTEGER, INTENT(IN) :: nmol_types 848 849 CHARACTER(len=*), PARAMETER :: routineN = 'mc_sim_par_create', & 850 routineP = moduleN//':'//routineN 851 852 NULLIFY (mc_par) 853 854 ALLOCATE (mc_par) 855 ALLOCATE (mc_par%pmtraion_mol(1:nmol_types)) 856 ALLOCATE (mc_par%pmtrans_mol(1:nmol_types)) 857 ALLOCATE (mc_par%pmrot_mol(1:nmol_types)) 858 ALLOCATE (mc_par%pmswap_mol(1:nmol_types)) 859 860 ALLOCATE (mc_par%eta(1:nmol_types)) 861 862 ALLOCATE (mc_par%rmbond(1:nmol_types)) 863 ALLOCATE (mc_par%rmangle(1:nmol_types)) 864 ALLOCATE (mc_par%rmdihedral(1:nmol_types)) 865 ALLOCATE (mc_par%rmtrans(1:nmol_types)) 866 ALLOCATE (mc_par%rmrot(1:nmol_types)) 867 868 ALLOCATE (mc_par%avbmc_atom(1:nmol_types)) 869 ALLOCATE (mc_par%avbmc_rmin(1:nmol_types)) 870 ALLOCATE (mc_par%avbmc_rmax(1:nmol_types)) 871 ALLOCATE (mc_par%pmavbmc_mol(1:nmol_types)) 872 ALLOCATE (mc_par%pbias(1:nmol_types)) 873 874 ALLOCATE (mc_par%mc_input_file) 875 ALLOCATE (mc_par%mc_bias_file) 876 877 END SUBROUTINE mc_sim_par_create 878 879! ************************************************************************************************** 880!> \brief destroys (deallocates) the mc_simulation_parameters type 881!> \param mc_par the structure that will be destroyed 882!> \author MJM 883! ************************************************************************************************** 884 SUBROUTINE mc_sim_par_destroy(mc_par) 885 886 TYPE(mc_simpar_type), POINTER :: mc_par 887 888 CHARACTER(len=*), PARAMETER :: routineN = 'mc_sim_par_destroy', & 889 routineP = moduleN//':'//routineN 890 891 DEALLOCATE (mc_par%mc_input_file) 892 DEALLOCATE (mc_par%mc_bias_file) 893 894 DEALLOCATE (mc_par%pmtraion_mol) 895 DEALLOCATE (mc_par%pmtrans_mol) 896 DEALLOCATE (mc_par%pmrot_mol) 897 DEALLOCATE (mc_par%pmswap_mol) 898 899 DEALLOCATE (mc_par%eta) 900 901 DEALLOCATE (mc_par%rmbond) 902 DEALLOCATE (mc_par%rmangle) 903 DEALLOCATE (mc_par%rmdihedral) 904 DEALLOCATE (mc_par%rmtrans) 905 DEALLOCATE (mc_par%rmrot) 906 907 DEALLOCATE (mc_par%avbmc_atom) 908 DEALLOCATE (mc_par%avbmc_rmin) 909 DEALLOCATE (mc_par%avbmc_rmax) 910 DEALLOCATE (mc_par%pmavbmc_mol) 911 DEALLOCATE (mc_par%pbias) 912 IF (mc_par%ensemble == "VIRIAL") THEN 913 DEALLOCATE (mc_par%virial_temps) 914 END IF 915 916 DEALLOCATE (mc_par) 917 NULLIFY (mc_par) 918 919 END SUBROUTINE mc_sim_par_destroy 920 921! ************************************************************************************************** 922!> \brief creates (allocates) the mc_input_file_type 923!> \param mc_input_file the structure that will be created (allocated) 924!> \param input_file_name the name of the file to read 925!> \param mc_molecule_info ... 926!> \param empty_coords ... 927!> \param lhmc ... 928!> \author MJM 929! ************************************************************************************************** 930 SUBROUTINE mc_input_file_create(mc_input_file, input_file_name, & 931 mc_molecule_info, empty_coords, lhmc) 932 933 TYPE(mc_input_file_type), POINTER :: mc_input_file 934 CHARACTER(LEN=*), INTENT(IN) :: input_file_name 935 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 936 REAL(dp), DIMENSION(:, :) :: empty_coords 937 LOGICAL, INTENT(IN) :: lhmc 938 939 CHARACTER(len=*), PARAMETER :: routineN = 'mc_input_file_create', & 940 routineP = moduleN//':'//routineN 941 942 CHARACTER(default_string_length) :: cdum, line, method_name 943 CHARACTER(default_string_length), & 944 DIMENSION(:, :), POINTER :: atom_names 945 INTEGER :: abc_column, abc_row, cell_column, cell_row, idum, iline, io_stat, irep, iunit, & 946 iw, line_number, nlines, nmol_types, nstart, row_number, unit 947 INTEGER, DIMENSION(:), POINTER :: nunits 948 949! some stuff in case we need to write error messages to the screen 950 951 iw = cp_logger_get_default_io_unit() 952 953! allocate the array we'll need in case we have an empty box 954 CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types, & 955 nunits=nunits, atom_names=atom_names) 956 ALLOCATE (mc_input_file%coordinates_empty(1:3, 1:nunits(1))) 957 ALLOCATE (mc_input_file%atom_names_empty(1:nunits(1))) 958 DO iunit = 1, nunits(1) 959 mc_input_file%atom_names_empty(iunit) = atom_names(iunit, 1) 960 mc_input_file%coordinates_empty(1:3, iunit) = empty_coords(1:3, iunit) 961 ENDDO 962 mc_input_file%nunits_empty = nunits(1) 963 964! allocate some arrays 965 ALLOCATE (mc_input_file%mol_set_nmol_row(1:nmol_types)) 966 ALLOCATE (mc_input_file%mol_set_nmol_column(1:nmol_types)) 967 968! now read in and store the input file, first finding out how many lines it is 969 CALL open_file(file_name=input_file_name, & 970 unit_number=unit, file_action='READ', file_status='OLD') 971 972 nlines = 0 973 DO 974 READ (unit, '(A)', IOSTAT=io_stat) line 975 IF (io_stat .NE. 0) EXIT 976 nlines = nlines + 1 977 ENDDO 978 979 ALLOCATE (mc_input_file%text(1:nlines)) 980 981 REWIND (unit) 982 DO iline = 1, nlines 983 READ (unit, '(A)') mc_input_file%text(iline) 984 ENDDO 985 986! close the input file 987 CALL close_file(unit_number=unit) 988 989! now we need to find the row and column numbers of a variety of things 990! first, Let's find the first END after GLOBAL 991 CALL mc_parse_text(mc_input_file%text, 1, nlines, "&GLOBAL", .TRUE., & 992 mc_input_file%global_row_end, idum) 993 IF (mc_input_file%global_row_end == 0) THEN 994 IF (iw > 0) THEN 995 WRITE (iw, *) 996 WRITE (iw, *) 'File name ', input_file_name 997 END IF 998 CALL cp_abort(__LOCATION__, & 999 'Could not find &END after &GLOBAL (make sure & is in the same column) and last line is blank') 1000 ENDIF 1001 1002! Let's find the first END after MOTION...we might need this whole section 1003! because of hybrid MD/MC moves, which requires the MD information to always 1004! be attached to the force env 1005 CALL mc_parse_text(mc_input_file%text, 1, nlines, "&MOTION", .TRUE., & 1006 mc_input_file%motion_row_end, idum, start_row_number=mc_input_file%motion_row_start) 1007 IF (mc_input_file%motion_row_end == 0) THEN 1008 IF (iw > 0) THEN 1009 WRITE (iw, *) 1010 WRITE (iw, *) 'File name ', input_file_name, mc_input_file%motion_row_start 1011 END IF 1012 CALL cp_abort(__LOCATION__, & 1013 'Could not find &END after &MOTION (make sure & is in the same column) and last line is blank') 1014 ENDIF 1015 1016! first, let's find the first END after &COORD, and the line of &COORD 1017 CALL mc_parse_text(mc_input_file%text, 1, nlines, "&COORD", .TRUE., & 1018 mc_input_file%coord_row_end, idum) 1019 CALL mc_parse_text(mc_input_file%text, 1, nlines, "&COORD", .FALSE., & 1020 mc_input_file%coord_row_start, idum) 1021 1022 IF (mc_input_file%coord_row_end == 0) THEN 1023 IF (iw > 0) THEN 1024 WRITE (iw, *) 1025 WRITE (iw, *) 'File name ', input_file_name 1026 END IF 1027 CALL cp_abort(__LOCATION__, & 1028 'Could not find &END after &COORD (make sure & is the first in the same column after &COORD)') 1029 ENDIF 1030 1031! now the RUN_TYPE 1032 CALL mc_parse_text(mc_input_file%text, 1, nlines, "RUN_TYPE", .FALSE., & 1033 mc_input_file%run_type_row, mc_input_file%run_type_column) 1034 mc_input_file%run_type_column = mc_input_file%run_type_column + 9 1035 IF (mc_input_file%run_type_row == 0) THEN 1036 IF (iw > 0) THEN 1037 WRITE (iw, *) 1038 WRITE (iw, *) 'File name ', input_file_name 1039 END IF 1040 CALL cp_abort(__LOCATION__, & 1041 'Could not find RUN_TYPE in the input file (should be in the &GLOBAL section).') 1042 ENDIF 1043 1044! now the CELL stuff..we don't care about REF_CELL since that should 1045! never change if it's there 1046 CALL mc_parse_text(mc_input_file%text, 1, nlines, "&CELL", .FALSE., & 1047 mc_input_file%cell_row, mc_input_file%cell_column) 1048! now find the ABC input line after CELL 1049 CALL mc_parse_text(mc_input_file%text, mc_input_file%cell_row + 1, nlines, & 1050 "ABC", .FALSE., abc_row, abc_column) 1051! is there a &CELL inbetween? If so, that ABC will be for the ref_cell 1052! and we need to find the next one 1053 CALL mc_parse_text(mc_input_file%text, mc_input_file%cell_row + 1, abc_row, & 1054 "&CELL", .FALSE., cell_row, cell_column) 1055 IF (cell_row == 0) THEN ! nothing in between...we found the correct ABC 1056 mc_input_file%cell_row = abc_row 1057 mc_input_file%cell_column = abc_column + 4 1058 ELSE 1059 CALL mc_parse_text(mc_input_file%text, abc_row + 1, nlines, & 1060 "ABC", .FALSE., mc_input_file%cell_row, mc_input_file%cell_column) 1061 ENDIF 1062 IF (mc_input_file%cell_row == 0) THEN 1063 IF (iw > 0) THEN 1064 WRITE (iw, *) 1065 WRITE (iw, *) 'File name ', input_file_name 1066 END IF 1067 CALL cp_abort(__LOCATION__, & 1068 'Could not find &CELL section in the input file. Or could not find the ABC line after it.') 1069 ENDIF 1070 1071! now we need all the MOL_SET NMOLS indices 1072 IF (.NOT. lhmc) THEN 1073 irep = 0 1074 nstart = 1 1075 DO 1076 CALL mc_parse_text(mc_input_file%text, nstart, nlines, "&MOLECULE", & 1077 .FALSE., row_number, idum) 1078 IF (row_number == 0) EXIT 1079 nstart = row_number + 1 1080 irep = irep + 1 1081 CALL mc_parse_text(mc_input_file%text, nstart, nlines, "NMOL", & 1082 .FALSE., mc_input_file%mol_set_nmol_row(irep), & 1083 mc_input_file%mol_set_nmol_column(irep)) 1084 mc_input_file%mol_set_nmol_column(irep) = mc_input_file%mol_set_nmol_column(irep) + 5 1085 1086 ENDDO 1087 IF (irep .NE. nmol_types) THEN 1088 IF (iw > 0) THEN 1089 WRITE (iw, *) 1090 WRITE (iw, *) 'File name ', input_file_name 1091 WRITE (iw, *) 'Number of molecule types ', nmol_types 1092 WRITE (iw, *) '&MOLECULE sections found ', irep 1093 END IF 1094 CALL cp_abort( & 1095 __LOCATION__, & 1096 'Did not find MOLECULE sections for every molecule in the simulation (make sure both input files have all types)') 1097 ENDIF 1098 ENDIF 1099 1100! now let's find the type of force_env...right now, I'm only concerned with 1101! QS, and Fist, though it should be trivial to add others 1102 CALL mc_parse_text(mc_input_file%text, 1, nlines, & 1103 "METHOD", .FALSE., line_number, idum) 1104 READ (mc_input_file%text(line_number), *) cdum, method_name 1105 CALL uppercase(method_name) 1106 SELECT CASE (TRIM(ADJUSTL(method_name))) 1107 CASE ("FIST") 1108 mc_input_file%in_use = use_fist_force 1109 CASE ("QS", "QUICKSTEP") 1110 mc_input_file%in_use = use_qs_force 1111 CASE default 1112 CPABORT('Cannot determine the type of force_env we are using (check METHOD)') 1113 END SELECT 1114 1115 END SUBROUTINE mc_input_file_create 1116 1117! ************************************************************************************************** 1118!> \brief destroys (deallocates) things in the mc_input_file_type 1119!> \param mc_input_file the structure that will be released (deallocated) 1120!> \author MJM 1121! ************************************************************************************************** 1122 SUBROUTINE mc_input_file_destroy(mc_input_file) 1123 1124 TYPE(mc_input_file_type), POINTER :: mc_input_file 1125 1126 CHARACTER(len=*), PARAMETER :: routineN = 'mc_input_file_destroy', & 1127 routineP = moduleN//':'//routineN 1128 1129 DEALLOCATE (mc_input_file%mol_set_nmol_row) 1130 DEALLOCATE (mc_input_file%mol_set_nmol_column) 1131 DEALLOCATE (mc_input_file%text) 1132 DEALLOCATE (mc_input_file%atom_names_empty) 1133 DEALLOCATE (mc_input_file%coordinates_empty) 1134 1135 END SUBROUTINE mc_input_file_destroy 1136 1137! ************************************************************************************************** 1138!> \brief a basic text parser used to find the row and column numbers of various strings 1139!> in the input file, to store as indices for when we create a dat_file... 1140!> returns 0 for the row if nothing is found 1141!> \param text the text to parse 1142!> \param nstart the line to start searching from 1143!> \param nend the line to end searching 1144!> \param string_search the text we're looking for 1145!> \param lend if .TRUE., find the &END that comes after string_search...assumes that 1146!> the row is the first row with & in the same column as the search string 1147!> \param row_number the row the string is first found on 1148!> \param column_number the column number that the string first appears on 1149!> \param start_row_number ... 1150!> \author MJM 1151! ************************************************************************************************** 1152 SUBROUTINE mc_parse_text(text, nstart, nend, string_search, lend, & 1153 row_number, column_number, start_row_number) 1154 1155 CHARACTER(LEN=*), DIMENSION(:), INTENT(IN) :: text 1156 INTEGER, INTENT(IN) :: nstart, nend 1157 CHARACTER(LEN=*), INTENT(IN) :: string_search 1158 LOGICAL, INTENT(IN) :: lend 1159 INTEGER, INTENT(OUT) :: row_number, column_number 1160 INTEGER, INTENT(OUT), OPTIONAL :: start_row_number 1161 1162 CHARACTER(default_string_length) :: text_temp 1163 INTEGER :: iline, index_string, jline, string_size 1164 1165! did not see any string utilities in the code to help, so here I go 1166 1167 row_number = 0 1168 column_number = 0 1169 1170! how long is our string to search for? 1171 string_size = LEN_TRIM(string_search) 1172 whole_file: DO iline = nstart, nend 1173 1174 index_string = -1 1175 index_string = INDEX(text(iline), string_search(1:string_size)) 1176 1177 IF (index_string .GT. 0) THEN ! we found one 1178 IF (.NOT. lend) THEN 1179 row_number = iline 1180 column_number = index_string 1181 ELSE 1182 IF (PRESENT(start_row_number)) start_row_number = iline 1183 column_number = index_string 1184 DO jline = iline + 1, nend 1185! now we find the &END that matches up with this one... 1186! I need proper indentation because I'm not very smart 1187 text_temp = text(jline) 1188 IF (text_temp(column_number:column_number) .EQ. "&") THEN 1189 row_number = jline 1190 EXIT 1191 ENDIF 1192 ENDDO 1193 ENDIF 1194 EXIT whole_file 1195 ENDIF 1196 ENDDO whole_file 1197 1198 END SUBROUTINE mc_parse_text 1199 1200! ************************************************************************************************** 1201!> \brief reads in the Monte Carlo simulation parameters from an input file 1202!> \param mc_par the structure that will store the parameters 1203!> \param para_env ... 1204!> \param globenv the global environment for the simulation 1205!> \param input_file_name the name of the input_file 1206!> \param input_file the structure that contains all the keywords in the input file 1207!> \param force_env_section used to grab the type of force_env 1208!> \author MJM 1209! ************************************************************************************************** 1210 SUBROUTINE read_mc_section(mc_par, para_env, globenv, input_file_name, input_file, force_env_section) 1211 1212 TYPE(mc_simpar_type), POINTER :: mc_par 1213 TYPE(cp_para_env_type), POINTER :: para_env 1214 TYPE(global_environment_type), POINTER :: globenv 1215 CHARACTER(LEN=*), INTENT(IN) :: input_file_name 1216 TYPE(section_vals_type), POINTER :: input_file, force_env_section 1217 1218 CHARACTER(len=*), PARAMETER :: routineN = 'read_mc_section', & 1219 routineP = moduleN//':'//routineN 1220 1221 CHARACTER(LEN=5) :: box_string, molecule_string, & 1222 tab_box_string, tab_string, & 1223 tab_string_int 1224 CHARACTER(LEN=default_string_length) :: format_box_string, format_string, & 1225 format_string_int 1226 INTEGER :: handle, ia, ie, imol, itype, iw, & 1227 method_name_id, nmol_types, stop_num 1228 INTEGER, DIMENSION(:), POINTER :: temp_i_array 1229 REAL(dp), DIMENSION(:), POINTER :: temp_r_array 1230 TYPE(section_vals_type), POINTER :: mc_section 1231 1232! begin the timing of the subroutine 1233 1234 CPASSERT(ASSOCIATED(input_file)) 1235 1236 CALL timeset(routineN, handle) 1237 1238 NULLIFY (mc_section) 1239 mc_section => section_vals_get_subs_vals(input_file, & 1240 "MOTION%MC") 1241 1242! need the input file sturcutre that we're reading from for when we make 1243! dat files 1244 mc_par%input_file => input_file 1245 1246! set the ionode and mepos 1247 mc_par%ionode = para_env%ionode 1248 mc_par%group = para_env%group 1249 mc_par%source = para_env%source 1250 1251! for convenience 1252 nmol_types = mc_par%mc_molecule_info%nmol_types 1253 1254!..defaults...most are set in input_cp2k_motion 1255 mc_par%nstart = 0 1256 CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id) 1257 SELECT CASE (method_name_id) 1258 CASE (do_fist) 1259 mc_par%iprint = 100 1260 CASE (do_qs) 1261 mc_par%iprint = 1 1262 END SELECT 1263 1264 iw = cp_logger_get_default_io_unit() 1265 IF (iw > 0) WRITE (iw, *) 1266 1267!..filenames 1268 mc_par%program = input_file_name 1269 CALL xstring(mc_par%program, ia, ie) 1270 mc_par%coords_file = mc_par%program(ia:ie)//'.coordinates' 1271 mc_par%molecules_file = mc_par%program(ia:ie)//'.molecules' 1272 mc_par%moves_file = mc_par%program(ia:ie)//'.moves' 1273 mc_par%energy_file = mc_par%program(ia:ie)//'.energy' 1274 mc_par%cell_file = mc_par%program(ia:ie)//'.cell' 1275 mc_par%displacement_file = mc_par%program(ia:ie) & 1276 //'.max_displacements' 1277 mc_par%data_file = mc_par%program(ia:ie)//'.data' 1278 stop_num = ie - 3 1279 mc_par%dat_file = mc_par%program(ia:stop_num)//'dat' 1280 1281! set them into the input parameter structure as the new defaults 1282 CALL section_vals_val_set(mc_section, "COORDINATE_FILE_NAME", & 1283 c_val=mc_par%coords_file) 1284 CALL section_vals_val_set(mc_section, "DATA_FILE_NAME", & 1285 c_val=mc_par%data_file) 1286 CALL section_vals_val_set(mc_section, "CELL_FILE_NAME", & 1287 c_val=mc_par%cell_file) 1288 CALL section_vals_val_set(mc_section, "MAX_DISP_FILE_NAME", & 1289 c_val=mc_par%displacement_file) 1290 CALL section_vals_val_set(mc_section, "MOVES_FILE_NAME", & 1291 c_val=mc_par%moves_file) 1292 CALL section_vals_val_set(mc_section, "MOLECULES_FILE_NAME", & 1293 c_val=mc_par%molecules_file) 1294 CALL section_vals_val_set(mc_section, "ENERGY_FILE_NAME", & 1295 c_val=mc_par%energy_file) 1296 1297! grab the FFT library name and print level...this is used for writing the dat file 1298! and hopefully will be changed 1299 mc_par%fft_lib = globenv%default_fft_library 1300 1301! get all the move probabilities first...if we are not doing certain types of moves, we don't care 1302! if the values for those move parameters are strange 1303 1304! find out if we're only doing HMC moves...we can ignore a lot of extra information 1305! then, which would ordinarily be cause for concern 1306 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMHMC", r_val=mc_par%pmhmc) 1307 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMSWAP", r_val=mc_par%pmswap) 1308 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMVOLUME", r_val=mc_par%pmvolume) 1309 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMAVBMC", r_val=mc_par%pmavbmc) 1310 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMTRANS", r_val=mc_par%pmtrans) 1311 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMTRAION", r_val=mc_par%pmtraion) 1312 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%PMCLTRANS", r_val=mc_par%pmcltrans) 1313 1314 ! first, grab all the integer values 1315 CALL section_vals_val_get(mc_section, "NSTEP", i_val=mc_par%nstep) 1316 CALL section_vals_val_get(mc_section, "NMOVES", i_val=mc_par%nmoves) 1317 CALL section_vals_val_get(mc_section, "NSWAPMOVES", i_val=mc_par%nswapmoves) 1318 CALL section_vals_val_get(mc_section, "MOVE_UPDATES%IUPVOLUME", & 1319 i_val=mc_par%iupvolume) 1320 CALL section_vals_val_get(mc_section, "NVIRIAL", i_val=mc_par%nvirial) 1321 CALL section_vals_val_get(mc_section, "MOVE_UPDATES%IUPTRANS", & 1322 i_val=mc_par%iuptrans) 1323 CALL section_vals_val_get(mc_section, "MOVE_UPDATES%IUPCLTRANS", & 1324 i_val=mc_par%iupcltrans) 1325 CALL section_vals_val_get(mc_section, "IPRINT", i_val=mc_par%iprint) 1326! now an integer array 1327 CALL section_vals_val_get(mc_section, "AVBMC%AVBMC_ATOM", i_vals=temp_i_array) 1328 1329 IF (mc_par%pmhmc - mc_par%pmswap >= 1.0_dp .AND. mc_par%pmhmc == 1.0_dp) THEN 1330 mc_par%lhmc = .TRUE. 1331 ELSE 1332 mc_par%lhmc = .FALSE. 1333 ENDIF 1334 1335 IF (.NOT. mc_par%lhmc) THEN 1336 IF (SIZE(temp_i_array) .NE. nmol_types) THEN 1337 CPABORT('AVBMC_ATOM must have as many values as there are molecule types.') 1338 ENDIF 1339 ENDIF 1340 DO imol = 1, SIZE(temp_i_array) 1341 mc_par%avbmc_atom(imol) = temp_i_array(imol) 1342 ENDDO 1343 1344! now the real values 1345 CALL section_vals_val_get(mc_section, "PRESSURE", r_val=mc_par%pressure) 1346 CALL section_vals_val_get(mc_section, "TEMPERATURE", r_val=mc_par%temperature) 1347 CALL section_vals_val_get(mc_section, "RCLUS", r_val=mc_par%rclus) 1348 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%BOX_PROBABILITIES%PMCLUS_BOX", & 1349 r_val=mc_par%pmclus_box) 1350 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%BOX_PROBABILITIES%PMHMC_BOX", & 1351 r_val=mc_par%pmhmc_box) 1352 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%BOX_PROBABILITIES%PMVOL_BOX", & 1353 r_val=mc_par%pmvol_box) 1354 CALL section_vals_val_get(mc_section, "DISCRETE_STEP", r_val=mc_par%discrete_step) 1355 CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%BOX_DISPLACEMENTS%RMVOLUME", & 1356 r_val=mc_par%rmvolume) 1357 CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%BOX_DISPLACEMENTS%RMCLTRANS", & 1358 r_val=mc_par%rmcltrans) 1359 1360! finally the character values 1361 CALL section_vals_val_get(mc_section, "ENSEMBLE", c_val=mc_par%ensemble) 1362 CALL section_vals_val_get(mc_section, "RESTART_FILE_NAME", c_val=mc_par%restart_file_name) 1363 CALL section_vals_val_get(mc_section, "COORDINATE_FILE_NAME", c_val=mc_par%coords_file) 1364 CALL section_vals_val_get(mc_section, "ENERGY_FILE_NAME", c_val=mc_par%energy_file) 1365 CALL section_vals_val_get(mc_section, "MOVES_FILE_NAME", c_val=mc_par%moves_file) 1366 CALL section_vals_val_get(mc_section, "MOLECULES_FILE_NAME", c_val=mc_par%molecules_file) 1367 CALL section_vals_val_get(mc_section, "CELL_FILE_NAME", c_val=mc_par%cell_file) 1368 CALL section_vals_val_get(mc_section, "DATA_FILE_NAME", c_val=mc_par%data_file) 1369 CALL section_vals_val_get(mc_section, "MAX_DISP_FILE_NAME", c_val=mc_par%displacement_file) 1370 CALL section_vals_val_get(mc_section, "BOX2_FILE_NAME", c_val=mc_par%box2_file) 1371 1372! set the values of the arrays...if we just point, we have problems when we start fooling around 1373! with releasing force_envs and wonky values get set (despite that these are private) 1374 IF (mc_par%ensemble == "VIRIAL") THEN 1375 CALL section_vals_val_get(mc_section, "VIRIAL_TEMPS", r_vals=temp_r_array) 1376! yes, I'm allocating here...I cannot find a better place to do it, though 1377 ALLOCATE (mc_par%virial_temps(1:SIZE(temp_r_array))) 1378 DO imol = 1, SIZE(temp_r_array) 1379 mc_par%virial_temps(imol) = temp_r_array(imol) 1380 ENDDO 1381 END IF 1382! all of these arrays should have one value for each type of molecule...so check that 1383 CALL section_vals_val_get(mc_section, "AVBMC%AVBMC_RMIN", r_vals=temp_r_array) 1384 IF (.NOT. mc_par%lhmc) THEN 1385 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1386 CPABORT('AVBMC_RMIN must have as many values as there are molecule types.') 1387 ENDIF 1388 ENDIF 1389 DO imol = 1, SIZE(temp_r_array) 1390 mc_par%avbmc_rmin(imol) = temp_r_array(imol) 1391 ENDDO 1392 CALL section_vals_val_get(mc_section, "AVBMC%AVBMC_RMAX", r_vals=temp_r_array) 1393 IF (.NOT. mc_par%lhmc) THEN 1394 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1395 CPABORT('AVBMC_RMAXL must have as many values as there are molecule types.') 1396 ENDIF 1397 ENDIF 1398 DO imol = 1, SIZE(temp_r_array) 1399 mc_par%avbmc_rmax(imol) = temp_r_array(imol) 1400 ENDDO 1401 CALL section_vals_val_get(mc_section, "AVBMC%PBIAS", r_vals=temp_r_array) 1402 IF (.NOT. mc_par%lhmc) THEN 1403 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1404 CPABORT('PBIAS must have as many values as there are molecule types.') 1405 ENDIF 1406 ENDIF 1407 DO imol = 1, SIZE(temp_r_array) 1408 mc_par%pbias(imol) = temp_r_array(imol) 1409 ENDDO 1410 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%MOL_PROBABILITIES%PMAVBMC_MOL", & 1411 r_vals=temp_r_array) 1412 IF (.NOT. mc_par%lhmc) THEN 1413 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1414 CPABORT('PMAVBMC_MOL must have as many values as there are molecule types.') 1415 ENDIF 1416 ENDIF 1417 DO imol = 1, SIZE(temp_r_array) 1418 mc_par%pmavbmc_mol(imol) = temp_r_array(imol) 1419 ENDDO 1420 CALL section_vals_val_get(mc_section, "ETA", r_vals=temp_r_array) 1421 IF (.NOT. mc_par%lhmc) THEN 1422 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1423 CPABORT('ETA must have as many values as there are molecule types.') 1424 ENDIF 1425 ENDIF 1426 DO imol = 1, SIZE(temp_r_array) 1427 mc_par%eta(imol) = temp_r_array(imol) 1428 ENDDO 1429 CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%MOL_DISPLACEMENTS%RMBOND", & 1430 r_vals=temp_r_array) 1431 IF (.NOT. mc_par%lhmc) THEN 1432 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1433 CPABORT('RMBOND must have as many values as there are molecule types.') 1434 ENDIF 1435 ENDIF 1436 DO imol = 1, SIZE(temp_r_array) 1437 mc_par%rmbond(imol) = temp_r_array(imol) 1438 ENDDO 1439 CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%MOL_DISPLACEMENTS%RMANGLE", & 1440 r_vals=temp_r_array) 1441 IF (.NOT. mc_par%lhmc) THEN 1442 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1443 CPABORT('RMANGLE must have as many values as there are molecule types.') 1444 ENDIF 1445 ENDIF 1446 DO imol = 1, SIZE(temp_r_array) 1447 mc_par%rmangle(imol) = temp_r_array(imol) 1448 ENDDO 1449 CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%MOL_DISPLACEMENTS%RMDIHEDRAL", & 1450 r_vals=temp_r_array) 1451 IF (.NOT. mc_par%lhmc) THEN 1452 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1453 CPABORT('RMDIHEDRAL must have as many values as there are molecule types.') 1454 ENDIF 1455 ENDIF 1456 DO imol = 1, SIZE(temp_r_array) 1457 mc_par%rmdihedral(imol) = temp_r_array(imol) 1458 ENDDO 1459 CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%MOL_DISPLACEMENTS%RMROT", & 1460 r_vals=temp_r_array) 1461 IF (.NOT. mc_par%lhmc) THEN 1462 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1463 CPABORT('RMROT must have as many values as there are molecule types.') 1464 ENDIF 1465 ENDIF 1466 DO imol = 1, SIZE(temp_r_array) 1467 mc_par%rmrot(imol) = temp_r_array(imol) 1468 ENDDO 1469 CALL section_vals_val_get(mc_section, "MAX_DISPLACEMENTS%MOL_DISPLACEMENTS%RMTRANS", & 1470 r_vals=temp_r_array) 1471 IF (.NOT. mc_par%lhmc) THEN 1472 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1473 CPABORT('RMTRANS must have as many values as there are molecule types.') 1474 ENDIF 1475 ENDIF 1476 DO imol = 1, SIZE(temp_r_array) 1477 mc_par%rmtrans(imol) = temp_r_array(imol) 1478 ENDDO 1479 1480 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%MOL_PROBABILITIES%PMTRAION_MOL", & 1481 r_vals=temp_r_array) 1482 IF (.NOT. mc_par%lhmc) THEN 1483 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1484 CPABORT('PMTRAION_MOL must have as many values as there are molecule types.') 1485 ENDIF 1486 ENDIF 1487 DO imol = 1, SIZE(temp_r_array) 1488 mc_par%pmtraion_mol(imol) = temp_r_array(imol) 1489 ENDDO 1490 1491 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%MOL_PROBABILITIES%PMTRANS_MOL", & 1492 r_vals=temp_r_array) 1493 IF (.NOT. mc_par%lhmc) THEN 1494 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1495 CPABORT('PMTRANS_MOL must have as many values as there are molecule types.') 1496 ENDIF 1497 ENDIF 1498 DO imol = 1, SIZE(temp_r_array) 1499 mc_par%pmtrans_mol(imol) = temp_r_array(imol) 1500 ENDDO 1501 1502 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%MOL_PROBABILITIES%PMROT_MOL", & 1503 r_vals=temp_r_array) 1504 IF (.NOT. mc_par%lhmc) THEN 1505 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1506 CPABORT('PMROT_MOL must have as many values as there are molecule types.') 1507 ENDIF 1508 ENDIF 1509 DO imol = 1, SIZE(temp_r_array) 1510 mc_par%pmrot_mol(imol) = temp_r_array(imol) 1511 ENDDO 1512 1513 CALL section_vals_val_get(mc_section, "MOVE_PROBABILITIES%MOL_PROBABILITIES%PMSWAP_MOL", & 1514 r_vals=temp_r_array) 1515 IF (.NOT. mc_par%lhmc) THEN 1516 IF (SIZE(temp_r_array) .NE. nmol_types) THEN 1517 CPABORT('PMSWAP_MOL must have as many values as there are molecule types.') 1518 ENDIF 1519 ENDIF 1520 DO imol = 1, SIZE(temp_r_array) 1521 mc_par%pmswap_mol(imol) = temp_r_array(imol) 1522 ENDDO 1523 1524! now some logical values 1525 CALL section_vals_val_get(mc_section, "LBIAS", l_val=mc_par%lbias) 1526 CALL section_vals_val_get(mc_section, "LDISCRETE", l_val=mc_par%ldiscrete) 1527 CALL section_vals_val_get(mc_section, "LSTOP", l_val=mc_par%lstop) 1528 CALL section_vals_val_get(mc_section, "RESTART", l_val=mc_par%lrestart) 1529 CALL section_vals_val_get(mc_section, "LBIAS", l_val=mc_par%lbias) 1530 1531!..end of parsing the input section 1532 1533!..write some information to output 1534 IF (iw > 0) THEN 1535 WRITE (box_string, '(I2)') mc_par%mc_molecule_info%nboxes 1536 WRITE (molecule_string, '(I2)') mc_par%mc_molecule_info%nmol_types 1537 WRITE (tab_string, '(I4)') 81 - 10*mc_par%mc_molecule_info%nmol_types 1538 WRITE (tab_box_string, '(I4)') 81 - 10*mc_par%mc_molecule_info%nboxes 1539 WRITE (tab_string_int, '(I4)') 81 - 5*mc_par%mc_molecule_info%nmol_types 1540 format_string = "(A,T"//TRIM(ADJUSTL(tab_string))//","//TRIM(ADJUSTL(molecule_string))//"(2X,F8.4))" 1541 format_box_string = "(A,T"//TRIM(ADJUSTL(tab_box_string))//","//TRIM(ADJUSTL(box_string))//"(2X,F8.4))" 1542 format_string_int = "(A,T"//TRIM(ADJUSTL(tab_string))//","//TRIM(ADJUSTL(molecule_string))//"(I3,2X))" 1543 WRITE (iw, '( A )') ' MC| Monte Carlo Protocol ' 1544 WRITE (iw, '( A,T71,I10 )') ' MC| total number of steps ', & 1545 mc_par%nstep 1546 WRITE (iw, '( A,T71,F10.3 )') ' MC| pmvolume ', & 1547 mc_par%pmvolume 1548 WRITE (iw, '( A,T71,F10.3 )') ' MC| pmvol_box ', & 1549 mc_par%pmvol_box 1550 WRITE (iw, '( A,T71,F10.3 )') ' MC| pmclus_box ', & 1551 mc_par%pmclus_box 1552 WRITE (iw, '( A,T71,F10.3 )') ' MC| pmhmc ', & 1553 mc_par%pmhmc 1554 WRITE (iw, '( A,T71,F10.3 )') ' MC| pmhmc_box ', & 1555 mc_par%pmhmc_box 1556 WRITE (iw, '( A,T71,F10.3 )') ' MC| pmswap ', & 1557 mc_par%pmswap 1558 WRITE (iw, '( A,T71,F10.3 )') ' MC| pmavbmc ', & 1559 mc_par%pmavbmc 1560 WRITE (iw, '( A,T71,F10.3 )') ' MC| pmtraion ', & 1561 mc_par%pmtraion 1562 WRITE (iw, '( A,T71,F10.3 )') ' MC| pmtrans ', & 1563 mc_par%pmtrans 1564 WRITE (iw, '( A,T71,F10.3 )') ' MC| pmcltrans ', & 1565 mc_par%pmcltrans 1566 WRITE (iw, '( A,T71,I10 )') ' MC| iupvolume ', & 1567 mc_par%iupvolume 1568 WRITE (iw, '( A,T71,I10 )') ' MC| nvirial ', & 1569 mc_par%nvirial 1570 WRITE (iw, '( A,T71,I10 )') ' MC| iuptrans ', & 1571 mc_par%iuptrans 1572 WRITE (iw, '( A,T71,I10 )') ' MC| iupcltrans ', & 1573 mc_par%iupcltrans 1574 WRITE (iw, '( A,T71,I10 )') ' MC| iprint ', & 1575 mc_par%iprint 1576 WRITE (iw, '( A,T61,A20 )') ' MC| ensemble ', & 1577 ADJUSTR(mc_par%ensemble) 1578! format string may not be valid if all the molecules are atomic, 1579! like in HMC 1580 IF (.NOT. mc_par%lhmc) THEN 1581 WRITE (iw, format_string) ' MC| pmswap_mol ', & 1582 mc_par%pmswap_mol 1583 WRITE (iw, format_string) ' MC| pmavbmc_mol ', & 1584 mc_par%pmavbmc_mol 1585 WRITE (iw, format_string) ' MC| pbias ', & 1586 mc_par%pbias 1587 WRITE (iw, format_string) ' MC| pmtraion_mol ', & 1588 mc_par%pmtraion_mol 1589 WRITE (iw, format_string) ' MC| pmtrans_mol ', & 1590 mc_par%pmtrans_mol 1591 WRITE (iw, format_string) ' MC| pmrot_mol ', & 1592 mc_par%pmrot_mol 1593 WRITE (iw, format_string) ' MC| eta [K]', & 1594 mc_par%eta 1595 WRITE (iw, format_string) ' MC| rmbond [angstroms]', & 1596 mc_par%rmbond 1597 WRITE (iw, format_string) ' MC| rmangle [degrees]', & 1598 mc_par%rmangle 1599 WRITE (iw, format_string) ' MC| rmdihedral [degrees]', & 1600 mc_par%rmdihedral 1601 WRITE (iw, format_string) ' MC| rmtrans [angstroms]', & 1602 mc_par%rmtrans 1603 WRITE (iw, format_string) ' MC| rmcltrans [angstroms]', & 1604 mc_par%rmcltrans 1605 WRITE (iw, format_string) ' MC| rmrot [degrees]', & 1606 mc_par%rmrot 1607 WRITE (iw, format_string_int) ' MC| AVBMC target atom ', & 1608 mc_par%avbmc_atom 1609 WRITE (iw, format_string) ' MC| AVBMC inner cutoff [ang]', & 1610 mc_par%avbmc_rmin 1611 WRITE (iw, format_string) ' MC| AVBMC outer cutoff [ang]', & 1612 mc_par%avbmc_rmax 1613 ENDIF 1614 IF (mc_par%ensemble .EQ. 'GEMC-NVT') THEN 1615 WRITE (iw, '( A,T58,A)') ' MC| Box 2 file', & 1616 TRIM(mc_par%box2_file) 1617 ENDIF 1618 WRITE (iw, '( A,T58,A )') ' MC| Name of restart file:', & 1619 TRIM(mc_par%restart_file_name) 1620 WRITE (iw, '( A,A,T44,A )') ' MC| Name of output ', & 1621 'coordinate file:', & 1622 TRIM(mc_par%coords_file) 1623 WRITE (iw, '( A,T44,A )') ' MC| Name of output data file:', & 1624 TRIM(mc_par%data_file) 1625 WRITE (iw, '( A,A,T44,A )') ' MC| Name of output ', & 1626 'molecules file:', & 1627 TRIM(mc_par%molecules_file) 1628 WRITE (iw, '( A,T44,A )') ' MC| Name of output moves file:', & 1629 TRIM(mc_par%moves_file) 1630 WRITE (iw, '( A,T44,A )') ' MC| Name of output energy file:', & 1631 TRIM(mc_par%energy_file) 1632 WRITE (iw, '( A,T44,A )') ' MC| Name of output cell file:', & 1633 TRIM(mc_par%cell_file) 1634 WRITE (iw, '( A,A,T44,A )') ' MC| Name of output', & 1635 ' displacement file:', & 1636 TRIM(mc_par%displacement_file) 1637 IF (mc_par%ldiscrete) THEN 1638 WRITE (iw, '(A,A,T71,F10.3)') ' MC| discrete step size', & 1639 '[angstroms]', & 1640 mc_par%discrete_step 1641 ELSE 1642 WRITE (iw, '( A,A,T71,F10.3 )') ' MC| rmvolume ', & 1643 '[cubic angstroms]', & 1644 mc_par%rmvolume 1645 ENDIF 1646 WRITE (iw, '( A,T71,F10.2 )') ' MC| Temperature [K] ', & 1647 mc_par%temperature 1648 WRITE (iw, '( A,T71,F10.5 )') ' MC| Pressure [bar] ', & 1649 mc_par%pressure 1650 WRITE (iw, '( A,T71,F10.5 )') ' MC| Rclus [ang] ', & 1651 mc_par%rclus 1652 IF (mc_par%lrestart) THEN 1653 WRITE (iw, '(A,A)') ' MC| Initial data will be read from a', & 1654 ' restart file.' 1655 ENDIF 1656 IF (mc_par%lbias) THEN 1657 WRITE (iw, '(A)') ' MC| The moves will be biased.' 1658 ELSE 1659 WRITE (iw, '(A)') ' MC| The moves will not be biased,' 1660 ENDIF 1661 IF (mc_par%nmoves .EQ. 1) THEN 1662 WRITE (iw, '(A,A)') ' MC| A full energy calculation ', & 1663 'will be done at every step.' 1664 ELSE 1665 WRITE (iw, '(A,I4,A,A)') ' MC| ', mc_par%nmoves, & 1666 ' moves will be attempted ', & 1667 'before a Quickstep energy calculation' 1668 WRITE (iw, '(A)') ' MC| takes place.' 1669 ENDIF 1670 1671 WRITE (iw, '(A,I4,A,A)') ' MC| ', mc_par%nswapmoves, & 1672 ' swap insertions will be attempted ', & 1673 'per molecular swap move' 1674 1675 IF (mc_par%lhmc) THEN 1676 WRITE (iw, '(A)') '********************************************************' 1677 WRITE (iw, '(A)') '************** ONLY DOING HYBRID MONTE CARLO MOVES **************************' 1678 WRITE (iw, '(A)') '********************************************************' 1679 1680 ENDIF 1681 END IF 1682 1683! figure out what beta (1/kT) is in atomic units (1/Hartree) 1684 mc_par%BETA = 1/mc_par%temperature/boltzmann*joule 1685 1686 ! 0.9_dp is to avoid overflow or underflow 1687 mc_par%exp_max_val = 0.9_dp*LOG(HUGE(0.0_dp)) 1688 mc_par%exp_min_val = 0.9_dp*LOG(TINY(0.0_dp)) 1689 mc_par%max_val = EXP(mc_par%exp_max_val) 1690 mc_par%min_val = 0.0_dp 1691 1692! convert from bar to a.u. 1693 mc_par%pressure = cp_unit_to_cp2k(mc_par%pressure, "bar") 1694 mc_par%rclus = cp_unit_to_cp2k(mc_par%rclus, "angstrom") 1695! convert from angstrom to a.u. 1696 DO itype = 1, mc_par%mc_molecule_info%nmol_types 1697! convert from Kelvin to a.u. 1698! mc_par % eta = mc_par%eta(itype) * boltzmann / joule 1699! convert from degrees to radians 1700 mc_par%rmrot(itype) = mc_par%rmrot(itype)/180.0e0_dp*pi 1701 mc_par%rmangle(itype) = mc_par%rmangle(itype)/180.0e0_dp*pi 1702 mc_par%rmdihedral(itype) = mc_par%rmdihedral(itype)/180.0e0_dp*pi 1703 mc_par%rmtrans(itype) = cp_unit_to_cp2k(mc_par%rmtrans(itype), "angstrom") 1704 mc_par%rmbond(itype) = cp_unit_to_cp2k(mc_par%rmbond(itype), "angstrom") 1705 mc_par%eta(itype) = cp_unit_to_cp2k(mc_par%eta(itype), "K_e") 1706 mc_par%avbmc_rmin(itype) = cp_unit_to_cp2k(mc_par%avbmc_rmin(itype), "angstrom") 1707 mc_par%avbmc_rmax(itype) = cp_unit_to_cp2k(mc_par%avbmc_rmax(itype), "angstrom") 1708 ENDDO 1709 mc_par%rmvolume = cp_unit_to_cp2k(mc_par%rmvolume, "angstrom^3") 1710 mc_par%rmcltrans = cp_unit_to_cp2k(mc_par%rmcltrans, "angstrom") 1711 mc_par%discrete_step = cp_unit_to_cp2k(mc_par%discrete_step, "angstrom") 1712 1713! end the timing 1714 CALL timestop(handle) 1715 1716 END SUBROUTINE read_mc_section 1717 1718! ************************************************************************************************** 1719!> \brief finds the largest interaction cutoff value in a classical simulation 1720!> so we know the smallest size we can make the box in a volume move 1721!> \param mc_par the structure that will store the parameters 1722!> \param force_env the force environment that we'll grab the rcut parameter 1723!> out of 1724!> \param lterminate set to .TRUE. if one of the sides of the box is 1725!> less than twice the cutoff 1726!> 1727!> Suitable for parallel. 1728!> \author MJM 1729! ************************************************************************************************** 1730 SUBROUTINE find_mc_rcut(mc_par, force_env, lterminate) 1731 1732 TYPE(mc_simpar_type), INTENT(INOUT) :: mc_par 1733 TYPE(force_env_type), POINTER :: force_env 1734 LOGICAL, INTENT(OUT) :: lterminate 1735 1736 CHARACTER(len=*), PARAMETER :: routineN = 'find_mc_rcut', routineP = moduleN//':'//routineN 1737 1738 INTEGER :: itype, jtype 1739 REAL(KIND=dp) :: rcutsq_max 1740 REAL(KIND=dp), DIMENSION(1:3) :: abc 1741 TYPE(cell_type), POINTER :: cell 1742 TYPE(fist_environment_type), POINTER :: fist_env 1743 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 1744 TYPE(pair_potential_pp_type), POINTER :: potparm 1745 1746 NULLIFY (cell, potparm, fist_nonbond_env, fist_env) 1747 1748 lterminate = .FALSE. 1749 CALL force_env_get(force_env, fist_env=fist_env) 1750 CALL fist_env_get(fist_env, cell=cell, & 1751 fist_nonbond_env=fist_nonbond_env) 1752 CALL fist_nonbond_env_get(fist_nonbond_env, potparm=potparm) 1753 CALL get_cell(cell, abc=abc) 1754 1755! find the largest value of rcutsq 1756 rcutsq_max = 0.0e0_dp 1757 DO itype = 1, SIZE(potparm%pot, 1) 1758 DO jtype = itype, SIZE(potparm%pot, 2) 1759 IF (potparm%pot(itype, jtype)%pot%rcutsq .GT. rcutsq_max) & 1760 rcutsq_max = potparm%pot(itype, jtype)%pot%rcutsq 1761 ENDDO 1762 ENDDO 1763 1764! check to make sure all box dimensions are greater than two times this 1765! value 1766 mc_par%rcut = rcutsq_max**0.5_dp 1767 DO itype = 1, 3 1768 IF (abc(itype) .LT. 2.0_dp*mc_par%rcut) THEN 1769 lterminate = .TRUE. 1770 ENDIF 1771 ENDDO 1772 1773 END SUBROUTINE find_mc_rcut 1774 1775! ************************************************************************************************** 1776!> \brief figures out the number of total molecules, the number of atoms in each 1777!> molecule, an array with the molecule types, etc...a lot of information 1778!> that we need. I did this because I use multiple force_env (simulation 1779!> boxes) for MC, and they don't know about each other. 1780!> \param force_env the pointer containing all the force environments in the 1781!> simulation 1782!> \param mc_molecule_info the structure that will hold all the information 1783!> for the molecule types 1784!> 1785!> Suitable for parallel. 1786!> \param box_number ... 1787!> \param coordinates_empty ... 1788!> \author MJM 1789! ************************************************************************************************** 1790 SUBROUTINE mc_determine_molecule_info(force_env, mc_molecule_info, box_number, & 1791 coordinates_empty) 1792 1793 TYPE(force_env_p_type), DIMENSION(:), POINTER :: force_env 1794 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 1795 INTEGER, INTENT(IN), OPTIONAL :: box_number 1796 REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: coordinates_empty 1797 1798 CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_determine_molecule_info', & 1799 routineP = moduleN//':'//routineN 1800 1801 CHARACTER(LEN=default_string_length) :: name 1802 CHARACTER(LEN=default_string_length), & 1803 ALLOCATABLE, DIMENSION(:) :: names_init 1804 INTEGER :: first_mol, iatom, ibox, imol, imolecule, itype, iunique, iunit, jtype, last_mol, & 1805 natom, natoms_large, nbend, nbond, nboxes, nmol_types, nmolecule, ntorsion, ntypes, & 1806 skip_box, this_molecule, total 1807 INTEGER, DIMENSION(:), POINTER :: molecule_list 1808 LOGICAL :: lnew_type 1809 TYPE(atom_type), DIMENSION(:), POINTER :: atom_list 1810 TYPE(atomic_kind_type), POINTER :: atomic_kind 1811 TYPE(cp_subsys_type), POINTER :: subsys 1812 TYPE(molecule_kind_list_p_type), DIMENSION(:), & 1813 POINTER :: molecule_kinds 1814 TYPE(molecule_kind_type), POINTER :: molecule_kind 1815 TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles 1816 1817! how many simulation boxes do we have? 1818 1819 nboxes = SIZE(force_env(:)) 1820 1821! allocate the pointer 1822 ALLOCATE (mc_molecule_info) 1823 mc_molecule_info%nboxes = nboxes 1824 1825! if box_number is present, that box is supposed to be empty, 1826! so we don't count it in anything 1827 IF (PRESENT(box_number)) THEN 1828 skip_box = box_number 1829 ELSE 1830 skip_box = 0 1831 ENDIF 1832 1833! we need to figure out how many different kinds of molecules we have... 1834! the fun stuff comes in when box 1 is missing a molecule type...I'll 1835! distinguish molecules based on names from the .psf files... 1836! this is only really a problem when we have more than one simulation 1837! box 1838 NULLIFY (subsys, molecule_kinds, molecule_kind) 1839 1840 ALLOCATE (particles(1:nboxes)) 1841 ALLOCATE (molecule_kinds(1:nboxes)) 1842 1843 ! find out how many types we have 1844 ntypes = 0 1845 DO ibox = 1, nboxes 1846 IF (ibox == skip_box) CYCLE 1847 CALL force_env_get(force_env(ibox)%force_env, & 1848 subsys=subsys) 1849 CALL cp_subsys_get(subsys, & 1850 molecule_kinds=molecule_kinds(ibox)%list, & 1851 particles=particles(ibox)%list) 1852 ntypes = ntypes + SIZE(molecule_kinds(ibox)%list%els(:)) 1853 ENDDO 1854 1855 ALLOCATE (names_init(1:ntypes)) 1856 1857! now get the names of all types so we can figure out how many distinct 1858! types we have 1859 itype = 1 1860 natoms_large = 0 1861 DO ibox = 1, nboxes 1862 IF (ibox == skip_box) CYCLE 1863 DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:)) 1864 molecule_kind => molecule_kinds(ibox)%list%els(imolecule) 1865 CALL get_molecule_kind(molecule_kind, name=names_init(itype), & 1866 natom=natom) 1867 IF (natom .GT. natoms_large) natoms_large = natom 1868 itype = itype + 1 1869 ENDDO 1870 ENDDO 1871 1872 nmol_types = 0 1873 DO itype = 1, ntypes 1874 lnew_type = .TRUE. 1875 DO jtype = 1, itype - 1 1876 IF (TRIM(names_init(itype)) .EQ. TRIM(names_init(jtype))) & 1877 lnew_type = .FALSE. 1878 ENDDO 1879 IF (lnew_type) THEN 1880 nmol_types = nmol_types + 1 1881 ELSE 1882 names_init(itype) = '' 1883 ENDIF 1884 ENDDO 1885 1886! allocate arrays for the names of the molecules, the number of atoms, and 1887! the conformational probabilities 1888 mc_molecule_info%nmol_types = nmol_types 1889 ALLOCATE (mc_molecule_info%names(1:nmol_types)) 1890 ALLOCATE (mc_molecule_info%nunits(1:nmol_types)) 1891 ALLOCATE (mc_molecule_info%conf_prob(1:3, 1:nmol_types)) 1892 ALLOCATE (mc_molecule_info%nchains(1:nmol_types, 1:nboxes)) 1893 ALLOCATE (mc_molecule_info%nunits_tot(1:nboxes)) 1894 ALLOCATE (mc_molecule_info%atom_names(1:natoms_large, 1:nmol_types)) 1895 ALLOCATE (mc_molecule_info%mass(1:natoms_large, 1:nmol_types)) 1896 1897! now record all the information for every molecule type 1898 iunique = 0 1899 itype = 0 1900 DO ibox = 1, nboxes 1901 IF (ibox == skip_box) CYCLE 1902 DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:)) 1903 itype = itype + 1 1904 IF (names_init(itype) .EQ. '') CYCLE 1905 iunique = iunique + 1 1906 mc_molecule_info%names(iunique) = names_init(itype) 1907 molecule_kind => molecule_kinds(ibox)%list%els(imolecule) 1908 CALL get_molecule_kind(molecule_kind, natom=mc_molecule_info%nunits(iunique), & 1909 nbond=nbond, nbend=nbend, ntorsion=ntorsion, atom_list=atom_list) 1910 1911 DO iatom = 1, mc_molecule_info%nunits(iunique) 1912 atomic_kind => atom_list(iatom)%atomic_kind 1913 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1914 name=mc_molecule_info%atom_names(iatom, iunique), & 1915 mass=mc_molecule_info%mass(iatom, iunique)) 1916 ENDDO 1917 1918! compute the probabilities of doing any particular kind of conformation change 1919 total = nbond + nbend + ntorsion 1920 1921 IF (total == 0) THEN 1922 mc_molecule_info%conf_prob(:, iunique) = 0.0e0_dp 1923 ELSE 1924 mc_molecule_info%conf_prob(1, iunique) = REAL(nbond, dp)/REAL(total, dp) 1925 mc_molecule_info%conf_prob(2, iunique) = REAL(nbend, dp)/REAL(total, dp) 1926 mc_molecule_info%conf_prob(3, iunique) = REAL(ntorsion, dp)/REAL(total, dp) 1927 ENDIF 1928 1929 ENDDO 1930 ENDDO 1931 1932! figure out the empty coordinates 1933 IF (PRESENT(coordinates_empty)) THEN 1934 ALLOCATE (coordinates_empty(1:3, 1:mc_molecule_info%nunits(1))) 1935 DO iunit = 1, mc_molecule_info%nunits(1) 1936 coordinates_empty(1:3, iunit) = particles(1)%list%els(iunit)%r(1:3) 1937 ENDDO 1938 ENDIF 1939 1940! now we need to figure out the total number of molecules of each kind in 1941! each box, and the total number of interaction sites in each box 1942 mc_molecule_info%nchains(:, :) = 0 1943 DO iunique = 1, nmol_types 1944 DO ibox = 1, nboxes 1945 IF (ibox == skip_box) CYCLE 1946 DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:)) 1947 molecule_kind => molecule_kinds(ibox)%list%els(imolecule) 1948 CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, & 1949 name=name) 1950 IF (TRIM(name) .NE. mc_molecule_info%names(iunique)) CYCLE 1951 mc_molecule_info%nchains(iunique, ibox) = mc_molecule_info%nchains(iunique, ibox) + nmolecule 1952 ENDDO 1953 ENDDO 1954 ENDDO 1955 mc_molecule_info%nchain_total = 0 1956 DO ibox = 1, nboxes 1957 mc_molecule_info%nunits_tot(ibox) = 0 1958 IF (ibox == skip_box) CYCLE 1959 DO iunique = 1, nmol_types 1960 mc_molecule_info%nunits_tot(ibox) = mc_molecule_info%nunits_tot(ibox) + & 1961 mc_molecule_info%nunits(iunique)*mc_molecule_info%nchains(iunique, ibox) 1962 ENDDO 1963 mc_molecule_info%nchain_total = mc_molecule_info%nchain_total + SUM(mc_molecule_info%nchains(:, ibox)) 1964 ENDDO 1965 1966! now we need to figure out which type every molecule is, 1967! and which force_env (box) it's in 1968 ALLOCATE (mc_molecule_info%mol_type(1:mc_molecule_info%nchain_total)) 1969 ALLOCATE (mc_molecule_info%in_box(1:mc_molecule_info%nchain_total)) 1970 1971 last_mol = 0 1972 DO ibox = 1, nboxes 1973 IF (ibox == skip_box) CYCLE 1974 first_mol = last_mol + 1 1975 last_mol = first_mol + SUM(mc_molecule_info%nchains(:, ibox)) - 1 1976 mc_molecule_info%in_box(first_mol:last_mol) = ibox 1977 DO imolecule = 1, SIZE(molecule_kinds(ibox)%list%els(:)) 1978 molecule_kind => molecule_kinds(ibox)%list%els(imolecule) 1979 CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, & 1980 name=name, molecule_list=molecule_list) 1981! figure out which molecule number this is 1982 this_molecule = 0 1983 DO iunique = 1, nmol_types 1984 IF (TRIM(name) .EQ. TRIM(mc_molecule_info%names(iunique))) THEN 1985 this_molecule = iunique 1986 ENDIF 1987 ENDDO 1988 DO imol = 1, SIZE(molecule_list(:)) 1989 mc_molecule_info%mol_type(first_mol + molecule_list(imol) - 1) = this_molecule 1990 ENDDO 1991 ENDDO 1992 ENDDO 1993 1994! get rid of stuff we don't need 1995 DEALLOCATE (names_init) 1996 DEALLOCATE (molecule_kinds) 1997 DEALLOCATE (particles) 1998 1999 END SUBROUTINE MC_DETERMINE_MOLECULE_INFO 2000 2001! ************************************************************************************************** 2002!> \brief deallocates all the arrays in the mc_molecule_info_type 2003!> \param mc_molecule_info the structure we wish to deallocate 2004!> 2005!> Suitable for parallel. 2006!> \author MJM 2007! ************************************************************************************************** 2008 SUBROUTINE mc_molecule_info_destroy(mc_molecule_info) 2009 2010 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 2011 2012 CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_molecule_info_destroy', & 2013 routineP = moduleN//':'//routineN 2014 2015 DEALLOCATE (mc_molecule_info%nchains) 2016 DEALLOCATE (mc_molecule_info%nunits) 2017 DEALLOCATE (mc_molecule_info%mol_type) 2018 DEALLOCATE (mc_molecule_info%in_box) 2019 DEALLOCATE (mc_molecule_info%names) 2020 DEALLOCATE (mc_molecule_info%atom_names) 2021 DEALLOCATE (mc_molecule_info%conf_prob) 2022 DEALLOCATE (mc_molecule_info%nunits_tot) 2023 DEALLOCATE (mc_molecule_info%mass) 2024 2025 DEALLOCATE (mc_molecule_info) 2026 NULLIFY (mc_molecule_info) 2027 2028 END SUBROUTINE mc_molecule_info_destroy 2029 2030! ************************************************************************************************** 2031!> \brief ... 2032!> \param mc_par ... 2033! ************************************************************************************************** 2034 SUBROUTINE mc_input_parameters_check(mc_par) 2035 2036! ************************************************************************************************** 2037!> \brief accesses the private elements of the mc_molecule_info_type 2038!> \param mc_molecule_info the structure you want the parameters for 2039!> \param nmol_types the number of molecule types in the simulation 2040!> 2041!> Suitable for parallel. 2042!> \author MJM 2043 TYPE(mc_simpar_type), POINTER :: mc_par 2044 2045 CHARACTER(len=*), PARAMETER :: routineN = 'mc_input_parameters_check', & 2046 routineP = moduleN//':'//routineN 2047 2048 INTEGER :: imol_type, nboxes, nchain_total, & 2049 nmol_types 2050 INTEGER, DIMENSION(:), POINTER :: nunits 2051 LOGICAL :: lhmc 2052 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 2053 2054 CALL get_mc_par(mc_par, mc_molecule_info=mc_molecule_info, lhmc=lhmc) 2055 2056 CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types, & 2057 nboxes=nboxes, nunits=nunits, nchain_total=nchain_total) 2058 2059! if we're only doing HMC, we can skip these checks 2060 IF (lhmc) RETURN 2061 2062! the final value of these arrays needs to be 1.0, otherwise there is 2063! a chance we won't select a molecule type for a move 2064 IF (mc_par%pmavbmc_mol(nmol_types) .LT. 0.9999E0_dp) THEN 2065 CPABORT('The last value of PMAVBMC_MOL needs to be 1.0') 2066 ENDIF 2067 IF (mc_par%pmswap_mol(nmol_types) .LT. 0.9999E0_dp) THEN 2068 CPABORT('The last value of PMSWAP_MOL needs to be 1.0') 2069 ENDIF 2070 IF (mc_par%pmtraion_mol(nmol_types) .LT. 0.9999E0_dp) THEN 2071 CPABORT('The last value of PMTRAION_MOL needs to be 1.0') 2072 ENDIF 2073 IF (mc_par%pmtrans_mol(nmol_types) .LT. 0.9999E0_dp) THEN 2074 CPABORT('The last value of PMTRANS_MOL needs to be 1.0') 2075 ENDIF 2076 IF (mc_par%pmrot_mol(nmol_types) .LT. 0.9999E0_dp) THEN 2077 CPABORT('The last value of PMROT_MOL needs to be 1.0') 2078 ENDIF 2079 2080! check to make sure we have all the desired values for various 2081 2082 ! some ensembles require multiple boxes or molecule types 2083 SELECT CASE (mc_par%ensemble) 2084 CASE ("GEMC_NPT") 2085 IF (nmol_types .LE. 1) & 2086 CPABORT('Cannot have GEMC-NPT simulation with only one molecule type') 2087 IF (nboxes .LE. 1) & 2088 CPABORT('Cannot have GEMC-NPT simulation with only one box') 2089 CASE ("GEMC_NVT") 2090 IF (nboxes .LE. 1) & 2091 CPABORT('Cannot have GEMC-NVT simulation with only one box') 2092 CASE ("TRADITIONAL") 2093 IF (mc_par%pmswap .GT. 0.0E0_dp) & 2094 CPABORT('You cannot do swap moves in a system with only one box') 2095 CASE ("VIRIAL") 2096 IF (nchain_total .NE. 2) & 2097 CPABORT('You need exactly two molecules in the box to compute the second virial.') 2098 END SELECT 2099 2100! can't choose an AVBMC target atom number higher than the number 2101! of units in that molecule 2102 DO imol_type = 1, nmol_types 2103 IF (mc_par%avbmc_atom(imol_type) .GT. nunits(imol_type)) THEN 2104 CPABORT('Cannot have avbmc_atom higher than the number of atoms for that molecule type!') 2105 ENDIF 2106 ENDDO 2107 2108 END SUBROUTINE mc_input_parameters_check 2109 2110END MODULE mc_types 2111 2112