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 : MODULE qs_tddfpt2_fprint
9 : USE atomic_kind_types, ONLY: atomic_kind_type,&
10 : get_atomic_kind_set
11 : USE cp_control_types, ONLY: dft_control_type
12 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
13 : USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
14 : USE cp_fm_struct, ONLY: cp_fm_struct_type
15 : USE cp_fm_types, ONLY: cp_fm_create,&
16 : cp_fm_get_info,&
17 : cp_fm_release,&
18 : cp_fm_to_fm,&
19 : cp_fm_type
20 : USE cp_log_handling, ONLY: cp_get_default_logger,&
21 : cp_logger_get_default_io_unit,&
22 : cp_logger_type
23 : USE cp_output_handling, ONLY: cp_p_file,&
24 : cp_print_key_finished_output,&
25 : cp_print_key_should_output,&
26 : cp_print_key_unit_nr
27 : USE efield_utils, ONLY: calculate_ecore_efield
28 : USE exstates_types, ONLY: excited_energy_type
29 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
30 : section_vals_type,&
31 : section_vals_val_get
32 : USE kinds, ONLY: dp
33 : USE lri_environment_types, ONLY: lri_environment_type
34 : USE message_passing, ONLY: mp_para_env_type
35 : USE physcon, ONLY: evolt
36 : USE qs_core_energies, ONLY: calculate_ecore_overlap,&
37 : calculate_ecore_self
38 : USE qs_core_hamiltonian, ONLY: build_core_hamiltonian_matrix
39 : USE qs_energy_init, ONLY: qs_energies_init
40 : USE qs_environment_methods, ONLY: qs_env_rebuild_pw_env
41 : USE qs_environment_types, ONLY: get_qs_env,&
42 : qs_environment_type,&
43 : set_qs_env
44 : USE qs_external_potential, ONLY: external_c_potential,&
45 : external_e_potential
46 : USE qs_force_types, ONLY: allocate_qs_force,&
47 : deallocate_qs_force,&
48 : qs_force_type,&
49 : replicate_qs_force,&
50 : sum_qs_force,&
51 : total_qs_force,&
52 : zero_qs_force
53 : USE qs_kernel_types, ONLY: kernel_env_type
54 : USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
55 : USE qs_ks_types, ONLY: qs_ks_env_type,&
56 : set_ks_env
57 : USE qs_matrix_w, ONLY: compute_matrix_w
58 : USE qs_p_env_types, ONLY: p_env_release,&
59 : qs_p_env_type
60 : USE qs_scf, ONLY: scf
61 : USE qs_tddfpt2_forces, ONLY: tddfpt_forces_main
62 : USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
63 : USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos,&
64 : tddfpt_work_matrices
65 : USE response_solver, ONLY: response_equation,&
66 : response_force,&
67 : response_force_xtb
68 : USE ri_environment_methods, ONLY: build_ri_matrices
69 : USE xtb_ks_matrix, ONLY: build_xtb_ks_matrix
70 : USE xtb_matrices, ONLY: build_xtb_matrices
71 : #include "./base/base_uses.f90"
72 :
73 : IMPLICIT NONE
74 :
75 : PRIVATE
76 :
77 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_fprint'
78 :
79 : PUBLIC :: tddfpt_print_forces
80 :
81 : ! **************************************************************************************************
82 :
83 : CONTAINS
84 :
85 : ! **************************************************************************************************
86 : !> \brief Calculate and print forces of selected excited states.
87 : !> \param qs_env Information on Kinds and Particles
88 : !> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
89 : !> SIZE(evects,2) -- number of excited states to print)
90 : !> \param evals TDDFPT eigenvalues
91 : !> \param ostrength ...
92 : !> \param print_section ...
93 : !> \param gs_mos molecular orbitals optimised for the ground state
94 : !> \param kernel_env ...
95 : !> \param sub_env ...
96 : !> \param work_matrices ...
97 : !> \par History
98 : !> * 10.2022 created [JGH]
99 : ! **************************************************************************************************
100 554 : SUBROUTINE tddfpt_print_forces(qs_env, evects, evals, ostrength, print_section, &
101 : gs_mos, kernel_env, sub_env, work_matrices)
102 : TYPE(qs_environment_type), POINTER :: qs_env
103 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
104 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals, ostrength
105 : TYPE(section_vals_type), POINTER :: print_section
106 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
107 : POINTER :: gs_mos
108 : TYPE(kernel_env_type) :: kernel_env
109 : TYPE(tddfpt_subgroup_env_type) :: sub_env
110 : TYPE(tddfpt_work_matrices) :: work_matrices
111 :
112 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_print_forces'
113 : LOGICAL, PARAMETER :: debug_forces = .FALSE.
114 :
115 : INTEGER :: handle, iounit, is, ispin, istate, iw, &
116 : iwunit, kstate, natom, nspins, nstates
117 1108 : INTEGER, ALLOCATABLE, DIMENSION(:) :: alist, natom_of_kind, state_list
118 : REAL(KIND=dp) :: eener, threshold
119 554 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
120 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
121 : TYPE(cp_logger_type), POINTER :: logger
122 : TYPE(dft_control_type), POINTER :: dft_control
123 : TYPE(excited_energy_type), POINTER :: ex_env
124 : TYPE(mp_para_env_type), POINTER :: para_env
125 1108 : TYPE(qs_force_type), DIMENSION(:), POINTER :: gs_force, ks_force, td_force
126 : TYPE(qs_p_env_type) :: p_env
127 : TYPE(section_vals_type), POINTER :: force_section, tdlr_section
128 :
129 554 : CALL timeset(routineN, handle)
130 :
131 554 : logger => cp_get_default_logger()
132 554 : iounit = cp_logger_get_default_io_unit(logger)
133 :
134 554 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_section, &
135 : "FORCES"), cp_p_file)) THEN
136 :
137 4 : IF (sub_env%is_split) THEN
138 : CALL cp_abort(__LOCATION__, "Excited state forces not possible when states"// &
139 0 : " are distributed to different CPU pools.")
140 : END IF
141 :
142 4 : nspins = SIZE(evects, 1)
143 4 : nstates = SIZE(evects, 2)
144 4 : IF (iounit > 0) THEN
145 2 : WRITE (iounit, "(1X,A)") "", &
146 2 : "-------------------------------------------------------------------------------", &
147 2 : "- TDDFPT PROPERTIES: Nuclear Forces -", &
148 4 : "-------------------------------------------------------------------------------"
149 : END IF
150 4 : force_section => section_vals_get_subs_vals(print_section, "FORCES")
151 4 : CALL section_vals_val_get(force_section, "THRESHOLD", r_val=threshold)
152 4 : tdlr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%LINRES")
153 12 : ALLOCATE (state_list(nstates))
154 4 : CALL build_state_list(force_section, state_list)
155 : ! screen with oscillator strength
156 : ! Warning: if oscillator strength are not calculated they are set to zero and forces are
157 : ! only calculated if threshold is also 0
158 24 : DO istate = 1, nstates
159 24 : IF (ostrength(istate) < threshold) state_list(istate) = 0
160 : END DO
161 4 : IF (iounit > 0) THEN
162 2 : WRITE (iounit, "(1X,A,T61,E20.8)") " Screening threshold for oscillator strength ", threshold
163 4 : ALLOCATE (alist(nstates))
164 2 : is = 0
165 12 : DO istate = 1, nstates
166 12 : IF (state_list(istate) == 1) THEN
167 8 : is = is + 1
168 8 : alist(is) = istate
169 : END IF
170 : END DO
171 2 : WRITE (iounit, "(1X,A,T71,I10)") " List of states requested for force calculation ", is
172 2 : WRITE (iounit, "(16I5)") alist(1:is)
173 : END IF
174 :
175 : iwunit = cp_print_key_unit_nr(logger, force_section, "", &
176 4 : extension=".tdfrc", file_status='REPLACE')
177 :
178 : ! prepare force array
179 4 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, natom=natom)
180 4 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
181 4 : NULLIFY (td_force, gs_force)
182 4 : CALL allocate_qs_force(gs_force, natom_of_kind)
183 4 : CALL allocate_qs_force(td_force, natom_of_kind)
184 : ! ground state forces
185 4 : CALL gs_forces(qs_env, gs_force)
186 : ! Swap force arrays
187 4 : CALL get_qs_env(qs_env, force=ks_force)
188 4 : CALL set_qs_env(qs_env, force=td_force)
189 :
190 4 : CALL get_qs_env(qs_env, exstate_env=ex_env, para_env=para_env)
191 4 : kstate = ex_env%state
192 4 : eener = ex_env%evalue
193 : ! Start force loop over states
194 24 : DO istate = 1, nstates
195 24 : IF (state_list(istate) == 1) THEN
196 16 : IF (iounit > 0) THEN
197 8 : WRITE (iounit, "(1X,A,I3,T30,F10.5,A,T50,A,T71,F10.7)") " STATE NR. ", istate, &
198 16 : evals(istate)*evolt, " eV", "Oscillator strength:", ostrength(istate)
199 : END IF
200 16 : IF (iwunit > 0) THEN
201 8 : WRITE (iwunit, "(1X,A,I3,T30,F10.5,A,T50,A,T71,F10.7)") " # STATE NR. ", istate, &
202 16 : evals(istate)*evolt, " eV", "Oscillator strength:", ostrength(istate)
203 : END IF
204 16 : ex_env%state = istate
205 16 : ex_env%evalue = evals(istate)
206 16 : CALL cp_fm_release(ex_env%evect)
207 64 : ALLOCATE (ex_env%evect(nspins))
208 32 : DO ispin = 1, nspins
209 16 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct)
210 16 : CALL cp_fm_create(ex_env%evect(ispin), matrix_struct)
211 32 : CALL cp_fm_to_fm(evects(ispin, istate), ex_env%evect(ispin))
212 : END DO
213 : ! force array
214 16 : CALL zero_qs_force(td_force)
215 : !
216 16 : CALL tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, sub_env, work_matrices)
217 : !
218 16 : IF (debug_forces) THEN
219 : iw = iounit
220 : ELSE
221 : iw = -1
222 : END IF
223 16 : CALL response_equation(qs_env, p_env, ex_env%cpmos, iw, tdlr_section)
224 : !
225 16 : CALL get_qs_env(qs_env, dft_control=dft_control)
226 16 : IF (dft_control%qs_control%semi_empirical) THEN
227 0 : CPABORT("Not available")
228 16 : ELSEIF (dft_control%qs_control%dftb) THEN
229 0 : CPABORT("Not available")
230 16 : ELSEIF (dft_control%qs_control%xtb) THEN
231 0 : CALL response_force_xtb(qs_env, p_env, ex_env%matrix_hz, ex_env, debug=debug_forces)
232 : ELSE
233 : CALL response_force(qs_env=qs_env, vh_rspace=ex_env%vh_rspace, &
234 : vxc_rspace=ex_env%vxc_rspace, vtau_rspace=ex_env%vtau_rspace, &
235 : vadmm_rspace=ex_env%vadmm_rspace, matrix_hz=ex_env%matrix_hz, &
236 : matrix_pz=ex_env%matrix_px1, matrix_pz_admm=p_env%p1_admm, &
237 16 : matrix_wz=p_env%w1, p_env=p_env, ex_env=ex_env, debug=debug_forces)
238 : END IF
239 16 : CALL p_env_release(p_env)
240 : !
241 16 : CALL replicate_qs_force(td_force, para_env)
242 16 : CALL sum_qs_force(td_force, gs_force)
243 16 : CALL pforce(iwunit, td_force, atomic_kind_set, natom)
244 : !
245 : ELSE
246 4 : IF (iwunit > 0) THEN
247 2 : WRITE (iwunit, "(1X,A,I3,T30,F10.5,A,T50,A,T71,F10.7)") " # STATE NR. ", istate, &
248 4 : evals(istate)*evolt, " eV", "Oscillator strength:", ostrength(istate)
249 : END IF
250 : END IF
251 : END DO
252 4 : CALL set_qs_env(qs_env, force=ks_force)
253 4 : CALL deallocate_qs_force(gs_force)
254 4 : CALL deallocate_qs_force(td_force)
255 4 : DEALLOCATE (state_list)
256 :
257 4 : ex_env%state = kstate
258 4 : ex_env%evalue = eener
259 :
260 4 : CALL cp_print_key_finished_output(iwunit, logger, force_section, "")
261 :
262 8 : IF (iounit > 0) THEN
263 : WRITE (iounit, "(1X,A)") &
264 2 : "-------------------------------------------------------------------------------"
265 : END IF
266 :
267 : END IF
268 :
269 554 : CALL timestop(handle)
270 :
271 3878 : END SUBROUTINE tddfpt_print_forces
272 :
273 : ! **************************************************************************************************
274 : !> \brief Calculate the Quickstep forces. Adapted from qs_forces
275 : !> \param qs_env ...
276 : !> \param gs_force ...
277 : !> \date 11.2022
278 : !> \author JGH
279 : !> \version 1.0
280 : ! **************************************************************************************************
281 4 : SUBROUTINE gs_forces(qs_env, gs_force)
282 :
283 : TYPE(qs_environment_type), POINTER :: qs_env
284 : TYPE(qs_force_type), DIMENSION(:), POINTER :: gs_force
285 :
286 : CHARACTER(len=*), PARAMETER :: routineN = 'gs_forces'
287 :
288 : INTEGER :: handle
289 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_w_b, matrix_w_kp
290 : TYPE(dft_control_type), POINTER :: dft_control
291 : TYPE(lri_environment_type), POINTER :: lri_env
292 : TYPE(mp_para_env_type), POINTER :: para_env
293 4 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
294 : TYPE(qs_ks_env_type), POINTER :: ks_env
295 :
296 4 : CALL timeset(routineN, handle)
297 :
298 : ! Swap froce arrays
299 4 : CALL get_qs_env(qs_env, force=force)
300 4 : CALL zero_qs_force(gs_force)
301 4 : CALL set_qs_env(qs_env, force=gs_force)
302 :
303 : ! Check for incompatible options
304 4 : CALL get_qs_env(qs_env, dft_control=dft_control)
305 4 : CPASSERT(.NOT. dft_control%qs_control%cdft)
306 4 : CPASSERT(.NOT. qs_env%run_rtp)
307 4 : CPASSERT(.NOT. dft_control%qs_control%mulliken_restraint)
308 4 : CPASSERT(.NOT. dft_control%dft_plus_u)
309 4 : CPASSERT(.NOT. qs_env%energy_correction)
310 4 : IF (ASSOCIATED(qs_env%mp2_env)) THEN
311 0 : CPABORT("TDDFPT| MP2 not available")
312 : END IF
313 :
314 : ! Save current W matrix
315 4 : CALL get_qs_env(qs_env=qs_env, matrix_w_kp=matrix_w_b)
316 4 : NULLIFY (matrix_w_kp)
317 4 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
318 4 : CALL set_ks_env(ks_env, matrix_w_kp=matrix_w_kp)
319 :
320 : ! recalculate energy with forces
321 4 : CPASSERT(.NOT. dft_control%qs_control%do_ls_scf)
322 4 : CPASSERT(.NOT. dft_control%qs_control%do_almo_scf)
323 4 : CALL qs_env_rebuild_pw_env(qs_env)
324 4 : CALL qs_energies_init(qs_env, .TRUE.)
325 4 : CALL scf(qs_env)
326 4 : CALL compute_matrix_w(qs_env, .TRUE.)
327 :
328 4 : NULLIFY (para_env)
329 4 : CALL get_qs_env(qs_env, para_env=para_env)
330 4 : IF (dft_control%qs_control%semi_empirical) THEN
331 0 : CPABORT("TDDFPT| SE not available")
332 4 : ELSEIF (dft_control%qs_control%dftb) THEN
333 0 : CPABORT("TDDFPT| DFTB not available")
334 4 : ELSEIF (dft_control%qs_control%xtb) THEN
335 0 : CALL build_xtb_matrices(qs_env=qs_env, calculate_forces=.TRUE.)
336 0 : CALL build_xtb_ks_matrix(qs_env, calculate_forces=.TRUE., just_energy=.FALSE.)
337 : ELSE
338 4 : CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.TRUE.)
339 4 : CALL calculate_ecore_self(qs_env)
340 4 : CALL calculate_ecore_overlap(qs_env, para_env, calculate_forces=.TRUE.)
341 4 : CALL calculate_ecore_efield(qs_env, calculate_forces=.TRUE.)
342 4 : CALL external_e_potential(qs_env)
343 4 : IF (.NOT. dft_control%qs_control%gapw) THEN
344 4 : CALL external_c_potential(qs_env, calculate_forces=.TRUE.)
345 : END IF
346 4 : IF (dft_control%qs_control%rigpw) THEN
347 0 : CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
348 0 : CALL build_ri_matrices(lri_env, qs_env, calculate_forces=.TRUE.)
349 : END IF
350 4 : CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.TRUE., just_energy=.FALSE.)
351 : END IF
352 :
353 : ! replicate forces (get current pointer)
354 4 : NULLIFY (gs_force)
355 4 : CALL get_qs_env(qs_env=qs_env, force=gs_force)
356 4 : CALL replicate_qs_force(gs_force, para_env)
357 : ! Swap back force array
358 4 : CALL set_qs_env(qs_env=qs_env, force=force)
359 :
360 : ! deallocate W Matrix and bring back saved one
361 4 : CALL get_qs_env(qs_env=qs_env, matrix_w_kp=matrix_w_kp)
362 4 : CALL dbcsr_deallocate_matrix_set(matrix_w_kp)
363 4 : NULLIFY (matrix_w_kp)
364 4 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
365 4 : CALL set_ks_env(ks_env, matrix_w_kp=matrix_w_b)
366 :
367 4 : CALL timestop(handle)
368 :
369 4 : END SUBROUTINE gs_forces
370 :
371 : ! **************************************************************************************************
372 : !> \brief building a state list
373 : !> \param section input section
374 : !> \param state_list ...
375 : ! **************************************************************************************************
376 4 : SUBROUTINE build_state_list(section, state_list)
377 :
378 : TYPE(section_vals_type), POINTER :: section
379 : INTEGER, DIMENSION(:), INTENT(OUT) :: state_list
380 :
381 : INTEGER :: i, is, k, n_rep, nstate
382 4 : INTEGER, DIMENSION(:), POINTER :: indexes
383 : LOGICAL :: explicit
384 :
385 4 : nstate = SIZE(state_list)
386 :
387 4 : CALL section_vals_val_get(section, "LIST", explicit=explicit, n_rep_val=n_rep)
388 4 : IF (explicit) THEN
389 24 : state_list = 0
390 8 : DO i = 1, n_rep
391 4 : CALL section_vals_val_get(section, "LIST", i_rep_val=i, i_vals=indexes)
392 26 : DO is = 1, SIZE(indexes)
393 18 : k = indexes(is)
394 18 : IF (k <= 0 .OR. k > nstate) THEN
395 0 : CALL cp_warn(__LOCATION__, "State List contains invalid state.")
396 0 : CPABORT("TDDFPT Print Forces: Invalid State")
397 : END IF
398 22 : state_list(k) = 1
399 : END DO
400 : END DO
401 : ELSE
402 0 : state_list = 1
403 : END IF
404 :
405 4 : END SUBROUTINE build_state_list
406 :
407 : ! **************************************************************************************************
408 : !> \brief ...
409 : !> \param iwunit ...
410 : !> \param td_force ...
411 : !> \param atomic_kind_set ...
412 : !> \param natom ...
413 : ! **************************************************************************************************
414 16 : SUBROUTINE pforce(iwunit, td_force, atomic_kind_set, natom)
415 : INTEGER, INTENT(IN) :: iwunit
416 : TYPE(qs_force_type), DIMENSION(:), POINTER :: td_force
417 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
418 : INTEGER, INTENT(IN) :: natom
419 :
420 : INTEGER :: iatom
421 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: force
422 :
423 48 : ALLOCATE (force(3, natom))
424 16 : CALL total_qs_force(force, td_force, atomic_kind_set)
425 :
426 16 : IF (iwunit > 0) THEN
427 8 : WRITE (iwunit, *) natom
428 8 : WRITE (iwunit, *)
429 32 : DO iatom = 1, natom
430 104 : WRITE (iwunit, "(3F24.14)") - force(1:3, iatom)
431 : END DO
432 : END IF
433 16 : DEALLOCATE (force)
434 :
435 16 : END SUBROUTINE pforce
436 :
437 : END MODULE qs_tddfpt2_fprint
|
