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