1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief contains miscellaneous subroutines used in the Monte Carlo runs,mostly 8!> geared towards changes in coordinates 9!> \author MJM 10! ************************************************************************************************** 11MODULE mc_coordinates 12 USE cell_types, ONLY: cell_type,& 13 get_cell 14 USE cp_subsys_types, ONLY: cp_subsys_get,& 15 cp_subsys_type 16 USE force_env_methods, ONLY: force_env_calc_energy_force 17 USE force_env_types, ONLY: force_env_get,& 18 force_env_type 19 USE kinds, ONLY: dp 20 USE mathconstants, ONLY: pi 21 USE mc_types, ONLY: get_mc_molecule_info,& 22 get_mc_par,& 23 mc_molecule_info_type,& 24 mc_simpar_type 25 USE message_passing, ONLY: mp_bcast 26 USE molecule_types, ONLY: molecule_type 27 USE parallel_rng_types, ONLY: rng_stream_type 28 USE particle_list_types, ONLY: particle_list_type 29 USE physcon, ONLY: angstrom 30#include "../../base/base_uses.f90" 31 32 IMPLICIT NONE 33 34 PRIVATE 35 36 PRIVATE :: generate_avbmc_insertion 37 38 PUBLIC :: generate_cbmc_swap_config, & 39 get_center_of_mass, mc_coordinate_fold, & 40 find_mc_test_molecule, & 41 create_discrete_array, & 42 check_for_overlap, & 43 rotate_molecule, & 44 cluster_search 45 46 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_coordinates' 47 48CONTAINS 49 50! ************************************************************************************************** 51!> \brief looks for overlaps (intermolecular distances less than rmin) 52!> \param force_env the force environment containing the coordinates 53!> \param nchains the number of molecules of each type in the box 54!> \param nunits the number of interaction sites for each molecule 55!> \param loverlap returns .TRUE. if atoms overlap 56!> \param mol_type an array that indicates the type of each molecule 57!> \param cell_length the length of the box...if none is specified, 58!> it uses the cell found in the force_env 59!> \param molecule_number if present, just look for overlaps with this 60!> molecule 61!> 62!> Suitable for parallel use. 63!> \author MJM 64! ************************************************************************************************** 65 SUBROUTINE check_for_overlap(force_env, nchains, nunits, loverlap, mol_type, & 66 cell_length, molecule_number) 67 68 TYPE(force_env_type), POINTER :: force_env 69 INTEGER, DIMENSION(:), INTENT(IN) :: nchains, nunits 70 LOGICAL, INTENT(OUT) :: loverlap 71 INTEGER, DIMENSION(:), INTENT(IN) :: mol_type 72 REAL(KIND=dp), DIMENSION(1:3), INTENT(IN), & 73 OPTIONAL :: cell_length 74 INTEGER, INTENT(IN), OPTIONAL :: molecule_number 75 76 CHARACTER(len=*), PARAMETER :: routineN = 'check_for_overlap', & 77 routineP = moduleN//':'//routineN 78 79 INTEGER :: handle, imol, iunit, jmol, jstart, & 80 junit, nend, nstart, nunit, nunits_i, & 81 nunits_j 82 LOGICAL :: lall 83 REAL(KIND=dp) :: dist, rmin 84 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r 85 REAL(KIND=dp), DIMENSION(1:3) :: abc, box_length, RIJ 86 TYPE(cell_type), POINTER :: cell 87 TYPE(cp_subsys_type), POINTER :: oldsys 88 TYPE(particle_list_type), POINTER :: particles 89 90! begin the timing of the subroutine 91 92 CALL timeset(routineN, handle) 93 94 NULLIFY (oldsys, particles) 95 96! initialize some stuff 97 loverlap = .FALSE. 98 rmin = 3.57106767_dp ! 1 angstrom squared 99 100! get the particle coordinates and the cell length 101 CALL force_env_get(force_env, cell=cell, subsys=oldsys) 102 CALL get_cell(cell, abc=abc) 103 CALL cp_subsys_get(oldsys, particles=particles) 104 105 ALLOCATE (r(1:3, 1:MAXVAL(nunits), 1:SUM(nchains))) 106 107 IF (PRESENT(cell_length)) THEN 108 box_length(1:3) = cell_length(1:3) 109 ELSE 110 box_length(1:3) = abc(1:3) 111 ENDIF 112 113! put the coordinates into an easier matrix to manipulate 114 junit = 0 115 DO imol = 1, SUM(nchains) 116 nunit = nunits(mol_type(imol)) 117 DO iunit = 1, nunit 118 junit = junit + 1 119 r(1:3, iunit, imol) = particles%els(junit)%r(1:3) 120 ENDDO 121 ENDDO 122 123! now let's find the LJ energy between all the oxygens and 124! the charge interactions 125 lall = .TRUE. 126 jstart = 1 127 IF (PRESENT(molecule_number)) THEN 128 lall = .FALSE. 129 nstart = molecule_number 130 nend = molecule_number 131 ELSE 132 nstart = 1 133 nend = SUM(nchains(:)) 134 ENDIF 135 DO imol = nstart, nend 136 IF (lall) jstart = imol + 1 137 nunits_i = nunits(mol_type(imol)) 138 DO jmol = jstart, SUM(nchains(:)) 139 IF (imol == jmol) CYCLE 140 nunits_j = nunits(mol_type(jmol)) 141 142 DO iunit = 1, nunits_i 143 DO junit = 1, nunits_j 144! find the minimum image distance 145 RIJ(1) = r(1, iunit, imol) - r(1, junit, jmol) - & 146 box_length(1)*ANINT( & 147 (r(1, iunit, imol) - r(1, junit, jmol))/box_length(1)) 148 RIJ(2) = r(2, iunit, imol) - r(2, junit, jmol) - & 149 box_length(2)*ANINT( & 150 (r(2, iunit, imol) - r(2, junit, jmol))/box_length(2)) 151 RIJ(3) = r(3, iunit, imol) - r(3, junit, jmol) - & 152 box_length(3)*ANINT( & 153 (r(3, iunit, imol) - r(3, junit, jmol))/box_length(3)) 154 155 dist = RIJ(1)**2 + RIJ(2)**2 + RIJ(3)**2 156 157 IF (dist < rmin) THEN 158 loverlap = .TRUE. 159 DEALLOCATE (r) 160 161 CALL timestop(handle) 162 RETURN 163 ENDIF 164 165 ENDDO 166 ENDDO 167 ENDDO 168 ENDDO 169 170 DEALLOCATE (r) 171 172! end the timing 173 CALL timestop(handle) 174 175 END SUBROUTINE check_for_overlap 176 177! ************************************************************************************************** 178!> \brief calculates the center of mass of a given molecule 179!> \param coordinates the coordinates of the atoms in the molecule 180!> \param natom the number of atoms in the molecule 181!> \param center_of_mass the coordinates of the center of mass 182!> \param mass the mass of the atoms in the molecule 183!> 184!> Designed for parallel use. 185!> \author MJM 186! ************************************************************************************************** 187 SUBROUTINE get_center_of_mass(coordinates, natom, center_of_mass, & 188 mass) 189 190 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: coordinates 191 INTEGER, INTENT(IN) :: natom 192 REAL(KIND=dp), DIMENSION(1:3), INTENT(OUT) :: center_of_mass 193 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: mass 194 195 CHARACTER(len=*), PARAMETER :: routineN = 'get_center_of_mass', & 196 routineP = moduleN//':'//routineN 197 198 INTEGER :: handle, i, iatom 199 REAL(KIND=dp) :: total_mass 200 201! begin the timing of the subroutine 202 203 CALL timeset(routineN, handle) 204 205 total_mass = SUM(mass(1:natom)) 206 center_of_mass(:) = 0.0E0_dp 207 208 DO iatom = 1, natom 209 DO i = 1, 3 210 center_of_mass(i) = center_of_mass(i) + & 211 mass(iatom)*coordinates(i, iatom) 212 ENDDO 213 ENDDO 214 215 center_of_mass(1:3) = center_of_mass(1:3)/total_mass 216 217! end the timing 218 CALL timestop(handle) 219 220 END SUBROUTINE get_center_of_mass 221 222! ************************************************************************************************** 223!> \brief folds all the coordinates into the center simulation box using 224!> a center of mass cutoff 225!> \param coordinates the coordinates of the atoms in the system 226!> \param nchains_tot the total number of molecules in the box 227!> \param mol_type an array that indicates the type of every molecule in the box 228!> \param mass the mass of every atom for all molecule types 229!> \param nunits the number of interaction sites for each molecule type 230!> \param box_length an array for the lengths of the simulation box sides 231!> 232!> Designed for parallel use. 233!> \author MJM 234! ************************************************************************************************** 235 SUBROUTINE mc_coordinate_fold(coordinates, nchains_tot, mol_type, mass, nunits, box_length) 236 237 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: coordinates 238 INTEGER, INTENT(IN) :: nchains_tot 239 INTEGER, DIMENSION(:), INTENT(IN) :: mol_type 240 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: mass 241 INTEGER, DIMENSION(:), INTENT(IN) :: nunits 242 REAL(KIND=dp), DIMENSION(1:3), INTENT(IN) :: box_length 243 244 CHARACTER(len=*), PARAMETER :: routineN = 'mc_coordinate_fold', & 245 routineP = moduleN//':'//routineN 246 247 INTEGER :: end_atom, handle, iatom, imolecule, & 248 jatom, molecule_type, natoms, & 249 start_atom 250 REAL(KIND=dp), DIMENSION(1:3) :: center_of_mass 251 252! begin the timing of the subroutine 253 254 CALL timeset(routineN, handle) 255 256! loop over all molecules 257 end_atom = 0 258 DO imolecule = 1, nchains_tot 259 molecule_type = mol_type(imolecule) 260 natoms = nunits(molecule_type) 261 start_atom = end_atom + 1 262 end_atom = start_atom + natoms - 1 263 CALL get_center_of_mass(coordinates(:, start_atom:end_atom), & 264 natoms, center_of_mass(:), mass(:, molecule_type)) 265 DO iatom = 1, natoms 266 jatom = iatom + start_atom - 1 267 coordinates(1, jatom) = coordinates(1, jatom) - & 268 box_length(1)*FLOOR(center_of_mass(1)/box_length(1)) 269 coordinates(2, jatom) = coordinates(2, jatom) - & 270 box_length(2)*FLOOR(center_of_mass(2)/box_length(2)) 271 coordinates(3, jatom) = coordinates(3, jatom) - & 272 box_length(3)*FLOOR(center_of_mass(3)/box_length(2)) 273 ENDDO 274 275 ENDDO 276 277! end the timing 278 CALL timestop(handle) 279 280 END SUBROUTINE mc_coordinate_fold 281 282! ************************************************************************************************** 283!> \brief takes the last molecule in a force environment and moves it around 284!> to different center of mass positions and orientations, selecting one 285!> based on the rosenbluth weight 286!> \param force_env the force environment containing the coordinates 287!> \param BETA the value of 1/kT for this simulations, in a.u. 288!> \param max_val ... 289!> \param min_val ... 290!> \param exp_max_val ... 291!> \param exp_min_val ... 292!> \param nswapmoves the number of desired trial configurations 293!> \param rosenbluth_weight the Rosenbluth weight for this set of configs 294!> \param start_atom the atom number that the molecule to be swapped starts on 295!> \param natoms_tot the total number of interaction sites in the box 296!> \param nunits the number of interaction sites for every molecule_type 297!> \param nunits_mol ... 298!> \param mass the mass for every interaction site of every molecule type 299!> \param loverlap the flag that indicates if all of the configs have an 300!> atomic overlap 301!> \param choosen_energy the energy of the chosen config 302!> \param old_energy the energy that we subtract from all of the trial 303!> energies to prevent numerical overflows 304!> \param ionode indicates if we're on the main CPU 305!> \param lremove is this the Rosenbluth weight for a removal box? 306!> \param mol_type an array that contains the molecule type for every atom in the box 307!> \param nchains the number of molecules of each type in this box 308!> \param source the MPI source value, for broadcasts 309!> \param group the MPI group value, for broadcasts 310!> \param rng_stream the random number stream that we draw from 311!> \param avbmc_atom ... 312!> \param rmin ... 313!> \param rmax ... 314!> \param move_type ... 315!> \param target_atom ... 316!> \par Optional Avbmc Flags 317!> - avbmc_atom: the atom number that serves for the target atom in each 318!> molecule (1 is the first atom in the molecule, etc.) 319!> - rmin: the minimum AVBMC radius for the shell around the target 320!> - rmax: the maximum AVBMC radius for the shell around the target 321!> - move_type: generate configs in the "in" or "out" volume 322!> - target_atom: the number of the avbmc atom in the target molecule 323!> \par 324!> Suitable for parallel. 325!> \author MJM 326! ************************************************************************************************** 327 SUBROUTINE generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, & 328 exp_min_val, nswapmoves, rosenbluth_weight, start_atom, natoms_tot, nunits, nunits_mol, & 329 mass, loverlap, choosen_energy, old_energy, ionode, lremove, mol_type, nchains, source, & 330 group, rng_stream, avbmc_atom, rmin, rmax, move_type, target_atom) 331 332 TYPE(force_env_type), POINTER :: force_env 333 REAL(KIND=dp), INTENT(IN) :: BETA, max_val, min_val, exp_max_val, & 334 exp_min_val 335 INTEGER, INTENT(IN) :: nswapmoves 336 REAL(KIND=dp), INTENT(OUT) :: rosenbluth_weight 337 INTEGER, INTENT(IN) :: start_atom, natoms_tot 338 INTEGER, DIMENSION(:), INTENT(IN) :: nunits 339 INTEGER, INTENT(IN) :: nunits_mol 340 REAL(dp), DIMENSION(1:nunits_mol), INTENT(IN) :: mass 341 LOGICAL, INTENT(OUT) :: loverlap 342 REAL(KIND=dp), INTENT(OUT) :: choosen_energy 343 REAL(KIND=dp), INTENT(IN) :: old_energy 344 LOGICAL, INTENT(IN) :: ionode, lremove 345 INTEGER, DIMENSION(:), INTENT(IN) :: mol_type, nchains 346 INTEGER, INTENT(IN) :: source, group 347 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 348 INTEGER, INTENT(IN), OPTIONAL :: avbmc_atom 349 REAL(KIND=dp), INTENT(IN), OPTIONAL :: rmin, rmax 350 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: move_type 351 INTEGER, INTENT(IN), OPTIONAL :: target_atom 352 353 CHARACTER(len=*), PARAMETER :: routineN = 'generate_cbmc_swap_config', & 354 routineP = moduleN//':'//routineN 355 356 INTEGER :: atom_number, choosen, end_atom, handle, & 357 i, iatom, imolecule, imove, & 358 molecule_number 359 LOGICAL :: all_overlaps 360 LOGICAL, ALLOCATABLE, DIMENSION(:) :: loverlap_array 361 REAL(KIND=dp) :: bias_energy, exponent, rand, & 362 total_running_weight 363 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: boltz_weights, trial_energy 364 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r_old 365 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r 366 REAL(KIND=dp), DIMENSION(1:3) :: abc, center_of_mass, diff, r_insert 367 TYPE(cell_type), POINTER :: cell 368 TYPE(cp_subsys_type), POINTER :: oldsys 369 TYPE(particle_list_type), POINTER :: particles 370 371! begin the timing of the subroutine 372 373 CALL timeset(routineN, handle) 374 375 NULLIFY (oldsys) 376! get the particle coordinates and the cell length 377 CALL force_env_get(force_env, cell=cell, subsys=oldsys) 378 CALL get_cell(cell, abc=abc) 379 CALL cp_subsys_get(oldsys, particles=particles) 380 381! do some checking to make sure we have all the data we need 382 IF (PRESENT(avbmc_atom)) THEN 383 IF (.NOT. PRESENT(rmin) .OR. .NOT. PRESENT(rmax) .OR. & 384 .NOT. PRESENT(move_type) .OR. .NOT. PRESENT(target_atom)) THEN 385 CPABORT('AVBMC swap move is missing information!') 386 ENDIF 387 ENDIF 388 389 ALLOCATE (r_old(1:3, 1:natoms_tot)) 390 ALLOCATE (r(1:3, 1:natoms_tot, 1:nswapmoves)) 391 ALLOCATE (trial_energy(1:nswapmoves)) 392 ALLOCATE (boltz_weights(1:nswapmoves)) 393 ALLOCATE (loverlap_array(1:nswapmoves)) 394 395! initialize the arrays that need it 396 loverlap_array(:) = .FALSE. 397 loverlap = .FALSE. 398 boltz_weights(:) = 0.0_dp 399 trial_energy(:) = 0.0_dp 400 r(:, :, :) = 0.0_dp 401 choosen_energy = 0.0_dp 402 rosenbluth_weight = 0.0_dp 403 404! save the positions of the molecules 405 DO imove = 1, nswapmoves 406 DO iatom = 1, natoms_tot 407 r(1:3, iatom, imove) = particles%els(iatom)%r(1:3) 408 ENDDO 409 ENDDO 410 411! save the remove coordinates 412 DO iatom = 1, natoms_tot 413 r_old(1:3, iatom) = r(1:3, iatom, 1) 414 ENDDO 415 416! figure out the numbers of the first and last atoms in the molecule 417 end_atom = start_atom + nunits_mol - 1 418! figure out which molecule number we're on 419 molecule_number = 0 420 atom_number = 1 421 DO imolecule = 1, SUM(nchains(:)) 422 IF (atom_number == start_atom) THEN 423 molecule_number = imolecule 424 EXIT 425 ENDIF 426 atom_number = atom_number + nunits(mol_type(imolecule)) 427 ENDDO 428 IF (molecule_number == 0) CALL cp_abort(__LOCATION__, & 429 'CBMC swap move cannot find which molecule number it needs') 430 431 IF (lremove) THEN 432 CALL check_for_overlap(force_env, nchains, nunits, loverlap_array(1), & 433 mol_type) 434 CALL mp_bcast(loverlap_array(1), source, group) 435 436 IF (loverlap_array(1)) THEN 437 IF (ionode) THEN 438 WRITE (*, *) start_atom, end_atom, natoms_tot 439 DO iatom = 1, natoms_tot 440 WRITE (*, *) r(1:3, iatom, 1) 441 ENDDO 442 ENDIF 443 CPABORT('CBMC swap move found an overlap in the old config') 444 ENDIF 445 ENDIF 446 447 DO imove = 1, nswapmoves 448 449! drop into serial 450 IF (ionode) THEN 451 452 IF (PRESENT(avbmc_atom)) THEN 453! find an AVBMC insertion point 454 CALL generate_avbmc_insertion(rmin, rmax, & 455 r_old(1:3, target_atom), & 456 move_type, r_insert(:), abc(:), rng_stream) 457 458 DO i = 1, 3 459 diff(i) = r_insert(i) - r_old(i, start_atom + avbmc_atom - 1) 460 ENDDO 461 462 ELSE 463! find a new insertion point somewhere in the box 464 DO i = 1, 3 465 rand = rng_stream%next() 466 r_insert(i) = rand*abc(i) 467 ENDDO 468 469! find the center of mass of the insertion molecule 470 CALL get_center_of_mass(r(:, start_atom:end_atom, imove), nunits_mol, & 471 center_of_mass(:), mass(:)) 472 473! move the molecule to the insertion point 474 475 DO i = 1, 3 476 diff(i) = r_insert(i) - center_of_mass(i) 477 ENDDO 478 479 ENDIF 480 481 DO iatom = start_atom, end_atom 482 r(1:3, iatom, imove) = r(1:3, iatom, imove) + diff(1:3) 483 ENDDO 484 485! rotate the molecule...this routine is only made for serial use 486 CALL rotate_molecule(r(:, start_atom:end_atom, imove), mass(:), & 487 nunits_mol, rng_stream) 488 489 IF (imove == 1 .AND. lremove) THEN 490 DO iatom = 1, natoms_tot 491 r(1:3, iatom, 1) = r_old(1:3, iatom) 492 ENDDO 493 ENDIF 494 495 ENDIF 496 497 CALL mp_bcast(r(:, :, imove), source, group) 498 499! calculate the energy and boltzman weight of the config 500 DO iatom = start_atom, end_atom 501 particles%els(iatom)%r(1:3) = r(1:3, iatom, imove) 502 ENDDO 503 504 CALL check_for_overlap(force_env, nchains, nunits, loverlap_array(imove), & 505 mol_type, molecule_number=molecule_number) 506 IF (loverlap_array(imove)) THEN 507 boltz_weights(imove) = 0.0_dp 508 CYCLE 509 ENDIF 510 511 CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.) 512 CALL force_env_get(force_env, & 513 potential_energy=bias_energy) 514 515 trial_energy(imove) = (bias_energy - old_energy) 516 exponent = -BETA*trial_energy(imove) 517 518 IF (exponent .GT. exp_max_val) THEN 519 boltz_weights(imove) = max_val 520 ELSEIF (exponent .LT. exp_min_val) THEN 521 boltz_weights(imove) = min_val 522 ELSE 523 boltz_weights(imove) = EXP(exponent) 524 ENDIF 525 526 ENDDO 527 528! now we need to pick a configuration based on the Rosenbluth weight, 529! which is just the sum of the Boltzmann weights 530 rosenbluth_weight = SUM(boltz_weights(:)) 531 IF (rosenbluth_weight .EQ. 0.0_dp .AND. lremove) THEN 532! should never have 0.0 for an old weight...causes a divide by zero 533! in the acceptance rule 534 IF (ionode) THEN 535 WRITE (*, *) boltz_weights(1:nswapmoves) 536 WRITE (*, *) start_atom, end_atom, lremove 537 WRITE (*, *) loverlap_array(1:nswapmoves) 538 WRITE (*, *) natoms_tot 539 WRITE (*, *) 540 DO iatom = 1, natoms_tot 541 WRITE (*, *) r(1:3, iatom, 1)*angstrom 542 ENDDO 543 ENDIF 544 CPABORT('CBMC swap move found a bad old weight') 545 ENDIF 546 all_overlaps = .TRUE. 547 total_running_weight = 0.0E0_dp 548 choosen = 0 549 IF (ionode) THEN 550 rand = rng_stream%next() 551! CALL random_number(rand) 552 ENDIF 553 CALL mp_bcast(rand, source, group) 554 DO imove = 1, nswapmoves 555 IF (loverlap_array(imove)) CYCLE 556 all_overlaps = .FALSE. 557 total_running_weight = total_running_weight + boltz_weights(imove) 558 IF (total_running_weight .GE. rand*rosenbluth_weight) THEN 559 choosen = imove 560 EXIT 561 ENDIF 562 ENDDO 563 564 IF (all_overlaps) THEN 565 loverlap = .TRUE. 566 567! if this is an old configuration, we always choose the first one... 568! this should never be the case, but I'm testing something 569 IF (lremove) THEN 570 IF (ionode) THEN 571 WRITE (*, *) boltz_weights(1:nswapmoves) 572 WRITE (*, *) start_atom, end_atom, lremove 573 WRITE (*, *) loverlap_array(1:nswapmoves) 574 DO iatom = 1, natoms_tot 575 WRITE (*, *) r(1:3, iatom, 1) 576 ENDDO 577 ENDIF 578 CPABORT('CBMC swap move found all overlaps for the remove config') 579 ENDIF 580 581 DEALLOCATE (r_old) 582 DEALLOCATE (r) 583 DEALLOCATE (trial_energy) 584 DEALLOCATE (boltz_weights) 585 DEALLOCATE (loverlap_array) 586 CALL timestop(handle) 587 RETURN 588 ENDIF 589 590! make sure a configuration was chosen 591 IF (choosen == 0) & 592 CPABORT('CBMC swap move failed to select config') 593 594! if this is an old configuration, we always choose the first one 595 IF (lremove) choosen = 1 596 597! set the energy for the configuration 598 choosen_energy = trial_energy(choosen) 599 600! copy the coordinates to the force environment 601 DO iatom = 1, natoms_tot 602 particles%els(iatom)%r(1:3) = r(1:3, iatom, choosen) 603 ENDDO 604 605 DEALLOCATE (r_old) 606 DEALLOCATE (r) 607 DEALLOCATE (trial_energy) 608 DEALLOCATE (boltz_weights) 609 DEALLOCATE (loverlap_array) 610 611! end the timing 612 CALL timestop(handle) 613 614 END SUBROUTINE generate_cbmc_swap_config 615 616! ************************************************************************************************** 617!> \brief rotates a molecule randomly around the center of mass, 618!> sequentially in x, y, and z directions 619!> \param r the coordinates of the molecule to rotate 620!> \param mass the mass of all the atoms in the molecule 621!> \param natoms the number of atoms in the molecule 622!> \param rng_stream the stream we pull random numbers from 623!> 624!> Use only in serial. 625!> \author MJM 626! ************************************************************************************************** 627 SUBROUTINE rotate_molecule(r, mass, natoms, rng_stream) 628 629 INTEGER, INTENT(IN) :: natoms 630 REAL(KIND=dp), DIMENSION(1:natoms), INTENT(IN) :: mass 631 REAL(KIND=dp), DIMENSION(1:3, 1:natoms), & 632 INTENT(INOUT) :: r 633 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 634 635 CHARACTER(len=*), PARAMETER :: routineN = 'rotate_molecule', & 636 routineP = moduleN//':'//routineN 637 638 INTEGER :: handle, iunit 639 REAL(KIND=dp) :: cosdg, dgamma, rand, rx, rxnew, ry, & 640 rynew, rz, rznew, sindg 641 REAL(KIND=dp), DIMENSION(1:3) :: center_of_mass 642 643! begin the timing of the subroutine 644 645 CALL timeset(routineN, handle) 646 647! find the center of mass of the molecule 648 CALL get_center_of_mass(r(:, :), natoms, center_of_mass(:), mass(:)) 649 650! call a random number to figure out how far we're moving 651 rand = rng_stream%next() 652 dgamma = pi*(rand - 0.5E0_dp)*2.0E0_dp 653 654! *** set up the rotation matrix *** 655 656 cosdg = COS(dgamma) 657 sindg = SIN(dgamma) 658 659! *** ROTATE UNITS OF I AROUND X-AXIS *** 660 661 DO iunit = 1, natoms 662 ry = r(2, iunit) - center_of_mass(2) 663 rz = r(3, iunit) - center_of_mass(3) 664 rynew = cosdg*ry + sindg*rz 665 rznew = cosdg*rz - sindg*ry 666 667 r(2, iunit) = rynew + center_of_mass(2) 668 r(3, iunit) = rznew + center_of_mass(3) 669 670 ENDDO 671 672! *** ROTATE UNITS OF I AROUND y-AXIS *** 673 674 DO iunit = 1, natoms 675 rx = r(1, iunit) - center_of_mass(1) 676 rz = r(3, iunit) - center_of_mass(3) 677 rxnew = cosdg*rx + sindg*rz 678 rznew = cosdg*rz - sindg*rx 679 680 r(1, iunit) = rxnew + center_of_mass(1) 681 r(3, iunit) = rznew + center_of_mass(3) 682 683 ENDDO 684 685! *** ROTATE UNITS OF I AROUND z-AXIS *** 686 687 DO iunit = 1, natoms 688 rx = r(1, iunit) - center_of_mass(1) 689 ry = r(2, iunit) - center_of_mass(2) 690 rxnew = cosdg*rx + sindg*ry 691 rynew = cosdg*ry - sindg*rx 692 693 r(1, iunit) = rxnew + center_of_mass(1) 694 r(2, iunit) = rynew + center_of_mass(2) 695 696 ENDDO 697 698! end the timing 699 CALL timestop(handle) 700 701 END SUBROUTINE rotate_molecule 702 703! ************************************************************************************************** 704!> \brief selects a molecule at random to perform a MC move on...you can specify 705!> the box the molecule should be in, its type, both, or neither 706!> \param mc_molecule_info the structure that contains some global molecule data 707!> \param start_atom the number of the first atom in the chosen molecule in relation 708!> to the force_env it's in 709!> \param box_number the box the chosen molecule is in 710!> \param molecule_type the type of molecule the chosen molecule is 711!> \param rng_stream the stream we pull random numbers from 712!> \param box if present, tells the routine which box to grab a molecule from 713!> \param molecule_type_old if present, tells the routine which molecule type to select from 714!> \author MJM 715! ************************************************************************************************** 716 SUBROUTINE find_mc_test_molecule(mc_molecule_info, start_atom, & 717 box_number, molecule_type, rng_stream, box, molecule_type_old) 718 719 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 720 INTEGER, INTENT(OUT) :: start_atom, box_number, molecule_type 721 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 722 INTEGER, INTENT(IN), OPTIONAL :: box, molecule_type_old 723 724 CHARACTER(LEN=*), PARAMETER :: routineN = 'find_mc_test_molecule', & 725 routineP = moduleN//':'//routineN 726 727 INTEGER :: handle, ibox, imol_type, imolecule, & 728 jbox, molecule_number, nchains_tot, & 729 start_mol 730 INTEGER, DIMENSION(:), POINTER :: mol_type, nunits 731 INTEGER, DIMENSION(:, :), POINTER :: nchains 732 REAL(KIND=dp) :: rand 733 734! begin the timing of the subroutine 735 736 CALL timeset(routineN, handle) 737 738 NULLIFY (nunits, mol_type, nchains) 739 CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, nunits=nunits, & 740 mol_type=mol_type) 741 742! initialize the outgoing variables 743 start_atom = 0 744 box_number = 0 745 molecule_type = 0 746 747 IF (PRESENT(box) .AND. PRESENT(molecule_type_old)) THEN 748! only need to find the atom number the molecule starts on 749 rand = rng_stream%next() 750 molecule_number = CEILING(rand*REAL(nchains(molecule_type_old, box), KIND=dp)) 751 752 start_mol = 1 753 DO jbox = 1, box - 1 754 start_mol = start_mol + SUM(nchains(:, jbox)) 755 ENDDO 756 757! adjust to take into account molecules of other types in the box 758 DO imol_type = 1, molecule_type_old - 1 759 molecule_number = molecule_number + nchains(imol_type, box) 760 ENDDO 761 762 start_atom = 1 763 DO imolecule = 1, molecule_number - 1 764 start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1)) 765 ENDDO 766 767 ELSEIF (PRESENT(box)) THEN 768! any molecule in box...need to find molecule type and start atom 769 rand = rng_stream%next() 770 molecule_number = CEILING(rand*REAL(SUM(nchains(:, box)), KIND=dp)) 771 772 start_mol = 1 773 DO jbox = 1, box - 1 774 start_mol = start_mol + SUM(nchains(:, jbox)) 775 ENDDO 776 777 molecule_type = mol_type(start_mol + molecule_number - 1) 778 779! now the starting atom 780 start_atom = 1 781 DO imolecule = 1, molecule_number - 1 782 start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1)) 783 ENDDO 784 785 ELSEIF (PRESENT(molecule_type_old)) THEN 786! any molecule of type molecule_type_old...need to find box number and start atom 787 rand = rng_stream%next() 788 molecule_number = CEILING(rand*REAL(SUM(nchains(molecule_type_old, :)), KIND=dp)) 789 790! find which box it's in 791 nchains_tot = 0 792 DO ibox = 1, SIZE(nchains(molecule_type_old, :)) 793 IF (molecule_number .LE. nchains(molecule_type_old, ibox)) THEN 794 box_number = ibox 795 EXIT 796 ENDIF 797 molecule_number = molecule_number - nchains(molecule_type_old, ibox) 798 ENDDO 799 800 start_mol = 1 801 DO jbox = 1, box_number - 1 802 start_mol = start_mol + SUM(nchains(:, jbox)) 803 ENDDO 804 805! now find the starting atom number 806 DO imol_type = 1, molecule_type_old - 1 807 molecule_number = molecule_number + nchains(imol_type, box_number) 808 ENDDO 809 start_atom = 1 810 DO imolecule = 1, molecule_number - 1 811 start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1)) 812 ENDDO 813 814 ELSE 815! no restrictions...need to find all pieces of data 816 nchains_tot = 0 817 DO ibox = 1, SIZE(nchains(1, :)) 818 nchains_tot = nchains_tot + SUM(nchains(:, ibox)) 819 ENDDO 820 rand = rng_stream%next() 821 molecule_number = CEILING(rand*REAL(nchains_tot, KIND=dp)) 822 823 molecule_type = mol_type(molecule_number) 824 825! now which box it's in 826 DO ibox = 1, SUM(nchains(1, :)) 827 IF (molecule_number .LE. SUM(nchains(:, ibox))) THEN 828 box_number = ibox 829 EXIT 830 ENDIF 831 molecule_number = molecule_number - SUM(nchains(:, ibox)) 832 ENDDO 833 834! now find the starting atom number 835 start_mol = 1 836 DO jbox = 1, box_number - 1 837 start_mol = start_mol + SUM(nchains(:, jbox)) 838 ENDDO 839 start_atom = 1 840 DO imolecule = 1, molecule_number - 1 841 start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1)) 842 ENDDO 843 844 ENDIF 845 846! make sure things are good 847 IF (PRESENT(box)) box_number = box 848 IF (PRESENT(molecule_type_old)) molecule_type = molecule_type_old 849 850 CPASSERT(start_atom /= 0) 851 CPASSERT(box_number /= 0) 852 CPASSERT(molecule_type /= 0) 853 854! end the timing 855 CALL timestop(handle) 856 857 END SUBROUTINE find_mc_test_molecule 858 859! ************************************************************************************************** 860!> \brief generates an array that tells us which sides of the simulation 861!> cell we can increase or decrease using a discrete volume move 862!> \param cell the lengths of the sides of the cell 863!> \param discrete_array the array that indicates which sides we can move 864!> \param step_size the size of the discrete volume move 865!> 866!> Suitable for parallel. 867!> \author MJM 868! ************************************************************************************************** 869 SUBROUTINE create_discrete_array(cell, discrete_array, step_size) 870 871! 1 is for increase, 2 is for decrease 872! 1 is for "yes, we can do the move", 0 is for no 873 874 REAL(dp), DIMENSION(1:3), INTENT(IN) :: cell 875 INTEGER, DIMENSION(1:3, 1:2), INTENT(OUT) :: discrete_array 876 REAL(dp), INTENT(IN) :: step_size 877 878 INTEGER :: iside 879 REAL(dp) :: high_value, length1, length2, low_value 880 881 discrete_array(:, :) = 0 882 883 length1 = ABS(cell(1) - cell(2)) 884 length2 = ABS(cell(2) - cell(3)) 885 886! now let's figure out all the different cases 887 IF (length1 .LT. 0.01_dp*step_size .AND. & 888 length2 .LT. 0.01_dp*step_size) THEN 889! all sides are equal, so we can move up or down 890 discrete_array(1:3, 1) = 1 891 discrete_array(1:3, 2) = 1 892 ELSE 893 894! find the low value and the high value 895 high_value = -1.0_dp 896 low_value = cell(1)*cell(2)*cell(3) 897 DO iside = 1, 3 898 IF (cell(iside) .LT. low_value) low_value = cell(iside) 899 IF (cell(iside) .GT. high_value) high_value = cell(iside) 900 ENDDO 901 DO iside = 1, 3 902! now we see if the value is a high value or a low value...it can only be 903! one of the two 904 IF (ABS(cell(iside) - low_value) .LT. 0.01_dp*step_size) THEN 905! low value, we can only increase the cell size 906 discrete_array(iside, 1) = 1 907 discrete_array(iside, 2) = 0 908 ELSE 909! high value, we can only decrease the cell size 910 discrete_array(iside, 1) = 0 911 discrete_array(iside, 2) = 1 912 ENDIF 913 ENDDO 914 ENDIF 915 916 END SUBROUTINE create_discrete_array 917 918! ************************************************************************************************** 919!> \brief generates an insertion point in either the "in" or the "out" volume 920!> of a target atom, where the "in" volume is a shell with inner radius 921!> rmin and outer radius rmax 922!> \param rmin the minimum AVBMC radius for the shell around the target 923!> \param rmax the maximum AVBMC radius for the shell around the target 924!> \param r_target the coordinates of the target atom 925!> \param move_type generate configs in the "in" or "out" volume 926!> \param r_insert the output insertion site 927!> \param abc the lengths of the sides of the simulation box 928!> \param rng_stream the random number stream that we draw from 929!> 930!> Use only in serial. 931!> \author MJM 932! ************************************************************************************************** 933 SUBROUTINE generate_avbmc_insertion(rmin, rmax, r_target, & 934 move_type, r_insert, abc, rng_stream) 935 936 REAL(KIND=dp), INTENT(IN) :: rmin, rmax 937 REAL(KIND=dp), DIMENSION(1:3), INTENT(IN) :: r_target 938 CHARACTER(LEN=*), INTENT(IN) :: move_type 939 REAL(KIND=dp), DIMENSION(1:3), INTENT(OUT) :: r_insert 940 REAL(KIND=dp), DIMENSION(1:3), INTENT(IN) :: abc 941 TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream 942 943 INTEGER :: i 944 REAL(dp) :: dist, eta_1, eta_2, eta_sq, rand 945 REAL(dp), DIMENSION(1:3) :: RIJ 946 947 r_insert(1:3) = 0.0_dp 948 949 IF (move_type == 'in') THEN 950! generate a random unit vector, from Allen and Tildesley 951 DO 952 eta_1 = rng_stream%next() 953 eta_2 = rng_stream%next() 954 eta_sq = eta_1**2 + eta_2**2 955 IF (eta_sq .LT. 1.0_dp) THEN 956 r_insert(1) = 2.0_dp*eta_1*SQRT(1.0_dp - eta_sq) 957 r_insert(2) = 2.0_dp*eta_2*SQRT(1.0_dp - eta_sq) 958 r_insert(3) = 1.0_dp - 2.0_dp*eta_sq 959 EXIT 960 ENDIF 961 ENDDO 962 963! now scale that vector to be within the "in" region 964 rand = rng_stream%next() 965 r_insert(1:3) = r_insert(1:3)*(rand*(rmax**3 - rmin**3) + rmin**3)** & 966 (1.0_dp/3.0_dp) 967 968 r_insert(1:3) = r_target(1:3) + r_insert(1:3) 969 ELSE 970 971! find a new insertion point somewhere in the box 972 DO 973 DO i = 1, 3 974 rand = rng_stream%next() 975 r_insert(i) = rand*abc(i) 976 ENDDO 977 978! make sure it's not in the "in" region 979 RIJ(1) = r_insert(1) - r_target(1) - abc(1)* & 980 ANINT((r_insert(1) - r_target(1))/abc(1)) 981 RIJ(2) = r_insert(2) - r_target(2) - abc(2)* & 982 ANINT((r_insert(2) - r_target(2))/abc(2)) 983 RIJ(3) = r_insert(3) - r_target(3) - abc(3)* & 984 ANINT((r_insert(3) - r_target(3))/abc(3)) 985 986 dist = RIJ(1)**2 + RIJ(2)**2 + RIJ(3)**2 987 988 IF (dist .LT. rmin**2 .OR. dist .GT. rmax**2) THEN 989 EXIT 990 ENDIF 991 992 ENDDO 993 ENDIF 994 995 END SUBROUTINE generate_avbmc_insertion 996 997! ***************************************************************************** 998! ************************************************************************************************** 999!> \brief determine the number of cluster present in the given configuration 1000!> based on the rclus value 1001!> \param mc_par the mc parameters for the force env 1002!> \param force_env the force environment containing the coordinates 1003!> \param cluster ... 1004!> \param nchains the number of molecules of each type in the box 1005!> \param nunits the number of interaction sites for each molecule 1006!> \param mol_type an array that indicates the type of each molecule 1007!> \param total_clus ... 1008!> \par 1009!> Original Multiparticle/Cluster Translation paper: 1010!> Orkoulas, Gerassimos, and Athanassios Z. Panagiotopoulos. Free energy and 1011!> phase equilibria for the restricted primitive model of ionic fluids from Monte 1012!> Carlo simulations. J. Chem. Phys. 1994,101.2,1452-1459. 1013!> \author Himanshu Goel 1014! ************************************************************************************************** 1015 1016 SUBROUTINE cluster_search(mc_par, force_env, cluster, nchains, nunits, mol_type, total_clus) 1017 1018 TYPE(mc_simpar_type), POINTER :: mc_par 1019 TYPE(force_env_type), POINTER :: force_env 1020 INTEGER, DIMENSION(:, :), INTENT(INOUT) :: cluster 1021 INTEGER, DIMENSION(:), INTENT(IN) :: nchains, nunits, mol_type 1022 INTEGER, INTENT(INOUT) :: total_clus 1023 1024 CHARACTER(len=*), PARAMETER :: routineN = 'cluster_search', routineP = moduleN//':'//routineN 1025 1026 INTEGER :: counter, handle, imol, iunit, jmol, & 1027 junit, nend, nstart, nunit, nunits_i, & 1028 nunits_j 1029 INTEGER, ALLOCATABLE, DIMENSION(:) :: clusmat, decision 1030 LOGICAL :: lclus 1031 REAL(KIND=dp) :: dx, dy, dz, rclus, rclussquare, rsquare 1032 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xcoord, ycoord, zcoord 1033 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r 1034 REAL(KIND=dp), DIMENSION(1:3) :: abc 1035 TYPE(cell_type), POINTER :: cell 1036 TYPE(cp_subsys_type), POINTER :: oldsys 1037 TYPE(particle_list_type), POINTER :: particles 1038 1039! begin the timing of the subroutine 1040 1041 CALL timeset(routineN, handle) 1042 1043 NULLIFY (oldsys, particles) 1044 1045! Getting Particle Coordinates 1046 1047 CALL force_env_get(force_env, cell=cell, subsys=oldsys) 1048 CALL get_cell(cell, abc=abc) 1049 CALL cp_subsys_get(oldsys, particles=particles) 1050 CALL get_mc_par(mc_par, rclus=rclus) 1051 1052 ALLOCATE (r(1:3, 1:MAXVAL(nunits), 1:SUM(nchains))) 1053 1054! Arranging particles coordinates into an easier matrix 1055 nend = SUM(nchains(:)) 1056 junit = 0 1057 DO imol = 1, nend 1058 nunit = nunits(mol_type(imol)) 1059 DO iunit = 1, nunit 1060 junit = junit + 1 1061 r(1:3, iunit, imol) = particles%els(junit)%r(1:3) 1062 ENDDO 1063 ENDDO 1064 1065 counter = 0 1066 1067! Allocating the size of matrix and decision matrix 1068 ALLOCATE (clusmat(nend), decision(nend)) 1069 1070!Initialize 1071 DO imol = 1, nend 1072 decision(imol) = 0 1073 clusmat(imol) = 0 1074 END DO 1075 1076 rclussquare = rclus*rclus 1077! Starting the cluster count loop 1078 DO WHILE (SUM(decision) .LT. nend) 1079 DO nstart = 1, nend 1080 IF (clusmat(nstart) .EQ. 0) THEN 1081 counter = counter + 1 1082 clusmat(nstart) = counter 1083 EXIT 1084 END IF 1085 END DO 1086 1087 lclus = .TRUE. 1088 DO WHILE (lclus .EQV. .TRUE.) 1089 DO imol = 1, nend 1090 nunits_i = nunits(mol_type(imol)) 1091! Allocating the xcoord,ycoord,zcoord based upon the size of molecule nunits 1092 lclus = .FALSE. 1093 IF (clusmat(imol) .EQ. counter .AND. decision(imol) .EQ. 0) THEN 1094 ALLOCATE (xcoord(nunits_i), ycoord(nunits_i), zcoord(nunits_i)) 1095 decision(imol) = 1 1096 lclus = .TRUE. 1097 DO iunit = 1, nunits_i 1098 xcoord(iunit) = r(1, iunit, imol) 1099 ycoord(iunit) = r(2, iunit, imol) 1100 zcoord(iunit) = r(3, iunit, imol) 1101 END DO 1102 EXIT 1103 END IF 1104 END DO 1105 IF (lclus .EQV. .TRUE.) THEN 1106 DO jmol = 1, nend 1107 nunits_j = nunits(mol_type(jmol)) 1108 IF (clusmat(jmol) .EQ. 0 .AND. decision(jmol) .EQ. 0) THEN 1109!Calculating the distance between atoms 1110 DO iunit = 1, nunits_i 1111 DO junit = 1, nunits_j 1112 dx = xcoord(iunit) - r(1, junit, jmol) 1113 dy = ycoord(iunit) - r(2, junit, jmol) 1114 dz = zcoord(iunit) - r(3, junit, jmol) 1115 dx = dx - abc(1)*ANINT(dx/abc(1)) 1116 dy = dy - abc(2)*ANINT(dy/abc(2)) 1117 dz = dz - abc(3)*ANINT(dz/abc(3)) 1118 rsquare = (dx*dx) + (dy*dy) + (dz*dz) 1119!Checking the distance based on rclus square(rclussq) 1120 IF (rsquare .LT. rclussquare) THEN 1121 clusmat(jmol) = counter 1122 END IF 1123 END DO 1124 END DO 1125 END IF 1126 END DO 1127 DEALLOCATE (xcoord, ycoord, zcoord) 1128 END IF 1129 END DO 1130 END DO 1131 1132!Putting cluster information in a cluster matrix 1133 total_clus = counter 1134 1135 DO imol = 1, counter 1136 DO jmol = 1, nend 1137 IF (imol .EQ. clusmat(jmol)) THEN 1138 cluster(imol, jmol) = jmol 1139 END IF 1140 END DO 1141 END DO 1142 DEALLOCATE (r) 1143 DEALLOCATE (decision) 1144 DEALLOCATE (clusmat) 1145 1146! end the timing 1147 CALL timestop(handle) 1148 1149 END SUBROUTINE cluster_search 1150 1151END MODULE mc_coordinates 1152 1153