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

            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 Used to run the bulk of the MC simulation, doing things like
      10              : !>      choosing move types and writing data to files
      11              : !> \author Matthew J. McGrath  (09.26.2003)
      12              : !>
      13              : !>    REVISIONS
      14              : !>      09.10.05  MJM combined the two subroutines in this module into one
      15              : ! **************************************************************************************************
      16              : MODULE mc_ensembles
      17              :    USE cell_types,                      ONLY: cell_p_type,&
      18              :                                               get_cell
      19              :    USE cp_external_control,             ONLY: external_control
      20              :    USE cp_files,                        ONLY: close_file,&
      21              :                                               open_file
      22              :    USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
      23              :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      24              :                                               cp_subsys_p_type,&
      25              :                                               cp_subsys_type
      26              :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      27              :    USE force_env_methods,               ONLY: force_env_calc_energy_force
      28              :    USE force_env_types,                 ONLY: force_env_get,&
      29              :                                               force_env_p_type,&
      30              :                                               force_env_release
      31              :    USE global_types,                    ONLY: global_environment_type
      32              :    USE input_constants,                 ONLY: dump_xmol
      33              :    USE input_section_types,             ONLY: section_type,&
      34              :                                               section_vals_type,&
      35              :                                               section_vals_val_get
      36              :    USE kinds,                           ONLY: default_string_length,&
      37              :                                               dp
      38              :    USE machine,                         ONLY: m_flush
      39              :    USE mathconstants,                   ONLY: pi
      40              :    USE mc_control,                      ONLY: mc_create_bias_force_env,&
      41              :                                               write_mc_restart
      42              :    USE mc_coordinates,                  ONLY: check_for_overlap,&
      43              :                                               create_discrete_array,&
      44              :                                               find_mc_test_molecule,&
      45              :                                               get_center_of_mass,&
      46              :                                               mc_coordinate_fold,&
      47              :                                               rotate_molecule
      48              :    USE mc_environment_types,            ONLY: get_mc_env,&
      49              :                                               mc_environment_p_type,&
      50              :                                               set_mc_env
      51              :    USE mc_ge_moves,                     ONLY: mc_ge_swap_move,&
      52              :                                               mc_ge_volume_move,&
      53              :                                               mc_quickstep_move
      54              :    USE mc_misc,                         ONLY: final_mc_write,&
      55              :                                               mc_averages_create,&
      56              :                                               mc_averages_release
      57              :    USE mc_move_control,                 ONLY: init_mc_moves,&
      58              :                                               mc_move_update,&
      59              :                                               mc_moves_release,&
      60              :                                               write_move_stats
      61              :    USE mc_moves,                        ONLY: mc_avbmc_move,&
      62              :                                               mc_cluster_translation,&
      63              :                                               mc_conformation_change,&
      64              :                                               mc_hmc_move,&
      65              :                                               mc_molecule_rotation,&
      66              :                                               mc_molecule_translation,&
      67              :                                               mc_volume_move
      68              :    USE mc_types,                        ONLY: get_mc_molecule_info,&
      69              :                                               get_mc_par,&
      70              :                                               mc_averages_p_type,&
      71              :                                               mc_input_file_type,&
      72              :                                               mc_molecule_info_type,&
      73              :                                               mc_moves_p_type,&
      74              :                                               mc_simulation_parameters_p_type,&
      75              :                                               set_mc_par
      76              :    USE message_passing,                 ONLY: mp_comm_type,&
      77              :                                               mp_para_env_type
      78              :    USE parallel_rng_types,              ONLY: rng_stream_type
      79              :    USE particle_list_types,             ONLY: particle_list_p_type,&
      80              :                                               particle_list_type
      81              :    USE particle_methods,                ONLY: write_particle_coordinates
      82              :    USE physcon,                         ONLY: angstrom,&
      83              :                                               boltzmann,&
      84              :                                               joule,&
      85              :                                               n_avogadro
      86              : #include "../../base/base_uses.f90"
      87              : 
      88              :    IMPLICIT NONE
      89              : 
      90              :    PRIVATE
      91              : 
      92              : ! *** Global parameters ***
      93              : 
      94              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_ensembles'
      95              :    LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
      96              : 
      97              :    PUBLIC :: mc_run_ensemble, mc_compute_virial
      98              : 
      99              : CONTAINS
     100              : 
     101              : ! **************************************************************************************************
     102              : !> \brief directs the program in running one or two box MC simulations
     103              : !> \param mc_env a pointer that contains all mc_env for all the simulation
     104              : !>          boxes
     105              : !> \param para_env ...
     106              : !> \param globenv the global environment for the simulation
     107              : !> \param input_declaration ...
     108              : !> \param nboxes the number of simulation boxes
     109              : !> \param rng_stream the stream we pull random numbers from
     110              : !>
     111              : !>    Suitable for parallel.
     112              : !> \author MJM
     113              : ! **************************************************************************************************
     114           18 :    SUBROUTINE mc_run_ensemble(mc_env, para_env, globenv, input_declaration, nboxes, rng_stream)
     115              : 
     116              :       TYPE(mc_environment_p_type), DIMENSION(:), POINTER :: mc_env
     117              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     118              :       TYPE(global_environment_type), POINTER             :: globenv
     119              :       TYPE(section_type), POINTER                        :: input_declaration
     120              :       INTEGER, INTENT(IN)                                :: nboxes
     121              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     122              : 
     123              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'mc_run_ensemble'
     124              : 
     125              :       CHARACTER(default_string_length), ALLOCATABLE, &
     126           18 :          DIMENSION(:)                                    :: atom_names_box
     127              :       CHARACTER(default_string_length), &
     128           18 :          DIMENSION(:, :), POINTER                        :: atom_names
     129              :       CHARACTER(LEN=20)                                  :: ensemble
     130              :       CHARACTER(LEN=40)                                  :: cbox, cstep, fft_lib, move_type, &
     131              :                                                             move_type_avbmc
     132           18 :       INTEGER, DIMENSION(:, :), POINTER                  :: nchains
     133           18 :       INTEGER, DIMENSION(:), POINTER                     :: avbmc_atom, mol_type, nchains_box, &
     134           18 :                                                             nunits, nunits_tot
     135           36 :       INTEGER, DIMENSION(1:nboxes)                       :: box_flag, cl, data_unit, diff, istep, &
     136           54 :                                                             move_unit, rm
     137              :       INTEGER, DIMENSION(1:3, 1:2)                       :: discrete_array
     138              :       INTEGER :: atom_number, box_number, cell_unit, com_crd, com_ene, com_mol, end_mol, handle, &
     139              :          ibox, idum, imol_type, imolecule, imove, iparticle, iprint, itype, iunit, iuptrans, &
     140              :          iupvolume, iw, jbox, jdum, molecule_type, molecule_type_swap, molecule_type_target, &
     141              :          nchain_total, nmol_types, nmoves, nnstep, nstart, nstep, source, start_atom, &
     142              :          start_atom_swap, start_atom_target, start_mol
     143              :       CHARACTER(LEN=default_string_length)               :: unit_str
     144           36 :       CHARACTER(LEN=40), DIMENSION(1:nboxes)             :: cell_file, coords_file, data_file, &
     145           36 :                                                             displacement_file, energy_file, &
     146           36 :                                                             molecules_file, moves_file
     147              :       LOGICAL                                            :: ionode, lbias, ldiscrete, lhmc, &
     148              :                                                             lnew_bias_env, loverlap, lreject, &
     149              :                                                             lstop, print_kind, should_stop
     150           18 :       REAL(dp), DIMENSION(:), POINTER                    :: pbias, pmavbmc_mol, pmclus_box, &
     151           18 :                                                             pmhmc_box, pmrot_mol, pmtraion_mol, &
     152           18 :                                                             pmtrans_mol, pmvol_box
     153           18 :       REAL(dp), DIMENSION(:, :), POINTER                 :: conf_prob, mass
     154              :       REAL(KIND=dp)                                      :: discrete_step, pmavbmc, pmcltrans, &
     155              :                                                             pmhmc, pmswap, pmtraion, pmtrans, &
     156              :                                                             pmvolume, rand, test_energy, unit_conv
     157           18 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: r_temp
     158           18 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: r_old
     159           36 :       REAL(KIND=dp), DIMENSION(1:3, 1:nboxes)            :: abc
     160           36 :       REAL(KIND=dp), DIMENSION(1:nboxes)                 :: bias_energy, energy_check, final_energy, &
     161           36 :                                                             initial_energy, last_bias_energy, &
     162           36 :                                                             old_energy
     163           18 :       TYPE(cell_p_type), DIMENSION(:), POINTER           :: cell
     164           18 :       TYPE(cp_subsys_p_type), DIMENSION(:), POINTER      :: oldsys
     165              :       TYPE(cp_subsys_type), POINTER                      :: biassys
     166           18 :       TYPE(force_env_p_type), DIMENSION(:), POINTER      :: bias_env, force_env
     167           18 :       TYPE(mc_averages_p_type), DIMENSION(:), POINTER    :: averages
     168              :       TYPE(mc_input_file_type), POINTER                  :: mc_bias_file
     169              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
     170           18 :       TYPE(mc_moves_p_type), DIMENSION(:), POINTER       :: test_moves
     171           18 :       TYPE(mc_moves_p_type), DIMENSION(:, :), POINTER    :: move_updates, moves
     172              :       TYPE(mc_simulation_parameters_p_type), &
     173           18 :          DIMENSION(:), POINTER                           :: mc_par
     174              :       TYPE(mp_comm_type)                                 :: group
     175           18 :       TYPE(particle_list_p_type), DIMENSION(:), POINTER  :: particles_old
     176              :       TYPE(particle_list_type), POINTER                  :: particles_bias
     177              :       TYPE(section_vals_type), POINTER                   :: root_section
     178              : 
     179           18 :       CALL timeset(routineN, handle)
     180              : 
     181              :       ! nullify some pointers
     182           18 :       NULLIFY (moves, move_updates, test_moves, root_section)
     183              : 
     184              :       ! allocate a whole bunch of stuff based on how many boxes we have
     185           96 :       ALLOCATE (force_env(1:nboxes))
     186           60 :       ALLOCATE (bias_env(1:nboxes))
     187           60 :       ALLOCATE (cell(1:nboxes))
     188           60 :       ALLOCATE (particles_old(1:nboxes))
     189           60 :       ALLOCATE (oldsys(1:nboxes))
     190           60 :       ALLOCATE (averages(1:nboxes))
     191           60 :       ALLOCATE (mc_par(1:nboxes))
     192           36 :       ALLOCATE (pmvol_box(1:nboxes))
     193           36 :       ALLOCATE (pmclus_box(1:nboxes))
     194           36 :       ALLOCATE (pmhmc_box(1:nboxes))
     195              : 
     196           42 :       DO ibox = 1, nboxes
     197              :          CALL get_mc_env(mc_env(ibox)%mc_env, &
     198              :                          mc_par=mc_par(ibox)%mc_par, &
     199           42 :                          force_env=force_env(ibox)%force_env)
     200              :       END DO
     201              : 
     202              :       ! Gather units of measure for output (if available)
     203           18 :       root_section => force_env(1)%force_env%root_section
     204              :       CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%UNIT", &
     205           18 :                                 c_val=unit_str)
     206           18 :       unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
     207              :       CALL section_vals_val_get(root_section, "MOTION%PRINT%TRAJECTORY%PRINT_ATOM_KIND", &
     208           18 :                                 l_val=print_kind)
     209              : 
     210              :       ! get some data out of mc_par
     211              :       CALL get_mc_par(mc_par(1)%mc_par, &
     212              :                       ionode=ionode, source=source, group=group, &
     213              :                       data_file=data_file(1), moves_file=moves_file(1), &
     214              :                       cell_file=cell_file(1), coords_file=coords_file(1), &
     215              :                       energy_file=energy_file(1), displacement_file=displacement_file(1), &
     216              :                       lstop=lstop, nstep=nstep, nstart=nstart, pmvolume=pmvolume, pmhmc=pmhmc, &
     217              :                       molecules_file=molecules_file(1), pmswap=pmswap, nmoves=nmoves, &
     218              :                       pmtraion=pmtraion, pmtrans=pmtrans, pmcltrans=pmcltrans, iuptrans=iuptrans, &
     219              :                       iupvolume=iupvolume, ldiscrete=ldiscrete, pmtraion_mol=pmtraion_mol, &
     220              :                       lbias=lbias, iprint=iprint, pmavbmc_mol=pmavbmc_mol, &
     221              :                       discrete_step=discrete_step, fft_lib=fft_lib, avbmc_atom=avbmc_atom, &
     222              :                       pmavbmc=pmavbmc, pbias=pbias, mc_molecule_info=mc_molecule_info, &
     223              :                       pmrot_mol=pmrot_mol, pmtrans_mol=pmtrans_mol, pmvol_box=pmvol_box(1), &
     224           18 :                       pmclus_box=pmclus_box(1), ensemble=ensemble, pmhmc_box=pmhmc_box(1), lhmc=lhmc)
     225              : 
     226              :       ! get some data from the molecule types
     227              :       CALL get_mc_molecule_info(mc_molecule_info, conf_prob=conf_prob, &
     228              :                                 nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
     229              :                                 mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
     230           18 :                                 atom_names=atom_names, mass=mass)
     231              : 
     232              :       ! allocate some stuff based on the number of molecule types we have
     233          136 :       ALLOCATE (moves(1:nmol_types, 1:nboxes))
     234          118 :       ALLOCATE (move_updates(1:nmol_types, 1:nboxes))
     235              : 
     236           18 :       IF (nboxes .GT. 1) THEN
     237           12 :          DO ibox = 2, nboxes
     238              :             CALL get_mc_par(mc_par(ibox)%mc_par, &
     239              :                             data_file=data_file(ibox), &
     240              :                             moves_file=moves_file(ibox), &
     241              :                             cell_file=cell_file(ibox), coords_file=coords_file(ibox), &
     242              :                             energy_file=energy_file(ibox), &
     243              :                             displacement_file=displacement_file(ibox), &
     244              :                             molecules_file=molecules_file(ibox), pmvol_box=pmvol_box(ibox), &
     245           12 :                             pmclus_box=pmclus_box(ibox), pmhmc_box=pmhmc_box(ibox))
     246              :          END DO
     247              :       END IF
     248              : 
     249              :       ! this is a check we can't do in the input checking
     250           18 :       IF (pmvol_box(nboxes) .LT. 1.0E0_dp) THEN
     251            0 :          CPABORT('The last value of PMVOL_BOX needs to be 1.0')
     252              :       END IF
     253           18 :       IF (pmclus_box(nboxes) .LT. 1.0E0_dp) THEN
     254            0 :          CPABORT('The last value of PMVOL_BOX needs to be 1.0')
     255              :       END IF
     256           18 :       IF (pmhmc_box(nboxes) .LT. 1.0E0_dp) THEN
     257            0 :          CPABORT('The last value of PMHMC_BOX needs to be 1.0')
     258              :       END IF
     259              : 
     260              :       ! allocate the particle positions array for broadcasting
     261           96 :       ALLOCATE (r_old(3, SUM(nunits_tot), 1:nboxes))
     262              : 
     263              :       ! figure out what the default write unit is
     264           18 :       iw = cp_logger_get_default_io_unit()
     265              : 
     266           18 :       IF (iw > 0) THEN
     267            9 :          WRITE (iw, *)
     268            9 :          WRITE (iw, *)
     269            9 :          WRITE (iw, *) 'Beginning the Monte Carlo calculation.'
     270            9 :          WRITE (iw, *)
     271            9 :          WRITE (iw, *)
     272              :       END IF
     273              : 
     274              :       ! initialize running average variables
     275           42 :       energy_check(:) = 0.0E0_dp
     276           42 :       box_flag(:) = 0
     277           42 :       istep(:) = 0
     278              : 
     279           42 :       DO ibox = 1, nboxes
     280              :          ! initialize the moves array, the arrays for updating maximum move
     281              :          ! displacements, and the averages array
     282           64 :          DO itype = 1, nmol_types
     283           40 :             CALL init_mc_moves(moves(itype, ibox)%moves)
     284           64 :             CALL init_mc_moves(move_updates(itype, ibox)%moves)
     285              :          END DO
     286           24 :          CALL mc_averages_create(averages(ibox)%averages)
     287              : 
     288              :          ! find the energy of the initial configuration
     289           64 :          IF (SUM(nchains(:, ibox)) .NE. 0) THEN
     290              :             CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
     291           24 :                                              calc_force=.FALSE.)
     292              :             CALL force_env_get(force_env(ibox)%force_env, &
     293           24 :                                potential_energy=old_energy(ibox))
     294              :          ELSE
     295            0 :             old_energy(ibox) = 0.0E0_dp
     296              :          END IF
     297           24 :          initial_energy(ibox) = old_energy(ibox)
     298              : 
     299              : ! don't care about overlaps if we're only doing HMC
     300              : 
     301           24 :          IF (.NOT. lhmc) THEN
     302              :             ! check for overlaps
     303              :             start_mol = 1
     304           28 :             DO jbox = 1, ibox - 1
     305           40 :                start_mol = start_mol + SUM(nchains(:, jbox))
     306              :             END DO
     307           60 :             end_mol = start_mol + SUM(nchains(:, ibox)) - 1
     308              :             CALL check_for_overlap(force_env(ibox)%force_env, nchains(:, ibox), &
     309           22 :                                    nunits, loverlap, mol_type(start_mol:end_mol))
     310           22 :             IF (loverlap) CPABORT("overlap in an initial configuration")
     311              :          END IF
     312              : 
     313              :          ! get the subsystems and the cell information
     314              :          CALL force_env_get(force_env(ibox)%force_env, &
     315           24 :                             subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
     316           24 :          CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
     317              :          CALL cp_subsys_get(oldsys(ibox)%subsys, &
     318           24 :                             particles=particles_old(ibox)%list)
     319              :          ! record the old coordinates, in case a move is rejected
     320         2656 :          DO iparticle = 1, nunits_tot(ibox)
     321              :             r_old(1:3, iparticle, ibox) = &
     322        10552 :                particles_old(ibox)%list%els(iparticle)%r(1:3)
     323              :          END DO
     324              : 
     325              :          ! find the bias energy of the initial run
     326           24 :          IF (lbias) THEN
     327              :             ! determine the atom names of every particle
     328           42 :             ALLOCATE (atom_names_box(1:nunits_tot(ibox)))
     329              : 
     330           14 :             atom_number = 1
     331          212 :             DO imolecule = 1, SUM(nchains(:, ibox))
     332          538 :                DO iunit = 1, nunits(mol_type(imolecule + start_mol - 1))
     333              :                   atom_names_box(atom_number) = &
     334          354 :                      atom_names(iunit, mol_type(imolecule + start_mol - 1))
     335          524 :                   atom_number = atom_number + 1
     336              :                END DO
     337              :             END DO
     338              : 
     339           14 :             CALL get_mc_par(mc_par(ibox)%mc_par, mc_bias_file=mc_bias_file)
     340           14 :             nchains_box => nchains(:, ibox)
     341              :             CALL mc_create_bias_force_env(bias_env(ibox)%force_env, &
     342              :                                           r_old(:, :, ibox), atom_names_box(:), nunits_tot(ibox), &
     343              :                                           para_env, abc(:, ibox), nchains_box, input_declaration, mc_bias_file, &
     344           14 :                                           ionode)
     345           42 :             IF (SUM(nchains(:, ibox)) .NE. 0) THEN
     346              :                CALL force_env_calc_energy_force(bias_env(ibox)%force_env, &
     347           14 :                                                 calc_force=.FALSE.)
     348              :                CALL force_env_get(bias_env(ibox)%force_env, &
     349           14 :                                   potential_energy=last_bias_energy(ibox))
     350              : 
     351              :             ELSE
     352            0 :                last_bias_energy(ibox) = 0.0E0_dp
     353              :             END IF
     354           14 :             bias_energy(ibox) = last_bias_energy(ibox)
     355           14 :             DEALLOCATE (atom_names_box)
     356              :          END IF
     357           42 :          lnew_bias_env = .FALSE.
     358              : 
     359              :       END DO
     360              : 
     361              :       ! back to seriel for a bunch of I/O stuff
     362           18 :       IF (ionode) THEN
     363              : 
     364              :          ! record the combined energies,coordinates, and cell lengths
     365              :          CALL open_file(file_name='mc_cell_length', &
     366              :                         unit_number=cell_unit, file_position='APPEND', &
     367            9 :                         file_action='WRITE', file_status='UNKNOWN')
     368              :          CALL open_file(file_name='mc_energies', &
     369              :                         unit_number=com_ene, file_position='APPEND', &
     370            9 :                         file_action='WRITE', file_status='UNKNOWN')
     371              :          CALL open_file(file_name='mc_coordinates', &
     372              :                         unit_number=com_crd, file_position='APPEND', &
     373            9 :                         file_action='WRITE', file_status='UNKNOWN')
     374              :          CALL open_file(file_name='mc_molecules', &
     375              :                         unit_number=com_mol, file_position='APPEND', &
     376            9 :                         file_action='WRITE', file_status='UNKNOWN')
     377            9 :          WRITE (com_ene, *) 'Initial Energies:       ', &
     378           18 :             old_energy(1:nboxes)
     379           21 :          DO ibox = 1, nboxes
     380           12 :             WRITE (com_mol, *) 'Initial Molecules:       ', &
     381           33 :                nchains(:, ibox)
     382              :          END DO
     383           21 :          DO ibox = 1, nboxes
     384           12 :             WRITE (cell_unit, *) 'Initial: ', &
     385           60 :                abc(1:3, ibox)*angstrom
     386           12 :             WRITE (cbox, '(I4)') ibox
     387              :             CALL open_file(file_name='energy_differences_box'// &
     388              :                            TRIM(ADJUSTL(cbox)), &
     389              :                            unit_number=diff(ibox), file_position='APPEND', &
     390           12 :                            file_action='WRITE', file_status='UNKNOWN')
     391           32 :             IF (SUM(nchains(:, ibox)) == 0) THEN
     392            0 :                WRITE (com_crd, *) ' 0'
     393            0 :                WRITE (com_crd, *) 'INITIAL BOX '//TRIM(ADJUSTL(cbox))
     394              :             ELSE
     395              :                CALL write_particle_coordinates(particles_old(ibox)%list%els, &
     396              :                                                com_crd, dump_xmol, 'POS', 'INITIAL BOX '//TRIM(ADJUSTL(cbox)), &
     397           12 :                                                unit_conv=unit_conv, print_kind=print_kind)
     398              :             END IF
     399              :             CALL open_file(file_name=data_file(ibox), &
     400              :                            unit_number=data_unit(ibox), file_position='APPEND', &
     401           12 :                            file_action='WRITE', file_status='UNKNOWN')
     402              :             CALL open_file(file_name=moves_file(ibox), &
     403              :                            unit_number=move_unit(ibox), file_position='APPEND', &
     404           12 :                            file_action='WRITE', file_status='UNKNOWN')
     405              :             CALL open_file(file_name=displacement_file(ibox), &
     406              :                            unit_number=rm(ibox), file_position='APPEND', &
     407           12 :                            file_action='WRITE', file_status='UNKNOWN')
     408              :             CALL open_file(file_name=cell_file(ibox), &
     409              :                            unit_number=cl(ibox), file_position='APPEND', &
     410           21 :                            file_action='WRITE', file_status='UNKNOWN')
     411              : 
     412              :          END DO
     413              : 
     414              :          ! back to parallel mode
     415              :       END IF
     416              : 
     417           42 :       DO ibox = 1, nboxes
     418           24 :          CALL group%bcast(cl(ibox), source)
     419           24 :          CALL group%bcast(rm(ibox), source)
     420           24 :          CALL group%bcast(diff(ibox), source)
     421              :          ! set all the units numbers that we just opened in the respective mc_par
     422              :          CALL set_mc_par(mc_par(ibox)%mc_par, cl=cl(ibox), rm=rm(ibox), &
     423           42 :                          diff=diff(ibox))
     424              :       END DO
     425              : 
     426              :       ! if we're doing a discrete volume move, we need to set up the array
     427              :       ! that keeps track of which direction we can move in
     428           18 :       IF (ldiscrete) THEN
     429            0 :          IF (nboxes .NE. 1) &
     430            0 :             CPABORT('ldiscrete=.true. ONLY for systems with 1 box')
     431              :          CALL create_discrete_array(abc(:, 1), discrete_array(:, :), &
     432            0 :                                     discrete_step)
     433              :       END IF
     434              : 
     435              :       ! find out how many steps we're doing...change the updates to be in cycles
     436              :       ! if the total number of steps is measured in cycles
     437           18 :       IF (.NOT. lstop) THEN
     438           10 :          nstep = nstep*nchain_total
     439           10 :          iuptrans = iuptrans*nchain_total
     440           10 :          iupvolume = iupvolume*nchain_total
     441              :       END IF
     442              : 
     443          486 :       DO nnstep = nstart + 1, nstart + nstep
     444              : 
     445          468 :          IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
     446           15 :             WRITE (iw, *)
     447           15 :             WRITE (iw, *) "------- On Monte Carlo Step ", nnstep
     448              :          END IF
     449              : 
     450          468 :          IF (ionode) rand = rng_stream%next()
     451              :          ! broadcast the random number, to make sure we're on the same move
     452          468 :          CALL group%bcast(rand, source)
     453              : 
     454          468 :          IF (rand .LT. pmvolume) THEN
     455              : 
     456           58 :             IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
     457            1 :                WRITE (iw, *) "Attempting a volume move"
     458            1 :                WRITE (iw, *)
     459              :             END IF
     460              : 
     461            8 :             SELECT CASE (ensemble)
     462              :             CASE ("TRADITIONAL")
     463              :                CALL mc_volume_move(mc_par(1)%mc_par, &
     464              :                                    force_env(1)%force_env, &
     465              :                                    moves(1, 1)%moves, move_updates(1, 1)%moves, &
     466              :                                    old_energy(1), 1, &
     467              :                                    energy_check(1), r_old(:, :, 1), iw, discrete_array(:, :), &
     468            8 :                                    rng_stream)
     469              :             CASE ("GEMC_NVT")
     470              :                CALL mc_ge_volume_move(mc_par, force_env, moves, &
     471              :                                       move_updates, nnstep, old_energy, energy_check, &
     472           24 :                                       r_old, rng_stream)
     473              :             CASE ("GEMC_NPT")
     474              :                ! we need to select a box based on the probability given in the input file
     475           26 :                IF (ionode) rand = rng_stream%next()
     476           26 :                CALL group%bcast(rand, source)
     477              : 
     478           38 :                DO ibox = 1, nboxes
     479           38 :                   IF (rand .LE. pmvol_box(ibox)) THEN
     480           26 :                      box_number = ibox
     481           26 :                      EXIT
     482              :                   END IF
     483              :                END DO
     484              : 
     485              :                CALL mc_volume_move(mc_par(box_number)%mc_par, &
     486              :                                    force_env(box_number)%force_env, &
     487              :                                    moves(1, box_number)%moves, &
     488              :                                    move_updates(1, box_number)%moves, &
     489              :                                    old_energy(box_number), box_number, &
     490              :                                    energy_check(box_number), r_old(:, :, box_number), iw, &
     491              :                                    discrete_array(:, :), &
     492           84 :                                    rng_stream)
     493              :             END SELECT
     494              : 
     495              : ! update all the pointers here, because otherwise we may pass wrong information when we're making a bias environment
     496          166 :             DO ibox = 1, nboxes
     497              :                CALL force_env_get(force_env(ibox)%force_env, &
     498          108 :                                   subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
     499          108 :                CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
     500              :                CALL cp_subsys_get(oldsys(ibox)%subsys, &
     501          166 :                                   particles=particles_old(ibox)%list)
     502              :             END DO
     503              : 
     504              :             ! we need a new biasing environment now, if we're into that sort of thing
     505           58 :             IF (lbias) THEN
     506          150 :                DO ibox = 1, nboxes
     507          100 :                   CALL force_env_release(bias_env(ibox)%force_env)
     508              :                   ! determine the atom names of every particle
     509          300 :                   ALLOCATE (atom_names_box(1:nunits_tot(ibox)))
     510          150 :                   start_mol = 1
     511          150 :                   DO jbox = 1, ibox - 1
     512          250 :                      start_mol = start_mol + SUM(nchains(:, jbox))
     513              :                   END DO
     514          300 :                   end_mol = start_mol + SUM(nchains(:, ibox)) - 1
     515          300 :                   atom_number = 1
     516         1500 :                   DO imolecule = 1, SUM(nchains(:, ibox))
     517         3800 :                      DO iunit = 1, nunits(mol_type(imolecule + start_mol - 1))
     518              :                         atom_names_box(atom_number) = &
     519         2500 :                            atom_names(iunit, mol_type(imolecule + start_mol - 1))
     520         3700 :                         atom_number = atom_number + 1
     521              :                      END DO
     522              :                   END DO
     523              : 
     524              : ! need to find out what the cell lengths are
     525              :                   CALL force_env_get(force_env(ibox)%force_env, &
     526          100 :                                      subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
     527          100 :                   CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
     528              : 
     529              :                   CALL get_mc_par(mc_par(ibox)%mc_par, &
     530          100 :                                   mc_bias_file=mc_bias_file)
     531          100 :                   nchains_box => nchains(:, ibox)
     532              : 
     533              :                   CALL mc_create_bias_force_env(bias_env(ibox)%force_env, &
     534              :                                                 r_old(:, :, ibox), atom_names_box(:), nunits_tot(ibox), &
     535              :                                                 para_env, abc(:, ibox), nchains_box, input_declaration, &
     536          100 :                                                 mc_bias_file, ionode)
     537              : 
     538          300 :                   IF (SUM(nchains(:, ibox)) .NE. 0) THEN
     539              :                      CALL force_env_calc_energy_force( &
     540              :                         bias_env(ibox)%force_env, &
     541          100 :                         calc_force=.FALSE.)
     542              :                      CALL force_env_get(bias_env(ibox)%force_env, &
     543          100 :                                         potential_energy=last_bias_energy(ibox))
     544              :                   ELSE
     545            0 :                      last_bias_energy(ibox) = 0.0E0_dp
     546              :                   END IF
     547          100 :                   bias_energy(ibox) = last_bias_energy(ibox)
     548          150 :                   DEALLOCATE (atom_names_box)
     549              :                END DO
     550              :             END IF
     551              : 
     552          410 :          ELSEIF (rand .LT. pmswap) THEN
     553              : 
     554              :             ! try a swap move
     555           22 :             IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
     556            0 :                WRITE (iw, *) "Attempting a swap move"
     557            0 :                WRITE (iw, *)
     558              :             END IF
     559              : 
     560              :             CALL mc_ge_swap_move(mc_par, force_env, bias_env, moves, &
     561              :                                  energy_check(:), r_old(:, :, :), old_energy(:), input_declaration, &
     562           22 :                                  para_env, bias_energy(:), last_bias_energy(:), rng_stream)
     563              : 
     564              :             ! the number of molecules may have changed, which deallocated the whole
     565              :             ! mc_molecule_info structure
     566           22 :             CALL get_mc_par(mc_par(1)%mc_par, mc_molecule_info=mc_molecule_info)
     567              :             CALL get_mc_molecule_info(mc_molecule_info, conf_prob=conf_prob, &
     568              :                                       nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
     569              :                                       mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
     570           22 :                                       atom_names=atom_names, mass=mass)
     571              : 
     572          388 :          ELSEIF (rand .LT. pmhmc) THEN
     573              : ! try hybrid Monte Carlo
     574           20 :             IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
     575            2 :                WRITE (iw, *) "Attempting a hybrid Monte Carlo move"
     576            2 :                WRITE (iw, *)
     577              :             END IF
     578              : 
     579              : ! pick a box at random
     580           20 :             IF (ionode) rand = rng_stream%next()
     581           20 :             CALL group%bcast(rand, source)
     582              : 
     583           20 :             DO ibox = 1, nboxes
     584           20 :                IF (rand .LE. pmhmc_box(ibox)) THEN
     585           20 :                   box_number = ibox
     586           20 :                   EXIT
     587              :                END IF
     588              :             END DO
     589              : 
     590              :             CALL mc_hmc_move(mc_par(box_number)%mc_par, &
     591              :                              force_env(box_number)%force_env, globenv, &
     592              :                              moves(1, box_number)%moves, &
     593              :                              move_updates(1, box_number)%moves, &
     594              :                              old_energy(box_number), box_number, &
     595              :                              energy_check(box_number), r_old(:, :, box_number), &
     596           20 :                              rng_stream)
     597              : 
     598          368 :          ELSEIF (rand .LT. pmavbmc) THEN
     599              :             ! try an AVBMC move
     600            0 :             IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
     601            0 :                WRITE (iw, *) "Attempting an AVBMC1 move"
     602            0 :                WRITE (iw, *)
     603              :             END IF
     604              : 
     605              :             ! first, pick a box to do it for
     606            0 :             IF (ionode) rand = rng_stream%next()
     607            0 :             CALL group%bcast(rand, source)
     608              : 
     609            0 :             IF (nboxes .EQ. 2) THEN
     610            0 :                IF (rand .LT. 0.1E0_dp) THEN
     611            0 :                   box_number = 1
     612              :                ELSE
     613            0 :                   box_number = 2
     614              :                END IF
     615              :             ELSE
     616            0 :                box_number = 1
     617              :             END IF
     618              : 
     619              :             ! now pick a molecule type to do it for
     620            0 :             IF (ionode) rand = rng_stream%next()
     621            0 :             CALL group%bcast(rand, source)
     622            0 :             molecule_type_swap = 0
     623            0 :             DO imol_type = 1, nmol_types
     624            0 :                IF (rand .LT. pmavbmc_mol(imol_type)) THEN
     625            0 :                   molecule_type_swap = imol_type
     626            0 :                   EXIT
     627              :                END IF
     628              :             END DO
     629            0 :             IF (molecule_type_swap == 0) &
     630            0 :                CPABORT('Did not choose a molecule type to swap...check AVBMC input')
     631              : 
     632              :             ! now pick a molecule, automatically rejecting the move if the
     633              :             ! box is empty or only has one molecule
     634            0 :             IF (SUM(nchains(:, box_number)) .LE. 1) THEN
     635              :                ! indicate that we tried a move
     636              :                moves(molecule_type_swap, box_number)%moves%empty_avbmc = &
     637            0 :                   moves(molecule_type_swap, box_number)%moves%empty_avbmc + 1
     638              :             ELSE
     639              : 
     640              :                ! pick a molecule to be swapped in the box
     641            0 :                IF (ionode) THEN
     642              :                   CALL find_mc_test_molecule(mc_molecule_info, &
     643              :                                              start_atom_swap, idum, jdum, rng_stream, &
     644            0 :                                              box=box_number, molecule_type_old=molecule_type_swap)
     645              : 
     646              :                   ! pick a molecule to act as the target in the box...we don't care what type
     647            0 :                   DO
     648              :                      CALL find_mc_test_molecule(mc_molecule_info, &
     649              :                                                 start_atom_target, idum, molecule_type_target, &
     650            0 :                                                 rng_stream, box=box_number)
     651            0 :                      IF (start_atom_swap .NE. start_atom_target) THEN
     652              :                         start_atom_target = start_atom_target + &
     653            0 :                                             avbmc_atom(molecule_type_target) - 1
     654              :                         EXIT
     655              :                      END IF
     656              :                   END DO
     657              : 
     658              :                   ! choose if we're swapping into the bonded region of mol_target, or
     659              :                   ! into the nonbonded region
     660            0 :                   rand = rng_stream%next()
     661              : 
     662              :                END IF
     663            0 :                CALL group%bcast(start_atom_swap, source)
     664            0 :                CALL group%bcast(box_number, source)
     665            0 :                CALL group%bcast(start_atom_target, source)
     666            0 :                CALL group%bcast(rand, source)
     667              : 
     668            0 :                IF (rand .LT. pbias(molecule_type_swap)) THEN
     669            0 :                   move_type_avbmc = 'in'
     670              :                ELSE
     671            0 :                   move_type_avbmc = 'out'
     672              :                END IF
     673              : 
     674              :                CALL mc_avbmc_move(mc_par(box_number)%mc_par, &
     675              :                                   force_env(box_number)%force_env, &
     676              :                                   bias_env(box_number)%force_env, &
     677              :                                   moves(molecule_type_swap, box_number)%moves, &
     678              :                                   energy_check(box_number), &
     679              :                                   r_old(:, :, box_number), old_energy(box_number), &
     680              :                                   start_atom_swap, start_atom_target, molecule_type_swap, &
     681              :                                   box_number, bias_energy(box_number), &
     682              :                                   last_bias_energy(box_number), &
     683            0 :                                   move_type_avbmc, rng_stream)
     684              : 
     685              :             END IF
     686              : 
     687              :          ELSE
     688              : 
     689          368 :             IF (MOD(nnstep, iprint) == 0 .AND. (iw > 0)) THEN
     690           12 :                WRITE (iw, *) "Attempting an inner move"
     691           12 :                WRITE (iw, *)
     692              :             END IF
     693              : 
     694         2010 :             DO imove = 1, nmoves
     695              : 
     696         1642 :                IF (ionode) rand = rng_stream%next()
     697         1642 :                CALL group%bcast(rand, source)
     698         2010 :                IF (rand .LT. pmtraion) THEN
     699              :                   ! change molecular conformation
     700              :                   ! first, pick a box to do it for
     701          506 :                   IF (ionode) rand = rng_stream%next()
     702          506 :                   CALL group%bcast(rand, source)
     703          506 :                   IF (nboxes .EQ. 2) THEN
     704            0 :                      IF (rand .LT. 0.75E0_dp) THEN
     705            0 :                         box_number = 1
     706              :                      ELSE
     707            0 :                         box_number = 2
     708              :                      END IF
     709              :                   ELSE
     710          506 :                      box_number = 1
     711              :                   END IF
     712              : 
     713              :                   ! figure out which molecule type we're looking for
     714          506 :                   IF (ionode) rand = rng_stream%next()
     715          506 :                   CALL group%bcast(rand, source)
     716          506 :                   molecule_type = 0
     717          506 :                   DO imol_type = 1, nmol_types
     718          506 :                      IF (rand .LT. pmtraion_mol(imol_type)) THEN
     719          506 :                         molecule_type = imol_type
     720          506 :                         EXIT
     721              :                      END IF
     722              :                   END DO
     723          506 :                   IF (molecule_type == 0) CALL cp_abort( &
     724              :                      __LOCATION__, &
     725            0 :                      'Did not choose a molecule type to conf change...PMTRAION_MOL should not be all 0.0')
     726              : 
     727              :                   ! now pick a molecule, automatically rejecting the move if the
     728              :                   ! box is empty
     729          506 :                   IF (nchains(molecule_type, box_number) == 0) THEN
     730              :                      ! indicate that we tried a move
     731              :                      moves(molecule_type, box_number)%moves%empty_conf = &
     732            0 :                         moves(molecule_type, box_number)%moves%empty_conf + 1
     733              :                   ELSE
     734              :                      ! pick a molecule in the box
     735          506 :                      IF (ionode) THEN
     736              :                         CALL find_mc_test_molecule(mc_molecule_info, &
     737              :                                                    start_atom, idum, &
     738              :                                                    jdum, rng_stream, &
     739          253 :                                                    box=box_number, molecule_type_old=molecule_type)
     740              : 
     741              :                         ! choose if we're changing a bond length or an angle
     742          253 :                         rand = rng_stream%next()
     743              :                      END IF
     744          506 :                      CALL group%bcast(rand, source)
     745          506 :                      CALL group%bcast(start_atom, source)
     746          506 :                      CALL group%bcast(box_number, source)
     747          506 :                      CALL group%bcast(molecule_type, source)
     748              : 
     749              :                      ! figure out what kind of move we're doing
     750          506 :                      IF (rand .LT. conf_prob(1, molecule_type)) THEN
     751          312 :                         move_type = 'bond'
     752          194 :                      ELSEIF (rand .LT. (conf_prob(1, molecule_type) + &
     753              :                                         conf_prob(2, molecule_type))) THEN
     754          194 :                         move_type = 'angle'
     755              :                      ELSE
     756            0 :                         move_type = 'dihedral'
     757              :                      END IF
     758          506 :                      box_flag(box_number) = 1
     759              :                      CALL mc_conformation_change(mc_par(box_number)%mc_par, &
     760              :                                                  force_env(box_number)%force_env, &
     761              :                                                  bias_env(box_number)%force_env, &
     762              :                                                  moves(molecule_type, box_number)%moves, &
     763              :                                                  move_updates(molecule_type, box_number)%moves, &
     764              :                                                  start_atom, molecule_type, box_number, &
     765              :                                                  bias_energy(box_number), &
     766          506 :                                                  move_type, lreject, rng_stream)
     767          506 :                      IF (lreject) EXIT
     768              :                   END IF
     769         1136 :                ELSEIF (rand .LT. pmtrans) THEN
     770              :                   ! translate a whole molecule in the system
     771              :                   ! pick a molecule type
     772          624 :                   IF (ionode) rand = rng_stream%next()
     773          624 :                   CALL group%bcast(rand, source)
     774          624 :                   molecule_type = 0
     775          924 :                   DO imol_type = 1, nmol_types
     776          924 :                      IF (rand .LT. pmtrans_mol(imol_type)) THEN
     777          624 :                         molecule_type = imol_type
     778          624 :                         EXIT
     779              :                      END IF
     780              :                   END DO
     781          624 :                   IF (molecule_type == 0) CALL cp_abort( &
     782              :                      __LOCATION__, &
     783            0 :                      'Did not choose a molecule type to translate...PMTRANS_MOL should not be all 0.0')
     784              : 
     785              :                   ! now pick a molecule of that type
     786          624 :                   IF (ionode) &
     787              :                      CALL find_mc_test_molecule(mc_molecule_info, &
     788              :                                                 start_atom, box_number, idum, rng_stream, &
     789          312 :                                                 molecule_type_old=molecule_type)
     790          624 :                   CALL group%bcast(start_atom, source)
     791          624 :                   CALL group%bcast(box_number, source)
     792          624 :                   box_flag(box_number) = 1
     793              :                   CALL mc_molecule_translation(mc_par(box_number)%mc_par, &
     794              :                                                force_env(box_number)%force_env, &
     795              :                                                bias_env(box_number)%force_env, &
     796              :                                                moves(molecule_type, box_number)%moves, &
     797              :                                                move_updates(molecule_type, box_number)%moves, &
     798              :                                                start_atom, box_number, bias_energy(box_number), &
     799          624 :                                                molecule_type, lreject, rng_stream)
     800          624 :                   IF (lreject) EXIT
     801          512 :                ELSEIF (rand .LT. pmcltrans) THEN
     802              :                   ! translate a whole cluster in the system
     803              :                   ! first, pick a box to do it for
     804           10 :                   IF (ionode) rand = rng_stream%next()
     805           10 :                   CALL group%bcast(rand, source)
     806              : 
     807           10 :                   DO ibox = 1, nboxes
     808           10 :                   IF (rand .LE. pmclus_box(ibox)) THEN
     809           10 :                      box_number = ibox
     810           10 :                      EXIT
     811              :                   END IF
     812              :                   END DO
     813           10 :                   box_flag(box_number) = 1
     814              :                   CALL mc_cluster_translation(mc_par(box_number)%mc_par, &
     815              :                                               force_env(box_number)%force_env, &
     816              :                                               bias_env(box_number)%force_env, &
     817              :                                               moves(1, box_number)%moves, &
     818              :                                               move_updates(1, box_number)%moves, &
     819              :                                               box_number, bias_energy(box_number), &
     820           10 :                                               lreject, rng_stream)
     821           10 :                   IF (lreject) EXIT
     822              :                ELSE
     823              :                   !     rotate a whole molecule in the system
     824              :                   ! pick a molecule type
     825          502 :                   IF (ionode) rand = rng_stream%next()
     826          502 :                   CALL group%bcast(rand, source)
     827          502 :                   molecule_type = 0
     828          502 :                   DO imol_type = 1, nmol_types
     829          502 :                      IF (rand .LT. pmrot_mol(imol_type)) THEN
     830          502 :                         molecule_type = imol_type
     831          502 :                         EXIT
     832              :                      END IF
     833              :                   END DO
     834          502 :                   IF (molecule_type == 0) CALL cp_abort( &
     835              :                      __LOCATION__, &
     836            0 :                      'Did not choose a molecule type to rotate...PMROT_MOL should not be all 0.0')
     837              : 
     838          502 :                   IF (ionode) &
     839              :                      CALL find_mc_test_molecule(mc_molecule_info, &
     840              :                                                 start_atom, box_number, idum, rng_stream, &
     841          251 :                                                 molecule_type_old=molecule_type)
     842          502 :                   CALL group%bcast(start_atom, source)
     843          502 :                   CALL group%bcast(box_number, source)
     844          502 :                   box_flag(box_number) = 1
     845              :                   CALL mc_molecule_rotation(mc_par(box_number)%mc_par, &
     846              :                                             force_env(box_number)%force_env, &
     847              :                                             bias_env(box_number)%force_env, &
     848              :                                             moves(molecule_type, box_number)%moves, &
     849              :                                             move_updates(molecule_type, box_number)%moves, &
     850              :                                             box_number, start_atom, &
     851              :                                             molecule_type, bias_energy(box_number), &
     852          502 :                                             lreject, rng_stream)
     853          502 :                   IF (lreject) EXIT
     854              :                END IF
     855              : 
     856              :             END DO
     857              : 
     858              :             ! now do a Quickstep calculation to see if we accept the sequence
     859              :             CALL mc_Quickstep_move(mc_par, force_env, bias_env, &
     860              :                                    moves, lreject, move_updates, energy_check(:), r_old(:, :, :), &
     861              :                                    nnstep, old_energy(:), bias_energy(:), last_bias_energy(:), &
     862              :                                    nboxes, box_flag(:), oldsys, particles_old, &
     863          368 :                                    rng_stream, unit_conv)
     864              : 
     865              :          END IF
     866              : 
     867              :          ! make sure the pointers are pointing correctly since the subsys may
     868              :          ! have changed
     869         1080 :          DO ibox = 1, nboxes
     870              :             CALL force_env_get(force_env(ibox)%force_env, &
     871          612 :                                subsys=oldsys(ibox)%subsys, cell=cell(ibox)%cell)
     872          612 :             CALL get_cell(cell(ibox)%cell, abc=abc(:, ibox))
     873              :             CALL cp_subsys_get(oldsys(ibox)%subsys, &
     874         1080 :                                particles=particles_old(ibox)%list)
     875              :          END DO
     876              : 
     877          468 :          IF (ionode) THEN
     878              : 
     879          234 :             IF (MOD(nnstep, iprint) == 0) THEN
     880           15 :                WRITE (com_ene, *) nnstep, old_energy(1:nboxes)
     881              : 
     882           33 :                DO ibox = 1, nboxes
     883              : 
     884              :                   ! write the molecule information
     885           18 :                   WRITE (com_mol, *) nnstep, nchains(:, ibox)
     886              : 
     887              :                   ! write the move statistics to file
     888           47 :                   DO itype = 1, nmol_types
     889              :                      CALL write_move_stats(moves(itype, ibox)%moves, &
     890           47 :                                            nnstep, move_unit(ibox))
     891              :                   END DO
     892              : 
     893              :                   ! write a restart file
     894              :                   CALL write_mc_restart(nnstep, mc_par(ibox)%mc_par, &
     895           18 :                                         nchains(:, ibox), force_env(ibox)%force_env)
     896              : 
     897              :                   ! write cell lengths
     898           72 :                   WRITE (cell_unit, *) nnstep, abc(1:3, ibox)*angstrom
     899              : 
     900              :                   ! write particle coordinates
     901           18 :                   WRITE (cbox, '(I4)') ibox
     902           18 :                   WRITE (cstep, '(I8)') nnstep
     903           62 :                   IF (SUM(nchains(:, ibox)) == 0) THEN
     904            0 :                      WRITE (com_crd, *) ' 0'
     905              :                      WRITE (com_crd, *) 'BOX '//TRIM(ADJUSTL(cbox))// &
     906            0 :                         ',  STEP '//TRIM(ADJUSTL(cstep))
     907              :                   ELSE
     908              :                      CALL write_particle_coordinates( &
     909              :                         particles_old(ibox)%list%els, &
     910              :                         com_crd, dump_xmol, 'POS', &
     911              :                         'BOX '//TRIM(ADJUSTL(cbox))// &
     912              :                         ',  STEP '//TRIM(ADJUSTL(cstep)), &
     913           18 :                         unit_conv=unit_conv)
     914              :                   END IF
     915              :                END DO
     916              :             END IF ! end the things we only do every iprint moves
     917              : 
     918          540 :             DO ibox = 1, nboxes
     919              :                ! compute some averages
     920              :                averages(ibox)%averages%ave_energy = &
     921              :                   averages(ibox)%averages%ave_energy*REAL(nnstep - &
     922              :                                                           nstart - 1, dp)/REAL(nnstep - nstart, dp) + &
     923          306 :                   old_energy(ibox)/REAL(nnstep - nstart, dp)
     924              :                averages(ibox)%averages%molecules = &
     925              :                   averages(ibox)%averages%molecules*REAL(nnstep - &
     926              :                                                          nstart - 1, dp)/REAL(nnstep - nstart, dp) + &
     927          899 :                   REAL(SUM(nchains(:, ibox)), dp)/REAL(nnstep - nstart, dp)
     928              :                averages(ibox)%averages%ave_volume = &
     929              :                   averages(ibox)%averages%ave_volume* &
     930              :                   REAL(nnstep - nstart - 1, dp)/REAL(nnstep - nstart, dp) + &
     931              :                   abc(1, ibox)*abc(2, ibox)*abc(3, ibox)/ &
     932          306 :                   REAL(nnstep - nstart, dp)
     933              : 
     934              :                ! flush the buffers to the files
     935          306 :                CALL m_flush(data_unit(ibox))
     936          306 :                CALL m_flush(diff(ibox))
     937          306 :                CALL m_flush(move_unit(ibox))
     938          306 :                CALL m_flush(cl(ibox))
     939          540 :                CALL m_flush(rm(ibox))
     940              : 
     941              :             END DO
     942              : 
     943              :             ! flush more buffers to the files
     944          234 :             CALL m_flush(cell_unit)
     945          234 :             CALL m_flush(com_ene)
     946          234 :             CALL m_flush(com_crd)
     947          234 :             CALL m_flush(com_mol)
     948              : 
     949              :          END IF
     950              : 
     951              :          ! reset the box flags
     952         1080 :          box_flag(:) = 0
     953              : 
     954              :          ! check to see if EXIT file exists...if so, end the calculation
     955          468 :          CALL external_control(should_stop, "MC", globenv=globenv)
     956          468 :          IF (should_stop) EXIT
     957              : 
     958              :          ! update the move displacements, if necessary
     959         1080 :          DO ibox = 1, nboxes
     960          612 :             IF (MOD(nnstep - nstart, iuptrans) == 0) THEN
     961            0 :                DO itype = 1, nmol_types
     962              :                   CALL mc_move_update(mc_par(ibox)%mc_par, &
     963              :                                       move_updates(itype, ibox)%moves, itype, &
     964            0 :                                       "trans", nnstep, ionode)
     965              :                END DO
     966              :             END IF
     967              : 
     968         1080 :             IF (MOD(nnstep - nstart, iupvolume) == 0) THEN
     969              :                CALL mc_move_update(mc_par(ibox)%mc_par, &
     970              :                                    move_updates(1, ibox)%moves, 1337, &
     971            0 :                                    "volume", nnstep, ionode)
     972              :             END IF
     973              :          END DO
     974              : 
     975              :          ! check to see if there are any overlaps in the boxes, and fold coordinates
     976              : ! don't care about overlaps if we're only doing HMC
     977          468 :          IF (.NOT. lhmc) THEN
     978         1040 :             DO ibox = 1, nboxes
     979         2206 :                IF (SUM(nchains(:, ibox)) .NE. 0) THEN
     980              :                   start_mol = 1
     981          736 :                   DO jbox = 1, ibox - 1
     982         1024 :                      start_mol = start_mol + SUM(nchains(:, jbox))
     983              :                   END DO
     984         1758 :                   end_mol = start_mol + SUM(nchains(:, ibox)) - 1
     985              :                   CALL check_for_overlap(force_env(ibox)%force_env, &
     986              :                                          nchains(:, ibox), nunits, loverlap, &
     987          592 :                                          mol_type(start_mol:end_mol))
     988          592 :                   IF (loverlap) THEN
     989            0 :                      IF (iw > 0) WRITE (iw, *) nnstep
     990            0 :                      CPABORT('coordinate overlap at the end of the above step')
     991              :                      ! now fold the coordinates...don't do this anywhere but here, because
     992              :                      ! we can get screwed up with the mc_molecule_info stuff (like in swap move)...
     993              :                      ! this is kind of ugly, with allocated and deallocating every time
     994            0 :                      ALLOCATE (r_temp(1:3, 1:nunits_tot(ibox)))
     995              : 
     996            0 :                      DO iunit = 1, nunits_tot(ibox)
     997              :                         r_temp(1:3, iunit) = &
     998            0 :                            particles_old(ibox)%list%els(iunit)%r(1:3)
     999              :                      END DO
    1000              : 
    1001              :                      CALL mc_coordinate_fold(r_temp(:, :), &
    1002              :                                              SUM(nchains(:, ibox)), mol_type(start_mol:end_mol), &
    1003            0 :                                              mass, nunits, abc(1:3, ibox))
    1004              : 
    1005              :                      ! save the folded coordinates
    1006            0 :                      DO iunit = 1, nunits_tot(ibox)
    1007            0 :                         r_old(1:3, iunit, ibox) = r_temp(1:3, iunit)
    1008              :                         particles_old(ibox)%list%els(iunit)%r(1:3) = &
    1009            0 :                            r_temp(1:3, iunit)
    1010              :                      END DO
    1011              : 
    1012              :                      ! if we're biasing, we need to do the same
    1013            0 :                      IF (lbias) THEN
    1014              :                         CALL force_env_get(bias_env(ibox)%force_env, &
    1015            0 :                                            subsys=biassys)
    1016              :                         CALL cp_subsys_get(biassys, &
    1017            0 :                                            particles=particles_bias)
    1018              : 
    1019            0 :                         DO iunit = 1, nunits_tot(ibox)
    1020              :                            particles_bias%els(iunit)%r(1:3) = &
    1021            0 :                               r_temp(1:3, iunit)
    1022              :                         END DO
    1023              :                      END IF
    1024              : 
    1025            0 :                      DEALLOCATE (r_temp)
    1026              :                   END IF
    1027              :                END IF
    1028              :             END DO
    1029              :          END IF
    1030              : 
    1031              :          !debug code
    1032          486 :          IF (debug_this_module) THEN
    1033              :             DO ibox = 1, nboxes
    1034              :                IF (SUM(nchains(:, ibox)) .NE. 0) THEN
    1035              :                   CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
    1036              :                                                    calc_force=.FALSE.)
    1037              :                   CALL force_env_get(force_env(ibox)%force_env, &
    1038              :                                      potential_energy=test_energy)
    1039              :                ELSE
    1040              :                   test_energy = 0.0E0_dp
    1041              :                END IF
    1042              : 
    1043              :                IF (ABS(initial_energy(ibox) + energy_check(ibox) - &
    1044              :                        test_energy) .GT. 0.0000001E0_dp) THEN
    1045              :                   IF (iw > 0) THEN
    1046              :                      WRITE (iw, *) '!!!!!!! We have an energy problem. !!!!!!!!'
    1047              :                      WRITE (iw, '(A,T64,F16.10)') 'Final Energy = ', test_energy
    1048              :                      WRITE (iw, '(A,T64,F16.10)') 'Initial Energy+energy_check=', &
    1049              :                         initial_energy(ibox) + energy_check(ibox)
    1050              :                      WRITE (iw, *) 'Box ', ibox
    1051              :                      WRITE (iw, *) 'nchains ', nchains(:, ibox)
    1052              :                   END IF
    1053              :                   CPABORT('!!!!!!! We have an energy problem. !!!!!!!!')
    1054              :                END IF
    1055              :             END DO
    1056              :          END IF
    1057              :       END DO
    1058              : 
    1059              :       ! write a restart file
    1060           18 :       IF (ionode) THEN
    1061           21 :          DO ibox = 1, nboxes
    1062              :             CALL write_mc_restart(nnstep, mc_par(ibox)%mc_par, &
    1063           21 :                                   nchains(:, ibox), force_env(ibox)%force_env)
    1064              :          END DO
    1065              :       END IF
    1066              : 
    1067              :       ! calculate the final energy
    1068           42 :       DO ibox = 1, nboxes
    1069           64 :          IF (SUM(nchains(:, ibox)) .NE. 0) THEN
    1070              :             CALL force_env_calc_energy_force(force_env(ibox)%force_env, &
    1071           24 :                                              calc_force=.FALSE.)
    1072              :             CALL force_env_get(force_env(ibox)%force_env, &
    1073           24 :                                potential_energy=final_energy(ibox))
    1074              :          ELSE
    1075            0 :             final_energy(ibox) = 0.0E0_dp
    1076              :          END IF
    1077           42 :          IF (lbias) THEN
    1078           14 :             CALL force_env_release(bias_env(ibox)%force_env)
    1079              :          END IF
    1080              :       END DO
    1081              : 
    1082              :       ! do some stuff in serial
    1083           18 :       IF (ionode .OR. (iw > 0)) THEN
    1084              : 
    1085            9 :          WRITE (com_ene, *) 'Final Energies:                      ', &
    1086           18 :             final_energy(1:nboxes)
    1087              : 
    1088           21 :          DO ibox = 1, nboxes
    1089           12 :             WRITE (cbox, '(I4)') ibox
    1090           32 :             IF (SUM(nchains(:, ibox)) == 0) THEN
    1091            0 :                WRITE (com_crd, *) ' 0'
    1092            0 :                WRITE (com_crd, *) 'BOX '//TRIM(ADJUSTL(cbox))
    1093              :             ELSE
    1094              :                CALL write_particle_coordinates( &
    1095              :                   particles_old(ibox)%list%els, &
    1096              :                   com_crd, dump_xmol, 'POS', &
    1097           12 :                   'FINAL BOX '//TRIM(ADJUSTL(cbox)), unit_conv=unit_conv)
    1098              :             END IF
    1099              : 
    1100              :             ! write a bunch of data to the screen
    1101              :             WRITE (iw, '(A)') &
    1102           12 :                '------------------------------------------------'
    1103              :             WRITE (iw, '(A,I1,A)') &
    1104           12 :                '|                   BOX ', ibox, &
    1105           24 :                '                      |'
    1106              :             WRITE (iw, '(A)') &
    1107           12 :                '------------------------------------------------'
    1108           12 :             test_moves => moves(:, ibox)
    1109              :             CALL final_mc_write(mc_par(ibox)%mc_par, test_moves, &
    1110              :                                 iw, energy_check(ibox), &
    1111              :                                 initial_energy(ibox), final_energy(ibox), &
    1112           12 :                                 averages(ibox)%averages)
    1113              : 
    1114              :             ! close any open files
    1115           12 :             CALL close_file(unit_number=diff(ibox))
    1116           12 :             CALL close_file(unit_number=data_unit(ibox))
    1117           12 :             CALL close_file(unit_number=move_unit(ibox))
    1118           12 :             CALL close_file(unit_number=cl(ibox))
    1119           21 :             CALL close_file(unit_number=rm(ibox))
    1120              :          END DO
    1121              : 
    1122              :          ! close some more files
    1123            9 :          CALL close_file(unit_number=cell_unit)
    1124            9 :          CALL close_file(unit_number=com_ene)
    1125            9 :          CALL close_file(unit_number=com_crd)
    1126            9 :          CALL close_file(unit_number=com_mol)
    1127              :       END IF
    1128              : 
    1129           42 :       DO ibox = 1, nboxes
    1130              :          CALL set_mc_env(mc_env(ibox)%mc_env, &
    1131              :                          mc_par=mc_par(ibox)%mc_par, &
    1132           24 :                          force_env=force_env(ibox)%force_env)
    1133              : 
    1134              :          ! deallocate some stuff
    1135           64 :          DO itype = 1, nmol_types
    1136           40 :             CALL mc_moves_release(move_updates(itype, ibox)%moves)
    1137           64 :             CALL mc_moves_release(moves(itype, ibox)%moves)
    1138              :          END DO
    1139           42 :          CALL mc_averages_release(averages(ibox)%averages)
    1140              :       END DO
    1141              : 
    1142           18 :       DEALLOCATE (pmhmc_box)
    1143           18 :       DEALLOCATE (pmvol_box)
    1144           18 :       DEALLOCATE (pmclus_box)
    1145           18 :       DEALLOCATE (r_old)
    1146           18 :       DEALLOCATE (force_env)
    1147           18 :       DEALLOCATE (bias_env)
    1148           18 :       DEALLOCATE (cell)
    1149           18 :       DEALLOCATE (particles_old)
    1150           18 :       DEALLOCATE (oldsys)
    1151           18 :       DEALLOCATE (averages)
    1152           18 :       DEALLOCATE (moves)
    1153           18 :       DEALLOCATE (move_updates)
    1154           18 :       DEALLOCATE (mc_par)
    1155              : 
    1156              :       ! end the timing
    1157           18 :       CALL timestop(handle)
    1158              : 
    1159           54 :    END SUBROUTINE mc_run_ensemble
    1160              : 
    1161              : ! **************************************************************************************************
    1162              : !> \brief Computes the second virial coefficient of a molecule by using the integral form
    1163              : !>      of the second virial coefficient found in McQuarrie "Statistical Thermodynamics",
    1164              : !>      B2(T) = -2Pi Int 0toInf [ Exp[-beta*u(r)] -1] r^2 dr     Eq. 15-25
    1165              : !>      I use trapazoidal integration with various step sizes
    1166              : !>      (the integral is broken up into three parts, currently, but that's easily
    1167              : !>      changed by the first variables found below).  It generates nvirial configurations,
    1168              : !>      doing the integration for each one, and then averages all the B2(T) to produce
    1169              : !>      the final answer.
    1170              : !> \param mc_env a pointer that contains all mc_env for all the simulation
    1171              : !>          boxes
    1172              : !> \param rng_stream the stream we pull random numbers from
    1173              : !>
    1174              : !>    Suitable for parallel.
    1175              : !> \author MJM
    1176              : ! **************************************************************************************************
    1177            2 :    SUBROUTINE mc_compute_virial(mc_env, rng_stream)
    1178              : 
    1179              :       TYPE(mc_environment_p_type), DIMENSION(:), POINTER :: mc_env
    1180              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
    1181              : 
    1182              :       INTEGER :: current_division, end_atom, ibin, idivision, iparticle, iprint, itemp, iunit, &
    1183              :          ivirial, iw, nbins, nchain_total, nintegral_divisions, nmol_types, nvirial, &
    1184              :          nvirial_temps, source, start_atom
    1185            2 :       INTEGER, DIMENSION(:), POINTER                     :: mol_type, nunits, nunits_tot
    1186            2 :       INTEGER, DIMENSION(:, :), POINTER                  :: nchains
    1187              :       LOGICAL                                            :: ionode, loverlap
    1188            2 :       REAL(dp), DIMENSION(:), POINTER                    :: BETA, virial_cutoffs, virial_stepsize, &
    1189            2 :                                                             virial_temps
    1190            2 :       REAL(dp), DIMENSION(:, :), POINTER                 :: mass, mayer, r_old
    1191              :       REAL(KIND=dp) :: ave_virial, current_value, distance, exp_max_val, exp_min_val, exponent, &
    1192              :          integral, previous_value, square_value, trial_energy, triangle_value
    1193              :       REAL(KIND=dp), DIMENSION(1:3)                      :: abc, center_of_mass
    1194            2 :       TYPE(cell_p_type), DIMENSION(:), POINTER           :: cell
    1195            2 :       TYPE(cp_subsys_p_type), DIMENSION(:), POINTER      :: subsys
    1196            2 :       TYPE(force_env_p_type), DIMENSION(:), POINTER      :: force_env
    1197              :       TYPE(mc_molecule_info_type), POINTER               :: mc_molecule_info
    1198              :       TYPE(mc_simulation_parameters_p_type), &
    1199            2 :          DIMENSION(:), POINTER                           :: mc_par
    1200              :       TYPE(mp_comm_type)                                 :: group
    1201            2 :       TYPE(particle_list_p_type), DIMENSION(:), POINTER  :: particles
    1202              : 
    1203              : ! these are current magic numbers for how we compute the virial...
    1204              : ! we break it up into three parts to integrate the function so provide
    1205              : ! better statistics
    1206              : 
    1207            2 :       nintegral_divisions = 3
    1208            0 :       ALLOCATE (virial_cutoffs(1:nintegral_divisions))
    1209            2 :       ALLOCATE (virial_stepsize(1:nintegral_divisions))
    1210            2 :       virial_cutoffs(1) = 8.0 ! first distance, in bohr
    1211            2 :       virial_cutoffs(2) = 13.0 ! second distance, in bohr
    1212            2 :       virial_cutoffs(3) = 22.0 ! maximum distance, in bohr
    1213            2 :       virial_stepsize(1) = 0.04 ! stepsize from 0 to virial_cutoffs(1)
    1214            2 :       virial_stepsize(2) = 0.1
    1215            2 :       virial_stepsize(3) = 0.2
    1216              : 
    1217              :       nbins = CEILING(virial_cutoffs(1)/virial_stepsize(1) + (virial_cutoffs(2) - virial_cutoffs(1))/ &
    1218            2 :                       virial_stepsize(2) + (virial_cutoffs(3) - virial_cutoffs(2))/virial_stepsize(3))
    1219              : 
    1220              :       ! figure out what the default write unit is
    1221            2 :       iw = cp_logger_get_default_io_unit()
    1222              : 
    1223              :       ! allocate a whole bunch of stuff based on how many boxes we have
    1224            4 :       ALLOCATE (force_env(1:1))
    1225            4 :       ALLOCATE (cell(1:1))
    1226            4 :       ALLOCATE (particles(1:1))
    1227            4 :       ALLOCATE (subsys(1:1))
    1228            4 :       ALLOCATE (mc_par(1:1))
    1229              : 
    1230              :       CALL get_mc_env(mc_env(1)%mc_env, &
    1231              :                       mc_par=mc_par(1)%mc_par, &
    1232            2 :                       force_env=force_env(1)%force_env)
    1233              : 
    1234              :       ! get some data out of mc_par
    1235              :       CALL get_mc_par(mc_par(1)%mc_par, &
    1236              :                       exp_max_val=exp_max_val, &
    1237              :                       exp_min_val=exp_min_val, nvirial=nvirial, &
    1238              :                       ionode=ionode, source=source, group=group, &
    1239            2 :                       mc_molecule_info=mc_molecule_info, virial_temps=virial_temps)
    1240              : 
    1241            2 :       IF (iw > 0) THEN
    1242            1 :          WRITE (iw, *)
    1243            1 :          WRITE (iw, *)
    1244            1 :          WRITE (iw, *) 'Beginning the calculation of the second virial coefficient'
    1245            1 :          WRITE (iw, *)
    1246            1 :          WRITE (iw, *)
    1247              :       END IF
    1248              : 
    1249              :       ! get some data from the molecule types
    1250              :       CALL get_mc_molecule_info(mc_molecule_info, &
    1251              :                                 nchains=nchains, nmol_types=nmol_types, nunits_tot=nunits_tot, &
    1252              :                                 mol_type=mol_type, nchain_total=nchain_total, nunits=nunits, &
    1253            2 :                                 mass=mass)
    1254              : 
    1255            2 :       nvirial_temps = SIZE(virial_temps)
    1256            6 :       ALLOCATE (BETA(1:nvirial_temps))
    1257              : 
    1258            6 :       DO itemp = 1, nvirial_temps
    1259            6 :          BETA(itemp) = 1/virial_temps(itemp)/boltzmann*joule
    1260              :       END DO
    1261              : 
    1262              :       ! get the subsystems and the cell information
    1263              :       CALL force_env_get(force_env(1)%force_env, &
    1264            2 :                          subsys=subsys(1)%subsys, cell=cell(1)%cell)
    1265            2 :       CALL get_cell(cell(1)%cell, abc=abc(:))
    1266              :       CALL cp_subsys_get(subsys(1)%subsys, &
    1267            2 :                          particles=particles(1)%list)
    1268              : 
    1269              :       ! check and make sure the box is big enough
    1270            2 :       IF (abc(1) .NE. abc(2) .OR. abc(2) .NE. abc(3)) THEN
    1271            0 :          CPABORT('The box needs to be cubic for a virial calculation (it is easiest).')
    1272              :       END IF
    1273            2 :       IF (virial_cutoffs(nintegral_divisions) .GT. abc(1)/2.0E0_dp) THEN
    1274            0 :          IF (iw > 0) &
    1275            0 :             WRITE (iw, *) "Box length ", abc(1)*angstrom, " virial cutoff ", &
    1276            0 :             virial_cutoffs(nintegral_divisions)*angstrom
    1277            0 :          CPABORT('You need a bigger box to deal with this virial cutoff (see above).')
    1278              :       END IF
    1279              : 
    1280              :       ! store the coordinates of the molecules in an array so we can work with it
    1281            6 :       ALLOCATE (r_old(1:3, 1:nunits_tot(1)))
    1282              : 
    1283           14 :       DO iparticle = 1, nunits_tot(1)
    1284              :          r_old(1:3, iparticle) = &
    1285           50 :             particles(1)%list%els(iparticle)%r(1:3)
    1286              :       END DO
    1287              : 
    1288              :       ! move the center of mass of molecule 1 to the origin
    1289            2 :       start_atom = 1
    1290            2 :       end_atom = nunits(mol_type(1))
    1291              :       CALL get_center_of_mass(r_old(:, start_atom:end_atom), nunits(mol_type(1)), &
    1292            2 :                               center_of_mass(:), mass(1:nunits(mol_type(1)), mol_type(1)))
    1293            8 :       DO iunit = start_atom, end_atom
    1294           26 :          r_old(:, iunit) = r_old(:, iunit) - center_of_mass(:)
    1295              :       END DO
    1296              :       ! set them in the force_env, so the first molecule is ready for the energy calc
    1297            8 :       DO iparticle = start_atom, end_atom
    1298           44 :          particles(1)%list%els(iparticle)%r(1:3) = r_old(1:3, iparticle)
    1299              :       END DO
    1300              : 
    1301              :       ! print out a notice every 1%
    1302            2 :       iprint = FLOOR(REAL(nvirial, KIND=dp)/100.0_dp)
    1303            2 :       IF (iprint == 0) iprint = 1
    1304              : 
    1305              :       ! we'll compute the average potential, and then integrate that, as opposed to
    1306              :       ! integrating every orientation and then averaging
    1307            8 :       ALLOCATE (mayer(1:nvirial_temps, 1:nbins))
    1308              : 
    1309         1778 :       mayer(:, :) = 0.0_dp
    1310              : 
    1311              :       ! loop over all nvirial random configurations
    1312           22 :       DO ivirial = 1, nvirial
    1313              : 
    1314              :          ! move molecule two back to the origin
    1315           20 :          start_atom = nunits(mol_type(1)) + 1
    1316           20 :          end_atom = nunits_tot(1)
    1317              :          CALL get_center_of_mass(r_old(:, start_atom:end_atom), nunits(mol_type(2)), &
    1318           20 :                                  center_of_mass(:), mass(1:nunits(mol_type(2)), mol_type(2)))
    1319           80 :          DO iunit = start_atom, end_atom
    1320          260 :             r_old(:, iunit) = r_old(:, iunit) - center_of_mass(:)
    1321              :          END DO
    1322              : 
    1323              :          ! now we need a random orientation for molecule 2...this routine is
    1324              :          ! only done in serial since it calls a random number
    1325           20 :          IF (ionode) THEN
    1326              :             CALL rotate_molecule(r_old(:, start_atom:end_atom), &
    1327              :                                  mass(1:nunits(mol_type(2)), mol_type(2)), &
    1328           10 :                                  nunits(mol_type(2)), rng_stream)
    1329              :          END IF
    1330          980 :          CALL group%bcast(r_old(:, :), source)
    1331              : 
    1332           20 :          distance = 0.0E0_dp
    1333           20 :          ibin = 1
    1334         5900 :          DO
    1335              :             ! find out what our stepsize is
    1336         5920 :             current_division = 0
    1337         8780 :             DO idivision = 1, nintegral_divisions
    1338         8780 :                IF (distance .LT. virial_cutoffs(idivision) - virial_stepsize(idivision)/2.0E0_dp) THEN
    1339              :                   current_division = idivision
    1340              :                   EXIT
    1341              :                END IF
    1342              :             END DO
    1343         5920 :             IF (current_division == 0) EXIT
    1344         5900 :             distance = distance + virial_stepsize(current_division)
    1345              : 
    1346              :             ! move the second molecule only along the x direction
    1347        23600 :             DO iparticle = start_atom, end_atom
    1348        17700 :                particles(1)%list%els(iparticle)%r(1) = r_old(1, iparticle) + distance
    1349        17700 :                particles(1)%list%els(iparticle)%r(2) = r_old(2, iparticle)
    1350        23600 :                particles(1)%list%els(iparticle)%r(3) = r_old(3, iparticle)
    1351              :             END DO
    1352              : 
    1353              :             ! check for overlaps
    1354         5900 :             CALL check_for_overlap(force_env(1)%force_env, nchains(:, 1), nunits, loverlap, mol_type)
    1355              : 
    1356              :             ! compute the energy if there is no overlap
    1357              :             ! exponent is exp(-beta*energy)-1, also called the Mayer term
    1358         5900 :             IF (loverlap) THEN
    1359         3834 :                DO itemp = 1, nvirial_temps
    1360         3834 :                   mayer(itemp, ibin) = mayer(itemp, ibin) - 1.0_dp
    1361              :                END DO
    1362              :             ELSE
    1363              :                CALL force_env_calc_energy_force(force_env(1)%force_env, &
    1364         4622 :                                                 calc_force=.FALSE.)
    1365              :                CALL force_env_get(force_env(1)%force_env, &
    1366         4622 :                                   potential_energy=trial_energy)
    1367              : 
    1368        13866 :                DO itemp = 1, nvirial_temps
    1369              : 
    1370         9244 :                   exponent = -BETA(itemp)*trial_energy
    1371              : 
    1372         9244 :                   IF (exponent .GT. exp_max_val) THEN
    1373              :                      exponent = exp_max_val
    1374         9244 :                   ELSEIF (exponent .LT. exp_min_val) THEN
    1375         1750 :                      exponent = exp_min_val
    1376              :                   END IF
    1377        13866 :                   mayer(itemp, ibin) = mayer(itemp, ibin) + EXP(exponent) - 1.0_dp
    1378              :                END DO
    1379              :             END IF
    1380              : 
    1381         5900 :             ibin = ibin + 1
    1382              :          END DO
    1383              :          ! write out some info that keeps track of where we are
    1384           22 :          IF (iw > 0) THEN
    1385           10 :             IF (MOD(ivirial, iprint) == 0) &
    1386           10 :                WRITE (iw, '(A,I6,A,I6)') ' Done with config ', ivirial, ' out of ', nvirial
    1387              :          END IF
    1388              :       END DO
    1389              : 
    1390              :       ! now we integrate this average potential
    1391         1778 :       mayer(:, :) = mayer(:, :)/REAL(nvirial, dp)
    1392              : 
    1393            6 :       DO itemp = 1, nvirial_temps
    1394              :          integral = 0.0_dp
    1395              :          previous_value = 0.0_dp
    1396              :          distance = 0.0E0_dp
    1397              :          ibin = 1
    1398         1180 :          DO
    1399         1184 :             current_division = 0
    1400         1756 :             DO idivision = 1, nintegral_divisions
    1401         1756 :                IF (distance .LT. virial_cutoffs(idivision) - virial_stepsize(idivision)/2.0E0_dp) THEN
    1402              :                   current_division = idivision
    1403              :                   EXIT
    1404              :                END IF
    1405              :             END DO
    1406         1184 :             IF (current_division == 0) EXIT
    1407         1180 :             distance = distance + virial_stepsize(current_division)
    1408              : 
    1409              :             ! now we need to integrate, using the trapazoidal method
    1410              :             ! first, find the value of the square
    1411         1180 :             current_value = mayer(itemp, ibin)*distance**2
    1412         1180 :             square_value = previous_value*virial_stepsize(current_division)
    1413              :             ! now the triangle that sits on top of it, which is half the size of this square...
    1414              :             ! notice this is negative if the current value is less than the previous value
    1415         1180 :             triangle_value = 0.5E0_dp*((current_value - previous_value)*virial_stepsize(current_division))
    1416              : 
    1417         1180 :             integral = integral + square_value + triangle_value
    1418         1180 :             previous_value = current_value
    1419         1180 :             ibin = ibin + 1
    1420              :          END DO
    1421              : 
    1422              :          ! now that the integration is done, compute the second virial that results
    1423            4 :          ave_virial = -2.0E0_dp*pi*integral
    1424              : 
    1425              :          ! convert from CP2K units to something else
    1426            4 :          ave_virial = ave_virial*n_avogadro*angstrom**3/1.0E8_dp**3
    1427              : 
    1428            6 :          IF (iw > 0) THEN
    1429            2 :             WRITE (iw, *)
    1430            2 :             WRITE (iw, *) '*********************************************************************'
    1431            2 :             WRITE (iw, '(A,F12.6,A)') ' ***                Temperature = ', virial_temps(itemp), &
    1432            4 :                '                     ***'
    1433            2 :             WRITE (iw, *) '***                                                               ***'
    1434            2 :             WRITE (iw, '(A,E12.6,A)') ' ***                  B2(T) = ', ave_virial, &
    1435            4 :                ' cm**3/mol               ***'
    1436            2 :             WRITE (iw, *) '*********************************************************************'
    1437            2 :             WRITE (iw, *)
    1438              :          END IF
    1439              :       END DO
    1440              : 
    1441              :       ! deallocate some stuff
    1442            2 :       DEALLOCATE (mc_par)
    1443            2 :       DEALLOCATE (subsys)
    1444            2 :       DEALLOCATE (force_env)
    1445            2 :       DEALLOCATE (particles)
    1446            2 :       DEALLOCATE (cell)
    1447            2 :       DEALLOCATE (virial_cutoffs)
    1448            2 :       DEALLOCATE (virial_stepsize)
    1449            2 :       DEALLOCATE (r_old)
    1450            2 :       DEALLOCATE (mayer)
    1451            2 :       DEALLOCATE (BETA)
    1452              : 
    1453            6 :    END SUBROUTINE mc_compute_virial
    1454              : 
    1455              : END MODULE mc_ensembles
    1456              : 
        

Generated by: LCOV version 2.0-1