LCOV - code coverage report
Current view: top level - src/motion/mc - mc_ge_moves.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 87.2 % 564 492
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 3 3

            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 contains the Monte Carlo moves that can handle more than one
      10              : !>      box, including the Quickstep move, a volume swap between boxes,
      11              : !>       and a particle swap between boxes
      12              : !> \par History
      13              : !>      MJM (07.28.2005): make the Quickstep move general, and changed
      14              : !>                        the swap and volume moves to work with the
      15              : !>                        CP2K classical routines
      16              : !> \author Matthew J. McGrath  (01.25.2004)
      17              : ! **************************************************************************************************
      18              : MODULE mc_ge_moves
      19              :    USE cell_methods,                    ONLY: cell_create
      20              :    USE cell_types,                      ONLY: cell_clone,&
      21              :                                               cell_p_type,&
      22              :                                               cell_release,&
      23              :                                               cell_type,&
      24              :                                               get_cell
      25              :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      26              :                                               cp_subsys_p_type,&
      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_p_type,&
      32              :                                               force_env_release
      33              :    USE input_constants,                 ONLY: dump_xmol
      34              :    USE input_section_types,             ONLY: section_type
      35              :    USE kinds,                           ONLY: default_string_length,&
      36              :                                               dp
      37              :    USE mc_control,                      ONLY: mc_create_force_env
      38              :    USE mc_coordinates,                  ONLY: check_for_overlap,&
      39              :                                               generate_cbmc_swap_config,&
      40              :                                               get_center_of_mass
      41              :    USE mc_misc,                         ONLY: mc_make_dat_file_new
      42              :    USE mc_move_control,                 ONLY: move_q_reinit,&
      43              :                                               q_move_accept
      44              :    USE mc_types,                        ONLY: &
      45              :         get_mc_molecule_info, get_mc_par, mc_determine_molecule_info, mc_input_file_type, &
      46              :         mc_molecule_info_destroy, mc_molecule_info_type, mc_moves_p_type, &
      47              :         mc_simulation_parameters_p_type, set_mc_par
      48              :    USE message_passing,                 ONLY: mp_comm_type,&
      49              :                                               mp_para_env_type
      50              :    USE parallel_rng_types,              ONLY: rng_stream_type
      51              :    USE particle_list_types,             ONLY: particle_list_p_type,&
      52              :                                               particle_list_type
      53              :    USE particle_methods,                ONLY: write_particle_coordinates
      54              :    USE physcon,                         ONLY: angstrom
      55              : #include "../../base/base_uses.f90"
      56              : 
      57              :    IMPLICIT NONE
      58              : 
      59              :    PRIVATE
      60              : 
      61              : ! *** Global parameters ***
      62              : 
      63              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_ge_moves'
      64              : 
      65              :    PUBLIC :: mc_ge_volume_move, mc_ge_swap_move, &
      66              :              mc_quickstep_move
      67              : 
      68              : CONTAINS
      69              : 
      70              : ! **************************************************************************************************
      71              : !> \brief computes the acceptance of a series of biased or unbiased moves
      72              : !>      (translation, rotation, conformational changes)
      73              : !> \param mc_par the mc parameters for the force envs of the boxes
      74              : !> \param force_env the force environments for the boxes
      75              : !> \param bias_env the force environments with the biasing potential for the boxes
      76              : !> \param moves the structure that keeps track of how many moves have been
      77              : !>               accepted/rejected for both boxes
      78              : !> \param lreject automatically rejects the move (used when an overlap occurs in
      79              : !>                the sequence of moves)
      80              : !> \param move_updates the structure that keeps track of how many moves have
      81              : !>               been accepted/rejected since the last time the displacements
      82              : !>               were updated for both boxes
      83              : !> \param energy_check the running total of how much the energy has changed
      84              : !>                      since the initial configuration
      85              : !> \param r_old the coordinates of the last accepted move before the sequence
      86              : !>         whose acceptance is determined by this call
      87              : !> \param nnstep the Monte Carlo step we're on
      88              : !> \param old_energy the energy of the last accepted move involving the full potential
      89              : !> \param bias_energy_new the energy of the current configuration involving the bias potential
      90              : !> \param last_bias_energy ...
      91              : !> \param nboxes the number of boxes (force environments) in the system
      92              : !> \param box_flag indicates if a move has been tried in a given box..if not, we don't
      93              : !>        recompute the energy
      94              : !> \param subsys the pointers for the particle subsystems of both boxes
      95              : !> \param particles the pointers for the particle sets
      96              : !> \param rng_stream the stream we pull random numbers from
      97              : !> \param unit_conv ...
      98              : !> \author MJM
      99              : ! **************************************************************************************************
     100          368 :    SUBROUTINE mc_Quickstep_move(mc_par, force_env, bias_env, moves, &
     101          368 :                                 lreject, move_updates, energy_check, r_old, &
     102          368 :                                 nnstep, old_energy, bias_energy_new, last_bias_energy, &
     103          368 :                                 nboxes, box_flag, subsys, particles, rng_stream, &
     104              :                                 unit_conv)
     105              : 
     106              :       TYPE(mc_simulation_parameters_p_type), &
     107              :          DIMENSION(:), POINTER                           :: mc_par
     108              :       TYPE(force_env_p_type), DIMENSION(:), POINTER      :: force_env, bias_env
     109              :       TYPE(mc_moves_p_type), DIMENSION(:, :), POINTER    :: moves
     110              :       LOGICAL, INTENT(IN)                                :: lreject
     111              :       TYPE(mc_moves_p_type), DIMENSION(:, :), POINTER    :: move_updates
     112              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: energy_check
     113              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: r_old
     114              :       INTEGER, INTENT(IN)                                :: nnstep
     115              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: old_energy, bias_energy_new, &
     116              :                                                             last_bias_energy
     117              :       INTEGER, INTENT(IN)                                :: nboxes
     118              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: box_flag
     119              :       TYPE(cp_subsys_p_type), DIMENSION(:), POINTER      :: subsys
     120              :       TYPE(particle_list_p_type), DIMENSION(:), POINTER  :: particles
     121              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     122              :       REAL(KIND=dp), INTENT(IN)                          :: unit_conv
     123              : 
     124              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'mc_Quickstep_move'
     125              : 
     126              :       INTEGER                                            :: end_mol, handle, ibox, iparticle, &
     127              :                                                             iprint, itype, jbox, nmol_types, &
     128              :                                                             source, start_mol
     129          368 :       INTEGER, DIMENSION(:, :), POINTER                  :: nchains
     130          368 :       INTEGER, DIMENSION(:), POINTER                     :: mol_type, nunits, nunits_tot
     131          736 :       INTEGER, DIMENSION(1:nboxes)                       :: diff
     132              :       LOGICAL                                            :: ionode, lbias, loverlap
     133              :       REAL(KIND=dp)                                      :: BETA, energies, rand, w
     134          736 :       REAL(KIND=dp), DIMENSION(1:nboxes)                 :: bias_energy_old, new_energy
     135          368 :       TYPE(cp_subsys_p_type), DIMENSION(:), POINTER      :: subsys_bias
     136              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
     137              :       TYPE(mp_comm_type)                                 :: group
     138          368 :       TYPE(particle_list_p_type), DIMENSION(:), POINTER  :: particles_bias
     139              : 
     140              : ! begin the timing of the subroutine
     141              : 
     142          368 :       CALL timeset(routineN, handle)
     143              : 
     144          368 :       NULLIFY (subsys_bias, particles_bias)
     145              : 
     146              : ! get a bunch of data from mc_par
     147              :       CALL get_mc_par(mc_par(1)%mc_par, ionode=ionode, lbias=lbias, &
     148              :                       BETA=BETA, diff=diff(1), source=source, group=group, &
     149              :                       iprint=iprint, &
     150          368 :                       mc_molecule_info=mc_molecule_info)
     151              :       CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types, &
     152          368 :                                 nchains=nchains, nunits_tot=nunits_tot, nunits=nunits, mol_type=mol_type)
     153              : 
     154          368 :       IF (nboxes .GT. 1) THEN
     155          144 :          DO ibox = 2, nboxes
     156          144 :             CALL get_mc_par(mc_par(ibox)%mc_par, diff=diff(ibox))
     157              :          END DO
     158              :       END IF
     159              : 
     160              : ! allocate some stuff
     161         1912 :       ALLOCATE (subsys_bias(1:nboxes))
     162         1176 :       ALLOCATE (particles_bias(1:nboxes))
     163              : 
     164              : ! record the attempt...we really only care about molecule type 1 and box
     165              : ! type 1, since the acceptance will be identical for all boxes and molecules
     166              :       moves(1, 1)%moves%Quickstep%attempts = &
     167          368 :          moves(1, 1)%moves%Quickstep%attempts + 1
     168              : 
     169              : ! grab the coordinates for the force_env
     170          808 :       DO ibox = 1, nboxes
     171              :          CALL force_env_get(force_env(ibox)%force_env, &
     172          440 :                             subsys=subsys(ibox)%subsys)
     173              :          CALL cp_subsys_get(subsys(ibox)%subsys, &
     174          808 :                             particles=particles(ibox)%list)
     175              :       END DO
     176              : 
     177              : ! calculate the new energy of the system...if we're biasing,
     178              : ! force_env hasn't changed but bias_env has
     179          808 :       DO ibox = 1, nboxes
     180          808 :          IF (box_flag(ibox) == 1) THEN
     181          368 :             IF (lbias) THEN
     182              : ! grab the coords from bias_env and put them into force_env
     183              :                CALL force_env_get(bias_env(ibox)%force_env, &
     184           98 :                                   subsys=subsys_bias(ibox)%subsys)
     185              :                CALL cp_subsys_get(subsys_bias(ibox)%subsys, &
     186           98 :                                   particles=particles_bias(ibox)%list)
     187              : 
     188         2606 :                DO iparticle = 1, nunits_tot(ibox)
     189              :                   particles(ibox)%list%els(iparticle)%r(1:3) = &
     190        17654 :                      particles_bias(ibox)%list%els(iparticle)%r(1:3)
     191              :                END DO
     192              : 
     193              :                CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
     194           98 :                                                 calc_force=.FALSE.)
     195              :                CALL force_env_get(force_env(ibox)%force_env, &
     196           98 :                                   potential_energy=new_energy(ibox))
     197              :             ELSE
     198          270 :                IF (.NOT. lreject) THEN
     199              :                   CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
     200          270 :                                                    calc_force=.FALSE.)
     201              :                   CALL force_env_get(force_env(ibox)%force_env, &
     202          270 :                                      potential_energy=new_energy(ibox))
     203              :                END IF
     204              :             END IF
     205              :          ELSE
     206           72 :             new_energy(ibox) = old_energy(ibox)
     207              :          END IF
     208              : 
     209              :       END DO
     210              : 
     211              : ! accept or reject the move based on Metropolis or the Iftimie rule
     212          368 :       IF (ionode) THEN
     213              : 
     214              : ! write them out in case something bad happens
     215          184 :          IF (MOD(nnstep, iprint) == 0) THEN
     216           26 :             DO ibox = 1, nboxes
     217           49 :                IF (SUM(nchains(:, ibox)) == 0) THEN
     218            0 :                   WRITE (diff(ibox), *) nnstep
     219            0 :                   WRITE (diff(ibox), *) nchains(:, ibox)
     220              :                ELSE
     221           14 :                   WRITE (diff(ibox), *) nnstep
     222              :                   CALL write_particle_coordinates( &
     223              :                      particles(ibox)%list%els, &
     224              :                      diff(ibox), dump_xmol, 'POS', 'TRIAL', &
     225           14 :                      unit_conv=unit_conv)
     226              :                END IF
     227              :             END DO
     228              :          END IF
     229              :       END IF
     230              : 
     231          368 :       IF (.NOT. lreject) THEN
     232          368 :          IF (lbias) THEN
     233              : 
     234          268 :             DO ibox = 1, nboxes
     235              : ! look for overlap
     236          510 :                IF (SUM(nchains(:, ibox)) .NE. 0) THEN
     237              : ! find the molecule bounds
     238              :                   start_mol = 1
     239          242 :                   DO jbox = 1, ibox - 1
     240          386 :                      start_mol = start_mol + SUM(nchains(:, jbox))
     241              :                   END DO
     242          510 :                   end_mol = start_mol + SUM(nchains(:, ibox)) - 1
     243              :                   CALL check_for_overlap(bias_env(ibox)%force_env, &
     244          170 :                                          nchains(:, ibox), nunits(:), loverlap, mol_type(start_mol:end_mol))
     245          170 :                   IF (loverlap) &
     246            0 :                      CPABORT('Quickstep move found an overlap in the old config')
     247              :                END IF
     248          268 :                bias_energy_old(ibox) = last_bias_energy(ibox)
     249              :             END DO
     250              : 
     251              :             energies = -BETA*((SUM(new_energy(:)) - SUM(bias_energy_new(:))) &
     252          778 :                               - (SUM(old_energy(:)) - SUM(bias_energy_old(:))))
     253              : 
     254              : ! used to prevent over and underflows
     255           98 :             IF (energies .GE. -1.0E-8) THEN
     256              :                w = 1.0_dp
     257           34 :             ELSEIF (energies .LE. -500.0_dp) THEN
     258              :                w = 0.0_dp
     259              :             ELSE
     260           34 :                w = EXP(energies)
     261              :             END IF
     262              : 
     263           98 :             IF (ionode) THEN
     264          134 :                DO ibox = 1, nboxes
     265           85 :                   WRITE (diff(ibox), *) nnstep, new_energy(ibox) - &
     266           85 :                      old_energy(ibox), &
     267          219 :                      bias_energy_new(ibox) - bias_energy_old(ibox)
     268              :                END DO
     269              :             END IF
     270              :          ELSE
     271          810 :             energies = -BETA*(SUM(new_energy(:)) - SUM(old_energy(:)))
     272              : ! used to prevent over and underflows
     273          270 :             IF (energies .GE. 0.0_dp) THEN
     274              :                w = 1.0_dp
     275          156 :             ELSEIF (energies .LE. -500.0_dp) THEN
     276              :                w = 0.0_dp
     277              :             ELSE
     278          156 :                w = EXP(energies)
     279              :             END IF
     280              :          END IF
     281              :       ELSE
     282              :          w = 0.0E0_dp
     283              :       END IF
     284          254 :       IF (w .GE. 1.0E0_dp) THEN
     285          178 :          w = 1.0E0_dp
     286          178 :          rand = 0.0E0_dp
     287              :       ELSE
     288          190 :          IF (ionode) rand = rng_stream%next()
     289          190 :          CALL group%bcast(rand, source)
     290              :       END IF
     291              : 
     292          368 :       IF (rand .LT. w) THEN
     293              : 
     294              : ! accept the move
     295              :          moves(1, 1)%moves%Quickstep%successes = &
     296          286 :             moves(1, 1)%moves%Quickstep%successes + 1
     297              : 
     298          644 :          DO ibox = 1, nboxes
     299              : ! remember what kind of move we did for lbias=.false.
     300          358 :             IF (.NOT. lbias) THEN
     301          554 :                DO itype = 1, nmol_types
     302          366 :                   CALL q_move_accept(moves(itype, ibox)%moves, .TRUE.)
     303          366 :                   CALL q_move_accept(move_updates(itype, ibox)%moves, .TRUE.)
     304              : 
     305              : ! reset the counters
     306          366 :                   CALL move_q_reinit(moves(itype, ibox)%moves, .TRUE.)
     307          554 :                   CALL move_q_reinit(move_updates(itype, ibox)%moves, .TRUE.)
     308              :                END DO
     309              :             END IF
     310              : 
     311         1064 :             DO itype = 1, nmol_types
     312              : ! we need to record all accepted moves since last Quickstep calculation
     313          706 :                CALL q_move_accept(moves(itype, ibox)%moves, .FALSE.)
     314          706 :                CALL q_move_accept(move_updates(itype, ibox)%moves, .FALSE.)
     315              : 
     316              : ! reset the counters
     317          706 :                CALL move_q_reinit(moves(itype, ibox)%moves, .FALSE.)
     318         1064 :                CALL move_q_reinit(move_updates(itype, ibox)%moves, .FALSE.)
     319              :             END DO
     320              : 
     321              : ! update energies
     322              :             energy_check(ibox) = energy_check(ibox) + &
     323          358 :                                  (new_energy(ibox) - old_energy(ibox))
     324          644 :             old_energy(ibox) = new_energy(ibox)
     325              : 
     326              :          END DO
     327              : 
     328          286 :          IF (lbias) THEN
     329          268 :             DO ibox = 1, nboxes
     330          268 :                last_bias_energy(ibox) = bias_energy_new(ibox)
     331              :             END DO
     332              :          END IF
     333              : 
     334              : ! update coordinates
     335          644 :          DO ibox = 1, nboxes
     336          644 :             IF (nunits_tot(ibox) .NE. 0) THEN
     337         9566 :                DO iparticle = 1, nunits_tot(ibox)
     338              :                   r_old(1:3, iparticle, ibox) = &
     339        37190 :                      particles(ibox)%list%els(iparticle)%r(1:3)
     340              :                END DO
     341              :             END IF
     342              :          END DO
     343              : 
     344              :       ELSE
     345              : 
     346              :          ! reject the move
     347          164 :          DO ibox = 1, nboxes
     348          328 :             DO itype = 1, nmol_types
     349          164 :                CALL move_q_reinit(moves(itype, ibox)%moves, .FALSE.)
     350          164 :                CALL move_q_reinit(move_updates(itype, ibox)%moves, .FALSE.)
     351          246 :                IF (.NOT. lbias) THEN
     352              : ! reset the counters
     353          164 :                   CALL move_q_reinit(moves(itype, ibox)%moves, .TRUE.)
     354          164 :                   CALL move_q_reinit(move_updates(itype, ibox)%moves, .TRUE.)
     355              :                END IF
     356              :             END DO
     357              : 
     358              :          END DO
     359              : 
     360         4551 :          IF (.NOT. ionode) r_old(:, :, :) = 0.0E0_dp
     361              : 
     362              : ! coodinates changed, so we need to broadcast those, even for the lbias
     363              : ! case since bias_env needs to have the same coords as force_env
     364        17958 :          CALL group%bcast(r_old, source)
     365              : 
     366          164 :          DO ibox = 1, nboxes
     367         2378 :             DO iparticle = 1, nunits_tot(ibox)
     368              :                particles(ibox)%list%els(iparticle)%r(1:3) = &
     369         8856 :                   r_old(1:3, iparticle, ibox)
     370         2214 :                IF (lbias .AND. box_flag(ibox) == 1) &
     371              :                   particles_bias(ibox)%list%els(iparticle)%r(1:3) = &
     372           82 :                   r_old(1:3, iparticle, ibox)
     373              :             END DO
     374              :          END DO
     375              : 
     376              : ! need to reset the energies of the biasing potential
     377           82 :          IF (lbias) THEN
     378            0 :             DO ibox = 1, nboxes
     379            0 :                bias_energy_new(ibox) = last_bias_energy(ibox)
     380              :             END DO
     381              :          END IF
     382              : 
     383              :       END IF
     384              : 
     385              : ! make sure the coordinates are transferred
     386          808 :       DO ibox = 1, nboxes
     387              :          CALL cp_subsys_set(subsys(ibox)%subsys, &
     388          440 :                             particles=particles(ibox)%list)
     389          440 :          IF (lbias .AND. box_flag(ibox) == 1) &
     390              :             CALL cp_subsys_set(subsys_bias(ibox)%subsys, &
     391          466 :                                particles=particles_bias(ibox)%list)
     392              :       END DO
     393              : 
     394              :       ! deallocate some stuff
     395          368 :       DEALLOCATE (subsys_bias)
     396          368 :       DEALLOCATE (particles_bias)
     397              : 
     398              : ! end the timing
     399          368 :       CALL timestop(handle)
     400              : 
     401          368 :    END SUBROUTINE mc_Quickstep_move
     402              : 
     403              : ! **************************************************************************************************
     404              : !> \brief attempts a swap move between two simulation boxes
     405              : !> \param mc_par the mc parameters for the force envs of the boxes
     406              : !> \param force_env the force environments for the boxes
     407              : !> \param bias_env the force environments used to bias moves for the boxes
     408              : !> \param moves the structure that keeps track of how many moves have been
     409              : !>               accepted/rejected for both boxes
     410              : !> \param energy_check the running total of how much the energy has changed
     411              : !>                      since the initial configuration
     412              : !> \param r_old the coordinates of the last accepted move involving a
     413              : !>               full potential calculation for both boxes
     414              : !> \param old_energy the energy of the last accepted move involving a
     415              : !>                    a full potential calculation
     416              : !> \param input_declaration ...
     417              : !> \param para_env the parallel environment for this simulation
     418              : !> \param bias_energy_old the energies of both boxes computed using the biasing
     419              : !>        potential
     420              : !> \param last_bias_energy the energy for the biased simulations
     421              : !> \param rng_stream the stream we pull random numbers from
     422              : !> \author MJM
     423              : ! **************************************************************************************************
     424           22 :    SUBROUTINE mc_ge_swap_move(mc_par, force_env, bias_env, moves, &
     425           22 :                               energy_check, r_old, old_energy, input_declaration, &
     426              :                               para_env, bias_energy_old, last_bias_energy, &
     427              :                               rng_stream)
     428              : 
     429              :       TYPE(mc_simulation_parameters_p_type), &
     430              :          DIMENSION(:), POINTER                           :: mc_par
     431              :       TYPE(force_env_p_type), DIMENSION(:), POINTER      :: force_env, bias_env
     432              :       TYPE(mc_moves_p_type), DIMENSION(:, :), POINTER    :: moves
     433              :       REAL(KIND=dp), DIMENSION(1:2), INTENT(INOUT)       :: energy_check
     434              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: r_old
     435              :       REAL(KIND=dp), DIMENSION(1:2), INTENT(INOUT)       :: old_energy
     436              :       TYPE(section_type), POINTER                        :: input_declaration
     437              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     438              :       REAL(KIND=dp), DIMENSION(1:2), INTENT(INOUT)       :: bias_energy_old, last_bias_energy
     439              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     440              : 
     441              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'mc_ge_swap_move'
     442              : 
     443              :       CHARACTER(default_string_length), ALLOCATABLE, &
     444           22 :          DIMENSION(:)                                    :: atom_names_insert, atom_names_remove
     445              :       CHARACTER(default_string_length), &
     446           22 :          DIMENSION(:, :), POINTER                        :: atom_names
     447              :       CHARACTER(LEN=200)                                 :: fft_lib
     448              :       CHARACTER(LEN=40), DIMENSION(1:2)                  :: dat_file
     449              :       INTEGER :: end_mol, handle, iatom, ibox, idim, iiatom, imolecule, ins_atoms, insert_box, &
     450              :          ipart, itype, jbox, molecule_type, nmol_types, nswapmoves, print_level, rem_atoms, &
     451              :          remove_box, source, start_atom_ins, start_atom_rem, start_mol
     452           22 :       INTEGER, DIMENSION(:), POINTER                     :: mol_type, mol_type_test, nunits, &
     453           22 :                                                             nunits_tot
     454           22 :       INTEGER, DIMENSION(:, :), POINTER                  :: nchains, nchains_test
     455              :       LOGICAL                                            :: ionode, lbias, loverlap, loverlap_ins, &
     456              :                                                             loverlap_rem
     457           22 :       REAL(dp), DIMENSION(:), POINTER                    :: eta_insert, eta_remove, pmswap_mol
     458           22 :       REAL(dp), DIMENSION(:, :), POINTER                 :: insert_coords, remove_coords
     459              :       REAL(KIND=dp) :: BETA, del_quickstep_energy, exp_max_val, exp_min_val, max_val, min_val, &
     460              :          prefactor, rand, rdum, vol_insert, vol_remove, w, weight_new, weight_old
     461           22 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: cbmc_energies, r_cbmc, r_insert_mol
     462              :       REAL(KIND=dp), DIMENSION(1:2)                      :: bias_energy_new, new_energy
     463              :       REAL(KIND=dp), DIMENSION(1:3)                      :: abc_insert, abc_remove, center_of_mass, &
     464              :                                                             displace_molecule, pos_insert
     465           22 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: mass
     466              :       TYPE(cell_type), POINTER                           :: cell_insert, cell_remove
     467           22 :       TYPE(cp_subsys_p_type), DIMENSION(:), POINTER      :: oldsys
     468              :       TYPE(cp_subsys_type), POINTER                      :: insert_sys, remove_sys
     469           22 :       TYPE(force_env_p_type), DIMENSION(:), POINTER      :: test_env, test_env_bias
     470              :       TYPE(mc_input_file_type), POINTER                  :: mc_bias_file, mc_input_file
     471              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info, mc_molecule_info_test
     472              :       TYPE(mp_comm_type)                                 :: group
     473           22 :       TYPE(particle_list_p_type), DIMENSION(:), POINTER  :: particles_old
     474              :       TYPE(particle_list_type), POINTER                  :: particles_insert, particles_remove
     475              : 
     476              : ! begin the timing of the subroutine
     477              : 
     478           22 :       CALL timeset(routineN, handle)
     479              : 
     480              : ! reset the overlap flag
     481           22 :       loverlap = .FALSE.
     482              : 
     483              : ! nullify some pointers
     484           22 :       NULLIFY (particles_old, mol_type, mol_type_test, mc_input_file, mc_bias_file)
     485           22 :       NULLIFY (oldsys, atom_names, pmswap_mol, insert_coords, remove_coords)
     486           22 :       NULLIFY (eta_insert, eta_remove)
     487              : 
     488              : ! grab some stuff from mc_par
     489              :       CALL get_mc_par(mc_par(1)%mc_par, ionode=ionode, BETA=BETA, &
     490              :                       max_val=max_val, min_val=min_val, exp_max_val=exp_max_val, &
     491              :                       exp_min_val=exp_min_val, nswapmoves=nswapmoves, group=group, source=source, &
     492              :                       lbias=lbias, dat_file=dat_file(1), fft_lib=fft_lib, &
     493           22 :                       mc_molecule_info=mc_molecule_info, pmswap_mol=pmswap_mol)
     494              :       CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, &
     495              :                                 nunits=nunits, nunits_tot=nunits_tot, nmol_types=nmol_types, &
     496           22 :                                 atom_names=atom_names, mass=mass, mol_type=mol_type)
     497              : 
     498           22 :       print_level = 1
     499              : 
     500           22 :       CALL get_mc_par(mc_par(2)%mc_par, dat_file=dat_file(2))
     501              : 
     502              : ! allocate some stuff
     503           66 :       ALLOCATE (oldsys(1:2))
     504           66 :       ALLOCATE (particles_old(1:2))
     505              : 
     506              : ! get the old coordinates
     507           66 :       DO ibox = 1, 2
     508              :          CALL force_env_get(force_env(ibox)%force_env, &
     509           44 :                             subsys=oldsys(ibox)%subsys)
     510              :          CALL cp_subsys_get(oldsys(ibox)%subsys, &
     511           66 :                             particles=particles_old(ibox)%list)
     512              :       END DO
     513              : 
     514              : !     choose a direction to swap
     515           22 :       IF (ionode) rand = rng_stream%next()
     516           22 :       CALL group%bcast(rand, source)
     517              : 
     518           22 :       IF (rand .LE. 0.50E0_dp) THEN
     519           12 :          remove_box = 1
     520           12 :          insert_box = 2
     521              :       ELSE
     522           10 :          remove_box = 2
     523           10 :          insert_box = 1
     524              :       END IF
     525              : 
     526              : ! now assign the eta values for the insert and remove boxes
     527           22 :       CALL get_mc_par(mc_par(remove_box)%mc_par, eta=eta_remove)
     528           22 :       CALL get_mc_par(mc_par(insert_box)%mc_par, eta=eta_insert)
     529              : 
     530              : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! testing
     531              : !      remove_box=2
     532              : !      insert_box=1
     533              : 
     534              : ! now choose a molecule type at random
     535           22 :       IF (ionode) rand = rng_stream%next()
     536           22 :       CALL group%bcast(rand, source)
     537           34 :       DO itype = 1, nmol_types
     538           34 :          IF (rand .LT. pmswap_mol(itype)) THEN
     539              :             molecule_type = itype
     540              :             EXIT
     541              :          END IF
     542              :       END DO
     543              : 
     544              : ! record the attempt for the box the particle is to be inserted into
     545              :       moves(molecule_type, insert_box)%moves%swap%attempts = &
     546           22 :          moves(molecule_type, insert_box)%moves%swap%attempts + 1
     547              : 
     548              : ! now choose a random molecule to remove from the removal box, checking
     549              : ! to make sure the box isn't empty
     550           22 :       IF (nchains(molecule_type, remove_box) == 0) THEN
     551            0 :          loverlap = .TRUE.
     552              :          moves(molecule_type, insert_box)%moves%empty = &
     553            0 :             moves(molecule_type, insert_box)%moves%empty + 1
     554              :       ELSE
     555              : 
     556           22 :          IF (ionode) rand = rng_stream%next()
     557           22 :          CALL group%bcast(rand, source)
     558           22 :          imolecule = CEILING(rand*nchains(molecule_type, remove_box))
     559              : ! figure out the atom number this molecule starts on
     560           22 :          start_atom_rem = 1
     561           34 :          DO itype = 1, nmol_types
     562           34 :             IF (itype == molecule_type) THEN
     563           22 :                start_atom_rem = start_atom_rem + (imolecule - 1)*nunits(itype)
     564           22 :                EXIT
     565              :             ELSE
     566           12 :                start_atom_rem = start_atom_rem + nchains(itype, remove_box)*nunits(itype)
     567              :             END IF
     568              :          END DO
     569              : 
     570              : ! check for overlap
     571           22 :          start_mol = 1
     572           32 :          DO jbox = 1, remove_box - 1
     573           52 :             start_mol = start_mol + SUM(nchains(:, jbox))
     574              :          END DO
     575           66 :          end_mol = start_mol + SUM(nchains(:, remove_box)) - 1
     576              :          CALL check_for_overlap(force_env(remove_box)%force_env, &
     577           22 :                                 nchains(:, remove_box), nunits, loverlap, mol_type(start_mol:end_mol))
     578           22 :          IF (loverlap) CALL cp_abort(__LOCATION__, &
     579            0 :                                      'CBMC swap move found an overlap in the old remove config')
     580           22 :          start_mol = 1
     581           34 :          DO jbox = 1, insert_box - 1
     582           58 :             start_mol = start_mol + SUM(nchains(:, jbox))
     583              :          END DO
     584           66 :          end_mol = start_mol + SUM(nchains(:, insert_box)) - 1
     585              :          CALL check_for_overlap(force_env(insert_box)%force_env, &
     586           22 :                                 nchains(:, insert_box), nunits, loverlap, mol_type(start_mol:end_mol))
     587           22 :          IF (loverlap) CALL cp_abort(__LOCATION__, &
     588            0 :                                      'CBMC swap move found an overlap in the old insert config')
     589              :       END IF
     590              : 
     591           22 :       IF (loverlap) THEN
     592            0 :          DEALLOCATE (oldsys)
     593            0 :          DEALLOCATE (particles_old)
     594            0 :          CALL timestop(handle)
     595            0 :          RETURN
     596              :       END IF
     597              : 
     598              : ! figure out how many atoms will be in each box after the move
     599           22 :       ins_atoms = nunits_tot(insert_box) + nunits(molecule_type)
     600           22 :       rem_atoms = nunits_tot(remove_box) - nunits(molecule_type)
     601              : ! now allocate the arrays that will hold the coordinates and the
     602              : ! atom name, for writing to the dat file
     603           22 :       IF (rem_atoms == 0) THEN
     604            0 :          ALLOCATE (remove_coords(1:3, 1:nunits(1)))
     605            0 :          ALLOCATE (atom_names_remove(1:nunits(1)))
     606              :       ELSE
     607           66 :          ALLOCATE (remove_coords(1:3, 1:rem_atoms))
     608           66 :          ALLOCATE (atom_names_remove(1:rem_atoms))
     609              :       END IF
     610           66 :       ALLOCATE (insert_coords(1:3, 1:ins_atoms))
     611           66 :       ALLOCATE (atom_names_insert(1:ins_atoms))
     612              : 
     613              : ! grab the cells for later...acceptance and insertion
     614           22 :       IF (lbias) THEN
     615              :          CALL force_env_get(bias_env(insert_box)%force_env, &
     616           22 :                             cell=cell_insert)
     617              :          CALL force_env_get(bias_env(remove_box)%force_env, &
     618           22 :                             cell=cell_remove)
     619              :       ELSE
     620              :          CALL force_env_get(force_env(insert_box)%force_env, &
     621            0 :                             cell=cell_insert)
     622              :          CALL force_env_get(force_env(remove_box)%force_env, &
     623            0 :                             cell=cell_remove)
     624              :       END IF
     625           22 :       CALL get_cell(cell_remove, abc=abc_remove, deth=vol_remove)
     626           22 :       CALL get_cell(cell_insert, abc=abc_insert, deth=vol_insert)
     627              : 
     628           22 :       IF (ionode) THEN
     629              : ! choose an insertion point
     630           44 :          DO idim = 1, 3
     631           33 :             rand = rng_stream%next()
     632           44 :             pos_insert(idim) = rand*abc_insert(idim)
     633              :          END DO
     634              :       END IF
     635           22 :       CALL group%bcast(pos_insert, source)
     636              : 
     637              : ! allocate some arrays we'll be using
     638           66 :       ALLOCATE (r_insert_mol(1:3, 1:nunits(molecule_type)))
     639              : 
     640           22 :       iiatom = 1
     641           64 :       DO iatom = start_atom_rem, start_atom_rem + nunits(molecule_type) - 1
     642              :          r_insert_mol(1:3, iiatom) = &
     643          168 :             particles_old(remove_box)%list%els(iatom)%r(1:3)
     644           64 :          iiatom = iiatom + 1
     645              :       END DO
     646              : 
     647              : !     find the center of mass of the molecule
     648              :       CALL get_center_of_mass(r_insert_mol(:, :), nunits(molecule_type), &
     649           22 :                               center_of_mass(:), mass(:, molecule_type))
     650              : 
     651              : !     move the center of mass to the insertion point
     652           88 :       displace_molecule(1:3) = pos_insert(1:3) - center_of_mass(1:3)
     653           64 :       DO iatom = 1, nunits(molecule_type)
     654              :          r_insert_mol(1:3, iatom) = r_insert_mol(1:3, iatom) + &
     655          190 :                                     displace_molecule(1:3)
     656              :       END DO
     657              : 
     658              : ! prepare the insertion coordinates to be written to the .dat file so
     659              : ! we can create a new force environment...remember there is still a particle
     660              : ! in the box even if nchain=0
     661           66 :       IF (SUM(nchains(:, insert_box)) == 0) THEN
     662            0 :          DO iatom = 1, nunits(molecule_type)
     663            0 :             insert_coords(1:3, iatom) = r_insert_mol(1:3, iatom)
     664              :             atom_names_insert(iatom) = &
     665            0 :                particles_old(remove_box)%list%els(start_atom_rem + iatom - 1)%atomic_kind%name
     666              :          END DO
     667            0 :          start_atom_ins = 1
     668              :       ELSE
     669              : ! the problem is I can't just tack the new molecule on to the end,
     670              : ! because of reading in the dat_file...the topology stuff will crash
     671              : ! if the molecules aren't all grouped together, so I have to insert it
     672              : ! at the end of the section of molecules with the same type, then
     673              : ! remember the start number for the CBMC stuff
     674           22 :          start_atom_ins = 1
     675           34 :          DO itype = 1, nmol_types
     676              :             start_atom_ins = start_atom_ins + &
     677           34 :                              nchains(itype, insert_box)*nunits(itype)
     678           34 :             IF (itype == molecule_type) EXIT
     679              :          END DO
     680              : 
     681          504 :          DO iatom = 1, start_atom_ins - 1
     682              :             insert_coords(1:3, iatom) = &
     683         1928 :                particles_old(insert_box)%list%els(iatom)%r(1:3)
     684              :             atom_names_insert(iatom) = &
     685          504 :                particles_old(insert_box)%list%els(iatom)%atomic_kind%name
     686              :          END DO
     687           22 :          iiatom = 1
     688           64 :          DO iatom = start_atom_ins, start_atom_ins + nunits(molecule_type) - 1
     689          168 :             insert_coords(1:3, iatom) = r_insert_mol(1:3, iiatom)
     690           42 :             atom_names_insert(iatom) = atom_names(iiatom, molecule_type)
     691           64 :             iiatom = iiatom + 1
     692              :          END DO
     693           76 :          DO iatom = start_atom_ins + nunits(molecule_type), ins_atoms
     694              :             insert_coords(1:3, iatom) = &
     695          216 :                particles_old(insert_box)%list%els(iatom - nunits(molecule_type))%r(1:3)
     696              :             atom_names_insert(iatom) = &
     697           76 :                particles_old(insert_box)%list%els(iatom - nunits(molecule_type))%atomic_kind%name
     698              :          END DO
     699              :       END IF
     700              : 
     701              : ! fold the coordinates into the box and check for overlaps
     702           22 :       start_mol = 1
     703           34 :       DO jbox = 1, insert_box - 1
     704           58 :          start_mol = start_mol + SUM(nchains(:, jbox))
     705              :       END DO
     706           66 :       end_mol = start_mol + SUM(nchains(:, insert_box)) - 1
     707              : 
     708              : ! make the .dat file
     709           22 :       IF (ionode) THEN
     710              : 
     711           11 :          nchains(molecule_type, insert_box) = nchains(molecule_type, insert_box) + 1
     712           11 :          IF (lbias) THEN
     713           11 :             CALL get_mc_par(mc_par(insert_box)%mc_par, mc_bias_file=mc_bias_file)
     714              :             CALL mc_make_dat_file_new(insert_coords(:, :), atom_names_insert(:), ins_atoms, &
     715              :                                       abc_insert(:), dat_file(insert_box), nchains(:, insert_box), &
     716           11 :                                       mc_bias_file)
     717              :          ELSE
     718            0 :             CALL get_mc_par(mc_par(insert_box)%mc_par, mc_input_file=mc_input_file)
     719              :             CALL mc_make_dat_file_new(insert_coords(:, :), atom_names_insert(:), ins_atoms, &
     720              :                                       abc_insert(:), dat_file(insert_box), nchains(:, insert_box), &
     721            0 :                                       mc_input_file)
     722              :          END IF
     723           11 :          nchains(molecule_type, insert_box) = nchains(molecule_type, insert_box) - 1
     724              : 
     725              :       END IF
     726              : 
     727              : ! now do the same for the removal box...be careful not to make an empty box
     728           22 :       IF (rem_atoms == 0) THEN
     729            0 :          DO iatom = 1, nunits(molecule_type)
     730            0 :             remove_coords(1:3, iatom) = r_insert_mol(1:3, iatom)
     731            0 :             atom_names_remove(iatom) = atom_names(iatom, molecule_type)
     732              :          END DO
     733              : 
     734              : ! need to adjust nchains, because otherwise if we are removing a molecule type
     735              : ! that is not the first molecule, the dat file will have two molecules in it but
     736              : ! only the coordinates for one
     737            0 :          nchains(molecule_type, remove_box) = nchains(molecule_type, remove_box) - 1
     738            0 :          IF (ionode) THEN
     739            0 :             IF (lbias) THEN
     740            0 :                CALL get_mc_par(mc_par(remove_box)%mc_par, mc_bias_file=mc_bias_file)
     741              :                CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
     742              :                                          abc_remove(:), dat_file(remove_box), nchains(:, remove_box), &
     743            0 :                                          mc_bias_file)
     744              :             ELSE
     745            0 :                CALL get_mc_par(mc_par(remove_box)%mc_par, mc_input_file=mc_input_file)
     746              :                CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
     747              :                                          abc_remove(:), dat_file(remove_box), nchains(:, remove_box), &
     748            0 :                                          mc_input_file)
     749              :             END IF
     750              : 
     751              :          END IF
     752            0 :          nchains(molecule_type, remove_box) = nchains(molecule_type, remove_box) + 1
     753              : 
     754              :       ELSE
     755          368 :          DO iatom = 1, start_atom_rem - 1
     756              :             remove_coords(1:3, iatom) = &
     757         1384 :                particles_old(remove_box)%list%els(iatom)%r(1:3)
     758              :             atom_names_remove(iatom) = &
     759          368 :                particles_old(remove_box)%list%els(iatom)%atomic_kind%name
     760              :          END DO
     761          198 :          DO iatom = start_atom_rem + nunits(molecule_type), nunits_tot(remove_box)
     762              :             remove_coords(1:3, iatom - nunits(molecule_type)) = &
     763          704 :                particles_old(remove_box)%list%els(iatom)%r(1:3)
     764              :             atom_names_remove(iatom - nunits(molecule_type)) = &
     765          198 :                particles_old(remove_box)%list%els(iatom)%atomic_kind%name
     766              :          END DO
     767              : 
     768              : ! make the .dat file
     769           22 :          IF (ionode) THEN
     770           11 :             nchains(molecule_type, remove_box) = nchains(molecule_type, remove_box) - 1
     771           11 :             IF (lbias) THEN
     772           11 :                CALL get_mc_par(mc_par(remove_box)%mc_par, mc_bias_file=mc_bias_file)
     773              :                CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
     774              :                                          abc_remove(:), dat_file(remove_box), nchains(:, remove_box), &
     775           11 :                                          mc_bias_file)
     776              :             ELSE
     777            0 :                CALL get_mc_par(mc_par(remove_box)%mc_par, mc_input_file=mc_input_file)
     778              :                CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
     779              :                                          abc_remove(:), dat_file(remove_box), nchains(:, remove_box), &
     780            0 :                                          mc_input_file)
     781              :             END IF
     782           11 :             nchains(molecule_type, remove_box) = nchains(molecule_type, remove_box) + 1
     783              : 
     784              :          END IF
     785              :       END IF
     786              : 
     787              : ! deallocate r_insert_mol
     788           22 :       DEALLOCATE (r_insert_mol)
     789              : 
     790              : ! now let's create the two new environments with the different number
     791              : ! of molecules
     792           66 :       ALLOCATE (test_env(1:2))
     793              :       CALL mc_create_force_env(test_env(insert_box)%force_env, input_declaration, &
     794           22 :                                para_env, dat_file(insert_box))
     795              :       CALL mc_create_force_env(test_env(remove_box)%force_env, input_declaration, &
     796           22 :                                para_env, dat_file(remove_box))
     797              : 
     798              : ! allocate an array we'll need
     799           44 :       ALLOCATE (r_cbmc(1:3, 1:ins_atoms))
     800           66 :       ALLOCATE (cbmc_energies(1:nswapmoves, 1:2))
     801              : 
     802           22 :       loverlap_ins = .FALSE.
     803           22 :       loverlap_rem = .FALSE.
     804              : 
     805              : ! compute the new molecule information...we need this for the CBMC part
     806           22 :       IF (rem_atoms == 0) THEN
     807              :          CALL mc_determine_molecule_info(test_env, mc_molecule_info_test, &
     808            0 :                                          box_number=remove_box)
     809              :       ELSE
     810           22 :          CALL mc_determine_molecule_info(test_env, mc_molecule_info_test)
     811              :       END IF
     812              :       CALL get_mc_molecule_info(mc_molecule_info_test, nchains=nchains_test, &
     813           22 :                                 mol_type=mol_type_test)
     814              : 
     815              : ! figure out the position of the molecule we're inserting, and the
     816              : ! Rosenbluth weight
     817           22 :       start_mol = 1
     818           34 :       DO jbox = 1, insert_box - 1
     819           58 :          start_mol = start_mol + SUM(nchains_test(:, jbox))
     820              :       END DO
     821           66 :       end_mol = start_mol + SUM(nchains_test(:, insert_box)) - 1
     822              : 
     823           22 :       IF (lbias) THEN
     824              :          CALL generate_cbmc_swap_config(test_env(insert_box)%force_env, &
     825              :                                         BETA, max_val, min_val, exp_max_val, &
     826              :                                         exp_min_val, nswapmoves, weight_new, start_atom_ins, ins_atoms, nunits(:), &
     827              :                                         nunits(molecule_type), mass(:, molecule_type), loverlap_ins, bias_energy_new(insert_box), &
     828              :                                         bias_energy_old(insert_box), ionode, .FALSE., mol_type_test(start_mol:end_mol), &
     829           22 :                                         nchains_test(:, insert_box), source, group, rng_stream)
     830              : 
     831              : ! the energy that comes out of the above routine is the difference...we want
     832              : ! the real energy for the acceptance rule...we don't do this for the
     833              : ! lbias=.false. case because it doesn't appear in the acceptance rule, and
     834              : ! we compensate in case of acceptance
     835              :          bias_energy_new(insert_box) = bias_energy_new(insert_box) + &
     836           22 :                                        bias_energy_old(insert_box)
     837              :       ELSE
     838              :          CALL generate_cbmc_swap_config(test_env(insert_box)%force_env, &
     839              :                                         BETA, max_val, min_val, exp_max_val, &
     840              :                                         exp_min_val, nswapmoves, weight_new, start_atom_ins, ins_atoms, nunits(:), &
     841              :                                         nunits(molecule_type), mass(:, molecule_type), loverlap_ins, new_energy(insert_box), &
     842              :                                         old_energy(insert_box), ionode, .FALSE., mol_type_test(start_mol:end_mol), &
     843            0 :                                         nchains_test(:, insert_box), source, group, rng_stream)
     844              :       END IF
     845              : 
     846              :       CALL force_env_get(test_env(insert_box)%force_env, &
     847           22 :                          subsys=insert_sys)
     848              :       CALL cp_subsys_get(insert_sys, &
     849           22 :                          particles=particles_insert)
     850              : 
     851          600 :       DO iatom = 1, ins_atoms
     852         2334 :          r_cbmc(1:3, iatom) = particles_insert%els(iatom)%r(1:3)
     853              :       END DO
     854              : 
     855              : ! make sure there is no overlap
     856              : 
     857           22 :       IF (loverlap_ins .OR. loverlap_rem) THEN
     858              : ! deallocate some stuff
     859            0 :          CALL mc_molecule_info_destroy(mc_molecule_info_test)
     860            0 :          CALL force_env_release(test_env(insert_box)%force_env)
     861            0 :          CALL force_env_release(test_env(remove_box)%force_env)
     862            0 :          DEALLOCATE (insert_coords)
     863            0 :          DEALLOCATE (remove_coords)
     864            0 :          DEALLOCATE (r_cbmc)
     865            0 :          DEALLOCATE (cbmc_energies)
     866            0 :          DEALLOCATE (oldsys)
     867            0 :          DEALLOCATE (particles_old)
     868            0 :          DEALLOCATE (test_env)
     869            0 :          CALL timestop(handle)
     870            0 :          RETURN
     871              :       END IF
     872              : 
     873              : ! broadcast the chosen coordinates to all processors
     874              : 
     875              :       CALL force_env_get(test_env(insert_box)%force_env, &
     876           22 :                          subsys=insert_sys)
     877              :       CALL cp_subsys_get(insert_sys, &
     878           22 :                          particles=particles_insert)
     879              : 
     880          600 :       DO iatom = 1, ins_atoms
     881              :          particles_insert%els(iatom)%r(1:3) = &
     882         2334 :             r_cbmc(1:3, iatom)
     883              :       END DO
     884              : 
     885              : ! if we made it this far, we have no overlaps
     886              :       moves(molecule_type, insert_box)%moves%grown = &
     887           22 :          moves(molecule_type, insert_box)%moves%grown + 1
     888              : 
     889              : ! if we're biasing, we need to make environments with the non-biasing
     890              : ! potentials, and calculate the energies
     891           22 :       IF (lbias) THEN
     892              : 
     893           66 :          ALLOCATE (test_env_bias(1:2))
     894              : 
     895              : ! first, the environment to which we added a molecule
     896           22 :          CALL get_mc_par(mc_par(insert_box)%mc_par, mc_input_file=mc_input_file)
     897           22 :          IF (ionode) CALL mc_make_dat_file_new(r_cbmc(:, :), atom_names_insert(:), ins_atoms, &
     898              :                                                abc_insert(:), dat_file(insert_box), nchains_test(:, insert_box), &
     899           11 :                                                mc_input_file)
     900           22 :          test_env_bias(insert_box)%force_env => test_env(insert_box)%force_env
     901           22 :          NULLIFY (test_env(insert_box)%force_env)
     902              :          CALL mc_create_force_env(test_env(insert_box)%force_env, input_declaration, &
     903           22 :                                   para_env, dat_file(insert_box))
     904              : 
     905              :          CALL force_env_calc_energy_force(test_env(insert_box)%force_env, &
     906           22 :                                           calc_force=.FALSE.)
     907              :          CALL force_env_get(test_env(insert_box)%force_env, &
     908           22 :                             potential_energy=new_energy(insert_box))
     909              : 
     910              : ! now the environment that has one less molecule
     911           66 :          IF (SUM(nchains_test(:, remove_box)) == 0) THEN
     912            0 :             CALL get_mc_par(mc_par(remove_box)%mc_par, mc_input_file=mc_input_file)
     913            0 :             IF (ionode) CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
     914              :                                                   abc_remove(:), dat_file(remove_box), nchains_test(:, remove_box), &
     915            0 :                                                   mc_input_file)
     916            0 :             test_env_bias(remove_box)%force_env => test_env(remove_box)%force_env
     917            0 :             NULLIFY (test_env(remove_box)%force_env)
     918              :             CALL mc_create_force_env(test_env(remove_box)%force_env, input_declaration, &
     919            0 :                                      para_env, dat_file(remove_box))
     920            0 :             new_energy(remove_box) = 0.0E0_dp
     921            0 :             bias_energy_new(remove_box) = 0.0E0_dp
     922              :          ELSE
     923           22 :             CALL get_mc_par(mc_par(remove_box)%mc_par, mc_input_file=mc_input_file)
     924           22 :             IF (ionode) CALL mc_make_dat_file_new(remove_coords(:, :), atom_names_remove(:), rem_atoms, &
     925              :                                                   abc_remove(:), dat_file(remove_box), nchains_test(:, remove_box), &
     926           11 :                                                   mc_input_file)
     927           22 :             test_env_bias(remove_box)%force_env => test_env(remove_box)%force_env
     928           22 :             NULLIFY (test_env(remove_box)%force_env)
     929              :             CALL mc_create_force_env(test_env(remove_box)%force_env, input_declaration, &
     930           22 :                                      para_env, dat_file(remove_box))
     931              :             CALL force_env_calc_energy_force(test_env(remove_box)%force_env, &
     932           22 :                                              calc_force=.FALSE.)
     933              :             CALL force_env_get(test_env(remove_box)%force_env, &
     934           22 :                                potential_energy=new_energy(remove_box))
     935              :             CALL force_env_calc_energy_force(test_env_bias(remove_box)%force_env, &
     936           22 :                                              calc_force=.FALSE.)
     937              :             CALL force_env_get(test_env_bias(remove_box)%force_env, &
     938           22 :                                potential_energy=bias_energy_new(remove_box))
     939              :          END IF
     940              :       ELSE
     941            0 :          IF (SUM(nchains_test(:, remove_box)) == 0) THEN
     942            0 :             new_energy(remove_box) = 0.0E0_dp
     943              :          ELSE
     944              :             CALL force_env_calc_energy_force(test_env(remove_box)%force_env, &
     945            0 :                                              calc_force=.FALSE.)
     946              :             CALL force_env_get(test_env(remove_box)%force_env, &
     947            0 :                                potential_energy=new_energy(remove_box))
     948              :          END IF
     949              :       END IF
     950              : 
     951              : ! now we need to figure out the rosenbluth weight for the old configuration...
     952              : ! we wait until now to do that because we need the energy of the box that
     953              : ! has had a molecule removed...notice we use the environment that has not
     954              : ! had a molecule removed for the CBMC configurations, and therefore nchains
     955              : ! and mol_type instead of nchains_test and mol_type_test
     956           22 :       start_mol = 1
     957           32 :       DO jbox = 1, remove_box - 1
     958           52 :          start_mol = start_mol + SUM(nchains(:, jbox))
     959              :       END DO
     960           66 :       end_mol = start_mol + SUM(nchains(:, remove_box)) - 1
     961           22 :       IF (lbias) THEN
     962              :          CALL generate_cbmc_swap_config(bias_env(remove_box)%force_env, &
     963              :                                         BETA, max_val, min_val, exp_max_val, &
     964              :                                         exp_min_val, nswapmoves, weight_old, start_atom_rem, nunits_tot(remove_box), &
     965              :                                         nunits(:), nunits(molecule_type), mass(:, molecule_type), loverlap_rem, rdum, &
     966              :                                         bias_energy_new(remove_box), ionode, .TRUE., mol_type(start_mol:end_mol), &
     967           22 :                                         nchains(:, remove_box), source, group, rng_stream)
     968              :       ELSE
     969              :          CALL generate_cbmc_swap_config(force_env(remove_box)%force_env, &
     970              :                                         BETA, max_val, min_val, exp_max_val, &
     971              :                                         exp_min_val, nswapmoves, weight_old, start_atom_rem, nunits_tot(remove_box), &
     972              :                                         nunits(:), nunits(molecule_type), mass(:, molecule_type), loverlap_rem, rdum, &
     973              :                                         new_energy(remove_box), ionode, .TRUE., mol_type(start_mol:end_mol), &
     974            0 :                                         nchains(:, remove_box), source, group, rng_stream)
     975              :       END IF
     976              : 
     977              : ! figure out the prefactor to the boltzmann weight in the acceptance
     978              : ! rule, based on numbers of particles and volumes
     979              : 
     980              :       prefactor = REAL(nchains(molecule_type, remove_box), dp)/ &
     981              :                   REAL(nchains(molecule_type, insert_box) + 1, dp)* &
     982           22 :                   vol_insert/vol_remove
     983              : 
     984           22 :       IF (lbias) THEN
     985              : 
     986              :          del_quickstep_energy = (-BETA)*(new_energy(insert_box) - &
     987              :                                          old_energy(insert_box) + new_energy(remove_box) - &
     988              :                                          old_energy(remove_box) - (bias_energy_new(insert_box) + &
     989              :                                                                    bias_energy_new(remove_box) - bias_energy_old(insert_box) &
     990           22 :                                                                    - bias_energy_old(remove_box)))
     991              : 
     992           22 :          IF (del_quickstep_energy .GT. exp_max_val) THEN
     993            0 :             del_quickstep_energy = max_val
     994           22 :          ELSEIF (del_quickstep_energy .LT. exp_min_val) THEN
     995            0 :             del_quickstep_energy = min_val
     996              :          ELSE
     997           22 :             del_quickstep_energy = EXP(del_quickstep_energy)
     998              :          END IF
     999              :          w = prefactor*del_quickstep_energy*weight_new/weight_old &
    1000           22 :              *EXP(BETA*(eta_remove(molecule_type) - eta_insert(molecule_type)))
    1001              : 
    1002              :       ELSE
    1003              :          w = prefactor*weight_new/weight_old &
    1004            0 :              *EXP(BETA*(eta_remove(molecule_type) - eta_insert(molecule_type)))
    1005              : 
    1006              :       END IF
    1007              : 
    1008              : ! check if the move is accepted
    1009           22 :       IF (w .GE. 1.0E0_dp) THEN
    1010            8 :          rand = 0.0E0_dp
    1011              :       ELSE
    1012           14 :          IF (ionode) rand = rng_stream%next()
    1013           14 :          CALL group%bcast(rand, source)
    1014              :       END IF
    1015              : 
    1016              : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1017           22 :       IF (rand .LT. w) THEN
    1018              : 
    1019              : ! accept the move
    1020              : 
    1021              : !     accept the move
    1022              :          moves(molecule_type, insert_box)%moves%swap%successes = &
    1023           14 :             moves(molecule_type, insert_box)%moves%swap%successes + 1
    1024              : 
    1025              : ! we need to compensate for the fact that we take the difference in
    1026              : ! generate_cbmc_config to keep the exponetials small
    1027           14 :          IF (.NOT. lbias) THEN
    1028              :             new_energy(insert_box) = new_energy(insert_box) + &
    1029            0 :                                      old_energy(insert_box)
    1030              :          END IF
    1031              : 
    1032           42 :          DO ibox = 1, 2
    1033              : ! update energies
    1034              :             energy_check(ibox) = energy_check(ibox) + (new_energy(ibox) - &
    1035           28 :                                                        old_energy(ibox))
    1036           28 :             old_energy(ibox) = new_energy(ibox)
    1037              : ! if we're biasing the update the biasing energy
    1038           42 :             IF (lbias) THEN
    1039           28 :                last_bias_energy(ibox) = bias_energy_new(ibox)
    1040           28 :                bias_energy_old(ibox) = bias_energy_new(ibox)
    1041              :             END IF
    1042              : 
    1043              :          END DO
    1044              : 
    1045              : ! change particle numbers...basically destroy the old mc_molecule_info and attach
    1046              : ! the new stuff to the mc_pars
    1047              : ! figure out the molecule information for the new environments
    1048           14 :          CALL mc_molecule_info_destroy(mc_molecule_info)
    1049           14 :          CALL set_mc_par(mc_par(insert_box)%mc_par, mc_molecule_info=mc_molecule_info_test)
    1050           14 :          CALL set_mc_par(mc_par(remove_box)%mc_par, mc_molecule_info=mc_molecule_info_test)
    1051              : 
    1052              : ! update coordinates
    1053              :          CALL force_env_get(test_env(insert_box)%force_env, &
    1054           14 :                             subsys=insert_sys)
    1055              :          CALL cp_subsys_get(insert_sys, &
    1056           14 :                             particles=particles_insert)
    1057          376 :          DO ipart = 1, ins_atoms
    1058         1462 :             r_old(1:3, ipart, insert_box) = particles_insert%els(ipart)%r(1:3)
    1059              :          END DO
    1060              :          CALL force_env_get(test_env(remove_box)%force_env, &
    1061           14 :                             subsys=remove_sys)
    1062              :          CALL cp_subsys_get(remove_sys, &
    1063           14 :                             particles=particles_remove)
    1064          352 :          DO ipart = 1, rem_atoms
    1065         1366 :             r_old(1:3, ipart, remove_box) = particles_remove%els(ipart)%r(1:3)
    1066              :          END DO
    1067              : 
    1068              :          ! insertion box
    1069           14 :          CALL force_env_release(force_env(insert_box)%force_env)
    1070           14 :          force_env(insert_box)%force_env => test_env(insert_box)%force_env
    1071              : 
    1072              :          ! removal box
    1073           14 :          CALL force_env_release(force_env(remove_box)%force_env)
    1074           14 :          force_env(remove_box)%force_env => test_env(remove_box)%force_env
    1075              : 
    1076              : ! if we're biasing, update the bias_env
    1077           14 :          IF (lbias) THEN
    1078           14 :             CALL force_env_release(bias_env(insert_box)%force_env)
    1079           14 :             bias_env(insert_box)%force_env => test_env_bias(insert_box)%force_env
    1080           14 :             CALL force_env_release(bias_env(remove_box)%force_env)
    1081           14 :             bias_env(remove_box)%force_env => test_env_bias(remove_box)%force_env
    1082           14 :             DEALLOCATE (test_env_bias)
    1083              :          END IF
    1084              : 
    1085              :       ELSE
    1086              : 
    1087              : ! reject the move
    1088            8 :          CALL mc_molecule_info_destroy(mc_molecule_info_test)
    1089            8 :          CALL force_env_release(test_env(insert_box)%force_env)
    1090            8 :          CALL force_env_release(test_env(remove_box)%force_env)
    1091            8 :          IF (lbias) THEN
    1092            8 :             CALL force_env_release(test_env_bias(insert_box)%force_env)
    1093            8 :             CALL force_env_release(test_env_bias(remove_box)%force_env)
    1094            8 :             DEALLOCATE (test_env_bias)
    1095              :          END IF
    1096              :       END IF
    1097              : 
    1098              : ! deallocate some stuff
    1099           22 :       DEALLOCATE (insert_coords)
    1100           22 :       DEALLOCATE (remove_coords)
    1101           22 :       DEALLOCATE (test_env)
    1102           22 :       DEALLOCATE (cbmc_energies)
    1103           22 :       DEALLOCATE (r_cbmc)
    1104           22 :       DEALLOCATE (oldsys)
    1105           22 :       DEALLOCATE (particles_old)
    1106              : 
    1107              : ! end the timing
    1108           22 :       CALL timestop(handle)
    1109              : 
    1110           66 :    END SUBROUTINE mc_ge_swap_move
    1111              : 
    1112              : ! **************************************************************************************************
    1113              : !> \brief performs a Monte Carlo move that alters the volume of the simulation boxes,
    1114              : !>      keeping the total volume of the two boxes the same
    1115              : !> \param mc_par the mc parameters for the force env
    1116              : !> \param force_env the force environments used in the move
    1117              : !> \param moves the structure that keeps track of how many moves have been
    1118              : !>               accepted/rejected
    1119              : !> \param move_updates the structure that keeps track of how many moves have
    1120              : !>               been accepted/rejected since the last time the displacements
    1121              : !>               were updated
    1122              : !> \param nnstep the total number of Monte Carlo moves already performed
    1123              : !> \param old_energy the energy of the last accepted move involving an
    1124              : !>                    unbiased potential calculation
    1125              : !> \param energy_check the running total of how much the energy has changed
    1126              : !>                      since the initial configuration
    1127              : !> \param r_old the coordinates of the last accepted move involving a
    1128              : !>               Quickstep calculation
    1129              : !> \param rng_stream the stream we pull random numbers from
    1130              : !> \author MJM
    1131              : ! **************************************************************************************************
    1132           24 :    SUBROUTINE mc_ge_volume_move(mc_par, force_env, moves, move_updates, &
    1133           24 :                                 nnstep, old_energy, energy_check, r_old, rng_stream)
    1134              : 
    1135              :       TYPE(mc_simulation_parameters_p_type), &
    1136              :          DIMENSION(:), POINTER                           :: mc_par
    1137              :       TYPE(force_env_p_type), DIMENSION(:), POINTER      :: force_env
    1138              :       TYPE(mc_moves_p_type), DIMENSION(:, :), POINTER    :: moves, move_updates
    1139              :       INTEGER, INTENT(IN)                                :: nnstep
    1140              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: old_energy, energy_check
    1141              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: r_old
    1142              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
    1143              : 
    1144              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'mc_ge_volume_move'
    1145              : 
    1146              :       CHARACTER(LEN=200)                                 :: fft_lib
    1147              :       CHARACTER(LEN=40), DIMENSION(1:2)                  :: dat_file
    1148              :       INTEGER :: cl, end_atom, end_mol, handle, iatom, ibox, imolecule, iside, j, jatom, jbox, &
    1149              :          max_atoms, molecule_index, molecule_type, print_level, source, start_atom, start_mol
    1150           24 :       INTEGER, DIMENSION(:), POINTER                     :: mol_type, nunits, nunits_tot
    1151           24 :       INTEGER, DIMENSION(:, :), POINTER                  :: nchains
    1152              :       LOGICAL                                            :: ionode
    1153           24 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: loverlap
    1154              :       LOGICAL, DIMENSION(1:2)                            :: lempty
    1155           24 :       REAL(dp), DIMENSION(:, :), POINTER                 :: mass
    1156              :       REAL(KIND=dp)                                      :: BETA, prefactor, rand, rmvolume, &
    1157              :                                                             vol_dis, w
    1158           24 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: r
    1159              :       REAL(KIND=dp), DIMENSION(1:2)                      :: new_energy, volume_new, volume_old
    1160              :       REAL(KIND=dp), DIMENSION(1:3)                      :: center_of_mass, center_of_mass_new, diff
    1161              :       REAL(KIND=dp), DIMENSION(1:3, 1:2)                 :: abc, new_cell_length, old_cell_length
    1162              :       REAL(KIND=dp), DIMENSION(1:3, 1:3, 1:2)            :: hmat_test
    1163           24 :       TYPE(cell_p_type), DIMENSION(:), POINTER           :: cell, cell_old, cell_test
    1164           24 :       TYPE(cp_subsys_p_type), DIMENSION(:), POINTER      :: oldsys
    1165              :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1166              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
    1167              :       TYPE(mp_comm_type)                                 :: group
    1168           24 :       TYPE(particle_list_p_type), DIMENSION(:), POINTER  :: particles_old
    1169              : 
    1170              : ! begin the timing of the subroutine
    1171              : 
    1172           24 :       CALL timeset(routineN, handle)
    1173              : 
    1174              : ! nullify some pointers
    1175           24 :       NULLIFY (particles_old, cell, oldsys, cell_old, cell_test, subsys)
    1176              : 
    1177              : ! get some data from mc_par
    1178              :       CALL get_mc_par(mc_par(1)%mc_par, ionode=ionode, source=source, &
    1179              :                       group=group, dat_file=dat_file(1), rmvolume=rmvolume, &
    1180              :                       BETA=BETA, cl=cl, fft_lib=fft_lib, &
    1181           24 :                       mc_molecule_info=mc_molecule_info)
    1182              :       CALL get_mc_molecule_info(mc_molecule_info, nunits_tot=nunits_tot, &
    1183           24 :                                 mass=mass, nchains=nchains, nunits=nunits, mol_type=mol_type)
    1184              : 
    1185           24 :       print_level = 1
    1186           24 :       CALL get_mc_par(mc_par(2)%mc_par, dat_file=dat_file(2))
    1187              : 
    1188              : ! allocate some stuff
    1189           24 :       max_atoms = MAX(nunits_tot(1), nunits_tot(2))
    1190           96 :       ALLOCATE (r(1:3, max_atoms, 1:2))
    1191           72 :       ALLOCATE (oldsys(1:2))
    1192           72 :       ALLOCATE (particles_old(1:2))
    1193           72 :       ALLOCATE (cell(1:2))
    1194           72 :       ALLOCATE (cell_test(1:2))
    1195           72 :       ALLOCATE (cell_old(1:2))
    1196           24 :       ALLOCATE (loverlap(1:2))
    1197              : 
    1198              : ! check for empty boxes...need to be careful because we can't build
    1199              : ! a force_env with no particles
    1200           72 :       DO ibox = 1, 2
    1201           48 :          lempty(ibox) = .FALSE.
    1202          168 :          IF (SUM(nchains(:, ibox)) == 0) THEN
    1203            0 :             lempty(ibox) = .TRUE.
    1204              :          END IF
    1205              :       END DO
    1206              : 
    1207              : ! record the attempt
    1208           72 :       DO ibox = 1, 2
    1209              :          moves(1, ibox)%moves%volume%attempts = &
    1210           48 :             moves(1, ibox)%moves%volume%attempts + 1
    1211              :          move_updates(1, ibox)%moves%volume%attempts = &
    1212           72 :             move_updates(1, ibox)%moves%volume%attempts + 1
    1213              :       END DO
    1214              : 
    1215              : ! now let's grab the cell length and particle positions
    1216           72 :       DO ibox = 1, 2
    1217              :          CALL force_env_get(force_env(ibox)%force_env, &
    1218           48 :                             subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
    1219           48 :          CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
    1220           48 :          NULLIFY (cell_old(ibox)%cell)
    1221           48 :          CALL cell_create(cell_old(ibox)%cell)
    1222           48 :          CALL cell_clone(cell(ibox)%cell, cell_old(ibox)%cell, tag="CELL_OLD")
    1223              :          CALL cp_subsys_get(oldsys(ibox)%subsys, &
    1224           48 :                             particles=particles_old(ibox)%list)
    1225              : 
    1226              : ! find the old cell length
    1227          216 :          old_cell_length(1:3, ibox) = abc(1:3, ibox)
    1228              : 
    1229              :       END DO
    1230              : 
    1231           72 :       DO ibox = 1, 2
    1232              : 
    1233              : ! save the old coordinates
    1234         1272 :          DO iatom = 1, nunits_tot(ibox)
    1235         4848 :             r(1:3, iatom, ibox) = particles_old(ibox)%list%els(iatom)%r(1:3)
    1236              :          END DO
    1237              : 
    1238              :       END DO
    1239              : 
    1240              : ! call a random number to figure out how far we're moving
    1241           24 :       IF (ionode) rand = rng_stream%next()
    1242           24 :       CALL group%bcast(rand, source)
    1243              : 
    1244           24 :       vol_dis = rmvolume*(rand - 0.5E0_dp)*2.0E0_dp
    1245              : 
    1246              : ! add to one box, subtract from the other
    1247           24 :       IF (old_cell_length(1, 1)*old_cell_length(2, 1)* &
    1248              :           old_cell_length(3, 1) + vol_dis .LE. (3.0E0_dp/angstrom)**3) &
    1249            0 :          CPABORT('GE_volume moves are trying to make box 1 smaller than 3')
    1250           24 :       IF (old_cell_length(1, 2)*old_cell_length(2, 2)* &
    1251              :           old_cell_length(3, 2) + vol_dis .LE. (3.0E0_dp/angstrom)**3) &
    1252            0 :          CPABORT('GE_volume moves are trying to make box 2 smaller than 3')
    1253              : 
    1254           96 :       DO iside = 1, 3
    1255              :          new_cell_length(iside, 1) = (old_cell_length(1, 1)**3 + &
    1256           72 :                                       vol_dis)**(1.0E0_dp/3.0E0_dp)
    1257              :          new_cell_length(iside, 2) = (old_cell_length(1, 2)**3 - &
    1258           96 :                                       vol_dis)**(1.0E0_dp/3.0E0_dp)
    1259              :       END DO
    1260              : 
    1261              : ! now we need to make the new cells
    1262           72 :       DO ibox = 1, 2
    1263          624 :          hmat_test(:, :, ibox) = 0.0e0_dp
    1264           48 :          hmat_test(1, 1, ibox) = new_cell_length(1, ibox)
    1265           48 :          hmat_test(2, 2, ibox) = new_cell_length(2, ibox)
    1266           48 :          hmat_test(3, 3, ibox) = new_cell_length(3, ibox)
    1267           48 :          NULLIFY (cell_test(ibox)%cell)
    1268              :          CALL cell_create(cell_test(ibox)%cell, hmat=hmat_test(:, :, ibox), &
    1269           48 :                           periodic=cell(ibox)%cell%perd)
    1270           48 :          CALL force_env_get(force_env(ibox)%force_env, subsys=subsys)
    1271           72 :          CALL cp_subsys_set(subsys, cell=cell_test(ibox)%cell)
    1272              :       END DO
    1273              : 
    1274           72 :       DO ibox = 1, 2
    1275              : 
    1276              : ! save the coords
    1277         1248 :          DO iatom = 1, nunits_tot(ibox)
    1278         4848 :             r(1:3, iatom, ibox) = particles_old(ibox)%list%els(iatom)%r(1:3)
    1279              :          END DO
    1280              : 
    1281              : ! now we need to scale the coordinates of all the molecules by the
    1282              : ! center of mass
    1283           72 :          start_atom = 1
    1284              :          molecule_index = 1
    1285           72 :          DO jbox = 1, ibox - 1
    1286              :             IF (jbox == ibox) EXIT
    1287          120 :             molecule_index = molecule_index + SUM(nchains(:, jbox))
    1288              :          END DO
    1289          720 :          DO imolecule = 1, SUM(nchains(:, ibox))
    1290          576 :             molecule_type = mol_type(imolecule + molecule_index - 1)
    1291          576 :             IF (imolecule .NE. 1) THEN
    1292          528 :                start_atom = start_atom + nunits(mol_type(imolecule + molecule_index - 2))
    1293              :             END IF
    1294          576 :             end_atom = start_atom + nunits(molecule_type) - 1
    1295              : 
    1296              : ! now find the center of mass
    1297              :             CALL get_center_of_mass(r(:, start_atom:end_atom, ibox), &
    1298          576 :                                     nunits(molecule_type), center_of_mass(:), mass(:, molecule_type))
    1299              : 
    1300              : ! scale the center of mass and determine the vector that points from the
    1301              : !    old COM to the new one
    1302              :             center_of_mass_new(1:3) = center_of_mass(1:3)* &
    1303         2304 :                                       new_cell_length(1:3, ibox)/old_cell_length(1:3, ibox)
    1304         2352 :             DO j = 1, 3
    1305         1728 :                diff(j) = center_of_mass_new(j) - center_of_mass(j)
    1306              : ! now change the particle positions
    1307         5904 :                DO jatom = start_atom, end_atom
    1308              :                   particles_old(ibox)%list%els(jatom)%r(j) = &
    1309         5328 :                      particles_old(ibox)%list%els(jatom)%r(j) + diff(j)
    1310              :                END DO
    1311              : 
    1312              :             END DO
    1313              :          END DO
    1314              : 
    1315              : ! check for any overlaps we might have
    1316              :          start_mol = 1
    1317           72 :          DO jbox = 1, ibox - 1
    1318          120 :             start_mol = start_mol + SUM(nchains(:, jbox))
    1319              :          END DO
    1320          144 :          end_mol = start_mol + SUM(nchains(:, ibox)) - 1
    1321              :          CALL check_for_overlap(force_env(ibox)%force_env, &
    1322              :                                 nchains(:, ibox), nunits, loverlap(ibox), mol_type(start_mol:end_mol), &
    1323           72 :                                 cell_length=new_cell_length(:, ibox))
    1324              : 
    1325              :       END DO
    1326              : 
    1327              : ! determine the overall energy difference
    1328              : 
    1329           72 :       DO ibox = 1, 2
    1330           48 :          IF (loverlap(ibox)) CYCLE
    1331              : ! remake the force environment and calculate the energy
    1332           72 :          IF (lempty(ibox)) THEN
    1333            0 :             new_energy(ibox) = 0.0E0_dp
    1334              :          ELSE
    1335              : 
    1336              :             CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
    1337           48 :                                              calc_force=.FALSE.)
    1338              :             CALL force_env_get(force_env(ibox)%force_env, &
    1339           48 :                                potential_energy=new_energy(ibox))
    1340              : 
    1341              :          END IF
    1342              :       END DO
    1343              : 
    1344              : ! accept or reject the move
    1345           72 :       DO ibox = 1, 2
    1346              :          volume_new(ibox) = new_cell_length(1, ibox)* &
    1347           48 :                             new_cell_length(2, ibox)*new_cell_length(3, ibox)
    1348              :          volume_old(ibox) = old_cell_length(1, ibox)* &
    1349           72 :                             old_cell_length(2, ibox)*old_cell_length(3, ibox)
    1350              :       END DO
    1351              :       prefactor = (volume_new(1)/volume_old(1))**(SUM(nchains(:, 1)))* &
    1352          120 :                   (volume_new(2)/volume_old(2))**(SUM(nchains(:, 2)))
    1353              : 
    1354           24 :       IF (loverlap(1) .OR. loverlap(2)) THEN
    1355            0 :          w = 0.0E0_dp
    1356              :       ELSE
    1357              :          w = prefactor*EXP(-BETA* &
    1358              :                            (new_energy(1) + new_energy(2) - &
    1359           24 :                             old_energy(1) - old_energy(2)))
    1360              : 
    1361              :       END IF
    1362              : 
    1363           24 :       IF (w .GE. 1.0E0_dp) THEN
    1364            6 :          w = 1.0E0_dp
    1365            6 :          rand = 0.0E0_dp
    1366              :       ELSE
    1367           18 :          IF (ionode) rand = rng_stream%next()
    1368           18 :          CALL group%bcast(rand, source)
    1369              :       END IF
    1370              : 
    1371           24 :       IF (rand .LT. w) THEN
    1372              : 
    1373              : ! write cell length, volume, density, and trial displacement to a file
    1374           18 :          IF (ionode) THEN
    1375              : 
    1376            9 :             WRITE (cl, *) nnstep, new_cell_length(1, 1)* &
    1377            9 :                angstrom, vol_dis*(angstrom)**3, new_cell_length(1, 2)* &
    1378           18 :                angstrom
    1379            9 :             WRITE (cl, *) nnstep, new_energy(1), &
    1380           18 :                old_energy(1), new_energy(2), old_energy(2)
    1381            9 :             WRITE (cl, *) prefactor, w
    1382              :          END IF
    1383              : 
    1384           54 :          DO ibox = 1, 2
    1385              : ! accept the move
    1386              :             moves(1, ibox)%moves%volume%successes = &
    1387           36 :                moves(1, ibox)%moves%volume%successes + 1
    1388              :             move_updates(1, ibox)%moves%volume%successes = &
    1389           36 :                move_updates(1, ibox)%moves%volume%successes + 1
    1390              : 
    1391              : ! update energies
    1392              :             energy_check(ibox) = energy_check(ibox) + (new_energy(ibox) - &
    1393           36 :                                                        old_energy(ibox))
    1394           36 :             old_energy(ibox) = new_energy(ibox)
    1395              : 
    1396              : ! and the new "old" coordinates
    1397          954 :             DO iatom = 1, nunits_tot(ibox)
    1398              :                r_old(1:3, iatom, ibox) = &
    1399         3636 :                   particles_old(ibox)%list%els(iatom)%r(1:3)
    1400              :             END DO
    1401              : 
    1402              :          END DO
    1403              :       ELSE
    1404              : 
    1405              : ! reject the move
    1406              : ! write cell length, volume, density, and trial displacement to a file
    1407            6 :          IF (ionode) THEN
    1408              : 
    1409            3 :             WRITE (cl, *) nnstep, new_cell_length(1, 1)* &
    1410            3 :                angstrom, vol_dis*(angstrom)**3, new_cell_length(1, 2)* &
    1411            6 :                angstrom
    1412            3 :             WRITE (cl, *) nnstep, new_energy(1), &
    1413            6 :                old_energy(1), new_energy(2), old_energy(2)
    1414            3 :             WRITE (cl, *) prefactor, w
    1415              : 
    1416              :          END IF
    1417              : 
    1418              : ! reset the cell and particle positions
    1419           18 :          DO ibox = 1, 2
    1420           12 :             CALL force_env_get(force_env(ibox)%force_env, subsys=subsys)
    1421           12 :             CALL cp_subsys_set(subsys, cell=cell_old(ibox)%cell)
    1422          318 :             DO iatom = 1, nunits_tot(ibox)
    1423         1212 :                particles_old(ibox)%list%els(iatom)%r(1:3) = r_old(1:3, iatom, ibox)
    1424              :             END DO
    1425              :          END DO
    1426              : 
    1427              :       END IF
    1428              : 
    1429              : ! free up some memory
    1430           72 :       DO ibox = 1, 2
    1431           48 :          CALL cell_release(cell_test(ibox)%cell)
    1432           72 :          CALL cell_release(cell_old(ibox)%cell)
    1433              :       END DO
    1434           24 :       DEALLOCATE (r)
    1435           24 :       DEALLOCATE (oldsys)
    1436           24 :       DEALLOCATE (particles_old)
    1437           24 :       DEALLOCATE (cell)
    1438           24 :       DEALLOCATE (cell_old)
    1439           24 :       DEALLOCATE (cell_test)
    1440           24 :       DEALLOCATE (loverlap)
    1441              : 
    1442              : ! end the timing
    1443           24 :       CALL timestop(handle)
    1444              : 
    1445           24 :    END SUBROUTINE mc_ge_volume_move
    1446              : 
    1447              : END MODULE mc_ge_moves
    1448              : 
        

Generated by: LCOV version 2.0-1