mc_coordinates.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief contains miscellaneous subroutines used in the Monte Carlo runs,mostly
!>      geared towards changes in coordinates
!> \author MJM
! **************************************************************************************************
MODULE mc_coordinates
   USE cell_types,                      ONLY: cell_type,&
                                              get_cell
   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
                                              cp_subsys_type
   USE force_env_methods,               ONLY: force_env_calc_energy_force
   USE force_env_types,                 ONLY: force_env_get,&
                                              force_env_type
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: pi
   USE mc_types,                        ONLY: get_mc_molecule_info,&
                                              get_mc_par,&
                                              mc_molecule_info_type,&
                                              mc_simpar_type
   USE message_passing,                 ONLY: mp_comm_type
   USE molecule_types,                  ONLY: molecule_type
   USE parallel_rng_types,              ONLY: rng_stream_type
   USE particle_list_types,             ONLY: particle_list_type
   USE physcon,                         ONLY: angstrom
#include "../../base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   PRIVATE :: generate_avbmc_insertion
 
   PUBLIC :: generate_cbmc_swap_config, &
             get_center_of_mass, mc_coordinate_fold, &
             find_mc_test_molecule, &
             create_discrete_array, &
             check_for_overlap, &
             rotate_molecule, &
             cluster_search
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_coordinates'
 
CONTAINS
 
! **************************************************************************************************
!> \brief looks for overlaps (intermolecular distances less than rmin)
!> \param force_env the force environment containing the coordinates
!> \param nchains the number of molecules of each type in the box
!> \param nunits the number of interaction sites for each molecule
!> \param loverlap returns .TRUE. if atoms overlap
!> \param mol_type an array that indicates the type of each molecule
!> \param cell_length the length of the box...if none is specified,
!>                     it uses the cell found in the force_env
!> \param molecule_number if present, just look for overlaps with this
!>            molecule
!>
!>      Suitable for parallel use.
!> \author MJM
! **************************************************************************************************
   SUBROUTINE check_for_overlap(force_env, nchains, nunits, loverlap, mol_type, &
                                cell_length, molecule_number)
 
      TYPE(force_env_type), POINTER                      :: force_env
      INTEGER, DIMENSION(:), INTENT(IN)                  :: nchains, nunits
      LOGICAL, INTENT(OUT)                               :: loverlap
      INTEGER, DIMENSION(:), INTENT(IN)                  :: mol_type
      REAL(kind=dp), DIMENSION(1:3), INTENT(IN), &
         OPTIONAL                                        :: cell_length
      INTEGER, INTENT(IN), OPTIONAL                      :: molecule_number
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'check_for_overlap'
 
      INTEGER                                            :: handle, imol, iunit, jmol, jstart, &
                                                            junit, nend, nstart, nunit, nunits_i, &
                                                            nunits_j
      LOGICAL                                            :: lall
      REAL(kind=dp)                                      :: dist, rmin
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: r
      REAL(kind=dp), DIMENSION(1:3)                      :: abc, box_length, rij
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_subsys_type), POINTER                      :: oldsys
      TYPE(particle_list_type), POINTER                  :: particles
 
! begin the timing of the subroutine
 
      CALL timeset(routinen, handle)
 
      NULLIFY (oldsys, particles)
 
! initialize some stuff
      loverlap = .false.
      rmin = 1.28558315_dp ! 0.6 angstrom squared
 
! get the particle coordinates and the cell length
      CALL force_env_get(force_env, cell=cell, subsys=oldsys)
      CALL get_cell(cell, abc=abc)
      CALL cp_subsys_get(oldsys, particles=particles)
 
      ALLOCATE (r(1:3, 1:maxval(nunits), 1:sum(nchains)))
 
      IF (PRESENT(cell_length)) THEN
         box_length(1:3) = cell_length(1:3)
      ELSE
         box_length(1:3) = abc(1:3)
      END IF
 
! put the coordinates into an easier matrix to manipulate
      junit = 0
      DO imol = 1, sum(nchains)
         nunit = nunits(mol_type(imol))
         DO iunit = 1, nunit
            junit = junit + 1
            r(1:3, iunit, imol) = particles%els(junit)%r(1:3)
         END DO
      END DO
 
! now let's find the LJ energy between all the oxygens and
! the charge interactions
      lall = .true.
      jstart = 1
      IF (PRESENT(molecule_number)) THEN
         lall = .false.
         nstart = molecule_number
         nend = molecule_number
      ELSE
         nstart = 1
         nend = sum(nchains(:))
      END IF
      DO imol = nstart, nend
         IF (lall) jstart = imol + 1
         nunits_i = nunits(mol_type(imol))
         DO jmol = jstart, sum(nchains(:))
            IF (imol == jmol) cycle
            nunits_j = nunits(mol_type(jmol))
 
            DO iunit = 1, nunits_i
               DO junit = 1, nunits_j
! find the minimum image distance
                  rij(1) = r(1, iunit, imol) - r(1, junit, jmol) - &
                           box_length(1)*anint( &
                           (r(1, iunit, imol) - r(1, junit, jmol))/box_length(1))
                  rij(2) = r(2, iunit, imol) - r(2, junit, jmol) - &
                           box_length(2)*anint( &
                           (r(2, iunit, imol) - r(2, junit, jmol))/box_length(2))
                  rij(3) = r(3, iunit, imol) - r(3, junit, jmol) - &
                           box_length(3)*anint( &
                           (r(3, iunit, imol) - r(3, junit, jmol))/box_length(3))
 
                  dist = rij(1)**2 + rij(2)**2 + rij(3)**2
 
                  IF (dist < rmin) THEN
                     loverlap = .true.
                     DEALLOCATE (r)
 
                     CALL timestop(handle)
                     RETURN
                  END IF
 
               END DO
            END DO
         END DO
      END DO
 
      DEALLOCATE (r)
 
