LCOV - code coverage report
Current view: top level - src/motion/thermostat - thermostat_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 619 793 78.1 %
Date: 2024-03-29 07:50:05 Functions: 20 23 87.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Utilities for thermostats
      10             : !> \author teo [tlaino] - University of Zurich - 10.2007
      11             : ! **************************************************************************************************
      12             : MODULE thermostat_utils
      13             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      14             :                                               get_atomic_kind
      15             :    USE cell_types,                      ONLY: cell_type
      16             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      17             :                                               cp_logger_get_default_io_unit,&
      18             :                                               cp_logger_type,&
      19             :                                               cp_to_string
      20             :    USE cp_output_handling,              ONLY: cp_p_file,&
      21             :                                               cp_print_key_finished_output,&
      22             :                                               cp_print_key_should_output,&
      23             :                                               cp_print_key_unit_nr
      24             :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      25             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      26             :    USE extended_system_types,           ONLY: lnhc_parameters_type,&
      27             :                                               map_info_type,&
      28             :                                               npt_info_type
      29             :    USE input_constants,                 ONLY: &
      30             :         do_constr_atomic, do_constr_molec, do_region_defined, do_region_global, do_region_massive, &
      31             :         do_region_molecule, do_thermo_al, do_thermo_communication, do_thermo_csvr, do_thermo_gle, &
      32             :         do_thermo_no_communication, do_thermo_nose, isokin_ensemble, langevin_ensemble, &
      33             :         npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
      34             :         npt_f_ensemble, npt_i_ensemble, npt_ia_ensemble, nve_ensemble, nvt_adiabatic_ensemble, &
      35             :         nvt_ensemble, reftraj_ensemble
      36             :    USE input_section_types,             ONLY: section_vals_get,&
      37             :                                               section_vals_get_subs_vals,&
      38             :                                               section_vals_type,&
      39             :                                               section_vals_val_get
      40             :    USE kinds,                           ONLY: default_string_length,&
      41             :                                               dp
      42             :    USE machine,                         ONLY: m_flush
      43             :    USE message_passing,                 ONLY: mp_comm_type,&
      44             :                                               mp_para_env_type
      45             :    USE molecule_kind_types,             ONLY: get_molecule_kind,&
      46             :                                               get_molecule_kind_set,&
      47             :                                               molecule_kind_type
      48             :    USE molecule_list_types,             ONLY: molecule_list_type
      49             :    USE molecule_types,                  ONLY: get_molecule,&
      50             :                                               global_constraint_type,&
      51             :                                               molecule_type
      52             :    USE motion_utils,                    ONLY: rot_ana
      53             :    USE particle_list_types,             ONLY: particle_list_type
      54             :    USE particle_types,                  ONLY: particle_type
      55             :    USE physcon,                         ONLY: femtoseconds
      56             :    USE qmmm_types,                      ONLY: qmmm_env_type
      57             :    USE shell_potential_types,           ONLY: shell_kind_type
      58             :    USE simpar_types,                    ONLY: simpar_type
      59             :    USE thermostat_types,                ONLY: thermostat_info_type,&
      60             :                                               thermostat_type,&
      61             :                                               thermostats_type
      62             : #include "../../base/base_uses.f90"
      63             : 
      64             :    IMPLICIT NONE
      65             : 
      66             :    PRIVATE
      67             :    PUBLIC :: compute_degrees_of_freedom, &
      68             :              compute_nfree, &
      69             :              setup_thermostat_info, &
      70             :              setup_adiabatic_thermostat_info, &
      71             :              ke_region_baro, &
      72             :              ke_region_particles, &
      73             :              ke_region_shells, &
      74             :              vel_rescale_baro, &
      75             :              vel_rescale_particles, &
      76             :              vel_rescale_shells, &
      77             :              get_thermostat_energies, &
      78             :              get_nhc_energies, &
      79             :              get_kin_energies, &
      80             :              communication_thermo_low2, &
      81             :              print_thermostats_status, &
      82             :              momentum_region_particles
      83             : 
      84             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'thermostat_utils'
      85             : 
      86             : CONTAINS
      87             : 
      88             : ! **************************************************************************************************
      89             : !> \brief ...
      90             : !> \param cell ...
      91             : !> \param simpar ...
      92             : !> \param molecule_kind_set ...
      93             : !> \param print_section ...
      94             : !> \param particles ...
      95             : !> \param gci ...
      96             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
      97             : ! **************************************************************************************************
      98           0 :    SUBROUTINE compute_nfree(cell, simpar, molecule_kind_set, &
      99             :                             print_section, particles, gci)
     100             : 
     101             :       TYPE(cell_type), POINTER                           :: cell
     102             :       TYPE(simpar_type), POINTER                         :: simpar
     103             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     104             :       TYPE(section_vals_type), POINTER                   :: print_section
     105             :       TYPE(particle_list_type), POINTER                  :: particles
     106             :       TYPE(global_constraint_type), POINTER              :: gci
     107             : 
     108             :       INTEGER                                            :: natom, nconstraint_ext, nconstraint_int, &
     109             :                                                             nrestraints_int, rot_dof, &
     110             :                                                             roto_trasl_dof
     111             : 
     112             : ! Retrieve information on number of atoms, constraints (external and internal)
     113             : 
     114             :       CALL get_molecule_kind_set(molecule_kind_set=molecule_kind_set, &
     115           0 :                                  natom=natom, nconstraint=nconstraint_int, nrestraints=nrestraints_int)
     116             : 
     117             :       ! Compute degrees of freedom
     118             :       CALL rot_ana(particles%els, dof=roto_trasl_dof, rot_dof=rot_dof, &
     119             :                    print_section=print_section, keep_rotations=.FALSE., &
     120           0 :                    mass_weighted=.TRUE., natoms=natom)
     121             : 
     122           0 :       roto_trasl_dof = roto_trasl_dof - MIN(SUM(cell%perd(1:3)), rot_dof)
     123             : 
     124             :       ! Saving this value of simpar preliminar to the real count of constraints..
     125           0 :       simpar%nfree_rot_transl = roto_trasl_dof
     126             : 
     127             :       ! compute the total number of degrees of freedom for temperature
     128           0 :       nconstraint_ext = gci%ntot - gci%nrestraint
     129           0 :       simpar%nfree = 3*natom - nconstraint_int - nconstraint_ext - roto_trasl_dof
     130             : 
     131           0 :    END SUBROUTINE compute_nfree
     132             : 
     133             : ! **************************************************************************************************
     134             : !> \brief ...
     135             : !> \param thermostats ...
     136             : !> \param cell ...
     137             : !> \param simpar ...
     138             : !> \param molecule_kind_set ...
     139             : !> \param local_molecules ...
     140             : !> \param molecules ...
     141             : !> \param particles ...
     142             : !> \param print_section ...
     143             : !> \param region_sections ...
     144             : !> \param gci ...
     145             : !> \param region ...
     146             : !> \param qmmm_env ...
     147             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
     148             : ! **************************************************************************************************
     149        1796 :    SUBROUTINE compute_degrees_of_freedom(thermostats, cell, simpar, molecule_kind_set, &
     150             :                                          local_molecules, molecules, particles, print_section, region_sections, gci, &
     151             :                                          region, qmmm_env)
     152             : 
     153             :       TYPE(thermostats_type), POINTER                    :: thermostats
     154             :       TYPE(cell_type), POINTER                           :: cell
     155             :       TYPE(simpar_type), POINTER                         :: simpar
     156             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     157             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     158             :       TYPE(molecule_list_type), POINTER                  :: molecules
     159             :       TYPE(particle_list_type), POINTER                  :: particles
     160             :       TYPE(section_vals_type), POINTER                   :: print_section, region_sections
     161             :       TYPE(global_constraint_type), POINTER              :: gci
     162             :       INTEGER, INTENT(IN)                                :: region
     163             :       TYPE(qmmm_env_type), POINTER                       :: qmmm_env
     164             : 
     165             :       INTEGER                                            :: iw, natom, nconstraint_ext, &
     166             :                                                             nconstraint_int, nrestraints_int, &
     167             :                                                             rot_dof, roto_trasl_dof
     168             :       TYPE(cp_logger_type), POINTER                      :: logger
     169             : 
     170             : ! Retrieve information on number of atoms, constraints (external and internal)
     171             : 
     172             :       CALL get_molecule_kind_set(molecule_kind_set=molecule_kind_set, &
     173        1796 :                                  natom=natom, nconstraint=nconstraint_int, nrestraints=nrestraints_int)
     174             : 
     175             :       ! Compute degrees of freedom
     176             :       CALL rot_ana(particles%els, dof=roto_trasl_dof, rot_dof=rot_dof, &
     177             :                    print_section=print_section, keep_rotations=.FALSE., &
     178        1796 :                    mass_weighted=.TRUE., natoms=natom)
     179             : 
     180        7184 :       roto_trasl_dof = roto_trasl_dof - MIN(SUM(cell%perd(1:3)), rot_dof)
     181             : 
     182             :       ! Collect info about thermostats
     183             :       CALL setup_thermostat_info(thermostats%thermostat_info_part, molecule_kind_set, &
     184             :                                  local_molecules, molecules, particles, region, simpar%ensemble, roto_trasl_dof, &
     185        1796 :                                  region_sections=region_sections, qmmm_env=qmmm_env)
     186             : 
     187             :       ! Saving this value of simpar preliminar to the real count of constraints..
     188        1796 :       simpar%nfree_rot_transl = roto_trasl_dof
     189             : 
     190             :       ! compute the total number of degrees of freedom for temperature
     191        1796 :       nconstraint_ext = gci%ntot - gci%nrestraint
     192        1796 :       simpar%nfree = 3*natom - nconstraint_int - nconstraint_ext - roto_trasl_dof
     193             : 
     194        1796 :       logger => cp_get_default_logger()
     195             :       iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", &
     196        1796 :                                 extension=".log")
     197        1796 :       IF (iw > 0) THEN
     198             :          WRITE (iw, '(/,T2,A)') &
     199         811 :             'DOF| Calculation of degrees of freedom'
     200             :          WRITE (iw, '(T2,A,T71,I10)') &
     201         811 :             'DOF| Number of atoms', natom, &
     202         811 :             'DOF| Number of intramolecular constraints', nconstraint_int, &
     203         811 :             'DOF| Number of intermolecular constraints', nconstraint_ext, &
     204         811 :             'DOF| Invariants (translations + rotations)', roto_trasl_dof, &
     205        1622 :             'DOF| Degrees of freedom', simpar%nfree
     206             :          WRITE (iw, '(/,T2,A)') &
     207         811 :             'DOF| Restraints information'
     208             :          WRITE (iw, '(T2,A,T71,I10)') &
     209         811 :             'DOF| Number of intramolecular restraints', nrestraints_int, &
     210        1622 :             'DOF| Number of intermolecular restraints', gci%nrestraint
     211             :       END IF
     212             :       CALL cp_print_key_finished_output(iw, logger, print_section, &
     213        1796 :                                         "PROGRAM_RUN_INFO")
     214             : 
     215        1796 :    END SUBROUTINE compute_degrees_of_freedom
     216             : 
     217             : ! **************************************************************************************************
     218             : !> \brief ...
     219             : !> \param thermostat_info ...
     220             : !> \param molecule_kind_set ...
     221             : !> \param local_molecules ...
     222             : !> \param molecules ...
     223             : !> \param particles ...
     224             : !> \param region ...
     225             : !> \param ensemble ...
     226             : !> \param nfree ...
     227             : !> \param shell ...
     228             : !> \param region_sections ...
     229             : !> \param qmmm_env ...
     230             : !> \author 10.2011  CJM - PNNL
     231             : ! **************************************************************************************************
     232           0 :    SUBROUTINE setup_adiabatic_thermostat_info(thermostat_info, molecule_kind_set, local_molecules, &
     233             :                                               molecules, particles, region, ensemble, nfree, shell, region_sections, qmmm_env)
     234             :       TYPE(thermostat_info_type), POINTER                :: thermostat_info
     235             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     236             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     237             :       TYPE(molecule_list_type), POINTER                  :: molecules
     238             :       TYPE(particle_list_type), POINTER                  :: particles
     239             :       INTEGER, INTENT(IN)                                :: region, ensemble
     240             :       INTEGER, INTENT(INOUT), OPTIONAL                   :: nfree
     241             :       LOGICAL, INTENT(IN), OPTIONAL                      :: shell
     242             :       TYPE(section_vals_type), POINTER                   :: region_sections
     243             :       TYPE(qmmm_env_type), POINTER                       :: qmmm_env
     244             : 
     245             :       INTEGER :: dis_type, first_atom, i, ikind, imol, imol_global, ipart, itherm, katom, &
     246             :          last_atom, natom, natom_local, nkind, nmol_local, nmol_per_kind, nmolecule, nshell, &
     247             :          number, stat, sum_of_thermostats
     248           0 :       INTEGER, POINTER                                   :: molecule_list(:), thermolist(:)
     249             :       LOGICAL                                            :: check, do_shell, nointer, on_therm
     250             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     251           0 :       TYPE(molecule_type), POINTER                       :: molecule, molecule_set(:)
     252             : 
     253           0 :       NULLIFY (molecule_kind, molecule, thermostat_info%map_loc_thermo_gen, thermolist)
     254           0 :       nkind = SIZE(molecule_kind_set)
     255           0 :       do_shell = .FALSE.
     256           0 :       IF (PRESENT(shell)) do_shell = shell
     257             :       ! Counting the global number of thermostats
     258           0 :       sum_of_thermostats = 0
     259             :       ! Variable to denote independent thermostats (no communication necessary)
     260           0 :       nointer = .TRUE.
     261           0 :       check = .TRUE.
     262           0 :       number = 0
     263           0 :       dis_type = do_thermo_no_communication
     264             : 
     265             :       CALL get_adiabatic_region_info(region_sections, sum_of_thermostats, &
     266             :                                      thermolist=thermolist, &
     267             :                                      molecule_kind_set=molecule_kind_set, &
     268           0 :                                      molecules=molecules, particles=particles, qmmm_env=qmmm_env)
     269             : 
     270             : !    map_loc_thermo_gen=>thermostat_info%map_loc_thermo_gen
     271           0 :       molecule_set => molecules%els
     272           0 :       SELECT CASE (ensemble)
     273             :       CASE DEFAULT
     274           0 :          CPABORT('Unknown ensemble')
     275             :       CASE (nvt_adiabatic_ensemble)
     276           0 :          SELECT CASE (region)
     277             :          CASE (do_region_global)
     278             :             ! Global Thermostat
     279           0 :             nointer = .FALSE.
     280           0 :             sum_of_thermostats = 1
     281             :          CASE (do_region_molecule)
     282             :             ! Molecular Thermostat
     283             :             itherm = 0
     284           0 :             DO ikind = 1, nkind
     285           0 :                molecule_kind => molecule_kind_set(ikind)
     286           0 :                nmol_per_kind = local_molecules%n_el(ikind)
     287             :                CALL get_molecule_kind(molecule_kind, natom=natom, &
     288           0 :                                       molecule_list=molecule_list)
     289             : ! use thermolist ( ipart ) to get global indexing correct
     290           0 :                DO imol_global = 1, SIZE(molecule_list)
     291           0 :                   molecule => molecule_set(molecule_list(imol_global))
     292             :                   CALL get_molecule(molecule, first_atom=first_atom, &
     293           0 :                                     last_atom=last_atom)
     294           0 :                   on_therm = .TRUE.
     295           0 :                   DO katom = first_atom, last_atom
     296           0 :                      IF (thermolist(katom) == HUGE(0)) THEN
     297             :                         on_therm = .FALSE.
     298             :                         EXIT
     299             :                      END IF
     300             :                   END DO
     301           0 :                   IF (on_therm) THEN
     302           0 :                      itherm = itherm + 1
     303           0 :                      DO katom = first_atom, last_atom
     304           0 :                         thermolist(katom) = itherm
     305             :                      END DO
     306             :                   END IF
     307             :                END DO
     308             :             END DO
     309           0 :             DO i = 1, nkind
     310           0 :                molecule_kind => molecule_kind_set(i)
     311           0 :                CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, nshell=nshell)
     312           0 :                IF ((do_shell) .AND. (nshell == 0)) nmolecule = 0
     313           0 :                sum_of_thermostats = sum_of_thermostats + nmolecule
     314             :             END DO
     315             :             ! If we have ONE kind and ONE molecule, then effectively we have a GLOBAL thermostat
     316             :             ! and the degrees of freedom will be computed correctly for this special case
     317           0 :             IF ((nmolecule == 1) .AND. (nkind == 1)) nointer = .FALSE.
     318             :          CASE (do_region_massive)
     319             :             ! Massive Thermostat
     320           0 :             DO i = 1, nkind
     321           0 :                molecule_kind => molecule_kind_set(i)
     322             :                CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, &
     323           0 :                                       natom=natom, nshell=nshell)
     324           0 :                IF (do_shell) natom = nshell
     325           0 :                sum_of_thermostats = sum_of_thermostats + 3*natom*nmolecule
     326             :             END DO
     327             :          END SELECT
     328             : 
     329           0 :          natom_local = 0
     330           0 :          DO ikind = 1, SIZE(molecule_kind_set)
     331           0 :             nmol_per_kind = local_molecules%n_el(ikind)
     332           0 :             DO imol = 1, nmol_per_kind
     333           0 :                i = local_molecules%list(ikind)%array(imol)
     334           0 :                molecule => molecule_set(i)
     335           0 :                CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
     336           0 :                DO ipart = first_atom, last_atom
     337           0 :                   natom_local = natom_local + 1
     338             :                END DO
     339             :             END DO
     340             :          END DO
     341             : 
     342             :          ! Now map the local atoms with the corresponding thermostat
     343           0 :          ALLOCATE (thermostat_info%map_loc_thermo_gen(natom_local), stat=stat)
     344           0 :          thermostat_info%map_loc_thermo_gen = HUGE(0)
     345           0 :          CPASSERT(stat == 0)
     346           0 :          natom_local = 0
     347           0 :          DO ikind = 1, SIZE(molecule_kind_set)
     348           0 :             nmol_per_kind = local_molecules%n_el(ikind)
     349           0 :             DO imol = 1, nmol_per_kind
     350           0 :                i = local_molecules%list(ikind)%array(imol)
     351           0 :                molecule => molecule_set(i)
     352           0 :                CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
     353           0 :                DO ipart = first_atom, last_atom
     354           0 :                   natom_local = natom_local + 1
     355             : ! only map the correct region to the thermostat
     356           0 :                   IF (thermolist(ipart) /= HUGE(0)) &
     357           0 :                      thermostat_info%map_loc_thermo_gen(natom_local) = thermolist(ipart)
     358             :                END DO
     359             :             END DO
     360             :          END DO
     361             :          ! Here we decide which parallel algorithm to use.
     362             :          ! if there are only massive and molecule type thermostats we can use
     363             :          ! a local scheme, in cases involving any combination with a
     364             :          ! global thermostat we assume a coupling of  degrees of freedom
     365             :          ! from different processors
     366           0 :          IF (nointer) THEN
     367             :             ! Distributed thermostats, no interaction
     368           0 :             dis_type = do_thermo_no_communication
     369             :             ! we only count thermostats on this processor
     370             :             number = 0
     371           0 :             DO ikind = 1, nkind
     372           0 :                nmol_local = local_molecules%n_el(ikind)
     373           0 :                molecule_kind => molecule_kind_set(ikind)
     374           0 :                CALL get_molecule_kind(molecule_kind, natom=natom, nshell=nshell)
     375           0 :                IF (do_shell) THEN
     376           0 :                   natom = nshell
     377           0 :                   IF (nshell == 0) nmol_local = 0
     378             :                END IF
     379           0 :                IF (region == do_region_molecule) THEN
     380           0 :                   number = number + nmol_local
     381           0 :                ELSE IF (region == do_region_massive) THEN
     382           0 :                   number = number + 3*nmol_local*natom
     383             :                ELSE
     384           0 :                   CPABORT('Invalid region setup')
     385             :                END IF
     386             :             END DO
     387             :          ELSE
     388             :             ! REPlicated thermostats, INTERacting via communication
     389           0 :             dis_type = do_thermo_communication
     390           0 :             IF ((region == do_region_global) .OR. (region == do_region_molecule)) number = 1
     391             :          END IF
     392             : 
     393           0 :          IF (PRESENT(nfree)) THEN
     394             :             ! re-initializing simpar%nfree to zero because of multiple thermostats in the adiabatic sampling
     395           0 :             nfree = 0
     396             :          END IF
     397             :       END SELECT
     398             : 
     399             :       ! Saving information about thermostats
     400           0 :       thermostat_info%sum_of_thermostats = sum_of_thermostats
     401           0 :       thermostat_info%number_of_thermostats = number
     402           0 :       thermostat_info%dis_type = dis_type
     403             : 
     404           0 :       DEALLOCATE (thermolist)
     405             : 
     406           0 :    END SUBROUTINE setup_adiabatic_thermostat_info
     407             : 
     408             : ! **************************************************************************************************
     409             : !> \brief ...
     410             : !> \param region_sections ...
     411             : !> \param sum_of_thermostats ...
     412             : !> \param thermolist ...
     413             : !> \param molecule_kind_set ...
     414             : !> \param molecules ...
     415             : !> \param particles ...
     416             : !> \param qmmm_env ...
     417             : !> \author 10.2011 CJM -PNNL
     418             : ! **************************************************************************************************
     419           0 :    SUBROUTINE get_adiabatic_region_info(region_sections, sum_of_thermostats, &
     420             :                                         thermolist, molecule_kind_set, molecules, particles, &
     421             :                                         qmmm_env)
     422             :       TYPE(section_vals_type), POINTER                   :: region_sections
     423             :       INTEGER, INTENT(INOUT), OPTIONAL                   :: sum_of_thermostats
     424             :       INTEGER, DIMENSION(:), POINTER                     :: thermolist(:)
     425             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     426             :       TYPE(molecule_list_type), POINTER                  :: molecules
     427             :       TYPE(particle_list_type), POINTER                  :: particles
     428             :       TYPE(qmmm_env_type), POINTER                       :: qmmm_env
     429             : 
     430             :       CHARACTER(LEN=default_string_length), &
     431           0 :          DIMENSION(:), POINTER                           :: tmpstringlist
     432             :       INTEGER                                            :: first_atom, i, ig, ikind, ilist, imol, &
     433             :                                                             ipart, itherm, jg, last_atom, &
     434             :                                                             mregions, n_rep, nregions, output_unit
     435           0 :       INTEGER, DIMENSION(:), POINTER                     :: tmplist
     436             :       TYPE(cp_logger_type), POINTER                      :: logger
     437             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     438             :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     439             :       TYPE(molecule_type), POINTER                       :: molecule
     440             : 
     441           0 :       NULLIFY (tmplist, tmpstringlist, thermolist, molecule_kind, molecule, molecule_set)
     442           0 :       NULLIFY (logger)
     443           0 :       logger => cp_get_default_logger()
     444           0 :       output_unit = cp_logger_get_default_io_unit(logger)
     445             :       ! CPASSERT(.NOT.(ASSOCIATED(map_loc_thermo_gen)))
     446           0 :       CALL section_vals_get(region_sections, n_repetition=nregions)
     447           0 :       ALLOCATE (thermolist(particles%n_els))
     448           0 :       thermolist = HUGE(0)
     449           0 :       molecule_set => molecules%els
     450           0 :       mregions = nregions
     451           0 :       itherm = 0
     452           0 :       DO ig = 1, mregions
     453           0 :          CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ig, n_rep_val=n_rep)
     454           0 :          DO jg = 1, n_rep
     455           0 :             CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ig, i_rep_val=jg, i_vals=tmplist)
     456           0 :             DO i = 1, SIZE(tmplist)
     457           0 :                ipart = tmplist(i)
     458           0 :                CPASSERT(((ipart > 0) .AND. (ipart <= particles%n_els)))
     459           0 :                IF (thermolist(ipart) == HUGE(0)) THEN
     460           0 :                   itherm = itherm + 1
     461           0 :                   thermolist(ipart) = itherm
     462             :                ELSE
     463           0 :                   CPABORT("")
     464             :                END IF
     465             :             END DO
     466             :          END DO
     467           0 :          CALL section_vals_val_get(region_sections, "MOLNAME", i_rep_section=ig, n_rep_val=n_rep)
     468           0 :          DO jg = 1, n_rep
     469           0 :             CALL section_vals_val_get(region_sections, "MOLNAME", i_rep_section=ig, i_rep_val=jg, c_vals=tmpstringlist)
     470           0 :             DO ilist = 1, SIZE(tmpstringlist)
     471           0 :                DO ikind = 1, SIZE(molecule_kind_set)
     472           0 :                   molecule_kind => molecule_kind_set(ikind)
     473           0 :                   IF (molecule_kind%name == tmpstringlist(ilist)) THEN
     474           0 :                      DO imol = 1, SIZE(molecule_kind%molecule_list)
     475           0 :                         molecule => molecule_set(molecule_kind%molecule_list(imol))
     476           0 :                         CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
     477           0 :                         DO ipart = first_atom, last_atom
     478           0 :                            IF (thermolist(ipart) == HUGE(0)) THEN
     479           0 :                               itherm = itherm + 1
     480           0 :                               thermolist(ipart) = itherm
     481             :                            ELSE
     482           0 :                               CPABORT("")
     483             :                            END IF
     484             :                         END DO
     485             :                      END DO
     486             :                   END IF
     487             :                END DO
     488             :             END DO
     489             :          END DO
     490             :          CALL setup_thermostat_subsys(region_sections, qmmm_env, thermolist, molecule_set, &
     491           0 :                                       subsys_qm=.FALSE., ig=ig, sum_of_thermostats=sum_of_thermostats, nregions=nregions)
     492             :          CALL setup_thermostat_subsys(region_sections, qmmm_env, thermolist, molecule_set, &
     493           0 :                                       subsys_qm=.TRUE., ig=ig, sum_of_thermostats=sum_of_thermostats, nregions=nregions)
     494             :       END DO
     495             : 
     496           0 :       CPASSERT(.NOT. ALL(thermolist == HUGE(0)))
     497             : 
     498             : !    natom_local = 0
     499             : !    DO ikind = 1, SIZE(molecule_kind_set)
     500             : !       nmol_per_kind = local_molecules%n_el(ikind)
     501             : !       DO imol = 1, nmol_per_kind
     502             : !          i = local_molecules%list(ikind)%array(imol)
     503             : !          molecule => molecule_set(i)
     504             : !          CALL get_molecule ( molecule, first_atom = first_atom, last_atom = last_atom )
     505             : !          DO ipart = first_atom, last_atom
     506             : !             natom_local = natom_local + 1
     507             : !          END DO
     508             : !       END DO
     509             : !    END DO
     510             : 
     511             :       ! Now map the local atoms with the corresponding thermostat
     512             : !    ALLOCATE(map_loc_thermo_gen(natom_local),stat=stat)
     513             : !    map_loc_thermo_gen  = HUGE ( 0 )
     514             : !    CPPostcondition(stat==0,cp_failure_level,routineP,failure)
     515             : !    natom_local = 0
     516             : !    DO ikind = 1, SIZE(molecule_kind_set)
     517             : !       nmol_per_kind = local_molecules%n_el(ikind)
     518             : !       DO imol = 1, nmol_per_kind
     519             : !          i = local_molecules%list(ikind)%array(imol)
     520             : !          molecule => molecule_set(i)
     521             : !          CALL get_molecule ( molecule, first_atom = first_atom, last_atom = last_atom )
     522             : !          DO ipart = first_atom, last_atom
     523             : !             natom_local = natom_local + 1
     524             : ! only map the correct region to the thermostat
     525             : !             IF ( thermolist (ipart ) /= HUGE ( 0 ) ) &
     526             : !             map_loc_thermo_gen(natom_local) = thermolist(ipart)
     527             : !          END DO
     528             : !       END DO
     529             : !    END DO
     530             : 
     531             : !    DEALLOCATE(thermolist, stat=stat)
     532             : !    CPPostcondition(stat==0,cp_failure_level,routineP,failure)
     533           0 :    END SUBROUTINE get_adiabatic_region_info
     534             : ! **************************************************************************************************
     535             : !> \brief ...
     536             : !> \param thermostat_info ...
     537             : !> \param molecule_kind_set ...
     538             : !> \param local_molecules ...
     539             : !> \param molecules ...
     540             : !> \param particles ...
     541             : !> \param region ...
     542             : !> \param ensemble ...
     543             : !> \param nfree ...
     544             : !> \param shell ...
     545             : !> \param region_sections ...
     546             : !> \param qmmm_env ...
     547             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
     548             : ! **************************************************************************************************
     549        1842 :    SUBROUTINE setup_thermostat_info(thermostat_info, molecule_kind_set, local_molecules, &
     550             :                                     molecules, particles, region, ensemble, nfree, shell, region_sections, qmmm_env)
     551             :       TYPE(thermostat_info_type), POINTER                :: thermostat_info
     552             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     553             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     554             :       TYPE(molecule_list_type), POINTER                  :: molecules
     555             :       TYPE(particle_list_type), POINTER                  :: particles
     556             :       INTEGER, INTENT(IN)                                :: region, ensemble
     557             :       INTEGER, INTENT(INOUT), OPTIONAL                   :: nfree
     558             :       LOGICAL, INTENT(IN), OPTIONAL                      :: shell
     559             :       TYPE(section_vals_type), POINTER                   :: region_sections
     560             :       TYPE(qmmm_env_type), POINTER                       :: qmmm_env
     561             : 
     562             :       INTEGER                                            :: dis_type, i, ikind, natom, nkind, &
     563             :                                                             nmol_local, nmolecule, nshell, number, &
     564             :                                                             sum_of_thermostats
     565             :       LOGICAL                                            :: check, do_shell, nointer
     566             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     567             : 
     568        1842 :       NULLIFY (molecule_kind)
     569        1842 :       nkind = SIZE(molecule_kind_set)
     570        1842 :       do_shell = .FALSE.
     571        1842 :       IF (PRESENT(shell)) do_shell = shell
     572             :       ! Counting the global number of thermostats
     573        1842 :       sum_of_thermostats = 0
     574             :       ! Variable to denote independent thermostats (no communication necessary)
     575        1842 :       nointer = .TRUE.
     576        1842 :       check = .TRUE.
     577        1842 :       number = 0
     578        1842 :       dis_type = do_thermo_no_communication
     579             : 
     580        1842 :       SELECT CASE (ensemble)
     581             :       CASE DEFAULT
     582           0 :          CPABORT('Unknown ensemble')
     583             :       CASE (isokin_ensemble, nph_uniaxial_ensemble, nph_uniaxial_damped_ensemble, &
     584             :             reftraj_ensemble, langevin_ensemble)
     585             :          ! Do Nothing
     586             :       CASE (nve_ensemble, nvt_ensemble, nvt_adiabatic_ensemble, npt_i_ensemble, &
     587             :             npt_f_ensemble, npe_i_ensemble, npe_f_ensemble, npt_ia_ensemble)
     588        1754 :          IF (ensemble == nve_ensemble) check = do_shell
     589        3016 :          IF (check) THEN
     590         584 :             SELECT CASE (region)
     591             :             CASE (do_region_global)
     592             :                ! Global Thermostat
     593         284 :                nointer = .FALSE.
     594         284 :                sum_of_thermostats = 1
     595             :             CASE (do_region_molecule)
     596             :                ! Molecular Thermostat
     597        6052 :                DO i = 1, nkind
     598        5916 :                   molecule_kind => molecule_kind_set(i)
     599        5916 :                   CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, nshell=nshell)
     600        5916 :                   IF ((do_shell) .AND. (nshell == 0)) nmolecule = 0
     601       11968 :                   sum_of_thermostats = sum_of_thermostats + nmolecule
     602             :                END DO
     603             :                ! If we have ONE kind and ONE molecule, then effectively we have a GLOBAL thermostat
     604             :                ! and the degrees of freedom will be computed correctly for this special case
     605         136 :                IF ((nmolecule == 1) .AND. (nkind == 1)) nointer = .FALSE.
     606             :             CASE (do_region_massive)
     607             :                ! Massive Thermostat
     608        8882 :                DO i = 1, nkind
     609        8750 :                   molecule_kind => molecule_kind_set(i)
     610             :                   CALL get_molecule_kind(molecule_kind, nmolecule=nmolecule, &
     611        8750 :                                          natom=natom, nshell=nshell)
     612        8750 :                   IF (do_shell) natom = nshell
     613       17632 :                   sum_of_thermostats = sum_of_thermostats + 3*natom*nmolecule
     614             :                END DO
     615             :             CASE (do_region_defined)
     616             :                ! User defined region to thermostat..
     617          32 :                nointer = .FALSE.
     618             :                ! Determine the number of thermostats defined in the input
     619          32 :                CALL section_vals_get(region_sections, n_repetition=sum_of_thermostats)
     620          32 :                IF (sum_of_thermostats < 1) &
     621             :                   CALL cp_abort(__LOCATION__, &
     622          32 :                                 "Provide at least 1 region (&DEFINE_REGION) when using the thermostat type DEFINED")
     623             :             END SELECT
     624             : 
     625             :             ! Here we decide which parallel algorithm to use.
     626             :             ! if there are only massive and molecule type thermostats we can use
     627             :             ! a local scheme, in cases involving any combination with a
     628             :             ! global thermostat we assume a coupling of  degrees of freedom
     629             :             ! from different processors
     630             :             IF (nointer) THEN
     631             :                ! Distributed thermostats, no interaction
     632       14926 :                dis_type = do_thermo_no_communication
     633             :                ! we only count thermostats on this processor
     634             :                number = 0
     635       14926 :                DO ikind = 1, nkind
     636       14662 :                   nmol_local = local_molecules%n_el(ikind)
     637       14662 :                   molecule_kind => molecule_kind_set(ikind)
     638       14662 :                   CALL get_molecule_kind(molecule_kind, natom=natom, nshell=nshell)
     639       14662 :                   IF (do_shell) THEN
     640          58 :                      natom = nshell
     641          58 :                      IF (nshell == 0) nmol_local = 0
     642             :                   END IF
     643       29588 :                   IF (region == do_region_molecule) THEN
     644        5912 :                      number = number + nmol_local
     645        8750 :                   ELSE IF (region == do_region_massive) THEN
     646        8750 :                      number = number + 3*nmol_local*natom
     647             :                   ELSE
     648           0 :                      CPABORT('Invalid region setup')
     649             :                   END IF
     650             :                END DO
     651             :             ELSE
     652             :                ! REPlicated thermostats, INTERacting via communication
     653         320 :                dis_type = do_thermo_communication
     654         320 :                IF ((region == do_region_global) .OR. (region == do_region_molecule)) THEN
     655         288 :                   number = 1
     656          32 :                ELSE IF (region == do_region_defined) THEN
     657             :                   CALL get_defined_region_info(region_sections, number, sum_of_thermostats, &
     658             :                                                map_loc_thermo_gen=thermostat_info%map_loc_thermo_gen, &
     659             :                                                local_molecules=local_molecules, molecule_kind_set=molecule_kind_set, &
     660          32 :                                                molecules=molecules, particles=particles, qmmm_env=qmmm_env)
     661             :                END IF
     662             :             END IF
     663             : 
     664         584 :             IF (PRESENT(nfree)) THEN
     665         538 :                IF ((sum_of_thermostats > 1) .OR. (dis_type == do_thermo_no_communication)) THEN
     666             :                   ! re-initializing simpar%nfree to zero because of multiple thermostats
     667         260 :                   nfree = 0
     668             :                END IF
     669             :             END IF
     670             :          END IF
     671             :       END SELECT
     672             : 
     673             :       ! Saving information about thermostats
     674        1842 :       thermostat_info%sum_of_thermostats = sum_of_thermostats
     675        1842 :       thermostat_info%number_of_thermostats = number
     676        1842 :       thermostat_info%dis_type = dis_type
     677        1842 :    END SUBROUTINE setup_thermostat_info
     678             : 
     679             : ! **************************************************************************************************
     680             : !> \brief ...
     681             : !> \param region_sections ...
     682             : !> \param number ...
     683             : !> \param sum_of_thermostats ...
     684             : !> \param map_loc_thermo_gen ...
     685             : !> \param local_molecules ...
     686             : !> \param molecule_kind_set ...
     687             : !> \param molecules ...
     688             : !> \param particles ...
     689             : !> \param qmmm_env ...
     690             : !> \author 11.2007 [tlaino] - Teodoro Laino - University of Zurich
     691             : ! **************************************************************************************************
     692          32 :    SUBROUTINE get_defined_region_info(region_sections, number, sum_of_thermostats, &
     693             :                                       map_loc_thermo_gen, local_molecules, molecule_kind_set, molecules, particles, &
     694             :                                       qmmm_env)
     695             :       TYPE(section_vals_type), POINTER                   :: region_sections
     696             :       INTEGER, INTENT(OUT), OPTIONAL                     :: number
     697             :       INTEGER, INTENT(INOUT), OPTIONAL                   :: sum_of_thermostats
     698             :       INTEGER, DIMENSION(:), POINTER                     :: map_loc_thermo_gen
     699             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     700             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     701             :       TYPE(molecule_list_type), POINTER                  :: molecules
     702             :       TYPE(particle_list_type), POINTER                  :: particles
     703             :       TYPE(qmmm_env_type), POINTER                       :: qmmm_env
     704             : 
     705             :       CHARACTER(LEN=default_string_length), &
     706          32 :          DIMENSION(:), POINTER                           :: tmpstringlist
     707             :       INTEGER :: first_atom, i, ig, ikind, ilist, imol, ipart, jg, last_atom, mregions, n_rep, &
     708             :          natom_local, nmol_per_kind, nregions, output_unit
     709          32 :       INTEGER, DIMENSION(:), POINTER                     :: thermolist, tmp, tmplist
     710             :       TYPE(cp_logger_type), POINTER                      :: logger
     711             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     712             :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     713             :       TYPE(molecule_type), POINTER                       :: molecule
     714             : 
     715          32 :       NULLIFY (tmplist, tmpstringlist, thermolist, molecule_kind, molecule, molecule_set)
     716          32 :       NULLIFY (logger)
     717          64 :       logger => cp_get_default_logger()
     718          32 :       output_unit = cp_logger_get_default_io_unit(logger)
     719          32 :       CPASSERT(.NOT. (ASSOCIATED(map_loc_thermo_gen)))
     720          32 :       CALL section_vals_get(region_sections, n_repetition=nregions)
     721          96 :       ALLOCATE (thermolist(particles%n_els))
     722       43970 :       thermolist = HUGE(0)
     723          32 :       molecule_set => molecules%els
     724          32 :       mregions = nregions
     725         102 :       DO ig = 1, mregions
     726          70 :          CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ig, n_rep_val=n_rep)
     727         184 :          DO jg = 1, n_rep
     728         114 :             CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ig, i_rep_val=jg, i_vals=tmplist)
     729        2428 :             DO i = 1, SIZE(tmplist)
     730        2244 :                ipart = tmplist(i)
     731        2244 :                CPASSERT(((ipart > 0) .AND. (ipart <= particles%n_els)))
     732        2358 :                IF (thermolist(ipart) == HUGE(0)) THEN
     733        2244 :                   thermolist(ipart) = ig
     734             :                ELSE
     735           0 :                   CPABORT("")
     736             :                END IF
     737             :             END DO
     738             :          END DO
     739          70 :          CALL section_vals_val_get(region_sections, "MOLNAME", i_rep_section=ig, n_rep_val=n_rep)
     740          74 :          DO jg = 1, n_rep
     741           4 :             CALL section_vals_val_get(region_sections, "MOLNAME", i_rep_section=ig, i_rep_val=jg, c_vals=tmpstringlist)
     742          78 :             DO ilist = 1, SIZE(tmpstringlist)
     743          20 :                DO ikind = 1, SIZE(molecule_kind_set)
     744          12 :                   molecule_kind => molecule_kind_set(ikind)
     745          16 :                   IF (molecule_kind%name == tmpstringlist(ilist)) THEN
     746          48 :                      DO imol = 1, SIZE(molecule_kind%molecule_list)
     747          44 :                         molecule => molecule_set(molecule_kind%molecule_list(imol))
     748          44 :                         CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
     749         180 :                         DO ipart = first_atom, last_atom
     750         176 :                            IF (thermolist(ipart) == HUGE(0)) THEN
     751         132 :                               thermolist(ipart) = ig
     752             :                            ELSE
     753           0 :                               CPABORT("")
     754             :                            END IF
     755             :                         END DO
     756             :                      END DO
     757             :                   END IF
     758             :                END DO
     759             :             END DO
     760             :          END DO
     761             :          CALL setup_thermostat_subsys(region_sections, qmmm_env, thermolist, molecule_set, &
     762          70 :                                       subsys_qm=.FALSE., ig=ig, sum_of_thermostats=sum_of_thermostats, nregions=nregions)
     763             :          CALL setup_thermostat_subsys(region_sections, qmmm_env, thermolist, molecule_set, &
     764         242 :                                       subsys_qm=.TRUE., ig=ig, sum_of_thermostats=sum_of_thermostats, nregions=nregions)
     765             :       END DO
     766             : 
     767             :       ! Dump IO warning for not thermalized particles
     768       29100 :       IF (ANY(thermolist == HUGE(0))) THEN
     769          14 :          nregions = nregions + 1
     770          14 :          sum_of_thermostats = sum_of_thermostats + 1
     771       15008 :          ALLOCATE (tmp(COUNT(thermolist == HUGE(0))))
     772          14 :          ilist = 0
     773       14980 :          DO i = 1, SIZE(thermolist)
     774       14980 :             IF (thermolist(i) == HUGE(0)) THEN
     775       13894 :                ilist = ilist + 1
     776       13894 :                tmp(ilist) = i
     777       13894 :                thermolist(i) = nregions
     778             :             END IF
     779             :          END DO
     780          14 :          IF (output_unit > 0) THEN
     781           7 :             WRITE (output_unit, '(A)') " WARNING| No thermostats defined for the following atoms:"
     782        6954 :             IF (ilist > 0) WRITE (output_unit, '(" WARNING|",9I8)') tmp
     783           7 :             WRITE (output_unit, '(A)') " WARNING| They will be included in a further unique thermostat!"
     784             :          END IF
     785          14 :          DEALLOCATE (tmp)
     786             :       END IF
     787       43970 :       CPASSERT(ALL(thermolist /= HUGE(0)))
     788             : 
     789             :       ! Now identify the local number of thermostats
     790          96 :       ALLOCATE (tmp(nregions))
     791         116 :       tmp = 0
     792          32 :       natom_local = 0
     793          98 :       DO ikind = 1, SIZE(molecule_kind_set)
     794          66 :          nmol_per_kind = local_molecules%n_el(ikind)
     795        7384 :          DO imol = 1, nmol_per_kind
     796        7286 :             i = local_molecules%list(ikind)%array(imol)
     797        7286 :             molecule => molecule_set(i)
     798        7286 :             CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
     799       29321 :             DO ipart = first_atom, last_atom
     800       21969 :                natom_local = natom_local + 1
     801       29255 :                tmp(thermolist(ipart)) = 1
     802             :             END DO
     803             :          END DO
     804             :       END DO
     805         116 :       number = SUM(tmp)
     806          32 :       DEALLOCATE (tmp)
     807             : 
     808             :       ! Now map the local atoms with the corresponding thermostat
     809          96 :       ALLOCATE (map_loc_thermo_gen(natom_local))
     810          32 :       natom_local = 0
     811          98 :       DO ikind = 1, SIZE(molecule_kind_set)
     812          66 :          nmol_per_kind = local_molecules%n_el(ikind)
     813        7384 :          DO imol = 1, nmol_per_kind
     814        7286 :             i = local_molecules%list(ikind)%array(imol)
     815        7286 :             molecule => molecule_set(i)
     816        7286 :             CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
     817       29321 :             DO ipart = first_atom, last_atom
     818       21969 :                natom_local = natom_local + 1
     819       29255 :                map_loc_thermo_gen(natom_local) = thermolist(ipart)
     820             :             END DO
     821             :          END DO
     822             :       END DO
     823             : 
     824          32 :       DEALLOCATE (thermolist)
     825          64 :    END SUBROUTINE get_defined_region_info
     826             : 
     827             : ! **************************************************************************************************
     828             : !> \brief ...
     829             : !> \param region_sections ...
     830             : !> \param qmmm_env ...
     831             : !> \param thermolist ...
     832             : !> \param molecule_set ...
     833             : !> \param subsys_qm ...
     834             : !> \param ig ...
     835             : !> \param sum_of_thermostats ...
     836             : !> \param nregions ...
     837             : !> \author 11.2007 [tlaino] - Teodoro Laino - University of Zurich
     838             : ! **************************************************************************************************
     839         140 :    SUBROUTINE setup_thermostat_subsys(region_sections, qmmm_env, thermolist, &
     840             :                                       molecule_set, subsys_qm, ig, sum_of_thermostats, nregions)
     841             :       TYPE(section_vals_type), POINTER                   :: region_sections
     842             :       TYPE(qmmm_env_type), POINTER                       :: qmmm_env
     843             :       INTEGER, DIMENSION(:), POINTER                     :: thermolist
     844             :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     845             :       LOGICAL, INTENT(IN)                                :: subsys_qm
     846             :       INTEGER, INTENT(IN)                                :: ig
     847             :       INTEGER, INTENT(INOUT)                             :: sum_of_thermostats, nregions
     848             : 
     849             :       CHARACTER(LEN=default_string_length)               :: label1, label2
     850             :       INTEGER                                            :: first_atom, i, imolecule, ipart, &
     851             :                                                             last_atom, nrep, thermo1
     852         140 :       INTEGER, DIMENSION(:), POINTER                     :: atom_index1
     853             :       LOGICAL                                            :: explicit
     854             :       TYPE(molecule_type), POINTER                       :: molecule
     855             : 
     856         140 :       label1 = "MM_SUBSYS"
     857         140 :       label2 = "QM_SUBSYS"
     858         140 :       IF (subsys_qm) THEN
     859          70 :          label1 = "QM_SUBSYS"
     860          70 :          label2 = "MM_SUBSYS"
     861             :       END IF
     862             :       CALL section_vals_val_get(region_sections, TRIM(label1), i_rep_section=ig, &
     863         140 :                                 n_rep_val=nrep, explicit=explicit)
     864         140 :       IF (nrep == 1 .AND. explicit) THEN
     865           8 :          IF (ASSOCIATED(qmmm_env)) THEN
     866           8 :             atom_index1 => qmmm_env%qm%mm_atom_index
     867           8 :             IF (subsys_qm) THEN
     868           4 :                atom_index1 => qmmm_env%qm%qm_atom_index
     869             :             END IF
     870           8 :             CALL section_vals_val_get(region_sections, TRIM(label1), i_val=thermo1, i_rep_section=ig)
     871           4 :             SELECT CASE (thermo1)
     872             :             CASE (do_constr_atomic)
     873       13820 :                DO i = 1, SIZE(atom_index1)
     874       13816 :                   ipart = atom_index1(i)
     875       13816 :                   IF (subsys_qm .AND. qmmm_env%qm%qmmm_link .AND. ASSOCIATED(qmmm_env%qm%mm_link_atoms)) THEN
     876          46 :                      IF (ANY(ipart == qmmm_env%qm%mm_link_atoms)) CYCLE
     877             :                   END IF
     878       13818 :                   IF (thermolist(ipart) == HUGE(0)) THEN
     879       13814 :                      thermolist(ipart) = ig
     880             :                   ELSE
     881             :                      CALL cp_abort(__LOCATION__, &
     882             :                                    'One atom ('//cp_to_string(ipart)//') of the '// &
     883             :                                    TRIM(label1)//' was already assigned to'// &
     884             :                                    ' the thermostatting region Nr.'//cp_to_string(thermolist(ipart))// &
     885           0 :                                    '. Please check the input for inconsistencies!')
     886             :                   END IF
     887             :                END DO
     888             :             CASE (do_constr_molec)
     889        9168 :                DO imolecule = 1, SIZE(molecule_set)
     890        9160 :                   molecule => molecule_set(imolecule)
     891        9160 :                   CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
     892    15908454 :                   IF (ANY(atom_index1 >= first_atom .AND. atom_index1 <= last_atom)) THEN
     893       18436 :                      DO ipart = first_atom, last_atom
     894       18436 :                         IF (thermolist(ipart) == HUGE(0)) THEN
     895       13854 :                            thermolist(ipart) = ig
     896             :                         ELSE
     897             :                            CALL cp_abort(__LOCATION__, &
     898             :                                          'One atom ('//cp_to_string(ipart)//') of the '// &
     899             :                                          TRIM(label1)//' was already assigned to'// &
     900             :                                          ' the thermostatting region Nr.'//cp_to_string(thermolist(ipart))// &
     901           0 :                                          '. Please check the input for inconsistencies!')
     902             :                         END IF
     903             :                      END DO
     904             :                   END IF
     905             :                END DO
     906             :             END SELECT
     907             :          ELSE
     908           0 :             sum_of_thermostats = sum_of_thermostats - 1
     909           0 :             nregions = nregions - 1
     910             :          END IF
     911             :       END IF
     912         140 :    END SUBROUTINE setup_thermostat_subsys
     913             : 
     914             : ! **************************************************************************************************
     915             : !> \brief ...
     916             : !> \param map_info ...
     917             : !> \param npt ...
     918             : !> \param group ...
     919             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
     920             : ! **************************************************************************************************
     921        5444 :    SUBROUTINE ke_region_baro(map_info, npt, group)
     922             :       TYPE(map_info_type), POINTER                       :: map_info
     923             :       TYPE(npt_info_type), DIMENSION(:, :), &
     924             :          INTENT(INOUT)                                   :: npt
     925             :       TYPE(mp_comm_type), INTENT(IN)                     :: group
     926             : 
     927             :       INTEGER                                            :: i, j, ncoef
     928             : 
     929       10888 :       map_info%v_scale = 1.0_dp
     930       10888 :       map_info%s_kin = 0.0_dp
     931        5444 :       ncoef = 0
     932       13672 :       DO i = 1, SIZE(npt, 1)
     933       30252 :          DO j = 1, SIZE(npt, 2)
     934       16580 :             ncoef = ncoef + 1
     935             :             map_info%p_kin(1, ncoef)%point = map_info%p_kin(1, ncoef)%point &
     936       24808 :                                              + npt(i, j)%mass*npt(i, j)%v**2
     937             :          END DO
     938             :       END DO
     939             : 
     940        5444 :       IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
     941             : 
     942        5444 :    END SUBROUTINE ke_region_baro
     943             : 
     944             : ! **************************************************************************************************
     945             : !> \brief ...
     946             : !> \param map_info ...
     947             : !> \param npt ...
     948             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
     949             : ! **************************************************************************************************
     950        4360 :    SUBROUTINE vel_rescale_baro(map_info, npt)
     951             :       TYPE(map_info_type), POINTER                       :: map_info
     952             :       TYPE(npt_info_type), DIMENSION(:, :), &
     953             :          INTENT(INOUT)                                   :: npt
     954             : 
     955             :       INTEGER                                            :: i, j, ncoef
     956             : 
     957        4360 :       ncoef = 0
     958       11344 :       DO i = 1, SIZE(npt, 1)
     959       26200 :          DO j = 1, SIZE(npt, 2)
     960       14856 :             ncoef = ncoef + 1
     961       21840 :             npt(i, j)%v = npt(i, j)%v*map_info%p_scale(1, ncoef)%point
     962             :          END DO
     963             :       END DO
     964             : 
     965        4360 :    END SUBROUTINE vel_rescale_baro
     966             : 
     967             : ! **************************************************************************************************
     968             : !> \brief ...
     969             : !> \param map_info ...
     970             : !> \param particle_set ...
     971             : !> \param molecule_kind_set ...
     972             : !> \param local_molecules ...
     973             : !> \param molecule_set ...
     974             : !> \param group ...
     975             : !> \param vel ...
     976             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
     977             : ! **************************************************************************************************
     978       24884 :    SUBROUTINE ke_region_particles(map_info, particle_set, molecule_kind_set, &
     979       24884 :                                   local_molecules, molecule_set, group, vel)
     980             : 
     981             :       TYPE(map_info_type), POINTER                       :: map_info
     982             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     983             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     984             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     985             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     986             :       TYPE(mp_comm_type), INTENT(IN)                     :: group
     987             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: vel(:, :)
     988             : 
     989             :       INTEGER                                            :: first_atom, ii, ikind, imol, imol_local, &
     990             :                                                             ipart, last_atom, nmol_local
     991             :       LOGICAL                                            :: present_vel
     992             :       REAL(KIND=dp)                                      :: mass
     993             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     994             :       TYPE(molecule_type), POINTER                       :: molecule
     995             : 
     996     1551928 :       map_info%v_scale = 1.0_dp
     997     1551928 :       map_info%s_kin = 0.0_dp
     998       24884 :       present_vel = PRESENT(vel)
     999       24884 :       ii = 0
    1000     1116680 :       DO ikind = 1, SIZE(molecule_kind_set)
    1001     1091796 :          nmol_local = local_molecules%n_el(ikind)
    1002     2425634 :          DO imol_local = 1, nmol_local
    1003     1308954 :             imol = local_molecules%list(ikind)%array(imol_local)
    1004     1308954 :             molecule => molecule_set(imol)
    1005     1308954 :             CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
    1006     5370552 :             DO ipart = first_atom, last_atom
    1007     2969802 :                ii = ii + 1
    1008     2969802 :                atomic_kind => particle_set(ipart)%atomic_kind
    1009     2969802 :                CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
    1010     4278756 :                IF (present_vel) THEN
    1011     1484901 :                   IF (ASSOCIATED(map_info%p_kin(1, ii)%point)) &
    1012     1484901 :                      map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mass*vel(1, ipart)**2
    1013     1484901 :                   IF (ASSOCIATED(map_info%p_kin(2, ii)%point)) &
    1014     1484901 :                      map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mass*vel(2, ipart)**2
    1015     1484901 :                   IF (ASSOCIATED(map_info%p_kin(3, ii)%point)) &
    1016     1484901 :                      map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mass*vel(3, ipart)**2
    1017             :                ELSE
    1018     1484901 :                   IF (ASSOCIATED(map_info%p_kin(1, ii)%point)) &
    1019     1484901 :                      map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mass*particle_set(ipart)%v(1)**2
    1020     1484901 :                   IF (ASSOCIATED(map_info%p_kin(2, ii)%point)) &
    1021     1484901 :                      map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mass*particle_set(ipart)%v(2)**2
    1022     1484901 :                   IF (ASSOCIATED(map_info%p_kin(3, ii)%point)) &
    1023     1484901 :                      map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mass*particle_set(ipart)%v(3)**2
    1024             :                END IF
    1025             :             END DO
    1026             :          END DO
    1027             :       END DO
    1028             : 
    1029       60260 :       IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
    1030             : 
    1031       24884 :    END SUBROUTINE ke_region_particles
    1032             : 
    1033             : ! **************************************************************************************************
    1034             : !> \brief ...
    1035             : !> \param map_info ...
    1036             : !> \param particle_set ...
    1037             : !> \param molecule_kind_set ...
    1038             : !> \param local_molecules ...
    1039             : !> \param molecule_set ...
    1040             : !> \param group ...
    1041             : !> \param vel ...
    1042             : !> \author 07.2009 MI
    1043             : ! **************************************************************************************************
    1044         800 :    SUBROUTINE momentum_region_particles(map_info, particle_set, molecule_kind_set, &
    1045         800 :                                         local_molecules, molecule_set, group, vel)
    1046             : 
    1047             :       TYPE(map_info_type), POINTER                       :: map_info
    1048             :       TYPE(particle_type), POINTER                       :: particle_set(:)
    1049             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
    1050             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
    1051             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
    1052             :       TYPE(mp_comm_type), INTENT(IN)                     :: group
    1053             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: vel(:, :)
    1054             : 
    1055             :       INTEGER                                            :: first_atom, ii, ikind, imol, imol_local, &
    1056             :                                                             ipart, last_atom, nmol_local
    1057             :       LOGICAL                                            :: present_vel
    1058             :       REAL(KIND=dp)                                      :: mass
    1059             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1060             :       TYPE(molecule_type), POINTER                       :: molecule
    1061             : 
    1062      130400 :       map_info%v_scale = 1.0_dp
    1063      130400 :       map_info%s_kin = 0.0_dp
    1064         800 :       present_vel = PRESENT(vel)
    1065         800 :       ii = 0
    1066       87200 :       DO ikind = 1, SIZE(molecule_kind_set)
    1067       86400 :          nmol_local = local_molecules%n_el(ikind)
    1068      130400 :          DO imol_local = 1, nmol_local
    1069       43200 :             imol = local_molecules%list(ikind)%array(imol_local)
    1070       43200 :             molecule => molecule_set(imol)
    1071       43200 :             CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
    1072      172800 :             DO ipart = first_atom, last_atom
    1073       43200 :                ii = ii + 1
    1074       43200 :                atomic_kind => particle_set(ipart)%atomic_kind
    1075       43200 :                CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
    1076       86400 :                IF (present_vel) THEN
    1077       21600 :                   map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + SQRT(mass)*vel(1, ipart)
    1078       21600 :                   map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + SQRT(mass)*vel(2, ipart)
    1079       21600 :                   map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + SQRT(mass)*vel(3, ipart)
    1080             :                ELSE
    1081       21600 :                   map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + SQRT(mass)*particle_set(ipart)%v(1)
    1082       21600 :                   map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + SQRT(mass)*particle_set(ipart)%v(2)
    1083       21600 :                   map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + SQRT(mass)*particle_set(ipart)%v(3)
    1084             :                END IF
    1085             :             END DO
    1086             :          END DO
    1087             :       END DO
    1088             : 
    1089         800 :       IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
    1090             : 
    1091         800 :    END SUBROUTINE momentum_region_particles
    1092             : 
    1093             : ! **************************************************************************************************
    1094             : !> \brief ...
    1095             : !> \param map_info ...
    1096             : !> \param molecule_kind_set ...
    1097             : !> \param molecule_set ...
    1098             : !> \param particle_set ...
    1099             : !> \param local_molecules ...
    1100             : !> \param shell_adiabatic ...
    1101             : !> \param shell_particle_set ...
    1102             : !> \param core_particle_set ...
    1103             : !> \param vel ...
    1104             : !> \param shell_vel ...
    1105             : !> \param core_vel ...
    1106             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
    1107             : ! **************************************************************************************************
    1108       19056 :    SUBROUTINE vel_rescale_particles(map_info, molecule_kind_set, molecule_set, &
    1109             :                                     particle_set, local_molecules, shell_adiabatic, shell_particle_set, &
    1110       19056 :                                     core_particle_set, vel, shell_vel, core_vel)
    1111             : 
    1112             :       TYPE(map_info_type), POINTER                       :: map_info
    1113             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
    1114             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
    1115             :       TYPE(particle_type), POINTER                       :: particle_set(:)
    1116             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
    1117             :       LOGICAL, INTENT(IN)                                :: shell_adiabatic
    1118             :       TYPE(particle_type), OPTIONAL, POINTER             :: shell_particle_set(:), &
    1119             :                                                             core_particle_set(:)
    1120             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: vel(:, :), shell_vel(:, :), &
    1121             :                                                             core_vel(:, :)
    1122             : 
    1123             :       INTEGER                                            :: first_atom, ii, ikind, imol, imol_local, &
    1124             :                                                             ipart, jj, last_atom, nmol_local, &
    1125             :                                                             shell_index
    1126             :       LOGICAL                                            :: present_vel
    1127             :       REAL(KIND=dp)                                      :: fac_massc, fac_masss, mass, vc(3), vs(3)
    1128             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1129             :       TYPE(molecule_type), POINTER                       :: molecule
    1130             :       TYPE(shell_kind_type), POINTER                     :: shell
    1131             : 
    1132       19056 :       ii = 0
    1133       19056 :       jj = 0
    1134       19056 :       present_vel = PRESENT(vel)
    1135             :       ! Just few checks for consistency
    1136       19056 :       IF (present_vel) THEN
    1137        9528 :          IF (shell_adiabatic) THEN
    1138        1410 :             CPASSERT(PRESENT(shell_vel))
    1139        1410 :             CPASSERT(PRESENT(core_vel))
    1140             :          END IF
    1141             :       ELSE
    1142        9528 :          IF (shell_adiabatic) THEN
    1143        1410 :             CPASSERT(PRESENT(shell_particle_set))
    1144        1410 :             CPASSERT(PRESENT(core_particle_set))
    1145             :          END IF
    1146             :       END IF
    1147      748740 :       Kind: DO ikind = 1, SIZE(molecule_kind_set)
    1148      729684 :          nmol_local = local_molecules%n_el(ikind)
    1149     1617340 :          Mol_local: DO imol_local = 1, nmol_local
    1150      868600 :             imol = local_molecules%list(ikind)%array(imol_local)
    1151      868600 :             molecule => molecule_set(imol)
    1152      868600 :             CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
    1153     3543700 :             Particle: DO ipart = first_atom, last_atom
    1154     1945416 :                ii = ii + 1
    1155     1945416 :                IF (present_vel) THEN
    1156      972708 :                   vel(1, ipart) = vel(1, ipart)*map_info%p_scale(1, ii)%point
    1157      972708 :                   vel(2, ipart) = vel(2, ipart)*map_info%p_scale(2, ii)%point
    1158      972708 :                   vel(3, ipart) = vel(3, ipart)*map_info%p_scale(3, ii)%point
    1159             :                ELSE
    1160      972708 :                   particle_set(ipart)%v(1) = particle_set(ipart)%v(1)*map_info%p_scale(1, ii)%point
    1161      972708 :                   particle_set(ipart)%v(2) = particle_set(ipart)%v(2)*map_info%p_scale(2, ii)%point
    1162      972708 :                   particle_set(ipart)%v(3) = particle_set(ipart)%v(3)*map_info%p_scale(3, ii)%point
    1163             :                END IF
    1164             :                ! If Shell Adiabatic then apply the NHC thermostat also to the Shells
    1165     2814016 :                IF (shell_adiabatic) THEN
    1166      152160 :                   shell_index = particle_set(ipart)%shell_index
    1167      152160 :                   IF (shell_index /= 0) THEN
    1168      150880 :                      jj = jj + 2
    1169      150880 :                      atomic_kind => particle_set(ipart)%atomic_kind
    1170      150880 :                      CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell=shell)
    1171      150880 :                      fac_masss = shell%mass_shell/mass
    1172      150880 :                      fac_massc = shell%mass_core/mass
    1173      150880 :                      IF (present_vel) THEN
    1174      301760 :                         vs(1:3) = shell_vel(1:3, shell_index)
    1175      301760 :                         vc(1:3) = core_vel(1:3, shell_index)
    1176       75440 :                         shell_vel(1, shell_index) = vel(1, ipart) + fac_massc*(vs(1) - vc(1))
    1177       75440 :                         shell_vel(2, shell_index) = vel(2, ipart) + fac_massc*(vs(2) - vc(2))
    1178       75440 :                         shell_vel(3, shell_index) = vel(3, ipart) + fac_massc*(vs(3) - vc(3))
    1179       75440 :                         core_vel(1, shell_index) = vel(1, ipart) + fac_masss*(vc(1) - vs(1))
    1180       75440 :                         core_vel(2, shell_index) = vel(2, ipart) + fac_masss*(vc(2) - vs(2))
    1181       75440 :                         core_vel(3, shell_index) = vel(3, ipart) + fac_masss*(vc(3) - vs(3))
    1182             :                      ELSE
    1183      301760 :                         vs(1:3) = shell_particle_set(shell_index)%v(1:3)
    1184      301760 :                         vc(1:3) = core_particle_set(shell_index)%v(1:3)
    1185       75440 :                         shell_particle_set(shell_index)%v(1) = particle_set(ipart)%v(1) + fac_massc*(vs(1) - vc(1))
    1186       75440 :                         shell_particle_set(shell_index)%v(2) = particle_set(ipart)%v(2) + fac_massc*(vs(2) - vc(2))
    1187       75440 :                         shell_particle_set(shell_index)%v(3) = particle_set(ipart)%v(3) + fac_massc*(vs(3) - vc(3))
    1188       75440 :                         core_particle_set(shell_index)%v(1) = particle_set(ipart)%v(1) + fac_masss*(vc(1) - vs(1))
    1189       75440 :                         core_particle_set(shell_index)%v(2) = particle_set(ipart)%v(2) + fac_masss*(vc(2) - vs(2))
    1190       75440 :                         core_particle_set(shell_index)%v(3) = particle_set(ipart)%v(3) + fac_masss*(vc(3) - vs(3))
    1191             :                      END IF
    1192             :                   END IF
    1193             :                END IF
    1194             :             END DO Particle
    1195             :          END DO Mol_local
    1196             :       END DO Kind
    1197             : 
    1198       19056 :    END SUBROUTINE vel_rescale_particles
    1199             : 
    1200             : ! **************************************************************************************************
    1201             : !> \brief ...
    1202             : !> \param map_info ...
    1203             : !> \param particle_set ...
    1204             : !> \param atomic_kind_set ...
    1205             : !> \param local_particles ...
    1206             : !> \param group ...
    1207             : !> \param core_particle_set ...
    1208             : !> \param shell_particle_set ...
    1209             : !> \param core_vel ...
    1210             : !> \param shell_vel ...
    1211             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
    1212             : ! **************************************************************************************************
    1213        1840 :    SUBROUTINE ke_region_shells(map_info, particle_set, atomic_kind_set, &
    1214             :                                local_particles, group, core_particle_set, shell_particle_set, &
    1215        1840 :                                core_vel, shell_vel)
    1216             : 
    1217             :       TYPE(map_info_type), POINTER                       :: map_info
    1218             :       TYPE(particle_type), POINTER                       :: particle_set(:)
    1219             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
    1220             :       TYPE(distribution_1d_type), POINTER                :: local_particles
    1221             :       TYPE(mp_comm_type), INTENT(IN)                     :: group
    1222             :       TYPE(particle_type), OPTIONAL, POINTER             :: core_particle_set(:), &
    1223             :                                                             shell_particle_set(:)
    1224             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: core_vel(:, :), shell_vel(:, :)
    1225             : 
    1226             :       INTEGER                                            :: ii, iparticle, iparticle_kind, &
    1227             :                                                             iparticle_local, nparticle_kind, &
    1228             :                                                             nparticle_local, shell_index
    1229             :       LOGICAL                                            :: is_shell, present_vel
    1230             :       REAL(dp)                                           :: mass, mu_mass, v_sc(3)
    1231             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1232             :       TYPE(shell_kind_type), POINTER                     :: shell
    1233             : 
    1234        1840 :       present_vel = PRESENT(shell_vel)
    1235             :       ! Preliminary checks for consistency usage
    1236        1840 :       IF (present_vel) THEN
    1237         920 :          CPASSERT(PRESENT(core_vel))
    1238             :       ELSE
    1239         920 :          CPASSERT(PRESENT(shell_particle_set))
    1240         920 :          CPASSERT(PRESENT(core_particle_set))
    1241             :       END IF
    1242             :       ! get force on first thermostat for all the chains in the system.
    1243      154040 :       map_info%v_scale = 1.0_dp
    1244      154040 :       map_info%s_kin = 0.0_dp
    1245        1840 :       ii = 0
    1246             : 
    1247        1840 :       nparticle_kind = SIZE(atomic_kind_set)
    1248        5520 :       DO iparticle_kind = 1, nparticle_kind
    1249        3680 :          atomic_kind => atomic_kind_set(iparticle_kind)
    1250        3680 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell_active=is_shell, shell=shell)
    1251        5520 :          IF (is_shell) THEN
    1252        3680 :             mu_mass = shell%mass_shell*shell%mass_core/mass
    1253        3680 :             nparticle_local = local_particles%n_el(iparticle_kind)
    1254       92000 :             DO iparticle_local = 1, nparticle_local
    1255       88320 :                iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    1256       88320 :                shell_index = particle_set(iparticle)%shell_index
    1257       88320 :                ii = ii + 1
    1258       92000 :                IF (present_vel) THEN
    1259       44160 :                   v_sc(1) = core_vel(1, shell_index) - shell_vel(1, shell_index)
    1260       44160 :                   v_sc(2) = core_vel(2, shell_index) - shell_vel(2, shell_index)
    1261       44160 :                   v_sc(3) = core_vel(3, shell_index) - shell_vel(3, shell_index)
    1262       44160 :                   map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mu_mass*v_sc(1)**2
    1263       44160 :                   map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mu_mass*v_sc(2)**2
    1264       44160 :                   map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mu_mass*v_sc(3)**2
    1265             :                ELSE
    1266       44160 :                   v_sc(1) = core_particle_set(shell_index)%v(1) - shell_particle_set(shell_index)%v(1)
    1267       44160 :                   v_sc(2) = core_particle_set(shell_index)%v(2) - shell_particle_set(shell_index)%v(2)
    1268       44160 :                   v_sc(3) = core_particle_set(shell_index)%v(3) - shell_particle_set(shell_index)%v(3)
    1269       44160 :                   map_info%p_kin(1, ii)%point = map_info%p_kin(1, ii)%point + mu_mass*v_sc(1)**2
    1270       44160 :                   map_info%p_kin(2, ii)%point = map_info%p_kin(2, ii)%point + mu_mass*v_sc(2)**2
    1271       44160 :                   map_info%p_kin(3, ii)%point = map_info%p_kin(3, ii)%point + mu_mass*v_sc(3)**2
    1272             :                END IF
    1273             :             END DO
    1274             :          END IF
    1275             :       END DO
    1276        2880 :       IF (map_info%dis_type == do_thermo_communication) CALL group%sum(map_info%s_kin)
    1277             : 
    1278        1840 :    END SUBROUTINE ke_region_shells
    1279             : 
    1280             : ! **************************************************************************************************
    1281             : !> \brief ...
    1282             : !> \param map_info ...
    1283             : !> \param atomic_kind_set ...
    1284             : !> \param particle_set ...
    1285             : !> \param local_particles ...
    1286             : !> \param shell_particle_set ...
    1287             : !> \param core_particle_set ...
    1288             : !> \param shell_vel ...
    1289             : !> \param core_vel ...
    1290             : !> \param vel ...
    1291             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
    1292             : ! **************************************************************************************************
    1293        1600 :    SUBROUTINE vel_rescale_shells(map_info, atomic_kind_set, particle_set, local_particles, &
    1294        1600 :                                  shell_particle_set, core_particle_set, shell_vel, core_vel, vel)
    1295             : 
    1296             :       TYPE(map_info_type), POINTER                       :: map_info
    1297             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
    1298             :       TYPE(particle_type), POINTER                       :: particle_set(:)
    1299             :       TYPE(distribution_1d_type), POINTER                :: local_particles
    1300             :       TYPE(particle_type), OPTIONAL, POINTER             :: shell_particle_set(:), &
    1301             :                                                             core_particle_set(:)
    1302             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: shell_vel(:, :), core_vel(:, :), &
    1303             :                                                             vel(:, :)
    1304             : 
    1305             :       INTEGER                                            :: ii, iparticle, iparticle_kind, &
    1306             :                                                             iparticle_local, nparticle_kind, &
    1307             :                                                             nparticle_local, shell_index
    1308             :       LOGICAL                                            :: is_shell, present_vel
    1309             :       REAL(dp)                                           :: mass, massc, masss, umass, v(3), vc(3), &
    1310             :                                                             vs(3)
    1311             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1312             :       TYPE(shell_kind_type), POINTER                     :: shell
    1313             : 
    1314        1600 :       present_vel = PRESENT(vel)
    1315             :       ! Preliminary checks for consistency usage
    1316        1600 :       IF (present_vel) THEN
    1317         800 :          CPASSERT(PRESENT(shell_vel))
    1318         800 :          CPASSERT(PRESENT(core_vel))
    1319             :       ELSE
    1320         800 :          CPASSERT(PRESENT(shell_particle_set))
    1321         800 :          CPASSERT(PRESENT(core_particle_set))
    1322             :       END IF
    1323        1600 :       ii = 0
    1324        1600 :       nparticle_kind = SIZE(atomic_kind_set)
    1325             :       ! now scale the core-shell velocities
    1326        4800 :       Kind: DO iparticle_kind = 1, nparticle_kind
    1327        3200 :          atomic_kind => atomic_kind_set(iparticle_kind)
    1328        3200 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell_active=is_shell, shell=shell)
    1329        4800 :          IF (is_shell) THEN
    1330        3200 :             umass = 1.0_dp/mass
    1331        3200 :             masss = shell%mass_shell*umass
    1332        3200 :             massc = shell%mass_core*umass
    1333             : 
    1334        3200 :             nparticle_local = local_particles%n_el(iparticle_kind)
    1335       80000 :             Particles: DO iparticle_local = 1, nparticle_local
    1336       76800 :                iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    1337       76800 :                shell_index = particle_set(iparticle)%shell_index
    1338       76800 :                ii = ii + 1
    1339       80000 :                IF (present_vel) THEN
    1340      153600 :                   vc(1:3) = core_vel(1:3, shell_index)
    1341      153600 :                   vs(1:3) = shell_vel(1:3, shell_index)
    1342      153600 :                   v(1:3) = vel(1:3, iparticle)
    1343       38400 :                   shell_vel(1, shell_index) = v(1) + map_info%p_scale(1, ii)%point*massc*(vs(1) - vc(1))
    1344       38400 :                   shell_vel(2, shell_index) = v(2) + map_info%p_scale(2, ii)%point*massc*(vs(2) - vc(2))
    1345       38400 :                   shell_vel(3, shell_index) = v(3) + map_info%p_scale(3, ii)%point*massc*(vs(3) - vc(3))
    1346       38400 :                   core_vel(1, shell_index) = v(1) + map_info%p_scale(1, ii)%point*masss*(vc(1) - vs(1))
    1347       38400 :                   core_vel(2, shell_index) = v(2) + map_info%p_scale(2, ii)%point*masss*(vc(2) - vs(2))
    1348       38400 :                   core_vel(3, shell_index) = v(3) + map_info%p_scale(3, ii)%point*masss*(vc(3) - vs(3))
    1349             :                ELSE
    1350      153600 :                   vc(1:3) = core_particle_set(shell_index)%v(1:3)
    1351      153600 :                   vs(1:3) = shell_particle_set(shell_index)%v(1:3)
    1352      153600 :                   v(1:3) = particle_set(iparticle)%v(1:3)
    1353       38400 :                   shell_particle_set(shell_index)%v(1) = v(1) + map_info%p_scale(1, ii)%point*massc*(vs(1) - vc(1))
    1354       38400 :                   shell_particle_set(shell_index)%v(2) = v(2) + map_info%p_scale(2, ii)%point*massc*(vs(2) - vc(2))
    1355       38400 :                   shell_particle_set(shell_index)%v(3) = v(3) + map_info%p_scale(3, ii)%point*massc*(vs(3) - vc(3))
    1356       38400 :                   core_particle_set(shell_index)%v(1) = v(1) + map_info%p_scale(1, ii)%point*masss*(vc(1) - vs(1))
    1357       38400 :                   core_particle_set(shell_index)%v(2) = v(2) + map_info%p_scale(2, ii)%point*masss*(vc(2) - vs(2))
    1358       38400 :                   core_particle_set(shell_index)%v(3) = v(3) + map_info%p_scale(3, ii)%point*masss*(vc(3) - vs(3))
    1359             :                END IF
    1360             :             END DO Particles
    1361             :          END IF
    1362             :       END DO Kind
    1363             : 
    1364        1600 :    END SUBROUTINE vel_rescale_shells
    1365             : 
    1366             : ! **************************************************************************************************
    1367             : !> \brief Calculates kinetic energy and potential energy of the nhc variables
    1368             : !> \param nhc ...
    1369             : !> \param nhc_pot ...
    1370             : !> \param nhc_kin ...
    1371             : !> \param para_env ...
    1372             : !> \param array_kin ...
    1373             : !> \param array_pot ...
    1374             : !> \par History
    1375             : !>      none
    1376             : !> \author CJM
    1377             : ! **************************************************************************************************
    1378       10650 :    SUBROUTINE get_nhc_energies(nhc, nhc_pot, nhc_kin, para_env, array_kin, array_pot)
    1379             :       TYPE(lnhc_parameters_type), POINTER                :: nhc
    1380             :       REAL(KIND=dp), INTENT(OUT)                         :: nhc_pot, nhc_kin
    1381             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1382             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: array_kin, array_pot
    1383             : 
    1384             :       INTEGER                                            :: imap, l, n, number
    1385             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: akin, vpot
    1386             : 
    1387       10650 :       number = nhc%glob_num_nhc
    1388       31950 :       ALLOCATE (akin(number))
    1389       31950 :       ALLOCATE (vpot(number))
    1390      774086 :       akin = 0.0_dp
    1391      774086 :       vpot = 0.0_dp
    1392      395893 :       DO n = 1, nhc%loc_num_nhc
    1393      385243 :          imap = nhc%map_info%index(n)
    1394     1796788 :          DO l = 1, nhc%nhc_len
    1395     1400895 :             akin(imap) = akin(imap) + 0.5_dp*nhc%nvt(l, n)%mass*nhc%nvt(l, n)%v**2
    1396     1786138 :             vpot(imap) = vpot(imap) + nhc%nvt(l, n)%nkt*nhc%nvt(l, n)%eta
    1397             :          END DO
    1398             :       END DO
    1399             : 
    1400             :       ! Handle the thermostat distribution
    1401       10650 :       IF (nhc%map_info%dis_type == do_thermo_no_communication) THEN
    1402        3726 :          CALL para_env%sum(akin)
    1403        3726 :          CALL para_env%sum(vpot)
    1404        6924 :       ELSE IF (nhc%map_info%dis_type == do_thermo_communication) THEN
    1405        4940 :          CALL communication_thermo_low1(akin, number, para_env)
    1406        4940 :          CALL communication_thermo_low1(vpot, number, para_env)
    1407             :       END IF
    1408      774086 :       nhc_kin = SUM(akin)
    1409      774086 :       nhc_pot = SUM(vpot)
    1410             : 
    1411             :       ! Possibly give back kinetic or potential energy arrays
    1412       10650 :       IF (PRESENT(array_pot)) THEN
    1413         274 :          IF (ASSOCIATED(array_pot)) THEN
    1414           0 :             CPASSERT(SIZE(array_pot) == number)
    1415             :          ELSE
    1416         822 :             ALLOCATE (array_pot(number))
    1417             :          END IF
    1418       35794 :          array_pot = vpot
    1419             :       END IF
    1420       10650 :       IF (PRESENT(array_kin)) THEN
    1421         274 :          IF (ASSOCIATED(array_kin)) THEN
    1422           0 :             CPASSERT(SIZE(array_kin) == number)
    1423             :          ELSE
    1424         822 :             ALLOCATE (array_kin(number))
    1425             :          END IF
    1426       35794 :          array_kin = akin
    1427             :       END IF
    1428       10650 :       DEALLOCATE (akin)
    1429       10650 :       DEALLOCATE (vpot)
    1430       10650 :    END SUBROUTINE get_nhc_energies
    1431             : 
    1432             : ! **************************************************************************************************
    1433             : !> \brief Calculates kinetic energy and potential energy
    1434             : !>      of the csvr  and gle thermostats
    1435             : !> \param map_info ...
    1436             : !> \param loc_num ...
    1437             : !> \param glob_num ...
    1438             : !> \param thermo_energy ...
    1439             : !> \param thermostat_kin ...
    1440             : !> \param para_env ...
    1441             : !> \param array_pot ...
    1442             : !> \param array_kin ...
    1443             : !> \par History generalized MI [07.2009]
    1444             : !> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
    1445             : ! **************************************************************************************************
    1446        4328 :    SUBROUTINE get_kin_energies(map_info, loc_num, glob_num, thermo_energy, thermostat_kin, &
    1447             :                                para_env, array_pot, array_kin)
    1448             : 
    1449             :       TYPE(map_info_type), POINTER                       :: map_info
    1450             :       INTEGER, INTENT(IN)                                :: loc_num, glob_num
    1451             :       REAL(dp), DIMENSION(:), INTENT(IN)                 :: thermo_energy
    1452             :       REAL(KIND=dp), INTENT(OUT)                         :: thermostat_kin
    1453             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1454             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: array_pot, array_kin
    1455             : 
    1456             :       INTEGER                                            :: imap, n, number
    1457             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: akin
    1458             : 
    1459        4328 :       number = glob_num
    1460       12984 :       ALLOCATE (akin(number))
    1461      278278 :       akin = 0.0_dp
    1462      142960 :       DO n = 1, loc_num
    1463      138632 :          imap = map_info%index(n)
    1464      142960 :          akin(imap) = thermo_energy(n)
    1465             :       END DO
    1466             : 
    1467             :       ! Handle the thermostat distribution
    1468        4328 :       IF (map_info%dis_type == do_thermo_no_communication) THEN
    1469        1306 :          CALL para_env%sum(akin)
    1470        3022 :       ELSE IF (map_info%dis_type == do_thermo_communication) THEN
    1471        2334 :          CALL communication_thermo_low1(akin, number, para_env)
    1472             :       END IF
    1473      278278 :       thermostat_kin = SUM(akin)
    1474             : 
    1475             :       ! Possibly give back kinetic or potential energy arrays
    1476        4328 :       IF (PRESENT(array_pot)) THEN
    1477          22 :          IF (ASSOCIATED(array_pot)) THEN
    1478           0 :             CPASSERT(SIZE(array_pot) == number)
    1479             :          ELSE
    1480          66 :             ALLOCATE (array_pot(number))
    1481             :          END IF
    1482          66 :          array_pot = 0.0_dp
    1483             :       END IF
    1484        4328 :       IF (PRESENT(array_kin)) THEN
    1485         440 :          IF (ASSOCIATED(array_kin)) THEN
    1486         418 :             CPASSERT(SIZE(array_kin) == number)
    1487             :          ELSE
    1488          66 :             ALLOCATE (array_kin(number))
    1489             :          END IF
    1490       22856 :          array_kin = akin
    1491             :       END IF
    1492        4328 :       DEALLOCATE (akin)
    1493        4328 :    END SUBROUTINE get_kin_energies
    1494             : 
    1495             : ! **************************************************************************************************
    1496             : !> \brief Calculates the temperatures of the regions when a thermostat is
    1497             : !>        applied
    1498             : !> \param map_info ...
    1499             : !> \param loc_num ...
    1500             : !> \param glob_num ...
    1501             : !> \param nkt ...
    1502             : !> \param dof ...
    1503             : !> \param para_env ...
    1504             : !> \param temp_tot ...
    1505             : !> \param array_temp ...
    1506             : !> \par History generalized MI [07.2009]
    1507             : !> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
    1508             : ! **************************************************************************************************
    1509         274 :    SUBROUTINE get_temperatures(map_info, loc_num, glob_num, nkt, dof, para_env, &
    1510             :                                temp_tot, array_temp)
    1511             :       TYPE(map_info_type), POINTER                       :: map_info
    1512             :       INTEGER, INTENT(IN)                                :: loc_num, glob_num
    1513             :       REAL(dp), DIMENSION(:), INTENT(IN)                 :: nkt, dof
    1514             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1515             :       REAL(KIND=dp), INTENT(OUT)                         :: temp_tot
    1516             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: array_temp
    1517             : 
    1518             :       INTEGER                                            :: i, imap, imap2, n, number
    1519             :       REAL(KIND=dp)                                      :: fdeg_of_free
    1520             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: akin, deg_of_free
    1521             : 
    1522         274 :       number = glob_num
    1523         822 :       ALLOCATE (akin(number))
    1524         822 :       ALLOCATE (deg_of_free(number))
    1525       35816 :       akin = 0.0_dp
    1526       35816 :       deg_of_free = 0.0_dp
    1527       18120 :       DO n = 1, loc_num
    1528       17846 :          imap = map_info%index(n)
    1529       17846 :          imap2 = map_info%map_index(n)
    1530       17846 :          IF (nkt(n) == 0.0_dp) CYCLE
    1531       17846 :          deg_of_free(imap) = REAL(dof(n), KIND=dp)
    1532       18120 :          akin(imap) = map_info%s_kin(imap2)
    1533             :       END DO
    1534             : 
    1535             :       ! Handle the thermostat distribution
    1536         274 :       IF (map_info%dis_type == do_thermo_no_communication) THEN
    1537         146 :          CALL para_env%sum(akin)
    1538         146 :          CALL para_env%sum(deg_of_free)
    1539         128 :       ELSE IF (map_info%dis_type == do_thermo_communication) THEN
    1540          22 :          CALL communication_thermo_low1(akin, number, para_env)
    1541          22 :          CALL communication_thermo_low1(deg_of_free, number, para_env)
    1542             :       END IF
    1543       35816 :       temp_tot = SUM(akin)
    1544       35816 :       fdeg_of_free = SUM(deg_of_free)
    1545             : 
    1546         274 :       temp_tot = temp_tot/fdeg_of_free
    1547         274 :       temp_tot = cp_unit_from_cp2k(temp_tot, "K_temp")
    1548             :       ! Possibly give back temperatures of the full set of regions
    1549         274 :       IF (PRESENT(array_temp)) THEN
    1550         274 :          IF (ASSOCIATED(array_temp)) THEN
    1551           0 :             CPASSERT(SIZE(array_temp) == number)
    1552             :          ELSE
    1553         822 :             ALLOCATE (array_temp(number))
    1554             :          END IF
    1555       35816 :          DO i = 1, number
    1556       35542 :             array_temp(i) = akin(i)/deg_of_free(i)
    1557       35816 :             array_temp(i) = cp_unit_from_cp2k(array_temp(i), "K_temp")
    1558             :          END DO
    1559             :       END IF
    1560         274 :       DEALLOCATE (akin)
    1561         274 :       DEALLOCATE (deg_of_free)
    1562         274 :    END SUBROUTINE get_temperatures
    1563             : 
    1564             : ! **************************************************************************************************
    1565             : !> \brief Calculates energy associated with a thermostat
    1566             : !> \param thermostat ...
    1567             : !> \param thermostat_pot ...
    1568             : !> \param thermostat_kin ...
    1569             : !> \param para_env ...
    1570             : !> \param array_pot ...
    1571             : !> \param array_kin ...
    1572             : !> \author Teodoro Laino [tlaino] - 10.2007 - University of Zurich
    1573             : ! **************************************************************************************************
    1574       57603 :    SUBROUTINE get_thermostat_energies(thermostat, thermostat_pot, thermostat_kin, para_env, &
    1575             :                                       array_pot, array_kin)
    1576             :       TYPE(thermostat_type), POINTER                     :: thermostat
    1577             :       REAL(KIND=dp), INTENT(OUT)                         :: thermostat_pot, thermostat_kin
    1578             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1579             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: array_pot, array_kin
    1580             : 
    1581             :       INTEGER                                            :: i
    1582       57603 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: thermo_energy
    1583             : 
    1584       57603 :       thermostat_pot = 0.0_dp
    1585       57603 :       thermostat_kin = 0.0_dp
    1586       57603 :       IF (ASSOCIATED(thermostat)) THEN
    1587       14248 :          IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
    1588             :             ! Energy associated with the Nose-Hoover thermostat
    1589       10322 :             CPASSERT(ASSOCIATED(thermostat%nhc))
    1590             :             CALL get_nhc_energies(thermostat%nhc, thermostat_pot, thermostat_kin, para_env, &
    1591       30418 :                                   array_pot, array_kin)
    1592        3926 :          ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
    1593             :             ! Energy associated with the CSVR thermostat
    1594        3502 :             CPASSERT(ASSOCIATED(thermostat%csvr))
    1595       10506 :             ALLOCATE (thermo_energy(thermostat%csvr%loc_num_csvr))
    1596       64700 :             DO i = 1, thermostat%csvr%loc_num_csvr
    1597       64700 :                thermo_energy(i) = thermostat%csvr%nvt(i)%thermostat_energy
    1598             :             END DO
    1599             :             CALL get_kin_energies(thermostat%csvr%map_info, thermostat%csvr%loc_num_csvr, &
    1600             :                                   thermostat%csvr%glob_num_csvr, thermo_energy, &
    1601       10462 :                                   thermostat_kin, para_env, array_pot, array_kin)
    1602        3502 :             DEALLOCATE (thermo_energy)
    1603             : 
    1604         424 :          ELSE IF (thermostat%type_of_thermostat == do_thermo_gle) THEN
    1605             :             ! Energy associated with the GLE thermostat
    1606         408 :             CPASSERT(ASSOCIATED(thermostat%gle))
    1607        1224 :             ALLOCATE (thermo_energy(thermostat%gle%loc_num_gle))
    1608       66504 :             DO i = 1, thermostat%gle%loc_num_gle
    1609       66504 :                thermo_energy(i) = thermostat%gle%nvt(i)%thermostat_energy
    1610             :             END DO
    1611             :             CALL get_kin_energies(thermostat%gle%map_info, thermostat%gle%loc_num_gle, &
    1612             :                                   thermostat%gle%glob_num_gle, thermo_energy, &
    1613        1224 :                                   thermostat_kin, para_env, array_pot, array_kin)
    1614         408 :             DEALLOCATE (thermo_energy)
    1615             : 
    1616             :             ![NB] nothing to do for Ad-Langevin?
    1617             : 
    1618             :          END IF
    1619             :       END IF
    1620             : 
    1621       57603 :    END SUBROUTINE get_thermostat_energies
    1622             : 
    1623             : ! **************************************************************************************************
    1624             : !> \brief Calculates the temperatures for each region associated to a thermostat
    1625             : !> \param thermostat ...
    1626             : !> \param tot_temperature ...
    1627             : !> \param para_env ...
    1628             : !> \param array_temp ...
    1629             : !> \author Teodoro Laino [tlaino] - 02.2008 - University of Zurich
    1630             : ! **************************************************************************************************
    1631         274 :    SUBROUTINE get_region_temperatures(thermostat, tot_temperature, para_env, array_temp)
    1632             :       TYPE(thermostat_type), POINTER                     :: thermostat
    1633             :       REAL(KIND=dp), INTENT(OUT)                         :: tot_temperature
    1634             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1635             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: array_temp
    1636             : 
    1637             :       INTEGER                                            :: i
    1638         274 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: dof, nkt
    1639             : 
    1640         274 :       IF (ASSOCIATED(thermostat)) THEN
    1641         274 :          IF (thermostat%type_of_thermostat == do_thermo_nose) THEN
    1642             :             ! Energy associated with the Nose-Hoover thermostat
    1643         252 :             CPASSERT(ASSOCIATED(thermostat%nhc))
    1644         756 :             ALLOCATE (nkt(thermostat%nhc%loc_num_nhc))
    1645         756 :             ALLOCATE (dof(thermostat%nhc%loc_num_nhc))
    1646       18054 :             DO i = 1, thermostat%nhc%loc_num_nhc
    1647       17802 :                nkt(i) = thermostat%nhc%nvt(1, i)%nkt
    1648       18054 :                dof(i) = REAL(thermostat%nhc%nvt(1, i)%degrees_of_freedom, KIND=dp)
    1649             :             END DO
    1650             :             CALL get_temperatures(thermostat%nhc%map_info, thermostat%nhc%loc_num_nhc, &
    1651         252 :                                   thermostat%nhc%glob_num_nhc, nkt, dof, para_env, tot_temperature, array_temp)
    1652         252 :             DEALLOCATE (nkt)
    1653         252 :             DEALLOCATE (dof)
    1654          22 :          ELSE IF (thermostat%type_of_thermostat == do_thermo_csvr) THEN
    1655             :             ! Energy associated with the CSVR thermostat
    1656          22 :             CPASSERT(ASSOCIATED(thermostat%csvr))
    1657             : 
    1658          66 :             ALLOCATE (nkt(thermostat%csvr%loc_num_csvr))
    1659          66 :             ALLOCATE (dof(thermostat%csvr%loc_num_csvr))
    1660          66 :             DO i = 1, thermostat%csvr%loc_num_csvr
    1661          44 :                nkt(i) = thermostat%csvr%nvt(i)%nkt
    1662          66 :                dof(i) = REAL(thermostat%csvr%nvt(i)%degrees_of_freedom, KIND=dp)
    1663             :             END DO
    1664             :             CALL get_temperatures(thermostat%csvr%map_info, thermostat%csvr%loc_num_csvr, &
    1665          22 :                                   thermostat%csvr%glob_num_csvr, nkt, dof, para_env, tot_temperature, array_temp)
    1666          22 :             DEALLOCATE (nkt)
    1667          22 :             DEALLOCATE (dof)
    1668           0 :          ELSE IF (thermostat%type_of_thermostat == do_thermo_al) THEN
    1669             :             ! Energy associated with the AD_LANGEVIN thermostat
    1670           0 :             CPASSERT(ASSOCIATED(thermostat%al))
    1671             : 
    1672           0 :             ALLOCATE (nkt(thermostat%al%loc_num_al))
    1673           0 :             ALLOCATE (dof(thermostat%al%loc_num_al))
    1674           0 :             DO i = 1, thermostat%al%loc_num_al
    1675           0 :                nkt(i) = thermostat%al%nvt(i)%nkt
    1676           0 :                dof(i) = REAL(thermostat%al%nvt(i)%degrees_of_freedom, KIND=dp)
    1677             :             END DO
    1678             :             CALL get_temperatures(thermostat%al%map_info, thermostat%al%loc_num_al, &
    1679           0 :                                   thermostat%al%glob_num_al, nkt, dof, para_env, tot_temperature, array_temp)
    1680           0 :             DEALLOCATE (nkt)
    1681           0 :             DEALLOCATE (dof)
    1682           0 :          ELSE IF (thermostat%type_of_thermostat == do_thermo_gle) THEN
    1683             :             ! Energy associated with the GLE thermostat
    1684           0 :             CPASSERT(ASSOCIATED(thermostat%gle))
    1685             : 
    1686           0 :             ALLOCATE (nkt(thermostat%gle%loc_num_gle))
    1687           0 :             ALLOCATE (dof(thermostat%gle%loc_num_gle))
    1688           0 :             DO i = 1, thermostat%gle%loc_num_gle
    1689           0 :                nkt(i) = thermostat%gle%nvt(i)%nkt
    1690           0 :                dof(i) = REAL(thermostat%gle%nvt(i)%degrees_of_freedom, KIND=dp)
    1691             :             END DO
    1692             :             CALL get_temperatures(thermostat%gle%map_info, thermostat%gle%loc_num_gle, &
    1693           0 :                                   thermostat%gle%glob_num_gle, nkt, dof, para_env, tot_temperature, array_temp)
    1694           0 :             DEALLOCATE (nkt)
    1695           0 :             DEALLOCATE (dof)
    1696             :          END IF
    1697             :       END IF
    1698             : 
    1699         274 :    END SUBROUTINE get_region_temperatures
    1700             : 
    1701             : ! **************************************************************************************************
    1702             : !> \brief Prints status of all thermostats during an MD run
    1703             : !> \param thermostats ...
    1704             : !> \param para_env ...
    1705             : !> \param my_pos ...
    1706             : !> \param my_act ...
    1707             : !> \param itimes ...
    1708             : !> \param time ...
    1709             : !> \author Teodoro Laino [tlaino] - 02.2008 - University of Zurich
    1710             : ! **************************************************************************************************
    1711       44645 :    SUBROUTINE print_thermostats_status(thermostats, para_env, my_pos, my_act, itimes, time)
    1712             :       TYPE(thermostats_type), POINTER                    :: thermostats
    1713             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1714             :       CHARACTER(LEN=default_string_length)               :: my_pos, my_act
    1715             :       INTEGER, INTENT(IN)                                :: itimes
    1716             :       REAL(KIND=dp), INTENT(IN)                          :: time
    1717             : 
    1718       44645 :       IF (ASSOCIATED(thermostats)) THEN
    1719       10264 :          IF (ASSOCIATED(thermostats%thermostat_part)) THEN
    1720        9930 :             CALL print_thermostat_status(thermostats%thermostat_part, para_env, my_pos, my_act, itimes, time)
    1721             :          END IF
    1722       10264 :          IF (ASSOCIATED(thermostats%thermostat_shell)) THEN
    1723         830 :             CALL print_thermostat_status(thermostats%thermostat_shell, para_env, my_pos, my_act, itimes, time)
    1724             :          END IF
    1725       10264 :          IF (ASSOCIATED(thermostats%thermostat_coef)) THEN
    1726           0 :             CALL print_thermostat_status(thermostats%thermostat_coef, para_env, my_pos, my_act, itimes, time)
    1727             :          END IF
    1728       10264 :          IF (ASSOCIATED(thermostats%thermostat_baro)) THEN
    1729        2302 :             CALL print_thermostat_status(thermostats%thermostat_baro, para_env, my_pos, my_act, itimes, time)
    1730             :          END IF
    1731             :       END IF
    1732       44645 :    END SUBROUTINE print_thermostats_status
    1733             : 
    1734             : ! **************************************************************************************************
    1735             : !> \brief Prints status of a specific thermostat
    1736             : !> \param thermostat ...
    1737             : !> \param para_env ...
    1738             : !> \param my_pos ...
    1739             : !> \param my_act ...
    1740             : !> \param itimes ...
    1741             : !> \param time ...
    1742             : !> \author Teodoro Laino [tlaino] - 02.2008 - University of Zurich
    1743             : ! **************************************************************************************************
    1744       13062 :    SUBROUTINE print_thermostat_status(thermostat, para_env, my_pos, my_act, itimes, time)
    1745             :       TYPE(thermostat_type), POINTER                     :: thermostat
    1746             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1747             :       CHARACTER(LEN=default_string_length)               :: my_pos, my_act
    1748             :       INTEGER, INTENT(IN)                                :: itimes
    1749             :       REAL(KIND=dp), INTENT(IN)                          :: time
    1750             : 
    1751             :       INTEGER                                            :: i, unit
    1752             :       LOGICAL                                            :: new_file
    1753             :       REAL(KIND=dp)                                      :: thermo_kin, thermo_pot, tot_temperature
    1754       13062 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: array_kin, array_pot, array_temp
    1755             :       TYPE(cp_logger_type), POINTER                      :: logger
    1756             :       TYPE(section_vals_type), POINTER                   :: print_key
    1757             : 
    1758       13062 :       NULLIFY (logger, print_key, array_pot, array_kin, array_temp)
    1759       26124 :       logger => cp_get_default_logger()
    1760             : 
    1761       13062 :       IF (ASSOCIATED(thermostat)) THEN
    1762             :          ! Print Energies
    1763       13062 :          print_key => section_vals_get_subs_vals(thermostat%section, "PRINT%ENERGY")
    1764       13062 :          IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    1765         296 :             CALL get_thermostat_energies(thermostat, thermo_pot, thermo_kin, para_env, array_pot, array_kin)
    1766             :             unit = cp_print_key_unit_nr(logger, thermostat%section, "PRINT%ENERGY", &
    1767             :                                         extension="."//TRIM(thermostat%label)//".tener", file_position=my_pos, &
    1768         296 :                                         file_action=my_act, is_new_file=new_file)
    1769         296 :             IF (unit > 0) THEN
    1770         148 :                IF (new_file) THEN
    1771          13 :                   WRITE (unit, '(A)') "# Thermostat Potential and Kinetic Energies - Total and per Region"
    1772          13 :                   WRITE (unit, '("#",3X,A,2X,A,13X,A,10X,A)') "Step Nr.", "Time[fs]", "Kin.[a.u.]", "Pot.[a.u.]"
    1773             :                END IF
    1774         148 :                WRITE (UNIT=unit, FMT="(I8, F12.3,6X,2F20.10)") itimes, time*femtoseconds, thermo_kin, thermo_pot
    1775         148 :                WRITE (unit, '(A,4F20.10)') "# KINETIC ENERGY REGIONS: ", array_kin(1:MIN(4, SIZE(array_kin)))
    1776        4499 :                DO i = 5, SIZE(array_kin), 4
    1777        4499 :                   WRITE (UNIT=unit, FMT='("#",25X,4F20.10)') array_kin(i:MIN(i + 3, SIZE(array_kin)))
    1778             :                END DO
    1779         148 :                WRITE (unit, '(A,4F20.10)') "# POTENT. ENERGY REGIONS: ", array_pot(1:MIN(4, SIZE(array_pot)))
    1780        4499 :                DO i = 5, SIZE(array_pot), 4
    1781        4499 :                   WRITE (UNIT=unit, FMT='("#",25X,4F20.10)') array_pot(i:MIN(i + 3, SIZE(array_pot)))
    1782             :                END DO
    1783         148 :                CALL m_flush(unit)
    1784             :             END IF
    1785         296 :             DEALLOCATE (array_kin)
    1786         296 :             DEALLOCATE (array_pot)
    1787         296 :             CALL cp_print_key_finished_output(unit, logger, thermostat%section, "PRINT%ENERGY")
    1788             :          END IF
    1789             :          ! Print Temperatures of the regions
    1790       13062 :          print_key => section_vals_get_subs_vals(thermostat%section, "PRINT%TEMPERATURE")
    1791       13062 :          IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    1792         274 :             CALL get_region_temperatures(thermostat, tot_temperature, para_env, array_temp)
    1793             :             unit = cp_print_key_unit_nr(logger, thermostat%section, "PRINT%TEMPERATURE", &
    1794             :                                         extension="."//TRIM(thermostat%label)//".temp", file_position=my_pos, &
    1795         274 :                                         file_action=my_act, is_new_file=new_file)
    1796         274 :             IF (unit > 0) THEN
    1797         137 :                IF (new_file) THEN
    1798          12 :                   WRITE (unit, '(A)') "# Temperature Total and per Region"
    1799          12 :                   WRITE (unit, '("#",3X,A,2X,A,10X,A)') "Step Nr.", "Time[fs]", "Temp.[K]"
    1800             :                END IF
    1801         137 :                WRITE (UNIT=unit, FMT="(I8, F12.3,3X,F20.10)") itimes, time*femtoseconds, tot_temperature
    1802         137 :                WRITE (unit, '(A,I10)') "# TEMPERATURE REGIONS: ", SIZE(array_temp)
    1803        4625 :                DO i = 1, SIZE(array_temp), 4
    1804        4625 :                   WRITE (UNIT=unit, FMT='("#",22X,4F20.10)') array_temp(i:MIN(i + 3, SIZE(array_temp)))
    1805             :                END DO
    1806         137 :                CALL m_flush(unit)
    1807             :             END IF
    1808         274 :             DEALLOCATE (array_temp)
    1809         274 :             CALL cp_print_key_finished_output(unit, logger, thermostat%section, "PRINT%TEMPERATURE")
    1810             :          END IF
    1811             :       END IF
    1812       13062 :    END SUBROUTINE print_thermostat_status
    1813             : 
    1814             : ! **************************************************************************************************
    1815             : !> \brief Handles the communication for thermostats (1D array)
    1816             : !> \param array ...
    1817             : !> \param number ...
    1818             : !> \param para_env ...
    1819             : !> \author Teodoro Laino [tlaino] - University of Zurich 11.2007
    1820             : ! **************************************************************************************************
    1821       12258 :    SUBROUTINE communication_thermo_low1(array, number, para_env)
    1822             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: array
    1823             :       INTEGER, INTENT(IN)                                :: number
    1824             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1825             : 
    1826             :       INTEGER                                            :: i, icheck, ncheck
    1827       12258 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: work, work2
    1828             : 
    1829       36774 :       ALLOCATE (work(para_env%num_pe))
    1830       25548 :       DO i = 1, number
    1831       39870 :          work = 0.0_dp
    1832       13290 :          work(para_env%mepos + 1) = array(i)
    1833       66450 :          CALL para_env%sum(work)
    1834       39870 :          ncheck = COUNT(work /= 0.0_dp)
    1835       13290 :          array(i) = 0.0_dp
    1836       25548 :          IF (ncheck /= 0) THEN
    1837       37512 :             ALLOCATE (work2(ncheck))
    1838       12504 :             ncheck = 0
    1839       37512 :             DO icheck = 1, para_env%num_pe
    1840       37512 :                IF (work(icheck) /= 0.0_dp) THEN
    1841       24600 :                   ncheck = ncheck + 1
    1842       24600 :                   work2(ncheck) = work(icheck)
    1843             :                END IF
    1844             :             END DO
    1845       12504 :             CPASSERT(ncheck == SIZE(work2))
    1846       37104 :             CPASSERT(ALL(work2 == work2(1)))
    1847             : 
    1848       12504 :             array(i) = work2(1)
    1849       12504 :             DEALLOCATE (work2)
    1850             :          END IF
    1851             :       END DO
    1852       12258 :       DEALLOCATE (work)
    1853       12258 :    END SUBROUTINE communication_thermo_low1
    1854             : 
    1855             : ! **************************************************************************************************
    1856             : !> \brief Handles the communication for thermostats (2D array)
    1857             : !> \param array ...
    1858             : !> \param number1 ...
    1859             : !> \param number2 ...
    1860             : !> \param para_env ...
    1861             : !> \author Teodoro Laino [tlaino] - University of Zurich 11.2007
    1862             : ! **************************************************************************************************
    1863         250 :    SUBROUTINE communication_thermo_low2(array, number1, number2, para_env)
    1864             :       INTEGER, DIMENSION(:, :), INTENT(INOUT)            :: array
    1865             :       INTEGER, INTENT(IN)                                :: number1, number2
    1866             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1867             : 
    1868             :       INTEGER                                            :: i, icheck, j, ncheck
    1869         250 :       INTEGER, DIMENSION(:, :), POINTER                  :: work, work2
    1870             : 
    1871        1000 :       ALLOCATE (work(number1, para_env%num_pe))
    1872         606 :       DO i = 1, number2
    1873      309364 :          work = 0
    1874      154504 :          work(:, para_env%mepos + 1) = array(:, i)
    1875      618372 :          CALL para_env%sum(work)
    1876         356 :          ncheck = 0
    1877        1068 :          DO j = 1, para_env%num_pe
    1878       23584 :             IF (ANY(work(:, j) /= 0)) THEN
    1879         660 :                ncheck = ncheck + 1
    1880             :             END IF
    1881             :          END DO
    1882      154504 :          array(:, i) = 0
    1883         606 :          IF (ncheck /= 0) THEN
    1884        1424 :             ALLOCATE (work2(number1, ncheck))
    1885         356 :             ncheck = 0
    1886        1068 :             DO icheck = 1, para_env%num_pe
    1887       23584 :                IF (ANY(work(:, icheck) /= 0)) THEN
    1888         660 :                   ncheck = ncheck + 1
    1889      572880 :                   work2(:, ncheck) = work(:, icheck)
    1890             :                END IF
    1891             :             END DO
    1892         356 :             CPASSERT(ncheck == SIZE(work2, 2))
    1893        1016 :             DO j = 1, ncheck
    1894      286796 :                CPASSERT(ALL(work2(:, j) == work2(:, 1)))
    1895             :             END DO
    1896      154504 :             array(:, i) = work2(:, 1)
    1897         356 :             DEALLOCATE (work2)
    1898             :          END IF
    1899             :       END DO
    1900         250 :       DEALLOCATE (work)
    1901         250 :    END SUBROUTINE communication_thermo_low2
    1902             : 
    1903             : END MODULE thermostat_utils

Generated by: LCOV version 1.15