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