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