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 Optional native tblite command-line diagnostics.
10 : !> \note This module has no dependency on the tblite library.
11 : ! **************************************************************************************************
12 :
13 : MODULE tblite_reference_cli
14 :
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind
17 : USE cell_types, ONLY: cell_type,&
18 : get_cell
19 : USE cp_control_types, ONLY: dft_control_type,&
20 : xtb_reference_cli_type
21 : USE cp_files, ONLY: close_file,&
22 : open_file
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_get_default_io_unit,&
25 : cp_logger_type
26 : USE cp_output_handling, ONLY: debug_print_level,&
27 : high_print_level,&
28 : silent_print_level
29 : USE cp_units, ONLY: cp_unit_from_cp2k
30 : USE input_constants, ONLY: &
31 : gfn1xtb, gfn2xtb, ipea1xtb, smear_energy_window, smear_fermi_dirac, smear_gaussian, &
32 : smear_list, smear_mp, smear_mv, tblite_cli_born_kernel_auto, tblite_cli_born_kernel_p16, &
33 : tblite_cli_born_kernel_still, tblite_cli_solution_state_bar1mol, &
34 : tblite_cli_solution_state_gsolv, tblite_cli_solution_state_reference, &
35 : tblite_cli_solvation_alpb, tblite_cli_solvation_cpcm, tblite_cli_solvation_gb, &
36 : tblite_cli_solvation_gbe, tblite_cli_solvation_gbsa, tblite_guess_ceh, tblite_guess_eeq, &
37 : tblite_guess_sad, tblite_mixer_damping_default, tblite_mixer_max_weight_default, &
38 : tblite_mixer_min_weight_default, tblite_mixer_omega0_default, &
39 : tblite_mixer_weight_factor_default, tblite_scc_mixer_cp2k, tblite_solver_gvd, &
40 : tblite_solver_gvr
41 : USE kinds, ONLY: default_path_length,&
42 : dp
43 : USE kpoint_types, ONLY: get_kpoint_info,&
44 : kpoint_type
45 : USE machine, ONLY: m_unlink
46 : USE message_passing, ONLY: mp_para_env_type
47 : USE particle_types, ONLY: particle_type
48 : USE qs_energy_types, ONLY: qs_energy_type
49 : USE qs_environment_types, ONLY: get_qs_env,&
50 : qs_environment_type
51 : USE qs_force_types, ONLY: qs_force_type,&
52 : total_qs_force
53 : USE scf_control_types, ONLY: scf_control_type
54 : USE virial_types, ONLY: virial_type
55 : #include "../base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 : PRIVATE
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tblite_reference_cli'
61 :
62 : PUBLIC :: tb_reference_cli_compare
63 :
64 : CONTAINS
65 :
66 : ! **************************************************************************************************
67 : !> \brief Run native tblite CLI and compare against CP2K/tblite.
68 : !> \param qs_env ...
69 : ! **************************************************************************************************
70 154 : SUBROUTINE tb_reference_cli_compare(qs_env)
71 :
72 : TYPE(qs_environment_type), POINTER :: qs_env
73 :
74 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tb_reference_cli_compare'
75 :
76 : CHARACTER(LEN=16) :: solvation_model_name
77 : CHARACTER(LEN=32) :: acc_str, charge_str, efield_x_str, efield_y_str, efield_z_str, &
78 : etemp_guess_val_str, etemp_str, iter_str, spin_str, spinpol_str
79 : CHARACTER(LEN=4*default_path_length+128) :: efield_str, etemp_guess_str, param_str, &
80 : post_processing_output_str, post_processing_str, restart_str, solvation_str, verbosity_str
81 : CHARACTER(LEN=64*default_path_length) :: command
82 : CHARACTER(LEN=8) :: guess, method, solver
83 : CHARACTER(LEN=default_path_length) :: file_base, gen_file, grad_file, &
84 : json_file, log_file, &
85 : post_processing_output_file
86 : INTEGER :: handle, iounit, n_periodic, natom, nkp, &
87 : reference_iterations, spin
88 : INTEGER, DIMENSION(3) :: periodic
89 : LOGICAL :: do_kpoints, have_energy, have_gradient, &
90 : have_virial, success, too_large, &
91 : unsupported_kpoints
92 : REAL(KIND=dp) :: cli_energy, cp_energy, ediff, etemp, &
93 : etemp_guess, fmax, fsum, vmax, vsum
94 154 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cli_gradient, cli_virial, cp_gradient
95 154 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
96 154 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
97 : TYPE(cell_type), POINTER :: cell
98 : TYPE(cp_logger_type), POINTER :: logger
99 : TYPE(dft_control_type), POINTER :: dft_control
100 : TYPE(kpoint_type), POINTER :: kpoints
101 : TYPE(mp_para_env_type), POINTER :: para_env
102 154 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
103 : TYPE(qs_energy_type), POINTER :: energy
104 154 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
105 : TYPE(scf_control_type), POINTER :: scf_control
106 : TYPE(virial_type), POINTER :: virial
107 : TYPE(xtb_reference_cli_type) :: ref
108 :
109 154 : CALL timeset(routineN, handle)
110 :
111 154 : NULLIFY (atomic_kind_set, cell, dft_control, energy, force, kpoints, logger, para_env, particle_set, &
112 154 : scf_control, virial, xkp)
113 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell, &
114 : dft_control=dft_control, energy=energy, force=force, &
115 : para_env=para_env, particle_set=particle_set, scf_control=scf_control, &
116 154 : virial=virial, do_kpoints=do_kpoints, kpoints=kpoints)
117 :
118 154 : ref = dft_control%qs_control%xtb_control%reference_cli
119 154 : IF (.NOT. ref%enabled) THEN
120 152 : CALL timestop(handle)
121 152 : RETURN
122 : END IF
123 2 : IF (.NOT. para_env%is_source()) THEN
124 1 : CALL timestop(handle)
125 1 : RETURN
126 : END IF
127 :
128 1 : logger => cp_get_default_logger()
129 1 : iounit = cp_logger_get_default_io_unit(logger)
130 1 : verbosity_str = ""
131 1 : IF (logger%iter_info%print_level == silent_print_level) THEN
132 0 : verbosity_str = " --silent"
133 1 : ELSE IF (logger%iter_info%print_level == high_print_level .OR. &
134 : logger%iter_info%print_level == debug_print_level) THEN
135 0 : verbosity_str = " --verbose"
136 : END IF
137 1 : IF (ref%solvation_active) THEN
138 1 : periodic = 0
139 1 : IF (ASSOCIATED(cell)) CALL get_cell(cell=cell, periodic=periodic)
140 4 : n_periodic = COUNT(periodic == 1)
141 1 : IF (n_periodic == 3) THEN
142 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
143 0 : "tblite reference CLI implicit solvation is not supported for PERIODIC XYZ."
144 : WRITE (UNIT=iounit, FMT="(T2,A)") &
145 0 : "Use PERIODIC NONE for molecular solvation diagnostics, or remove IMPLICIT_SOLVATION."
146 0 : CPABORT("REFERENCE_CLI implicit solvation is incompatible with PERIODIC XYZ")
147 1 : ELSE IF (n_periodic > 0) THEN
148 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
149 0 : "WARNING: tblite reference CLI implicit solvation with finite periodicity is diagnostic only."
150 : WRITE (UNIT=iounit, FMT="(T2,A,I0,A)") &
151 0 : "The generated native tblite reference geometry has ", n_periodic, &
152 0 : " periodic direction(s); continuum-solvation conventions are primarily molecular."
153 : END IF
154 : END IF
155 1 : unsupported_kpoints = .FALSE.
156 1 : IF (do_kpoints .AND. ASSOCIATED(kpoints)) THEN
157 0 : nkp = 0
158 0 : CALL get_kpoint_info(kpoint=kpoints, nkp=nkp, xkp=xkp)
159 0 : unsupported_kpoints = nkp > 1
160 0 : IF (nkp == 1 .AND. ASSOCIATED(xkp)) unsupported_kpoints = ANY(ABS(xkp(:, 1)) > 1.0E-12_dp)
161 : END IF
162 0 : IF (unsupported_kpoints) THEN
163 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
164 0 : "tblite reference CLI check skipped: CP2K KPOINTS are active."
165 : WRITE (UNIT=iounit, FMT="(T2,A)") &
166 0 : "The native tblite CLI reference path does not reproduce CP2K multi-k-point sampling."
167 0 : IF (ref%stop_on_error) THEN
168 0 : CPABORT("tblite reference CLI cannot check CP2K k-point calculations")
169 : END IF
170 0 : CALL timestop(handle)
171 0 : RETURN
172 : END IF
173 1 : IF (dft_control%qs_control%xtb_control%tblite_scc_mixer == tblite_scc_mixer_cp2k) THEN
174 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
175 0 : "WARNING: tblite reference CLI cannot reproduce XTB/SCC_MIXER CP2K."
176 : WRITE (UNIT=iounit, FMT="(T2,A)") &
177 0 : "The external native tblite run uses tblite's own SCC mixer; only the converged result is compared."
178 0 : IF (ref%stop_on_error) CPABORT("tblite reference CLI cannot reproduce SCC_MIXER CP2K")
179 : END IF
180 1 : natom = SIZE(particle_set)
181 1 : method = tb_reference_method_name(dft_control%qs_control%xtb_control%tblite_method)
182 1 : guess = tb_reference_guess_name(ref%guess)
183 1 : solver = tb_reference_solver_name(dft_control%qs_control%xtb_control%tblite_mixer_solver)
184 1 : file_base = tb_join_path(ref%work_directory, ref%prefix)
185 1 : gen_file = TRIM(file_base)//".gen"
186 1 : grad_file = TRIM(file_base)//".grad"
187 1 : json_file = TRIM(file_base)//".json"
188 1 : log_file = TRIM(file_base)//".log"
189 1 : post_processing_output_file = ""
190 1 : IF (LEN_TRIM(ref%grad_file) > 0) grad_file = ref%grad_file
191 1 : IF (LEN_TRIM(ref%json_file) > 0) json_file = ref%json_file
192 1 : IF (LEN_TRIM(ref%post_processing_output_file) > 0) THEN
193 1 : post_processing_output_file = ref%post_processing_output_file
194 : END IF
195 :
196 1 : WRITE (charge_str, "(I0)") dft_control%charge
197 1 : spin = MAX(0, dft_control%multiplicity - 1)
198 1 : WRITE (spin_str, "(I0)") spin
199 1 : WRITE (acc_str, "(ES16.8)") dft_control%qs_control%xtb_control%tblite_accuracy
200 1 : reference_iterations = dft_control%qs_control%xtb_control%tblite_mixer_iterations
201 1 : WRITE (iter_str, "(I0)") reference_iterations
202 1 : efield_str = ""
203 1 : IF (ref%efield_active) THEN
204 1 : WRITE (efield_x_str, "(ES16.8)") ref%efield(1)
205 1 : WRITE (efield_y_str, "(ES16.8)") ref%efield(2)
206 1 : WRITE (efield_z_str, "(ES16.8)") ref%efield(3)
207 : efield_str = " --efield "//TRIM(ADJUSTL(efield_x_str))//","// &
208 1 : TRIM(ADJUSTL(efield_y_str))//","//TRIM(ADJUSTL(efield_z_str))
209 : END IF
210 1 : solvation_str = ""
211 1 : solvation_model_name = ""
212 1 : IF (ref%solvation_active) THEN
213 2 : SELECT CASE (ref%solvation_model)
214 : CASE (tblite_cli_solvation_alpb)
215 1 : solvation_model_name = "ALPB"
216 1 : solvation_str = " --alpb "
217 : CASE (tblite_cli_solvation_gbsa)
218 0 : solvation_model_name = "GBSA"
219 0 : solvation_str = " --gbsa "
220 : CASE (tblite_cli_solvation_gbe)
221 0 : solvation_model_name = "GBE"
222 0 : solvation_str = " --gbe "
223 : CASE (tblite_cli_solvation_gb)
224 0 : solvation_model_name = "GB"
225 0 : solvation_str = " --gb "
226 : CASE (tblite_cli_solvation_cpcm)
227 0 : solvation_model_name = "CPCM"
228 0 : solvation_str = " --cpcm "
229 : CASE DEFAULT
230 1 : CPABORT("Unknown tblite reference CLI implicit-solvation model")
231 : END SELECT
232 1 : solvation_str = TRIM(solvation_str)//" "//TRIM(tb_shell_quote(ref%solvation_solvent))
233 2 : SELECT CASE (ref%solvation_born_kernel)
234 : CASE (tblite_cli_born_kernel_auto)
235 : CASE (tblite_cli_born_kernel_p16)
236 1 : solvation_str = TRIM(solvation_str)//" --born-kernel p16"
237 : CASE (tblite_cli_born_kernel_still)
238 0 : solvation_str = TRIM(solvation_str)//" --born-kernel still"
239 : CASE DEFAULT
240 1 : CPABORT("Unknown tblite reference CLI Born kernel")
241 : END SELECT
242 2 : SELECT CASE (ref%solvation_state)
243 : CASE (tblite_cli_solution_state_gsolv)
244 : CASE (tblite_cli_solution_state_bar1mol)
245 1 : solvation_str = TRIM(solvation_str)//" --solv-state bar1mol"
246 : CASE (tblite_cli_solution_state_reference)
247 0 : solvation_str = TRIM(solvation_str)//" --solv-state reference"
248 : CASE DEFAULT
249 1 : CPABORT("Unknown tblite reference CLI solution state")
250 : END SELECT
251 : END IF
252 1 : IF (dft_control%qs_control%xtb_control%tblite_mixer_memory /= reference_iterations) THEN
253 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
254 0 : "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/MEMORY."
255 : WRITE (UNIT=iounit, FMT="(T2,A,I0,A,I0,A)") &
256 0 : "The native reference run uses tblite's internal mixer memory tied to --iterations (", &
257 0 : reference_iterations, "), while CP2K uses MEMORY ", &
258 0 : dft_control%qs_control%xtb_control%tblite_mixer_memory, "."
259 0 : IF (ref%stop_on_error) THEN
260 0 : CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/MEMORY")
261 : END IF
262 : END IF
263 1 : IF (ABS(dft_control%qs_control%xtb_control%tblite_mixer_damping - &
264 : tblite_mixer_damping_default) > 10.0_dp*EPSILON(1.0_dp)) THEN
265 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
266 0 : "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/DAMPING."
267 : WRITE (UNIT=iounit, FMT="(T2,A,F8.4,A,F8.4,A)") &
268 0 : "The native reference run uses tblite's library default ", tblite_mixer_damping_default, &
269 0 : ", while CP2K uses DAMPING ", dft_control%qs_control%xtb_control%tblite_mixer_damping, "."
270 0 : IF (ref%stop_on_error) THEN
271 0 : CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/DAMPING")
272 : END IF
273 : END IF
274 1 : IF (ABS(dft_control%qs_control%xtb_control%tblite_mixer_omega0 - &
275 : tblite_mixer_omega0_default) > 10.0_dp*EPSILON(1.0_dp)) THEN
276 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
277 0 : "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/OMEGA0."
278 : WRITE (UNIT=iounit, FMT="(T2,A,ES12.4,A,ES12.4,A)") &
279 0 : "The native reference run uses tblite's library default ", tblite_mixer_omega0_default, &
280 0 : ", while CP2K uses OMEGA0 ", dft_control%qs_control%xtb_control%tblite_mixer_omega0, "."
281 0 : IF (ref%stop_on_error) THEN
282 0 : CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/OMEGA0")
283 : END IF
284 : END IF
285 1 : IF (ABS(dft_control%qs_control%xtb_control%tblite_mixer_min_weight - &
286 : tblite_mixer_min_weight_default) > 10.0_dp*EPSILON(1.0_dp)) THEN
287 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
288 0 : "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/MIN_WEIGHT."
289 : WRITE (UNIT=iounit, FMT="(T2,A,ES12.4,A,ES12.4,A)") &
290 0 : "The native reference run uses tblite's library default ", tblite_mixer_min_weight_default, &
291 0 : ", while CP2K uses MIN_WEIGHT ", dft_control%qs_control%xtb_control%tblite_mixer_min_weight, "."
292 0 : IF (ref%stop_on_error) THEN
293 0 : CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/MIN_WEIGHT")
294 : END IF
295 : END IF
296 1 : IF (ABS(dft_control%qs_control%xtb_control%tblite_mixer_max_weight - &
297 : tblite_mixer_max_weight_default) > 10.0_dp*EPSILON(1.0_dp)) THEN
298 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
299 0 : "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/MAX_WEIGHT."
300 : WRITE (UNIT=iounit, FMT="(T2,A,ES12.4,A,ES12.4,A)") &
301 0 : "The native reference run uses tblite's library default ", tblite_mixer_max_weight_default, &
302 0 : ", while CP2K uses MAX_WEIGHT ", dft_control%qs_control%xtb_control%tblite_mixer_max_weight, "."
303 0 : IF (ref%stop_on_error) THEN
304 0 : CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/MAX_WEIGHT")
305 : END IF
306 : END IF
307 1 : IF (ABS(dft_control%qs_control%xtb_control%tblite_mixer_weight_factor - &
308 : tblite_mixer_weight_factor_default) > 10.0_dp*EPSILON(1.0_dp)) THEN
309 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
310 0 : "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/WEIGHT_FACTOR."
311 : WRITE (UNIT=iounit, FMT="(T2,A,ES12.4,A,ES12.4,A)") &
312 0 : "The native reference run uses tblite's library default ", tblite_mixer_weight_factor_default, &
313 0 : ", while CP2K uses WEIGHT_FACTOR ", &
314 0 : dft_control%qs_control%xtb_control%tblite_mixer_weight_factor, "."
315 0 : IF (ref%stop_on_error) THEN
316 0 : CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/WEIGHT_FACTOR")
317 : END IF
318 : END IF
319 1 : etemp = 300.0_dp
320 1 : IF (ASSOCIATED(scf_control)) THEN
321 1 : IF (ASSOCIATED(scf_control%smear)) THEN
322 1 : IF (scf_control%smear%do_smear) THEN
323 1 : etemp = cp_unit_from_cp2k(scf_control%smear%electronic_temperature, "K")
324 1 : IF (scf_control%smear%method /= smear_fermi_dirac) THEN
325 : WRITE (UNIT=iounit, FMT="(/,T2,A,A,A)") &
326 0 : "WARNING: tblite reference CLI cannot reproduce CP2K smearing method ", &
327 0 : TRIM(tb_reference_smear_method_name(scf_control%smear%method)), "."
328 : WRITE (UNIT=iounit, FMT="(T2,A,F12.3,A)") &
329 0 : "The native reference run uses Fermi-Dirac electronic temperature ", etemp, " K instead."
330 : END IF
331 : END IF
332 : END IF
333 : END IF
334 1 : WRITE (etemp_str, "(ES16.8)") etemp
335 1 : etemp_guess = 0.0_dp
336 1 : etemp_guess_str = ""
337 1 : IF (ref%electronic_temperature_guess > 0.0_dp) THEN
338 0 : etemp_guess = cp_unit_from_cp2k(ref%electronic_temperature_guess, "K")
339 0 : WRITE (etemp_guess_val_str, "(ES16.8)") etemp_guess
340 0 : etemp_guess_str = " --etemp-guess "//TRIM(ADJUSTL(etemp_guess_val_str))
341 : END IF
342 1 : param_str = ""
343 1 : IF (LEN_TRIM(dft_control%qs_control%xtb_control%tblite_param_file) > 0) THEN
344 0 : param_str = " --param "//TRIM(tb_shell_quote(dft_control%qs_control%xtb_control%tblite_param_file))
345 : END IF
346 1 : spinpol_str = ""
347 1 : IF (dft_control%lsd) spinpol_str = " --spin-polarized"
348 1 : post_processing_str = ""
349 1 : IF (LEN_TRIM(ref%post_processing) > 0) THEN
350 1 : post_processing_str = " --post-processing "//TRIM(tb_shell_quote(ref%post_processing))
351 : END IF
352 1 : post_processing_output_str = ""
353 1 : IF (LEN_TRIM(post_processing_output_file) > 0) THEN
354 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
355 1 : "WARNING: tblite reference CLI POST_PROCESSING_OUTPUT was requested explicitly."
356 : WRITE (UNIT=iounit, FMT="(T2,A)") &
357 1 : "Some tblite 0.5.0 command-line builds document --post-processing-output but do not parse it."
358 : post_processing_output_str = " --post-processing-output "// &
359 1 : TRIM(tb_shell_quote(post_processing_output_file))
360 : END IF
361 1 : restart_str = " --no-restart"
362 1 : IF (LEN_TRIM(ref%restart_file) > 0) THEN
363 0 : restart_str = " --restart "//TRIM(tb_shell_quote(ref%restart_file))
364 : END IF
365 :
366 1 : CALL tb_write_reference_gen(qs_env, TRIM(gen_file))
367 :
368 : command = TRIM(tb_shell_quote(ref%program_name))//" run --method "//TRIM(method)// &
369 : TRIM(param_str)// &
370 : TRIM(spinpol_str)// &
371 : " --charge "//TRIM(ADJUSTL(charge_str))// &
372 : " --spin "//TRIM(ADJUSTL(spin_str))// &
373 : " --acc "//TRIM(ADJUSTL(acc_str))// &
374 : " --guess "//TRIM(guess)// &
375 : " --solver "//TRIM(solver)// &
376 : " --iterations "//TRIM(ADJUSTL(iter_str))// &
377 : " --etemp "//TRIM(ADJUSTL(etemp_str))// &
378 : TRIM(etemp_guess_str)// &
379 : TRIM(efield_str)// &
380 : TRIM(solvation_str)// &
381 : TRIM(post_processing_str)// &
382 : TRIM(post_processing_output_str)// &
383 : TRIM(restart_str)// &
384 : TRIM(verbosity_str)//" --input "//TRIM(tb_shell_quote(ref%input_format))// &
385 : " --grad "//TRIM(tb_shell_quote(grad_file))// &
386 : " --json "//TRIM(tb_shell_quote(json_file))//" "//TRIM(tb_shell_quote(gen_file))// &
387 1 : " > "//TRIM(tb_shell_quote(log_file))//" 2>&1"
388 :
389 1 : CALL tb_reference_cli_execute(ref, "run", command, log_file, iounit, success)
390 1 : IF (.NOT. success) THEN
391 : ! Retain the geometry and diagnostic output on failure, even without KEEP_FILES.
392 1 : CALL timestop(handle)
393 1 : RETURN
394 : END IF
395 :
396 0 : ALLOCATE (cli_gradient(3, natom), cli_virial(3, 3))
397 : CALL tb_read_reference_grad(TRIM(grad_file), natom, cli_energy, cli_gradient, cli_virial, &
398 0 : have_energy, have_gradient, have_virial)
399 :
400 0 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "tblite reference CLI check"
401 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Executable: ", TRIM(ref%program_name)
402 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Method: ", TRIM(method)
403 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Guess: ", TRIM(guess)
404 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Solver: ", TRIM(solver)
405 0 : WRITE (UNIT=iounit, FMT="(T2,A,L1)") "Spin-pol.: ", dft_control%lsd
406 0 : IF (ref%efield_active) THEN
407 0 : WRITE (UNIT=iounit, FMT="(T2,A,3ES16.8,A)") "Efield: ", ref%efield, " V/Angstrom"
408 : END IF
409 0 : IF (ref%solvation_active) THEN
410 0 : WRITE (UNIT=iounit, FMT="(T2,A,A,1X,A)") "Solvation: ", TRIM(solvation_model_name), &
411 0 : TRIM(ref%solvation_solvent)
412 : END IF
413 0 : IF (LEN_TRIM(ref%post_processing) > 0) THEN
414 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Post proc.: ", TRIM(ref%post_processing)
415 : END IF
416 0 : IF (LEN_TRIM(post_processing_output_file) > 0) THEN
417 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "PP output: ", TRIM(post_processing_output_file)
418 : END IF
419 0 : IF (ref%electronic_temperature_guess > 0.0_dp) THEN
420 0 : WRITE (UNIT=iounit, FMT="(T2,A,F12.3,A)") "Guess etemp:", etemp_guess, " K"
421 : END IF
422 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Grad file: ", TRIM(grad_file)
423 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "JSON file: ", TRIM(json_file)
424 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Log file: ", TRIM(log_file)
425 :
426 0 : too_large = .FALSE.
427 0 : IF (ref%check_energy) THEN
428 0 : IF (have_energy) THEN
429 0 : cp_energy = energy%total
430 0 : ediff = ABS(cp_energy - cli_energy)
431 : WRITE (UNIT=iounit, FMT="(T2,A,3ES22.12)") &
432 0 : "Energy CP2K/CLI/absdiff:", cp_energy, cli_energy, ediff
433 0 : too_large = too_large .OR. ediff > ref%error_limit
434 : ELSE
435 0 : WRITE (UNIT=iounit, FMT="(T2,A)") "Energy check skipped: no CLI energy found."
436 : END IF
437 : END IF
438 :
439 0 : IF (ref%check_forces) THEN
440 0 : IF (have_gradient .AND. ASSOCIATED(force)) THEN
441 0 : ALLOCATE (cp_gradient(3, natom))
442 0 : CALL total_qs_force(cp_gradient, force, atomic_kind_set)
443 0 : fsum = SUM(ABS(cp_gradient - cli_gradient))
444 0 : fmax = MAXVAL(ABS(cp_gradient - cli_gradient))
445 0 : WRITE (UNIT=iounit, FMT="(T2,A,2ES22.12)") "Gradient diff sum/max:", fsum, fmax
446 0 : too_large = too_large .OR. fmax > ref%error_limit
447 0 : DEALLOCATE (cp_gradient)
448 : ELSE
449 0 : WRITE (UNIT=iounit, FMT="(T2,A)") "Gradient check skipped: no CLI gradient or CP2K force found."
450 : END IF
451 : END IF
452 :
453 0 : IF (ref%check_virial) THEN
454 0 : IF (have_virial .AND. ASSOCIATED(virial)) THEN
455 : ! Native tblite prints the positive cell derivative; CP2K stores the PV virial
456 : ! with the opposite sign.
457 0 : vsum = SUM(ABS(-virial%pv_virial - cli_virial))
458 0 : vmax = MAXVAL(ABS(-virial%pv_virial - cli_virial))
459 0 : WRITE (UNIT=iounit, FMT="(T2,A,2ES22.12)") "Virial diff sum/max:", vsum, vmax
460 0 : too_large = too_large .OR. vmax > ref%error_limit
461 : ELSE
462 0 : WRITE (UNIT=iounit, FMT="(T2,A)") "Virial check skipped: no CLI virial or CP2K virial found."
463 : END IF
464 : END IF
465 :
466 0 : IF (too_large) THEN
467 : WRITE (UNIT=iounit, FMT="(T2,A,ES12.4)") &
468 0 : "tblite reference CLI deviation exceeded ERROR_LIMIT = ", ref%error_limit
469 0 : IF (ref%stop_on_error) CPABORT("tblite reference CLI deviation exceeded ERROR_LIMIT")
470 : END IF
471 :
472 0 : CALL tb_reference_cli_aux_commands(ref, dft_control, gen_file, file_base, verbosity_str, iounit, success)
473 :
474 0 : IF (success .AND. .NOT. too_large) THEN
475 0 : CALL tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file, post_processing_output_file)
476 : END IF
477 0 : DEALLOCATE (cli_gradient, cli_virial)
478 :
479 0 : CALL timestop(handle)
480 :
481 1233 : END SUBROUTINE tb_reference_cli_compare
482 :
483 : ! **************************************************************************************************
484 : !> \brief Map CP2K tblite method id to native tblite CLI method name.
485 : !> \param method_id ...
486 : !> \return ...
487 : ! **************************************************************************************************
488 1 : FUNCTION tb_reference_method_name(method_id) RESULT(method)
489 : INTEGER, INTENT(IN) :: method_id
490 : CHARACTER(LEN=8) :: method
491 :
492 1 : SELECT CASE (method_id)
493 : CASE (gfn1xtb)
494 0 : method = "gfn1"
495 : CASE (gfn2xtb)
496 1 : method = "gfn2"
497 : CASE (ipea1xtb)
498 0 : method = "ipea1"
499 : CASE DEFAULT
500 1 : CPABORT("Unknown tblite reference CLI method")
501 : END SELECT
502 :
503 1 : END FUNCTION tb_reference_method_name
504 :
505 : ! **************************************************************************************************
506 : !> \brief Map CP2K tblite reference CLI guess id to native tblite CLI guess name.
507 : !> \param guess_id ...
508 : !> \return ...
509 : ! **************************************************************************************************
510 1 : FUNCTION tb_reference_guess_name(guess_id) RESULT(guess)
511 : INTEGER, INTENT(IN) :: guess_id
512 : CHARACTER(LEN=8) :: guess
513 :
514 2 : SELECT CASE (guess_id)
515 : CASE (tblite_guess_sad)
516 1 : guess = "sad"
517 : CASE (tblite_guess_eeq)
518 0 : guess = "eeq"
519 : CASE (tblite_guess_ceh)
520 0 : guess = "ceh"
521 : CASE DEFAULT
522 1 : CPABORT("Unknown tblite reference CLI guess")
523 : END SELECT
524 :
525 1 : END FUNCTION tb_reference_guess_name
526 :
527 : ! **************************************************************************************************
528 : !> \brief Map CP2K tblite solver id to native tblite CLI solver name.
529 : !> \param solver_id ...
530 : !> \return ...
531 : ! **************************************************************************************************
532 1 : FUNCTION tb_reference_solver_name(solver_id) RESULT(solver)
533 : INTEGER, INTENT(IN) :: solver_id
534 : CHARACTER(LEN=8) :: solver
535 :
536 2 : SELECT CASE (solver_id)
537 : CASE (tblite_solver_gvd)
538 1 : solver = "gvd"
539 : CASE (tblite_solver_gvr)
540 0 : solver = "gvr"
541 : CASE DEFAULT
542 1 : CPABORT("Unknown tblite reference CLI solver")
543 : END SELECT
544 :
545 1 : END FUNCTION tb_reference_solver_name
546 :
547 : ! **************************************************************************************************
548 : !> \brief Map CP2K smearing method id to a diagnostic label.
549 : !> \param method_id ...
550 : !> \return ...
551 : ! **************************************************************************************************
552 0 : FUNCTION tb_reference_smear_method_name(method_id) RESULT(method)
553 : INTEGER, INTENT(IN) :: method_id
554 : CHARACTER(LEN=24) :: method
555 :
556 0 : SELECT CASE (method_id)
557 : CASE (smear_fermi_dirac)
558 0 : method = "FERMI_DIRAC"
559 : CASE (smear_energy_window)
560 0 : method = "ENERGY_WINDOW"
561 : CASE (smear_list)
562 0 : method = "LIST"
563 : CASE (smear_gaussian)
564 0 : method = "GAUSSIAN"
565 : CASE (smear_mp)
566 0 : method = "METHFESSEL_PAXTON"
567 : CASE (smear_mv)
568 0 : method = "MARZARI_VANDERBILT"
569 : CASE DEFAULT
570 0 : method = "UNKNOWN"
571 : END SELECT
572 :
573 0 : END FUNCTION tb_reference_smear_method_name
574 :
575 : ! **************************************************************************************************
576 : !> \brief Run optional native tblite auxiliary subcommands.
577 : !> \param ref ...
578 : !> \param dft_control ...
579 : !> \param gen_file ...
580 : !> \param file_base ...
581 : !> \param verbosity_str ...
582 : !> \param iounit ...
583 : !> \param success Whether all requested commands completed successfully.
584 : ! **************************************************************************************************
585 0 : SUBROUTINE tb_reference_cli_aux_commands(ref, dft_control, gen_file, file_base, verbosity_str, iounit, success)
586 :
587 : TYPE(xtb_reference_cli_type), INTENT(IN) :: ref
588 : TYPE(dft_control_type), INTENT(IN) :: dft_control
589 : CHARACTER(LEN=*), INTENT(IN) :: gen_file, file_base, verbosity_str
590 : INTEGER, INTENT(IN) :: iounit
591 : LOGICAL, INTENT(OUT) :: success
592 :
593 : CHARACTER(LEN=32) :: charge_str, efield_x_str, efield_y_str, &
594 : efield_z_str, etemp_guess_val_str, &
595 : spin_str
596 : CHARACTER(LEN=4*default_path_length+128) :: copy_str, dry_run_str, efield_str, &
597 : etemp_guess_str, grad_str, json_str, &
598 : method_str, output_str
599 : CHARACTER(LEN=64*default_path_length) :: command
600 : CHARACTER(LEN=default_path_length) :: guess_input, log_file
601 : INTEGER :: spin
602 : LOGICAL :: command_success
603 : REAL(KIND=dp) :: etemp_guess
604 :
605 0 : success = .TRUE.
606 0 : WRITE (charge_str, "(I0)") dft_control%charge
607 0 : spin = MAX(0, dft_control%multiplicity - 1)
608 0 : WRITE (spin_str, "(I0)") spin
609 :
610 0 : IF (ref%guess_cli%enabled) THEN
611 0 : guess_input = ref%guess_cli%input_file
612 0 : IF (LEN_TRIM(guess_input) == 0) guess_input = gen_file
613 0 : etemp_guess_str = ""
614 0 : IF (ref%guess_cli%electronic_temperature_guess > 0.0_dp) THEN
615 0 : etemp_guess = cp_unit_from_cp2k(ref%guess_cli%electronic_temperature_guess, "K")
616 0 : WRITE (etemp_guess_val_str, "(ES16.8)") etemp_guess
617 0 : etemp_guess_str = " --etemp-guess "//TRIM(ADJUSTL(etemp_guess_val_str))
618 : END IF
619 0 : efield_str = ""
620 0 : IF (ref%guess_cli%efield_active) THEN
621 0 : WRITE (efield_x_str, "(ES16.8)") ref%guess_cli%efield(1)
622 0 : WRITE (efield_y_str, "(ES16.8)") ref%guess_cli%efield(2)
623 0 : WRITE (efield_z_str, "(ES16.8)") ref%guess_cli%efield(3)
624 : efield_str = " --efield "//TRIM(ADJUSTL(efield_x_str))//","// &
625 0 : TRIM(ADJUSTL(efield_y_str))//","//TRIM(ADJUSTL(efield_z_str))
626 : END IF
627 0 : grad_str = ""
628 0 : IF (ref%guess_cli%grad) grad_str = " --grad"
629 0 : json_str = ""
630 0 : IF (LEN_TRIM(ref%guess_cli%json_file) > 0) THEN
631 0 : json_str = " --json "//TRIM(tb_shell_quote(ref%guess_cli%json_file))
632 : END IF
633 0 : log_file = TRIM(file_base)//".guess.log"
634 : command = TRIM(tb_shell_quote(ref%program_name))// &
635 : " guess --charge "//TRIM(ADJUSTL(charge_str))// &
636 : " --spin "//TRIM(ADJUSTL(spin_str))// &
637 : " --method "//TRIM(tb_reference_guess_name(ref%guess_cli%method))// &
638 : " --solver "//TRIM(tb_reference_solver_name(ref%guess_cli%solver))// &
639 : TRIM(etemp_guess_str)// &
640 : TRIM(efield_str)// &
641 : TRIM(grad_str)// &
642 : TRIM(json_str)// &
643 : TRIM(verbosity_str)//" --input "//TRIM(tb_shell_quote(ref%guess_cli%input_format))// &
644 : " "//TRIM(tb_shell_quote(guess_input))// &
645 0 : " > "//TRIM(tb_shell_quote(log_file))//" 2>&1"
646 0 : CALL tb_reference_cli_execute(ref, "guess", command, log_file, iounit, command_success)
647 0 : success = success .AND. command_success
648 0 : IF (command_success .AND. .NOT. ref%keep_files) CALL tb_delete_file(log_file)
649 0 : IF (command_success .AND. .NOT. ref%keep_files .AND. LEN_TRIM(ref%guess_cli%json_file) > 0) THEN
650 0 : CALL tb_delete_file(ref%guess_cli%json_file)
651 : END IF
652 : END IF
653 :
654 0 : IF (ref%param_cli%enabled) THEN
655 0 : method_str = ""
656 0 : IF (ref%param_cli%method_explicit .OR. LEN_TRIM(ref%param_cli%input_file) == 0) THEN
657 : method_str = " --method "// &
658 : TRIM(tb_reference_method_name(MERGE(ref%param_cli%method, &
659 : dft_control%qs_control%xtb_control%tblite_method, &
660 0 : ref%param_cli%method_explicit)))
661 : END IF
662 0 : output_str = ""
663 0 : IF (LEN_TRIM(ref%param_cli%output_file) > 0) THEN
664 0 : output_str = " --output "//TRIM(tb_shell_quote(ref%param_cli%output_file))
665 : END IF
666 0 : log_file = TRIM(file_base)//".param.log"
667 : command = TRIM(tb_shell_quote(ref%program_name))//" param"// &
668 : TRIM(method_str)// &
669 0 : TRIM(output_str)
670 0 : IF (LEN_TRIM(ref%param_cli%input_file) > 0) THEN
671 0 : command = TRIM(command)//" "//TRIM(tb_shell_quote(ref%param_cli%input_file))
672 : END IF
673 0 : command = TRIM(command)//" > "//TRIM(tb_shell_quote(log_file))//" 2>&1"
674 0 : CALL tb_reference_cli_execute(ref, "param", command, log_file, iounit, command_success)
675 0 : success = success .AND. command_success
676 0 : IF (command_success .AND. .NOT. ref%keep_files) CALL tb_delete_file(log_file)
677 0 : IF (command_success .AND. .NOT. ref%keep_files .AND. LEN_TRIM(ref%param_cli%output_file) > 0) THEN
678 0 : CALL tb_delete_file(ref%param_cli%output_file)
679 : END IF
680 : END IF
681 :
682 0 : IF (ref%fit_cli%enabled) THEN
683 0 : dry_run_str = ""
684 0 : IF (ref%fit_cli%dry_run) dry_run_str = " --dry-run"
685 0 : copy_str = ""
686 0 : IF (LEN_TRIM(ref%fit_cli%copy_file) > 0) THEN
687 0 : copy_str = " --copy "//TRIM(tb_shell_quote(ref%fit_cli%copy_file))
688 : END IF
689 0 : log_file = TRIM(file_base)//".fit.log"
690 : command = TRIM(tb_shell_quote(ref%program_name))//" fit"// &
691 : TRIM(dry_run_str)// &
692 : TRIM(copy_str)// &
693 : TRIM(verbosity_str)//" "//TRIM(tb_shell_quote(ref%fit_cli%param_file))// &
694 : " "//TRIM(tb_shell_quote(ref%fit_cli%input_file))// &
695 0 : " > "//TRIM(tb_shell_quote(log_file))//" 2>&1"
696 0 : CALL tb_reference_cli_execute(ref, "fit", command, log_file, iounit, command_success)
697 0 : success = success .AND. command_success
698 0 : IF (command_success .AND. .NOT. ref%keep_files) CALL tb_delete_file(log_file)
699 0 : IF (command_success .AND. .NOT. ref%keep_files .AND. LEN_TRIM(ref%fit_cli%copy_file) > 0) THEN
700 0 : CALL tb_delete_file(ref%fit_cli%copy_file)
701 : END IF
702 : END IF
703 :
704 0 : IF (ref%tagdiff_cli%enabled) THEN
705 0 : method_str = ""
706 0 : IF (ref%tagdiff_cli%fit) method_str = " --fit"
707 0 : log_file = TRIM(file_base)//".tagdiff.log"
708 : command = TRIM(tb_shell_quote(ref%program_name))//" tagdiff"// &
709 : TRIM(method_str)//" "//TRIM(tb_shell_quote(ref%tagdiff_cli%actual_file))// &
710 : " "//TRIM(tb_shell_quote(ref%tagdiff_cli%reference_file))// &
711 0 : " > "//TRIM(tb_shell_quote(log_file))//" 2>&1"
712 0 : CALL tb_reference_cli_execute(ref, "tagdiff", command, log_file, iounit, command_success)
713 0 : success = success .AND. command_success
714 0 : IF (command_success .AND. .NOT. ref%keep_files) CALL tb_delete_file(log_file)
715 : END IF
716 :
717 0 : END SUBROUTINE tb_reference_cli_aux_commands
718 :
719 : ! **************************************************************************************************
720 : !> \brief Execute one native tblite REFERENCE_CLI command.
721 : !> \param ref ...
722 : !> \param label ...
723 : !> \param command ...
724 : !> \param log_file ...
725 : !> \param iounit ...
726 : !> \param success Whether all requested commands completed successfully.
727 : !> \note Keep the external-process exception confined to this optional diagnostic path.
728 : ! **************************************************************************************************
729 1 : SUBROUTINE tb_reference_cli_execute(ref, label, command, log_file, iounit, success)
730 :
731 : TYPE(xtb_reference_cli_type), INTENT(IN) :: ref
732 : CHARACTER(LEN=*), INTENT(IN) :: label, command, log_file
733 : INTEGER, INTENT(IN) :: iounit
734 : LOGICAL, INTENT(OUT) :: success
735 :
736 : CHARACTER(LEN=default_path_length) :: cmdmsg
737 : INTEGER :: cmdstat, exitstat
738 :
739 1 : cmdmsg = ""
740 1 : cmdstat = 0
741 1 : exitstat = 0
742 1 : CALL EXECUTE_COMMAND_LINE(TRIM(command), EXITSTAT=exitstat, CMDSTAT=cmdstat, CMDMSG=cmdmsg)
743 1 : success = cmdstat == 0 .AND. exitstat == 0
744 1 : IF (.NOT. success) THEN
745 1 : WRITE (UNIT=iounit, FMT="(/,T2,A)") "tblite reference CLI check failed to run."
746 1 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Command: ", TRIM(command)
747 1 : WRITE (UNIT=iounit, FMT="(T2,A,I0,T32,A,I0)") "cmdstat:", cmdstat, "exitstat:", exitstat
748 1 : IF (cmdstat /= 0) WRITE (UNIT=iounit, FMT="(T2,A,A)") "Execution error: ", TRIM(cmdmsg)
749 1 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Diagnostic files retained; log: ", TRIM(log_file)
750 1 : IF (ref%stop_on_error) CPABORT("tblite reference CLI command failed")
751 : ELSE
752 0 : WRITE (UNIT=iounit, FMT="(/,T2,A,A)") "tblite reference CLI command completed: ", &
753 0 : TRIM(label)
754 0 : WRITE (UNIT=iounit, FMT="(T2,A,A)") "Log file: ", TRIM(log_file)
755 : END IF
756 :
757 1 : END SUBROUTINE tb_reference_cli_execute
758 :
759 : ! **************************************************************************************************
760 : !> \brief Join directory and filename.
761 : !> \param directory ...
762 : !> \param filename ...
763 : !> \return ...
764 : ! **************************************************************************************************
765 1 : FUNCTION tb_join_path(directory, filename) RESULT(path)
766 : CHARACTER(LEN=*), INTENT(IN) :: directory, filename
767 : CHARACTER(LEN=default_path_length) :: path
768 :
769 1 : IF (LEN_TRIM(directory) == 0 .OR. TRIM(directory) == ".") THEN
770 1 : path = TRIM(filename)
771 0 : ELSE IF (directory(LEN_TRIM(directory):LEN_TRIM(directory)) == "/") THEN
772 0 : path = TRIM(directory)//TRIM(filename)
773 : ELSE
774 0 : path = TRIM(directory)//"/"//TRIM(filename)
775 : END IF
776 :
777 1 : END FUNCTION tb_join_path
778 :
779 : ! **************************************************************************************************
780 : !> \brief Shell-quote a filename or executable path.
781 : !> \param text ...
782 : !> \return ...
783 : ! **************************************************************************************************
784 9 : FUNCTION tb_shell_quote(text) RESULT(quoted)
785 : CHARACTER(LEN=*), INTENT(IN) :: text
786 : CHARACTER(LEN=4*LEN_TRIM(text)+2) :: quoted
787 :
788 : INTEGER :: i, pos
789 :
790 9 : quoted = "'"
791 9 : pos = 2
792 150 : DO i = 1, LEN_TRIM(text)
793 150 : IF (text(i:i) == "'") THEN
794 : ! End the quoted part, escape the apostrophe, then resume quoting.
795 : ! ACHAR avoids dependence on the compiler's backslash-literal mode.
796 0 : quoted(pos:pos + 3) = "'"//ACHAR(92)//"''"
797 0 : pos = pos + 4
798 : ELSE
799 141 : quoted(pos:pos) = text(i:i)
800 141 : pos = pos + 1
801 : END IF
802 : END DO
803 9 : quoted(pos:pos) = "'"
804 :
805 9 : END FUNCTION tb_shell_quote
806 :
807 : ! **************************************************************************************************
808 : !> \brief Write current CP2K geometry as DFTB+ gen format for native tblite.
809 : !> \param qs_env ...
810 : !> \param filename ...
811 : ! **************************************************************************************************
812 1 : SUBROUTINE tb_write_reference_gen(qs_env, filename)
813 :
814 : TYPE(qs_environment_type), POINTER :: qs_env
815 : CHARACTER(LEN=*), INTENT(IN) :: filename
816 :
817 1 : CHARACTER(LEN=2), ALLOCATABLE, DIMENSION(:) :: symbols, unique_symbols
818 : INTEGER :: iatom, ikind, natom, nuniq, unit_nr
819 1 : INTEGER, ALLOCATABLE, DIMENSION(:) :: species
820 : INTEGER, DIMENSION(3) :: periodic
821 : LOGICAL :: found
822 : REAL(KIND=dp) :: to_angstrom
823 : TYPE(cell_type), POINTER :: cell
824 1 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
825 :
826 1 : NULLIFY (cell, particle_set)
827 1 : CALL get_qs_env(qs_env=qs_env, cell=cell, particle_set=particle_set)
828 :
829 1 : natom = SIZE(particle_set)
830 1 : to_angstrom = cp_unit_from_cp2k(1.0_dp, "angstrom")
831 5 : ALLOCATE (symbols(natom), unique_symbols(natom), species(natom))
832 1 : nuniq = 0
833 5 : DO iatom = 1, natom
834 4 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=symbols(iatom))
835 4 : found = .FALSE.
836 9 : DO ikind = 1, nuniq
837 9 : IF (TRIM(unique_symbols(ikind)) == TRIM(symbols(iatom))) THEN
838 : found = .TRUE.
839 : EXIT
840 : END IF
841 : END DO
842 4 : IF (.NOT. found) THEN
843 3 : nuniq = nuniq + 1
844 3 : unique_symbols(nuniq) = symbols(iatom)
845 3 : ikind = nuniq
846 : END IF
847 5 : species(iatom) = ikind
848 : END DO
849 :
850 : CALL open_file(file_name=filename, unit_number=unit_nr, file_status="REPLACE", &
851 1 : file_action="WRITE")
852 :
853 1 : CALL get_cell(cell=cell, periodic=periodic)
854 4 : IF (ANY(periodic == 1)) THEN
855 0 : WRITE (UNIT=unit_nr, FMT="(I0,1X,A)") natom, "S"
856 : ELSE
857 1 : WRITE (UNIT=unit_nr, FMT="(I0,1X,A)") natom, "C"
858 : END IF
859 4 : WRITE (UNIT=unit_nr, FMT="(*(A,1X))") (TRIM(unique_symbols(ikind)), ikind=1, nuniq)
860 5 : DO iatom = 1, natom
861 : WRITE (UNIT=unit_nr, FMT="(I0,1X,I0,3(1X,ES24.16))") &
862 17 : iatom, species(iatom), particle_set(iatom)%r(:)*to_angstrom
863 : END DO
864 4 : IF (ANY(periodic == 1)) THEN
865 0 : WRITE (UNIT=unit_nr, FMT="(3(1X,ES24.16))") 0.0_dp, 0.0_dp, 0.0_dp
866 0 : DO ikind = 1, 3
867 0 : WRITE (UNIT=unit_nr, FMT="(3(1X,ES24.16))") cell%hmat(:, ikind)*to_angstrom
868 : END DO
869 : END IF
870 1 : CALL close_file(unit_nr)
871 :
872 1 : DEALLOCATE (symbols, unique_symbols, species)
873 :
874 2 : END SUBROUTINE tb_write_reference_gen
875 :
876 : ! **************************************************************************************************
877 : !> \brief Read native tblite gradient file.
878 : !> \param filename ...
879 : !> \param natom ...
880 : !> \param energy ...
881 : !> \param gradient ...
882 : !> \param virial ...
883 : !> \param have_energy ...
884 : !> \param have_gradient ...
885 : !> \param have_virial ...
886 : ! **************************************************************************************************
887 0 : SUBROUTINE tb_read_reference_grad(filename, natom, energy, gradient, virial, &
888 : have_energy, have_gradient, have_virial)
889 :
890 : CHARACTER(LEN=*), INTENT(IN) :: filename
891 : INTEGER, INTENT(IN) :: natom
892 : REAL(KIND=dp), INTENT(OUT) :: energy
893 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: gradient, virial
894 : LOGICAL, INTENT(OUT) :: have_energy, have_gradient, have_virial
895 :
896 : CHARACTER(LEN=1024) :: line
897 : INTEGER :: ios, nread, unit_nr
898 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: values
899 :
900 0 : have_energy = .FALSE.
901 0 : have_gradient = .FALSE.
902 0 : have_virial = .FALSE.
903 0 : energy = 0.0_dp
904 0 : gradient = 0.0_dp
905 0 : virial = 0.0_dp
906 :
907 : ! Keep this diagnostic read nonfatal and confined to the requested output path.
908 : OPEN (NEWUNIT=unit_nr, FILE=TRIM(filename), STATUS="OLD", ACTION="READ", &
909 0 : FORM="FORMATTED", IOSTAT=ios)
910 0 : IF (ios /= 0) RETURN
911 :
912 : DO
913 0 : READ (UNIT=unit_nr, FMT="(A)", IOSTAT=ios) line
914 0 : IF (ios /= 0) EXIT
915 0 : IF (INDEX(line, "energy :real:0:") > 0) THEN
916 0 : READ (UNIT=unit_nr, FMT="(A)", IOSTAT=ios) line
917 0 : IF (ios == 0) THEN
918 0 : READ (line, *, IOSTAT=ios) energy
919 0 : have_energy = ios == 0
920 : END IF
921 0 : ELSE IF (INDEX(line, "gradient :real:2:3,") > 0) THEN
922 0 : ALLOCATE (values(3*natom))
923 0 : CALL tb_read_real_values(unit_nr, values, nread)
924 0 : IF (nread == 3*natom) THEN
925 0 : CALL tb_values_to_matrix(values, gradient)
926 0 : have_gradient = .TRUE.
927 : END IF
928 0 : DEALLOCATE (values)
929 0 : ELSE IF (INDEX(line, "virial :real:2:3,3") > 0) THEN
930 0 : ALLOCATE (values(9))
931 0 : CALL tb_read_real_values(unit_nr, values, nread)
932 0 : IF (nread == 9) THEN
933 0 : CALL tb_values_to_matrix(values, virial)
934 0 : have_virial = .TRUE.
935 : END IF
936 0 : DEALLOCATE (values)
937 : END IF
938 : END DO
939 0 : CALL close_file(unit_nr)
940 :
941 0 : END SUBROUTINE tb_read_reference_grad
942 :
943 : ! **************************************************************************************************
944 : !> \brief Read a fixed number of real values from following lines.
945 : !> \param unit_nr ...
946 : !> \param values ...
947 : !> \param nread ...
948 : ! **************************************************************************************************
949 0 : SUBROUTINE tb_read_real_values(unit_nr, values, nread)
950 :
951 : INTEGER, INTENT(IN) :: unit_nr
952 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: values
953 : INTEGER, INTENT(OUT) :: nread
954 :
955 : CHARACTER(LEN=1024) :: line
956 : INTEGER :: ios
957 :
958 0 : nread = 0
959 0 : DO WHILE (nread < SIZE(values))
960 0 : READ (UNIT=unit_nr, FMT="(A)", IOSTAT=ios) line
961 0 : IF (ios /= 0) EXIT
962 0 : CALL tb_parse_real_line(line, values, nread)
963 : END DO
964 :
965 0 : END SUBROUTINE tb_read_real_values
966 :
967 : ! **************************************************************************************************
968 : !> \brief Parse real values from one text line.
969 : !> \param line ...
970 : !> \param values ...
971 : !> \param nread ...
972 : ! **************************************************************************************************
973 0 : SUBROUTINE tb_parse_real_line(line, values, nread)
974 :
975 : CHARACTER(LEN=*), INTENT(IN) :: line
976 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: values
977 : INTEGER, INTENT(INOUT) :: nread
978 :
979 : CHARACTER(LEN=128) :: token
980 : INTEGER :: first, ios, last, pos
981 :
982 0 : pos = 1
983 0 : DO WHILE (pos <= LEN_TRIM(line) .AND. nread < SIZE(values))
984 0 : DO WHILE (pos <= LEN_TRIM(line) .AND. INDEX(" ,[]", line(pos:pos)) > 0)
985 0 : pos = pos + 1
986 : END DO
987 0 : IF (pos > LEN_TRIM(line)) EXIT
988 : first = pos
989 0 : DO WHILE (pos <= LEN_TRIM(line) .AND. INDEX(" ,[]", line(pos:pos)) == 0)
990 0 : pos = pos + 1
991 : END DO
992 0 : last = pos - 1
993 0 : token = line(first:last)
994 0 : READ (token, *, IOSTAT=ios) values(nread + 1)
995 0 : IF (ios == 0) nread = nread + 1
996 : END DO
997 :
998 0 : END SUBROUTINE tb_parse_real_line
999 :
1000 : ! **************************************************************************************************
1001 : !> \brief Convert flat native tblite values to CP2K atom-major matrix layout.
1002 : !> \param values ...
1003 : !> \param matrix ...
1004 : ! **************************************************************************************************
1005 0 : SUBROUTINE tb_values_to_matrix(values, matrix)
1006 :
1007 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: values
1008 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: matrix
1009 :
1010 : INTEGER :: i, j, n
1011 :
1012 0 : n = 0
1013 0 : DO j = 1, SIZE(matrix, 2)
1014 0 : DO i = 1, SIZE(matrix, 1)
1015 0 : n = n + 1
1016 0 : matrix(i, j) = values(n)
1017 : END DO
1018 : END DO
1019 :
1020 0 : END SUBROUTINE tb_values_to_matrix
1021 :
1022 : ! **************************************************************************************************
1023 : !> \brief Delete temporary files unless requested otherwise.
1024 : !> \param ref ...
1025 : !> \param gen_file ...
1026 : !> \param grad_file ...
1027 : !> \param json_file ...
1028 : !> \param log_file ...
1029 : !> \param post_processing_output_file ...
1030 : ! **************************************************************************************************
1031 0 : SUBROUTINE tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file, post_processing_output_file)
1032 :
1033 : TYPE(xtb_reference_cli_type), INTENT(IN) :: ref
1034 : CHARACTER(LEN=*), INTENT(IN) :: gen_file, grad_file, json_file, &
1035 : log_file, post_processing_output_file
1036 :
1037 0 : IF (ref%keep_files) RETURN
1038 0 : CALL tb_delete_file(gen_file)
1039 0 : CALL tb_delete_file(grad_file)
1040 0 : CALL tb_delete_file(json_file)
1041 0 : CALL tb_delete_file(log_file)
1042 0 : IF (LEN_TRIM(post_processing_output_file) > 0) THEN
1043 0 : CALL tb_delete_file(post_processing_output_file)
1044 : END IF
1045 :
1046 : END SUBROUTINE tb_reference_cleanup
1047 :
1048 : ! **************************************************************************************************
1049 : !> \brief Delete a file if it exists.
1050 : !> \param filename ...
1051 : ! **************************************************************************************************
1052 0 : SUBROUTINE tb_delete_file(filename)
1053 :
1054 : CHARACTER(LEN=*), INTENT(IN) :: filename
1055 :
1056 : INTEGER :: istat
1057 :
1058 : ! Cleanup is best effort, including files the reference CLI did not create.
1059 0 : istat = m_unlink(filename)
1060 : IF (istat /= 0) RETURN
1061 :
1062 : END SUBROUTINE tb_delete_file
1063 :
1064 : END MODULE tblite_reference_cli
|