! end the timing
      CALL timestop(handle)
 
   END SUBROUTINE check_for_overlap
 
! **************************************************************************************************
!> \brief calculates the center of mass of a given molecule
!> \param coordinates the coordinates of the atoms in the molecule
!> \param natom the number of atoms in the molecule
!> \param center_of_mass the coordinates of the center of mass
!> \param mass the mass of the atoms in the molecule
!>
!>    Designed for parallel use.
!> \author MJM
! **************************************************************************************************
   SUBROUTINE get_center_of_mass(coordinates, natom, center_of_mass, &
                                 mass)
 
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: coordinates
      INTEGER, INTENT(IN)                                :: natom
      REAL(kind=dp), DIMENSION(1:3), INTENT(OUT)         :: center_of_mass
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: mass
 
      CHARACTER(len=*), PARAMETER :: routinen = 'get_center_of_mass'
 
      INTEGER                                            :: handle, i, iatom
      REAL(kind=dp)                                      :: total_mass
 
! begin the timing of the subroutine
 
      CALL timeset(routinen, handle)
 
      total_mass = sum(mass(1:natom))
      center_of_mass(:) = 0.0e0_dp
 
      DO iatom = 1, natom
         DO i = 1, 3
            center_of_mass(i) = center_of_mass(i) + &
                                mass(iatom)*coordinates(i, iatom)
         END DO
      END DO
 
      center_of_mass(1:3) = center_of_mass(1:3)/total_mass
 
! end the timing
      CALL timestop(handle)
 
   END SUBROUTINE get_center_of_mass
 
! **************************************************************************************************
!> \brief folds all the coordinates into the center simulation box using
!>      a center of mass cutoff
!> \param coordinates the coordinates of the atoms in the system
!> \param nchains_tot the total number of molecules in the box
!> \param mol_type an array that indicates the type of every molecule in the box
!> \param mass the mass of every atom for all molecule types
!> \param nunits the number of interaction sites for each molecule type
!> \param box_length an array for the lengths of the simulation box sides
!>
!>      Designed for parallel use.
!> \author MJM
! **************************************************************************************************
   SUBROUTINE mc_coordinate_fold(coordinates, nchains_tot, mol_type, mass, nunits, box_length)
 
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: coordinates
      INTEGER, INTENT(IN)                                :: nchains_tot
      INTEGER, DIMENSION(:), INTENT(IN)                  :: mol_type
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: mass
      INTEGER, DIMENSION(:), INTENT(IN)                  :: nunits
      REAL(kind=dp), DIMENSION(1:3), INTENT(IN)          :: box_length
 
      CHARACTER(len=*), PARAMETER :: routinen = 'mc_coordinate_fold'
 
      INTEGER                                            :: end_atom, handle, iatom, imolecule, &
                                                            jatom, molecule_type, natoms, &
                                                            start_atom
      REAL(kind=dp), DIMENSION(1:3)                      :: center_of_mass
 
! begin the timing of the subroutine
 
      CALL timeset(routinen, handle)
 
! loop over all molecules
      end_atom = 0
      DO imolecule = 1, nchains_tot
         molecule_type = mol_type(imolecule)
         natoms = nunits(molecule_type)
         start_atom = end_atom + 1
         end_atom = start_atom + natoms - 1
         CALL get_center_of_mass(coordinates(:, start_atom:end_atom), &
                                 natoms, center_of_mass(:), mass(:, molecule_type))
         DO iatom = 1, natoms
            jatom = iatom + start_atom - 1
            coordinates(1, jatom) = coordinates(1, jatom) - &
                                    box_length(1)*floor(center_of_mass(1)/box_length(1))
            coordinates(2, jatom) = coordinates(2, jatom) - &
                                    box_length(2)*floor(center_of_mass(2)/box_length(2))
            coordinates(3, jatom) = coordinates(3, jatom) - &
                                    box_length(3)*floor(center_of_mass(3)/box_length(2))
         END DO
 
      END DO
 
! end the timing
      CALL timestop(handle)
 
   END SUBROUTINE mc_coordinate_fold
 
