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 : MODULE qs_tddfpt2_methods
9 : USE admm_methods, ONLY: admm_fit_mo_coeffs
10 : USE admm_types, ONLY: admm_type,&
11 : get_admm_env
12 : USE atomic_kind_types, ONLY: atomic_kind_type
13 : USE bibliography, ONLY: Grimme2013,&
14 : Grimme2016,&
15 : Hernandez2025,&
16 : Iannuzzi2005,&
17 : cite_reference
18 : USE cell_types, ONLY: cell_type
19 : USE cp_blacs_env, ONLY: cp_blacs_env_type
20 : USE cp_control_types, ONLY: dft_control_type,&
21 : rixs_control_type,&
22 : tddfpt2_control_type
23 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
24 : USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
25 : USE cp_fm_pool_types, ONLY: fm_pool_create_fm
26 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
27 : cp_fm_struct_release,&
28 : cp_fm_struct_type
29 : USE cp_fm_types, ONLY: cp_fm_create,&
30 : cp_fm_get_info,&
31 : cp_fm_release,&
32 : cp_fm_to_fm,&
33 : cp_fm_to_fm_submat,&
34 : cp_fm_type
35 : USE cp_log_handling, ONLY: cp_get_default_logger,&
36 : cp_logger_get_default_io_unit,&
37 : cp_logger_type
38 : USE cp_output_handling, ONLY: cp_add_iter_level,&
39 : cp_iterate,&
40 : cp_print_key_finished_output,&
41 : cp_print_key_unit_nr,&
42 : cp_rm_iter_level
43 : USE exstates_types, ONLY: excited_energy_type
44 : USE header, ONLY: tddfpt_header,&
45 : tddfpt_soc_header
46 : USE hfx_admm_utils, ONLY: aux_admm_init
47 : USE hfx_types, ONLY: compare_hfx_sections,&
48 : hfx_create
49 : USE input_constants, ONLY: &
50 : do_admm_aux_exch_func_none, do_admm_basis_projection, do_admm_exch_scaling_none, &
51 : do_admm_purify_none, do_potential_truncated, no_sf_tddfpt, tddfpt_dipole_velocity, &
52 : tddfpt_kernel_full, tddfpt_kernel_none, tddfpt_kernel_stda, tddfpt_sf_col, tddfpt_sf_noncol
53 : USE input_section_types, ONLY: section_vals_get,&
54 : section_vals_get_subs_vals,&
55 : section_vals_type,&
56 : section_vals_val_get,&
57 : section_vals_val_set
58 : USE kinds, ONLY: dp
59 : USE lri_environment_methods, ONLY: lri_print_stat
60 : USE lri_environment_types, ONLY: lri_density_release,&
61 : lri_env_release
62 : USE machine, ONLY: m_flush
63 : USE message_passing, ONLY: mp_para_env_type
64 : USE min_basis_set, ONLY: create_minbas_set
65 : USE particle_types, ONLY: particle_type
66 : USE physcon, ONLY: evolt
67 : USE qs_environment_types, ONLY: get_qs_env,&
68 : qs_environment_type
69 : USE qs_kernel_methods, ONLY: create_fxc_kernel,&
70 : create_kernel_env
71 : USE qs_kernel_types, ONLY: full_kernel_env_type,&
72 : kernel_env_type,&
73 : release_kernel_env
74 : USE qs_kind_types, ONLY: qs_kind_type
75 : USE qs_mo_types, ONLY: mo_set_type
76 : USE qs_scf_methods, ONLY: eigensolver
77 : USE qs_scf_types, ONLY: qs_scf_env_type
78 : USE qs_tddfpt2_assign, ONLY: assign_state
79 : USE qs_tddfpt2_densities, ONLY: tddfpt_construct_aux_fit_density,&
80 : tddfpt_construct_ground_state_orb_density
81 : USE qs_tddfpt2_eigensolver, ONLY: tddfpt_davidson_solver,&
82 : tddfpt_orthogonalize_psi1_psi0,&
83 : tddfpt_orthonormalize_psi1_psi1
84 : USE qs_tddfpt2_forces, ONLY: tddfpt_forces_main
85 : USE qs_tddfpt2_fprint, ONLY: tddfpt_print_forces
86 : USE qs_tddfpt2_lri_utils, ONLY: tddfpt2_lri_init
87 : USE qs_tddfpt2_properties, ONLY: tddfpt_dipole_operator,&
88 : tddfpt_print_excitation_analysis,&
89 : tddfpt_print_exciton_descriptors,&
90 : tddfpt_print_nto_analysis,&
91 : tddfpt_print_summary
92 : USE qs_tddfpt2_restart, ONLY: tddfpt_read_restart,&
93 : tddfpt_write_newtonx_output,&
94 : tddfpt_write_restart
95 : USE qs_tddfpt2_smearing_methods, ONLY: tddfpt_smeared_occupation
96 : USE qs_tddfpt2_soc, ONLY: tddfpt_soc
97 : USE qs_tddfpt2_stda_types, ONLY: allocate_stda_env,&
98 : deallocate_stda_env,&
99 : stda_env_type,&
100 : stda_init_param
101 : USE qs_tddfpt2_stda_utils, ONLY: get_lowdin_mo_coefficients,&
102 : stda_init_matrices
103 : USE qs_tddfpt2_subgroups, ONLY: tddfpt_sub_env_init,&
104 : tddfpt_sub_env_release,&
105 : tddfpt_subgroup_env_type
106 : USE qs_tddfpt2_types, ONLY: hfxsr_create_work_matrices,&
107 : stda_create_work_matrices,&
108 : tddfpt_create_work_matrices,&
109 : tddfpt_ground_state_mos,&
110 : tddfpt_release_work_matrices,&
111 : tddfpt_work_matrices
112 : USE qs_tddfpt2_utils, ONLY: tddfpt_guess_vectors,&
113 : tddfpt_init_mos,&
114 : tddfpt_oecorr,&
115 : tddfpt_release_ground_state_mos
116 : USE rixs_types, ONLY: rixs_env_type,&
117 : tddfpt2_valence_type
118 : USE string_utilities, ONLY: integer_to_string
119 : USE xc_write_output, ONLY: xc_write
120 : #include "./base/base_uses.f90"
121 :
122 : IMPLICIT NONE
123 :
124 : PRIVATE
125 :
126 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_methods'
127 :
128 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
129 : ! number of first derivative components (3: d/dx, d/dy, d/dz)
130 : INTEGER, PARAMETER, PRIVATE :: nderivs = 3
131 : INTEGER, PARAMETER, PRIVATE :: maxspins = 2
132 :
133 : PUBLIC :: tddfpt, tddfpt_energies, tddfpt_input
134 :
135 : ! **************************************************************************************************
136 :
137 : CONTAINS
138 :
139 : ! **************************************************************************************************
140 : !> \brief Perform TDDFPT calculation. If calc_forces then it also builds the response vector for the
141 : !> Z-vector method and calculates some contributions to the force
142 : !> \param qs_env Quickstep environment
143 : !> \param calc_forces ...
144 : !> \param rixs_env ...
145 : !> \par History
146 : !> * 05.2016 created [Sergey Chulkov]
147 : !> * 06.2016 refactored to be used with Davidson eigensolver [Sergey Chulkov]
148 : !> * 03.2017 cleaned and refactored [Sergey Chulkov]
149 : !> \note Based on the subroutines tddfpt_env_init(), and tddfpt_env_deallocate().
150 : ! **************************************************************************************************
151 1342 : SUBROUTINE tddfpt(qs_env, calc_forces, rixs_env)
152 : TYPE(qs_environment_type), POINTER :: qs_env
153 : LOGICAL, INTENT(IN) :: calc_forces
154 : TYPE(rixs_env_type), OPTIONAL, POINTER :: rixs_env
155 :
156 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt'
157 :
158 : INTEGER :: handle, ispin, istate, log_unit, mult, &
159 : my_state, nao, nocc, nspins, &
160 : nstate_max, nstates, nvirt, old_state
161 : INTEGER, DIMENSION(maxspins) :: nactive
162 : LOGICAL :: do_admm, do_exck, do_hfx, do_hfxlr, &
163 : do_hfxsr, do_rixs, do_sf, do_soc, &
164 : lmult_tmp, state_change
165 : REAL(kind=dp) :: gsmin, gsval, xsval
166 1342 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals, ostrength
167 1342 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
168 : TYPE(cell_type), POINTER :: cell
169 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
170 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
171 1342 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: my_active, my_mos
172 1342 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: dipole_op_mos_occ, evects, S_evects
173 : TYPE(cp_logger_type), POINTER :: logger
174 1342 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_oep, matrix_s, &
175 1342 : matrix_s_aux_fit, &
176 1342 : matrix_s_aux_fit_vs_orb
177 : TYPE(dft_control_type), POINTER :: dft_control
178 : TYPE(excited_energy_type), POINTER :: ex_env
179 : TYPE(full_kernel_env_type), TARGET :: full_kernel_env, kernel_env_admm_aux
180 : TYPE(kernel_env_type) :: kernel_env
181 1342 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
182 : TYPE(mp_para_env_type), POINTER :: para_env
183 1342 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
184 1342 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
185 : TYPE(qs_scf_env_type), POINTER :: scf_env
186 : TYPE(rixs_control_type), POINTER :: rixs_control
187 : TYPE(section_vals_type), POINTER :: hfxsr_section, kernel_section, &
188 : lri_section, soc_section, &
189 : tddfpt_print_section, tddfpt_section, &
190 : xc_section
191 : TYPE(stda_env_type), TARGET :: stda_kernel
192 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
193 : TYPE(tddfpt2_valence_type), POINTER :: valence_state
194 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
195 1342 : POINTER :: gs_mos
196 1342 : TYPE(tddfpt_subgroup_env_type) :: sub_env
197 1342 : TYPE(tddfpt_work_matrices) :: work_matrices
198 :
199 1342 : CALL timeset(routineN, handle)
200 :
201 1342 : NULLIFY (logger)
202 1342 : logger => cp_get_default_logger()
203 :
204 1342 : NULLIFY (tddfpt_section, tddfpt_control)
205 :
206 : CALL get_qs_env(qs_env, &
207 : blacs_env=blacs_env, &
208 : cell=cell, &
209 : dft_control=dft_control, &
210 : matrix_ks=matrix_ks, &
211 : matrix_s=matrix_s, &
212 : mos=mos, &
213 : scf_env=scf_env, &
214 1342 : do_rixs=do_rixs)
215 :
216 1342 : IF (do_rixs) THEN
217 14 : tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%RIXS%TDDFPT")
218 14 : NULLIFY (rixs_control, valence_state)
219 14 : rixs_control => dft_control%rixs_control
220 14 : tddfpt_control => rixs_control%tddfpt2_control
221 14 : valence_state => rixs_env%valence_state
222 : ELSE
223 1328 : tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT")
224 1328 : tddfpt_control => dft_control%tddfpt2_control
225 : END IF
226 :
227 : ! input section print/xc
228 : CALL tddfpt_input(qs_env, tddfpt_section, tddfpt_control, do_hfx, do_admm, do_exck, &
229 : do_hfxsr, do_hfxlr, xc_section, tddfpt_print_section, &
230 1342 : lri_section, hfxsr_section)
231 :
232 : log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, "PROGRAM_BANNER", &
233 1342 : extension=".tddfptLog")
234 :
235 1342 : tddfpt_control%do_hfx = do_hfx
236 1342 : tddfpt_control%do_admm = do_admm
237 1342 : tddfpt_control%do_hfxsr = do_hfxsr
238 1342 : tddfpt_control%hfxsr_primbas = 0
239 1342 : tddfpt_control%hfxsr_re_int = .TRUE.
240 1342 : tddfpt_control%do_hfxlr = do_hfxlr
241 1342 : tddfpt_control%do_exck = do_exck
242 1342 : do_sf = tddfpt_control%spinflip /= no_sf_tddfpt
243 1342 : IF (do_sf) CALL cite_reference(Hernandez2025)
244 1342 : IF (tddfpt_control%do_hfxlr) THEN
245 6 : kernel_section => section_vals_get_subs_vals(tddfpt_section, "XC%HFX_KERNEL%HFXLR")
246 6 : CALL section_vals_val_get(kernel_section, "RCUT", r_val=tddfpt_control%hfxlr_rcut)
247 6 : CALL section_vals_val_get(kernel_section, "SCALE", r_val=tddfpt_control%hfxlr_scale)
248 : END IF
249 :
250 1342 : soc_section => section_vals_get_subs_vals(tddfpt_section, "SOC")
251 1342 : CALL section_vals_get(soc_section, explicit=do_soc)
252 :
253 1342 : IF (do_soc) THEN
254 : ! start with multiplicity that is not specified in input
255 : ! so that excited-state gradient is for multiplicity given in input
256 10 : lmult_tmp = tddfpt_control%rks_triplets
257 10 : tddfpt_control%rks_triplets = .NOT. (tddfpt_control%rks_triplets)
258 : END IF
259 :
260 1342 : CALL cite_reference(Iannuzzi2005)
261 1342 : IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
262 408 : CALL cite_reference(Grimme2013)
263 408 : CALL cite_reference(Grimme2016)
264 : END IF
265 :
266 1342 : CALL tddfpt_header(log_unit)
267 1342 : CALL kernel_info(log_unit, dft_control, tddfpt_control, xc_section)
268 :
269 : ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals
270 1342 : NULLIFY (gs_mos)
271 1342 : CALL tddfpt_init_mos(qs_env, gs_mos, log_unit)
272 :
273 : ! obtain smeared occupation numbers
274 1342 : IF (tddfpt_control%do_smearing) THEN
275 2 : CALL tddfpt_smeared_occupation(qs_env, gs_mos, log_unit)
276 : END IF
277 :
278 : ! obtain corrected KS-matrix
279 1342 : CALL tddfpt_oecorr(qs_env, gs_mos, matrix_ks_oep)
280 :
281 1342 : IF ((tddfpt_control%do_lrigpw) .AND. &
282 : (tddfpt_control%kernel /= tddfpt_kernel_full)) THEN
283 0 : CALL cp_abort(__LOCATION__, "LRI only implemented for full kernel")
284 : END IF
285 :
286 1342 : IF (ASSOCIATED(matrix_ks_oep)) matrix_ks => matrix_ks_oep
287 :
288 : ! determine active orbitals
289 : ! default is all occupied MOs
290 1342 : CALL init_res_method(qs_env, gs_mos, tddfpt_control, tddfpt_section, log_unit)
291 :
292 : ! components of the dipole operator
293 : CALL tddfpt_dipole_operator(dipole_op_mos_occ, &
294 : tddfpt_control, &
295 : gs_mos, &
296 1342 : qs_env)
297 :
298 1342 : nspins = SIZE(gs_mos)
299 : ! multiplicity of molecular system
300 1342 : IF (nspins > 1) THEN
301 152 : mult = ABS(SIZE(gs_mos(1)%evals_occ) - SIZE(gs_mos(2)%evals_occ)) + 1
302 152 : IF (mult > 2) &
303 22 : CALL cp_warn(__LOCATION__, "There is a convergence issue for multiplicity >= 3")
304 : ELSE
305 1190 : IF (tddfpt_control%rks_triplets) THEN
306 204 : mult = 3
307 : ELSE
308 986 : mult = 1
309 : END IF
310 : END IF
311 :
312 : ! split mpi communicator
313 8356 : ALLOCATE (my_mos(nspins), my_active(nspins))
314 2836 : DO ispin = 1, nspins
315 1494 : my_mos(ispin) = gs_mos(ispin)%mos_occ
316 2836 : my_active(ispin) = gs_mos(ispin)%mos_active
317 : END DO
318 : CALL tddfpt_sub_env_init(sub_env, qs_env, &
319 : mos_occ=my_mos(:), mos_active=my_active(:), &
320 1342 : kernel=tddfpt_control%kernel)
321 1342 : DEALLOCATE (my_mos, my_active)
322 :
323 1342 : IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
324 : ! create environment for Full Kernel
325 810 : IF (dft_control%qs_control%xtb) THEN
326 0 : CPABORT("TDDFPT: xTB only works with sTDA Kernel")
327 : END IF
328 :
329 810 : IF (tddfpt_control%do_hfxsr) THEN
330 4 : kernel_section => section_vals_get_subs_vals(tddfpt_section, "XC%HFX_KERNEL")
331 : CALL section_vals_val_get(kernel_section, "HFXSR_PRIMBAS", &
332 4 : i_val=tddfpt_control%hfxsr_primbas)
333 : ! basis set
334 : CALL create_minbas_set(qs_env, log_unit, basis_type="TDA_HFX", &
335 4 : primitive=tddfpt_control%hfxsr_primbas)
336 : ! admm control
337 20 : ALLOCATE (full_kernel_env%admm_control)
338 4 : full_kernel_env%admm_control%purification_method = do_admm_purify_none
339 : full_kernel_env%admm_control%method = do_admm_basis_projection
340 : full_kernel_env%admm_control%scaling_model = do_admm_exch_scaling_none
341 4 : full_kernel_env%admm_control%aux_exch_func = do_admm_aux_exch_func_none
342 : ! hfx section
343 4 : full_kernel_env%hfxsr_section => hfxsr_section
344 : !
345 : CALL aux_admm_init(qs_env, mos, full_kernel_env%admm_env, &
346 4 : full_kernel_env%admm_control, "TDA_HFX")
347 : CALL get_admm_env(full_kernel_env%admm_env, mos_aux_fit=mos_aux_fit, &
348 : matrix_s_aux_fit=matrix_s_aux_fit, &
349 4 : matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb)
350 : CALL admm_fit_mo_coeffs(full_kernel_env%admm_env, matrix_s_aux_fit, &
351 4 : matrix_s_aux_fit_vs_orb, mos, mos_aux_fit, .TRUE.)
352 : ! x_data
353 : CALL get_qs_env(qs_env, cell=cell, atomic_kind_set=atomic_kind_set, &
354 : qs_kind_set=qs_kind_set, particle_set=particle_set, &
355 4 : para_env=para_env)
356 : CALL hfx_create(full_kernel_env%x_data, para_env, hfxsr_section, atomic_kind_set, &
357 4 : qs_kind_set, particle_set, dft_control, cell, orb_basis="TDA_HFX")
358 : END IF
359 :
360 : ! allocate pools and work matrices
361 810 : nstates = tddfpt_control%nstates
362 : !! Too many states can lead to Problems
363 : !! You should be warned if there are more states
364 : !! than occ-virt Combinations!!
365 810 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=nocc)
366 810 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
367 788 : CALL cp_fm_get_info(gs_mos(1)%mos_virt, ncol_global=nvirt)
368 : ELSE
369 22 : CALL cp_fm_get_info(gs_mos(2)%mos_virt, ncol_global=nvirt)
370 : END IF
371 810 : nstate_max = nocc*nvirt
372 810 : IF (nstates > nstate_max) THEN
373 0 : CPWARN("NUMBER OF EXCITED STATES COULD LEAD TO PROBLEMS!")
374 0 : CPWARN("Experimental: CHANGED NSTATES TO ITS MAXIMUM VALUE!")
375 0 : nstates = nstate_max
376 0 : tddfpt_control%nstates = nstate_max
377 : END IF
378 : CALL tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, &
379 810 : do_hfx, do_admm, do_hfxlr, do_exck, do_sf, qs_env, sub_env)
380 :
381 : ! create full_kernel and admm_kernel within tddfpt_energies
382 810 : kernel_env%full_kernel => full_kernel_env
383 810 : kernel_env%admm_kernel => kernel_env_admm_aux
384 810 : NULLIFY (kernel_env%stda_kernel)
385 810 : IF (do_hfxsr) THEN
386 : ! work matrices for SR HFX
387 4 : CALL hfxsr_create_work_matrices(work_matrices, qs_env, full_kernel_env%admm_env)
388 : END IF
389 810 : IF (do_hfxlr) THEN
390 : ! calculate S_half and Lowdin MO coefficients
391 6 : CALL get_lowdin_mo_coefficients(qs_env, sub_env, work_matrices)
392 : END IF
393 532 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
394 : ! setup for kernel_stda outside tddfpt_energies
395 408 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
396 1224 : nactive = tddfpt_control%nactive
397 408 : CALL allocate_stda_env(qs_env, stda_kernel, nao, nactive)
398 : ! sTDA parameters
399 408 : CALL stda_init_param(qs_env, stda_kernel, tddfpt_control%stda_control)
400 : ! allocate pools and work matrices
401 408 : nstates = tddfpt_control%nstates
402 408 : CALL stda_create_work_matrices(work_matrices, gs_mos, nstates, qs_env, sub_env)
403 : !
404 : CALL stda_init_matrices(qs_env, stda_kernel, sub_env, &
405 408 : work_matrices, tddfpt_control)
406 : !
407 408 : kernel_env%stda_kernel => stda_kernel
408 408 : NULLIFY (kernel_env%full_kernel)
409 408 : NULLIFY (kernel_env%admm_kernel)
410 124 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
411 : ! allocate pools and work matrices
412 124 : nstates = tddfpt_control%nstates
413 124 : CALL stda_create_work_matrices(work_matrices, gs_mos, nstates, qs_env, sub_env)
414 124 : NULLIFY (kernel_env%full_kernel)
415 124 : NULLIFY (kernel_env%admm_kernel)
416 124 : NULLIFY (kernel_env%stda_kernel)
417 : END IF
418 :
419 1342 : IF (do_sf) THEN
420 : ! only alpha -> beta excitations are considered in spin-flip TDDFT
421 246 : ALLOCATE (evects(1, nstates))
422 : ELSE
423 12172 : ALLOCATE (evects(nspins, nstates))
424 : END IF
425 4026 : ALLOCATE (evals(nstates))
426 12440 : ALLOCATE (S_evects(SIZE(evects, 1), nstates))
427 :
428 4676 : DO istate = 1, nstates
429 8414 : DO ispin = 1, SIZE(evects, 1)
430 : CALL fm_pool_create_fm( &
431 : work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
432 7072 : S_evects(ispin, istate))
433 : END DO
434 : END DO
435 :
436 1342 : IF (.NOT. do_soc) THEN
437 : ! compute tddfpt excitation energies of multiplicity mult
438 : CALL tddfpt_energies(qs_env, nstates, nspins, work_matrices, &
439 : tddfpt_control, logger, tddfpt_print_section, evects, evals, &
440 : gs_mos, tddfpt_section, S_evects, matrix_s, kernel_env, matrix_ks, &
441 : sub_env, ostrength, dipole_op_mos_occ, mult, xc_section, full_kernel_env, &
442 1332 : kernel_env_admm_aux)
443 : ELSE
444 : CALL tddfpt_soc_energies(qs_env, nstates, work_matrices, &
445 : tddfpt_control, logger, tddfpt_print_section, &
446 : evects, evals, ostrength, &
447 : gs_mos, tddfpt_section, S_evects, matrix_s, kernel_env, matrix_ks, &
448 : sub_env, dipole_op_mos_occ, lmult_tmp, xc_section, full_kernel_env, &
449 10 : kernel_env_admm_aux)
450 : END IF
451 :
452 : !print forces for selected states
453 1342 : IF (calc_forces) THEN
454 : CALL tddfpt_print_forces(qs_env, evects, evals, ostrength, &
455 : tddfpt_print_section, gs_mos, &
456 638 : kernel_env, sub_env, work_matrices)
457 : END IF
458 :
459 : ! excited state potential energy surface
460 1342 : IF (qs_env%excited_state) THEN
461 1132 : IF (sub_env%is_split) THEN
462 : CALL cp_abort(__LOCATION__, &
463 : "Excited state forces not possible when states"// &
464 0 : " are distributed to different CPU pools.")
465 : END IF
466 : ! for gradients unshifted KS matrix
467 1132 : IF (ASSOCIATED(matrix_ks_oep)) CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
468 1132 : CALL get_qs_env(qs_env, exstate_env=ex_env)
469 1132 : state_change = .FALSE.
470 1132 : IF (ex_env%state > 0) THEN
471 1124 : my_state = ex_env%state
472 8 : ELSEIF (ex_env%state < 0) THEN
473 : ! state following
474 32 : ALLOCATE (my_mos(nspins))
475 16 : DO ispin = 1, nspins
476 16 : my_mos(ispin) = gs_mos(ispin)%mos_occ
477 : END DO
478 8 : my_state = ABS(ex_env%state)
479 8 : CALL assign_state(qs_env, matrix_s, evects, my_mos, ex_env%wfn_history, my_state)
480 8 : DEALLOCATE (my_mos)
481 8 : IF (my_state /= ABS(ex_env%state)) THEN
482 0 : state_change = .TRUE.
483 0 : old_state = ABS(ex_env%state)
484 : END IF
485 8 : ex_env%state = -my_state
486 : ELSE
487 : CALL cp_warn(__LOCATION__, &
488 0 : "Active excited state not assigned. Use the first state.")
489 0 : my_state = 1
490 : END IF
491 1132 : CPASSERT(my_state > 0)
492 1132 : IF (my_state > nstates) THEN
493 : CALL cp_warn(__LOCATION__, &
494 0 : "There were not enough excited states calculated.")
495 0 : CPABORT("excited state potential energy surface")
496 : END IF
497 : !
498 : ! energy
499 1132 : ex_env%evalue = evals(my_state)
500 : ! excitation vector
501 1132 : CALL cp_fm_release(ex_env%evect)
502 4624 : ALLOCATE (ex_env%evect(SIZE(evects, 1)))
503 2360 : DO ispin = 1, SIZE(evects, 1)
504 : CALL cp_fm_get_info(matrix=evects(ispin, 1), &
505 1228 : matrix_struct=matrix_struct)
506 1228 : CALL cp_fm_create(ex_env%evect(ispin), matrix_struct)
507 2360 : CALL cp_fm_to_fm(evects(ispin, my_state), ex_env%evect(ispin))
508 : END DO
509 :
510 1132 : IF (log_unit > 0) THEN
511 566 : gsval = ex_env%wfn_history%gsval
512 566 : gsmin = ex_env%wfn_history%gsmin
513 566 : xsval = ex_env%wfn_history%xsval
514 566 : WRITE (log_unit, "(1X,A,T40,F10.6,A,T62,F10.6,A)") "Ground state orbital alignment:", &
515 1132 : gsmin, "[MinVal]", gsval, "[Average]"
516 566 : WRITE (log_unit, "(1X,A,T71,F10.6)") "Excitation vector alignment:", xsval
517 566 : IF (state_change) THEN
518 : WRITE (log_unit, "(1X,A,I5,T60,A14,T76,I5)") &
519 0 : "Target state has been changed from state ", &
520 0 : old_state, " to new state ", my_state
521 : END IF
522 566 : WRITE (log_unit, "(1X,A,I4,A,F12.5,A)") "Calculate properties for state:", &
523 1132 : my_state, " with excitation energy ", ex_env%evalue*evolt, " eV"
524 : END IF
525 :
526 : ! Calculate response vector
527 1132 : IF (calc_forces) THEN
528 : CALL tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, &
529 636 : sub_env, work_matrices)
530 : END IF
531 : END IF
532 :
533 : ! share evals, evects and mo_coefs with rixs
534 1342 : IF (do_rixs) THEN
535 : ! copy evals
536 14 : valence_state%nstates = nstates
537 42 : ALLOCATE (valence_state%evals(SIZE(evals)))
538 62 : valence_state%evals(:) = evals(:)
539 :
540 172 : ALLOCATE (valence_state%evects(nspins, nstates))
541 60 : ALLOCATE (valence_state%mos_active(nspins))
542 32 : DO ispin = 1, nspins
543 : ! copy evects
544 86 : DO istate = 1, nstates
545 : CALL cp_fm_get_info(matrix=evects(ispin, istate), &
546 68 : matrix_struct=matrix_struct)
547 68 : CALL cp_fm_create(valence_state%evects(ispin, istate), matrix_struct)
548 86 : CALL cp_fm_to_fm(evects(ispin, istate), valence_state%evects(ispin, istate))
549 : END DO
550 : ! copy mos_occ
551 : CALL cp_fm_get_info(matrix=gs_mos(ispin)%mos_active, &
552 18 : matrix_struct=matrix_struct)
553 18 : CALL cp_fm_create(valence_state%mos_active(ispin), matrix_struct)
554 32 : CALL cp_fm_to_fm(gs_mos(ispin)%mos_active, valence_state%mos_active(ispin))
555 : END DO
556 : END IF
557 :
558 : ! clean up
559 1342 : CALL cp_fm_release(evects)
560 1342 : CALL cp_fm_release(S_evects)
561 :
562 : CALL cp_print_key_finished_output(log_unit, &
563 : logger, &
564 : tddfpt_print_section, &
565 1342 : "PROGRAM_BANNER")
566 :
567 1342 : DEALLOCATE (evals, ostrength)
568 :
569 1342 : IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
570 810 : IF (do_admm) CALL release_kernel_env(kernel_env%admm_kernel)
571 810 : IF (tddfpt_control%do_lrigpw) THEN
572 10 : CALL lri_env_release(kernel_env%full_kernel%lri_env)
573 10 : DEALLOCATE (kernel_env%full_kernel%lri_env)
574 10 : CALL lri_density_release(kernel_env%full_kernel%lri_density)
575 10 : DEALLOCATE (kernel_env%full_kernel%lri_density)
576 : END IF
577 810 : CALL release_kernel_env(kernel_env%full_kernel)
578 532 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
579 408 : CALL deallocate_stda_env(stda_kernel)
580 124 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
581 : !
582 : ELSE
583 0 : CPABORT('Unknown kernel type')
584 : END IF
585 1342 : CALL tddfpt_release_work_matrices(work_matrices, sub_env)
586 1342 : CALL tddfpt_sub_env_release(sub_env)
587 :
588 1342 : CALL cp_fm_release(dipole_op_mos_occ)
589 :
590 2836 : DO ispin = nspins, 1, -1
591 2836 : CALL tddfpt_release_ground_state_mos(gs_mos(ispin))
592 : END DO
593 1342 : DEALLOCATE (gs_mos)
594 :
595 1342 : IF (ASSOCIATED(matrix_ks_oep)) &
596 30 : CALL dbcsr_deallocate_matrix_set(matrix_ks_oep)
597 :
598 1342 : CALL timestop(handle)
599 :
600 10736 : END SUBROUTINE tddfpt
601 :
602 : ! **************************************************************************************************
603 : !> \brief TDDFPT input
604 : !> \param qs_env Quickstep environment
605 : !> \param tddfpt_section ...
606 : !> \param tddfpt_control ...
607 : !> \param do_hfx ...
608 : !> \param do_admm ...
609 : !> \param do_exck ...
610 : !> \param do_hfxsr ...
611 : !> \param do_hfxlr ...
612 : !> \param xc_section ...
613 : !> \param tddfpt_print_section ...
614 : !> \param lri_section ...
615 : !> \param hfxsr_section ...
616 : ! **************************************************************************************************
617 1342 : SUBROUTINE tddfpt_input(qs_env, tddfpt_section, tddfpt_control, do_hfx, do_admm, do_exck, &
618 : do_hfxsr, do_hfxlr, xc_section, tddfpt_print_section, lri_section, &
619 : hfxsr_section)
620 : TYPE(qs_environment_type), POINTER :: qs_env
621 : TYPE(section_vals_type), POINTER :: tddfpt_section
622 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
623 : LOGICAL, INTENT(INOUT) :: do_hfx, do_admm, do_exck, do_hfxsr, &
624 : do_hfxlr
625 : TYPE(section_vals_type), POINTER :: xc_section, tddfpt_print_section, &
626 : lri_section, hfxsr_section
627 :
628 : CHARACTER(len=20) :: nstates_str
629 : LOGICAL :: exar, exf, exgcp, exhf, exhfxk, exk, &
630 : explicit_root, expot, exvdw, exwfn, &
631 : found, same_hfx
632 : REAL(kind=dp) :: C_hf
633 : TYPE(dft_control_type), POINTER :: dft_control
634 : TYPE(section_vals_type), POINTER :: hfx_section, hfx_section_gs, input, &
635 : xc_root, xc_sub
636 :
637 1342 : NULLIFY (dft_control, input)
638 1342 : CALL get_qs_env(qs_env, dft_control=dft_control, input=input)
639 :
640 : ! no k-points
641 1342 : IF (dft_control%nimages > 1) CPABORT("k-points not implemented for TDDFPT")
642 :
643 1342 : IF (tddfpt_control%nstates <= 0) THEN
644 0 : CALL integer_to_string(tddfpt_control%nstates, nstates_str)
645 : CALL cp_warn(__LOCATION__, "TDDFPT calculation was requested for "// &
646 0 : TRIM(nstates_str)//" excited states: nothing to do.")
647 0 : RETURN
648 : END IF
649 :
650 1342 : NULLIFY (tddfpt_print_section)
651 1342 : tddfpt_print_section => section_vals_get_subs_vals(tddfpt_section, "PRINT")
652 :
653 1342 : IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
654 810 : NULLIFY (xc_root)
655 810 : xc_root => section_vals_get_subs_vals(tddfpt_section, "XC")
656 810 : CALL section_vals_get(xc_root, explicit=explicit_root)
657 810 : NULLIFY (xc_section)
658 810 : IF (explicit_root) THEN
659 : ! No ADIABATIC_RESCALING option possible
660 476 : NULLIFY (xc_sub)
661 476 : xc_sub => section_vals_get_subs_vals(xc_root, "ADIABATIC_RESCALING")
662 476 : CALL section_vals_get(xc_sub, explicit=exar)
663 476 : IF (exar) THEN
664 0 : CALL cp_warn(__LOCATION__, "TDDFPT Kernel with ADIABATIC_RESCALING not possible.")
665 0 : CPABORT("TDDFPT Input")
666 : END IF
667 : ! No GCP_POTENTIAL option possible
668 476 : NULLIFY (xc_sub)
669 476 : xc_sub => section_vals_get_subs_vals(xc_root, "GCP_POTENTIAL")
670 476 : CALL section_vals_get(xc_sub, explicit=exgcp)
671 476 : IF (exgcp) THEN
672 0 : CALL cp_warn(__LOCATION__, "TDDFPT Kernel with GCP_POTENTIAL not possible.")
673 0 : CPABORT("TDDFPT Input")
674 : END IF
675 : ! No VDW_POTENTIAL option possible
676 476 : NULLIFY (xc_sub)
677 476 : xc_sub => section_vals_get_subs_vals(xc_root, "VDW_POTENTIAL")
678 476 : CALL section_vals_get(xc_sub, explicit=exvdw)
679 476 : IF (exvdw) THEN
680 0 : CALL cp_warn(__LOCATION__, "TDDFPT Kernel with VDW_POTENTIAL not possible.")
681 0 : CPABORT("TDDFPT Input")
682 : END IF
683 : ! No WF_CORRELATION option possible
684 476 : NULLIFY (xc_sub)
685 476 : xc_sub => section_vals_get_subs_vals(xc_root, "WF_CORRELATION")
686 476 : CALL section_vals_get(xc_sub, explicit=exwfn)
687 476 : IF (exwfn) THEN
688 0 : CALL cp_warn(__LOCATION__, "TDDFPT Kernel with WF_CORRELATION not possible.")
689 0 : CPABORT("TDDFPT Input")
690 : END IF
691 : ! No XC_POTENTIAL option possible
692 476 : NULLIFY (xc_sub)
693 476 : xc_sub => section_vals_get_subs_vals(xc_root, "XC_POTENTIAL")
694 476 : CALL section_vals_get(xc_sub, explicit=expot)
695 476 : IF (expot) THEN
696 0 : CALL cp_warn(__LOCATION__, "TDDFPT Kernel with XC_POTENTIAL not possible.")
697 0 : CPABORT("TDDFPT Input")
698 : END IF
699 : !
700 476 : NULLIFY (xc_sub)
701 476 : xc_sub => section_vals_get_subs_vals(xc_root, "XC_FUNCTIONAL")
702 476 : CALL section_vals_get(xc_sub, explicit=exf)
703 476 : NULLIFY (xc_sub)
704 476 : xc_sub => section_vals_get_subs_vals(xc_root, "XC_KERNEL")
705 476 : CALL section_vals_get(xc_sub, explicit=exk)
706 476 : IF ((exf .AND. exk) .OR. .NOT. (exf .OR. exk)) THEN
707 0 : CALL cp_warn(__LOCATION__, "TDDFPT Kernel needs XC_FUNCTIONAL or XC_KERNEL section.")
708 0 : CPABORT("TDDFPT Input")
709 : END IF
710 476 : NULLIFY (xc_sub)
711 476 : xc_sub => section_vals_get_subs_vals(xc_root, "HF")
712 476 : CALL section_vals_get(xc_sub, explicit=exhf)
713 476 : NULLIFY (xc_sub)
714 476 : xc_sub => section_vals_get_subs_vals(xc_root, "HFX_KERNEL")
715 476 : CALL section_vals_get(xc_sub, explicit=exhfxk)
716 : !
717 476 : xc_section => xc_root
718 476 : hfx_section => section_vals_get_subs_vals(xc_section, "HF")
719 476 : CALL section_vals_get(hfx_section, explicit=do_hfx)
720 476 : IF (do_hfx) THEN
721 24 : CALL section_vals_val_get(hfx_section, "FRACTION", r_val=C_hf)
722 24 : do_hfx = (C_hf /= 0.0_dp)
723 : END IF
724 : !TDDFPT only works if the kernel has the same HF section as the DFT%XC one
725 476 : IF (do_hfx) THEN
726 24 : hfx_section_gs => section_vals_get_subs_vals(input, "DFT%XC%HF")
727 24 : CALL compare_hfx_sections(hfx_section, hfx_section_gs, same_hfx)
728 24 : IF (.NOT. same_hfx) THEN
729 0 : CPABORT("TDDFPT Kernel must use the same HF section as DFT%XC or no HF at all.")
730 : END IF
731 : END IF
732 :
733 476 : do_admm = do_hfx .AND. dft_control%do_admm
734 476 : IF (do_admm) THEN
735 : ! 'admm_env%xc_section_primary' and 'admm_env%xc_section_aux' need to be redefined
736 : CALL cp_abort(__LOCATION__, &
737 : "ADMM is not implemented for a TDDFT kernel XC-functional which is different from "// &
738 0 : "the one used for the ground-state calculation. A ground-state 'admm_env' cannot be reused.")
739 : END IF
740 : ! SET HFX_KERNEL and/or XC_KERNEL
741 476 : IF (exk) THEN
742 12 : do_exck = .TRUE.
743 : ELSE
744 464 : do_exck = .FALSE.
745 : END IF
746 476 : IF (exhfxk) THEN
747 6 : xc_sub => section_vals_get_subs_vals(xc_root, "HFX_KERNEL")
748 6 : CALL section_vals_val_get(xc_sub, "DO_HFXSR", l_val=do_hfxsr)
749 6 : xc_sub => section_vals_get_subs_vals(xc_root, "HFX_KERNEL%HFXLR")
750 6 : CALL section_vals_get(xc_sub, explicit=do_hfxlr)
751 : ELSE
752 470 : do_hfxsr = .FALSE.
753 470 : do_hfxlr = .FALSE.
754 : END IF
755 : ELSE
756 334 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
757 334 : hfx_section => section_vals_get_subs_vals(xc_section, "HF")
758 334 : CALL section_vals_get(hfx_section, explicit=do_hfx)
759 334 : IF (do_hfx) THEN
760 266 : CALL section_vals_val_get(hfx_section, "FRACTION", r_val=C_hf)
761 266 : do_hfx = (C_hf /= 0.0_dp)
762 : END IF
763 334 : do_admm = do_hfx .AND. dft_control%do_admm
764 334 : do_exck = .FALSE.
765 334 : do_hfxsr = .FALSE.
766 334 : do_hfxlr = .FALSE.
767 : END IF
768 : ELSE
769 532 : do_hfx = .FALSE.
770 532 : do_admm = .FALSE.
771 532 : do_exck = .FALSE.
772 532 : do_hfxsr = .FALSE.
773 532 : do_hfxlr = .FALSE.
774 : END IF
775 :
776 : ! reset rks_triplets if UKS is in use
777 1342 : IF (tddfpt_control%rks_triplets .AND. dft_control%nspins > 1) THEN
778 8 : tddfpt_control%rks_triplets = .FALSE.
779 8 : CALL cp_warn(__LOCATION__, "Keyword RKS_TRIPLETS has been ignored for spin-polarised calculations")
780 : END IF
781 :
782 : ! lri input
783 1342 : IF (tddfpt_control%do_lrigpw) THEN
784 10 : lri_section => section_vals_get_subs_vals(tddfpt_section, "LRIGPW")
785 : END IF
786 :
787 : ! set defaults for short range HFX
788 1342 : NULLIFY (hfxsr_section)
789 1342 : IF (do_hfxsr) THEN
790 4 : hfxsr_section => section_vals_get_subs_vals(tddfpt_section, "XC%HFX_KERNEL%HF")
791 4 : CALL section_vals_get(hfxsr_section, explicit=found)
792 4 : IF (.NOT. found) THEN
793 0 : CPABORT("HFXSR option needs &HF section defined")
794 : END IF
795 4 : CALL section_vals_val_get(hfxsr_section, "INTERACTION_POTENTIAL%POTENTIAL_TYPE", explicit=found)
796 4 : IF (.NOT. found) THEN
797 : CALL section_vals_val_set(hfxsr_section, "INTERACTION_POTENTIAL%POTENTIAL_TYPE", &
798 4 : i_val=do_potential_truncated)
799 : END IF
800 4 : CALL section_vals_val_get(hfxsr_section, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", explicit=found)
801 4 : IF (.NOT. found) THEN
802 4 : CALL section_vals_val_set(hfxsr_section, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", r_val=7.5589_dp)
803 : END IF
804 4 : CALL section_vals_val_get(hfxsr_section, "RI%_SECTION_PARAMETERS_", l_val=found)
805 4 : IF (found) THEN
806 0 : CALL cp_abort(__LOCATION__, "Short range TDA kernel with RI not possible")
807 : END IF
808 : END IF
809 :
810 : END SUBROUTINE tddfpt_input
811 :
812 : ! **************************************************************************************************
813 : !> \brief ...
814 : !> \param log_unit ...
815 : !> \param dft_control ...
816 : !> \param tddfpt_control ...
817 : !> \param xc_section ...
818 : ! **************************************************************************************************
819 1342 : SUBROUTINE kernel_info(log_unit, dft_control, tddfpt_control, xc_section)
820 : INTEGER, INTENT(IN) :: log_unit
821 : TYPE(dft_control_type), POINTER :: dft_control
822 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
823 : TYPE(section_vals_type), POINTER :: xc_section
824 :
825 : CHARACTER(LEN=4) :: ktype
826 : LOGICAL :: lsd
827 :
828 1342 : lsd = (dft_control%nspins > 1)
829 1342 : IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
830 810 : ktype = "FULL"
831 810 : IF (log_unit > 0) THEN
832 405 : WRITE (log_unit, "(T2,A,T77,A4)") "KERNEL|", TRIM(ktype)
833 405 : CALL xc_write(log_unit, xc_section, lsd)
834 405 : IF (tddfpt_control%do_hfx) THEN
835 145 : IF (tddfpt_control%do_admm) THEN
836 84 : WRITE (log_unit, "(T2,A,T62,A19)") "KERNEL|", "ADMM Exact Exchange"
837 84 : IF (tddfpt_control%admm_xc_correction) THEN
838 64 : WRITE (log_unit, "(T2,A,T60,A21)") "KERNEL|", "Apply ADMM Kernel XC Correction"
839 : END IF
840 84 : IF (tddfpt_control%admm_symm) THEN
841 84 : WRITE (log_unit, "(T2,A,T60,A21)") "KERNEL|", "Symmetric ADMM Kernel"
842 : END IF
843 : ELSE
844 61 : WRITE (log_unit, "(T2,A,T67,A14)") "KERNEL|", "Exact Exchange"
845 : END IF
846 : END IF
847 405 : IF (tddfpt_control%do_hfxsr) THEN
848 2 : WRITE (log_unit, "(T2,A,T43,A38)") "KERNEL|", "Short range HFX approximation"
849 : END IF
850 405 : IF (tddfpt_control%do_hfxlr) THEN
851 3 : WRITE (log_unit, "(T2,A,T43,A38)") "KERNEL|", "Long range HFX approximation"
852 : END IF
853 405 : IF (tddfpt_control%do_lrigpw) THEN
854 5 : WRITE (log_unit, "(T2,A,T42,A39)") "KERNEL|", "LRI approximation of transition density"
855 : END IF
856 : END IF
857 532 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
858 408 : ktype = "sTDA"
859 408 : IF (log_unit > 0) THEN
860 204 : WRITE (log_unit, "(T2,A,T77,A4)") "KERNEL|", TRIM(ktype)
861 204 : IF (tddfpt_control%stda_control%do_ewald) THEN
862 47 : WRITE (log_unit, "(T2,A,T78,A3)") "KERNEL| Coulomb term uses Ewald summation"
863 : ELSE
864 157 : WRITE (log_unit, "(T2,A,T78,A3)") "KERNEL| Coulomb term uses direct summation (MIC)"
865 : END IF
866 204 : IF (tddfpt_control%stda_control%do_exchange) THEN
867 188 : WRITE (log_unit, "(T2,A,T78,A3)") "KERNEL| Exact exchange term", "YES"
868 188 : WRITE (log_unit, "(T2,A,T71,F10.3)") "KERNEL| Short range HFX fraction:", &
869 376 : tddfpt_control%stda_control%hfx_fraction
870 : ELSE
871 16 : WRITE (log_unit, "(T2,A,T79,A2)") "KERNEL| Exact exchange term", "NO"
872 : END IF
873 204 : WRITE (log_unit, "(T2,A,T66,E15.3)") "KERNEL| Transition density filter", &
874 408 : tddfpt_control%stda_control%eps_td_filter
875 : END IF
876 124 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
877 124 : ktype = "NONE"
878 124 : IF (log_unit > 0) THEN
879 62 : WRITE (log_unit, "(T2,A,T77,A4)") "KERNEL|", TRIM(ktype)
880 : END IF
881 : ELSE
882 : !CPABORT("Unknown kernel")
883 : END IF
884 : !
885 1342 : IF (log_unit > 0) THEN
886 671 : IF (tddfpt_control%rks_triplets) THEN
887 102 : WRITE (log_unit, "(T2,A,T74,A7)") "KERNEL| Spin symmetry of excitations", "Triplet"
888 569 : ELSE IF (lsd) THEN
889 : ! Spin-conserving excitations where requested
890 76 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
891 65 : WRITE (log_unit, "(T2,A,T69,A12)") "KERNEL| Spin symmetry of excitations", "Unrestricted"
892 : ! Spin-flip excitations with collinear exchange-correlation kernel requested
893 11 : ELSE IF (tddfpt_control%spinflip == tddfpt_sf_col) THEN
894 5 : WRITE (log_unit, "(T2,A,T72,A9)") "KERNEL| Spin flip", "Collinear"
895 : ! Spin-flip excitations with noncollinear exchange-correlation kernel requested
896 6 : ELSE IF (tddfpt_control%spinflip == tddfpt_sf_noncol) THEN
897 6 : WRITE (log_unit, "(T2,A,T69,A12)") "KERNEL| Spin flip", "Noncollinear"
898 : END IF
899 : ELSE
900 493 : WRITE (log_unit, "(T2,A,T74,A7)") "KERNEL| Spin symmetry of excitations", "Singlet"
901 : END IF
902 671 : WRITE (log_unit, "(T2,A,T73,I8)") "TDDFPT| Number of states calculated", tddfpt_control%nstates
903 671 : WRITE (log_unit, "(T2,A,T73,I8)") "TDDFPT| Number of Davidson iterations", tddfpt_control%niters
904 671 : WRITE (log_unit, "(T2,A,T66,E15.3)") "TDDFPT| Davidson iteration convergence", tddfpt_control%conv
905 671 : WRITE (log_unit, "(T2,A,T73,I8)") "TDDFPT| Max. number of Krylov space vectors", tddfpt_control%nkvs
906 : END IF
907 :
908 1342 : END SUBROUTINE kernel_info
909 :
910 : ! **************************************************************************************************
911 : !> \brief The energy calculation has been moved to its own subroutine
912 : !> \param qs_env ...
913 : !> \param nstates ...
914 : !> \param nspins ...
915 : !> \param work_matrices ...
916 : !> \param tddfpt_control ...
917 : !> \param logger ...
918 : !> \param tddfpt_print_section ...
919 : !> \param evects ...
920 : !> \param evals ...
921 : !> \param gs_mos ...
922 : !> \param tddfpt_section ...
923 : !> \param S_evects ...
924 : !> \param matrix_s ...
925 : !> \param kernel_env ...
926 : !> \param matrix_ks ...
927 : !> \param sub_env ...
928 : !> \param ostrength ...
929 : !> \param dipole_op_mos_occ ...
930 : !> \param mult ...
931 : !> \param xc_section ...
932 : !> \param full_kernel_env ...
933 : !> \param kernel_env_admm_aux ...
934 : ! **************************************************************************************************
935 1352 : SUBROUTINE tddfpt_energies(qs_env, nstates, nspins, work_matrices, &
936 : tddfpt_control, logger, tddfpt_print_section, evects, evals, &
937 : gs_mos, tddfpt_section, S_evects, matrix_s, kernel_env, matrix_ks, &
938 : sub_env, ostrength, dipole_op_mos_occ, mult, xc_section, full_kernel_env, &
939 : kernel_env_admm_aux)
940 :
941 : TYPE(qs_environment_type), POINTER :: qs_env
942 : INTEGER :: nstates, nspins
943 : TYPE(tddfpt_work_matrices) :: work_matrices
944 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
945 : TYPE(cp_logger_type), POINTER :: logger
946 : TYPE(section_vals_type), POINTER :: tddfpt_print_section
947 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: evects
948 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
949 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
950 : POINTER :: gs_mos
951 : TYPE(section_vals_type), POINTER :: tddfpt_section
952 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: S_evects
953 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
954 : TYPE(kernel_env_type) :: kernel_env
955 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
956 : TYPE(tddfpt_subgroup_env_type) :: sub_env
957 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ostrength
958 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: dipole_op_mos_occ
959 : INTEGER :: mult
960 : TYPE(section_vals_type), POINTER :: xc_section
961 : TYPE(full_kernel_env_type), TARGET :: full_kernel_env, kernel_env_admm_aux
962 :
963 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_energies'
964 :
965 : CHARACTER(len=20) :: nstates_str
966 : INTEGER :: energy_unit, handle, iter, log_unit, &
967 : niters, nocc, nstate_max, &
968 : nstates_read, nvirt
969 : LOGICAL :: do_admm, do_exck, do_soc, explicit
970 : REAL(kind=dp) :: conv
971 : TYPE(admm_type), POINTER :: admm_env
972 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
973 1352 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_oep
974 : TYPE(section_vals_type), POINTER :: lri_section, namd_print_section, &
975 : soc_section
976 :
977 1352 : CALL timeset(routineN, handle)
978 :
979 1352 : NULLIFY (admm_env, matrix_ks_oep)
980 1352 : do_admm = tddfpt_control%do_admm
981 1352 : IF (do_admm) CALL get_qs_env(qs_env, admm_env=admm_env)
982 :
983 : ! setup for full_kernel and admm_kernel within tddfpt_energies due to dependence on multiplicity
984 1352 : IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
985 :
986 : CALL tddfpt_construct_ground_state_orb_density( &
987 : rho_orb_struct=work_matrices%rho_orb_struct_sub, &
988 : rho_xc_struct=work_matrices%rho_xc_struct_sub, &
989 : is_rks_triplets=tddfpt_control%rks_triplets, &
990 : qs_env=qs_env, sub_env=sub_env, &
991 820 : wfm_rho_orb=work_matrices%rho_ao_orb_fm_sub)
992 :
993 820 : IF (do_admm) THEN
994 : ! Full kernel with ADMM
995 168 : IF (tddfpt_control%admm_xc_correction) THEN
996 : CALL create_kernel_env(kernel_env=full_kernel_env, &
997 : rho_struct_sub=work_matrices%rho_orb_struct_sub, &
998 : xc_section=admm_env%xc_section_primary, &
999 : is_rks_triplets=tddfpt_control%rks_triplets, &
1000 128 : sub_env=sub_env)
1001 : ELSE
1002 : CALL create_kernel_env(kernel_env=full_kernel_env, &
1003 : rho_struct_sub=work_matrices%rho_orb_struct_sub, &
1004 : xc_section=xc_section, &
1005 : is_rks_triplets=tddfpt_control%rks_triplets, &
1006 40 : sub_env=sub_env)
1007 : END IF
1008 :
1009 : CALL tddfpt_construct_aux_fit_density( &
1010 : rho_orb_struct=work_matrices%rho_orb_struct_sub, &
1011 : rho_aux_fit_struct=work_matrices%rho_aux_fit_struct_sub, &
1012 : local_rho_set=sub_env%local_rho_set_admm, &
1013 : qs_env=qs_env, sub_env=sub_env, &
1014 : wfm_rho_orb=work_matrices%rho_ao_orb_fm_sub, &
1015 : wfm_rho_aux_fit=work_matrices%rho_ao_aux_fit_fm_sub, &
1016 168 : wfm_aux_orb=work_matrices%wfm_aux_orb_sub)
1017 :
1018 : CALL create_kernel_env(kernel_env=kernel_env_admm_aux, &
1019 : rho_struct_sub=work_matrices%rho_aux_fit_struct_sub, &
1020 : xc_section=admm_env%xc_section_aux, &
1021 : is_rks_triplets=tddfpt_control%rks_triplets, &
1022 168 : sub_env=sub_env)
1023 168 : kernel_env%full_kernel => full_kernel_env
1024 168 : kernel_env%admm_kernel => kernel_env_admm_aux
1025 : ELSE
1026 : ! Full kernel
1027 : CALL create_kernel_env(kernel_env=full_kernel_env, &
1028 : rho_struct_sub=work_matrices%rho_orb_struct_sub, &
1029 : xc_section=xc_section, &
1030 : is_rks_triplets=tddfpt_control%rks_triplets, &
1031 652 : sub_env=sub_env)
1032 652 : kernel_env%full_kernel => full_kernel_env
1033 652 : NULLIFY (kernel_env%admm_kernel)
1034 : END IF
1035 : ! Fxc from kernel definition
1036 820 : do_exck = tddfpt_control%do_exck
1037 820 : kernel_env%full_kernel%do_exck = do_exck
1038 : ! initilize xc kernel
1039 820 : IF (do_exck) THEN
1040 : CALL create_fxc_kernel(work_matrices%rho_orb_struct_sub, work_matrices%fxc_rspace_sub, &
1041 12 : xc_section, tddfpt_control%rks_triplets, sub_env, qs_env)
1042 : END IF
1043 : END IF
1044 :
1045 : ! lri input
1046 1352 : IF (tddfpt_control%do_lrigpw) THEN
1047 10 : lri_section => section_vals_get_subs_vals(tddfpt_section, "LRIGPW")
1048 : CALL tddfpt2_lri_init(qs_env, kernel_env, lri_section, &
1049 10 : tddfpt_print_section)
1050 : END IF
1051 :
1052 : !! Too many states can lead to Problems
1053 : !! You should be warned if there are more states
1054 : !! than occ-virt Combinations!!
1055 1352 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=nocc)
1056 1352 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1057 1330 : CALL cp_fm_get_info(gs_mos(1)%mos_virt, ncol_global=nvirt)
1058 : ELSE
1059 22 : CALL cp_fm_get_info(gs_mos(2)%mos_virt, ncol_global=nvirt)
1060 : END IF
1061 1352 : nstate_max = nocc*nvirt
1062 1352 : IF ((SIZE(gs_mos, 1) == 2) .AND. (tddfpt_control%spinflip == no_sf_tddfpt)) THEN
1063 130 : CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nocc)
1064 130 : CALL cp_fm_get_info(gs_mos(2)%mos_virt, ncol_global=nvirt)
1065 130 : nstate_max = nocc*nvirt + nstate_max
1066 : END IF
1067 1352 : IF (nstates > nstate_max) THEN
1068 0 : CPWARN("NUMBER OF EXCITED STATES COULD LEAD TO PROBLEMS!")
1069 0 : CPWARN("Experimental: CHANGED NSTATES TO ITS MAXIMUM VALUE!")
1070 0 : nstates = nstate_max
1071 : END IF
1072 :
1073 1352 : soc_section => section_vals_get_subs_vals(tddfpt_section, "SOC")
1074 1352 : CALL section_vals_get(soc_section, explicit=do_soc)
1075 :
1076 : ! reuse Ritz vectors from the previous calculation if available
1077 1352 : IF (tddfpt_control%is_restart .AND. .NOT. do_soc) THEN
1078 6 : CALL get_qs_env(qs_env, blacs_env=blacs_env)
1079 :
1080 : nstates_read = tddfpt_read_restart( &
1081 : evects=evects, &
1082 : evals=evals, &
1083 : gs_mos=gs_mos, &
1084 : logger=logger, &
1085 : tddfpt_section=tddfpt_section, &
1086 : tddfpt_print_section=tddfpt_print_section, &
1087 : fm_pool_ao_mo_active=work_matrices%fm_pool_ao_mo_active, &
1088 6 : blacs_env_global=blacs_env)
1089 : ELSE
1090 : nstates_read = 0
1091 : END IF
1092 :
1093 : ! build the list of missed singly excited states and sort them in ascending order
1094 : ! according to their excitation energies
1095 : log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, &
1096 1352 : "GUESS_VECTORS", extension=".tddfptLog")
1097 : CALL tddfpt_guess_vectors(evects=evects, evals=evals, &
1098 : gs_mos=gs_mos, log_unit=log_unit, tddfpt_control=tddfpt_control, &
1099 : fm_pool_ao_mo_active=work_matrices%fm_pool_ao_mo_active, &
1100 1352 : qs_env=qs_env, nspins=nspins)
1101 : CALL cp_print_key_finished_output(log_unit, logger, &
1102 1352 : tddfpt_print_section, "GUESS_VECTORS")
1103 :
1104 : CALL tddfpt_orthogonalize_psi1_psi0(evects, work_matrices%S_C0_C0T, qs_env, &
1105 1352 : gs_mos, evals, tddfpt_control, work_matrices%S_C0)
1106 1352 : CALL tddfpt_orthonormalize_psi1_psi1(evects, nstates, S_evects, matrix_s(1)%matrix)
1107 :
1108 1352 : niters = tddfpt_control%niters
1109 1352 : IF (niters > 0) THEN
1110 : log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, &
1111 1352 : "ITERATION_INFO", extension=".tddfptLog")
1112 : energy_unit = cp_print_key_unit_nr(logger, &
1113 : tddfpt_print_section, &
1114 : "DETAILED_ENERGY", &
1115 1352 : extension=".tddfptLog")
1116 :
1117 1352 : IF (log_unit > 0) THEN
1118 676 : WRITE (log_unit, "(1X,A)") "", &
1119 676 : "-------------------------------------------------------------------------------", &
1120 676 : "- TDDFPT WAVEFUNCTION OPTIMIZATION -", &
1121 1352 : "-------------------------------------------------------------------------------"
1122 :
1123 676 : WRITE (log_unit, '(/,T11,A,T27,A,T40,A,T62,A)') "Step", "Time", "Convergence", "Conv. states"
1124 676 : WRITE (log_unit, '(1X,79("-"))')
1125 : END IF
1126 :
1127 1352 : CALL cp_add_iter_level(logger%iter_info, "TDDFT_SCF")
1128 :
1129 : DO
1130 : ! *** perform Davidson iterations ***
1131 : conv = tddfpt_davidson_solver( &
1132 : evects=evects, &
1133 : evals=evals, &
1134 : S_evects=S_evects, &
1135 : gs_mos=gs_mos, &
1136 : tddfpt_control=tddfpt_control, &
1137 : matrix_ks=matrix_ks, &
1138 : qs_env=qs_env, &
1139 : kernel_env=kernel_env, &
1140 : sub_env=sub_env, &
1141 : logger=logger, &
1142 : iter_unit=log_unit, &
1143 : energy_unit=energy_unit, &
1144 : tddfpt_print_section=tddfpt_print_section, &
1145 1440 : work_matrices=work_matrices)
1146 :
1147 : ! at this point at least one of the following conditions are met:
1148 : ! a) convergence criteria has been achieved;
1149 : ! b) maximum number of iterations has been reached;
1150 : ! c) Davidson iterations must be restarted due to lack of Krylov vectors
1151 :
1152 1440 : CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter)
1153 : ! terminate the loop if either (a) or (b) is true ...
1154 1440 : IF ((conv <= tddfpt_control%conv) .OR. iter >= niters) EXIT
1155 :
1156 : ! ... otherwise restart Davidson iterations
1157 522 : evals = 0.0_dp
1158 1440 : IF (log_unit > 0) THEN
1159 44 : WRITE (log_unit, '(1X,25("-"),1X,A,1X,25("-"))') "Restart Davidson iterations"
1160 44 : CALL m_flush(log_unit)
1161 : END IF
1162 : END DO
1163 :
1164 : ! write TDDFPT restart file at the last iteration if requested to do so
1165 1352 : CALL cp_iterate(logger%iter_info, increment=0, last=.TRUE.)
1166 : CALL tddfpt_write_restart(evects=evects, &
1167 : evals=evals, &
1168 : gs_mos=gs_mos, &
1169 : logger=logger, &
1170 1352 : tddfpt_print_section=tddfpt_print_section)
1171 :
1172 1352 : CALL cp_rm_iter_level(logger%iter_info, "TDDFT_SCF")
1173 :
1174 : ! print convergence summary
1175 1352 : IF (log_unit > 0) THEN
1176 676 : CALL integer_to_string(iter, nstates_str)
1177 676 : IF (conv <= tddfpt_control%conv) THEN
1178 676 : WRITE (log_unit, "(1X,A)") "", &
1179 676 : "-------------------------------------------------------------------------------", &
1180 676 : "- TDDFPT run converged in "//TRIM(nstates_str)//" iteration(s) ", &
1181 1352 : "-------------------------------------------------------------------------------"
1182 : ELSE
1183 0 : WRITE (log_unit, "(1X,A)") "", &
1184 0 : "-------------------------------------------------------------------------------", &
1185 0 : "- TDDFPT run did NOT converge after "//TRIM(nstates_str)//" iteration(s) ", &
1186 0 : "-------------------------------------------------------------------------------"
1187 : END IF
1188 : END IF
1189 :
1190 : CALL cp_print_key_finished_output(energy_unit, logger, &
1191 1352 : tddfpt_print_section, "DETAILED_ENERGY")
1192 : CALL cp_print_key_finished_output(log_unit, logger, &
1193 1352 : tddfpt_print_section, "ITERATION_INFO")
1194 : ELSE
1195 : CALL cp_warn(__LOCATION__, &
1196 0 : "Skipping TDDFPT wavefunction optimization")
1197 : END IF
1198 :
1199 : IF (ASSOCIATED(matrix_ks_oep)) THEN
1200 : IF (tddfpt_control%dipole_form == tddfpt_dipole_velocity) THEN
1201 : CALL cp_warn(__LOCATION__, &
1202 : "Transition dipole moments and oscillator strengths are likely to be incorrect "// &
1203 : "when computed using an orbital energy correction XC-potential together with "// &
1204 : "the velocity form of dipole transition integrals")
1205 : END IF
1206 : END IF
1207 :
1208 : ! *** print summary information ***
1209 1352 : log_unit = cp_logger_get_default_io_unit(logger)
1210 :
1211 : namd_print_section => section_vals_get_subs_vals( &
1212 : tddfpt_print_section, &
1213 1352 : "NAMD_PRINT")
1214 1352 : CALL section_vals_get(namd_print_section, explicit=explicit)
1215 1352 : IF (explicit) THEN
1216 : CALL tddfpt_write_newtonx_output(evects, &
1217 : evals, &
1218 : gs_mos, &
1219 : logger, &
1220 : tddfpt_print_section, &
1221 : matrix_s(1)%matrix, &
1222 : S_evects, &
1223 2 : sub_env)
1224 : END IF
1225 4056 : ALLOCATE (ostrength(nstates))
1226 4714 : ostrength = 0.0_dp
1227 : CALL tddfpt_print_summary(log_unit, &
1228 : evects, &
1229 : evals, &
1230 : gs_mos, &
1231 : ostrength, &
1232 : mult, &
1233 : dipole_op_mos_occ, &
1234 1352 : tddfpt_control%dipole_form)
1235 : CALL tddfpt_print_excitation_analysis( &
1236 : log_unit, &
1237 : evects, &
1238 : evals, &
1239 : gs_mos, &
1240 : matrix_s(1)%matrix, &
1241 : tddfpt_control%spinflip, &
1242 1352 : min_amplitude=tddfpt_control%min_excitation_amplitude)
1243 : CALL tddfpt_print_nto_analysis(qs_env, &
1244 : evects, evals, &
1245 : ostrength, &
1246 : gs_mos, &
1247 : matrix_s(1)%matrix, &
1248 1352 : tddfpt_print_section)
1249 1352 : IF (tddfpt_control%do_exciton_descriptors) THEN
1250 : CALL tddfpt_print_exciton_descriptors( &
1251 : log_unit, &
1252 : evects, &
1253 : gs_mos, &
1254 : matrix_s(1)%matrix, &
1255 : tddfpt_control%do_directional_exciton_descriptors, &
1256 2 : qs_env)
1257 : END IF
1258 :
1259 1352 : IF (tddfpt_control%do_lrigpw) THEN
1260 : CALL lri_print_stat(qs_env, &
1261 : ltddfpt=.TRUE., &
1262 10 : tddfpt_lri_env=kernel_env%full_kernel%lri_env)
1263 : END IF
1264 :
1265 1352 : CALL timestop(handle)
1266 5408 : END SUBROUTINE tddfpt_energies
1267 :
1268 : ! **************************************************************************************************
1269 : !> \brief Perform singlet and triplet computations for subsequent TDDFPT-SOC calculation.
1270 : !> \param qs_env Quickstep environment
1271 : !> \param nstates number of requested exited states
1272 : !> \param work_matrices ...
1273 : !> \param tddfpt_control ...
1274 : !> \param logger ...
1275 : !> \param tddfpt_print_section ...
1276 : !> \param evects Eigenvector of the requested multiplicity
1277 : !> \param evals Eigenvalue of the requested multiplicity
1278 : !> \param ostrength Oscillatorstrength
1279 : !> \param gs_mos ...
1280 : !> \param tddfpt_section ...
1281 : !> \param S_evects ...
1282 : !> \param matrix_s ...
1283 : !> \param kernel_env ...
1284 : !> \param matrix_ks ...
1285 : !> \param sub_env ...
1286 : !> \param dipole_op_mos_occ ...
1287 : !> \param lmult_tmp ...
1288 : !> \param xc_section ...
1289 : !> \param full_kernel_env ...
1290 : !> \param kernel_env_admm_aux ...
1291 : !> \par History
1292 : !> * 02.2023 created [Jan-Robert Vogt]
1293 : !> \note Based on tddfpt2_methods and xas_tdp_utils.
1294 : !> \note only the values of one multiplicity will be passed back for force calculations!
1295 : ! **************************************************************************************************
1296 :
1297 10 : SUBROUTINE tddfpt_soc_energies(qs_env, nstates, work_matrices, &
1298 : tddfpt_control, logger, tddfpt_print_section, &
1299 : evects, evals, ostrength, &
1300 : gs_mos, tddfpt_section, S_evects, matrix_s, kernel_env, matrix_ks, &
1301 : sub_env, dipole_op_mos_occ, lmult_tmp, xc_section, full_kernel_env, &
1302 : kernel_env_admm_aux)
1303 :
1304 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1305 : INTEGER, INTENT(in) :: nstates
1306 : TYPE(tddfpt_work_matrices) :: work_matrices
1307 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
1308 : TYPE(cp_logger_type), POINTER :: logger
1309 : TYPE(section_vals_type), POINTER :: tddfpt_print_section
1310 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: evects
1311 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals, ostrength
1312 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1313 : POINTER :: gs_mos
1314 : TYPE(section_vals_type), POINTER :: tddfpt_section
1315 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: S_evects
1316 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1317 : TYPE(kernel_env_type) :: kernel_env
1318 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1319 : TYPE(tddfpt_subgroup_env_type) :: sub_env
1320 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: dipole_op_mos_occ
1321 : LOGICAL, INTENT(in) :: lmult_tmp
1322 : TYPE(section_vals_type), POINTER :: xc_section
1323 : TYPE(full_kernel_env_type), TARGET :: full_kernel_env, kernel_env_admm_aux
1324 :
1325 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_soc_energies'
1326 :
1327 : INTEGER :: handle, ispin, istate, log_unit, mult, &
1328 : nspins
1329 : LOGICAL :: do_sf
1330 10 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_mult, ostrength_mult
1331 10 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: evects_mult
1332 :
1333 10 : CALL timeset(routineN, handle)
1334 :
1335 : log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, &
1336 : "PROGRAM_BANNER", &
1337 10 : extension=".tddfptLog")
1338 10 : CALL tddfpt_soc_header(log_unit)
1339 :
1340 10 : nspins = SIZE(gs_mos)
1341 96 : ALLOCATE (evects_mult(nspins, nstates))
1342 30 : ALLOCATE (evals_mult(nstates))
1343 10 : do_sf = tddfpt_control%spinflip /= no_sf_tddfpt
1344 :
1345 : ! First multiplicity
1346 10 : IF (lmult_tmp) THEN
1347 2 : IF (log_unit > 0) THEN
1348 1 : WRITE (log_unit, "(1X,A)") "", &
1349 1 : "-------------------------------------------------------------------------------", &
1350 1 : "- TDDFPT SINGLET ENERGIES -", &
1351 2 : "-------------------------------------------------------------------------------"
1352 : END IF
1353 2 : mult = 1
1354 : ELSE
1355 8 : IF (log_unit > 0) THEN
1356 4 : WRITE (log_unit, "(1X,A)") "", &
1357 4 : "-------------------------------------------------------------------------------", &
1358 4 : "- TDDFPT TRIPLET ENERGIES -", &
1359 8 : "-------------------------------------------------------------------------------"
1360 : END IF
1361 8 : mult = 3
1362 : END IF
1363 :
1364 : CALL tddfpt_energies(qs_env, nstates, nspins, work_matrices, tddfpt_control, logger, &
1365 : tddfpt_print_section, evects_mult, evals_mult, &
1366 : gs_mos, tddfpt_section, S_evects, matrix_s, &
1367 : kernel_env, matrix_ks, sub_env, ostrength_mult, &
1368 : dipole_op_mos_occ, mult, xc_section, full_kernel_env, &
1369 10 : kernel_env_admm_aux)
1370 :
1371 : ! Clean up in between for full kernel
1372 10 : IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
1373 10 : IF (tddfpt_control%do_admm) CALL release_kernel_env(kernel_env%admm_kernel)
1374 10 : CALL release_kernel_env(kernel_env%full_kernel)
1375 10 : CALL tddfpt_release_work_matrices(work_matrices, sub_env)
1376 : CALL tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, &
1377 : tddfpt_control%do_hfx, &
1378 : tddfpt_control%do_admm, tddfpt_control%do_hfxlr, &
1379 10 : tddfpt_control%do_exck, do_sf, qs_env, sub_env)
1380 : END IF
1381 :
1382 38 : DO istate = 1, nstates
1383 66 : DO ispin = 1, nspins
1384 56 : CALL cp_fm_release(S_evects(ispin, istate))
1385 : END DO
1386 : END DO
1387 :
1388 38 : DO istate = 1, nstates
1389 66 : DO ispin = 1, nspins
1390 : CALL fm_pool_create_fm( &
1391 : work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
1392 56 : S_evects(ispin, istate))
1393 : END DO
1394 : END DO
1395 :
1396 10 : tddfpt_control%rks_triplets = lmult_tmp
1397 :
1398 : ! Second multiplicity
1399 10 : IF (lmult_tmp) THEN
1400 2 : IF (log_unit > 0) THEN
1401 1 : WRITE (log_unit, "(1X,A)") "", &
1402 1 : " singlet excitations finished ", &
1403 1 : " ", &
1404 1 : "-------------------------------------------------------------------------------", &
1405 1 : "- TDDFPT TRIPLET ENERGIES -", &
1406 2 : "-------------------------------------------------------------------------------"
1407 : END IF !log_unit
1408 2 : mult = 3
1409 : ELSE
1410 8 : IF (log_unit > 0) THEN
1411 4 : WRITE (log_unit, "(1X,A)") "", &
1412 4 : " triplet excitations finished ", &
1413 4 : " ", &
1414 4 : "-------------------------------------------------------------------------------", &
1415 4 : "- TDDFPT SINGLET ENERGIES -", &
1416 8 : "-------------------------------------------------------------------------------"
1417 : END IF !log_unit
1418 8 : mult = 1
1419 : END IF
1420 :
1421 : CALL tddfpt_energies(qs_env, nstates, nspins, work_matrices, tddfpt_control, logger, &
1422 : tddfpt_print_section, evects, evals, &
1423 : gs_mos, tddfpt_section, S_evects, matrix_s, &
1424 : kernel_env, matrix_ks, sub_env, ostrength, &
1425 : dipole_op_mos_occ, mult, xc_section, full_kernel_env, &
1426 10 : kernel_env_admm_aux)
1427 :
1428 : ! Compute perturbative SOC correction
1429 : ! Order should always be singlet triplet in tddfpt_soc
1430 10 : IF (lmult_tmp) THEN
1431 2 : CALL tddfpt_soc(qs_env, evals_mult, evals, evects_mult, evects, gs_mos) !mult=singlet
1432 : ELSE
1433 8 : CALL tddfpt_soc(qs_env, evals, evals_mult, evects, evects_mult, gs_mos) !mult=triplet
1434 : END IF
1435 :
1436 : ! deallocate the additional multiplicity
1437 20 : DO ispin = 1, SIZE(evects_mult, 1)
1438 48 : DO istate = 1, SIZE(evects_mult, 2)
1439 38 : CALL cp_fm_release(evects_mult(ispin, istate))
1440 : END DO
1441 : END DO
1442 10 : DEALLOCATE (evects_mult, evals_mult, ostrength_mult)
1443 :
1444 10 : CALL timestop(handle)
1445 :
1446 20 : END SUBROUTINE tddfpt_soc_energies
1447 :
1448 : ! **************************************************************************************************
1449 : !> \brief ...
1450 : !> \param qs_env ...
1451 : !> \param gs_mos ...
1452 : !> \param tddfpt_control ...
1453 : !> \param tddfpt_section ...
1454 : !> \param iounit ...
1455 : ! **************************************************************************************************
1456 1342 : SUBROUTINE init_res_method(qs_env, gs_mos, tddfpt_control, tddfpt_section, iounit)
1457 :
1458 : TYPE(qs_environment_type), POINTER :: qs_env
1459 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1460 : POINTER :: gs_mos
1461 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
1462 : TYPE(section_vals_type), POINTER :: tddfpt_section
1463 : INTEGER, INTENT(IN) :: iounit
1464 :
1465 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_res_method'
1466 :
1467 : INTEGER :: handle, i, io, ispin, nao, nmo, nmol, &
1468 : nspins
1469 1342 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: orblist
1470 1342 : INTEGER, DIMENSION(:), POINTER :: mollist
1471 : LOGICAL :: do_res, do_sf, ew1, ew2, ew3, ewcut, lms
1472 : REAL(KIND=dp) :: eclow, ecup, eint, emo
1473 1342 : REAL(KIND=dp), DIMENSION(:), POINTER :: rvint
1474 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1475 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1476 : TYPE(section_vals_type), POINTER :: res_section
1477 :
1478 1342 : CALL timeset(routineN, handle)
1479 :
1480 1342 : res_section => section_vals_get_subs_vals(tddfpt_section, "REDUCED_EXCITATION_SPACE")
1481 1342 : CALL section_vals_val_get(res_section, "_SECTION_PARAMETERS_", l_val=do_res)
1482 :
1483 : ! spin flip TDA
1484 1342 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1485 : do_sf = .FALSE.
1486 : ELSE
1487 22 : do_sf = .TRUE.
1488 : END IF
1489 :
1490 1342 : nspins = SIZE(gs_mos)
1491 1342 : IF (.NOT. do_res) THEN
1492 2808 : DO ispin = 1, nspins
1493 1480 : nmo = gs_mos(ispin)%nmo_occ
1494 1480 : tddfpt_control%nactive(ispin) = nmo
1495 1480 : gs_mos(ispin)%nmo_active = nmo
1496 4440 : ALLOCATE (gs_mos(ispin)%index_active(nmo))
1497 10434 : DO i = 1, nmo
1498 9106 : gs_mos(ispin)%index_active(i) = i
1499 : END DO
1500 : END DO
1501 : ELSE
1502 14 : IF (iounit > 0) THEN
1503 7 : WRITE (iounit, "(/,1X,27('='),A,26('='))") ' REDUCED EXCITATION SPACE '
1504 : END IF
1505 14 : CALL section_vals_val_get(res_section, "ENERGY_WINDOW", explicit=ew1)
1506 14 : CALL section_vals_val_get(res_section, "UPPER_ENERGY_CUTOFF", explicit=ew2)
1507 14 : CALL section_vals_val_get(res_section, "LOWER_ENERGY_CUTOFF", explicit=ew3)
1508 14 : ewcut = (ew1 .OR. ew2 .OR. ew3)
1509 14 : CALL section_vals_val_get(res_section, "MOLECULE_LIST", explicit=lms)
1510 :
1511 14 : CALL section_vals_val_get(res_section, "ENERGY_WINDOW", r_vals=rvint)
1512 14 : CPASSERT(SIZE(rvint) == 2)
1513 14 : eclow = rvint(1)
1514 14 : ecup = rvint(2)
1515 14 : CALL section_vals_val_get(res_section, "UPPER_ENERGY_CUTOFF", r_val=eint)
1516 14 : ecup = MIN(ecup, eint)
1517 14 : CALL section_vals_val_get(res_section, "LOWER_ENERGY_CUTOFF", r_val=eint)
1518 14 : eclow = MAX(eclow, eint)
1519 14 : IF (ewcut .AND. (iounit > 0)) THEN
1520 7 : IF (eclow < -1.E8_dp .AND. ecup > 1.E8_dp) THEN
1521 : WRITE (iounit, "(1X,A,T51,A10,T71,A10)") &
1522 0 : 'Orbital Energy Window [eV]', " -Inf", " Inf"
1523 : ELSE IF (eclow < -1.E8_dp) THEN
1524 : WRITE (iounit, "(1X,A,T51,A10,T71,F10.4)") &
1525 3 : 'Orbital Energy Window [eV]', " -Inf", evolt*ecup
1526 4 : ELSE IF (ecup > 1.E8_dp) THEN
1527 : WRITE (iounit, "(1X,A,T51,F10.4,T71,A10)") &
1528 1 : 'Orbital Energy Window [eV]', evolt*eclow, " Inf"
1529 : ELSE
1530 : WRITE (iounit, "(1X,A,T51,F10.4,T71,F10.4)") &
1531 3 : 'Orbital Energy Window [eV]', evolt*eclow, evolt*ecup
1532 : END IF
1533 : END IF
1534 :
1535 14 : nmol = 0
1536 14 : IF (lms) THEN
1537 0 : CALL section_vals_val_get(res_section, "MOLECULE_LIST", i_vals=mollist)
1538 0 : nmol = SIZE(mollist)
1539 0 : WRITE (iounit, "(1X,A)") 'List of Selected Molecules'
1540 0 : WRITE (iounit, "(1X,15(I5))") mollist(1:nmol)
1541 : END IF
1542 :
1543 28 : DO ispin = 1, nspins
1544 28 : tddfpt_control%nactive(ispin) = gs_mos(ispin)%nmo_occ
1545 : END DO
1546 42 : nmo = MAXVAL(tddfpt_control%nactive)
1547 56 : ALLOCATE (orblist(nmo, nspins))
1548 104 : orblist = 0
1549 :
1550 14 : IF (lms) THEN
1551 : ! ignore for now
1552 0 : orblist = 1
1553 0 : DO ispin = 1, nspins
1554 0 : CPASSERT(.NOT. ASSOCIATED(gs_mos(ispin)%evals_occ_matrix))
1555 : END DO
1556 14 : ELSEIF (ewcut) THEN
1557 : ! Filter orbitals wrt energy window
1558 28 : DO ispin = 1, nspins
1559 104 : DO i = 1, gs_mos(ispin)%nmo_occ
1560 76 : emo = gs_mos(ispin)%evals_occ(i)
1561 90 : IF (emo > eclow .AND. emo < ecup) orblist(i, ispin) = 1
1562 : END DO
1563 : END DO
1564 : ELSE
1565 0 : orblist = 1
1566 : END IF
1567 :
1568 : ! count active orbitals
1569 104 : nmo = SUM(orblist)
1570 14 : IF (nmo == 0) THEN
1571 0 : CPABORT("RSE TDA: no active orbitals selected.")
1572 : END IF
1573 28 : DO ispin = 1, nspins
1574 90 : nmo = SUM(orblist(:, ispin))
1575 14 : tddfpt_control%nactive(ispin) = nmo
1576 14 : gs_mos(ispin)%nmo_active = nmo
1577 42 : ALLOCATE (gs_mos(ispin)%index_active(nmo))
1578 14 : io = 0
1579 104 : DO i = 1, SIZE(ORBLIST, 1)
1580 90 : IF (orblist(i, ispin) == 1) THEN
1581 24 : io = io + 1
1582 24 : gs_mos(ispin)%index_active(io) = i
1583 : END IF
1584 : END DO
1585 : END DO
1586 14 : DEALLOCATE (orblist)
1587 :
1588 14 : IF (lms) THEN
1589 : ! output information
1590 : ELSE
1591 14 : IF (iounit > 0) THEN
1592 7 : WRITE (iounit, "(1X,A)") 'List of Selected States'
1593 7 : IF (nspins == 1) THEN
1594 7 : WRITE (iounit, "(A,T67,A)") ' Active State Orbital', 'Orbital Energy'
1595 19 : DO i = 1, gs_mos(1)%nmo_active
1596 12 : io = gs_mos(1)%index_active(i)
1597 19 : WRITE (iounit, "(T8,I6,T21,I6,T61,F20.4)") i, io, gs_mos(1)%evals_occ(io)*evolt
1598 : END DO
1599 : ELSE
1600 0 : DO ispin = 1, nspins
1601 0 : WRITE (iounit, "(1X,A,I2)") 'Spin ', ispin
1602 0 : WRITE (iounit, "(A,T67,A)") ' Active State Orbital', 'Orbital Energy'
1603 0 : DO i = 1, gs_mos(ispin)%nmo_active
1604 0 : io = gs_mos(ispin)%index_active(i)
1605 0 : WRITE (iounit, "(T8,I6,T21,I6,T61,F20.4)") i, io, gs_mos(ispin)%evals_occ(io)*evolt
1606 : END DO
1607 : END DO
1608 : END IF
1609 : END IF
1610 : END IF
1611 :
1612 14 : IF (do_sf) THEN
1613 0 : CPABORT("Restricted Active Space with spin flip TDA NYA")
1614 : END IF
1615 :
1616 56 : IF (iounit > 0) THEN
1617 7 : WRITE (iounit, "(1X,79('='))")
1618 : END IF
1619 : END IF
1620 :
1621 : ! Allocate mos_active
1622 2836 : DO ispin = 1, nspins
1623 1494 : CALL get_qs_env(qs_env, blacs_env=blacs_env)
1624 1494 : CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, nrow_global=nao)
1625 1494 : nmo = gs_mos(ispin)%nmo_active
1626 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=gs_mos(ispin)%mos_occ%matrix_struct, &
1627 1494 : ncol_global=nmo, context=blacs_env)
1628 1494 : NULLIFY (gs_mos(ispin)%mos_active)
1629 1494 : ALLOCATE (gs_mos(ispin)%mos_active)
1630 1494 : CALL cp_fm_create(gs_mos(ispin)%mos_active, fm_struct)
1631 1494 : CALL cp_fm_struct_release(fm_struct)
1632 : ! copy the active orbitals
1633 4330 : IF (gs_mos(ispin)%nmo_active == gs_mos(ispin)%nmo_occ) THEN
1634 9106 : DO i = 1, gs_mos(ispin)%nmo_active
1635 9106 : CPASSERT(i == gs_mos(ispin)%index_active(i))
1636 : END DO
1637 : CALL cp_fm_to_fm_submat(gs_mos(ispin)%mos_occ, gs_mos(ispin)%mos_active, &
1638 1480 : nao, nmo, 1, 1, 1, 1)
1639 : ELSE
1640 38 : DO i = 1, gs_mos(ispin)%nmo_active
1641 24 : io = gs_mos(ispin)%index_active(i)
1642 : CALL cp_fm_to_fm_submat(gs_mos(ispin)%mos_occ, gs_mos(ispin)%mos_active, &
1643 38 : nao, 1, 1, io, 1, i)
1644 : END DO
1645 : END IF
1646 : END DO
1647 :
1648 1342 : CALL timestop(handle)
1649 :
1650 1342 : END SUBROUTINE init_res_method
1651 :
1652 : END MODULE qs_tddfpt2_methods
|
