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