! **************************************************************************************************
!> \brief takes the last molecule in a force environment and moves it around
!>      to different center of mass positions and orientations, selecting one
!>      based on the rosenbluth weight
!> \param force_env the force environment containing the coordinates
!> \param BETA the value of 1/kT for this simulations, in a.u.
!> \param max_val ...
!> \param min_val ...
!> \param exp_max_val ...
!> \param exp_min_val ...
!> \param nswapmoves the number of desired trial configurations
!> \param rosenbluth_weight the Rosenbluth weight for this set of configs
!> \param start_atom the atom number that the molecule to be swapped starts on
!> \param natoms_tot the total number of interaction sites in the box
!> \param nunits the number of interaction sites for every molecule_type
!> \param nunits_mol ...
!> \param mass the mass for every interaction site of every molecule type
!> \param loverlap the flag that indicates if all of the configs have an
!>                  atomic overlap
!> \param choosen_energy the energy of the chosen config
!> \param old_energy the energy that we subtract from all of the trial
!>        energies to prevent numerical overflows
!> \param ionode indicates if we're on the main CPU
!> \param lremove is this the Rosenbluth weight for a removal box?
!> \param mol_type an array that contains the molecule type for every atom in the box
!> \param nchains the number of molecules of each type in this box
!> \param source the MPI source value, for broadcasts
!> \param group the MPI group value, for broadcasts
!> \param rng_stream the random number stream that we draw from
!> \param avbmc_atom ...
!> \param rmin ...
!> \param rmax ...
!> \param move_type ...
!> \param target_atom ...
!> \par Optional Avbmc Flags
!>      - avbmc_atom: the atom number that serves for the target atom in each
!>        molecule (1 is the first atom in the molecule, etc.)
!>      - rmin: the minimum AVBMC radius for the shell around the target
!>      - rmax: the maximum AVBMC radius for the shell around the target
!>      - move_type: generate configs in the "in" or "out" volume
!>      - target_atom: the number of the avbmc atom in the target molecule
!> \par
!>      Suitable for parallel.
!> \author MJM
! **************************************************************************************************
   SUBROUTINE generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, &
                                        exp_min_val, nswapmoves, rosenbluth_weight, start_atom, natoms_tot, nunits, nunits_mol, &
                                        mass, loverlap, choosen_energy, old_energy, ionode, lremove, mol_type, nchains, source, &
                                        group, rng_stream, avbmc_atom, rmin, rmax, move_type, target_atom)
 
      TYPE(force_env_type), POINTER                      :: force_env
      REAL(kind=dp), INTENT(IN)                          :: beta, max_val, min_val, exp_max_val, &
                                                            exp_min_val
      INTEGER, INTENT(IN)                                :: nswapmoves
      REAL(kind=dp), INTENT(OUT)                         :: rosenbluth_weight
      INTEGER, INTENT(IN)                                :: start_atom, natoms_tot
      INTEGER, DIMENSION(:), INTENT(IN)                  :: nunits
      INTEGER, INTENT(IN)                                :: nunits_mol
      REAL(dp), DIMENSION(1:nunits_mol), INTENT(IN)      :: mass
      LOGICAL, INTENT(OUT)                               :: loverlap
      REAL(kind=dp), INTENT(OUT)                         :: choosen_energy
      REAL(kind=dp), INTENT(IN)                          :: old_energy
      LOGICAL, INTENT(IN)                                :: ionode, lremove
      INTEGER, DIMENSION(:), INTENT(IN)                  :: mol_type, nchains
      INTEGER, INTENT(IN)                                :: source
      TYPE(mp_comm_type)                                 :: group
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
      INTEGER, INTENT(IN), OPTIONAL                      :: avbmc_atom
      REAL(kind=dp), INTENT(IN), OPTIONAL                :: rmin, rmax
      CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: move_type
      INTEGER, INTENT(IN), OPTIONAL                      :: target_atom
 
      CHARACTER(len=*), PARAMETER :: routinen = 'generate_cbmc_swap_config'
 
      INTEGER                                            :: atom_number, choosen, end_atom, handle, &
                                                            i, iatom, imolecule, imove, &
                                                            molecule_number
      LOGICAL                                            :: all_overlaps
      LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: loverlap_array
      REAL(kind=dp)                                      :: bias_energy, exponent, rand, &
                                                            total_running_weight
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: boltz_weights, trial_energy
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: r_old
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: r
      REAL(kind=dp), DIMENSION(1:3)                      :: abc, center_of_mass, diff, r_insert
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_subsys_type), POINTER                      :: oldsys
      TYPE(particle_list_type), POINTER                  :: particles
 
! begin the timing of the subroutine
 
      CALL timeset(routinen, handle)
 
      NULLIFY (oldsys)
! get the particle coordinates and the cell length
      CALL force_env_get(force_env, cell=cell, subsys=oldsys)
      CALL get_cell(cell, abc=abc)
      CALL cp_subsys_get(oldsys, particles=particles)
 
! do some checking to make sure we have all the data we need
      IF (PRESENT(avbmc_atom)) THEN
         IF (.NOT. PRESENT(rmin) .OR. .NOT. PRESENT(rmax) .OR. &
             .NOT. PRESENT(move_type) .OR. .NOT. PRESENT(target_atom)) THEN
            cpabort('AVBMC swap move is missing information!')
         END IF
      END IF
 
      ALLOCATE (r_old(1:3, 1:natoms_tot))
      ALLOCATE (r(1:3, 1:natoms_tot, 1:nswapmoves))
      ALLOCATE (trial_energy(1:nswapmoves))
      ALLOCATE (boltz_weights(1:nswapmoves))
      ALLOCATE (loverlap_array(1:nswapmoves))
 
! initialize the arrays that need it
      loverlap_array(:) = .false.
      loverlap = .false.
      boltz_weights(:) = 0.0_dp
      trial_energy(:) = 0.0_dp
      r(:, :, :) = 0.0_dp
      choosen_energy = 0.0_dp
      rosenbluth_weight = 0.0_dp
 
! save the positions of the molecules
      DO imove = 1, nswapmoves
         DO iatom = 1, natoms_tot
            r(1:3, iatom, imove) = particles%els(iatom)%r(1:3)
         END DO
      END DO
 
