(git:fd1133a)
Loading...
Searching...
No Matches
qs_tddfpt2_methods.F
Go to the documentation of this file.
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
10 USE admm_types, ONLY: admm_type,&
13 USE bibliography, ONLY: grimme2013,&
17 cite_reference
18 USE cell_types, ONLY: cell_type
23 USE cp_dbcsr_api, ONLY: dbcsr_p_type
29 USE cp_fm_types, ONLY: cp_fm_create,&
44 USE header, ONLY: tddfpt_header,&
49 USE input_constants, ONLY: &
58 USE kinds, ONLY: dp
62 USE machine, ONLY: m_flush
66 USE physcon, ONLY: evolt
75 USE qs_mo_types, ONLY: mo_set_type
96 USE qs_tddfpt2_soc, ONLY: tddfpt_soc
116 USE rixs_types, ONLY: rixs_env_type,&
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
134
135! **************************************************************************************************
136
137CONTAINS
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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals, ostrength
167 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 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: my_active, my_mos
172 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: dipole_op_mos_occ, evects, s_evects
173 TYPE(cp_logger_type), POINTER :: logger
174 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_oep, matrix_s, &
175 matrix_s_aux_fit, &
176 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 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
182 TYPE(mp_para_env_type), POINTER :: para_env
183 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
184 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 POINTER :: gs_mos
196 TYPE(tddfpt_subgroup_env_type) :: sub_env
197 TYPE(tddfpt_work_matrices) :: work_matrices
198
199 CALL timeset(routinen, handle)
200
201 NULLIFY (logger)
202 logger => cp_get_default_logger()
203
204 ! input section print/xc
205 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 lri_section, hfxsr_section)
210
211 log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, "PROGRAM_BANNER", &
212 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 do_rixs=do_rixs)
223
224 IF (do_rixs) THEN
225 tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%RIXS%TDDFPT")
226 NULLIFY (rixs_control, valence_state)
227 rixs_control => dft_control%rixs_control
228 tddfpt_control => rixs_control%tddfpt2_control
229 valence_state => rixs_env%valence_state
230 ELSE
231 tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT")
232 tddfpt_control => dft_control%tddfpt2_control
233 END IF
234
235 tddfpt_control%do_hfx = do_hfx
236 tddfpt_control%do_admm = do_admm
237 tddfpt_control%do_hfxsr = do_hfxsr
238 tddfpt_control%hfxsr_primbas = 0
239 tddfpt_control%hfxsr_re_int = .true.
240 tddfpt_control%do_hfxlr = do_hfxlr
241 tddfpt_control%do_exck = do_exck
242 do_sf = tddfpt_control%spinflip /= no_sf_tddfpt
243 IF (do_sf) CALL cite_reference(hernandez2025)
244 IF (tddfpt_control%do_hfxlr) THEN
245 kernel_section => section_vals_get_subs_vals(tddfpt_section, "XC%HFX_KERNEL%HFXLR")
246 CALL section_vals_val_get(kernel_section, "RCUT", r_val=tddfpt_control%hfxlr_rcut)
247 CALL section_vals_val_get(kernel_section, "SCALE", r_val=tddfpt_control%hfxlr_scale)
248 END IF
249
250 soc_section => section_vals_get_subs_vals(tddfpt_section, "SOC")
251 CALL section_vals_get(soc_section, explicit=do_soc)
252
253 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 lmult_tmp = tddfpt_control%rks_triplets
257 tddfpt_control%rks_triplets = .NOT. (tddfpt_control%rks_triplets)
258 END IF
259
260 CALL cite_reference(iannuzzi2005)
261 IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
262 CALL cite_reference(grimme2013)
263 CALL cite_reference(grimme2016)
264 END IF
265
266 CALL tddfpt_header(log_unit)
267 CALL kernel_info(log_unit, dft_control, tddfpt_control, xc_section)
268
269 ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals
270 NULLIFY (gs_mos)
271 CALL tddfpt_init_mos(qs_env, gs_mos, log_unit)
272
273 ! obtain smeared occupation numbers
274 IF (tddfpt_control%do_smearing) THEN
275 CALL tddfpt_smeared_occupation(qs_env, gs_mos, log_unit)
276 END IF
277
278 ! obtain corrected KS-matrix
279 CALL tddfpt_oecorr(qs_env, gs_mos, matrix_ks_oep)
280
281 IF ((tddfpt_control%do_lrigpw) .AND. &
282 (tddfpt_control%kernel /= tddfpt_kernel_full)) THEN
283 CALL cp_abort(__location__, "LRI only implemented for full kernel")
284 END IF
285
286 IF (ASSOCIATED(matrix_ks_oep)) matrix_ks => matrix_ks_oep
287
288 ! determine active orbitals
289 ! default is all occupied MOs
290 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 qs_env)
297
298 nspins = SIZE(gs_mos)
299 ! multiplicity of molecular system
300 IF (nspins > 1) THEN
301 mult = abs(SIZE(gs_mos(1)%evals_occ) - SIZE(gs_mos(2)%evals_occ)) + 1
302 IF (mult > 2) &
303 CALL cp_warn(__location__, "There is a convergence issue for multiplicity >= 3")
304 ELSE
305 IF (tddfpt_control%rks_triplets) THEN
306 mult = 3
307 ELSE
308 mult = 1
309 END IF
310 END IF
311
312 ! split mpi communicator
313 ALLOCATE (my_mos(nspins), my_active(nspins))
314 DO ispin = 1, nspins
315 my_mos(ispin) = gs_mos(ispin)%mos_occ
316 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 kernel=tddfpt_control%kernel)
321 DEALLOCATE (my_mos, my_active)
322
323 IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
324 ! create environment for Full Kernel
325 IF (dft_control%qs_control%xtb) THEN
326 cpabort("TDDFPT: xTB only works with sTDA Kernel")
327 END IF
328
329 IF (tddfpt_control%do_hfxsr) THEN
330 kernel_section => section_vals_get_subs_vals(tddfpt_section, "XC%HFX_KERNEL")
331 CALL section_vals_val_get(kernel_section, "HFXSR_PRIMBAS", &
332 i_val=tddfpt_control%hfxsr_primbas)
333 ! basis set
334 CALL create_minbas_set(qs_env, log_unit, basis_type="TDA_HFX", &
335 primitive=tddfpt_control%hfxsr_primbas)
336 ! admm control
337 ALLOCATE (full_kernel_env%admm_control)
338 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 full_kernel_env%admm_control%aux_exch_func = do_admm_aux_exch_func_none
342 ! hfx section
343 full_kernel_env%hfxsr_section => hfxsr_section
344 !
345 CALL aux_admm_init(qs_env, mos, full_kernel_env%admm_env, &
346 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 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 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 para_env=para_env)
356 CALL hfx_create(full_kernel_env%x_data, para_env, hfxsr_section, atomic_kind_set, &
357 qs_kind_set, particle_set, dft_control, cell, orb_basis="TDA_HFX")
358 END IF
359
360 ! allocate pools and work matrices
361 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 CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=nocc)
366 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
367 CALL cp_fm_get_info(gs_mos(1)%mos_virt, ncol_global=nvirt)
368 ELSE
369 CALL cp_fm_get_info(gs_mos(2)%mos_virt, ncol_global=nvirt)
370 END IF
371 nstate_max = nocc*nvirt
372 IF (nstates > nstate_max) THEN
373 cpwarn("NUMBER OF EXCITED STATES COULD LEAD TO PROBLEMS!")
374 cpwarn("Experimental: CHANGED NSTATES TO ITS MAXIMUM VALUE!")
375 nstates = nstate_max
376 tddfpt_control%nstates = nstate_max
377 END IF
378 CALL tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, &
379 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 kernel_env%full_kernel => full_kernel_env
383 kernel_env%admm_kernel => kernel_env_admm_aux
384 NULLIFY (kernel_env%stda_kernel)
385 IF (do_hfxsr) THEN
386 ! work matrices for SR HFX
387 CALL hfxsr_create_work_matrices(work_matrices, qs_env, full_kernel_env%admm_env)
388 END IF
389 IF (do_hfxlr) THEN
390 ! calculate S_half and Lowdin MO coefficients
391 CALL get_lowdin_mo_coefficients(qs_env, sub_env, work_matrices)
392 END IF
393 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
394 ! setup for kernel_stda outside tddfpt_energies
395 CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
396 nactive = tddfpt_control%nactive
397 CALL allocate_stda_env(qs_env, stda_kernel, nao, nactive)
398 ! sTDA parameters
399 CALL stda_init_param(qs_env, stda_kernel, tddfpt_control%stda_control)
400 ! allocate pools and work matrices
401 nstates = tddfpt_control%nstates
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 work_matrices, tddfpt_control)
406 !
407 kernel_env%stda_kernel => stda_kernel
408 NULLIFY (kernel_env%full_kernel)
409 NULLIFY (kernel_env%admm_kernel)
410 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
411 ! allocate pools and work matrices
412 nstates = tddfpt_control%nstates
413 CALL stda_create_work_matrices(work_matrices, gs_mos, nstates, qs_env, sub_env)
414 NULLIFY (kernel_env%full_kernel)
415 NULLIFY (kernel_env%admm_kernel)
416 NULLIFY (kernel_env%stda_kernel)
417 END IF
418
419 IF (do_sf) THEN
420 ! only alpha -> beta excitations are considered in spin-flip TDDFT
421 ALLOCATE (evects(1, nstates))
422 ELSE
423 ALLOCATE (evects(nspins, nstates))
424 END IF
425 ALLOCATE (evals(nstates))
426 ALLOCATE (s_evects(SIZE(evects, 1), nstates))
427
428 DO istate = 1, nstates
429 DO ispin = 1, SIZE(evects, 1)
430 CALL fm_pool_create_fm( &
431 work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
432 s_evects(ispin, istate))
433 END DO
434 END DO
435
436 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 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 kernel_env_admm_aux)
450 END IF
451
452 !print forces for selected states
453 IF (calc_forces) THEN
454 CALL tddfpt_print_forces(qs_env, evects, evals, ostrength, &
455 tddfpt_print_section, gs_mos, &
456 kernel_env, sub_env, work_matrices)
457 END IF
458
459 ! excited state potential energy surface
460 IF (qs_env%excited_state) THEN
461 IF (sub_env%is_split) THEN
462 CALL cp_abort(__location__, &
463 "Excited state forces not possible when states"// &
464 " are distributed to different CPU pools.")
465 END IF
466 ! for gradients unshifted KS matrix
467 IF (ASSOCIATED(matrix_ks_oep)) CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
468 CALL get_qs_env(qs_env, exstate_env=ex_env)
469 state_change = .false.
470 IF (ex_env%state > 0) THEN
471 my_state = ex_env%state
472 ELSEIF (ex_env%state < 0) THEN
473 ! state following
474 ALLOCATE (my_mos(nspins))
475 DO ispin = 1, nspins
476 my_mos(ispin) = gs_mos(ispin)%mos_occ
477 END DO
478 my_state = abs(ex_env%state)
479 CALL assign_state(qs_env, matrix_s, evects, my_mos, ex_env%wfn_history, my_state)
480 DEALLOCATE (my_mos)
481 IF (my_state /= abs(ex_env%state)) THEN
482 state_change = .true.
483 old_state = abs(ex_env%state)
484 END IF
485 ex_env%state = -my_state
486 ELSE
487 CALL cp_warn(__location__, &
488 "Active excited state not assigned. Use the first state.")
489 my_state = 1
490 END IF
491 cpassert(my_state > 0)
492 IF (my_state > nstates) THEN
493 CALL cp_warn(__location__, &
494 "There were not enough excited states calculated.")
495 cpabort("excited state potential energy surface")
496 END IF
497 !
498 ! energy
499 ex_env%evalue = evals(my_state)
500 ! excitation vector
501 CALL cp_fm_release(ex_env%evect)
502 ALLOCATE (ex_env%evect(SIZE(evects, 1)))
503 DO ispin = 1, SIZE(evects, 1)
504 CALL cp_fm_get_info(matrix=evects(ispin, 1), &
505 matrix_struct=matrix_struct)
506 CALL cp_fm_create(ex_env%evect(ispin), matrix_struct)
507 CALL cp_fm_to_fm(evects(ispin, my_state), ex_env%evect(ispin))
508 END DO
509
510 IF (log_unit > 0) THEN
511 gsval = ex_env%wfn_history%gsval
512 gsmin = ex_env%wfn_history%gsmin
513 xsval = ex_env%wfn_history%xsval
514 WRITE (log_unit, "(1X,A,T40,F10.6,A,T62,F10.6,A)") "Ground state orbital alignment:", &
515 gsmin, "[MinVal]", gsval, "[Average]"
516 WRITE (log_unit, "(1X,A,T71,F10.6)") "Excitation vector alignment:", xsval
517 IF (state_change) THEN
518 WRITE (log_unit, "(1X,A,I5,T60,A14,T76,I5)") &
519 "Target state has been changed from state ", &
520 old_state, " to new state ", my_state
521 END IF
522 WRITE (log_unit, "(1X,A,I4,A,F12.5,A)") "Calculate properties for state:", &
523 my_state, " with excitation energy ", ex_env%evalue*evolt, " eV"
524 END IF
525
526 ! Calculate response vector
527 IF (calc_forces) THEN
528 CALL tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, &
529 sub_env, work_matrices)
530 END IF
531 END IF
532
533 ! share evals, evects and mo_coefs with rixs
534 IF (do_rixs) THEN
535 ! copy evals
536 valence_state%nstates = nstates
537 ALLOCATE (valence_state%evals(SIZE(evals)))
538 valence_state%evals(:) = evals(:)
539
540 ALLOCATE (valence_state%evects(nspins, nstates))
541 ALLOCATE (valence_state%mos_occ(nspins))
542 DO ispin = 1, nspins
543 ! copy evects
544 DO istate = 1, nstates
545 CALL cp_fm_get_info(matrix=evects(ispin, istate), &
546 matrix_struct=matrix_struct)
547 CALL cp_fm_create(valence_state%evects(ispin, istate), matrix_struct)
548 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 matrix_struct=matrix_struct)
553 CALL cp_fm_create(valence_state%mos_occ(ispin), matrix_struct)
554 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 CALL cp_fm_release(evects)
560 CALL cp_fm_release(s_evects)
561
562 CALL cp_print_key_finished_output(log_unit, &
563 logger, &
564 tddfpt_print_section, &
565 "PROGRAM_BANNER")
566
567 DEALLOCATE (evals, ostrength)
568
569 IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
570 IF (do_admm) CALL release_kernel_env(kernel_env%admm_kernel)
571 IF (tddfpt_control%do_lrigpw) THEN
572 CALL lri_env_release(kernel_env%full_kernel%lri_env)
573 DEALLOCATE (kernel_env%full_kernel%lri_env)
574 CALL lri_density_release(kernel_env%full_kernel%lri_density)
575 DEALLOCATE (kernel_env%full_kernel%lri_density)
576 END IF
577 CALL release_kernel_env(kernel_env%full_kernel)
578 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
579 CALL deallocate_stda_env(stda_kernel)
580 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
581 !
582 ELSE
583 cpabort('Unknown kernel type')
584 END IF
585 CALL tddfpt_release_work_matrices(work_matrices, sub_env)
586 CALL tddfpt_sub_env_release(sub_env)
587
588 CALL cp_fm_release(dipole_op_mos_occ)
589
590 DO ispin = nspins, 1, -1
591 CALL tddfpt_release_ground_state_mos(gs_mos(ispin))
592 END DO
593 DEALLOCATE (gs_mos)
594
595 IF (ASSOCIATED(matrix_ks_oep)) &
596 CALL dbcsr_deallocate_matrix_set(matrix_ks_oep)
597
598 CALL timestop(handle)
599
600 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 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 NULLIFY (dft_control, input)
634 CALL get_qs_env(qs_env, dft_control=dft_control, input=input)
635 tddfpt_control => dft_control%tddfpt2_control
636
637 ! no k-points
638 IF (dft_control%nimages > 1) cpabort("k-points not implemented for TDDFPT")
639
640 IF (tddfpt_control%nstates <= 0) THEN
641 CALL integer_to_string(tddfpt_control%nstates, nstates_str)
642 CALL cp_warn(__location__, "TDDFPT calculation was requested for "// &
643 trim(nstates_str)//" excited states: nothing to do.")
644 RETURN
645 END IF
646
647 NULLIFY (tddfpt_section, tddfpt_print_section)
648 tddfpt_section => section_vals_get_subs_vals(input, "PROPERTIES%TDDFPT")
649 tddfpt_print_section => section_vals_get_subs_vals(tddfpt_section, "PRINT")
650
651 IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
652 NULLIFY (xc_root)
653 xc_root => section_vals_get_subs_vals(tddfpt_section, "XC")
654 CALL section_vals_get(xc_root, explicit=explicit_root)
655 NULLIFY (xc_section)
656 IF (explicit_root) THEN
657 ! No ADIABATIC_RESCALING option possible
658 NULLIFY (xc_sub)
659 xc_sub => section_vals_get_subs_vals(xc_root, "ADIABATIC_RESCALING")
660 CALL section_vals_get(xc_sub, explicit=exar)
661 IF (exar) THEN
662 CALL cp_warn(__location__, "TDDFPT Kernel with ADIABATIC_RESCALING not possible.")
663 cpabort("TDDFPT Input")
664 END IF
665 ! No GCP_POTENTIAL option possible
666 NULLIFY (xc_sub)
667 xc_sub => section_vals_get_subs_vals(xc_root, "GCP_POTENTIAL")
668 CALL section_vals_get(xc_sub, explicit=exgcp)
669 IF (exgcp) THEN
670 CALL cp_warn(__location__, "TDDFPT Kernel with GCP_POTENTIAL not possible.")
671 cpabort("TDDFPT Input")
672 END IF
673 ! No VDW_POTENTIAL option possible
674 NULLIFY (xc_sub)
675 xc_sub => section_vals_get_subs_vals(xc_root, "VDW_POTENTIAL")
676 CALL section_vals_get(xc_sub, explicit=exvdw)
677 IF (exvdw) THEN
678 CALL cp_warn(__location__, "TDDFPT Kernel with VDW_POTENTIAL not possible.")
679 cpabort("TDDFPT Input")
680 END IF
681 ! No WF_CORRELATION option possible
682 NULLIFY (xc_sub)
683 xc_sub => section_vals_get_subs_vals(xc_root, "WF_CORRELATION")
684 CALL section_vals_get(xc_sub, explicit=exwfn)
685 IF (exwfn) THEN
686 CALL cp_warn(__location__, "TDDFPT Kernel with WF_CORRELATION not possible.")
687 cpabort("TDDFPT Input")
688 END IF
689 ! No XC_POTENTIAL option possible
690 NULLIFY (xc_sub)
691 xc_sub => section_vals_get_subs_vals(xc_root, "XC_POTENTIAL")
692 CALL section_vals_get(xc_sub, explicit=expot)
693 IF (expot) THEN
694 CALL cp_warn(__location__, "TDDFPT Kernel with XC_POTENTIAL not possible.")
695 cpabort("TDDFPT Input")
696 END IF
697 !
698 NULLIFY (xc_sub)
699 xc_sub => section_vals_get_subs_vals(xc_root, "XC_FUNCTIONAL")
700 CALL section_vals_get(xc_sub, explicit=exf)
701 NULLIFY (xc_sub)
702 xc_sub => section_vals_get_subs_vals(xc_root, "XC_KERNEL")
703 CALL section_vals_get(xc_sub, explicit=exk)
704 IF ((exf .AND. exk) .OR. .NOT. (exf .OR. exk)) THEN
705 CALL cp_warn(__location__, "TDDFPT Kernel needs XC_FUNCTIONAL or XC_KERNEL section.")
706 cpabort("TDDFPT Input")
707 END IF
708 NULLIFY (xc_sub)
709 xc_sub => section_vals_get_subs_vals(xc_root, "HF")
710 CALL section_vals_get(xc_sub, explicit=exhf)
711 NULLIFY (xc_sub)
712 xc_sub => section_vals_get_subs_vals(xc_root, "HFX_KERNEL")
713 CALL section_vals_get(xc_sub, explicit=exhfxk)
714 !
715 xc_section => xc_root
716 hfx_section => section_vals_get_subs_vals(xc_section, "HF")
717 CALL section_vals_get(hfx_section, explicit=do_hfx)
718 IF (do_hfx) THEN
719 CALL section_vals_val_get(hfx_section, "FRACTION", r_val=c_hf)
720 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 IF (do_hfx) THEN
724 hfx_section_gs => section_vals_get_subs_vals(input, "DFT%XC%HF")
725 CALL compare_hfx_sections(hfx_section, hfx_section_gs, same_hfx)
726 IF (.NOT. same_hfx) THEN
727 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 do_admm = do_hfx .AND. dft_control%do_admm
732 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 "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 IF (exk) THEN
740 do_exck = .true.
741 ELSE
742 do_exck = .false.
743 END IF
744 IF (exhfxk) THEN
745 xc_sub => section_vals_get_subs_vals(xc_root, "HFX_KERNEL")
746 CALL section_vals_val_get(xc_sub, "DO_HFXSR", l_val=do_hfxsr)
747 xc_sub => section_vals_get_subs_vals(xc_root, "HFX_KERNEL%HFXLR")
748 CALL section_vals_get(xc_sub, explicit=do_hfxlr)
749 ELSE
750 do_hfxsr = .false.
751 do_hfxlr = .false.
752 END IF
753 ELSE
754 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
755 hfx_section => section_vals_get_subs_vals(xc_section, "HF")
756 CALL section_vals_get(hfx_section, explicit=do_hfx)
757 IF (do_hfx) THEN
758 CALL section_vals_val_get(hfx_section, "FRACTION", r_val=c_hf)
759 do_hfx = (c_hf /= 0.0_dp)
760 END IF
761 do_admm = do_hfx .AND. dft_control%do_admm
762 do_exck = .false.
763 do_hfxsr = .false.
764 do_hfxlr = .false.
765 END IF
766 ELSE
767 do_hfx = .false.
768 do_admm = .false.
769 do_exck = .false.
770 do_hfxsr = .false.
771 do_hfxlr = .false.
772 END IF
773
774 ! reset rks_triplets if UKS is in use
775 IF (tddfpt_control%rks_triplets .AND. dft_control%nspins > 1) THEN
776 tddfpt_control%rks_triplets = .false.
777 CALL cp_warn(__location__, "Keyword RKS_TRIPLETS has been ignored for spin-polarised calculations")
778 END IF
779
780 ! lri input
781 IF (tddfpt_control%do_lrigpw) THEN
782 lri_section => section_vals_get_subs_vals(tddfpt_section, "LRIGPW")
783 END IF
784
785 ! set defaults for short range HFX
786 NULLIFY (hfxsr_section)
787 IF (do_hfxsr) THEN
788 hfxsr_section => section_vals_get_subs_vals(tddfpt_section, "XC%HFX_KERNEL%HF")
789 CALL section_vals_get(hfxsr_section, explicit=found)
790 IF (.NOT. found) THEN
791 cpabort("HFXSR option needs &HF section defined")
792 END IF
793 CALL section_vals_val_get(hfxsr_section, "INTERACTION_POTENTIAL%POTENTIAL_TYPE", explicit=found)
794 IF (.NOT. found) THEN
795 CALL section_vals_val_set(hfxsr_section, "INTERACTION_POTENTIAL%POTENTIAL_TYPE", &
797 END IF
798 CALL section_vals_val_get(hfxsr_section, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", explicit=found)
799 IF (.NOT. found) THEN
800 CALL section_vals_val_set(hfxsr_section, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", r_val=7.5589_dp)
801 END IF
802 CALL section_vals_val_get(hfxsr_section, "RI%_SECTION_PARAMETERS_", l_val=found)
803 IF (found) THEN
804 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 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 lsd = (dft_control%nspins > 1)
827 IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
828 ktype = "FULL"
829 IF (log_unit > 0) THEN
830 WRITE (log_unit, "(T2,A,T77,A4)") "KERNEL|", trim(ktype)
831 CALL xc_write(log_unit, xc_section, lsd)
832 IF (tddfpt_control%do_hfx) THEN
833 IF (tddfpt_control%do_admm) THEN
834 WRITE (log_unit, "(T2,A,T62,A19)") "KERNEL|", "ADMM Exact Exchange"
835 IF (tddfpt_control%admm_xc_correction) THEN
836 WRITE (log_unit, "(T2,A,T60,A21)") "KERNEL|", "Apply ADMM Kernel XC Correction"
837 END IF
838 IF (tddfpt_control%admm_symm) THEN
839 WRITE (log_unit, "(T2,A,T60,A21)") "KERNEL|", "Symmetric ADMM Kernel"
840 END IF
841 ELSE
842 WRITE (log_unit, "(T2,A,T67,A14)") "KERNEL|", "Exact Exchange"
843 END IF
844 END IF
845 IF (tddfpt_control%do_hfxsr) THEN
846 WRITE (log_unit, "(T2,A,T43,A38)") "KERNEL|", "Short range HFX approximation"
847 END IF
848 IF (tddfpt_control%do_hfxlr) THEN
849 WRITE (log_unit, "(T2,A,T43,A38)") "KERNEL|", "Long range HFX approximation"
850 END IF
851 IF (tddfpt_control%do_lrigpw) THEN
852 WRITE (log_unit, "(T2,A,T42,A39)") "KERNEL|", "LRI approximation of transition density"
853 END IF
854 END IF
855 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
856 ktype = "sTDA"
857 IF (log_unit > 0) THEN
858 WRITE (log_unit, "(T2,A,T77,A4)") "KERNEL|", trim(ktype)
859 IF (tddfpt_control%stda_control%do_ewald) THEN
860 WRITE (log_unit, "(T2,A,T78,A3)") "KERNEL| Coulomb term uses Ewald summation"
861 ELSE
862 WRITE (log_unit, "(T2,A,T78,A3)") "KERNEL| Coulomb term uses direct summation (MIC)"
863 END IF
864 IF (tddfpt_control%stda_control%do_exchange) THEN
865 WRITE (log_unit, "(T2,A,T78,A3)") "KERNEL| Exact exchange term", "YES"
866 WRITE (log_unit, "(T2,A,T71,F10.3)") "KERNEL| Short range HFX fraction:", &
867 tddfpt_control%stda_control%hfx_fraction
868 ELSE
869 WRITE (log_unit, "(T2,A,T79,A2)") "KERNEL| Exact exchange term", "NO"
870 END IF
871 WRITE (log_unit, "(T2,A,T66,E15.3)") "KERNEL| Transition density filter", &
872 tddfpt_control%stda_control%eps_td_filter
873 END IF
874 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
875 ktype = "NONE"
876 IF (log_unit > 0) THEN
877 WRITE (log_unit, "(T2,A,T77,A4)") "KERNEL|", trim(ktype)
878 END IF
879 ELSE
880 !CPABORT("Unknown kernel")
881 END IF
882 !
883 IF (log_unit > 0) THEN
884 IF (tddfpt_control%rks_triplets) THEN
885 WRITE (log_unit, "(T2,A,T74,A7)") "KERNEL| Spin symmetry of excitations", "Triplet"
886 ELSE IF (lsd) THEN
887 ! Spin-conserving excitations where requested
888 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
889 WRITE (log_unit, "(T2,A,T69,A12)") "KERNEL| Spin symmetry of excitations", "Unrestricted"
890 ! Spin-flip excitations with collinear exchange-correlation kernel requested
891 ELSE IF (tddfpt_control%spinflip == tddfpt_sf_col) THEN
892 WRITE (log_unit, "(T2,A,T72,A9)") "KERNEL| Spin flip", "Collinear"
893 ! Spin-flip excitations with noncollinear exchange-correlation kernel requested
894 ELSE IF (tddfpt_control%spinflip == tddfpt_sf_noncol) THEN
895 WRITE (log_unit, "(T2,A,T69,A12)") "KERNEL| Spin flip", "Noncollinear"
896 END IF
897 ELSE
898 WRITE (log_unit, "(T2,A,T74,A7)") "KERNEL| Spin symmetry of excitations", "Singlet"
899 END IF
900 WRITE (log_unit, "(T2,A,T73,I8)") "TDDFPT| Number of states calculated", tddfpt_control%nstates
901 WRITE (log_unit, "(T2,A,T73,I8)") "TDDFPT| Number of Davidson iterations", tddfpt_control%niters
902 WRITE (log_unit, "(T2,A,T66,E15.3)") "TDDFPT| Davidson iteration convergence", tddfpt_control%conv
903 WRITE (log_unit, "(T2,A,T73,I8)") "TDDFPT| Max. number of Krylov space vectors", tddfpt_control%nkvs
904 END IF
905
906 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 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 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 CALL timeset(routinen, handle)
976
977 NULLIFY (admm_env, matrix_ks_oep)
978 do_admm = tddfpt_control%do_admm
979 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 IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
983
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 wfm_rho_orb=work_matrices%rho_ao_orb_fm_sub)
990
991 IF (do_admm) THEN
992 ! Full kernel with ADMM
993 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 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 sub_env=sub_env)
1005 END IF
1006
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 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 sub_env=sub_env)
1021 kernel_env%full_kernel => full_kernel_env
1022 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 sub_env=sub_env)
1030 kernel_env%full_kernel => full_kernel_env
1031 NULLIFY (kernel_env%admm_kernel)
1032 END IF
1033 ! Fxc from kernel definition
1034 do_exck = tddfpt_control%do_exck
1035 kernel_env%full_kernel%do_exck = do_exck
1036 ! initilize xc kernel
1037 IF (do_exck) THEN
1038 CALL create_fxc_kernel(work_matrices%rho_orb_struct_sub, work_matrices%fxc_rspace_sub, &
1039 xc_section, tddfpt_control%rks_triplets, sub_env, qs_env)
1040 END IF
1041 END IF
1042
1043 ! lri input
1044 IF (tddfpt_control%do_lrigpw) THEN
1045 lri_section => section_vals_get_subs_vals(tddfpt_section, "LRIGPW")
1046 CALL tddfpt2_lri_init(qs_env, kernel_env, lri_section, &
1047 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 CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=nocc)
1054 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1055 CALL cp_fm_get_info(gs_mos(1)%mos_virt, ncol_global=nvirt)
1056 ELSE
1057 CALL cp_fm_get_info(gs_mos(2)%mos_virt, ncol_global=nvirt)
1058 END IF
1059 nstate_max = nocc*nvirt
1060 IF ((SIZE(gs_mos, 1) == 2) .AND. (tddfpt_control%spinflip == no_sf_tddfpt)) THEN
1061 CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nocc)
1062 CALL cp_fm_get_info(gs_mos(2)%mos_virt, ncol_global=nvirt)
1063 nstate_max = nocc*nvirt + nstate_max
1064 END IF
1065 IF (nstates > nstate_max) THEN
1066 cpwarn("NUMBER OF EXCITED STATES COULD LEAD TO PROBLEMS!")
1067 cpwarn("Experimental: CHANGED NSTATES TO ITS MAXIMUM VALUE!")
1068 nstates = nstate_max
1069 END IF
1070
1071 soc_section => section_vals_get_subs_vals(tddfpt_section, "SOC")
1072 CALL section_vals_get(soc_section, explicit=do_soc)
1073
1074 ! reuse Ritz vectors from the previous calculation if available
1075 IF (tddfpt_control%is_restart .AND. .NOT. do_soc) THEN
1076 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 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 "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 qs_env=qs_env, nspins=nspins)
1099 CALL cp_print_key_finished_output(log_unit, logger, &
1100 tddfpt_print_section, "GUESS_VECTORS")
1101
1102 CALL tddfpt_orthogonalize_psi1_psi0(evects, work_matrices%S_C0_C0T, qs_env, &
1103 gs_mos, evals, tddfpt_control, work_matrices%S_C0)
1104 CALL tddfpt_orthonormalize_psi1_psi1(evects, nstates, s_evects, matrix_s(1)%matrix)
1105
1106 niters = tddfpt_control%niters
1107 IF (niters > 0) THEN
1108 log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, &
1109 "ITERATION_INFO", extension=".tddfptLog")
1110 energy_unit = cp_print_key_unit_nr(logger, &
1111 tddfpt_print_section, &
1112 "DETAILED_ENERGY", &
1113 extension=".tddfptLog")
1114
1115 IF (log_unit > 0) THEN
1116 WRITE (log_unit, "(1X,A)") "", &
1117 "-------------------------------------------------------------------------------", &
1118 "- TDDFPT WAVEFUNCTION OPTIMIZATION -", &
1119 "-------------------------------------------------------------------------------"
1120
1121 WRITE (log_unit, '(/,T11,A,T27,A,T40,A,T62,A)') "Step", "Time", "Convergence", "Conv. states"
1122 WRITE (log_unit, '(1X,79("-"))')
1123 END IF
1124
1125 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 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 CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter)
1151 ! terminate the loop if either (a) or (b) is true ...
1152 IF ((conv <= tddfpt_control%conv) .OR. iter >= niters) EXIT
1153
1154 ! ... otherwise restart Davidson iterations
1155 evals = 0.0_dp
1156 IF (log_unit > 0) THEN
1157 WRITE (log_unit, '(1X,25("-"),1X,A,1X,25("-"))') "Restart Davidson iterations"
1158 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 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 tddfpt_print_section=tddfpt_print_section)
1169
1170 CALL cp_rm_iter_level(logger%iter_info, "TDDFT_SCF")
1171
1172 ! print convergence summary
1173 IF (log_unit > 0) THEN
1174 CALL integer_to_string(iter, nstates_str)
1175 IF (conv <= tddfpt_control%conv) THEN
1176 WRITE (log_unit, "(1X,A)") "", &
1177 "-------------------------------------------------------------------------------", &
1178 "- TDDFPT run converged in "//trim(nstates_str)//" iteration(s) ", &
1179 "-------------------------------------------------------------------------------"
1180 ELSE
1181 WRITE (log_unit, "(1X,A)") "", &
1182 "-------------------------------------------------------------------------------", &
1183 "- TDDFPT run did NOT converge after "//trim(nstates_str)//" iteration(s) ", &
1184 "-------------------------------------------------------------------------------"
1185 END IF
1186 END IF
1187
1188 CALL cp_print_key_finished_output(energy_unit, logger, &
1189 tddfpt_print_section, "DETAILED_ENERGY")
1190 CALL cp_print_key_finished_output(log_unit, logger, &
1191 tddfpt_print_section, "ITERATION_INFO")
1192 ELSE
1193 CALL cp_warn(__location__, &
1194 "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 log_unit = cp_logger_get_default_io_unit(logger)
1208
1209 namd_print_section => section_vals_get_subs_vals( &
1210 tddfpt_print_section, &
1211 "NAMD_PRINT")
1212 CALL section_vals_get(namd_print_section, explicit=explicit)
1213 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 sub_env)
1222 END IF
1223 ALLOCATE (ostrength(nstates))
1224 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 tddfpt_control%dipole_form)
1234 log_unit, &
1235 evects, &
1236 evals, &
1237 gs_mos, &
1238 matrix_s(1)%matrix, &
1239 tddfpt_control%spinflip, &
1240 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 tddfpt_print_section)
1247 IF (tddfpt_control%do_exciton_descriptors) THEN
1249 log_unit, &
1250 evects, &
1251 gs_mos, &
1252 matrix_s(1)%matrix, &
1253 tddfpt_control%do_directional_exciton_descriptors, &
1254 qs_env)
1255 END IF
1256
1257 IF (tddfpt_control%do_lrigpw) THEN
1258 CALL lri_print_stat(qs_env, &
1259 ltddfpt=.true., &
1260 tddfpt_lri_env=kernel_env%full_kernel%lri_env)
1261 END IF
1262
1263 CALL timestop(handle)
1264 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_mult, ostrength_mult
1329 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: evects_mult
1330
1331 CALL timeset(routinen, handle)
1332
1333 log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, &
1334 "PROGRAM_BANNER", &
1335 extension=".tddfptLog")
1336 CALL tddfpt_soc_header(log_unit)
1337
1338 nspins = SIZE(gs_mos)
1339 ALLOCATE (evects_mult(nspins, nstates))
1340 ALLOCATE (evals_mult(nstates))
1341 do_sf = tddfpt_control%spinflip /= no_sf_tddfpt
1342
1343 ! First multiplicity
1344 IF (lmult_tmp) THEN
1345 IF (log_unit > 0) THEN
1346 WRITE (log_unit, "(1X,A)") "", &
1347 "-------------------------------------------------------------------------------", &
1348 "- TDDFPT SINGLET ENERGIES -", &
1349 "-------------------------------------------------------------------------------"
1350 END IF
1351 mult = 1
1352 ELSE
1353 IF (log_unit > 0) THEN
1354 WRITE (log_unit, "(1X,A)") "", &
1355 "-------------------------------------------------------------------------------", &
1356 "- TDDFPT TRIPLET ENERGIES -", &
1357 "-------------------------------------------------------------------------------"
1358 END IF
1359 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 kernel_env_admm_aux)
1368
1369 ! Clean up in between for full kernel
1370 IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
1371 IF (tddfpt_control%do_admm) CALL release_kernel_env(kernel_env%admm_kernel)
1372 CALL release_kernel_env(kernel_env%full_kernel)
1373 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 tddfpt_control%do_exck, do_sf, qs_env, sub_env)
1378 END IF
1379
1380 DO istate = 1, nstates
1381 DO ispin = 1, nspins
1382 CALL cp_fm_release(s_evects(ispin, istate))
1383 END DO
1384 END DO
1385
1386 DO istate = 1, nstates
1387 DO ispin = 1, nspins
1388 CALL fm_pool_create_fm( &
1389 work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
1390 s_evects(ispin, istate))
1391 END DO
1392 END DO
1393
1394 tddfpt_control%rks_triplets = lmult_tmp
1395
1396 ! Second multiplicity
1397 IF (lmult_tmp) THEN
1398 IF (log_unit > 0) THEN
1399 WRITE (log_unit, "(1X,A)") "", &
1400 " singlet excitations finished ", &
1401 " ", &
1402 "-------------------------------------------------------------------------------", &
1403 "- TDDFPT TRIPLET ENERGIES -", &
1404 "-------------------------------------------------------------------------------"
1405 END IF !log_unit
1406 mult = 3
1407 ELSE
1408 IF (log_unit > 0) THEN
1409 WRITE (log_unit, "(1X,A)") "", &
1410 " triplet excitations finished ", &
1411 " ", &
1412 "-------------------------------------------------------------------------------", &
1413 "- TDDFPT SINGLET ENERGIES -", &
1414 "-------------------------------------------------------------------------------"
1415 END IF !log_unit
1416 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 kernel_env_admm_aux)
1425
1426 ! Compute perturbative SOC correction
1427 ! Order should always be singlet triplet in tddfpt_soc
1428 IF (lmult_tmp) THEN
1429 CALL tddfpt_soc(qs_env, evals_mult, evals, evects_mult, evects, gs_mos) !mult=singlet
1430 ELSE
1431 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 DO ispin = 1, SIZE(evects_mult, 1)
1436 DO istate = 1, SIZE(evects_mult, 2)
1437 CALL cp_fm_release(evects_mult(ispin, istate))
1438 END DO
1439 END DO
1440 DEALLOCATE (evects_mult, evals_mult, ostrength_mult)
1441
1442 CALL timestop(handle)
1443
1444 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 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 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: orblist
1468 INTEGER, DIMENSION(:), POINTER :: mollist
1469 LOGICAL :: do_res, do_sf, ew1, ew2, ew3, ewcut, lms
1470 REAL(kind=dp) :: eclow, ecup, eint, emo
1471 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 CALL timeset(routinen, handle)
1477
1478 res_section => section_vals_get_subs_vals(tddfpt_section, "REDUCED_EXCITATION_SPACE")
1479 CALL section_vals_val_get(res_section, "_SECTION_PARAMETERS_", l_val=do_res)
1480
1481 ! spin flip TDA
1482 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1483 do_sf = .false.
1484 ELSE
1485 do_sf = .true.
1486 END IF
1487
1488 nspins = SIZE(gs_mos)
1489 IF (.NOT. do_res) THEN
1490 DO ispin = 1, nspins
1491 nmo = gs_mos(ispin)%nmo_occ
1492 tddfpt_control%nactive(ispin) = nmo
1493 gs_mos(ispin)%nmo_active = nmo
1494 ALLOCATE (gs_mos(ispin)%index_active(nmo))
1495 DO i = 1, nmo
1496 gs_mos(ispin)%index_active(i) = i
1497 END DO
1498 END DO
1499 ELSE
1500 IF (iounit > 0) THEN
1501 WRITE (iounit, "(/,1X,27('='),A,26('='))") ' REDUCED EXCITATION SPACE '
1502 END IF
1503 CALL section_vals_val_get(res_section, "ENERGY_WINDOW", explicit=ew1)
1504 CALL section_vals_val_get(res_section, "UPPER_ENERGY_CUTOFF", explicit=ew2)
1505 CALL section_vals_val_get(res_section, "LOWER_ENERGY_CUTOFF", explicit=ew3)
1506 ewcut = (ew1 .OR. ew2 .OR. ew3)
1507 CALL section_vals_val_get(res_section, "MOLECULE_LIST", explicit=lms)
1508
1509 CALL section_vals_val_get(res_section, "ENERGY_WINDOW", r_vals=rvint)
1510 cpassert(SIZE(rvint) == 2)
1511 eclow = rvint(1)
1512 ecup = rvint(2)
1513 CALL section_vals_val_get(res_section, "UPPER_ENERGY_CUTOFF", r_val=eint)
1514 ecup = min(ecup, eint)
1515 CALL section_vals_val_get(res_section, "LOWER_ENERGY_CUTOFF", r_val=eint)
1516 eclow = max(eclow, eint)
1517 IF (ewcut .AND. (iounit > 0)) THEN
1518 IF (eclow < -1.e8_dp .AND. ecup > 1.e8_dp) THEN
1519 WRITE (iounit, "(1X,A,T51,A10,T71,A10)") &
1520 '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 'Orbital Energy Window [eV]', " -Inf", evolt*ecup
1524 ELSE IF (ecup > 1.e8_dp) THEN
1525 WRITE (iounit, "(1X,A,T51,F10.4,T71,A10)") &
1526 'Orbital Energy Window [eV]', evolt*eclow, " Inf"
1527 ELSE
1528 WRITE (iounit, "(1X,A,T51,F10.4,T71,F10.4)") &
1529 'Orbital Energy Window [eV]', evolt*eclow, evolt*ecup
1530 END IF
1531 END IF
1532
1533 nmol = 0
1534 IF (lms) THEN
1535 CALL section_vals_val_get(res_section, "MOLECULE_LIST", i_vals=mollist)
1536 nmol = SIZE(mollist)
1537 WRITE (iounit, "(1X,A)") 'List of Selected Molecules'
1538 WRITE (iounit, "(1X,15(I5))") mollist(1:nmol)
1539 END IF
1540
1541 DO ispin = 1, nspins
1542 tddfpt_control%nactive(ispin) = gs_mos(ispin)%nmo_occ
1543 END DO
1544 nmo = maxval(tddfpt_control%nactive)
1545 ALLOCATE (orblist(nmo, nspins))
1546 orblist = 0
1547
1548 IF (lms) THEN
1549 ! ignore for now
1550 orblist = 1
1551 DO ispin = 1, nspins
1552 cpassert(.NOT. ASSOCIATED(gs_mos(ispin)%evals_occ_matrix))
1553 END DO
1554 ELSEIF (ewcut) THEN
1555 ! Filter orbitals wrt energy window
1556 DO ispin = 1, nspins
1557 DO i = 1, gs_mos(ispin)%nmo_occ
1558 emo = gs_mos(ispin)%evals_occ(i)
1559 IF (emo > eclow .AND. emo < ecup) orblist(i, ispin) = 1
1560 END DO
1561 END DO
1562 ELSE
1563 orblist = 1
1564 END IF
1565
1566 ! count active orbitals
1567 nmo = sum(orblist)
1568 IF (nmo == 0) THEN
1569 cpabort("RSE TDA: no active orbitals selected.")
1570 END IF
1571 DO ispin = 1, nspins
1572 nmo = sum(orblist(:, ispin))
1573 tddfpt_control%nactive(ispin) = nmo
1574 gs_mos(ispin)%nmo_active = nmo
1575 ALLOCATE (gs_mos(ispin)%index_active(nmo))
1576 io = 0
1577 DO i = 1, SIZE(orblist, 1)
1578 IF (orblist(i, ispin) == 1) THEN
1579 io = io + 1
1580 gs_mos(ispin)%index_active(io) = i
1581 END IF
1582 END DO
1583 END DO
1584 DEALLOCATE (orblist)
1585
1586 IF (lms) THEN
1587 ! output information
1588 ELSE
1589 IF (iounit > 0) THEN
1590 WRITE (iounit, "(1X,A)") 'List of Selected States'
1591 IF (nspins == 1) THEN
1592 WRITE (iounit, "(A,T67,A)") ' Active State Orbital', 'Orbital Energy'
1593 DO i = 1, gs_mos(1)%nmo_active
1594 io = gs_mos(1)%index_active(i)
1595 WRITE (iounit, "(T8,I6,T21,I6,T61,F20.4)") i, io, gs_mos(1)%evals_occ(io)*evolt
1596 END DO
1597 ELSE
1598 DO ispin = 1, nspins
1599 WRITE (iounit, "(1X,A,I2)") 'Spin ', ispin
1600 WRITE (iounit, "(A,T67,A)") ' Active State Orbital', 'Orbital Energy'
1601 DO i = 1, gs_mos(ispin)%nmo_active
1602 io = gs_mos(ispin)%index_active(i)
1603 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 IF (do_sf) THEN
1611 cpabort("Restricted Active Space with spin flip TDA NYA")
1612 END IF
1613
1614 IF (iounit > 0) THEN
1615 WRITE (iounit, "(1X,79('='))")
1616 END IF
1617 END IF
1618
1619 ! Allocate mos_active
1620 DO ispin = 1, nspins
1621 CALL get_qs_env(qs_env, blacs_env=blacs_env)
1622 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, nrow_global=nao)
1623 nmo = gs_mos(ispin)%nmo_active
1624 CALL cp_fm_struct_create(fm_struct, template_fmstruct=gs_mos(ispin)%mos_occ%matrix_struct, &
1625 ncol_global=nmo, context=blacs_env)
1626 NULLIFY (gs_mos(ispin)%mos_active)
1627 ALLOCATE (gs_mos(ispin)%mos_active)
1628 CALL cp_fm_create(gs_mos(ispin)%mos_active, fm_struct)
1629 CALL cp_fm_struct_release(fm_struct)
1630 ! copy the active orbitals
1631 IF (gs_mos(ispin)%nmo_active == gs_mos(ispin)%nmo_occ) THEN
1632 DO i = 1, gs_mos(ispin)%nmo_active
1633 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 nao, nmo, 1, 1, 1, 1)
1637 ELSE
1638 DO i = 1, gs_mos(ispin)%nmo_active
1639 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 nao, 1, 1, io, 1, i)
1642 END DO
1643 END IF
1644 END DO
1645
1646 CALL timestop(handle)
1647
1648 END SUBROUTINE init_res_method
1649
1650END MODULE qs_tddfpt2_methods
Contains ADMM methods which require molecular orbitals.
subroutine, public admm_fit_mo_coeffs(admm_env, matrix_s_aux_fit, matrix_s_mixed, mos, mos_aux_fit, geometry_did_change)
...
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition admm_types.F:593
Define the atomic kind types and their sub types.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public grimme2016
integer, save, public iannuzzi2005
integer, save, public hernandez2025
integer, save, public grimme2013
Handles all functions related to the CELL.
Definition cell_types.F:15
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
pool for for elements that are retained and released
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
subroutine, public cp_rm_iter_level(iteration_info, level_name, n_rlevel_att)
Removes an iteration level.
subroutine, public cp_add_iter_level(iteration_info, level_name, n_rlevel_new)
Adds an iteration level.
Types for excited states potential energies.
subroutine, public tddfpt_header(iw)
...
Definition header.F:122
subroutine, public tddfpt_soc_header(iw)
...
Definition header.F:143
Utilities for hfx and admm methods.
subroutine, public aux_admm_init(qs_env, mos, admm_env, admm_control, basis_type)
Minimal setup routine for admm_env No forces No k-points No DFT correction terms.
Types and set/get functions for HFX.
Definition hfx_types.F:15
subroutine, public hfx_create(x_data, para_env, hfx_section, atomic_kind_set, qs_kind_set, particle_set, dft_control, cell, orb_basis, ri_basis, nelectron_total, nkp_grid)
This routine allocates and initializes all types in hfx_data
Definition hfx_types.F:595
subroutine, public compare_hfx_sections(hfx_section1, hfx_section2, is_identical, same_except_frac)
Compares the non-technical parts of two HFX input section and check whether they are the same Ignore ...
Definition hfx_types.F:2957
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_admm_purify_none
integer, parameter, public tddfpt_sf_col
integer, parameter, public tddfpt_kernel_none
integer, parameter, public do_admm_basis_projection
integer, parameter, public do_admm_aux_exch_func_none
integer, parameter, public tddfpt_dipole_velocity
integer, parameter, public do_potential_truncated
integer, parameter, public tddfpt_kernel_full
integer, parameter, public tddfpt_sf_noncol
integer, parameter, public tddfpt_kernel_stda
integer, parameter, public no_sf_tddfpt
integer, parameter, public do_admm_exch_scaling_none
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Calculates integral matrices for LRIGPW method lri : local resolution of the identity.
subroutine, public lri_print_stat(qs_env, ltddfpt, tddfpt_lri_env)
...
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
subroutine, public lri_density_release(lri_density)
releases the given lri_density
subroutine, public lri_env_release(lri_env)
releases the given lri_env
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:136
Interface to the message passing library MPI.
generate or use from input minimal basis set
subroutine, public create_minbas_set(qs_env, unit_nr, basis_type, primitive)
...
Define the data structure for the particle information.
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public evolt
Definition physcon.F:183
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
subroutine, public create_fxc_kernel(rho_struct, fxc_rspace, xc_section, is_rks_triplets, sub_env, qs_env)
Create the xc kernel potential for the approximate Fxc kernel model.
subroutine, public create_kernel_env(kernel_env, xc_section, is_rks_triplets, rho_struct_sub, lsd_singlets, do_excitations, sub_env, qs_env)
Create kernel environment.
subroutine, public release_kernel_env(kernel_env)
Release kernel environment.
Define the quickstep kind type and their sub types.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
subroutine, public eigensolver(matrix_ks_fm, mo_set, ortho, work, cholesky_method, do_level_shift, level_shift, matrix_u_fm, use_jacobi)
Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- vectors and MO eigenvalues....
module that contains the definitions of the scf types
subroutine, public assign_state(qs_env, matrix_s, evects, psi0, wfn_history, my_state)
...
subroutine, public tddfpt_construct_aux_fit_density(rho_orb_struct, rho_aux_fit_struct, local_rho_set, qs_env, sub_env, wfm_rho_orb, wfm_rho_aux_fit, wfm_aux_orb)
Project a charge density expressed in primary basis set into the auxiliary basis set.
subroutine, public tddfpt_construct_ground_state_orb_density(rho_orb_struct, rho_xc_struct, is_rks_triplets, qs_env, sub_env, wfm_rho_orb)
Compute the ground-state charge density expressed in primary basis set.
subroutine, public tddfpt_orthonormalize_psi1_psi1(evects, nvects_new, s_evects, matrix_s)
Make new TDDFPT trial vectors orthonormal to all previous trial vectors.
subroutine, public tddfpt_orthogonalize_psi1_psi0(evects, s_c0_c0t, qs_env, gs_mos, evals, tddfpt_control, s_c0)
Make TDDFPT trial vectors orthogonal to all occupied molecular orbitals.
real(kind=dp) function, public tddfpt_davidson_solver(evects, evals, s_evects, gs_mos, tddfpt_control, matrix_ks, qs_env, kernel_env, sub_env, logger, iter_unit, energy_unit, tddfpt_print_section, work_matrices)
Perform Davidson iterations.
subroutine, public tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, sub_env, work_matrices)
Perform TDDFPT gradient calculation. This routine calculates the response vector R of Eq....
subroutine, public tddfpt_print_forces(qs_env, evects, evals, ostrength, print_section, gs_mos, kernel_env, sub_env, work_matrices)
Calculate and print forces of selected excited states.
subroutine, public tddfpt2_lri_init(qs_env, kernel_env, lri_section, tddfpt_print_section)
Initialize LRI environment, basis, neighborlists and matrices.
subroutine, public tddfpt(qs_env, calc_forces, rixs_env)
Perform TDDFPT calculation. If calc_forces then it also builds the response vector for the Z-vector m...
subroutine, public tddfpt_input(qs_env, do_hfx, do_admm, do_exck, do_hfxsr, do_hfxlr, xc_section, tddfpt_print_section, lri_section, hfxsr_section)
TDDFPT input.
subroutine, public tddfpt_energies(qs_env, nstates, nspins, work_matrices, tddfpt_control, logger, tddfpt_print_section, evects, evals, gs_mos, tddfpt_section, s_evects, matrix_s, kernel_env, matrix_ks, sub_env, ostrength, dipole_op_mos_occ, mult, xc_section, full_kernel_env, kernel_env_admm_aux)
The energy calculation has been moved to its own subroutine.
subroutine, public tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_env)
Compute the action of the dipole operator on the ground state wave function.
subroutine, public tddfpt_print_nto_analysis(qs_env, evects, evals, ostrength, gs_mos, matrix_s, print_section)
Print natural transition orbital analysis.
subroutine, public tddfpt_print_excitation_analysis(log_unit, evects, evals, gs_mos, matrix_s, spinflip, min_amplitude)
Print excitation analysis.
subroutine, public tddfpt_print_exciton_descriptors(log_unit, evects, gs_mos, matrix_s, do_directional_exciton_descriptors, qs_env)
Print exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018)
subroutine, public tddfpt_print_summary(log_unit, evects, evals, gs_mos, ostrength, mult, dipole_op_mos_occ, dipole_form)
Print final TDDFPT excitation energies and oscillator strengths.
subroutine, public tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, matrix_s, s_evects, sub_env)
Write Ritz vectors to a binary restart file.
subroutine, public tddfpt_write_restart(evects, evals, gs_mos, logger, tddfpt_print_section)
Write Ritz vectors to a binary restart file.
integer function, public tddfpt_read_restart(evects, evals, gs_mos, logger, tddfpt_section, tddfpt_print_section, fm_pool_ao_mo_active, blacs_env_global)
Initialise initial guess vectors by reading (un-normalised) Ritz vectors from a binary restart file.
subroutine, public tddfpt_smeared_occupation(qs_env, gs_mos, log_unit)
...
subroutine, public tddfpt_soc(qs_env, evals_a, evals_b, evects_a, evects_b, gs_mos)
Perform TDDFPT-SOC calculation.
Simplified Tamm Dancoff approach (sTDA).
subroutine, public stda_init_param(qs_env, stda_kernel, stda_control)
Get the parameters needed for an sTDA calculation.
subroutine, public deallocate_stda_env(stda_kernel)
Deallocate the sTDA environment.
subroutine, public allocate_stda_env(qs_env, stda_kernel, n_ao, nactive)
Allocate the sTDA environment.
Simplified Tamm Dancoff approach (sTDA).
subroutine, public get_lowdin_mo_coefficients(qs_env, sub_env, work)
Calculate Lowdin MO coefficients.
subroutine, public stda_init_matrices(qs_env, stda_kernel, sub_env, work, tddfpt_control)
Calculate sTDA matrices.
subroutine, public tddfpt_sub_env_init(sub_env, qs_env, mos_occ, mos_active, kernel)
Split MPI communicator to create a set of parallel (sub)groups.
subroutine, public tddfpt_sub_env_release(sub_env)
Release parallel group environment.
subroutine, public hfxsr_create_work_matrices(work_matrices, qs_env, admm_env)
Allocate work matrices for hfxsr.
subroutine, public tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, do_hfx, do_admm, do_hfxlr, do_exck, do_sf, qs_env, sub_env)
Allocate work matrices for full kernel.
subroutine, public tddfpt_release_work_matrices(work_matrices, sub_env)
Release work matrices.
subroutine, public stda_create_work_matrices(work_matrices, gs_mos, nstates, qs_env, sub_env)
Allocate work matrices for sTDA kernel.
subroutine, public tddfpt_release_ground_state_mos(gs_mos)
Release molecular orbitals.
subroutine, public tddfpt_oecorr(qs_env, gs_mos, matrix_ks_oep)
Callculate orbital corrected KS matrix for TDDFPT.
subroutine, public tddfpt_guess_vectors(evects, evals, gs_mos, log_unit, tddfpt_control, fm_pool_ao_mo_active, qs_env, nspins)
Generate missed guess vectors.
subroutine, public tddfpt_init_mos(qs_env, gs_mos, iounit)
Prepare MOs for TDDFPT Calculations.
Define Resonant Inelastic XRAY Scattering (RIXS) control type and associated create,...
Definition rixs_types.F:13
Utilities for string manipulations.
subroutine, public integer_to_string(inumber, string)
Converts an integer number to a string. The WRITE statement will return an error message,...
Writes information on XC functionals to output.
subroutine, public xc_write(iounit, xc_section, lsd)
...
stores some data used in wavefunction fitting
Definition admm_types.F:120
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Contains information on the excited states energy.
stores all the informations relevant to an mpi environment
Collection of variables required to evaluate adiabatic TDDFPT kernel.
Type to hold environments for the different kernels.
Provides all information about a quickstep kind.
Ground state molecular orbitals.
Set of temporary ("work") matrices.
Valence state coming from the qs_tddfpt2 routine.
Definition rixs_types.F:40