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