! save the remove coordinates
      DO iatom = 1, natoms_tot
         r_old(1:3, iatom) = r(1:3, iatom, 1)
      END DO
 
! figure out the numbers of the first and last atoms in the molecule
      end_atom = start_atom + nunits_mol - 1
! figure out which molecule number we're on
      molecule_number = 0
      atom_number = 1
      DO imolecule = 1, sum(nchains(:))
         IF (atom_number == start_atom) THEN
            molecule_number = imolecule
            EXIT
         END IF
         atom_number = atom_number + nunits(mol_type(imolecule))
      END DO
      IF (molecule_number == 0) CALL cp_abort(__location__, &
                                              'CBMC swap move cannot find which molecule number it needs')
 
      IF (lremove) THEN
         CALL check_for_overlap(force_env, nchains, nunits, loverlap_array(1), &
                                mol_type)
         CALL group%bcast(loverlap_array(1), source)
 
         IF (loverlap_array(1)) THEN
            IF (ionode) THEN
               WRITE (*, *) start_atom, end_atom, natoms_tot
               DO iatom = 1, natoms_tot
                  WRITE (*, *) r(1:3, iatom, 1)
               END DO
            END IF
            cpabort('CBMC swap move found an overlap in the old config')
         END IF
      END IF
 
      DO imove = 1, nswapmoves
 
! drop into serial
         IF (ionode) THEN
 
            IF (PRESENT(avbmc_atom)) THEN
! find an AVBMC insertion point
               CALL generate_avbmc_insertion(rmin, rmax, &
                                             r_old(1:3, target_atom), &
                                             move_type, r_insert(:), abc(:), rng_stream)
 
               DO i = 1, 3
                  diff(i) = r_insert(i) - r_old(i, start_atom + avbmc_atom - 1)
               END DO
 
            ELSE
! find a new insertion point somewhere in the box
               DO i = 1, 3
                  rand = rng_stream%next()
                  r_insert(i) = rand*abc(i)
               END DO
 
! find the center of mass of the insertion molecule
               CALL get_center_of_mass(r(:, start_atom:end_atom, imove), nunits_mol, &
                                       center_of_mass(:), mass(:))
 
! move the molecule to the insertion point
 
               DO i = 1, 3
                  diff(i) = r_insert(i) - center_of_mass(i)
               END DO
 
            END IF
 
            DO iatom = start_atom, end_atom
               r(1:3, iatom, imove) = r(1:3, iatom, imove) + diff(1:3)
            END DO
 
! rotate the molecule...this routine is only made for serial use
            CALL rotate_molecule(r(:, start_atom:end_atom, imove), mass(:), &
                                 nunits_mol, rng_stream)
 
            IF (imove == 1 .AND. lremove) THEN
               DO iatom = 1, natoms_tot
                  r(1:3, iatom, 1) = r_old(1:3, iatom)
               END DO
            END IF
 
         END IF
 
         CALL group%bcast(r(:, :, imove), source)
 
! calculate the energy and boltzman weight of the config
         DO iatom = start_atom, end_atom
            particles%els(iatom)%r(1:3) = r(1:3, iatom, imove)
         END DO
 
         CALL check_for_overlap(force_env, nchains, nunits, loverlap_array(imove), &
                                mol_type, molecule_number=molecule_number)
         IF (loverlap_array(imove)) THEN
            boltz_weights(imove) = 0.0_dp
            cycle
         END IF
 
         CALL force_env_calc_energy_force(force_env, calc_force=.false.)
         CALL force_env_get(force_env, &
                            potential_energy=bias_energy)
 
         trial_energy(imove) = (bias_energy - old_energy)
         exponent = -beta*trial_energy(imove)
 
         IF (exponent > exp_max_val) THEN
            boltz_weights(imove) = max_val
         ELSE IF (exponent < exp_min_val) THEN
            boltz_weights(imove) = min_val
         ELSE
            boltz_weights(imove) = exp(exponent)
         END IF
 
      END DO
 
! now we need to pick a configuration based on the Rosenbluth weight,
! which is just the sum of the Boltzmann weights
      rosenbluth_weight = sum(boltz_weights(:))
      IF (rosenbluth_weight == 0.0_dp .AND. lremove) THEN
! should never have 0.0 for an old weight...causes a divide by zero
! in the acceptance rule
         IF (ionode) THEN
            WRITE (*, *) boltz_weights(1:nswapmoves)
            WRITE (*, *) start_atom, end_atom, lremove
            WRITE (*, *) loverlap_array(1:nswapmoves)
            WRITE (*, *) natoms_tot
            WRITE (*, *)
            DO iatom = 1, natoms_tot
               WRITE (*, *) r(1:3, iatom, 1)*angstrom
            END DO
         END IF
         cpabort('CBMC swap move found a bad old weight')
      END IF
      all_overlaps = .true.
      total_running_weight = 0.0e0_dp
      choosen = 0
      IF (ionode) THEN
         rand = rng_stream%next()
!         CALL random_number(rand)
      END IF
      CALL group%bcast(rand, source)
      DO imove = 1, nswapmoves
         IF (loverlap_array(imove)) cycle
         all_overlaps = .false.
         total_running_weight = total_running_weight + boltz_weights(imove)
         IF (total_running_weight >= rand*rosenbluth_weight) THEN
            choosen = imove
            EXIT
         END IF
      END DO
 
      IF (all_overlaps) THEN
         loverlap = .true.
 
