Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Utilities for Molecular Dynamics
10 : !> \author Teodoro Laino [tlaino] - University of Zurich - 09.2007
11 : ! **************************************************************************************************
12 : MODULE md_util
13 :
14 : USE cp_files, ONLY: close_file,&
15 : open_file
16 : USE cp_log_handling, ONLY: cp_get_default_logger,&
17 : cp_logger_type,&
18 : cp_to_string
19 : USE cp_output_handling, ONLY: cp_print_key_generate_filename
20 : USE cp_units, ONLY: cp_unit_from_cp2k,&
21 : cp_unit_to_cp2k
22 : USE fparser, ONLY: EvalErrType,&
23 : evalf,&
24 : finalizef,&
25 : initf,&
26 : parsef
27 : USE input_cp2k_restarts, ONLY: write_restart
28 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
29 : section_vals_type,&
30 : section_vals_val_get
31 : USE kinds, ONLY: default_path_length,&
32 : default_string_length,&
33 : dp
34 : USE md_energies, ONLY: md_write_output
35 : USE md_environment_types, ONLY: get_md_env,&
36 : md_environment_type
37 : USE message_passing, ONLY: mp_para_env_type
38 : USE simpar_types, ONLY: simpar_type
39 : USE string_utilities, ONLY: integer_to_string
40 : USE thermal_region_types, ONLY: thermal_region_type,&
41 : thermal_regions_type
42 : #include "../base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 :
46 : PRIVATE
47 :
48 : ! *** Global parameters ***
49 :
50 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_util'
51 :
52 : PUBLIC :: md_output, &
53 : read_vib_eigs_unformatted, &
54 : update_expected_temperature
55 :
56 : CONTAINS
57 :
58 : ! **************************************************************************************************
59 : !> \brief update_expected_temperature
60 : !> \param md_env ...
61 : !> \date 03.2026
62 : !> \par Note on units: t_region%temperature is updated by
63 : !> SUBROUTINE scale_velocity_region in md_vel_utils.F
64 : !> which is in Kelvin, but t_region%temp_expected is
65 : !> in cp2k internal units. See "Note the unit conversion
66 : !> here" in SUBROUTINE print_thermal_regions_temperature
67 : !> in thermal_region_utils.F for output processing.
68 : !> \author HE Zilong
69 : ! **************************************************************************************************
70 42347 : SUBROUTINE update_expected_temperature(md_env)
71 : TYPE(md_environment_type), POINTER :: md_env
72 :
73 : CHARACTER(LEN=default_string_length) :: my_itimes, my_region
74 : INTEGER :: ireg
75 : INTEGER, POINTER :: itimes
76 : TYPE(simpar_type), POINTER :: simpar
77 : TYPE(thermal_region_type), POINTER :: t_region
78 : TYPE(thermal_regions_type), POINTER :: thermal_regions
79 :
80 42347 : NULLIFY (itimes, simpar)
81 42347 : CALL get_md_env(md_env, itimes=itimes, simpar=simpar)
82 :
83 : ! Update thermal regions
84 42347 : IF (simpar%do_thermal_region) THEN
85 96 : NULLIFY (thermal_regions)
86 96 : CALL get_md_env(md_env, thermal_regions=thermal_regions)
87 256 : DO ireg = 1, thermal_regions%nregions
88 160 : NULLIFY (t_region)
89 160 : t_region => thermal_regions%thermal_region(ireg)
90 256 : IF (LEN_TRIM(t_region%temperature_function) > 0) THEN
91 0 : CALL initf(1)
92 0 : CALL parsef(1, TRIM(t_region%temperature_function), t_region%temperature_function_parameters)
93 0 : t_region%temperature_function_values(1) = itimes*1.0_dp
94 0 : t_region%temperature_function_values(2) = t_region%temperature
95 0 : t_region%temp_expected = cp_unit_to_cp2k(evalf(1, t_region%temperature_function_values), "K")
96 0 : CPASSERT(EvalErrType <= 0)
97 0 : CALL finalizef()
98 0 : IF (t_region%temp_expected < 0.0_dp) THEN
99 0 : CALL integer_to_string(ireg, my_region)
100 0 : CALL integer_to_string(itimes, my_itimes)
101 : CALL cp_abort(__LOCATION__, &
102 : "At step "//TRIM(my_itimes)//" , the temperature "// &
103 : "evaluated for thermal region "//TRIM(my_region)//" is "// &
104 : TRIM(ADJUSTL(cp_to_string(cp_unit_from_cp2k(t_region%temp_expected, "K"))))// &
105 0 : " K, which is negative and unphysical!")
106 : END IF
107 : END IF
108 : END DO
109 : END IF
110 :
111 : ! Update thermostat regions: WIP
112 :
113 42347 : END SUBROUTINE update_expected_temperature
114 :
115 : ! **************************************************************************************************
116 : !> \brief collects the part of the MD that, basically, does the output
117 : !> \param md_env ...
118 : !> \param md_section ...
119 : !> \param root_section ...
120 : !> \param forced_io ...
121 : !> \par History
122 : !> 03.2006 created [Joost VandeVondele]
123 : ! **************************************************************************************************
124 42347 : SUBROUTINE md_output(md_env, md_section, root_section, forced_io)
125 : TYPE(md_environment_type), POINTER :: md_env
126 : TYPE(section_vals_type), POINTER :: md_section, root_section
127 : LOGICAL, INTENT(IN) :: forced_io
128 :
129 : CHARACTER(LEN=*), PARAMETER :: routineN = 'md_output'
130 :
131 : INTEGER :: handle
132 : LOGICAL :: do_print
133 : TYPE(section_vals_type), POINTER :: print_section
134 :
135 42347 : CALL timeset(routineN, handle)
136 42347 : do_print = .TRUE.
137 42347 : IF (forced_io) THEN
138 46 : print_section => section_vals_get_subs_vals(md_section, "PRINT")
139 46 : CALL section_vals_val_get(print_section, "FORCE_LAST", l_val=do_print)
140 : END IF
141 42347 : IF (do_print) THEN
142 : ! Dumps all files related to the MD run
143 42343 : CALL md_write_output(md_env)
144 42343 : CALL write_restart(md_env=md_env, root_section=root_section)
145 : END IF
146 42347 : CALL timestop(handle)
147 :
148 42347 : END SUBROUTINE md_output
149 :
150 : ! **************************************************************************************************
151 : !> \brief read eigenvalues and eigenvectors of Hessian from vibrational analysis results, for use
152 : !> of initialising MD simulations. Expects to read an unformatted binary file
153 : !> \param md_section : input section object containing MD subsections and keywords. This should
154 : !> provide the filename to read vib analysis eigenvalues and eigenvectors.
155 : !> If the filename is not explicitly specified by the user in the input, then
156 : !> it will use the default CARTESIAN_EIGS print key filename defined in the
157 : !> vibrational analysis input section as the filename.
158 : !> \param vib_section : input section object containing vibrational analysis subsections
159 : !> and keywords
160 : !> \param para_env : cp2k mpi environment object, needed for IO in parallel computations
161 : !> \param dof : outputs the total number of eigenvalues (no. degrees of freedom) read from the file
162 : !> \param eigenvalues : outputs the eigenvalues (Cartesian frequencies) read from the file
163 : !> \param eigenvectors : outputs the corresponding eigenvectors read from the file
164 : !> \author Lianheng Tong, lianheng.tong@kcl.ac.uk
165 : ! **************************************************************************************************
166 4 : SUBROUTINE read_vib_eigs_unformatted(md_section, &
167 : vib_section, &
168 : para_env, &
169 : dof, &
170 2 : eigenvalues, &
171 2 : eigenvectors)
172 : TYPE(section_vals_type), POINTER :: md_section, vib_section
173 : TYPE(mp_para_env_type), POINTER :: para_env
174 : INTEGER, INTENT(OUT) :: dof
175 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
176 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eigenvectors
177 :
178 : CHARACTER(LEN=default_path_length) :: filename
179 : INTEGER :: jj, n_rep_val, unit_nr
180 : LOGICAL :: exist
181 : TYPE(cp_logger_type), POINTER :: logger
182 : TYPE(section_vals_type), POINTER :: print_key
183 :
184 2 : logger => cp_get_default_logger()
185 2 : dof = 0
186 20 : eigenvalues = 0.0_dp
187 182 : eigenvectors = 0.0_dp
188 : ! obtain file name
189 : CALL section_vals_val_get(md_section, "INITIAL_VIBRATION%VIB_EIGS_FILE_NAME", &
190 2 : n_rep_val=n_rep_val)
191 2 : IF (n_rep_val > 0) THEN
192 2 : CALL section_vals_val_get(md_section, "INITIAL_VIBRATION%VIB_EIGS_FILE_NAME", c_val=filename)
193 : ELSE
194 0 : print_key => section_vals_get_subs_vals(vib_section, "PRINT%CARTESIAN_EIGS")
195 : filename = cp_print_key_generate_filename(logger, print_key, extension="eig", &
196 0 : my_local=.FALSE.)
197 : END IF
198 : ! read file
199 2 : IF (para_env%is_source()) THEN
200 1 : INQUIRE (FILE=filename, exist=exist)
201 1 : IF (.NOT. exist) THEN
202 0 : CPABORT("File "//filename//" is not found.")
203 : END IF
204 : CALL open_file(file_name=filename, &
205 : file_action="READ", &
206 : file_form="UNFORMATTED", &
207 : file_status="OLD", &
208 1 : unit_number=unit_nr)
209 : ! the first record contains one integer giving degrees of freedom
210 1 : READ (unit_nr) dof
211 1 : IF (dof > SIZE(eigenvalues)) THEN
212 0 : CPABORT("Too many DoFs found in "//filename)
213 : END IF
214 : ! the second record contains the eigenvalues
215 1 : READ (unit_nr) eigenvalues(1:dof)
216 : ! the rest of the records contain the eigenvectors
217 10 : DO jj = 1, dof
218 10 : READ (unit_nr) eigenvectors(1:dof, jj)
219 : END DO
220 : END IF
221 : ! broadcast to all compulational nodes. note that it is assumed
222 : ! that source is the ionode
223 2 : CALL para_env%bcast(dof)
224 38 : CALL para_env%bcast(eigenvalues)
225 362 : CALL para_env%bcast(eigenvectors)
226 : ! close file
227 2 : IF (para_env%is_source()) THEN
228 1 : CALL close_file(unit_number=unit_nr)
229 : END IF
230 2 : END SUBROUTINE read_vib_eigs_unformatted
231 :
232 : END MODULE md_util
|
