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 43107 : 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 43107 : NULLIFY (logger)
116 43107 : logger => cp_get_default_logger()
117 :
118 43107 : 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 43107 : END SUBROUTINE print_barostat_status
137 :
138 : END MODULE barostat_utils
|