! if this is an old configuration, we always choose the first one...
! this should never be the case, but I'm testing something
         IF (lremove) THEN
            IF (ionode) THEN
               WRITE (*, *) boltz_weights(1:nswapmoves)
               WRITE (*, *) start_atom, end_atom, lremove
               WRITE (*, *) loverlap_array(1:nswapmoves)
               DO iatom = 1, natoms_tot
                  WRITE (*, *) r(1:3, iatom, 1)
               END DO
            END IF
            cpabort('CBMC swap move found all overlaps for the remove config')
         END IF
 
         DEALLOCATE (r_old)
         DEALLOCATE (r)
         DEALLOCATE (trial_energy)
         DEALLOCATE (boltz_weights)
         DEALLOCATE (loverlap_array)
         CALL timestop(handle)
         RETURN
      END IF
 
! make sure a configuration was chosen
      IF (choosen == 0) THEN
         cpabort('CBMC swap move failed to select config')
      END IF
 
! if this is an old configuration, we always choose the first one
      IF (lremove) choosen = 1
 
! set the energy for the configuration
      choosen_energy = trial_energy(choosen)
 
! copy the coordinates to the force environment
      DO iatom = 1, natoms_tot
         particles%els(iatom)%r(1:3) = r(1:3, iatom, choosen)
      END DO
 
      DEALLOCATE (r_old)
      DEALLOCATE (r)
      DEALLOCATE (trial_energy)
      DEALLOCATE (boltz_weights)
      DEALLOCATE (loverlap_array)
 
! end the timing
      CALL timestop(handle)
 
   END SUBROUTINE generate_cbmc_swap_config
 
! **************************************************************************************************
!> \brief rotates a molecule randomly around the center of mass,
!>      sequentially in x, y, and z directions
!> \param r the coordinates of the molecule to rotate
!> \param mass the mass of all the atoms in the molecule
!> \param natoms the number of atoms in the molecule
!> \param rng_stream the stream we pull random numbers from
!>
!>    Use only in serial.
!> \author MJM
! **************************************************************************************************
   SUBROUTINE rotate_molecule(r, mass, natoms, rng_stream)
 
      INTEGER, INTENT(IN)                                :: natoms
      REAL(kind=dp), DIMENSION(1:natoms), INTENT(IN)     :: mass
      REAL(kind=dp), DIMENSION(1:3, 1:natoms), &
         INTENT(INOUT)                                   :: r
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'rotate_molecule'
 
      INTEGER                                            :: handle, iunit
      REAL(kind=dp)                                      :: cosdg, dgamma, rand, rx, rxnew, ry, &
                                                            rynew, rz, rznew, sindg
      REAL(kind=dp), DIMENSION(1:3)                      :: center_of_mass
 
! begin the timing of the subroutine
 
      CALL timeset(routinen, handle)
 
! find the center of mass of the molecule
      CALL get_center_of_mass(r(:, :), natoms, center_of_mass(:), mass(:))
 
! call a random number to figure out how far we're moving
      rand = rng_stream%next()
      dgamma = pi*(rand - 0.5e0_dp)*2.0e0_dp
 
! *** set up the rotation matrix ***
 
      cosdg = cos(dgamma)
      sindg = sin(dgamma)
 
! ***    ROTATE UNITS OF I AROUND X-AXIS ***
 
      DO iunit = 1, natoms
         ry = r(2, iunit) - center_of_mass(2)
         rz = r(3, iunit) - center_of_mass(3)
         rynew = cosdg*ry + sindg*rz
         rznew = cosdg*rz - sindg*ry
 
         r(2, iunit) = rynew + center_of_mass(2)
         r(3, iunit) = rznew + center_of_mass(3)
 
      END DO
 
! ***    ROTATE UNITS OF I AROUND y-AXIS ***
 
      DO iunit = 1, natoms
         rx = r(1, iunit) - center_of_mass(1)
         rz = r(3, iunit) - center_of_mass(3)
         rxnew = cosdg*rx + sindg*rz
         rznew = cosdg*rz - sindg*rx
 
         r(1, iunit) = rxnew + center_of_mass(1)
         r(3, iunit) = rznew + center_of_mass(3)
 
      END DO
 
! ***    ROTATE UNITS OF I AROUND z-AXIS ***
 
      DO iunit = 1, natoms
         rx = r(1, iunit) - center_of_mass(1)
         ry = r(2, iunit) - center_of_mass(2)
         rxnew = cosdg*rx + sindg*ry
         rynew = cosdg*ry - sindg*rx
 
         r(1, iunit) = rxnew + center_of_mass(1)
         r(2, iunit) = rynew + center_of_mass(2)
 
      END DO
 
! end the timing
      CALL timestop(handle)
 
   END SUBROUTINE rotate_molecule
 
