LCOV - code coverage report
Current view: top level - src/motion - thermal_region_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 125 134 93.3 %
Date: 2024-04-20 06:29:22 Functions: 3 3 100.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 Setup of regions with different temperature
      10             : !> \par History
      11             : !>   - Added support for Langevin regions (2014/01/08, LT)
      12             : !>   - Added print subroutine for langevin regions (2014/02/04, LT)
      13             : !>   - Changed print_thermal_regions to print_thermal_regions_temperature
      14             : !>     (2014/02/04, LT)
      15             : !> \author MI
      16             : ! **************************************************************************************************
      17             : MODULE thermal_region_utils
      18             : 
      19             :    USE bibliography,                    ONLY: Kantorovich2008,&
      20             :                                               Kantorovich2008a,&
      21             :                                               cite_reference
      22             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      23             :                                               cp_logger_type,&
      24             :                                               cp_to_string
      25             :    USE cp_output_handling,              ONLY: cp_p_file,&
      26             :                                               cp_print_key_finished_output,&
      27             :                                               cp_print_key_should_output,&
      28             :                                               cp_print_key_unit_nr
      29             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      30             :                                               cp_subsys_type
      31             :    USE force_env_types,                 ONLY: force_env_get,&
      32             :                                               force_env_type
      33             :    USE input_constants,                 ONLY: langevin_ensemble,&
      34             :                                               npt_f_ensemble,&
      35             :                                               npt_i_ensemble,&
      36             :                                               npt_ia_ensemble,&
      37             :                                               nvt_ensemble
      38             :    USE input_section_types,             ONLY: section_vals_get,&
      39             :                                               section_vals_get_subs_vals,&
      40             :                                               section_vals_type,&
      41             :                                               section_vals_val_get
      42             :    USE kinds,                           ONLY: default_string_length,&
      43             :                                               dp
      44             :    USE memory_utilities,                ONLY: reallocate
      45             :    USE particle_list_types,             ONLY: particle_list_type
      46             :    USE physcon,                         ONLY: femtoseconds,&
      47             :                                               kelvin
      48             :    USE simpar_types,                    ONLY: simpar_type
      49             :    USE string_utilities,                ONLY: integer_to_string
      50             :    USE thermal_region_types,            ONLY: allocate_thermal_regions,&
      51             :                                               release_thermal_regions,&
      52             :                                               thermal_region_type,&
      53             :                                               thermal_regions_type
      54             : #include "../base/base_uses.f90"
      55             : 
      56             :    IMPLICIT NONE
      57             : 
      58             :    PRIVATE
      59             :    PUBLIC :: create_thermal_regions, &
      60             :              print_thermal_regions_temperature, &
      61             :              print_thermal_regions_langevin
      62             : 
      63             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'thermal_region_utils'
      64             : 
      65             : CONTAINS
      66             : 
      67             : ! **************************************************************************************************
      68             : !> \brief create thermal_regions
      69             : !> \param thermal_regions ...
      70             : !> \param md_section ...
      71             : !> \param simpar ...
      72             : !> \param force_env ...
      73             : !> \par History
      74             : !>   - Added support for Langevin regions (2014/01/08, LT)
      75             : !> \author
      76             : ! **************************************************************************************************
      77        3592 :    SUBROUTINE create_thermal_regions(thermal_regions, md_section, simpar, force_env)
      78             :       TYPE(thermal_regions_type), POINTER                :: thermal_regions
      79             :       TYPE(section_vals_type), POINTER                   :: md_section
      80             :       TYPE(simpar_type), POINTER                         :: simpar
      81             :       TYPE(force_env_type), POINTER                      :: force_env
      82             : 
      83             :       CHARACTER(LEN=default_string_length)               :: my_region
      84             :       INTEGER                                            :: i, il, ipart, ireg, nlist, nregions
      85        1796 :       INTEGER, DIMENSION(:), POINTER                     :: tmplist
      86             :       LOGICAL                                            :: apply_thermostat, do_langevin, &
      87             :                                                             do_langevin_default, do_read_ngr, &
      88             :                                                             explicit
      89             :       REAL(KIND=dp)                                      :: temp, temp_tol
      90             :       TYPE(cp_subsys_type), POINTER                      :: subsys
      91             :       TYPE(particle_list_type), POINTER                  :: particles
      92             :       TYPE(section_vals_type), POINTER                   :: region_sections, thermal_region_section
      93             :       TYPE(thermal_region_type), POINTER                 :: t_region
      94             : 
      95        1796 :       NULLIFY (region_sections, t_region, thermal_region_section, particles, subsys, tmplist)
      96           0 :       ALLOCATE (thermal_regions)
      97        1796 :       CALL allocate_thermal_regions(thermal_regions)
      98        1796 :       thermal_region_section => section_vals_get_subs_vals(md_section, "THERMAL_REGION")
      99        1796 :       CALL section_vals_get(thermal_region_section, explicit=explicit)
     100        1796 :       IF (explicit) THEN
     101             :          apply_thermostat = (simpar%ensemble == nvt_ensemble) .OR. &
     102             :                             (simpar%ensemble == npt_f_ensemble) .OR. &
     103             :                             (simpar%ensemble == npt_ia_ensemble) .OR. &
     104          14 :                             (simpar%ensemble == npt_i_ensemble)
     105             :          IF (apply_thermostat) THEN
     106             :             CALL cp_warn(__LOCATION__, &
     107             :                          "With the chosen ensemble the temperature is "// &
     108             :                          "controlled by thermostats. The definition of different thermal "// &
     109           0 :                          "regions might result inconsistent with the presence of thermostats.")
     110             :          END IF
     111          14 :          IF (simpar%temp_tol > 0.0_dp) THEN
     112             :             CALL cp_warn(__LOCATION__, &
     113             :                          "Control of the global temperature by rescaling of the velocity "// &
     114             :                          "is not consistent with the presence of different thermal regions. "// &
     115           0 :                          "The temperature of different regions is rescaled separatedly.")
     116             :          END IF
     117             :          CALL section_vals_val_get(thermal_region_section, "FORCE_RESCALING", &
     118          14 :                                    l_val=thermal_regions%force_rescaling)
     119             :          region_sections => section_vals_get_subs_vals(thermal_region_section, &
     120          14 :                                                        "DEFINE_REGION")
     121          14 :          CALL section_vals_get(region_sections, n_repetition=nregions)
     122          14 :          IF (nregions > 0) THEN
     123          14 :             thermal_regions%nregions = nregions
     124          14 :             thermal_regions%section => thermal_region_section
     125          62 :             ALLOCATE (thermal_regions%thermal_region(nregions))
     126          14 :             CALL force_env_get(force_env, subsys=subsys)
     127          14 :             CALL cp_subsys_get(subsys, particles=particles)
     128          14 :             IF (simpar%ensemble == langevin_ensemble) THEN
     129          12 :                CALL cite_reference(Kantorovich2008)
     130          12 :                CALL cite_reference(Kantorovich2008a)
     131             :                CALL section_vals_val_get(thermal_region_section, "DO_LANGEVIN_DEFAULT", &
     132          12 :                                          l_val=do_langevin_default)
     133          36 :                ALLOCATE (thermal_regions%do_langevin(particles%n_els))
     134          96 :                thermal_regions%do_langevin = do_langevin_default
     135             :             END IF
     136          34 :             DO ireg = 1, nregions
     137          20 :                NULLIFY (t_region)
     138          20 :                t_region => thermal_regions%thermal_region(ireg)
     139          20 :                t_region%region_index = ireg
     140             :                CALL section_vals_val_get(region_sections, "LIST", &
     141          20 :                                          i_rep_section=ireg, n_rep_val=nlist)
     142          20 :                NULLIFY (t_region%part_index)
     143          20 :                t_region%npart = 0
     144          20 :                IF (simpar%ensemble == langevin_ensemble) THEN
     145             :                   CALL section_vals_val_get(region_sections, "DO_LANGEVIN", &
     146          16 :                                             i_rep_section=ireg, l_val=do_langevin)
     147             :                END IF
     148          42 :                DO il = 1, nlist
     149             :                   CALL section_vals_val_get(region_sections, "LIST", i_rep_section=ireg, &
     150          22 :                                             i_rep_val=il, i_vals=tmplist)
     151          22 :                   CALL reallocate(t_region%part_index, 1, t_region%npart + SIZE(tmplist))
     152         162 :                   DO i = 1, SIZE(tmplist)
     153         120 :                      ipart = tmplist(i)
     154         120 :                      CPASSERT(((ipart > 0) .AND. (ipart <= particles%n_els)))
     155         120 :                      t_region%npart = t_region%npart + 1
     156         120 :                      t_region%part_index(t_region%npart) = ipart
     157         120 :                      particles%els(ipart)%t_region_index = ireg
     158         142 :                      IF (simpar%ensemble == langevin_ensemble) THEN
     159          44 :                         thermal_regions%do_langevin(ipart) = do_langevin
     160             :                      END IF
     161             :                   END DO
     162             :                END DO
     163             :                CALL section_vals_val_get(region_sections, "TEMPERATURE", i_rep_section=ireg, &
     164          20 :                                          r_val=temp)
     165          20 :                t_region%temp_expected = temp
     166             :                CALL section_vals_val_get(region_sections, "TEMP_TOL", i_rep_section=ireg, &
     167          20 :                                          r_val=temp_tol)
     168          20 :                t_region%temp_tol = temp_tol
     169          20 :                CALL section_vals_val_get(region_sections, "NOISY_GAMMA_REGION", i_rep_section=ireg, explicit=do_read_ngr)
     170          54 :                IF (do_read_ngr) THEN
     171             :                   CALL section_vals_val_get(region_sections, "NOISY_GAMMA_REGION", i_rep_section=ireg, &
     172           6 :                                             r_val=t_region%noisy_gamma_region)
     173           6 :                   IF (simpar%ensemble == langevin_ensemble) THEN
     174           6 :                      IF (.NOT. do_langevin) THEN
     175           0 :                         CALL integer_to_string(ireg, my_region)
     176             :                         CALL cp_warn(__LOCATION__, &
     177             :                                      "You provided NOISY_GAMMA_REGION but atoms in thermal region "//TRIM(my_region)// &
     178             :                                      " will not undergo Langevin MD. "// &
     179           0 :                                      "NOISY_GAMMA_REGION will be ignored and its value discarded!")
     180             :                      END IF
     181             :                   ELSE
     182             :                      CALL cp_warn(__LOCATION__, &
     183             :                                   "You provided NOISY_GAMMA_REGION but the Langevin Ensamble is not selected "// &
     184           0 :                                   "NOISY_GAMMA_REGION will be ignored and its value discarded!")
     185             :                   END IF
     186             :                ELSE
     187          14 :                   t_region%noisy_gamma_region = simpar%noisy_gamma
     188             :                END IF
     189             :             END DO
     190          14 :             simpar%do_thermal_region = .TRUE.
     191             :          ELSE
     192           0 :             CALL release_thermal_regions(thermal_regions)
     193           0 :             DEALLOCATE (thermal_regions)
     194           0 :             simpar%do_thermal_region = .FALSE.
     195             :          END IF
     196             :       ELSE
     197        1782 :          CALL release_thermal_regions(thermal_regions)
     198        1782 :          DEALLOCATE (thermal_regions)
     199        1782 :          simpar%do_thermal_region = .FALSE.
     200             :       END IF
     201             : 
     202        1796 :    END SUBROUTINE create_thermal_regions
     203             : 
     204             : ! **************************************************************************************************
     205             : !> \brief print_thermal_regions_temperature
     206             : !> \param thermal_regions : thermal regions type contains information
     207             : !>                          about the regions
     208             : !> \param itimes          : iteration number of the time step
     209             : !> \param time            : simulation time of the time step
     210             : !> \param pos             : file position
     211             : !> \param act             : file action
     212             : !> \par History
     213             : !>   - added doxygen header and changed subroutine name from
     214             : !>     print_thermal_regions to print_thermal_regions_temperature
     215             : !>     (2014/02/04, LT)
     216             : !> \author
     217             : ! **************************************************************************************************
     218       44645 :    SUBROUTINE print_thermal_regions_temperature(thermal_regions, itimes, time, pos, act)
     219             :       TYPE(thermal_regions_type), POINTER                :: thermal_regions
     220             :       INTEGER, INTENT(IN)                                :: itimes
     221             :       REAL(KIND=dp), INTENT(IN)                          :: time
     222             :       CHARACTER(LEN=default_string_length)               :: pos, act
     223             : 
     224             :       CHARACTER(LEN=default_string_length)               :: fmd
     225             :       INTEGER                                            :: ireg, nregions, unit
     226             :       LOGICAL                                            :: new_file
     227       44645 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: temp
     228             :       TYPE(cp_logger_type), POINTER                      :: logger
     229             :       TYPE(section_vals_type), POINTER                   :: print_key
     230             : 
     231       44645 :       NULLIFY (logger)
     232       44645 :       logger => cp_get_default_logger()
     233             : 
     234       44645 :       IF (ASSOCIATED(thermal_regions)) THEN
     235         110 :          print_key => section_vals_get_subs_vals(thermal_regions%section, "PRINT%TEMPERATURE")
     236         110 :          IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     237             :             unit = cp_print_key_unit_nr(logger, thermal_regions%section, "PRINT%TEMPERATURE", &
     238             :                                         extension=".tregion", file_position=pos, &
     239          42 :                                         file_action=act, is_new_file=new_file)
     240          42 :             IF (unit > 0) THEN
     241          21 :                IF (new_file) THEN
     242           1 :                   WRITE (unit, '(A)') "# Temperature per Region"
     243           1 :                   WRITE (unit, '("#",3X,A,2X,A,13X,A)') "Step Nr.", "Time[fs]", "Temp.[K] ...."
     244             :                END IF
     245          21 :                nregions = thermal_regions%nregions
     246          63 :                ALLOCATE (temp(0:nregions))
     247          84 :                temp = 0.0_dp
     248          21 :                temp(0) = thermal_regions%temp_reg0
     249          63 :                DO ireg = 1, nregions
     250          63 :                   temp(ireg) = thermal_regions%thermal_region(ireg)%temperature
     251             :                END DO
     252          21 :                fmd = "(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(nregions + 1)))//"F20.6)"
     253             :                fmd = TRIM(fmd)
     254          21 :                WRITE (UNIT=unit, FMT=fmd) itimes, time, temp(0:nregions)
     255          21 :                DEALLOCATE (temp)
     256             :             END IF
     257          42 :             CALL cp_print_key_finished_output(unit, logger, thermal_regions%section, "PRINT%TEMPERATURE")
     258             :          END IF
     259             :       END IF
     260       44645 :    END SUBROUTINE print_thermal_regions_temperature
     261             : 
     262             : ! **************************************************************************************************
     263             : !> \brief print out information regarding to langevin regions defined in
     264             : !>        thermal_regions section
     265             : !> \param thermal_regions : thermal regions type containing the relevant
     266             : !>                          langevin regions information
     267             : !> \param simpar          : wrapper for simulation parameters
     268             : !> \param pos             : file position
     269             : !> \param act             : file action
     270             : !> \par History
     271             : !>   - created (2014/02/02, LT)
     272             : !> \author Lianheng Tong [LT] (tonglianheng@gmail.com)
     273             : ! **************************************************************************************************
     274          42 :    SUBROUTINE print_thermal_regions_langevin(thermal_regions, simpar, pos, act)
     275             :       TYPE(thermal_regions_type), POINTER                :: thermal_regions
     276             :       TYPE(simpar_type), POINTER                         :: simpar
     277             :       CHARACTER(LEN=default_string_length)               :: pos, act
     278             : 
     279             :       INTEGER                                            :: ipart, ipart_reg, ireg, natoms, &
     280             :                                                             print_unit
     281          42 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: region_id
     282             :       LOGICAL                                            :: new_file
     283          42 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: noisy_gamma_region, temperature
     284             :       TYPE(cp_logger_type), POINTER                      :: logger
     285             :       TYPE(section_vals_type), POINTER                   :: print_key
     286             : 
     287          42 :       NULLIFY (logger)
     288          42 :       logger => cp_get_default_logger()
     289             : 
     290          42 :       IF (ASSOCIATED(thermal_regions)) THEN
     291          12 :          IF (ASSOCIATED(thermal_regions%do_langevin)) THEN
     292             :             print_key => section_vals_get_subs_vals(thermal_regions%section, &
     293          12 :                                                     "PRINT%LANGEVIN_REGIONS")
     294          12 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
     295             :                       cp_p_file)) THEN
     296             :                print_unit = cp_print_key_unit_nr(logger, thermal_regions%section, &
     297             :                                                  "PRINT%LANGEVIN_REGIONS", &
     298             :                                                  extension=".lgv_regions", &
     299             :                                                  file_position=pos, file_action=act, &
     300          12 :                                                  is_new_file=new_file)
     301          12 :                IF (print_unit > 0) THEN
     302           6 :                   IF (new_file) THEN
     303           6 :                      WRITE (print_unit, '(A)') "# Atoms Undergoing Langevin MD"
     304             :                      WRITE (print_unit, '(A,3X,A,3X,A,3X,A,3X,A,3X,A)') &
     305           6 :                         "#", "Atom_ID", "Region_ID", "Langevin(L)/NVE(N)", "Expected_T[K]", "[NoisyGamma]"
     306             :                   END IF
     307           6 :                   natoms = SIZE(thermal_regions%do_langevin)
     308          18 :                   ALLOCATE (temperature(natoms))
     309          18 :                   ALLOCATE (region_id(natoms))
     310          18 :                   ALLOCATE (noisy_gamma_region(natoms))
     311          42 :                   temperature(:) = simpar%temp_ext
     312          42 :                   region_id(:) = 0
     313          42 :                   noisy_gamma_region(:) = simpar%noisy_gamma
     314          14 :                   DO ireg = 1, thermal_regions%nregions
     315          36 :                      DO ipart_reg = 1, thermal_regions%thermal_region(ireg)%npart
     316          22 :                         ipart = thermal_regions%thermal_region(ireg)%part_index(ipart_reg)
     317          22 :                         temperature(ipart) = thermal_regions%thermal_region(ireg)%temp_expected
     318          22 :                         region_id(ipart) = thermal_regions%thermal_region(ireg)%region_index
     319          30 :                         noisy_gamma_region(ipart) = thermal_regions%thermal_region(ireg)%noisy_gamma_region
     320             :                      END DO
     321             :                   END DO
     322          42 :                   DO ipart = 1, natoms
     323          36 :                      WRITE (print_unit, '(1X,I10,2X)', advance='no') ipart
     324          36 :                      WRITE (print_unit, '(I10,3X)', advance='no') region_id(ipart)
     325          42 :                      IF (thermal_regions%do_langevin(ipart)) THEN
     326          24 :                         WRITE (print_unit, '(A,3X)', advance='no') "L"
     327          24 :                         IF (noisy_gamma_region(ipart) > 0._dp) THEN
     328          10 :                            WRITE (print_unit, '(10X,F20.3,3X,ES9.3)') temperature(ipart)*kelvin, &
     329          20 :                               noisy_gamma_region(ipart)/femtoseconds
     330             :                         ELSE
     331          14 :                            WRITE (print_unit, '(10X,F20.3)') temperature(ipart)*kelvin
     332             :                         END IF
     333             :                      ELSE
     334          12 :                         WRITE (print_unit, '(A,3X)', advance='no') "N"
     335          12 :                         WRITE (print_unit, '(18X,A)') "--"
     336             :                      END IF
     337             :                   END DO
     338           6 :                   DEALLOCATE (region_id)
     339           6 :                   DEALLOCATE (temperature)
     340           6 :                   DEALLOCATE (noisy_gamma_region)
     341             :                END IF
     342             :                CALL cp_print_key_finished_output(print_unit, logger, thermal_regions%section, &
     343          12 :                                                  "PRINT%LANGEVIN_REGIONS")
     344             :             END IF
     345             :          END IF
     346             :       END IF
     347          42 :    END SUBROUTINE print_thermal_regions_langevin
     348             : 
     349             : END MODULE thermal_region_utils

Generated by: LCOV version 1.15