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

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief 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         5307 :    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         1769 :       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         1769 :       NULLIFY (region_sections, t_region, thermal_region_section, particles, subsys, tmplist)
      96            0 :       ALLOCATE (thermal_regions)
      97         1769 :       CALL allocate_thermal_regions(thermal_regions)
      98         1769 :       thermal_region_section => section_vals_get_subs_vals(md_section, "THERMAL_REGION")
      99         1769 :       CALL section_vals_get(thermal_region_section, explicit=explicit)
     100         1769 :       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         1755 :          CALL release_thermal_regions(thermal_regions)
     198         1755 :          DEALLOCATE (thermal_regions)
     199         1755 :          simpar%do_thermal_region = .FALSE.
     200              :       END IF
     201              : 
     202         1769 :    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        42149 :    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        42149 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: temp
     228              :       TYPE(cp_logger_type), POINTER                      :: logger
     229              :       TYPE(section_vals_type), POINTER                   :: print_key
     230              : 
     231        42149 :       NULLIFY (logger)
     232        42149 :       logger => cp_get_default_logger()
     233              : 
     234        42149 :       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        42149 :    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           12 :                   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 2.0-1