! **************************************************************************************************
!> \brief selects a molecule at random to perform a MC move on...you can specify
!>      the box the molecule should be in, its type, both, or neither
!> \param mc_molecule_info the structure that contains some global molecule data
!> \param start_atom the number of the first atom in the chosen molecule in relation
!>        to the force_env it's in
!> \param box_number the box the chosen molecule is in
!> \param molecule_type the type of molecule the chosen molecule is
!> \param rng_stream the stream we pull random numbers from
!> \param box if present, tells the routine which box to grab a molecule from
!> \param molecule_type_old if present, tells the routine which molecule type to select from
!> \author MJM
! **************************************************************************************************
   SUBROUTINE find_mc_test_molecule(mc_molecule_info, start_atom, &
                                    box_number, molecule_type, rng_stream, box, molecule_type_old)
 
      TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
      INTEGER, INTENT(OUT)                               :: start_atom, box_number, molecule_type
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
      INTEGER, INTENT(IN), OPTIONAL                      :: box, molecule_type_old
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'find_mc_test_molecule'
 
      INTEGER                                            :: handle, ibox, imol_type, imolecule, &
                                                            jbox, molecule_number, nchains_tot, &
                                                            start_mol
      INTEGER, DIMENSION(:), POINTER                     :: mol_type, nunits
      INTEGER, DIMENSION(:, :), POINTER                  :: nchains
      REAL(kind=dp)                                      :: rand
 
! begin the timing of the subroutine
 
      CALL timeset(routinen, handle)
 
      NULLIFY (nunits, mol_type, nchains)
      CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, nunits=nunits, &
                                mol_type=mol_type)
 
! initialize the outgoing variables
      start_atom = 0
      box_number = 0
      molecule_type = 0
 
      IF (PRESENT(box) .AND. PRESENT(molecule_type_old)) THEN
! only need to find the atom number the molecule starts on
         rand = rng_stream%next()
         molecule_number = ceiling(rand*real(nchains(molecule_type_old, box), kind=dp))
 
         start_mol = 1
         DO jbox = 1, box - 1
            start_mol = start_mol + sum(nchains(:, jbox))
         END DO
 
! adjust to take into account molecules of other types in the box
         DO imol_type = 1, molecule_type_old - 1
            molecule_number = molecule_number + nchains(imol_type, box)
         END DO
 
         start_atom = 1
         DO imolecule = 1, molecule_number - 1
            start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
         END DO
 
      ELSE IF (PRESENT(box)) THEN
! any molecule in box...need to find molecule type and start atom
         rand = rng_stream%next()
         molecule_number = ceiling(rand*real(sum(nchains(:, box)), kind=dp))
 
         start_mol = 1
         DO jbox = 1, box - 1
            start_mol = start_mol + sum(nchains(:, jbox))
         END DO
 
         molecule_type = mol_type(start_mol + molecule_number - 1)
 
! now the starting atom
         start_atom = 1
         DO imolecule = 1, molecule_number - 1
            start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
         END DO
 
      ELSE IF (PRESENT(molecule_type_old)) THEN
! any molecule of type molecule_type_old...need to find box number and start atom
         rand = rng_stream%next()
         molecule_number = ceiling(rand*real(sum(nchains(molecule_type_old, :)), kind=dp))
 
! find which box it's in
         nchains_tot = 0
         DO ibox = 1, SIZE(nchains(molecule_type_old, :))
            IF (molecule_number <= nchains(molecule_type_old, ibox)) THEN
               box_number = ibox
               EXIT
            END IF
            molecule_number = molecule_number - nchains(molecule_type_old, ibox)
         END DO
 
         start_mol = 1
         DO jbox = 1, box_number - 1
            start_mol = start_mol + sum(nchains(:, jbox))
         END DO
 
! now find the starting atom number
         DO imol_type = 1, molecule_type_old - 1
            molecule_number = molecule_number + nchains(imol_type, box_number)
         END DO
         start_atom = 1
         DO imolecule = 1, molecule_number - 1
            start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
         END DO
 
      ELSE
! no restrictions...need to find all pieces of data
         nchains_tot = 0
         DO ibox = 1, SIZE(nchains(1, :))
            nchains_tot = nchains_tot + sum(nchains(:, ibox))
         END DO
         rand = rng_stream%next()
         molecule_number = ceiling(rand*real(nchains_tot, kind=dp))
 
         molecule_type = mol_type(molecule_number)
 
! now which box it's in
         DO ibox = 1, sum(nchains(1, :))
            IF (molecule_number <= sum(nchains(:, ibox))) THEN
               box_number = ibox
               EXIT
            END IF
            molecule_number = molecule_number - sum(nchains(:, ibox))
         END DO
 
! now find the starting atom number
         start_mol = 1
         DO jbox = 1, box_number - 1
            start_mol = start_mol + sum(nchains(:, jbox))
         END DO
         start_atom = 1
         DO imolecule = 1, molecule_number - 1
            start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
         END DO
 
      END IF
 
! make sure things are good
      IF (PRESENT(box)) box_number = box
      IF (PRESENT(molecule_type_old)) molecule_type = molecule_type_old
 
      cpassert(start_atom /= 0)
      cpassert(box_number /= 0)
      cpassert(molecule_type /= 0)
 
! end the timing
      CALL timestop(handle)
 
   END SUBROUTINE find_mc_test_molecule
 
! **************************************************************************************************
!> \brief generates an array that tells us which sides of the simulation
!>      cell we can increase or decrease using a discrete volume move
!> \param cell the lengths of the sides of the cell
!> \param discrete_array the array that indicates which sides we can move
!> \param step_size the size of the discrete volume move
!>
!>    Suitable for parallel.
!> \author MJM
! **************************************************************************************************
   SUBROUTINE create_discrete_array(cell, discrete_array, step_size)
 
