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 Barostat utils
      10             : !> \author teo [tlaino] - University of Zurich - 02.2008
      11             : ! **************************************************************************************************
      12             : MODULE barostat_utils
      13             :    USE barostat_types,                  ONLY: barostat_type
      14             :    USE cell_types,                      ONLY: cell_type
      15             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      16             :                                               cp_logger_type
      17             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      18             :                                               cp_print_key_unit_nr
      19             :    USE extended_system_types,           ONLY: npt_info_type
      20             :    USE input_constants,                 ONLY: npe_f_ensemble,&
      21             :                                               npe_i_ensemble,&
      22             :                                               nph_uniaxial_damped_ensemble,&
      23             :                                               nph_uniaxial_ensemble,&
      24             :                                               npt_f_ensemble,&
      25             :                                               npt_i_ensemble,&
      26             :                                               npt_ia_ensemble
      27             :    USE kinds,                           ONLY: default_string_length,&
      28             :                                               dp
      29             :    USE machine,                         ONLY: m_flush
      30             :    USE physcon,                         ONLY: angstrom,&
      31             :                                               femtoseconds,&
      32             :                                               kelvin
      33             :    USE simpar_types,                    ONLY: simpar_type
      34             : #include "../../base/base_uses.f90"
      35             : 
      36             :    IMPLICIT NONE
      37             : 
      38             :    PRIVATE
      39             :    PUBLIC ::  get_baro_energies, print_barostat_status
      40             : 
      41             : ! *** Global parameters ***
      42             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'barostat_utils'
      43             : 
      44             : CONTAINS
      45             : ! **************************************************************************************************
      46             : !> \brief Calculates kinetic energy and potential of barostat
      47             : !> \param cell ...
      48             : !> \param simpar ...
      49             : !> \param npt ...
      50             : !> \param baro_kin ...
      51             : !> \param baro_pot ...
      52             : !> \par History
      53             : !>      none
      54             : !> \author CJM
      55             : ! **************************************************************************************************
      56        5352 :    SUBROUTINE get_baro_energies(cell, simpar, npt, baro_kin, baro_pot)
      57             : 
      58             :       TYPE(cell_type), POINTER                           :: cell
      59             :       TYPE(simpar_type), INTENT(IN)                      :: simpar
      60             :       TYPE(npt_info_type), DIMENSION(:, :), INTENT(IN)   :: npt
      61             :       REAL(KIND=dp), INTENT(OUT)                         :: baro_kin, baro_pot
      62             : 
      63             :       INTEGER                                            :: i, j
      64             :       REAL(dp)                                           :: iv0, v0, v_shock
      65             : 
      66             :       IF (simpar%ensemble == npt_i_ensemble .OR. simpar%ensemble == npe_i_ensemble &
      67        5352 :           .OR. simpar%ensemble == npt_ia_ensemble) THEN
      68        3332 :          baro_pot = simpar%p_ext*cell%deth
      69        3332 :          baro_kin = 0.5_dp*npt(1, 1)%v**2*npt(1, 1)%mass
      70             :       ELSE IF (simpar%ensemble == npt_f_ensemble .OR. simpar%ensemble == npe_f_ensemble) THEN
      71        1888 :          baro_pot = simpar%p_ext*cell%deth
      72        1888 :          baro_kin = 0.0_dp
      73        7552 :          DO i = 1, 3
      74       24544 :             DO j = 1, 3
      75       22656 :                baro_kin = baro_kin + 0.5_dp*npt(i, j)%v**2*npt(i, j)%mass
      76             :             END DO
      77             :          END DO
      78             :       ELSEIF (simpar%ensemble == nph_uniaxial_ensemble .OR. simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
      79         132 :          v0 = simpar%v0
      80         132 :          iv0 = 1._dp/v0
      81         132 :          v_shock = simpar%v_shock
      82             : 
      83             :          ! Valid only for orthorhombic cell
      84         132 :          baro_pot = -0.5_dp*v_shock*v_shock*(1._dp - cell%deth*iv0)**2 - simpar%p0*(v0 - cell%deth)
      85             :          ! Valid only for orthorhombic cell
      86         132 :          baro_kin = 0.5_dp*npt(1, 1)%v*npt(1, 1)%v*npt(1, 1)%mass
      87             :       END IF
      88             : 
      89        5352 :    END SUBROUTINE get_baro_energies
      90             : 
      91             : ! **************************************************************************************************
      92             : !> \brief Prints status of barostat during an MD run
      93             : !> \param barostat ...
      94             : !> \param simpar ...
      95             : !> \param my_pos ...
      96             : !> \param my_act ...
      97             : !> \param cell ...
      98             : !> \param itimes ...
      99             : !> \param time ...
     100             : !> \author Teodoro Laino [tlaino] - 02.2008 - University of Zurich
     101             : ! **************************************************************************************************
     102       41773 :    SUBROUTINE print_barostat_status(barostat, simpar, my_pos, my_act, cell, itimes, time)
     103             :       TYPE(barostat_type), POINTER                       :: barostat
     104             :       TYPE(simpar_type), INTENT(IN)                      :: simpar
     105             :       CHARACTER(LEN=default_string_length)               :: my_pos, my_act
     106             :       TYPE(cell_type), POINTER                           :: cell
     107             :       INTEGER, INTENT(IN)                                :: itimes
     108             :       REAL(KIND=dp), INTENT(IN)                          :: time
     109             : 
     110             :       INTEGER                                            :: baro, nfree
     111             :       LOGICAL                                            :: new_file
     112             :       REAL(KIND=dp)                                      :: baro_kin, baro_pot, temp
     113             :       TYPE(cp_logger_type), POINTER                      :: logger
     114             : 
     115       41773 :       NULLIFY (logger)
     116       41773 :       logger => cp_get_default_logger()
     117             : 
     118       41773 :       IF (ASSOCIATED(barostat)) THEN
     119             :          baro = cp_print_key_unit_nr(logger, barostat%section, "PRINT%ENERGY", &
     120        2660 :                                      extension=".bener", file_position=my_pos, file_action=my_act, is_new_file=new_file)
     121        2660 :          CALL get_baro_energies(cell, simpar, barostat%npt, baro_kin, baro_pot)
     122        2660 :          nfree = SIZE(barostat%npt, 1)*SIZE(barostat%npt, 2)
     123        2660 :          temp = 2.0_dp*baro_kin/REAL(nfree, dp)*kelvin
     124        2660 :          IF (baro > 0) THEN
     125         129 :             IF (new_file) THEN
     126          11 :                WRITE (baro, '("#",3X,A,10X,A,8X,3(5X,A,5X),3X,A)') "Step Nr.", "Time[fs]", "Kin.[a.u.]", &
     127          22 :                   "Temp.[K]", "Pot.[a.u.]", "Vol[Ang.^3]"
     128             :             END IF
     129         129 :             WRITE (UNIT=baro, FMT="(I10, F20.3,4F20.10)") itimes, time*femtoseconds, &
     130         258 :                baro_kin, temp, baro_pot, cell%deth*angstrom*angstrom*angstrom
     131         129 :             CALL m_flush(baro)
     132             :          END IF
     133        2660 :          CALL cp_print_key_finished_output(baro, logger, barostat%section, "PRINT%ENERGY")
     134             :       END IF
     135             : 
     136       41773 :    END SUBROUTINE print_barostat_status
     137             : 
     138             : END MODULE barostat_utils