LCOV - code coverage report
Current view: top level - src/motion/mc - mc_moves.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 69.6 % 1001 697
Test Date: 2025-07-25 12:55:17 Functions: 81.8 % 11 9

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief the various moves in Monte Carlo (MC) simulations, including
      10              : !>      change of internal conformation, translation of a molecule, rotation
      11              : !>      of a molecule, and changing the size of the simulation box
      12              : !> \par History
      13              : !>      none
      14              : !> \author Matthew J. McGrath  (10.16.2003)
      15              : ! **************************************************************************************************
      16              : MODULE mc_moves
      17              :    USE atomic_kind_types,               ONLY: get_atomic_kind
      18              :    USE cell_methods,                    ONLY: cell_create
      19              :    USE cell_types,                      ONLY: cell_clone,&
      20              :                                               cell_release,&
      21              :                                               cell_type,&
      22              :                                               get_cell
      23              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      24              :                                               cp_logger_get_default_io_unit,&
      25              :                                               cp_logger_type
      26              :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      27              :                                               cp_subsys_set,&
      28              :                                               cp_subsys_type
      29              :    USE force_env_methods,               ONLY: force_env_calc_energy_force
      30              :    USE force_env_types,                 ONLY: force_env_get,&
      31              :                                               force_env_type,&
      32              :                                               use_fist_force
      33              :    USE global_types,                    ONLY: global_environment_type
      34              :    USE kinds,                           ONLY: default_string_length,&
      35              :                                               dp
      36              :    USE mathconstants,                   ONLY: pi
      37              :    USE mc_coordinates,                  ONLY: check_for_overlap,&
      38              :                                               cluster_search,&
      39              :                                               create_discrete_array,&
      40              :                                               generate_cbmc_swap_config,&
      41              :                                               get_center_of_mass
      42              :    USE mc_types,                        ONLY: get_mc_molecule_info,&
      43              :                                               get_mc_par,&
      44              :                                               mc_ekin_type,&
      45              :                                               mc_molecule_info_type,&
      46              :                                               mc_moves_type,&
      47              :                                               mc_simpar_type
      48              :    USE md_run,                          ONLY: qs_mol_dyn
      49              :    USE message_passing,                 ONLY: mp_comm_type
      50              :    USE molecule_kind_list_types,        ONLY: molecule_kind_list_type
      51              :    USE molecule_kind_types,             ONLY: bend_type,&
      52              :                                               bond_type,&
      53              :                                               get_molecule_kind,&
      54              :                                               molecule_kind_type,&
      55              :                                               torsion_type
      56              :    USE parallel_rng_types,              ONLY: rng_stream_type
      57              :    USE particle_list_types,             ONLY: particle_list_type
      58              :    USE physcon,                         ONLY: angstrom
      59              : #include "../../base/base_uses.f90"
      60              : 
      61              :    IMPLICIT NONE
      62              : 
      63              :    PRIVATE
      64              : 
      65              :    PRIVATE  :: change_bond_angle, change_bond_length, depth_first_search, &
      66              :                change_dihedral
      67              : 
      68              : ! *** Global parameters ***
      69              : 
      70              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_moves'
      71              : 
      72              :    PUBLIC :: mc_conformation_change, mc_molecule_translation, &
      73              :              mc_molecule_rotation, mc_volume_move, mc_avbmc_move, &
      74              :              mc_hmc_move, mc_cluster_translation
      75              : 
      76              : CONTAINS
      77              : 
      78              : ! **************************************************************************************************
      79              : !> \brief essentially performs a depth-first search of the molecule structure
      80              : !>      to find all atoms connected to a specific atom excluding one branch...
      81              : !>      for instance, if water is labelled 1-2-3 for O-H-H, calling this
      82              : !>      routine with current_atom=1,avoid_atom=2 returns the array
      83              : !>      atom=(0,0,1)
      84              : !> \param current_atom the atom whose connections we're looking at
      85              : !> \param avoid_atom the atom whose direction the search is not supposed to go
      86              : !> \param connectivity an array telling us the neighbors of all atoms
      87              : !> \param atom the array that tells us if one can get to a given atom by
      88              : !>        starting at current_atom and not going through avoid_atom...0 is no,
      89              : !>        1 is yes
      90              : !> \author MJM
      91              : ! **************************************************************************************************
      92         1324 :    RECURSIVE SUBROUTINE depth_first_search(current_atom, avoid_atom, &
      93         1324 :                                            connectivity, atom)
      94              : 
      95              :       INTEGER, INTENT(IN)                                :: current_atom, avoid_atom
      96              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: connectivity
      97              :       INTEGER, DIMENSION(:), INTENT(INOUT)               :: atom
      98              : 
      99              :       INTEGER                                            :: iatom
     100              : 
     101         2960 :       DO iatom = 1, 6
     102         2960 :          IF (connectivity(iatom, current_atom) .NE. 0) THEN
     103         1636 :             IF (connectivity(iatom, current_atom) .NE. avoid_atom) THEN
     104          312 :                atom(connectivity(iatom, current_atom)) = 1
     105              :                CALL depth_first_search(connectivity(iatom, current_atom), &
     106          312 :                                        current_atom, connectivity, atom)
     107              :             END IF
     108              :          ELSE
     109              :             RETURN
     110              :          END IF
     111              :       END DO
     112              : 
     113              :    END SUBROUTINE depth_first_search
     114              : 
     115              : ! **************************************************************************************************
     116              : !> \brief performs either a bond or angle change move for a given molecule
     117              : !> \param mc_par the mc parameters for the force env
     118              : !> \param force_env the force environment used in the move
     119              : !> \param bias_env the force environment used to bias the move, if any (it may
     120              : !>            be null if lbias=.false. in mc_par)
     121              : !> \param moves the structure that keeps track of how many moves have been
     122              : !>               accepted/rejected
     123              : !> \param move_updates the structure that keeps track of how many moves have
     124              : !>               been accepted/rejected since the last time the displacements
     125              : !>               were updated
     126              : !> \param start_atom the number of the molecule's first atom, assuming the rest
     127              : !>        of the atoms follow sequentially
     128              : !> \param molecule_type the type of the molecule we're moving
     129              : !> \param box_number the box the molecule is in
     130              : !> \param bias_energy the biased energy of the system before the move
     131              : !> \param move_type dictates what kind of conformational change we do
     132              : !> \param lreject set to .true. if there is an overlap
     133              : !> \param rng_stream the random number stream that we draw from
     134              : !> \author MJM
     135              : ! **************************************************************************************************
     136          506 :    SUBROUTINE mc_conformation_change(mc_par, force_env, bias_env, moves, &
     137              :                                      move_updates, start_atom, molecule_type, box_number, &
     138              :                                      bias_energy, move_type, lreject, &
     139              :                                      rng_stream)
     140              : 
     141              :       TYPE(mc_simpar_type), POINTER                      :: mc_par
     142              :       TYPE(force_env_type), POINTER                      :: force_env, bias_env
     143              :       TYPE(mc_moves_type), POINTER                       :: moves, move_updates
     144              :       INTEGER, INTENT(IN)                                :: start_atom, molecule_type, box_number
     145              :       REAL(KIND=dp), INTENT(INOUT)                       :: bias_energy
     146              :       CHARACTER(LEN=*), INTENT(IN)                       :: move_type
     147              :       LOGICAL, INTENT(OUT)                               :: lreject
     148              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     149              : 
     150              :       CHARACTER(len=*), PARAMETER :: routineN = 'mc_conformation_change'
     151              : 
     152              :       CHARACTER(default_string_length)                   :: name
     153              :       CHARACTER(default_string_length), DIMENSION(:), &
     154          506 :          POINTER                                         :: names
     155              :       INTEGER :: atom_number, end_atom, end_mol, handle, imol_type, imolecule, ipart, jbox, &
     156              :          molecule_number, nunits_mol, source, start_mol
     157          506 :       INTEGER, DIMENSION(:), POINTER                     :: mol_type, nunits
     158          506 :       INTEGER, DIMENSION(:, :), POINTER                  :: nchains
     159              :       LOGICAL                                            :: ionode, lbias, loverlap
     160              :       REAL(KIND=dp)                                      :: BETA, bias_energy_new, bias_energy_old, &
     161              :                                                             dis_length, exp_max_val, exp_min_val, &
     162              :                                                             rand, value, w
     163              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: r_new, r_old
     164              :       TYPE(cp_subsys_type), POINTER                      :: subsys
     165              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
     166              :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
     167              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind, molecule_kind_test
     168              :       TYPE(mp_comm_type)                                 :: group
     169              :       TYPE(particle_list_type), POINTER                  :: particles
     170              : 
     171              : ! begin the timing of the subroutine
     172              : 
     173          506 :       CALL timeset(routineN, handle)
     174              : 
     175              : ! nullify some pointers
     176          506 :       NULLIFY (particles, subsys, molecule_kinds, molecule_kind, &
     177          506 :                molecule_kind_test)
     178              : 
     179              : ! get a bunch of stuff from mc_par
     180              :       CALL get_mc_par(mc_par, lbias=lbias, mc_molecule_info=mc_molecule_info, &
     181              :                       BETA=BETA, exp_max_val=exp_max_val, &
     182          506 :                       exp_min_val=exp_min_val, group=group, source=source, ionode=ionode)
     183              :       CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, nunits=nunits, &
     184          506 :                                 mol_type=mol_type, names=names)
     185              : 
     186              : ! do some allocation
     187          506 :       nunits_mol = nunits(molecule_type)
     188         1518 :       ALLOCATE (r_old(1:3, 1:nunits_mol))
     189         1012 :       ALLOCATE (r_new(1:3, 1:nunits_mol))
     190              : 
     191              : ! find out some bounds for mol_type
     192          506 :       start_mol = 1
     193          506 :       DO jbox = 1, box_number - 1
     194          506 :          start_mol = start_mol + SUM(nchains(:, jbox))
     195              :       END DO
     196         1518 :       end_mol = start_mol + SUM(nchains(:, box_number)) - 1
     197              : 
     198              : ! figure out which molecule number we are
     199          506 :       end_atom = start_atom + nunits_mol - 1
     200          506 :       molecule_number = 0
     201          506 :       atom_number = 1
     202         2974 :       DO imolecule = 1, SUM(nchains(:, box_number))
     203         1962 :          IF (atom_number == start_atom) THEN
     204          506 :             molecule_number = imolecule
     205          506 :             EXIT
     206              :          END IF
     207         1456 :          atom_number = atom_number + nunits(mol_type(imolecule + start_mol - 1))
     208              :       END DO
     209          506 :       IF (molecule_number == 0) CPABORT('Cannot find the molecule number')
     210              : 
     211              : ! are we biasing this move?
     212          506 :       IF (lbias) THEN
     213              : 
     214              : ! grab the coordinates
     215          420 :          CALL force_env_get(bias_env, subsys=subsys)
     216              : ! save the energy
     217          420 :          bias_energy_old = bias_energy
     218              : 
     219              :       ELSE
     220              : 
     221              : ! grab the coordinates
     222           86 :          CALL force_env_get(force_env, subsys=subsys)
     223              :       END IF
     224              : 
     225              : ! now find the molecule type associated with this guy
     226              :       CALL cp_subsys_get(subsys, &
     227          506 :                          particles=particles, molecule_kinds=molecule_kinds)
     228          506 :       DO imol_type = 1, SIZE(molecule_kinds%els(:))
     229          506 :          molecule_kind_test => molecule_kinds%els(imol_type)
     230          506 :          CALL get_molecule_kind(molecule_kind_test, name=name)
     231          506 :          IF (TRIM(ADJUSTL(name)) == TRIM(ADJUSTL(names(molecule_type)))) THEN
     232          506 :             molecule_kind => molecule_kinds%els(imol_type)
     233          506 :             EXIT
     234              :          END IF
     235              :       END DO
     236              : 
     237              : ! save the coordinates
     238         2024 :       DO ipart = start_atom, end_atom
     239         6578 :          r_old(1:3, ipart - start_atom + 1) = particles%els(ipart)%r(1:3)
     240              :       END DO
     241              : 
     242          506 :       IF (.NOT. ASSOCIATED(molecule_kind)) CPABORT('Cannot find the molecule type')
     243              : ! do the move
     244          506 :       IF (move_type == 'bond') THEN
     245              : 
     246              : ! record the attempt
     247          312 :          moves%bond%attempts = moves%bond%attempts + 1
     248          312 :          move_updates%bond%attempts = move_updates%bond%attempts + 1
     249          312 :          moves%bias_bond%attempts = moves%bias_bond%attempts + 1
     250          312 :          move_updates%bias_bond%attempts = move_updates%bias_bond%attempts + 1
     251          312 :          IF (.NOT. lbias) THEN
     252           48 :             moves%bond%qsuccesses = moves%bond%qsuccesses + 1
     253              :             move_updates%bond%qsuccesses = &
     254           48 :                move_updates%bond%qsuccesses + 1
     255           48 :             moves%bias_bond%qsuccesses = moves%bias_bond%qsuccesses + 1
     256              :             move_updates%bias_bond%qsuccesses = &
     257           48 :                move_updates%bias_bond%qsuccesses + 1
     258              :          END IF
     259              : 
     260              : ! do the move
     261              :          CALL change_bond_length(r_old, r_new, mc_par, molecule_type, &
     262          312 :                                  molecule_kind, dis_length, particles, rng_stream)
     263              : 
     264          194 :       ELSEIF (move_type == 'angle') THEN
     265              : 
     266              : ! record the attempt
     267          194 :          moves%angle%attempts = moves%angle%attempts + 1
     268          194 :          move_updates%angle%attempts = move_updates%angle%attempts + 1
     269          194 :          moves%bias_angle%attempts = moves%bias_angle%attempts + 1
     270          194 :          move_updates%bias_angle%attempts = move_updates%bias_angle%attempts + 1
     271          194 :          IF (.NOT. lbias) THEN
     272           38 :             moves%angle%qsuccesses = moves%angle%qsuccesses + 1
     273              :             move_updates%angle%qsuccesses = &
     274           38 :                move_updates%angle%qsuccesses + 1
     275           38 :             moves%bias_angle%qsuccesses = moves%bias_angle%qsuccesses + 1
     276              :             move_updates%bias_angle%qsuccesses = &
     277           38 :                move_updates%bias_angle%qsuccesses + 1
     278              :          END IF
     279              : 
     280              : ! do the move
     281              :          CALL change_bond_angle(r_old, r_new, mc_par, molecule_type, &
     282          194 :                                 molecule_kind, particles, rng_stream)
     283          194 :          dis_length = 1.0E0_dp
     284              :       ELSE
     285              : ! record the attempt
     286            0 :          moves%dihedral%attempts = moves%dihedral%attempts + 1
     287            0 :          move_updates%dihedral%attempts = move_updates%dihedral%attempts + 1
     288            0 :          moves%bias_dihedral%attempts = moves%bias_dihedral%attempts + 1
     289            0 :          move_updates%bias_dihedral%attempts = move_updates%bias_dihedral%attempts + 1
     290            0 :          IF (.NOT. lbias) THEN
     291            0 :             moves%dihedral%qsuccesses = moves%dihedral%qsuccesses + 1
     292              :             move_updates%dihedral%qsuccesses = &
     293            0 :                move_updates%dihedral%qsuccesses + 1
     294            0 :             moves%bias_dihedral%qsuccesses = moves%bias_dihedral%qsuccesses + 1
     295              :             move_updates%bias_dihedral%qsuccesses = &
     296            0 :                move_updates%bias_dihedral%qsuccesses + 1
     297              :          END IF
     298              : 
     299              : ! do the move
     300              :          CALL change_dihedral(r_old, r_new, mc_par, molecule_type, &
     301            0 :                               molecule_kind, particles, rng_stream)
     302            0 :          dis_length = 1.0E0_dp
     303              : 
     304              :       END IF
     305              : 
     306              : ! set the coordinates
     307         2024 :       DO ipart = start_atom, end_atom
     308         6578 :          particles%els(ipart)%r(1:3) = r_new(1:3, ipart - start_atom + 1)
     309              :       END DO
     310              : 
     311              : ! check for overlap
     312          506 :       lreject = .FALSE.
     313          506 :       IF (lbias) THEN
     314              :          CALL check_for_overlap(bias_env, nchains(:, box_number), &
     315              :                                 nunits(:), loverlap, mol_type(start_mol:end_mol), &
     316          420 :                                 molecule_number=molecule_number)
     317              :       ELSE
     318              :          CALL check_for_overlap(force_env, nchains(:, box_number), &
     319              :                                 nunits(:), loverlap, mol_type(start_mol:end_mol), &
     320           86 :                                 molecule_number=molecule_number)
     321           86 :          IF (loverlap) lreject = .TRUE.
     322              :       END IF
     323              : 
     324              : ! if we're biasing classical, check for acceptance
     325          506 :       IF (lbias) THEN
     326              : 
     327              : ! here's where we bias the moves
     328              : 
     329          420 :          IF (loverlap) THEN
     330              :             w = 0.0E0_dp
     331              :          ELSE
     332          420 :             CALL force_env_calc_energy_force(bias_env, calc_force=.FALSE.)
     333              :             CALL force_env_get(bias_env, &
     334          420 :                                potential_energy=bias_energy_new)
     335              : ! accept or reject the move based on the Metropolis rule with a
     336              : ! correction factor for the change in phase space...dis_length is
     337              : ! made unitless in change_bond_length
     338          420 :             value = -BETA*(bias_energy_new - bias_energy_old)
     339          420 :             IF (value .GT. exp_max_val) THEN
     340              :                w = 10.0_dp
     341          420 :             ELSEIF (value .LT. exp_min_val) THEN
     342              :                w = 0.0_dp
     343              :             ELSE
     344          420 :                w = EXP(value)*dis_length**2
     345              :             END IF
     346              : 
     347              :          END IF
     348              : 
     349          420 :          IF (w .GE. 1.0E0_dp) THEN
     350          194 :             w = 1.0E0_dp
     351          194 :             rand = 0.0E0_dp
     352              :          ELSE
     353          226 :             IF (ionode) THEN
     354          113 :                rand = rng_stream%next()
     355              :             END IF
     356          226 :             CALL group%bcast(rand, source)
     357              :          END IF
     358              : 
     359          420 :          IF (rand .LT. w) THEN
     360              : 
     361              : ! accept the move
     362          252 :             IF (move_type == 'bond') THEN
     363          140 :                moves%bond%qsuccesses = moves%bond%qsuccesses + 1
     364              :                move_updates%bond%successes = &
     365          140 :                   move_updates%bond%successes + 1
     366          140 :                moves%bias_bond%successes = moves%bias_bond%successes + 1
     367              :                move_updates%bias_bond%successes = &
     368          140 :                   move_updates%bias_bond%successes + 1
     369          112 :             ELSEIF (move_type == 'angle') THEN
     370          112 :                moves%angle%qsuccesses = moves%angle%qsuccesses + 1
     371              :                move_updates%angle%successes = &
     372          112 :                   move_updates%angle%successes + 1
     373          112 :                moves%bias_angle%successes = moves%bias_angle%successes + 1
     374              :                move_updates%bias_angle%successes = &
     375          112 :                   move_updates%bias_angle%successes + 1
     376              :             ELSE
     377            0 :                moves%dihedral%qsuccesses = moves%dihedral%qsuccesses + 1
     378              :                move_updates%dihedral%successes = &
     379            0 :                   move_updates%dihedral%successes + 1
     380            0 :                moves%bias_dihedral%successes = moves%bias_dihedral%successes + 1
     381              :                move_updates%bias_dihedral%successes = &
     382            0 :                   move_updates%bias_dihedral%successes + 1
     383              :             END IF
     384              : 
     385              :             bias_energy = bias_energy + bias_energy_new - &
     386          252 :                           bias_energy_old
     387              : 
     388              :          ELSE
     389              : 
     390              : ! reject the move
     391              : ! restore the coordinates
     392          168 :             CALL force_env_get(bias_env, subsys=subsys)
     393          168 :             CALL cp_subsys_get(subsys, particles=particles)
     394          672 :             DO ipart = start_atom, end_atom
     395         2184 :                particles%els(ipart)%r(1:3) = r_old(1:3, ipart - start_atom + 1)
     396              :             END DO
     397          168 :             CALL cp_subsys_set(subsys, particles=particles)
     398              : 
     399              :          END IF
     400              : 
     401              :       END IF
     402              : 
     403              : ! deallocate some stuff
     404          506 :       DEALLOCATE (r_old)
     405          506 :       DEALLOCATE (r_new)
     406              : 
     407              : ! end the timing
     408          506 :       CALL timestop(handle)
     409              : 
     410          506 :    END SUBROUTINE mc_conformation_change
     411              : 
     412              : ! **************************************************************************************************
     413              : !> \brief translates the given molecule randomly in either the x,y, or z direction
     414              : !> \param mc_par the mc parameters for the force env
     415              : !> \param force_env the force environment used in the move
     416              : !> \param bias_env the force environment used to bias the move, if any (it may
     417              : !>            be null if lbias=.false. in mc_par)
     418              : !> \param moves the structure that keeps track of how many moves have been
     419              : !>               accepted/rejected
     420              : !> \param move_updates the structure that keeps track of how many moves have
     421              : !>               been accepted/rejected since the last time the displacements
     422              : !>               were updated
     423              : !> \param start_atom the number of the molecule's first atom, assuming the rest of
     424              : !>        the atoms follow sequentially
     425              : !> \param box_number the box the molecule is in
     426              : !> \param bias_energy the biased energy of the system before the move
     427              : !> \param molecule_type the type of molecule we're moving
     428              : !> \param lreject set to .true. if there is an overlap
     429              : !> \param rng_stream the random number stream that we draw from
     430              : !> \author MJM
     431              : ! **************************************************************************************************
     432          624 :    SUBROUTINE mc_molecule_translation(mc_par, force_env, bias_env, moves, &
     433              :                                       move_updates, start_atom, box_number, &
     434              :                                       bias_energy, molecule_type, &
     435              :                                       lreject, rng_stream)
     436              : 
     437              :       TYPE(mc_simpar_type), POINTER                      :: mc_par
     438              :       TYPE(force_env_type), POINTER                      :: force_env, bias_env
     439              :       TYPE(mc_moves_type), POINTER                       :: moves, move_updates
     440              :       INTEGER, INTENT(IN)                                :: start_atom, box_number
     441              :       REAL(KIND=dp), INTENT(INOUT)                       :: bias_energy
     442              :       INTEGER, INTENT(IN)                                :: molecule_type
     443              :       LOGICAL, INTENT(OUT)                               :: lreject
     444              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     445              : 
     446              :       CHARACTER(len=*), PARAMETER :: routineN = 'mc_molecule_translation'
     447              : 
     448              :       INTEGER :: atom_number, end_atom, end_mol, handle, imolecule, ipart, iparticle, jbox, &
     449              :          molecule_number, move_direction, nunits_mol, source, start_mol
     450          624 :       INTEGER, DIMENSION(:), POINTER                     :: mol_type, nunits, nunits_tot
     451          624 :       INTEGER, DIMENSION(:, :), POINTER                  :: nchains
     452              :       LOGICAL                                            :: ionode, lbias, loverlap
     453          624 :       REAL(dp), DIMENSION(:), POINTER                    :: rmtrans
     454              :       REAL(KIND=dp)                                      :: BETA, bias_energy_new, bias_energy_old, &
     455              :                                                             dis_mol, exp_max_val, exp_min_val, &
     456              :                                                             rand, value, w
     457          624 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: r_old
     458              :       TYPE(cp_subsys_type), POINTER                      :: subsys
     459              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
     460              :       TYPE(mp_comm_type)                                 :: group
     461              :       TYPE(particle_list_type), POINTER                  :: particles
     462              : 
     463              : !   *** Local Counters ***
     464              : ! begin the timing of the subroutine
     465              : 
     466          624 :       CALL timeset(routineN, handle)
     467              : 
     468              : ! nullify some pointers
     469          624 :       NULLIFY (particles, subsys)
     470              : 
     471              : ! get a bunch of stuff from mc_par
     472              :       CALL get_mc_par(mc_par, lbias=lbias, &
     473              :                       BETA=BETA, exp_max_val=exp_max_val, &
     474              :                       exp_min_val=exp_min_val, rmtrans=rmtrans, ionode=ionode, source=source, &
     475          624 :                       group=group, mc_molecule_info=mc_molecule_info)
     476              :       CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot, &
     477          624 :                                 nchains=nchains, nunits=nunits, mol_type=mol_type)
     478              : 
     479              : ! find out some bounds for mol_type
     480          624 :       start_mol = 1
     481          640 :       DO jbox = 1, box_number - 1
     482          672 :          start_mol = start_mol + SUM(nchains(:, jbox))
     483              :       END DO
     484         1872 :       end_mol = start_mol + SUM(nchains(:, box_number)) - 1
     485              : 
     486              : ! do some allocation
     487         1872 :       ALLOCATE (r_old(1:3, 1:nunits_tot(box_number)))
     488              : 
     489              : ! find the index of the last atom of this molecule, and the molecule number
     490          624 :       nunits_mol = nunits(molecule_type)
     491          624 :       end_atom = start_atom + nunits_mol - 1
     492          624 :       molecule_number = 0
     493          624 :       atom_number = 1
     494         5770 :       DO imolecule = 1, SUM(nchains(:, box_number))
     495         4522 :          IF (atom_number == start_atom) THEN
     496          624 :             molecule_number = imolecule
     497          624 :             EXIT
     498              :          END IF
     499         3898 :          atom_number = atom_number + nunits(mol_type(imolecule + start_mol - 1))
     500              :       END DO
     501          624 :       IF (molecule_number == 0) CPABORT('Cannot find the molecule number')
     502              : 
     503              : ! are we biasing this move?
     504          624 :       IF (lbias) THEN
     505              : 
     506              : ! grab the coordinates
     507          528 :          CALL force_env_get(bias_env, subsys=subsys)
     508          528 :          CALL cp_subsys_get(subsys, particles=particles)
     509              : 
     510              : ! save the coordinates
     511        14710 :          DO ipart = 1, nunits_tot(box_number)
     512        57256 :             r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
     513              :          END DO
     514              : 
     515              : ! save the energy
     516          528 :          bias_energy_old = bias_energy
     517              : 
     518              :       ELSE
     519              : 
     520              : ! grab the coordinates
     521           96 :          CALL force_env_get(force_env, subsys=subsys)
     522           96 :          CALL cp_subsys_get(subsys, particles=particles)
     523              :       END IF
     524              : 
     525              : ! record the attempt
     526          624 :       moves%trans%attempts = moves%trans%attempts + 1
     527          624 :       move_updates%trans%attempts = move_updates%trans%attempts + 1
     528          624 :       moves%bias_trans%attempts = moves%bias_trans%attempts + 1
     529          624 :       move_updates%bias_trans%attempts = move_updates%bias_trans%attempts + 1
     530          624 :       IF (.NOT. lbias) THEN
     531           96 :          moves%trans%qsuccesses = moves%trans%qsuccesses + 1
     532           96 :          move_updates%trans%qsuccesses = move_updates%trans%qsuccesses + 1
     533           96 :          moves%bias_trans%qsuccesses = moves%bias_trans%qsuccesses + 1
     534           96 :          move_updates%bias_trans%qsuccesses = move_updates%bias_trans%qsuccesses + 1
     535              :       END IF
     536              : 
     537              : ! move one molecule in the system
     538              : 
     539              : ! call a random number to figure out which direction we're moving
     540          624 :       IF (ionode) rand = rng_stream%next()
     541          624 :       CALL group%bcast(rand, source)
     542              :       ! 1,2,3 with equal prob
     543          624 :       move_direction = INT(3*rand) + 1
     544              : 
     545              : ! call a random number to figure out how far we're moving
     546          624 :       IF (ionode) rand = rng_stream%next()
     547          624 :       CALL group%bcast(rand, source)
     548          624 :       dis_mol = rmtrans(molecule_type)*(rand - 0.5E0_dp)*2.0E0_dp
     549              : 
     550              : ! do the move
     551         1896 :       DO iparticle = start_atom, end_atom
     552              :          particles%els(iparticle)%r(move_direction) = &
     553         1896 :             particles%els(iparticle)%r(move_direction) + dis_mol
     554              :       END DO
     555          624 :       CALL cp_subsys_set(subsys, particles=particles)
     556              : 
     557              : ! figure out if there is any overlap...need the number of the molecule
     558          624 :       lreject = .FALSE.
     559          624 :       IF (lbias) THEN
     560              :          CALL check_for_overlap(bias_env, nchains(:, box_number), &
     561              :                                 nunits(:), loverlap, mol_type(start_mol:end_mol), &
     562          528 :                                 molecule_number=molecule_number)
     563              :       ELSE
     564              :          CALL check_for_overlap(force_env, nchains(:, box_number), &
     565              :                                 nunits(:), loverlap, mol_type(start_mol:end_mol), &
     566           96 :                                 molecule_number=molecule_number)
     567           96 :          IF (loverlap) lreject = .TRUE.
     568              :       END IF
     569              : 
     570              : ! if we're biasing with a cheaper potential, check for acceptance
     571          624 :       IF (lbias) THEN
     572              : 
     573              : ! here's where we bias the moves
     574          528 :          IF (loverlap) THEN
     575              :             w = 0.0E0_dp
     576              :          ELSE
     577          528 :             CALL force_env_calc_energy_force(bias_env, calc_force=.FALSE.)
     578              :             CALL force_env_get(bias_env, &
     579          528 :                                potential_energy=bias_energy_new)
     580              : ! accept or reject the move based on the Metropolis rule
     581          528 :             value = -BETA*(bias_energy_new - bias_energy_old)
     582          528 :             IF (value .GT. exp_max_val) THEN
     583              :                w = 10.0_dp
     584          528 :             ELSEIF (value .LT. exp_min_val) THEN
     585              :                w = 0.0_dp
     586              :             ELSE
     587          528 :                w = EXP(value)
     588              :             END IF
     589              : 
     590              :          END IF
     591              : 
     592          528 :          IF (w .GE. 1.0E0_dp) THEN
     593          258 :             w = 1.0E0_dp
     594          258 :             rand = 0.0E0_dp
     595              :          ELSE
     596          270 :             IF (ionode) rand = rng_stream%next()
     597          270 :             CALL group%bcast(rand, source)
     598              :          END IF
     599              : 
     600          528 :          IF (rand .LT. w) THEN
     601              : 
     602              : ! accept the move
     603          454 :             moves%bias_trans%successes = moves%bias_trans%successes + 1
     604          454 :             move_updates%bias_trans%successes = move_updates%bias_trans%successes + 1
     605          454 :             moves%trans%qsuccesses = moves%trans%qsuccesses + 1
     606              :             move_updates%trans%successes = &
     607          454 :                move_updates%trans%successes + 1
     608          454 :             moves%qtrans_dis = moves%qtrans_dis + ABS(dis_mol)
     609              :             bias_energy = bias_energy + bias_energy_new - &
     610          454 :                           bias_energy_old
     611              : 
     612              :          ELSE
     613              : 
     614              : ! reject the move
     615              : ! restore the coordinates
     616           74 :             CALL force_env_get(bias_env, subsys=subsys)
     617           74 :             CALL cp_subsys_get(subsys, particles=particles)
     618         2072 :             DO ipart = 1, nunits_tot(box_number)
     619         8066 :                particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
     620              :             END DO
     621           74 :             CALL cp_subsys_set(subsys, particles=particles)
     622              : 
     623              :          END IF
     624              : 
     625              :       END IF
     626              : 
     627              : ! deallocate some stuff
     628          624 :       DEALLOCATE (r_old)
     629              : 
     630              : ! end the timing
     631          624 :       CALL timestop(handle)
     632              : 
     633          624 :    END SUBROUTINE mc_molecule_translation
     634              : 
     635              : ! **************************************************************************************************
     636              : !> \brief rotates the given molecule randomly around the x,y, or z axis...
     637              : !>      only works for water at the moment
     638              : !> \param mc_par the mc parameters for the force env
     639              : !> \param force_env the force environment used in the move
     640              : !> \param bias_env the force environment used to bias the move, if any (it may
     641              : !>            be null if lbias=.false. in mc_par)
     642              : !> \param moves the structure that keeps track of how many moves have been
     643              : !>               accepted/rejected
     644              : !> \param move_updates the structure that keeps track of how many moves have
     645              : !>               been accepted/rejected since the last time the displacements
     646              : !>               were updated
     647              : !> \param box_number the box the molecule is in
     648              : !> \param start_atom the number of the molecule's first atom, assuming the rest of
     649              : !>        the atoms follow sequentially
     650              : !> \param molecule_type the type of molecule we're moving
     651              : !> \param bias_energy the biased energy of the system before the move
     652              : !> \param lreject set to .true. if there is an overlap
     653              : !> \param rng_stream the random number stream that we draw from
     654              : !> \author MJM
     655              : ! **************************************************************************************************
     656          502 :    SUBROUTINE mc_molecule_rotation(mc_par, force_env, bias_env, moves, &
     657              :                                    move_updates, box_number, &
     658              :                                    start_atom, molecule_type, bias_energy, lreject, &
     659              :                                    rng_stream)
     660              : 
     661              :       TYPE(mc_simpar_type), POINTER                      :: mc_par
     662              :       TYPE(force_env_type), POINTER                      :: force_env, bias_env
     663              :       TYPE(mc_moves_type), POINTER                       :: moves, move_updates
     664              :       INTEGER, INTENT(IN)                                :: box_number, start_atom, molecule_type
     665              :       REAL(KIND=dp), INTENT(INOUT)                       :: bias_energy
     666              :       LOGICAL, INTENT(OUT)                               :: lreject
     667              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     668              : 
     669              :       CHARACTER(len=*), PARAMETER :: routineN = 'mc_molecule_rotation'
     670              : 
     671              :       INTEGER :: atom_number, dir, end_atom, end_mol, handle, ii, imolecule, ipart, iunit, jbox, &
     672              :          molecule_number, nunits_mol, source, start_mol
     673          502 :       INTEGER, DIMENSION(:), POINTER                     :: mol_type, nunits, nunits_tot
     674          502 :       INTEGER, DIMENSION(:, :), POINTER                  :: nchains
     675              :       LOGICAL                                            :: ionode, lbias, loverlap, lx, ly
     676          502 :       REAL(dp), DIMENSION(:), POINTER                    :: rmrot
     677          502 :       REAL(dp), DIMENSION(:, :), POINTER                 :: mass
     678              :       REAL(KIND=dp) :: BETA, bias_energy_new, bias_energy_old, cosdg, dgamma, exp_max_val, &
     679              :          exp_min_val, masstot, nxcm, nycm, nzcm, rand, rx, rxnew, ry, rynew, rz, rznew, sindg, &
     680              :          value, w
     681          502 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: r_old
     682              :       TYPE(cp_subsys_type), POINTER                      :: subsys
     683              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
     684              :       TYPE(mp_comm_type)                                 :: group
     685              :       TYPE(particle_list_type), POINTER                  :: particles
     686              : 
     687              : ! begin the timing of the subroutine
     688              : 
     689          502 :       CALL timeset(routineN, handle)
     690              : 
     691          502 :       NULLIFY (rmrot, subsys, particles)
     692              : 
     693              : ! get a bunch of stuff from mc_par
     694              :       CALL get_mc_par(mc_par, lbias=lbias, &
     695              :                       BETA=BETA, exp_max_val=exp_max_val, &
     696              :                       exp_min_val=exp_min_val, rmrot=rmrot, mc_molecule_info=mc_molecule_info, &
     697          502 :                       ionode=ionode, group=group, source=source)
     698              :       CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits, &
     699              :                                 nunits_tot=nunits_tot, nchains=nchains, mass=mass, &
     700          502 :                                 mol_type=mol_type)
     701              : 
     702              : ! figure out some bounds for mol_type
     703          502 :       start_mol = 1
     704          516 :       DO jbox = 1, box_number - 1
     705          544 :          start_mol = start_mol + SUM(nchains(:, jbox))
     706              :       END DO
     707         1506 :       end_mol = start_mol + SUM(nchains(:, box_number)) - 1
     708              : 
     709          502 :       nunits_mol = nunits(molecule_type)
     710              : 
     711              : ! nullify some pointers
     712          502 :       NULLIFY (particles, subsys)
     713              : 
     714              : ! do some allocation
     715         1506 :       ALLOCATE (r_old(1:3, 1:nunits_tot(box_number)))
     716              : 
     717              : ! initialize some stuff
     718          502 :       lx = .FALSE.
     719          502 :       ly = .FALSE.
     720              : 
     721              : ! determine what the final atom in the molecule is numbered, and which
     722              : ! molecule number this is
     723          502 :       end_atom = start_atom + nunits_mol - 1
     724          502 :       molecule_number = 0
     725          502 :       atom_number = 1
     726         2956 :       DO imolecule = 1, SUM(nchains(:, box_number))
     727         1952 :          IF (atom_number == start_atom) THEN
     728          502 :             molecule_number = imolecule
     729          502 :             EXIT
     730              :          END IF
     731         1450 :          atom_number = atom_number + nunits(mol_type(imolecule + start_mol - 1))
     732              :       END DO
     733          502 :       IF (molecule_number == 0) CPABORT('Cannot find the molecule number')
     734              : 
     735              : ! are we biasing this move?
     736          502 :       IF (lbias) THEN
     737              : 
     738              : ! grab the coordinates
     739          424 :          CALL force_env_get(bias_env, subsys=subsys)
     740          424 :          CALL cp_subsys_get(subsys, particles=particles)
     741              : 
     742              : ! save the coordinates
     743        11808 :          DO ipart = 1, nunits_tot(box_number)
     744        45960 :             r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
     745              :          END DO
     746              : 
     747              : ! save the energy
     748          424 :          bias_energy_old = bias_energy
     749              : 
     750              :       ELSE
     751              : 
     752              : ! grab the coordinates
     753           78 :          CALL force_env_get(force_env, subsys=subsys)
     754           78 :          CALL cp_subsys_get(subsys, particles=particles)
     755              :       END IF
     756              : 
     757              : ! grab the masses
     758         2008 :       masstot = SUM(mass(1:nunits(molecule_type), molecule_type))
     759              : 
     760              : ! record the attempt
     761          502 :       moves%bias_rot%attempts = moves%bias_rot%attempts + 1
     762          502 :       move_updates%bias_rot%attempts = move_updates%bias_rot%attempts + 1
     763          502 :       moves%rot%attempts = moves%rot%attempts + 1
     764          502 :       move_updates%rot%attempts = move_updates%rot%attempts + 1
     765          502 :       IF (.NOT. lbias) THEN
     766           78 :          moves%rot%qsuccesses = moves%rot%qsuccesses + 1
     767           78 :          move_updates%rot%qsuccesses = move_updates%rot%qsuccesses + 1
     768           78 :          moves%bias_rot%qsuccesses = moves%bias_rot%qsuccesses + 1
     769           78 :          move_updates%bias_rot%qsuccesses = move_updates%bias_rot%qsuccesses + 1
     770              :       END IF
     771              : 
     772              : ! rotate one molecule in the system
     773              : 
     774              : ! call a random number to figure out which direction we're moving
     775          502 :       IF (ionode) rand = rng_stream%next()
     776              : !      CALL RANDOM_NUMBER(rand)
     777          502 :       CALL group%bcast(rand, source)
     778              :       ! 1,2,3 with equal prob
     779          502 :       dir = INT(3*rand) + 1
     780              : 
     781          502 :       IF (dir .EQ. 1) THEN
     782              :          lx = .TRUE.
     783          334 :       ELSEIF (dir .EQ. 2) THEN
     784          176 :          ly = .TRUE.
     785              :       END IF
     786              : 
     787              : ! Determine new center of mass for chain i by finding the sum
     788              : ! of m*r for each unit, then dividing by the total mass of the chain
     789          502 :       nxcm = 0.0E0_dp
     790          502 :       nycm = 0.0E0_dp
     791          502 :       nzcm = 0.0E0_dp
     792         2008 :       DO ii = 1, nunits_mol
     793         1506 :          nxcm = nxcm + particles%els(start_atom - 1 + ii)%r(1)*mass(ii, molecule_type)
     794         1506 :          nycm = nycm + particles%els(start_atom - 1 + ii)%r(2)*mass(ii, molecule_type)
     795         2008 :          nzcm = nzcm + particles%els(start_atom - 1 + ii)%r(3)*mass(ii, molecule_type)
     796              :       END DO
     797          502 :       nxcm = nxcm/masstot
     798          502 :       nycm = nycm/masstot
     799          502 :       nzcm = nzcm/masstot
     800              : 
     801              : ! call a random number to figure out how far we're moving
     802          502 :       IF (ionode) rand = rng_stream%next()
     803          502 :       CALL group%bcast(rand, source)
     804          502 :       dgamma = rmrot(molecule_type)*(rand - 0.5E0_dp)*2.0E0_dp
     805              : 
     806              : ! *** set up the rotation matrix ***
     807              : 
     808          502 :       cosdg = COS(dgamma)
     809          502 :       sindg = SIN(dgamma)
     810              : 
     811          502 :       IF (lx) THEN
     812              : 
     813              : ! ***    ROTATE UNITS OF I AROUND X-AXIS ***
     814              : 
     815          672 :          DO iunit = start_atom, end_atom
     816          504 :             ry = particles%els(iunit)%r(2) - nycm
     817          504 :             rz = particles%els(iunit)%r(3) - nzcm
     818          504 :             rynew = cosdg*ry - sindg*rz
     819          504 :             rznew = cosdg*rz + sindg*ry
     820              : 
     821          504 :             particles%els(iunit)%r(2) = rynew + nycm
     822          672 :             particles%els(iunit)%r(3) = rznew + nzcm
     823              : 
     824              :          END DO
     825          334 :       ELSEIF (ly) THEN
     826              : 
     827              : ! ***    ROTATE UNITS OF I AROUND y-AXIS ***
     828              : 
     829          704 :          DO iunit = start_atom, end_atom
     830          528 :             rx = particles%els(iunit)%r(1) - nxcm
     831          528 :             rz = particles%els(iunit)%r(3) - nzcm
     832          528 :             rxnew = cosdg*rx + sindg*rz
     833          528 :             rznew = cosdg*rz - sindg*rx
     834              : 
     835          528 :             particles%els(iunit)%r(1) = rxnew + nxcm
     836          704 :             particles%els(iunit)%r(3) = rznew + nzcm
     837              : 
     838              :          END DO
     839              : 
     840              :       ELSE
     841              : 
     842              : ! ***    ROTATE UNITS OF I AROUND z-AXIS ***
     843              : 
     844          632 :          DO iunit = start_atom, end_atom
     845          474 :             rx = particles%els(iunit)%r(1) - nxcm
     846          474 :             ry = particles%els(iunit)%r(2) - nycm
     847              : 
     848          474 :             rxnew = cosdg*rx - sindg*ry
     849          474 :             rynew = cosdg*ry + sindg*rx
     850              : 
     851          474 :             particles%els(iunit)%r(1) = rxnew + nxcm
     852          632 :             particles%els(iunit)%r(2) = rynew + nycm
     853              : 
     854              :          END DO
     855              : 
     856              :       END IF
     857          502 :       CALL cp_subsys_set(subsys, particles=particles)
     858              : 
     859              : ! check for overlap
     860          502 :       lreject = .FALSE.
     861          502 :       IF (lbias) THEN
     862              :          CALL check_for_overlap(bias_env, nchains(:, box_number), &
     863              :                                 nunits(:), loverlap, mol_type(start_mol:end_mol), &
     864          424 :                                 molecule_number=molecule_number)
     865              :       ELSE
     866              :          CALL check_for_overlap(force_env, nchains(:, box_number), &
     867              :                                 nunits(:), loverlap, mol_type(start_mol:end_mol), &
     868           78 :                                 molecule_number=molecule_number)
     869           78 :          IF (loverlap) lreject = .TRUE.
     870              :       END IF
     871              : 
     872              : ! if we're biasing classical, check for acceptance
     873          502 :       IF (lbias) THEN
     874              : 
     875              : ! here's where we bias the moves
     876              : 
     877          424 :          IF (loverlap) THEN
     878              :             w = 0.0E0_dp
     879              :          ELSE
     880          424 :             CALL force_env_calc_energy_force(bias_env, calc_force=.FALSE.)
     881              :             CALL force_env_get(bias_env, &
     882          424 :                                potential_energy=bias_energy_new)
     883              : ! accept or reject the move based on the Metropolis rule
     884          424 :             value = -BETA*(bias_energy_new - bias_energy_old)
     885          424 :             IF (value .GT. exp_max_val) THEN
     886              :                w = 10.0_dp
     887          424 :             ELSEIF (value .LT. exp_min_val) THEN
     888              :                w = 0.0_dp
     889              :             ELSE
     890          424 :                w = EXP(value)
     891              :             END IF
     892              : 
     893              :          END IF
     894              : 
     895          424 :          IF (w .GE. 1.0E0_dp) THEN
     896          180 :             w = 1.0E0_dp
     897          180 :             rand = 0.0E0_dp
     898              :          ELSE
     899          244 :             IF (ionode) rand = rng_stream%next()
     900          244 :             CALL group%bcast(rand, source)
     901              :          END IF
     902              : 
     903          424 :          IF (rand .LT. w) THEN
     904              : 
     905              : ! accept the move
     906          340 :             moves%bias_rot%successes = moves%bias_rot%successes + 1
     907          340 :             move_updates%bias_rot%successes = move_updates%bias_rot%successes + 1
     908          340 :             moves%rot%qsuccesses = moves%rot%qsuccesses + 1
     909          340 :             move_updates%rot%successes = move_updates%rot%successes + 1
     910              :             bias_energy = bias_energy + bias_energy_new - &
     911          340 :                           bias_energy_old
     912              : 
     913              :          ELSE
     914              : 
     915              : ! reject the move
     916              : ! restore the coordinates
     917           84 :             CALL force_env_get(bias_env, subsys=subsys)
     918           84 :             CALL cp_subsys_get(subsys, particles=particles)
     919         2316 :             DO ipart = 1, nunits_tot(box_number)
     920         9012 :                particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
     921              :             END DO
     922           84 :             CALL cp_subsys_set(subsys, particles=particles)
     923              : 
     924              :          END IF
     925              : 
     926              :       END IF
     927              : 
     928              : ! deallocate some stuff
     929          502 :       DEALLOCATE (r_old)
     930              : 
     931              : ! end the timing
     932          502 :       CALL timestop(handle)
     933              : 
     934          502 :    END SUBROUTINE mc_molecule_rotation
     935              : 
     936              : ! **************************************************************************************************
     937              : !> \brief performs a Monte Carlo move that alters the volume of the simulation box
     938              : !> \param mc_par the mc parameters for the force env
     939              : !> \param force_env the force environment whose cell we're changing
     940              : !> \param moves the structure that keeps track of how many moves have been
     941              : !>               accepted/rejected
     942              : !> \param move_updates the structure that keeps track of how many moves have
     943              : !>               been accepted/rejected since the last time the displacements
     944              : !>               were updated
     945              : !> \param old_energy the energy of the last accepted move involving an
     946              : !>                    unbiased calculation
     947              : !> \param box_number the box we're changing the volume of
     948              : !> \param energy_check the running total of how much the energy has changed
     949              : !>                      since the initial configuration
     950              : !> \param r_old the coordinates of the last accepted move involving an
     951              : !>               unbiased calculation
     952              : !> \param iw the unit number that writes to the screen
     953              : !> \param discrete_array tells use which volumes we can do for the discrete
     954              : !>            case
     955              : !> \param rng_stream the random number stream that we draw from
     956              : !> \author MJM
     957              : !> \note     Designed for parallel use.
     958              : ! **************************************************************************************************
     959           34 :    SUBROUTINE mc_volume_move(mc_par, force_env, moves, move_updates, &
     960              :                              old_energy, box_number, &
     961           34 :                              energy_check, r_old, iw, discrete_array, rng_stream)
     962              : 
     963              :       TYPE(mc_simpar_type), POINTER                      :: mc_par
     964              :       TYPE(force_env_type), POINTER                      :: force_env
     965              :       TYPE(mc_moves_type), POINTER                       :: moves, move_updates
     966              :       REAL(KIND=dp), INTENT(INOUT)                       :: old_energy
     967              :       INTEGER, INTENT(IN)                                :: box_number
     968              :       REAL(KIND=dp), INTENT(INOUT)                       :: energy_check
     969              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: r_old
     970              :       INTEGER, INTENT(IN)                                :: iw
     971              :       INTEGER, DIMENSION(1:3, 1:2), INTENT(INOUT)        :: discrete_array
     972              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     973              : 
     974              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'mc_volume_move'
     975              : 
     976              :       CHARACTER(LEN=200)                                 :: fft_lib
     977              :       CHARACTER(LEN=40)                                  :: dat_file
     978              :       INTEGER :: cl, end_atom, end_mol, handle, iatom, idim, imolecule, iside, iside_change, &
     979              :          iunit, jbox, nunits_mol, output_unit, print_level, source, start_atom, start_mol
     980           34 :       INTEGER, DIMENSION(:), POINTER                     :: mol_type, nunits, nunits_tot
     981           34 :       INTEGER, DIMENSION(:, :), POINTER                  :: nchains
     982              :       LOGICAL                                            :: ionode, ldiscrete, lincrease, loverlap, &
     983              :                                                             ltoo_small
     984           34 :       REAL(dp), DIMENSION(:, :), POINTER                 :: mass
     985              :       REAL(KIND=dp) :: BETA, discrete_step, energy_term, exp_max_val, exp_min_val, new_energy, &
     986              :          pressure, pressure_term, rand, rcut, rmvolume, temp_var, value, vol_dis, volume_term, w
     987           34 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: r
     988              :       REAL(KIND=dp), DIMENSION(1:3)                      :: abc, center_of_mass, center_of_mass_new, &
     989              :                                                             diff, new_cell_length, old_cell_length
     990              :       REAL(KIND=dp), DIMENSION(1:3, 1:3)                 :: hmat_test
     991              :       TYPE(cell_type), POINTER                           :: cell, cell_old, cell_test
     992              :       TYPE(cp_logger_type), POINTER                      :: logger
     993              :       TYPE(cp_subsys_type), POINTER                      :: oldsys
     994              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
     995              :       TYPE(mp_comm_type)                                 :: group
     996              :       TYPE(particle_list_type), POINTER                  :: particles_old
     997              : 
     998              : ! begin the timing of the subroutine
     999              : 
    1000           34 :       CALL timeset(routineN, handle)
    1001              : 
    1002              : ! get a bunch of stuff from mc_par
    1003              :       CALL get_mc_par(mc_par, ionode=ionode, &
    1004              :                       BETA=BETA, exp_max_val=exp_max_val, &
    1005              :                       exp_min_val=exp_min_val, source=source, group=group, &
    1006              :                       dat_file=dat_file, rmvolume=rmvolume, pressure=pressure, cl=cl, &
    1007              :                       fft_lib=fft_lib, discrete_step=discrete_step, &
    1008           34 :                       ldiscrete=ldiscrete, mc_molecule_info=mc_molecule_info)
    1009              :       CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, &
    1010              :                                 nunits=nunits, nunits_tot=nunits_tot, mol_type=mol_type, &
    1011           34 :                                 mass=mass)
    1012              : ! figure out some bounds for mol_type
    1013           34 :       start_mol = 1
    1014           46 :       DO jbox = 1, box_number - 1
    1015           70 :          start_mol = start_mol + SUM(nchains(:, jbox))
    1016              :       END DO
    1017           94 :       end_mol = start_mol + SUM(nchains(:, box_number)) - 1
    1018              : 
    1019           34 :       print_level = 1 ! hack, printlevel is for print_keys
    1020              : 
    1021              : ! nullify some pointers
    1022           34 :       NULLIFY (particles_old, cell_old, oldsys, cell_test, cell)
    1023              : 
    1024              : ! do some allocation
    1025          102 :       ALLOCATE (r(1:3, 1:nunits_tot(box_number)))
    1026              : 
    1027              : ! record the attempt
    1028           34 :       moves%volume%attempts = moves%volume%attempts + 1
    1029           34 :       move_updates%volume%attempts = move_updates%volume%attempts + 1
    1030              : 
    1031              : ! now let's grab the cell length and particle positions
    1032           34 :       CALL force_env_get(force_env, subsys=oldsys, cell=cell)
    1033           34 :       CALL get_cell(cell, abc=abc)
    1034           34 :       CALL cell_create(cell_old)
    1035           34 :       CALL cell_clone(cell, cell_old, tag="CELL_OLD")
    1036           34 :       CALL cp_subsys_get(oldsys, particles=particles_old)
    1037              : 
    1038              : ! find the old cell length
    1039           34 :       old_cell_length(1) = abc(1)
    1040           34 :       old_cell_length(2) = abc(2)
    1041           34 :       old_cell_length(3) = abc(3)
    1042              : 
    1043              : ! save the old coordinates
    1044          760 :       DO iatom = 1, nunits_tot(box_number)
    1045         2938 :          r(1:3, iatom) = particles_old%els(iatom)%r(1:3)
    1046              :       END DO
    1047              : 
    1048              : ! now do the move
    1049              : 
    1050              : ! call a random number to figure out how far we're moving
    1051           34 :       IF (ionode) rand = rng_stream%next()
    1052           34 :       CALL group%bcast(rand, source)
    1053              : 
    1054              : ! find the test cell lengths for the discrete volume move
    1055           34 :       IF (ldiscrete) THEN
    1056            0 :          IF (rand .LT. 0.5_dp) THEN
    1057              :             lincrease = .TRUE.
    1058              :          ELSE
    1059            0 :             lincrease = .FALSE.
    1060              :          END IF
    1061              : 
    1062            0 :          new_cell_length(1:3) = old_cell_length(1:3)
    1063              : 
    1064              : ! if we're increasing the volume, we need to find a side we can increase
    1065            0 :          IF (lincrease) THEN
    1066              :             DO
    1067            0 :                IF (ionode) rand = rng_stream%next()
    1068            0 :                CALL group%bcast(rand, source)
    1069            0 :                iside_change = CEILING(3.0_dp*rand)
    1070            0 :                IF (discrete_array(iside_change, 1) .EQ. 1) THEN
    1071              :                   new_cell_length(iside_change) = &
    1072            0 :                      new_cell_length(iside_change) + discrete_step
    1073              :                   EXIT
    1074              :                END IF
    1075              :             END DO
    1076              :          ELSE
    1077              :             DO
    1078            0 :                IF (ionode) rand = rng_stream%next()
    1079            0 :                CALL group%bcast(rand, source)
    1080            0 :                iside_change = CEILING(3.0_dp*rand)
    1081            0 :                IF (discrete_array(iside_change, 2) .EQ. 1) THEN
    1082              :                   new_cell_length(iside_change) = &
    1083            0 :                      new_cell_length(iside_change) - discrete_step
    1084            0 :                   EXIT
    1085              :                END IF
    1086              :             END DO
    1087              :          END IF
    1088              :          vol_dis = (new_cell_length(1)*new_cell_length(2)*new_cell_length(3)) &
    1089            0 :                    - old_cell_length(1)*old_cell_length(2)*old_cell_length(3)
    1090              :       ELSE
    1091              : ! now for the not discrete volume move
    1092              : !!!!!!!!!!!!!!!! for E_V curves
    1093           34 :          vol_dis = rmvolume*(rand - 0.5E0_dp)*2.0E0_dp
    1094              : !         WRITE(output_unit,*) '************************ be sure to change back!',&
    1095              : !                 old_cell_length(1),14.64_dp/angstrom
    1096              : !         vol_dis=-56.423592_dp/angstrom**3
    1097              : !         IF(old_cell_length(1) .LE. 14.64_dp/angstrom) THEN
    1098              : !            vol_dis=0.0_dp
    1099              : !            WRITE(output_unit,*) 'Found the correct box length!'
    1100              : !         ENDIF
    1101              : 
    1102              :          temp_var = vol_dis + &
    1103              :                     old_cell_length(1)*old_cell_length(2)* &
    1104           34 :                     old_cell_length(3)
    1105              : 
    1106           34 :          IF (temp_var .LE. 0.0E0_dp) THEN
    1107            0 :             loverlap = .TRUE. ! cannot have a negative volume
    1108              :          ELSE
    1109           34 :             new_cell_length(1) = (temp_var)**(1.0E0_dp/3.0E0_dp)
    1110           34 :             new_cell_length(2) = new_cell_length(1)
    1111           34 :             new_cell_length(3) = new_cell_length(1)
    1112           34 :             loverlap = .FALSE.
    1113              :          END IF
    1114              :       END IF
    1115           34 :       CALL group%bcast(loverlap, source)
    1116              : 
    1117           34 :       IF (loverlap) THEN
    1118              : ! deallocate some stuff
    1119            0 :          DEALLOCATE (r)
    1120            0 :          logger => cp_get_default_logger()
    1121            0 :          output_unit = cp_logger_get_default_io_unit(logger)
    1122            0 :          IF (output_unit > 0) WRITE (output_unit, *) &
    1123            0 :             "Volume move rejected because we tried to make too small of box.", vol_dis
    1124              : !     end the timing
    1125            0 :          CALL timestop(handle)
    1126            0 :          RETURN
    1127              :       END IF
    1128              : 
    1129              : ! now we need to make the new cell
    1130           34 :       hmat_test(:, :) = 0.0e0_dp
    1131           34 :       hmat_test(1, 1) = new_cell_length(1)
    1132           34 :       hmat_test(2, 2) = new_cell_length(2)
    1133           34 :       hmat_test(3, 3) = new_cell_length(3)
    1134           34 :       CALL cell_create(cell_test, hmat=hmat_test(:, :), periodic=cell%perd)
    1135           34 :       CALL cp_subsys_set(oldsys, cell=cell_test)
    1136              : 
    1137              : ! now we need to scale the coordinates of all the molecules by the
    1138              : ! center of mass, using the minimum image (not all molecules are in
    1139              : ! the central box)
    1140              : 
    1141              : ! now we need to scale the coordinates of all the molecules by the
    1142              : ! center of mass
    1143           34 :       end_atom = 0
    1144          432 :       DO imolecule = 1, SUM(nchains(:, box_number))
    1145          338 :          nunits_mol = nunits(mol_type(imolecule + start_mol - 1))
    1146          338 :          start_atom = end_atom + 1
    1147          338 :          end_atom = start_atom + nunits_mol - 1
    1148              : ! now find the center of mass
    1149              :          CALL get_center_of_mass(r(:, start_atom:end_atom), nunits_mol, &
    1150          338 :                                  center_of_mass(:), mass(:, mol_type(imolecule + start_mol - 1)))
    1151              : 
    1152              : ! scale the center of mass and determine the vector that points from the
    1153              : !    old COM to the new one
    1154         1352 :          DO iside = 1, 3
    1155              :             center_of_mass_new(iside) = center_of_mass(iside)* &
    1156         1352 :                                         new_cell_length(iside)/old_cell_length(iside)
    1157              :          END DO
    1158              : 
    1159         1386 :          DO idim = 1, 3
    1160         1014 :             diff(idim) = center_of_mass_new(idim) - center_of_mass(idim)
    1161              : ! now change the particle positions
    1162         3530 :             DO iunit = start_atom, end_atom
    1163              :                particles_old%els(iunit)%r(idim) = &
    1164         3192 :                   particles_old%els(iunit)%r(idim) + diff(idim)
    1165              :             END DO
    1166              :          END DO
    1167              :       END DO
    1168              : 
    1169              : ! check for overlap
    1170              :       CALL check_for_overlap(force_env, nchains(:, box_number), &
    1171              :                              nunits(:), loverlap, mol_type(start_mol:end_mol), &
    1172           34 :                              cell_length=new_cell_length)
    1173              : 
    1174              : ! figure out if we have overlap problems
    1175           34 :       CALL group%bcast(loverlap, source)
    1176           34 :       IF (loverlap) THEN
    1177              : ! deallocate some stuff
    1178            0 :          DEALLOCATE (r)
    1179              : 
    1180            0 :          logger => cp_get_default_logger()
    1181            0 :          output_unit = cp_logger_get_default_io_unit(logger)
    1182            0 :          IF (output_unit > 0) WRITE (output_unit, *) &
    1183            0 :             "Volume move rejected due to overlap.", vol_dis
    1184              : !     end the timing
    1185            0 :          CALL timestop(handle)
    1186              : ! reset the cell and particle positions
    1187            0 :          CALL cp_subsys_set(oldsys, cell=cell_old)
    1188            0 :          DO iatom = 1, nunits_tot(box_number)
    1189            0 :             particles_old%els(iatom)%r(1:3) = r_old(1:3, iatom)
    1190              :          END DO
    1191              :          RETURN
    1192              :       END IF
    1193              : 
    1194              : ! stop if we're trying to change a box to a boxlength smaller than rcut
    1195           34 :       IF (ionode) THEN
    1196           17 :          ltoo_small = .FALSE.
    1197           17 :          IF (force_env%in_use .EQ. use_fist_force) THEN
    1198           13 :             CALL get_mc_par(mc_par, rcut=rcut)
    1199           13 :             IF (new_cell_length(1) .LT. 2.0_dp*rcut) ltoo_small = .TRUE.
    1200           13 :             IF (new_cell_length(2) .LT. 2.0_dp*rcut) ltoo_small = .TRUE.
    1201           13 :             IF (new_cell_length(3) .LT. 2.0_dp*rcut) ltoo_small = .TRUE.
    1202              : 
    1203           13 :             IF (ltoo_small) THEN
    1204            0 :                WRITE (iw, *) 'new_cell_lengths ', &
    1205            0 :                   new_cell_length(1:3)/angstrom
    1206            0 :                WRITE (iw, *) 'rcut ', rcut/angstrom
    1207              :             END IF
    1208              :          END IF
    1209              :       END IF
    1210           34 :       CALL group%bcast(ltoo_small, source)
    1211           34 :       IF (ltoo_small) &
    1212            0 :          CPABORT("Attempted a volume move where box size got too small.")
    1213              : 
    1214              : ! now compute the energy
    1215           34 :       CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
    1216              :       CALL force_env_get(force_env, &
    1217           34 :                          potential_energy=new_energy)
    1218              : 
    1219              : ! accept or reject the move
    1220              : ! to prevent overflows
    1221           34 :       energy_term = new_energy - old_energy
    1222              :       volume_term = -REAL(SUM(nchains(:, box_number)), dp)/BETA* &
    1223              :                     LOG(new_cell_length(1)*new_cell_length(2)*new_cell_length(3)/ &
    1224           94 :                         (old_cell_length(1)*old_cell_length(2)*old_cell_length(3)))
    1225           34 :       pressure_term = pressure*vol_dis
    1226              : 
    1227           34 :       value = -BETA*(energy_term + volume_term + pressure_term)
    1228           34 :       IF (value .GT. exp_max_val) THEN
    1229              :          w = 10.0_dp
    1230           34 :       ELSEIF (value .LT. exp_min_val) THEN
    1231              :          w = 0.0_dp
    1232              :       ELSE
    1233           34 :          w = EXP(value)
    1234              :       END IF
    1235              : 
    1236              : !!!!!!!!!!!!!!!! for E_V curves
    1237              : !         w=1.0E0_dp
    1238              : !         w=0.0E0_dp
    1239              : 
    1240           34 :       IF (w .GE. 1.0E0_dp) THEN
    1241           18 :          w = 1.0E0_dp
    1242           18 :          rand = 0.0E0_dp
    1243              :       ELSE
    1244           16 :          IF (ionode) rand = rng_stream%next()
    1245           16 :          CALL group%bcast(rand, source)
    1246              :       END IF
    1247              : 
    1248           34 :       IF (rand .LT. w) THEN
    1249              : 
    1250              : ! accept the move
    1251           30 :          moves%volume%successes = moves%volume%successes + 1
    1252           30 :          move_updates%volume%successes = move_updates%volume%successes + 1
    1253              : 
    1254              : ! update energies
    1255           30 :          energy_check = energy_check + (new_energy - old_energy)
    1256           30 :          old_energy = new_energy
    1257              : 
    1258          720 :          DO iatom = 1, nunits_tot(box_number)
    1259         2790 :             r_old(1:3, iatom) = particles_old%els(iatom)%r(1:3)
    1260              :          END DO
    1261              : 
    1262              : ! update discrete_array if we're doing a discrete volume move
    1263           30 :          IF (ldiscrete) THEN
    1264              :             CALL create_discrete_array(new_cell_length(:), &
    1265            0 :                                        discrete_array(:, :), discrete_step)
    1266              :          END IF
    1267              : 
    1268              :       ELSE
    1269              : 
    1270              : ! reset the cell and particle positions
    1271            4 :          CALL cp_subsys_set(oldsys, cell=cell_old)
    1272           40 :          DO iatom = 1, nunits_tot(box_number)
    1273          148 :             particles_old%els(iatom)%r(1:3) = r_old(1:3, iatom)
    1274              :          END DO
    1275              : 
    1276              :       END IF
    1277              : 
    1278              : ! deallocate some stuff
    1279           34 :       DEALLOCATE (r)
    1280           34 :       CALL cell_release(cell_test)
    1281           34 :       CALL cell_release(cell_old)
    1282              : 
    1283              : ! end the timing
    1284           34 :       CALL timestop(handle)
    1285              : 
    1286           68 :    END SUBROUTINE mc_volume_move
    1287              : 
    1288              : ! **************************************************************************************************
    1289              : !> \brief alters the length of a random bond for the given molecule, using
    1290              : !>      a mass weighted scheme so the lightest atoms move the most
    1291              : !> \param r_old the initial coordinates of all molecules in the system
    1292              : !> \param r_new the new coordinates of all molecules in the system
    1293              : !> \param mc_par the mc parameters for the force env
    1294              : !> \param molecule_type the molecule type that we're moving
    1295              : !> \param molecule_kind the structure containing the molecule information
    1296              : !> \param dis_length the ratio of the new bond length to the old bond length,
    1297              : !>                    used in the acceptance rule
    1298              : !> \param particles the particle_list_type for all particles in the force_env..
    1299              : !>             used to grab the mass of each atom
    1300              : !> \param rng_stream the random number stream that we draw from
    1301              : !>
    1302              : !>    This subroutine is written to be parallel.
    1303              : !> \author MJM
    1304              : ! **************************************************************************************************
    1305          312 :    SUBROUTINE change_bond_length(r_old, r_new, mc_par, molecule_type, molecule_kind, &
    1306              :                                  dis_length, particles, rng_stream)
    1307              : 
    1308              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_old
    1309              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: r_new
    1310              :       TYPE(mc_simpar_type), POINTER                      :: mc_par
    1311              :       INTEGER, INTENT(IN)                                :: molecule_type
    1312              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
    1313              :       REAL(KIND=dp), INTENT(OUT)                         :: dis_length
    1314              :       TYPE(particle_list_type), POINTER                  :: particles
    1315              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
    1316              : 
    1317              :       CHARACTER(len=*), PARAMETER :: routineN = 'change_bond_length'
    1318              : 
    1319              :       INTEGER                                            :: bond_number, handle, i, iatom, ibond, &
    1320              :                                                             ipart, natom, nbond, source
    1321          312 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_a, atom_b, counter
    1322          312 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: connection, connectivity
    1323          312 :       INTEGER, DIMENSION(:), POINTER                     :: nunits
    1324              :       LOGICAL                                            :: ionode
    1325          312 :       REAL(dp), DIMENSION(:), POINTER                    :: rmbond
    1326              :       REAL(KIND=dp)                                      :: atom_mass, mass_a, mass_b, new_length_a, &
    1327              :                                                             new_length_b, old_length, rand
    1328              :       REAL(KIND=dp), DIMENSION(1:3)                      :: bond_a, bond_b
    1329          312 :       TYPE(bond_type), DIMENSION(:), POINTER             :: bond_list
    1330              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
    1331              :       TYPE(mp_comm_type)                                 :: group
    1332              : 
    1333              : ! begin the timing of the subroutine
    1334              : 
    1335          312 :       CALL timeset(routineN, handle)
    1336              : 
    1337          312 :       NULLIFY (rmbond, mc_molecule_info)
    1338              : 
    1339              : ! get some stuff from mc_par
    1340              :       CALL get_mc_par(mc_par, mc_molecule_info=mc_molecule_info, source=source, &
    1341          312 :                       group=group, rmbond=rmbond, ionode=ionode)
    1342          312 :       CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
    1343              : 
    1344              : ! copy the incoming coordinates so we can change them
    1345         1248 :       DO ipart = 1, nunits(molecule_type)
    1346         4056 :          r_new(1:3, ipart) = r_old(1:3, ipart)
    1347              :       END DO
    1348              : 
    1349              : ! pick which bond in the molecule at random
    1350          312 :       IF (ionode) THEN
    1351          156 :          rand = rng_stream%next()
    1352              :       END IF
    1353          312 :       CALL group%bcast(rand, source)
    1354              :       CALL get_molecule_kind(molecule_kind, natom=natom, nbond=nbond, &
    1355          312 :                              bond_list=bond_list)
    1356          312 :       bond_number = CEILING(rand*REAL(nbond, dp))
    1357              : 
    1358          936 :       ALLOCATE (connection(1:natom, 1:2))
    1359              : ! assume at most six bonds per atom
    1360          936 :       ALLOCATE (connectivity(1:6, 1:natom))
    1361          936 :       ALLOCATE (counter(1:natom))
    1362          624 :       ALLOCATE (atom_a(1:natom))
    1363          624 :       ALLOCATE (atom_b(1:natom))
    1364         2808 :       connection(:, :) = 0
    1365         6864 :       connectivity(:, :) = 0
    1366         1248 :       counter(:) = 0
    1367         1248 :       atom_a(:) = 0
    1368         1248 :       atom_b(:) = 0
    1369              : 
    1370              : ! now we need to find a list of atoms that each atom in this bond is connected
    1371              : ! to
    1372         1248 :       DO iatom = 1, natom
    1373         3120 :          DO ibond = 1, nbond
    1374         2808 :             IF (bond_list(ibond)%a == iatom) THEN
    1375          624 :                counter(iatom) = counter(iatom) + 1
    1376          624 :                connectivity(counter(iatom), iatom) = bond_list(ibond)%b
    1377         1248 :             ELSEIF (bond_list(ibond)%b == iatom) THEN
    1378          624 :                counter(iatom) = counter(iatom) + 1
    1379          624 :                connectivity(counter(iatom), iatom) = bond_list(ibond)%a
    1380              :             END IF
    1381              :          END DO
    1382              :       END DO
    1383              : 
    1384              : ! now I need to do a depth first search to figure out which atoms are on atom a's
    1385              : ! side and which are on atom b's
    1386         1248 :       atom_a(:) = 0
    1387          312 :       atom_a(bond_list(bond_number)%a) = 1
    1388              :       CALL depth_first_search(bond_list(bond_number)%a, bond_list(bond_number)%b, &
    1389          312 :                               connectivity(:, :), atom_a(:))
    1390         1248 :       atom_b(:) = 0
    1391          312 :       atom_b(bond_list(bond_number)%b) = 1
    1392              :       CALL depth_first_search(bond_list(bond_number)%b, bond_list(bond_number)%a, &
    1393          312 :                               connectivity(:, :), atom_b(:))
    1394              : 
    1395              : ! now figure out the masses of the various sides, so we can weight how far we move each
    1396              : ! group of atoms
    1397          312 :       mass_a = 0.0_dp
    1398          312 :       mass_b = 0.0_dp
    1399         1248 :       DO iatom = 1, natom
    1400              :          CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
    1401          936 :                               mass=atom_mass)
    1402         1248 :          IF (atom_a(iatom) == 1) THEN
    1403          624 :             mass_a = mass_a + atom_mass
    1404              :          ELSE
    1405          312 :             mass_b = mass_b + atom_mass
    1406              :          END IF
    1407              :       END DO
    1408              : 
    1409              : ! choose a displacement
    1410          312 :       IF (ionode) rand = rng_stream%next()
    1411          312 :       CALL group%bcast(rand, source)
    1412              : 
    1413          312 :       dis_length = rmbond(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
    1414              : 
    1415              : ! find the bond vector that atom a will be moving
    1416         1248 :       DO i = 1, 3
    1417              :          bond_a(i) = r_new(i, bond_list(bond_number)%a) - &
    1418          936 :                      r_new(i, bond_list(bond_number)%b)
    1419         1248 :          bond_b(i) = -bond_a(i)
    1420              :       END DO
    1421              : 
    1422              : ! notice we weight by the opposite masses...therefore lighter segments
    1423              : ! will move further
    1424         1248 :       old_length = SQRT(DOT_PRODUCT(bond_a, bond_a))
    1425          312 :       new_length_a = dis_length*mass_b/(mass_a + mass_b)
    1426          312 :       new_length_b = dis_length*mass_a/(mass_a + mass_b)
    1427              : 
    1428         1248 :       DO i = 1, 3
    1429          936 :          bond_a(i) = bond_a(i)/old_length*new_length_a
    1430         1248 :          bond_b(i) = bond_b(i)/old_length*new_length_b
    1431              :       END DO
    1432              : 
    1433         1248 :       DO iatom = 1, natom
    1434         1248 :          IF (atom_a(iatom) == 1) THEN
    1435          624 :             r_new(1, iatom) = r_new(1, iatom) + bond_a(1)
    1436          624 :             r_new(2, iatom) = r_new(2, iatom) + bond_a(2)
    1437          624 :             r_new(3, iatom) = r_new(3, iatom) + bond_a(3)
    1438              :          ELSE
    1439          312 :             r_new(1, iatom) = r_new(1, iatom) + bond_b(1)
    1440          312 :             r_new(2, iatom) = r_new(2, iatom) + bond_b(2)
    1441          312 :             r_new(3, iatom) = r_new(3, iatom) + bond_b(3)
    1442              :          END IF
    1443              :       END DO
    1444              : 
    1445              : ! correct the value of dis_length for the acceptance rule
    1446          312 :       dis_length = (old_length + dis_length)/old_length
    1447              : 
    1448          312 :       DEALLOCATE (connection)
    1449          312 :       DEALLOCATE (connectivity)
    1450          312 :       DEALLOCATE (counter)
    1451          312 :       DEALLOCATE (atom_a)
    1452          312 :       DEALLOCATE (atom_b)
    1453              : ! end the timing
    1454          312 :       CALL timestop(handle)
    1455              : 
    1456          624 :    END SUBROUTINE change_bond_length
    1457              : 
    1458              : ! **************************************************************************************************
    1459              : !> \brief Alters the magnitude of a random angle in a molecule centered on atom C
    1460              : !>      (connected to atoms A and B).  Atoms A and B are moved amounts related
    1461              : !>      to their masses (and masses of all connecting atoms), so that heavier
    1462              : !>      segments are moved less.
    1463              : !> \param r_old the initial coordinates of all molecules in the system
    1464              : !> \param r_new the new coordinates of all molecules in the system
    1465              : !> \param mc_par the mc parameters for the force env
    1466              : !> \param molecule_type the type of molecule we're playing with
    1467              : !> \param molecule_kind the structure containing the molecule information
    1468              : !> \param particles the particle_list_type for all particles in the force_env...
    1469              : !>             used to grab the mass of each atom
    1470              : !> \param rng_stream the random number stream that we draw from
    1471              : !> \author MJM
    1472              : ! **************************************************************************************************
    1473          194 :    SUBROUTINE change_bond_angle(r_old, r_new, mc_par, molecule_type, molecule_kind, &
    1474              :                                 particles, rng_stream)
    1475              : 
    1476              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_old
    1477              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: r_new
    1478              :       TYPE(mc_simpar_type), POINTER                      :: mc_par
    1479              :       INTEGER, INTENT(IN)                                :: molecule_type
    1480              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
    1481              :       TYPE(particle_list_type), POINTER                  :: particles
    1482              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
    1483              : 
    1484              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'change_bond_angle'
    1485              : 
    1486              :       INTEGER                                            :: bend_number, handle, i, iatom, ibond, &
    1487              :                                                             ipart, natom, nbend, nbond, source
    1488          194 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_a, atom_c, counter
    1489          194 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: connection, connectivity
    1490          194 :       INTEGER, DIMENSION(:), POINTER                     :: nunits
    1491              :       LOGICAL                                            :: ionode
    1492          194 :       REAL(dp), DIMENSION(:), POINTER                    :: rmangle
    1493              :       REAL(KIND=dp) :: atom_mass, bis_length, dis_angle, dis_angle_a, dis_angle_c, mass_a, mass_c, &
    1494              :          new_angle_a, new_angle_c, old_angle, old_length_a, old_length_c, rand, temp_length
    1495              :       REAL(KIND=dp), DIMENSION(1:3)                      :: bisector, bond_a, bond_c, cross_prod, &
    1496              :                                                             cross_prod_plane, temp
    1497          194 :       TYPE(bend_type), DIMENSION(:), POINTER             :: bend_list
    1498          194 :       TYPE(bond_type), DIMENSION(:), POINTER             :: bond_list
    1499              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
    1500              :       TYPE(mp_comm_type)                                 :: group
    1501              : 
    1502              : ! begin the timing of the subroutine
    1503              : 
    1504          194 :       CALL timeset(routineN, handle)
    1505              : 
    1506          194 :       NULLIFY (bend_list, bond_list, rmangle, mc_molecule_info)
    1507              : 
    1508              : ! get some stuff from mc_par
    1509              :       CALL get_mc_par(mc_par, rmangle=rmangle, source=source, &
    1510          194 :                       group=group, ionode=ionode, mc_molecule_info=mc_molecule_info)
    1511          194 :       CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
    1512              : 
    1513              : ! copy the incoming coordinates so we can change them
    1514          776 :       DO ipart = 1, nunits(molecule_type)
    1515         2522 :          r_new(1:3, ipart) = r_old(1:3, ipart)
    1516              :       END DO
    1517              : 
    1518              : ! pick which bond in the molecule at random
    1519          194 :       IF (ionode) THEN
    1520           97 :          rand = rng_stream%next()
    1521              :       END IF
    1522          194 :       CALL group%bcast(rand, source)
    1523              :       CALL get_molecule_kind(molecule_kind, natom=natom, nbend=nbend, &
    1524          194 :                              bend_list=bend_list, bond_list=bond_list, nbond=nbond)
    1525          194 :       bend_number = CEILING(rand*REAL(nbend, dp))
    1526              : 
    1527          582 :       ALLOCATE (connection(1:natom, 1:2))
    1528              : ! assume at most six bonds per atom
    1529          582 :       ALLOCATE (connectivity(1:6, 1:natom))
    1530          582 :       ALLOCATE (counter(1:natom))
    1531          388 :       ALLOCATE (atom_a(1:natom))
    1532          388 :       ALLOCATE (atom_c(1:natom))
    1533         1746 :       connection(:, :) = 0
    1534         4268 :       connectivity(:, :) = 0
    1535          776 :       counter(:) = 0
    1536          776 :       atom_a(:) = 0
    1537          776 :       atom_c(:) = 0
    1538              : 
    1539              : ! now we need to find a list of atoms that each atom in this bond is connected
    1540              : ! to
    1541          776 :       DO iatom = 1, natom
    1542         1940 :          DO ibond = 1, nbond
    1543         1746 :             IF (bond_list(ibond)%a == iatom) THEN
    1544          388 :                counter(iatom) = counter(iatom) + 1
    1545          388 :                connectivity(counter(iatom), iatom) = bond_list(ibond)%b
    1546          776 :             ELSEIF (bond_list(ibond)%b == iatom) THEN
    1547          388 :                counter(iatom) = counter(iatom) + 1
    1548          388 :                connectivity(counter(iatom), iatom) = bond_list(ibond)%a
    1549              :             END IF
    1550              :          END DO
    1551              :       END DO
    1552              : 
    1553              : ! now I need to do a depth first search to figure out which atoms are on atom a's
    1554              : ! side and which are on atom c's
    1555          776 :       atom_a(:) = 0
    1556          194 :       atom_a(bend_list(bend_number)%a) = 1
    1557              :       CALL depth_first_search(bend_list(bend_number)%a, bend_list(bend_number)%b, &
    1558          194 :                               connectivity(:, :), atom_a(:))
    1559          776 :       atom_c(:) = 0
    1560          194 :       atom_c(bend_list(bend_number)%c) = 1
    1561              :       CALL depth_first_search(bend_list(bend_number)%c, bend_list(bend_number)%b, &
    1562          194 :                               connectivity(:, :), atom_c(:))
    1563              : 
    1564              : ! now figure out the masses of the various sides, so we can weight how far we move each
    1565              : ! group of atoms
    1566          194 :       mass_a = 0.0_dp
    1567          194 :       mass_c = 0.0_dp
    1568          776 :       DO iatom = 1, natom
    1569              :          CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
    1570          582 :                               mass=atom_mass)
    1571          582 :          IF (atom_a(iatom) == 1) mass_a = mass_a + atom_mass
    1572         1358 :          IF (atom_c(iatom) == 1) mass_c = mass_c + atom_mass
    1573              :       END DO
    1574              : 
    1575              : ! choose a displacement
    1576          194 :       IF (ionode) rand = rng_stream%next()
    1577          194 :       CALL group%bcast(rand, source)
    1578              : 
    1579          194 :       dis_angle = rmangle(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
    1580              : 
    1581              : ! need to find the A-B-C bisector
    1582              : 
    1583              : ! this going to be tough...we need to find the plane of the A-B-C bond and only shift
    1584              : ! that component for all atoms connected to A and C...otherwise we change other
    1585              : ! internal degrees of freedom
    1586              : 
    1587              : ! find the bond vectors
    1588          776 :       DO i = 1, 3
    1589              :          bond_a(i) = r_new(i, bend_list(bend_number)%a) - &
    1590          582 :                      r_new(i, bend_list(bend_number)%b)
    1591              :          bond_c(i) = r_new(i, bend_list(bend_number)%c) - &
    1592          776 :                      r_new(i, bend_list(bend_number)%b)
    1593              :       END DO
    1594          776 :       old_length_a = SQRT(DOT_PRODUCT(bond_a, bond_a))
    1595          776 :       old_length_c = SQRT(DOT_PRODUCT(bond_c, bond_c))
    1596          776 :       old_angle = ACOS(DOT_PRODUCT(bond_a, bond_c)/(old_length_a*old_length_c))
    1597              : 
    1598          776 :       DO i = 1, 3
    1599              :          bisector(i) = bond_a(i)/old_length_a + & ! not yet normalized
    1600          776 :                        bond_c(i)/old_length_c
    1601              :       END DO
    1602          776 :       bis_length = SQRT(DOT_PRODUCT(bisector, bisector))
    1603          776 :       bisector(1:3) = bisector(1:3)/bis_length
    1604              : 
    1605              : ! now we need to find the cross product of the B-A and B-C vectors and normalize
    1606              : ! it, so we have a vector that defines the bend plane
    1607          194 :       cross_prod(1) = bond_a(2)*bond_c(3) - bond_a(3)*bond_c(2)
    1608          194 :       cross_prod(2) = bond_a(3)*bond_c(1) - bond_a(1)*bond_c(3)
    1609          194 :       cross_prod(3) = bond_a(1)*bond_c(2) - bond_a(2)*bond_c(1)
    1610         1358 :       cross_prod(1:3) = cross_prod(1:3)/SQRT(DOT_PRODUCT(cross_prod, cross_prod))
    1611              : 
    1612              : ! we have two axis of a coordinate system...let's get the third
    1613          194 :       cross_prod_plane(1) = cross_prod(2)*bisector(3) - cross_prod(3)*bisector(2)
    1614          194 :       cross_prod_plane(2) = cross_prod(3)*bisector(1) - cross_prod(1)*bisector(3)
    1615          194 :       cross_prod_plane(3) = cross_prod(1)*bisector(2) - cross_prod(2)*bisector(1)
    1616              :       cross_prod_plane(1:3) = cross_prod_plane(1:3)/ &
    1617         1358 :                               SQRT(DOT_PRODUCT(cross_prod_plane, cross_prod_plane))
    1618              : 
    1619              : ! now bisector is x, cross_prod_plane is the y vector (pointing towards c),
    1620              : ! and cross_prod is z
    1621              : ! shift the molecule so that atom b is at the origin
    1622          776 :       DO iatom = 1, natom
    1623              :          r_new(1:3, iatom) = r_new(1:3, iatom) - &
    1624         2522 :                              r_old(1:3, bend_list(bend_number)%b)
    1625              :       END DO
    1626              : 
    1627              : ! figure out how much we move each side, since we're mass-weighting, by the
    1628              : ! opposite masses, so lighter moves farther..this angle is the angle between
    1629              : ! the bond vector BA or BC and the bisector
    1630          194 :       dis_angle_a = dis_angle*mass_c/(mass_a + mass_c)
    1631          194 :       dis_angle_c = dis_angle*mass_a/(mass_a + mass_c)
    1632              : 
    1633              : ! now loop through all the atoms, moving the ones that are connected to a or c
    1634          776 :       DO iatom = 1, natom
    1635              : ! subtract out the z component (perpendicular to the angle plane)
    1636              :          temp(1:3) = r_new(1:3, iatom) - &
    1637              :                      DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
    1638         4074 :                      cross_prod(1:3)
    1639         2328 :          temp_length = SQRT(DOT_PRODUCT(temp, temp))
    1640              : 
    1641              : ! we can now compute all three components of the new bond vector along the
    1642              : ! axis defined above
    1643          776 :          IF (atom_a(iatom) == 1) THEN
    1644              : 
    1645              : ! if the y-coordinate is less than zero, we need to switch the sign when we make the vector,
    1646              : ! as the angle computed by the dot product can't distinguish between that
    1647          776 :             IF (DOT_PRODUCT(cross_prod_plane(1:3), r_new(1:3, iatom)) &
    1648              :                 .LT. 0.0_dp) THEN
    1649              : 
    1650              : ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is
    1651              :                new_angle_a = ACOS(DOT_PRODUCT(bisector, temp(1:3))/ &
    1652          776 :                                   (temp_length)) + dis_angle_a
    1653              : 
    1654              :                r_new(1:3, iatom) = COS(new_angle_a)*temp_length*bisector(1:3) - &
    1655              :                                    SIN(new_angle_a)*temp_length*cross_prod_plane(1:3) + &
    1656              :                                    DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
    1657         1358 :                                    cross_prod(1:3)
    1658              :             ELSE
    1659              : 
    1660              : ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is
    1661              :                new_angle_a = ACOS(DOT_PRODUCT(bisector, temp(1:3))/ &
    1662            0 :                                   (temp_length)) - dis_angle_a
    1663              : 
    1664              :                r_new(1:3, iatom) = COS(new_angle_a)*temp_length*bisector(1:3) + &
    1665              :                                    SIN(new_angle_a)*temp_length*cross_prod_plane(1:3) + &
    1666              :                                    DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
    1667            0 :                                    cross_prod(1:3)
    1668              :             END IF
    1669              : 
    1670          388 :          ELSEIF (atom_c(iatom) == 1) THEN
    1671              : 
    1672              : ! if the y-coordinate is less than zero, we need to switch the sign when we make the vector,
    1673              : ! as the angle computed by the dot product can't distinguish between that
    1674          776 :             IF (DOT_PRODUCT(cross_prod_plane(1:3), r_new(1:3, iatom)) &
    1675              :                 .LT. 0.0_dp) THEN
    1676              : ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is
    1677              :                new_angle_c = ACOS(DOT_PRODUCT(bisector(1:3), temp(1:3))/ &
    1678            0 :                                   (temp_length)) - dis_angle_c
    1679              : 
    1680              :                r_new(1:3, iatom) = COS(new_angle_c)*temp_length*bisector(1:3) - &
    1681              :                                    SIN(new_angle_c)*temp_length*cross_prod_plane(1:3) + &
    1682              :                                    DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
    1683            0 :                                    cross_prod(1:3)
    1684              :             ELSE
    1685              :                new_angle_c = ACOS(DOT_PRODUCT(bisector(1:3), temp(1:3))/ &
    1686          776 :                                   (temp_length)) + dis_angle_c
    1687              : 
    1688              :                r_new(1:3, iatom) = COS(new_angle_c)*temp_length*bisector(1:3) + &
    1689              :                                    SIN(new_angle_c)*temp_length*cross_prod_plane(1:3) + &
    1690              :                                    DOT_PRODUCT(cross_prod(1:3), r_new(1:3, iatom))* &
    1691         1358 :                                    cross_prod(1:3)
    1692              :             END IF
    1693              :          END IF
    1694              : 
    1695              :       END DO
    1696              : 
    1697          776 :       DO iatom = 1, natom
    1698              :          r_new(1:3, iatom) = r_new(1:3, iatom) + &
    1699         2522 :                              r_old(1:3, bend_list(bend_number)%b)
    1700              :       END DO
    1701              : 
    1702              : ! deallocate some stuff
    1703          194 :       DEALLOCATE (connection)
    1704          194 :       DEALLOCATE (connectivity)
    1705          194 :       DEALLOCATE (counter)
    1706          194 :       DEALLOCATE (atom_a)
    1707          194 :       DEALLOCATE (atom_c)
    1708              : 
    1709              : ! end the timing
    1710          194 :       CALL timestop(handle)
    1711              : 
    1712          388 :    END SUBROUTINE change_bond_angle
    1713              : 
    1714              : ! **************************************************************************************************
    1715              : !> \brief Alters a dihedral (A-B-C-D) in the molecule so that all other internal
    1716              : !>      degrees of freedom remain the same.  If other dihedrals are centered
    1717              : !>      on B-C, they rotate as well to keep the relationship between the
    1718              : !>      dihedrals the same.  Atoms A and D are moved amounts related to their
    1719              : !>      masses (and masses of all connecting atoms), so that heavier segments
    1720              : !>      are moved less.  All atoms except B and C are rotated around the
    1721              : !>      B-C bond vector (B and C are not moved).
    1722              : !> \param r_old the initial coordinates of all molecules in the system
    1723              : !> \param r_new the new coordinates of all molecules in the system
    1724              : !> \param mc_par the mc parameters for the force env
    1725              : !> \param molecule_type the type of molecule we're playing with
    1726              : !> \param molecule_kind the structure containing the molecule information
    1727              : !> \param particles the particle_list_type for all particles in the force_env..
    1728              : !>             used to grab the mass of each atom
    1729              : !> \param rng_stream the random number stream that we draw from
    1730              : !> \author MJM
    1731              : ! **************************************************************************************************
    1732            0 :    SUBROUTINE change_dihedral(r_old, r_new, mc_par, molecule_type, molecule_kind, &
    1733              :                               particles, rng_stream)
    1734              : 
    1735              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_old
    1736              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: r_new
    1737              :       TYPE(mc_simpar_type), POINTER                      :: mc_par
    1738              :       INTEGER, INTENT(IN)                                :: molecule_type
    1739              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
    1740              :       TYPE(particle_list_type), POINTER                  :: particles
    1741              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
    1742              : 
    1743              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'change_dihedral'
    1744              : 
    1745              :       INTEGER                                            :: handle, i, iatom, ibond, ipart, natom, &
    1746              :                                                             nbond, ntorsion, source, torsion_number
    1747            0 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_a, atom_d, counter
    1748            0 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: connection, connectivity
    1749            0 :       INTEGER, DIMENSION(:), POINTER                     :: nunits
    1750              :       LOGICAL                                            :: ionode
    1751            0 :       REAL(dp), DIMENSION(:), POINTER                    :: rmdihedral
    1752              :       REAL(KIND=dp)                                      :: atom_mass, dis_angle, dis_angle_a, &
    1753              :                                                             dis_angle_d, mass_a, mass_d, &
    1754              :                                                             old_length_a, rand, u, v, w, x, y, z
    1755              :       REAL(KIND=dp), DIMENSION(1:3)                      :: bond_a, temp
    1756            0 :       TYPE(bond_type), DIMENSION(:), POINTER             :: bond_list
    1757              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
    1758              :       TYPE(mp_comm_type)                                 :: group
    1759            0 :       TYPE(torsion_type), DIMENSION(:), POINTER          :: torsion_list
    1760              : 
    1761              : ! begin the timing of the subroutine
    1762              : 
    1763            0 :       CALL timeset(routineN, handle)
    1764              : 
    1765            0 :       NULLIFY (rmdihedral, torsion_list, bond_list, mc_molecule_info)
    1766              : 
    1767              : ! get some stuff from mc_par
    1768              :       CALL get_mc_par(mc_par, rmdihedral=rmdihedral, &
    1769              :                       source=source, group=group, ionode=ionode, &
    1770            0 :                       mc_molecule_info=mc_molecule_info)
    1771            0 :       CALL get_mc_molecule_info(mc_molecule_info, nunits=nunits)
    1772              : 
    1773              : ! copy the incoming coordinates so we can change them
    1774            0 :       DO ipart = 1, nunits(molecule_type)
    1775            0 :          r_new(1:3, ipart) = r_old(1:3, ipart)
    1776              :       END DO
    1777              : 
    1778              : ! pick which bond in the molecule at random
    1779            0 :       IF (ionode) THEN
    1780            0 :          rand = rng_stream%next()
    1781              : !      CALL RANDOM_NUMBER(rand)
    1782              :       END IF
    1783            0 :       CALL group%bcast(rand, source)
    1784              :       CALL get_molecule_kind(molecule_kind, natom=natom, &
    1785              :                              bond_list=bond_list, nbond=nbond, &
    1786            0 :                              ntorsion=ntorsion, torsion_list=torsion_list)
    1787            0 :       torsion_number = CEILING(rand*REAL(ntorsion, dp))
    1788              : 
    1789            0 :       ALLOCATE (connection(1:natom, 1:2))
    1790              : ! assume at most six bonds per atom
    1791            0 :       ALLOCATE (connectivity(1:6, 1:natom))
    1792            0 :       ALLOCATE (counter(1:natom))
    1793            0 :       ALLOCATE (atom_a(1:natom))
    1794            0 :       ALLOCATE (atom_d(1:natom))
    1795            0 :       connection(:, :) = 0
    1796            0 :       connectivity(:, :) = 0
    1797            0 :       counter(:) = 0
    1798            0 :       atom_a(:) = 0
    1799            0 :       atom_d(:) = 0
    1800              : 
    1801              : ! now we need to find a list of atoms that each atom in this bond is connected
    1802              : ! to
    1803            0 :       DO iatom = 1, natom
    1804            0 :          DO ibond = 1, nbond
    1805            0 :             IF (bond_list(ibond)%a == iatom) THEN
    1806            0 :                counter(iatom) = counter(iatom) + 1
    1807            0 :                connectivity(counter(iatom), iatom) = bond_list(ibond)%b
    1808            0 :             ELSEIF (bond_list(ibond)%b == iatom) THEN
    1809            0 :                counter(iatom) = counter(iatom) + 1
    1810            0 :                connectivity(counter(iatom), iatom) = bond_list(ibond)%a
    1811              :             END IF
    1812              :          END DO
    1813              :       END DO
    1814              : 
    1815              : ! now I need to do a depth first search to figure out which atoms are on atom
    1816              : ! a's side and which are on atom d's, but remember we're moving all atoms on a's
    1817              : ! side of b, including atoms not in a's branch
    1818            0 :       atom_a(:) = 0
    1819            0 :       atom_a(torsion_list(torsion_number)%a) = 1
    1820              :       CALL depth_first_search(torsion_list(torsion_number)%b, &
    1821            0 :                               torsion_list(torsion_number)%c, connectivity(:, :), atom_a(:))
    1822            0 :       atom_d(:) = 0
    1823            0 :       atom_d(torsion_list(torsion_number)%d) = 1
    1824              :       CALL depth_first_search(torsion_list(torsion_number)%c, &
    1825            0 :                               torsion_list(torsion_number)%b, connectivity(:, :), atom_d(:))
    1826              : 
    1827              : ! now figure out the masses of the various sides, so we can weight how far we
    1828              : ! move each group of atoms
    1829            0 :       mass_a = 0.0_dp
    1830            0 :       mass_d = 0.0_dp
    1831            0 :       DO iatom = 1, natom
    1832              :          CALL get_atomic_kind(particles%els(iatom)%atomic_kind, &
    1833            0 :                               mass=atom_mass)
    1834            0 :          IF (atom_a(iatom) == 1) mass_a = mass_a + atom_mass
    1835            0 :          IF (atom_d(iatom) == 1) mass_d = mass_d + atom_mass
    1836              :       END DO
    1837              : 
    1838              : ! choose a displacement
    1839            0 :       IF (ionode) rand = rng_stream%next()
    1840            0 :       CALL group%bcast(rand, source)
    1841              : 
    1842            0 :       dis_angle = rmdihedral(molecule_type)*2.0E0_dp*(rand - 0.5E0_dp)
    1843              : 
    1844              : ! find the bond vectors, B-C, so we know what to rotate around
    1845            0 :       DO i = 1, 3
    1846              :          bond_a(i) = r_new(i, torsion_list(torsion_number)%c) - &
    1847            0 :                      r_new(i, torsion_list(torsion_number)%b)
    1848              :       END DO
    1849            0 :       old_length_a = SQRT(DOT_PRODUCT(bond_a, bond_a))
    1850            0 :       bond_a(1:3) = bond_a(1:3)/old_length_a
    1851              : 
    1852              : ! figure out how much we move each side, since we're mass-weighting, by the
    1853              : ! opposite masses, so lighter moves farther...we take the opposite sign of d
    1854              : ! so we're not rotating both angles in the same direction
    1855            0 :       dis_angle_a = dis_angle*mass_d/(mass_a + mass_d)
    1856            0 :       dis_angle_d = -dis_angle*mass_a/(mass_a + mass_d)
    1857              : 
    1858            0 :       DO iatom = 1, natom
    1859              : 
    1860            0 :          IF (atom_a(iatom) == 1) THEN
    1861              : ! shift the coords so b is at the origin
    1862              :             r_new(1:3, iatom) = r_new(1:3, iatom) - &
    1863            0 :                                 r_new(1:3, torsion_list(torsion_number)%b)
    1864              : 
    1865              : ! multiply by the rotation matrix
    1866            0 :             u = bond_a(1)
    1867            0 :             v = bond_a(2)
    1868            0 :             w = bond_a(3)
    1869            0 :             x = r_new(1, iatom)
    1870            0 :             y = r_new(2, iatom)
    1871            0 :             z = r_new(3, iatom)
    1872              :             temp(1) = (u*(u*x + v*y + w*z) + (x*(v**2 + w**2) - u*(v*y + w*z))*COS(dis_angle_a) + &
    1873            0 :                        SQRT(u**2 + v**2 + w**2)*(v*z - w*y)*SIN(dis_angle_a))/(u**2 + v**2 + w**2)
    1874              :             temp(2) = (v*(u*x + v*y + w*z) + (y*(u**2 + w**2) - v*(u*x + w*z))*COS(dis_angle_a) + &
    1875            0 :                        SQRT(u**2 + v**2 + w**2)*(w*x - u*z)*SIN(dis_angle_a))/(u**2 + v**2 + w**2)
    1876              :             temp(3) = (w*(u*x + v*y + w*z) + (z*(v**2 + u**2) - w*(u*x + v*y))*COS(dis_angle_a) + &
    1877            0 :                        SQRT(u**2 + v**2 + w**2)*(u*y - v*x)*SIN(dis_angle_a))/(u**2 + v**2 + w**2)
    1878              : 
    1879              : ! shift back to the original position
    1880            0 :             temp(1:3) = temp(1:3) + r_new(1:3, torsion_list(torsion_number)%b)
    1881            0 :             r_new(1:3, iatom) = temp(1:3)
    1882              : 
    1883            0 :          ELSEIF (atom_d(iatom) == 1) THEN
    1884              : 
    1885              : ! shift the coords so c is at the origin
    1886              :             r_new(1:3, iatom) = r_new(1:3, iatom) - &
    1887            0 :                                 r_new(1:3, torsion_list(torsion_number)%c)
    1888              : 
    1889              : ! multiply by the rotation matrix
    1890            0 :             u = bond_a(1)
    1891            0 :             v = bond_a(2)
    1892            0 :             w = bond_a(3)
    1893            0 :             x = r_new(1, iatom)
    1894            0 :             y = r_new(2, iatom)
    1895            0 :             z = r_new(3, iatom)
    1896              :             temp(1) = (u*(u*x + v*y + w*z) + (x*(v**2 + w**2) - u*(v*y + w*z))*COS(dis_angle_d) + &
    1897            0 :                        SQRT(u**2 + v**2 + w**2)*(v*z - w*y)*SIN(dis_angle_d))/(u**2 + v**2 + w**2)
    1898              :             temp(2) = (v*(u*x + v*y + w*z) + (y*(u**2 + w**2) - v*(u*x + w*z))*COS(dis_angle_d) + &
    1899            0 :                        SQRT(u**2 + v**2 + w**2)*(w*x - u*z)*SIN(dis_angle_d))/(u**2 + v**2 + w**2)
    1900              :             temp(3) = (w*(u*x + v*y + w*z) + (z*(v**2 + u**2) - w*(u*x + v*y))*COS(dis_angle_d) + &
    1901            0 :                        SQRT(u**2 + v**2 + w**2)*(u*y - v*x)*SIN(dis_angle_d))/(u**2 + v**2 + w**2)
    1902              : 
    1903              : ! shift back to the original position
    1904            0 :             temp(1:3) = temp(1:3) + r_new(1:3, torsion_list(torsion_number)%c)
    1905            0 :             r_new(1:3, iatom) = temp(1:3)
    1906              :          END IF
    1907              :       END DO
    1908              : 
    1909              : ! deallocate some stuff
    1910            0 :       DEALLOCATE (connection)
    1911            0 :       DEALLOCATE (connectivity)
    1912            0 :       DEALLOCATE (counter)
    1913            0 :       DEALLOCATE (atom_a)
    1914            0 :       DEALLOCATE (atom_d)
    1915              : 
    1916              : ! end the timing
    1917            0 :       CALL timestop(handle)
    1918              : 
    1919            0 :    END SUBROUTINE change_dihedral
    1920              : 
    1921              : ! **************************************************************************************************
    1922              : !> \brief performs either a bond or angle change move for a given molecule
    1923              : !> \param mc_par the mc parameters for the force env
    1924              : !> \param force_env the force environment used in the move
    1925              : !> \param bias_env the force environment used to bias the move, if any (it may
    1926              : !>            be null if lbias=.false. in mc_par)
    1927              : !> \param moves the structure that keeps track of how many moves have been
    1928              : !>               accepted/rejected
    1929              : !> \param energy_check the running energy difference between now and the initial
    1930              : !>        energy
    1931              : !> \param r_old the coordinates of force_env before the move
    1932              : !> \param old_energy the energy of the force_env before the move
    1933              : !> \param start_atom_swap the number of the swap molecule's first atom, assuming the rest of
    1934              : !>        the atoms follow sequentially
    1935              : !> \param target_atom the number of the target atom for swapping
    1936              : !> \param molecule_type the molecule type for the atom we're swapping
    1937              : !> \param box_number the number of the box we're doing this move in
    1938              : !> \param bias_energy_old the biased energy of the system before the move
    1939              : !> \param last_bias_energy the last biased energy of the system
    1940              : !> \param move_type dictates if we're moving to an "in" or "out" region
    1941              : !> \param rng_stream the random number stream that we draw from
    1942              : !> \author MJM
    1943              : !> \note     Designed for parallel.
    1944              : ! **************************************************************************************************
    1945            0 :    SUBROUTINE mc_avbmc_move(mc_par, force_env, bias_env, moves, &
    1946            0 :                             energy_check, r_old, old_energy, start_atom_swap, &
    1947              :                             target_atom, &
    1948              :                             molecule_type, box_number, bias_energy_old, last_bias_energy, &
    1949              :                             move_type, rng_stream)
    1950              : 
    1951              :       TYPE(mc_simpar_type), POINTER                      :: mc_par
    1952              :       TYPE(force_env_type), POINTER                      :: force_env, bias_env
    1953              :       TYPE(mc_moves_type), POINTER                       :: moves
    1954              :       REAL(KIND=dp), INTENT(INOUT)                       :: energy_check
    1955              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: r_old
    1956              :       REAL(KIND=dp), INTENT(INOUT)                       :: old_energy
    1957              :       INTEGER, INTENT(IN)                                :: start_atom_swap, target_atom, &
    1958              :                                                             molecule_type, box_number
    1959              :       REAL(KIND=dp), INTENT(INOUT)                       :: bias_energy_old, last_bias_energy
    1960              :       CHARACTER(LEN=*), INTENT(IN)                       :: move_type
    1961              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
    1962              : 
    1963              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'mc_avbmc_move'
    1964              : 
    1965              :       INTEGER                                            :: end_mol, handle, ipart, jbox, natom, &
    1966              :                                                             nswapmoves, source, start_mol
    1967            0 :       INTEGER, DIMENSION(:), POINTER                     :: avbmc_atom, mol_type, nunits, nunits_tot
    1968            0 :       INTEGER, DIMENSION(:, :), POINTER                  :: nchains
    1969              :       LOGICAL                                            :: ionode, lbias, ldum, lin, loverlap
    1970            0 :       REAL(dp), DIMENSION(:), POINTER                    :: avbmc_rmax, avbmc_rmin, pbias
    1971            0 :       REAL(dp), DIMENSION(:, :), POINTER                 :: mass
    1972              :       REAL(KIND=dp) :: BETA, bias_energy_new, del_quickstep_energy, distance, exp_max_val, &
    1973              :          exp_min_val, max_val, min_val, new_energy, prefactor, rand, rdum, volume_in, volume_out, &
    1974              :          w, weight_new, weight_old
    1975            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: r_new
    1976              :       REAL(KIND=dp), DIMENSION(1:3)                      :: abc, RIJ
    1977              :       TYPE(cell_type), POINTER                           :: cell
    1978              :       TYPE(cp_subsys_type), POINTER                      :: subsys, subsys_force
    1979              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
    1980              :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
    1981              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
    1982              :       TYPE(mp_comm_type)                                 :: group
    1983              :       TYPE(particle_list_type), POINTER                  :: particles, particles_force
    1984              : 
    1985            0 :       rdum = 1.0_dp
    1986              : 
    1987              : ! begin the timing of the subroutine
    1988            0 :       CALL timeset(routineN, handle)
    1989              : 
    1990              : ! get a bunch of stuff from mc_par
    1991              :       CALL get_mc_par(mc_par, lbias=lbias, &
    1992              :                       BETA=BETA, max_val=max_val, min_val=min_val, exp_max_val=exp_max_val, &
    1993              :                       exp_min_val=exp_min_val, avbmc_atom=avbmc_atom, &
    1994              :                       avbmc_rmin=avbmc_rmin, avbmc_rmax=avbmc_rmax, &
    1995              :                       nswapmoves=nswapmoves, ionode=ionode, source=source, &
    1996            0 :                       group=group, pbias=pbias, mc_molecule_info=mc_molecule_info)
    1997              :       CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, &
    1998            0 :                                 mass=mass, nunits=nunits, nunits_tot=nunits_tot, mol_type=mol_type)
    1999              : ! figure out some bounds for mol_type
    2000            0 :       start_mol = 1
    2001            0 :       DO jbox = 1, box_number - 1
    2002            0 :          start_mol = start_mol + SUM(nchains(:, jbox))
    2003              :       END DO
    2004            0 :       end_mol = start_mol + SUM(nchains(:, box_number)) - 1
    2005              : 
    2006              : ! nullify some pointers
    2007            0 :       NULLIFY (particles, subsys, molecule_kinds, molecule_kind, &
    2008            0 :                particles_force, subsys_force)
    2009              : 
    2010              : ! do some allocation
    2011            0 :       ALLOCATE (r_new(1:3, 1:nunits_tot(box_number)))
    2012              : 
    2013              : ! now we need to grab and save coordinates, in case we reject
    2014              : ! are we biasing this move?
    2015            0 :       IF (lbias) THEN
    2016              : 
    2017              : ! grab the coordinates
    2018            0 :          CALL force_env_get(bias_env, cell=cell, subsys=subsys)
    2019              :          CALL cp_subsys_get(subsys, &
    2020            0 :                             particles=particles, molecule_kinds=molecule_kinds)
    2021            0 :          molecule_kind => molecule_kinds%els(1)
    2022            0 :          CALL get_molecule_kind(molecule_kind, natom=natom)
    2023            0 :          CALL get_cell(cell, abc=abc)
    2024              : 
    2025              : ! save the energy
    2026              : !         bias_energy_old=bias_energy
    2027              : 
    2028              :       ELSE
    2029              : 
    2030              : ! grab the coordinates
    2031            0 :          CALL force_env_get(force_env, cell=cell, subsys=subsys)
    2032              :          CALL cp_subsys_get(subsys, &
    2033            0 :                             particles=particles, molecule_kinds=molecule_kinds)
    2034            0 :          molecule_kind => molecule_kinds%els(1)
    2035            0 :          CALL get_molecule_kind(molecule_kind, natom=natom)
    2036            0 :          CALL get_cell(cell, abc=abc)
    2037              : 
    2038              :       END IF
    2039              : 
    2040              : ! let's determine if the molecule to be moved is in the "in" region or the
    2041              : ! "out" region of the target
    2042              :       RIJ(1) = particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(1) - &
    2043              :                particles%els(target_atom)%r(1) - abc(1)*ANINT( &
    2044              :                (particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(1) - &
    2045            0 :                 particles%els(target_atom)%r(1))/abc(1))
    2046              :       RIJ(2) = particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(2) - &
    2047              :                particles%els(target_atom)%r(2) - abc(2)*ANINT( &
    2048              :                (particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(2) - &
    2049            0 :                 particles%els(target_atom)%r(2))/abc(2))
    2050              :       RIJ(3) = particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(3) - &
    2051              :                particles%els(target_atom)%r(3) - abc(3)*ANINT( &
    2052              :                (particles%els(start_atom_swap + avbmc_atom(molecule_type) - 1)%r(3) - &
    2053            0 :                 particles%els(target_atom)%r(3))/abc(3))
    2054            0 :       distance = SQRT(RIJ(1)**2 + RIJ(2)**2 + RIJ(3)**2)
    2055            0 :       IF (distance .LE. avbmc_rmax(molecule_type) .AND. distance .GE. avbmc_rmin(molecule_type)) THEN
    2056              :          lin = .TRUE.
    2057              :       ELSE
    2058              :          lin = .FALSE.
    2059              :       END IF
    2060              : 
    2061              : ! increment the counter of the particular move we've done
    2062              : !     swapping into the "in" region of mol_target
    2063              :       IF (lin) THEN
    2064            0 :          IF (move_type == 'in') THEN
    2065              :             moves%avbmc_inin%attempts = &
    2066            0 :                moves%avbmc_inin%attempts + 1
    2067              :          ELSE
    2068              :             moves%avbmc_inout%attempts = &
    2069            0 :                moves%avbmc_inout%attempts + 1
    2070              :          END IF
    2071              :       ELSE
    2072            0 :          IF (move_type == 'in') THEN
    2073              :             moves%avbmc_outin%attempts = &
    2074            0 :                moves%avbmc_outin%attempts + 1
    2075              :          ELSE
    2076              :             moves%avbmc_outout%attempts = &
    2077            0 :                moves%avbmc_outout%attempts + 1
    2078              :          END IF
    2079              :       END IF
    2080              : 
    2081            0 :       IF (lbias) THEN
    2082              : 
    2083            0 :          IF (move_type == 'in') THEN
    2084              : 
    2085              : ! do CBMC for the old config
    2086              :             CALL generate_cbmc_swap_config(bias_env, BETA, max_val, min_val, exp_max_val, &
    2087              :                                            exp_min_val, nswapmoves, &
    2088              :                                            weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
    2089              :                                            mass(:, molecule_type), ldum, rdum, &
    2090              :                                            bias_energy_old, ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
    2091              :                                            source, group, rng_stream, &
    2092              :                                            avbmc_atom=avbmc_atom(molecule_type), &
    2093              :                                            rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='out', &
    2094            0 :                                            target_atom=target_atom)
    2095              : 
    2096              :          ELSE
    2097              : 
    2098              : ! do CBMC for the old config
    2099              :             CALL generate_cbmc_swap_config(bias_env, BETA, max_val, min_val, exp_max_val, &
    2100              :                                            exp_min_val, nswapmoves, &
    2101              :                                            weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
    2102              :                                            mass(:, molecule_type), ldum, rdum, &
    2103              :                                            bias_energy_old, ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
    2104              :                                            source, group, rng_stream, &
    2105              :                                            avbmc_atom=avbmc_atom(molecule_type), &
    2106              :                                            rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='in', &
    2107            0 :                                            target_atom=target_atom)
    2108              : 
    2109              :          END IF
    2110              : 
    2111              : ! generate the new config
    2112              :          CALL generate_cbmc_swap_config(bias_env, BETA, max_val, min_val, exp_max_val, &
    2113              :                                         exp_min_val, nswapmoves, &
    2114              :                                         weight_new, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
    2115              :                                         mass(:, molecule_type), loverlap, bias_energy_new, &
    2116              :                                         bias_energy_old, ionode, .FALSE., mol_type(start_mol:end_mol), nchains(:, box_number), &
    2117              :                                         source, group, rng_stream, &
    2118              :                                         avbmc_atom=avbmc_atom(molecule_type), &
    2119              :                                         rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type=move_type, &
    2120            0 :                                         target_atom=target_atom)
    2121              : 
    2122              : ! the energy that comes out of the above routine is the difference...we want
    2123              : ! the real energy for the acceptance rule...we don't do this for the
    2124              : ! lbias=.false. case because it doesn't appear in the acceptance rule, and
    2125              : ! we compensate in case of acceptance
    2126            0 :          bias_energy_new = bias_energy_new + bias_energy_old
    2127              : 
    2128              :       ELSE
    2129              : 
    2130            0 :          IF (move_type == 'in') THEN
    2131              : 
    2132              : ! find the weight of the old config
    2133              :             CALL generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, &
    2134              :                                            exp_min_val, nswapmoves, &
    2135              :                                            weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
    2136              :                                            mass(:, molecule_type), ldum, rdum, old_energy, &
    2137              :                                            ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
    2138              :                                            source, group, rng_stream, &
    2139              :                                            avbmc_atom=avbmc_atom(molecule_type), &
    2140              :                                            rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='out', &
    2141            0 :                                            target_atom=target_atom)
    2142              : 
    2143              :          ELSE
    2144              : 
    2145              : ! find the weight of the old config
    2146              :             CALL generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, &
    2147              :                                            exp_min_val, nswapmoves, &
    2148              :                                            weight_old, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
    2149              :                                            mass(:, molecule_type), ldum, rdum, old_energy, &
    2150              :                                            ionode, .TRUE., mol_type(start_mol:end_mol), nchains(:, box_number), &
    2151              :                                            source, group, rng_stream, &
    2152              :                                            avbmc_atom=avbmc_atom(molecule_type), &
    2153              :                                            rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type='in', &
    2154            0 :                                            target_atom=target_atom)
    2155              : 
    2156              :          END IF
    2157              : 
    2158              :          ! generate the new config...do this after, because it changes the force_env
    2159              :          CALL generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, &
    2160              :                                         exp_min_val, nswapmoves, &
    2161              :                                         weight_new, start_atom_swap, nunits_tot(box_number), nunits, nunits(molecule_type), &
    2162              :                                         mass(:, molecule_type), loverlap, new_energy, old_energy, &
    2163              :                                         ionode, .FALSE., mol_type(start_mol:end_mol), nchains(:, box_number), &
    2164              :                                         source, group, rng_stream, &
    2165              :                                         avbmc_atom=avbmc_atom(molecule_type), &
    2166              :                                         rmin=avbmc_rmin(molecule_type), rmax=avbmc_rmax(molecule_type), move_type=move_type, &
    2167            0 :                                         target_atom=target_atom)
    2168              : 
    2169              :       END IF
    2170              : 
    2171            0 :       IF (loverlap) THEN
    2172            0 :          DEALLOCATE (r_new)
    2173              : 
    2174              : ! need to reset the old coordinates
    2175            0 :          IF (lbias) THEN
    2176            0 :             CALL force_env_get(bias_env, subsys=subsys)
    2177            0 :             CALL cp_subsys_get(subsys, particles=particles)
    2178              :          ELSE
    2179            0 :             CALL force_env_get(force_env, subsys=subsys)
    2180            0 :             CALL cp_subsys_get(subsys, particles=particles)
    2181              :          END IF
    2182            0 :          DO ipart = 1, nunits_tot(box_number)
    2183            0 :             particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
    2184              :          END DO
    2185              : 
    2186            0 :          CALL timestop(handle)
    2187              : 
    2188              :          RETURN
    2189              :       END IF
    2190              : 
    2191              : ! if we're biasing, we need to compute the new energy with the full
    2192              : ! potential
    2193            0 :       IF (lbias) THEN
    2194              : ! need to give the force_env the coords from the bias_env
    2195            0 :          CALL force_env_get(force_env, subsys=subsys_force)
    2196            0 :          CALL cp_subsys_get(subsys_force, particles=particles_force)
    2197            0 :          CALL force_env_get(bias_env, subsys=subsys)
    2198            0 :          CALL cp_subsys_get(subsys, particles=particles)
    2199            0 :          DO ipart = 1, nunits_tot(box_number)
    2200            0 :             particles_force%els(ipart)%r(1:3) = particles%els(ipart)%r(1:3)
    2201              :          END DO
    2202              : 
    2203              :          CALL force_env_calc_energy_force(force_env, &
    2204            0 :                                           calc_force=.FALSE.)
    2205              :          CALL force_env_get(force_env, &
    2206            0 :                             potential_energy=new_energy)
    2207              : 
    2208              :       END IF
    2209              : 
    2210            0 :       volume_in = 4.0_dp/3.0_dp*pi*(avbmc_rmax(molecule_type)**3 - avbmc_rmin(molecule_type)**3)
    2211            0 :       volume_out = abc(1)*abc(2)*abc(3) - volume_in
    2212              : 
    2213            0 :       IF (lin .AND. move_type == 'in' .OR. &
    2214              :           .NOT. lin .AND. move_type == 'out') THEN
    2215              : ! standard Metropolis rule
    2216              :          prefactor = 1.0_dp
    2217            0 :       ELSEIF (.NOT. lin .AND. move_type == 'in') THEN
    2218            0 :          prefactor = (1.0_dp - pbias(molecule_type))*volume_in/(pbias(molecule_type)*volume_out)
    2219              :       ELSE
    2220            0 :          prefactor = pbias(molecule_type)*volume_out/((1.0_dp - pbias(molecule_type))*volume_in)
    2221              :       END IF
    2222              : 
    2223            0 :       IF (lbias) THEN
    2224              : ! AVBMC with CBMC and a biasing potential...notice that if the biasing
    2225              : ! potential equals the quickstep potential, this cancels out to the
    2226              : ! acceptance below
    2227              :          del_quickstep_energy = (-BETA)*(new_energy - old_energy - &
    2228            0 :                                          (bias_energy_new - bias_energy_old))
    2229              : 
    2230            0 :          IF (del_quickstep_energy .GT. exp_max_val) THEN
    2231            0 :             del_quickstep_energy = max_val
    2232            0 :          ELSEIF (del_quickstep_energy .LT. exp_min_val) THEN
    2233              :             del_quickstep_energy = 0.0_dp
    2234              :          ELSE
    2235            0 :             del_quickstep_energy = EXP(del_quickstep_energy)
    2236              :          END IF
    2237              : 
    2238            0 :          w = prefactor*del_quickstep_energy*weight_new/weight_old
    2239              : 
    2240              :       ELSE
    2241              : 
    2242              : ! AVBMC with CBMC
    2243            0 :          w = prefactor*weight_new/weight_old
    2244              :       END IF
    2245              : 
    2246              : ! check if the move is accepted
    2247            0 :       IF (w .GE. 1.0E0_dp) THEN
    2248            0 :          rand = 0.0E0_dp
    2249              :       ELSE
    2250            0 :          IF (ionode) rand = rng_stream%next()
    2251            0 :          CALL group%bcast(rand, source)
    2252              :       END IF
    2253              : 
    2254            0 :       IF (rand .LT. w) THEN
    2255              : 
    2256              : ! accept the move
    2257              : 
    2258            0 :          IF (lin) THEN
    2259            0 :             IF (move_type == 'in') THEN
    2260              :                moves%avbmc_inin%successes = &
    2261            0 :                   moves%avbmc_inin%successes + 1
    2262              :             ELSE
    2263              :                moves%avbmc_inout%successes = &
    2264            0 :                   moves%avbmc_inout%successes + 1
    2265              :             END IF
    2266              :          ELSE
    2267            0 :             IF (move_type == 'in') THEN
    2268              :                moves%avbmc_outin%successes = &
    2269            0 :                   moves%avbmc_outin%successes + 1
    2270              :             ELSE
    2271              :                moves%avbmc_outout%successes = &
    2272            0 :                   moves%avbmc_outout%successes + 1
    2273              :             END IF
    2274              :          END IF
    2275              : 
    2276              : ! we need to compensate for the fact that we take the difference in
    2277              : ! generate_cbmc_config to keep the exponetials small
    2278            0 :          IF (.NOT. lbias) THEN
    2279            0 :             new_energy = new_energy + old_energy
    2280              :          END IF
    2281              : 
    2282              : ! update energies
    2283            0 :          energy_check = energy_check + (new_energy - old_energy)
    2284            0 :          old_energy = new_energy
    2285              : 
    2286              : ! if we're biasing the update the biasing energy
    2287            0 :          IF (lbias) THEN
    2288              : ! need to do this outside of the routine
    2289            0 :             last_bias_energy = bias_energy_new
    2290            0 :             bias_energy_old = bias_energy_new
    2291              :          END IF
    2292              : 
    2293              : ! update coordinates
    2294            0 :          CALL force_env_get(force_env, subsys=subsys)
    2295            0 :          CALL cp_subsys_get(subsys, particles=particles)
    2296            0 :          DO ipart = 1, nunits_tot(box_number)
    2297            0 :             r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
    2298              :          END DO
    2299              :       ELSE
    2300              : ! reject the move...need to restore the old coordinates
    2301            0 :          IF (lbias) THEN
    2302            0 :             CALL force_env_get(bias_env, subsys=subsys)
    2303            0 :             CALL cp_subsys_get(subsys, particles=particles)
    2304            0 :             DO ipart = 1, nunits_tot(box_number)
    2305            0 :                particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
    2306              :             END DO
    2307            0 :             CALL cp_subsys_set(subsys, particles=particles)
    2308              :          END IF
    2309            0 :          CALL force_env_get(force_env, subsys=subsys)
    2310            0 :          CALL cp_subsys_get(subsys, particles=particles)
    2311            0 :          DO ipart = 1, nunits_tot(box_number)
    2312            0 :             particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
    2313              :          END DO
    2314            0 :          CALL cp_subsys_set(subsys, particles=particles)
    2315              : 
    2316              :       END IF
    2317              : 
    2318              : ! deallocate some stuff
    2319            0 :       DEALLOCATE (r_new)
    2320              : ! end the timing
    2321            0 :       CALL timestop(handle)
    2322              : 
    2323            0 :    END SUBROUTINE mc_avbmc_move
    2324              : 
    2325              : ! **************************************************************************************************
    2326              : !> \brief performs a hybrid Monte Carlo move that runs a short MD sequence
    2327              : !> \param mc_par the mc parameters for the force env
    2328              : !> \param force_env the force environment whose cell we're changing
    2329              : !> \param globenv ...
    2330              : !> \param moves the structure that keeps track of how many moves have been
    2331              : !>               accepted/rejected
    2332              : !> \param move_updates the structure that keeps track of how many moves have
    2333              : !>               been accepted/rejected since the last time the displacements
    2334              : !>               were updated
    2335              : !> \param old_energy the energy of the last accepted move involving an
    2336              : !>                    unbiased calculation
    2337              : !> \param box_number the box we're changing the volume of
    2338              : !> \param energy_check the running total of how much the energy has changed
    2339              : !>                      since the initial configuration
    2340              : !> \param r_old the coordinates of the last accepted move involving an
    2341              : !>               unbiased calculation
    2342              : !> \param rng_stream the random number stream that we draw from
    2343              : !> \author MJM
    2344              : !> \note     Designed for parallel use.
    2345              : ! **************************************************************************************************
    2346           20 :    SUBROUTINE mc_hmc_move(mc_par, force_env, globenv, moves, move_updates, &
    2347              :                           old_energy, box_number, &
    2348           20 :                           energy_check, r_old, rng_stream)
    2349              : 
    2350              :       TYPE(mc_simpar_type), POINTER                      :: mc_par
    2351              :       TYPE(force_env_type), POINTER                      :: force_env
    2352              :       TYPE(global_environment_type), POINTER             :: globenv
    2353              :       TYPE(mc_moves_type), POINTER                       :: moves, move_updates
    2354              :       REAL(KIND=dp), INTENT(INOUT)                       :: old_energy
    2355              :       INTEGER, INTENT(IN)                                :: box_number
    2356              :       REAL(KIND=dp), INTENT(INOUT)                       :: energy_check
    2357              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: r_old
    2358              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
    2359              : 
    2360              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'mc_hmc_move'
    2361              : 
    2362              :       INTEGER                                            :: handle, iatom, source
    2363           20 :       INTEGER, DIMENSION(:), POINTER                     :: nunits_tot
    2364              :       LOGICAL                                            :: ionode
    2365              :       REAL(KIND=dp)                                      :: BETA, energy_term, exp_max_val, &
    2366              :                                                             exp_min_val, new_energy, rand, value, w
    2367           20 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: r
    2368              :       TYPE(cp_subsys_type), POINTER                      :: oldsys
    2369              :       TYPE(mc_ekin_type), POINTER                        :: hmc_ekin
    2370              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
    2371              :       TYPE(mp_comm_type)                                 :: group
    2372              :       TYPE(particle_list_type), POINTER                  :: particles_old
    2373              : 
    2374              : ! begin the timing of the subroutine
    2375              : 
    2376           20 :       CALL timeset(routineN, handle)
    2377              : 
    2378              : ! get a bunch of stuff from mc_par
    2379              :       CALL get_mc_par(mc_par, ionode=ionode, &
    2380              :                       BETA=BETA, exp_max_val=exp_max_val, &
    2381              :                       exp_min_val=exp_min_val, source=source, group=group, &
    2382           20 :                       mc_molecule_info=mc_molecule_info)
    2383           20 :       CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot)
    2384              : 
    2385              : ! nullify some pointers
    2386           20 :       NULLIFY (particles_old, oldsys, hmc_ekin)
    2387              : 
    2388              : ! do some allocation
    2389           60 :       ALLOCATE (r(1:3, 1:nunits_tot(box_number)))
    2390           20 :       ALLOCATE (hmc_ekin)
    2391              : 
    2392              : ! record the attempt
    2393           20 :       moves%hmc%attempts = moves%hmc%attempts + 1
    2394           20 :       move_updates%hmc%attempts = move_updates%hmc%attempts + 1
    2395              : 
    2396              : ! now let's grab the particle positions
    2397           20 :       CALL force_env_get(force_env, subsys=oldsys)
    2398           20 :       CALL cp_subsys_get(oldsys, particles=particles_old)
    2399              : 
    2400              : ! save the old coordinates
    2401        21700 :       DO iatom = 1, nunits_tot(box_number)
    2402        86740 :          r(1:3, iatom) = particles_old%els(iatom)%r(1:3)
    2403              :       END DO
    2404              : 
    2405              : ! now run the MD simulation
    2406           20 :       CALL qs_mol_dyn(force_env, globenv, hmc_e_initial=hmc_ekin%initial_ekin, hmc_e_final=hmc_ekin%final_ekin)
    2407              : 
    2408              : ! get the energy
    2409              :       CALL force_env_get(force_env, &
    2410           20 :                          potential_energy=new_energy)
    2411              : 
    2412              : ! accept or reject the move
    2413              : ! to prevent overflows
    2414           20 :       energy_term = new_energy + hmc_ekin%final_ekin - old_energy - hmc_ekin%initial_ekin
    2415              : 
    2416           20 :       value = -BETA*(energy_term)
    2417           20 :       IF (value .GT. exp_max_val) THEN
    2418              :          w = 10.0_dp
    2419           20 :       ELSEIF (value .LT. exp_min_val) THEN
    2420              :          w = 0.0_dp
    2421              :       ELSE
    2422           20 :          w = EXP(value)
    2423              :       END IF
    2424              : 
    2425           20 :       IF (w .GE. 1.0E0_dp) THEN
    2426           10 :          w = 1.0E0_dp
    2427           10 :          rand = 0.0E0_dp
    2428              :       ELSE
    2429           10 :          IF (ionode) rand = rng_stream%next()
    2430           10 :          CALL group%bcast(rand, source)
    2431              :       END IF
    2432              : 
    2433           20 :       IF (rand .LT. w) THEN
    2434              : 
    2435              : ! accept the move
    2436           14 :          moves%hmc%successes = moves%hmc%successes + 1
    2437           14 :          move_updates%hmc%successes = move_updates%hmc%successes + 1
    2438              : 
    2439              : ! update energies
    2440           14 :          energy_check = energy_check + (new_energy - old_energy)
    2441           14 :          old_energy = new_energy
    2442              : 
    2443        15190 :          DO iatom = 1, nunits_tot(box_number)
    2444        60718 :             r_old(1:3, iatom) = particles_old%els(iatom)%r(1:3)
    2445              :          END DO
    2446              : 
    2447              :       ELSE
    2448              : 
    2449              : ! reset the cell and particle positions
    2450         6510 :          DO iatom = 1, nunits_tot(box_number)
    2451        26022 :             particles_old%els(iatom)%r(1:3) = r_old(1:3, iatom)
    2452              :          END DO
    2453              : 
    2454              :       END IF
    2455              : 
    2456              : ! deallocate some stuff
    2457           20 :       DEALLOCATE (r)
    2458           20 :       DEALLOCATE (hmc_ekin)
    2459              : 
    2460              : ! end the timing
    2461           20 :       CALL timestop(handle)
    2462              : 
    2463           40 :    END SUBROUTINE mc_hmc_move
    2464              : 
    2465              : ! *****************************************************************************
    2466              : !> \brief translates the cluster randomly in either the x,y, or z
    2467              : !>direction
    2468              : !> \param mc_par the mc parameters for the force env
    2469              : !> \param force_env the force environment used in the move
    2470              : !> \param bias_env the force environment used to bias the move, if any (it may
    2471              : !>            be null if lbias=.false. in mc_par)
    2472              : !> \param moves the structure that keeps track of how many moves have been
    2473              : !>               accepted/rejected
    2474              : !> \param move_updates the structure that keeps track of how many moves have
    2475              : !>               been accepted/rejected since the last time the displacements
    2476              : !>               were updated
    2477              : !> \param box_number ...
    2478              : !> \param bias_energy the biased energy of the system before the move
    2479              : !> \param lreject set to .true. if there is an overlap
    2480              : !> \param rng_stream the random number stream that we draw from
    2481              : !> \author Himanshu Goel
    2482              : !> \note     Designed for parallel use.
    2483              : ! **************************************************************************************************
    2484              : 
    2485           10 :    SUBROUTINE mc_cluster_translation(mc_par, force_env, bias_env, moves, &
    2486              :                                      move_updates, box_number, bias_energy, lreject, rng_stream)
    2487              : 
    2488              :       TYPE(mc_simpar_type), POINTER                      :: mc_par
    2489              :       TYPE(force_env_type), POINTER                      :: force_env, bias_env
    2490              :       TYPE(mc_moves_type), POINTER                       :: moves, move_updates
    2491              :       INTEGER, INTENT(IN)                                :: box_number
    2492              :       REAL(KIND=dp), INTENT(INOUT)                       :: bias_energy
    2493              :       LOGICAL, INTENT(OUT)                               :: lreject
    2494              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
    2495              : 
    2496              :       CHARACTER(len=*), PARAMETER :: routineN = 'mc_cluster_translation'
    2497              : 
    2498              :       INTEGER :: cstart, end_mol, handle, imol, ipart, iparticle, iunit, jbox, jpart, junit, &
    2499              :          move_direction, nend, nunit, source, start_mol, total_clus, total_clusafmo
    2500              :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: cluster
    2501           10 :       INTEGER, DIMENSION(:), POINTER                     :: mol_type, nunits, nunits_tot
    2502           10 :       INTEGER, DIMENSION(:, :), POINTER                  :: nchains
    2503              :       LOGICAL                                            :: ionode, lbias, loverlap
    2504              :       REAL(KIND=dp)                                      :: BETA, bias_energy_new, bias_energy_old, &
    2505              :                                                             dis_mol, exp_max_val, exp_min_val, &
    2506              :                                                             rand, rmcltrans, value, w
    2507           10 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: r_old
    2508              :       TYPE(cp_subsys_type), POINTER                      :: subsys
    2509              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
    2510              :       TYPE(mp_comm_type)                                 :: group
    2511              :       TYPE(particle_list_type), POINTER                  :: particles
    2512              : 
    2513              : !   *** Local Counters ***
    2514              : ! begin the timing of the subroutine
    2515              : 
    2516           10 :       CALL timeset(routineN, handle)
    2517              : 
    2518              : ! nullify some pointers
    2519           10 :       NULLIFY (particles, subsys)
    2520              : 
    2521              : ! get a bunch of stuff from mc_par
    2522              :       CALL get_mc_par(mc_par, lbias=lbias, &
    2523              :                       BETA=BETA, exp_max_val=exp_max_val, &
    2524              :                       exp_min_val=exp_min_val, rmcltrans=rmcltrans, ionode=ionode, source=source, &
    2525           10 :                       group=group, mc_molecule_info=mc_molecule_info)
    2526              :       CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot, &
    2527           10 :                                 nchains=nchains, nunits=nunits, mol_type=mol_type)
    2528              : 
    2529              : ! find out some bounds for mol_type
    2530           10 :       start_mol = 1
    2531           10 :       DO jbox = 1, box_number - 1
    2532           10 :          start_mol = start_mol + SUM(nchains(:, jbox))
    2533              :       END DO
    2534           20 :       end_mol = start_mol + SUM(nchains(:, box_number)) - 1
    2535              : 
    2536              : ! do some allocation
    2537           30 :       ALLOCATE (r_old(1:3, 1:nunits_tot(box_number)))
    2538              : 
    2539              : ! Allocating cluster matrix size
    2540           20 :       nend = SUM(nchains(:, box_number))
    2541           40 :       ALLOCATE (cluster(nend, nend))
    2542           60 :       DO ipart = 1, nend
    2543          310 :          DO jpart = 1, nend
    2544          300 :             cluster(ipart, jpart) = 0
    2545              :          END DO
    2546              :       END DO
    2547              : 
    2548              : ! Get cluster information in cluster matrix from cluster_search subroutine
    2549           10 :       IF (lbias) THEN
    2550              :          CALL cluster_search(mc_par, bias_env, cluster, nchains(:, box_number), &
    2551            0 :                              nunits, mol_type(start_mol:end_mol), total_clus)
    2552              :       ELSE
    2553              :          CALL cluster_search(mc_par, force_env, cluster, nchains(:, box_number), &
    2554           10 :                              nunits, mol_type(start_mol:end_mol), total_clus)
    2555              :       END IF
    2556              : 
    2557           10 :       IF (lbias) THEN
    2558              : 
    2559              : ! grab the coordinates
    2560            0 :          CALL force_env_get(bias_env, subsys=subsys)
    2561            0 :          CALL cp_subsys_get(subsys, particles=particles)
    2562              : 
    2563              : ! save the coordinates
    2564            0 :          DO ipart = 1, nunits_tot(box_number)
    2565            0 :             r_old(1:3, ipart) = particles%els(ipart)%r(1:3)
    2566              :          END DO
    2567              : 
    2568              : ! save the energy
    2569            0 :          bias_energy_old = bias_energy
    2570              :       ELSE
    2571              : 
    2572              : ! grab the coordinates
    2573           10 :          CALL force_env_get(force_env, subsys=subsys)
    2574           10 :          CALL cp_subsys_get(subsys, particles=particles)
    2575              :       END IF
    2576              : 
    2577              : ! record the attempt
    2578           10 :       moves%cltrans%attempts = moves%cltrans%attempts + 1
    2579           10 :       move_updates%cltrans%attempts = move_updates%cltrans%attempts + 1
    2580           10 :       moves%bias_cltrans%attempts = moves%bias_cltrans%attempts + 1
    2581           10 :       move_updates%bias_cltrans%attempts = move_updates%bias_cltrans%attempts + 1
    2582           10 :       IF (.NOT. lbias) THEN
    2583           10 :          moves%cltrans%qsuccesses = moves%cltrans%qsuccesses + 1
    2584           10 :          move_updates%cltrans%qsuccesses = move_updates%cltrans%qsuccesses + 1
    2585           10 :          moves%bias_cltrans%qsuccesses = moves%bias_cltrans%qsuccesses + 1
    2586           10 :          move_updates%bias_cltrans%qsuccesses = move_updates%bias_cltrans%qsuccesses + 1
    2587              :       END IF
    2588              : 
    2589              : ! call a random number to figure out which direction we're moving
    2590           10 :       IF (ionode) rand = rng_stream%next()
    2591           10 :       CALL group%bcast(rand, source)
    2592           10 :       move_direction = INT(3*rand) + 1
    2593              : 
    2594              : ! call a random number to figure out how far we're moving
    2595           10 :       IF (ionode) rand = rng_stream%next()
    2596           10 :       CALL group%bcast(rand, source)
    2597           10 :       dis_mol = rmcltrans*(rand - 0.5E0_dp)*2.0E0_dp
    2598              : 
    2599              : ! choosing cluster
    2600           10 :       IF (ionode) rand = rng_stream%next()
    2601           10 :       CALL group%bcast(rand, source)
    2602           10 :       jpart = INT(1 + rand*total_clus)
    2603              : 
    2604              : ! do the cluster move
    2605           60 :       DO cstart = 1, nend
    2606           50 :          imol = 0
    2607           60 :          IF (cluster(jpart, cstart) .NE. 0) THEN
    2608           34 :             imol = cluster(jpart, cstart)
    2609              :             iunit = 1
    2610           34 :             DO ipart = 1, imol - 1
    2611           24 :                nunit = nunits(mol_type(ipart + start_mol - 1))
    2612           34 :                iunit = iunit + nunit
    2613              :             END DO
    2614           10 :             nunit = nunits(mol_type(imol + start_mol - 1))
    2615           10 :             junit = iunit + nunit - 1
    2616           30 :             DO iparticle = iunit, junit
    2617              :                particles%els(iparticle)%r(move_direction) = &
    2618           30 :                   particles%els(iparticle)%r(move_direction) + dis_mol
    2619              :             END DO
    2620              :          END IF
    2621              :       END DO
    2622           10 :       CALL cp_subsys_set(subsys, particles=particles)
    2623              : 
    2624              : !Make cluster matrix null
    2625           60 :       DO ipart = 1, nend
    2626          310 :          DO jpart = 1, nend
    2627          300 :             cluster(ipart, jpart) = 0
    2628              :          END DO
    2629              :       END DO
    2630              : 
    2631              : ! checking the number of cluster are same or got changed after cluster translation move
    2632           10 :       IF (lbias) THEN
    2633              :          CALL cluster_search(mc_par, bias_env, cluster, nchains(:, box_number), &
    2634            0 :                              nunits, mol_type(start_mol:end_mol), total_clusafmo)
    2635              :       ELSE
    2636              :          CALL cluster_search(mc_par, force_env, cluster, nchains(:, box_number), &
    2637           10 :                              nunits, mol_type(start_mol:end_mol), total_clusafmo)
    2638              :       END IF
    2639              : 
    2640              : ! figure out if there is any overlap...need the number of the molecule
    2641           10 :       lreject = .FALSE.
    2642           10 :       IF (lbias) THEN
    2643              :          CALL check_for_overlap(bias_env, nchains(:, box_number), &
    2644            0 :                                 nunits(:), loverlap, mol_type(start_mol:end_mol))
    2645              :       ELSE
    2646              :          CALL check_for_overlap(force_env, nchains(:, box_number), &
    2647           10 :                                 nunits(:), loverlap, mol_type(start_mol:end_mol))
    2648           10 :          IF (loverlap) lreject = .TRUE.
    2649              :       END IF
    2650              : 
    2651              : ! check if cluster size changes then reject the move
    2652           10 :       IF (lbias) THEN
    2653            0 :          IF (total_clusafmo .NE. total_clus) THEN
    2654            0 :             loverlap = .TRUE.
    2655              :          END IF
    2656              :       ELSE
    2657           10 :          IF (total_clusafmo .NE. total_clus) THEN
    2658            0 :             loverlap = .TRUE.
    2659            0 :             lreject = .TRUE.
    2660              :          END IF
    2661              :       END IF
    2662              : 
    2663              : ! if we're biasing with a cheaper potential, check for acceptance
    2664           10 :       IF (lbias) THEN
    2665              : 
    2666              : ! here's where we bias the moves
    2667            0 :          IF (loverlap) THEN
    2668              :             w = 0.0E0_dp
    2669              :          ELSE
    2670            0 :             CALL force_env_calc_energy_force(bias_env, calc_force=.FALSE.)
    2671              :             CALL force_env_get(bias_env, &
    2672            0 :                                potential_energy=bias_energy_new)
    2673              : ! accept or reject the move based on the Metropolis rule
    2674            0 :             value = -BETA*(bias_energy_new - bias_energy_old)
    2675            0 :             IF (value .GT. exp_max_val) THEN
    2676              :                w = 10.0_dp
    2677            0 :             ELSEIF (value .LT. exp_min_val) THEN
    2678              :                w = 0.0_dp
    2679              :             ELSE
    2680            0 :                w = EXP(value)
    2681              :             END IF
    2682              : 
    2683              :          END IF
    2684              : 
    2685            0 :          IF (w .GE. 1.0E0_dp) THEN
    2686            0 :             w = 1.0E0_dp
    2687            0 :             rand = 0.0E0_dp
    2688              :          ELSE
    2689            0 :             IF (ionode) rand = rng_stream%next()
    2690            0 :             CALL group%bcast(rand, source)
    2691              :          END IF
    2692            0 :          IF (rand .LT. w) THEN
    2693              : 
    2694              : ! accept the move
    2695            0 :             moves%bias_cltrans%successes = moves%bias_cltrans%successes + 1
    2696            0 :             move_updates%bias_cltrans%successes = move_updates%bias_cltrans%successes + 1
    2697            0 :             moves%cltrans%qsuccesses = moves%cltrans%qsuccesses + 1
    2698              :             move_updates%cltrans%successes = &
    2699            0 :                move_updates%cltrans%successes + 1
    2700            0 :             moves%qcltrans_dis = moves%qcltrans_dis + ABS(dis_mol)
    2701              :             bias_energy = bias_energy + bias_energy_new - &
    2702            0 :                           bias_energy_old
    2703              : 
    2704              :          ELSE
    2705              : 
    2706              : ! reject the move
    2707              : ! restore the coordinates
    2708            0 :             CALL force_env_get(bias_env, subsys=subsys)
    2709            0 :             CALL cp_subsys_get(subsys, particles=particles)
    2710            0 :             DO ipart = 1, nunits_tot(box_number)
    2711            0 :                particles%els(ipart)%r(1:3) = r_old(1:3, ipart)
    2712              :             END DO
    2713            0 :             CALL cp_subsys_set(subsys, particles=particles)
    2714              : 
    2715              :          END IF
    2716              : 
    2717              :       END IF
    2718              : 
    2719              : ! deallocate some stuff
    2720           10 :       DEALLOCATE (cluster)
    2721           10 :       DEALLOCATE (r_old)
    2722              : 
    2723              : ! end the timing
    2724           10 :       CALL timestop(handle)
    2725              : 
    2726           10 :    END SUBROUTINE mc_cluster_translation
    2727              : 
    2728              : END MODULE mc_moves
        

Generated by: LCOV version 2.0-1