! 1 is for increase, 2 is for decrease
! 1 is for "yes, we can do the move", 0 is for no
 
      REAL(dp), DIMENSION(1:3), INTENT(IN)               :: cell
      INTEGER, DIMENSION(1:3, 1:2), INTENT(OUT)          :: discrete_array
      REAL(dp), INTENT(IN)                               :: step_size
 
      INTEGER                                            :: iside
      REAL(dp)                                           :: high_value, length1, length2, low_value
 
      discrete_array(:, :) = 0
 
      length1 = abs(cell(1) - cell(2))
      length2 = abs(cell(2) - cell(3))
 
! now let's figure out all the different cases
      IF (length1 < 0.01_dp*step_size .AND. &
          length2 < 0.01_dp*step_size) THEN
! all sides are equal, so we can move up or down
         discrete_array(1:3, 1) = 1
         discrete_array(1:3, 2) = 1
      ELSE
 
! find the low value and the high value
         high_value = -1.0_dp
         low_value = cell(1)*cell(2)*cell(3)
         DO iside = 1, 3
            IF (cell(iside) < low_value) low_value = cell(iside)
            IF (cell(iside) > high_value) high_value = cell(iside)
         END DO
         DO iside = 1, 3
! now we see if the value is a high value or a low value...it can only be
! one of the two
            IF (abs(cell(iside) - low_value) < 0.01_dp*step_size) THEN
! low value, we can only increase the cell size
               discrete_array(iside, 1) = 1
               discrete_array(iside, 2) = 0
            ELSE
! high value, we can only decrease the cell size
               discrete_array(iside, 1) = 0
               discrete_array(iside, 2) = 1
            END IF
         END DO
      END IF
 
   END SUBROUTINE create_discrete_array
 
! **************************************************************************************************
!> \brief generates an insertion point in either the "in" or the "out" volume
!>      of a target atom, where the "in" volume is a shell with inner radius
!>       rmin and outer radius rmax
!> \param rmin the minimum AVBMC radius for the shell around the target
!> \param rmax the maximum AVBMC radius for the shell around the target
!> \param r_target the coordinates of the target atom
!> \param move_type generate configs in the "in" or "out" volume
!> \param r_insert the output insertion site
!> \param abc the lengths of the sides of the simulation box
!> \param rng_stream the random number stream that we draw from
!>
!>      Use only in serial.
!> \author MJM
! **************************************************************************************************
   SUBROUTINE generate_avbmc_insertion(rmin, rmax, r_target, &
                                       move_type, r_insert, abc, rng_stream)
 
      REAL(kind=dp), INTENT(IN)                          :: rmin, rmax
      REAL(kind=dp), DIMENSION(1:3), INTENT(IN)          :: r_target
      CHARACTER(LEN=*), INTENT(IN)                       :: move_type
      REAL(kind=dp), DIMENSION(1:3), INTENT(OUT)         :: r_insert
      REAL(kind=dp), DIMENSION(1:3), INTENT(IN)          :: abc
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
      INTEGER                                            :: i
      REAL(dp)                                           :: dist, eta_1, eta_2, eta_sq, rand
      REAL(dp), DIMENSION(1:3)                           :: rij
 
      r_insert(1:3) = 0.0_dp
 
      IF (move_type == 'in') THEN
! generate a random unit vector, from Allen and Tildesley
         DO
            eta_1 = rng_stream%next()
            eta_2 = rng_stream%next()
            eta_sq = eta_1**2 + eta_2**2
            IF (eta_sq < 1.0_dp) THEN
               r_insert(1) = 2.0_dp*eta_1*sqrt(1.0_dp - eta_sq)
               r_insert(2) = 2.0_dp*eta_2*sqrt(1.0_dp - eta_sq)
               r_insert(3) = 1.0_dp - 2.0_dp*eta_sq
               EXIT
            END IF
         END DO
 
! now scale that vector to be within the "in" region
         rand = rng_stream%next()
         r_insert(1:3) = r_insert(1:3)*(rand*(rmax**3 - rmin**3) + rmin**3)** &
                         (1.0_dp/3.0_dp)
 
         r_insert(1:3) = r_target(1:3) + r_insert(1:3)
      ELSE
 
! find a new insertion point somewhere in the box
         DO
            DO i = 1, 3
               rand = rng_stream%next()
               r_insert(i) = rand*abc(i)
            END DO
 
! make sure it's not in the "in" region
            rij(1) = r_insert(1) - r_target(1) - abc(1)* &
                     anint((r_insert(1) - r_target(1))/abc(1))
            rij(2) = r_insert(2) - r_target(2) - abc(2)* &
                     anint((r_insert(2) - r_target(2))/abc(2))
            rij(3) = r_insert(3) - r_target(3) - abc(3)* &
                     anint((r_insert(3) - r_target(3))/abc(3))
 
            dist = rij(1)**2 + rij(2)**2 + rij(3)**2
 
            IF (dist < rmin**2 .OR. dist > rmax**2) THEN
               EXIT
            END IF
 
         END DO
      END IF
 
   END SUBROUTINE generate_avbmc_insertion
 
! *****************************************************************************
! **************************************************************************************************
!> \brief determine the number of cluster present in the given configuration
!>   based on the rclus value
!> \param mc_par the mc parameters for the force env
!> \param force_env the force environment containing the coordinates
!> \param cluster ...
!> \param nchains the number of molecules of each type in the box
!> \param nunits the number of interaction sites for each molecule
!> \param mol_type an array that indicates the type of each molecule
!> \param total_clus ...
!> \par
!>    Original Multiparticle/Cluster Translation paper:
!> Orkoulas, Gerassimos, and Athanassios Z. Panagiotopoulos. Free energy and
!> phase equilibria for the restricted primitive model of ionic fluids from Monte
!> Carlo simulations. J. Chem. Phys. 1994,101.2,1452-1459.
!> \author Himanshu Goel
! **************************************************************************************************
 
   SUBROUTINE cluster_search(mc_par, force_env, cluster, nchains, nunits, mol_type, total_clus)
 
      TYPE(mc_simpar_type), POINTER                      :: mc_par
      TYPE(force_env_type), POINTER                      :: force_env
      INTEGER, DIMENSION(:, :), INTENT(INOUT)            :: cluster
      INTEGER, DIMENSION(:), INTENT(IN)                  :: nchains, nunits, mol_type
      INTEGER, INTENT(INOUT)                             :: total_clus
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'cluster_search'
 
      INTEGER                                            :: counter, handle, imol, iunit, jmol, &
                                                            junit, nend, nstart, nunit, nunits_i, &
                                                            nunits_j
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: clusmat, decision
      LOGICAL                                            :: lclus
      REAL(kind=dp)                                      :: dx, dy, dz, rclus, rclussquare, rsquare
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: xcoord, ycoord, zcoord
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: r
      REAL(kind=dp), DIMENSION(1:3)                      :: abc
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_subsys_type), POINTER                      :: oldsys
      TYPE(particle_list_type), POINTER                  :: particles
 
! begin the timing of the subroutine
 
      CALL timeset(routinen, handle)
 
      NULLIFY (oldsys, particles)
 
! Getting Particle Coordinates
 
      CALL force_env_get(force_env, cell=cell, subsys=oldsys)
      CALL get_cell(cell, abc=abc)
      CALL cp_subsys_get(oldsys, particles=particles)
      CALL get_mc_par(mc_par, rclus=rclus)
 
      ALLOCATE (r(1:3, 1:maxval(nunits), 1:sum(nchains)))
 
! Arranging particles coordinates into an easier matrix
      nend = sum(nchains(:))
      junit = 0
      DO imol = 1, nend
         nunit = nunits(mol_type(imol))
         DO iunit = 1, nunit
            junit = junit + 1
            r(1:3, iunit, imol) = particles%els(junit)%r(1:3)
         END DO
      END DO
 
      counter = 0
 
! Allocating the size of matrix and decision matrix
      ALLOCATE (clusmat(nend), decision(nend))
 
!Initialize
      DO imol = 1, nend
         decision(imol) = 0
         clusmat(imol) = 0
      END DO
 
      rclussquare = rclus*rclus
! Starting the cluster count loop
      DO WHILE (sum(decision) < nend)
         DO nstart = 1, nend
            IF (clusmat(nstart) == 0) THEN
               counter = counter + 1
               clusmat(nstart) = counter
               EXIT
            END IF
         END DO
 
         lclus = .true.
         DO WHILE (lclus .EQV. .true.)
            DO imol = 1, nend
               nunits_i = nunits(mol_type(imol))
! Allocating the xcoord,ycoord,zcoord based upon the size of molecule nunits
               lclus = .false.
               IF (clusmat(imol) == counter .AND. decision(imol) == 0) THEN
                  ALLOCATE (xcoord(nunits_i), ycoord(nunits_i), zcoord(nunits_i))
                  decision(imol) = 1
                  lclus = .true.
                  DO iunit = 1, nunits_i
                     xcoord(iunit) = r(1, iunit, imol)
                     ycoord(iunit) = r(2, iunit, imol)
                     zcoord(iunit) = r(3, iunit, imol)
                  END DO
                  EXIT
               END IF
            END DO
            IF (lclus .EQV. .true.) THEN
               DO jmol = 1, nend
                  nunits_j = nunits(mol_type(jmol))
                  IF (clusmat(jmol) == 0 .AND. decision(jmol) == 0) THEN
!Calculating the distance between atoms
                     DO iunit = 1, nunits_i
                        DO junit = 1, nunits_j
                           dx = xcoord(iunit) - r(1, junit, jmol)
                           dy = ycoord(iunit) - r(2, junit, jmol)
                           dz = zcoord(iunit) - r(3, junit, jmol)
                           dx = dx - abc(1)*anint(dx/abc(1))
                           dy = dy - abc(2)*anint(dy/abc(2))
                           dz = dz - abc(3)*anint(dz/abc(3))
                           rsquare = (dx*dx) + (dy*dy) + (dz*dz)
!Checking the distance based on rclus square(rclussq)
                           IF (rsquare < rclussquare) THEN
                              clusmat(jmol) = counter
                           END IF
                        END DO
                     END DO
                  END IF
               END DO
               DEALLOCATE (xcoord, ycoord, zcoord)
            END IF
         END DO
      END DO
 
!Putting cluster information in a cluster matrix
      total_clus = counter
 
      DO imol = 1, counter
         DO jmol = 1, nend
            IF (imol == clusmat(jmol)) THEN
               cluster(imol, jmol) = jmol
            END IF
         END DO
      END DO
      DEALLOCATE (r)
      DEALLOCATE (decision)
      DEALLOCATE (clusmat)
 
! end the timing
      CALL timestop(handle)
 
   END SUBROUTINE cluster_search
 
END MODULE mc_coordinates