(git:64e569e)
Loading...
Searching...
No Matches
qs_tddfpt2_properties.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
10 USE bibliography, ONLY: martin2003,&
11 cite_reference
16 USE cell_types, ONLY: cell_type
19 USE cp_cfm_types, ONLY: cp_cfm_create,&
27 USE cp_dbcsr_api, ONLY: &
44 USE cp_fm_types, ONLY: cp_fm_create,&
55 USE cp_output_handling, ONLY: cp_p_file,&
60 USE input_constants, ONLY: no_sf_tddfpt,&
69 USE kinds, ONLY: default_path_length,&
70 dp,&
71 int_8
72 USE mathconstants, ONLY: twopi,&
73 z_one,&
74 z_zero
75 USE message_passing, ONLY: mp_comm_type,&
83 USE physcon, ONLY: evolt
84 USE pw_env_types, ONLY: pw_env_get,&
87 USE pw_pool_types, ONLY: pw_pool_p_type,&
89 USE pw_types, ONLY: pw_c1d_gs_type,&
96 USE qs_mo_types, ONLY: allocate_mo_set,&
104 USE qs_operators_ao, ONLY: rrc_xyz_ao
106 USE qs_subsys_types, ONLY: qs_subsys_get,&
110 USE util, ONLY: sort
111#include "./base/base_uses.f90"
112
113 IMPLICIT NONE
114
115 PRIVATE
116
117 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_properties'
118
119 ! number of first derivative components (3: d/dx, d/dy, d/dz)
120 INTEGER, PARAMETER, PRIVATE :: nderivs = 3
121 INTEGER, PARAMETER, PRIVATE :: maxspins = 2
122
125
126! **************************************************************************************************
127
128CONTAINS
129
130! **************************************************************************************************
131!> \brief Compute the action of the dipole operator on the ground state wave function.
132!> \param dipole_op_mos_occ 2-D array [x,y,z ; spin] of matrices where to put the computed quantity
133!> (allocated and initialised on exit)
134!> \param tddfpt_control TDDFPT control parameters
135!> \param gs_mos molecular orbitals optimised for the ground state
136!> \param qs_env Quickstep environment
137!> \par History
138!> * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov]
139!> * 06.2018 dipole operator based on the Berry-phase formula [Sergey Chulkov]
140!> * 08.2018 splited of from 'tddfpt_print_summary' and merged with code from 'tddfpt'
141!> [Sergey Chulkov]
142!> \note \parblock
143!> Adapted version of the subroutine find_contributions() which was originally created
144!> by Thomas Chassaing on 02.2005.
145!>
146!> The relation between dipole integrals in velocity and length forms are the following:
147!> \f[<\psi_i|\nabla|\psi_a> = <\psi_i|\vec{r}|\hat{H}\psi_a> - <\hat{H}\psi_i|\vec{r}|\psi_a>
148!> = (\epsilon_a - \epsilon_i) <\psi_i|\vec{r}|\psi_a> .\f],
149!> due to the commutation identity:
150!> \f[\vec{r}\hat{H} - \hat{H}\vec{r} = [\vec{r},\hat{H}] = [\vec{r},-1/2 \nabla^2] = \nabla\f] .
151!> \endparblock
152! **************************************************************************************************
153 SUBROUTINE tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_env)
154 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :), &
155 INTENT(inout) :: dipole_op_mos_occ
156 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
157 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
158 INTENT(in) :: gs_mos
159 TYPE(qs_environment_type), POINTER :: qs_env
160
161 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_dipole_operator'
162
163 INTEGER :: handle, i_cos_sin, icol, ideriv, irow, &
164 ispin, jderiv, nao, ncols_local, &
165 ndim_periodic, nrows_local, nspins
166 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
167 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt
168 REAL(kind=dp) :: eval_occ
169 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
170 POINTER :: local_data_ediff, local_data_wfm
171 REAL(kind=dp), DIMENSION(3) :: kvec, reference_point
172 TYPE(cell_type), POINTER :: cell
173 TYPE(cp_blacs_env_type), POINTER :: blacs_env
174 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: gamma_00, gamma_inv_00
175 TYPE(cp_fm_struct_type), POINTER :: fm_struct
176 TYPE(cp_fm_type) :: ediff_inv, wfm_ao_ao, wfm_mo_virt_mo_occ
177 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: s_mos_virt
178 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: dberry_mos_occ, gamma_real_imag, opvec
179 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: berry_cossin_xyz, matrix_s, rrc_xyz, scrm
180 TYPE(dft_control_type), POINTER :: dft_control
181 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
182 POINTER :: sab_orb
183 TYPE(pw_env_type), POINTER :: pw_env
184 TYPE(pw_poisson_type), POINTER :: poisson_env
185 TYPE(qs_ks_env_type), POINTER :: ks_env
186
187 CALL timeset(routinen, handle)
188
189 NULLIFY (blacs_env, cell, matrix_s, pw_env)
190 CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, matrix_s=matrix_s, pw_env=pw_env)
191
192 nspins = SIZE(gs_mos)
193 CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
194 DO ispin = 1, nspins
195 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
196 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
197 END DO
198
199 ! +++ allocate dipole operator matrices (must be deallocated elsewhere)
200 ALLOCATE (dipole_op_mos_occ(nderivs, nspins))
201 DO ispin = 1, nspins
202 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct)
203
204 DO ideriv = 1, nderivs
205 CALL cp_fm_create(dipole_op_mos_occ(ideriv, ispin), fm_struct)
206 END DO
207 END DO
208
209 ! +++ allocate work matrices
210 ALLOCATE (s_mos_virt(nspins))
211 DO ispin = 1, nspins
212 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct)
213 CALL cp_fm_create(s_mos_virt(ispin), fm_struct)
214 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
215 gs_mos(ispin)%mos_virt, &
216 s_mos_virt(ispin), &
217 ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
218 END DO
219
220 ! check that the chosen dipole operator is consistent with the periodic boundary conditions used
221 CALL pw_env_get(pw_env, poisson_env=poisson_env)
222 ndim_periodic = count(poisson_env%parameters%periodic == 1)
223
224 ! select default for dipole form
225 IF (tddfpt_control%dipole_form == 0) THEN
226 CALL get_qs_env(qs_env, dft_control=dft_control)
227 IF (dft_control%qs_control%xtb) THEN
228 IF (ndim_periodic == 0) THEN
229 tddfpt_control%dipole_form = tddfpt_dipole_length
230 ELSE
231 tddfpt_control%dipole_form = tddfpt_dipole_velocity
232 END IF
233 ELSE
234 tddfpt_control%dipole_form = tddfpt_dipole_velocity
235 END IF
236 END IF
237
238 SELECT CASE (tddfpt_control%dipole_form)
240 IF (ndim_periodic /= 3) THEN
241 CALL cp_warn(__location__, &
242 "Fully periodic Poisson solver (PERIODIC xyz) "// &
243 "or a large supercell in non-periodic directions is needed "// &
244 "for oscillator strengths based on the Berry phase formula")
245 END IF
246
247 NULLIFY (berry_cossin_xyz)
248 ! index: 1 = Re[exp(-i * G_t * t)],
249 ! 2 = Im[exp(-i * G_t * t)];
250 ! t = x,y,z
251 CALL dbcsr_allocate_matrix_set(berry_cossin_xyz, 2)
252
253 DO i_cos_sin = 1, 2
254 CALL dbcsr_init_p(berry_cossin_xyz(i_cos_sin)%matrix)
255 CALL dbcsr_copy(berry_cossin_xyz(i_cos_sin)%matrix, matrix_s(1)%matrix)
256 END DO
257
258 ! +++ allocate berry-phase-related work matrices
259 ALLOCATE (gamma_00(nspins), gamma_inv_00(nspins), gamma_real_imag(2, nspins), opvec(2, nspins))
260 ALLOCATE (dberry_mos_occ(nderivs, nspins))
261 DO ispin = 1, nspins
262 NULLIFY (fm_struct)
263 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_occ(ispin), &
264 ncol_global=nmo_occ(ispin), context=blacs_env)
265
266 CALL cp_cfm_create(gamma_00(ispin), fm_struct)
267 CALL cp_cfm_create(gamma_inv_00(ispin), fm_struct)
268
269 DO i_cos_sin = 1, 2
270 CALL cp_fm_create(gamma_real_imag(i_cos_sin, ispin), fm_struct)
271 END DO
272 CALL cp_fm_struct_release(fm_struct)
273
274 ! G_real C_0, G_imag C_0
275 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct)
276 DO i_cos_sin = 1, 2
277 CALL cp_fm_create(opvec(i_cos_sin, ispin), fm_struct)
278 END DO
279
280 ! dBerry * C_0
281 DO ideriv = 1, nderivs
282 CALL cp_fm_create(dberry_mos_occ(ideriv, ispin), fm_struct)
283 CALL cp_fm_set_all(dberry_mos_occ(ideriv, ispin), 0.0_dp)
284 END DO
285 END DO
286
287 DO ideriv = 1, nderivs
288 kvec(:) = twopi*cell%h_inv(ideriv, :)
289 DO i_cos_sin = 1, 2
290 CALL dbcsr_set(berry_cossin_xyz(i_cos_sin)%matrix, 0.0_dp)
291 END DO
292 CALL build_berry_moment_matrix(qs_env, berry_cossin_xyz(1)%matrix, &
293 berry_cossin_xyz(2)%matrix, kvec)
294
295 DO ispin = 1, nspins
296 ! i_cos_sin = 1: cos (real) component; opvec(1) = gamma_real C_0
297 ! i_cos_sin = 2: sin (imaginary) component; opvec(2) = gamma_imag C_0
298 DO i_cos_sin = 1, 2
299 CALL cp_dbcsr_sm_fm_multiply(berry_cossin_xyz(i_cos_sin)%matrix, &
300 gs_mos(ispin)%mos_occ, &
301 opvec(i_cos_sin, ispin), &
302 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
303 END DO
304
305 CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, &
306 1.0_dp, gs_mos(ispin)%mos_occ, opvec(1, ispin), &
307 0.0_dp, gamma_real_imag(1, ispin))
308
309 CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, &
310 -1.0_dp, gs_mos(ispin)%mos_occ, opvec(2, ispin), &
311 0.0_dp, gamma_real_imag(2, ispin))
312
313 CALL cp_fm_to_cfm(msourcer=gamma_real_imag(1, ispin), &
314 msourcei=gamma_real_imag(2, ispin), &
315 mtarget=gamma_00(ispin))
316
317 ! gamma_inv_00 = Q = [C_0^T (gamma_real - i gamma_imag) C_0] ^ {-1}
318 CALL cp_cfm_set_all(gamma_inv_00(ispin), z_zero, z_one)
319 CALL cp_cfm_solve(gamma_00(ispin), gamma_inv_00(ispin))
320
321 CALL cp_cfm_to_fm(msource=gamma_inv_00(ispin), &
322 mtargetr=gamma_real_imag(1, ispin), &
323 mtargeti=gamma_real_imag(2, ispin))
324
325 ! dBerry_mos_occ is identical to dBerry_psi0 from qs_linres_op % polar_operators()
326 CALL parallel_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), &
327 1.0_dp, opvec(1, ispin), gamma_real_imag(2, ispin), &
328 0.0_dp, dipole_op_mos_occ(1, ispin))
329 CALL parallel_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), &
330 -1.0_dp, opvec(2, ispin), gamma_real_imag(1, ispin), &
331 1.0_dp, dipole_op_mos_occ(1, ispin))
332
333 DO jderiv = 1, nderivs
334 CALL cp_fm_scale_and_add(1.0_dp, dberry_mos_occ(jderiv, ispin), &
335 cell%hmat(jderiv, ideriv), dipole_op_mos_occ(1, ispin))
336 END DO
337 END DO
338 END DO
339
340 ! --- release berry-phase-related work matrices
341 CALL cp_fm_release(opvec)
342 CALL cp_fm_release(gamma_real_imag)
343 DO ispin = nspins, 1, -1
344 CALL cp_cfm_release(gamma_inv_00(ispin))
345 CALL cp_cfm_release(gamma_00(ispin))
346 END DO
347 DEALLOCATE (gamma_00, gamma_inv_00)
348 CALL dbcsr_deallocate_matrix_set(berry_cossin_xyz)
349
350 ! trans_dipole = 2|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00) +
351 ! 2|e|/|G_mu| * Tr Imag(C_0^T * (gamma_real - i gamma_imag) * evects * gamma_inv_00) ,
352 !
353 ! Taking into account the symmetry of the matrices 'gamma_real' and 'gamma_imag' and the fact
354 ! that the response wave-function is a real-valued function, the above expression can be simplified as
355 ! trans_dipole = 4|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00)
356 !
357 ! 1/|G_mu| = |lattice_vector_mu| / (2*pi) .
358 DO ispin = 1, nspins
359
360 DO ideriv = 1, nderivs
361 CALL cp_fm_to_fm(dberry_mos_occ(ideriv, ispin), dipole_op_mos_occ(ideriv, ispin))
362 END DO
363 END DO
364
365 CALL cp_fm_release(wfm_ao_ao)
366 CALL cp_fm_release(dberry_mos_occ)
367
369 IF (ndim_periodic /= 0) THEN
370 CALL cp_warn(__location__, &
371 "Non-periodic Poisson solver (PERIODIC none) "// &
372 "or a large supercell approach is needed "// &
373 "for oscillator strengths based on the length operator")
374 END IF
375
376 ! compute components of the dipole operator in the length form
377 NULLIFY (rrc_xyz)
378 CALL dbcsr_allocate_matrix_set(rrc_xyz, nderivs)
379
380 DO ideriv = 1, nderivs
381 CALL dbcsr_init_p(rrc_xyz(ideriv)%matrix)
382 CALL dbcsr_copy(rrc_xyz(ideriv)%matrix, matrix_s(1)%matrix)
383 END DO
384
385 CALL get_reference_point(reference_point, qs_env=qs_env, &
386 reference=tddfpt_control%dipole_reference, &
387 ref_point=tddfpt_control%dipole_ref_point)
388
389 CALL rrc_xyz_ao(op=rrc_xyz, qs_env=qs_env, rc=reference_point, order=1, &
390 minimum_image=.false., soft=.false.)
391
392 DO ispin = 1, nspins
393
394 DO ideriv = 1, nderivs
395 CALL cp_dbcsr_sm_fm_multiply(rrc_xyz(ideriv)%matrix, &
396 gs_mos(ispin)%mos_occ, &
397 dipole_op_mos_occ(ideriv, ispin), &
398 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
399 END DO
400
401 END DO
402
403 CALL dbcsr_deallocate_matrix_set(rrc_xyz)
404
406 ! generate overlap derivatives
407 CALL get_qs_env(qs_env, ks_env=ks_env, sab_orb=sab_orb)
408 NULLIFY (scrm)
409 CALL build_overlap_matrix(ks_env, matrix_s=scrm, nderivative=1, &
410 basis_type_a="ORB", basis_type_b="ORB", &
411 sab_nl=sab_orb)
412
413 DO ispin = 1, nspins
414 DO ideriv = 1, nderivs
415 CALL cp_dbcsr_sm_fm_multiply(scrm(ideriv + 1)%matrix, &
416 gs_mos(ispin)%mos_occ, &
417 dipole_op_mos_occ(ideriv, ispin), &
418 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
419 END DO
420
421 CALL cp_fm_release(wfm_mo_virt_mo_occ)
422 END DO
424
426 ! generate overlap derivatives
427 CALL get_qs_env(qs_env, ks_env=ks_env, sab_orb=sab_orb)
428 NULLIFY (scrm)
429 CALL build_overlap_matrix(ks_env, matrix_s=scrm, nderivative=1, &
430 basis_type_a="ORB", basis_type_b="ORB", &
431 sab_nl=sab_orb)
432
433 DO ispin = 1, nspins
434 NULLIFY (fm_struct)
435 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), &
436 ncol_global=nmo_occ(ispin), context=blacs_env)
437 CALL cp_fm_create(ediff_inv, fm_struct)
438 CALL cp_fm_create(wfm_mo_virt_mo_occ, fm_struct)
439 CALL cp_fm_struct_release(fm_struct)
440
441 CALL cp_fm_get_info(ediff_inv, nrow_local=nrows_local, ncol_local=ncols_local, &
442 row_indices=row_indices, col_indices=col_indices, local_data=local_data_ediff)
443 CALL cp_fm_get_info(wfm_mo_virt_mo_occ, local_data=local_data_wfm)
444
445!$OMP PARALLEL DO DEFAULT(NONE), &
446!$OMP PRIVATE(eval_occ, icol, irow), &
447!$OMP SHARED(col_indices, gs_mos, ispin, local_data_ediff, ncols_local, nrows_local, row_indices)
448 DO icol = 1, ncols_local
449 ! E_occ_i ; imo_occ = col_indices(icol)
450 eval_occ = gs_mos(ispin)%evals_occ(col_indices(icol))
451
452 DO irow = 1, nrows_local
453 ! ediff_inv_weights(a, i) = 1.0 / (E_virt_a - E_occ_i)
454 ! imo_virt = row_indices(irow)
455 local_data_ediff(irow, icol) = 1.0_dp/(gs_mos(ispin)%evals_virt(row_indices(irow)) - eval_occ)
456 END DO
457 END DO
458!$OMP END PARALLEL DO
459
460 DO ideriv = 1, nderivs
461 CALL cp_dbcsr_sm_fm_multiply(scrm(ideriv + 1)%matrix, &
462 gs_mos(ispin)%mos_occ, &
463 dipole_op_mos_occ(ideriv, ispin), &
464 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
465
466 CALL parallel_gemm('T', 'N', nmo_virt(ispin), nmo_occ(ispin), nao, &
467 1.0_dp, gs_mos(ispin)%mos_virt, dipole_op_mos_occ(ideriv, ispin), &
468 0.0_dp, wfm_mo_virt_mo_occ)
469
470 ! in-place element-wise (Schur) product;
471 ! avoid allocation of a temporary [nmo_virt x nmo_occ] matrix which is needed
472 ! for cp_fm_schur_product() subroutine call
473
474!$OMP PARALLEL DO DEFAULT(NONE), &
475!$OMP PRIVATE(icol, irow), &
476!$OMP SHARED(ispin, local_data_ediff, local_data_wfm, ncols_local, nrows_local)
477 DO icol = 1, ncols_local
478 DO irow = 1, nrows_local
479 local_data_wfm(irow, icol) = local_data_wfm(irow, icol)*local_data_ediff(irow, icol)
480 END DO
481 END DO
482!$OMP END PARALLEL DO
483
484 CALL parallel_gemm('N', 'N', nao, nmo_occ(ispin), nmo_virt(ispin), &
485 1.0_dp, s_mos_virt(ispin), wfm_mo_virt_mo_occ, &
486 0.0_dp, dipole_op_mos_occ(ideriv, ispin))
487 END DO
488
489 CALL cp_fm_release(wfm_mo_virt_mo_occ)
490 CALL cp_fm_release(ediff_inv)
491 END DO
493
494 CASE DEFAULT
495 cpabort("Unimplemented form of the dipole operator")
496 END SELECT
497
498 ! --- release work matrices
499 CALL cp_fm_release(s_mos_virt)
500
501 CALL timestop(handle)
502 END SUBROUTINE tddfpt_dipole_operator
503
504! **************************************************************************************************
505!> \brief Print final TDDFPT excitation energies and oscillator strengths.
506!> \param log_unit output unit
507!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
508!> SIZE(evects,2) -- number of excited states to print)
509!> \param evals TDDFPT eigenvalues
510!> \param gs_mos ...
511!> \param ostrength TDDFPT oscillator strength
512!> \param mult multiplicity
513!> \param dipole_op_mos_occ action of the dipole operator on the ground state wave function
514!> [x,y,z ; spin]
515!> \param dipole_form ...
516!> \par History
517!> * 05.2016 created [Sergey Chulkov]
518!> * 06.2016 transition dipole moments and oscillator strengths [Sergey Chulkov]
519!> * 07.2016 spin-unpolarised electron density [Sergey Chulkov]
520!> * 08.2018 compute 'dipole_op_mos_occ' in a separate subroutine [Sergey Chulkov]
521!> \note \parblock
522!> Adapted version of the subroutine find_contributions() which was originally created
523!> by Thomas Chassaing on 02.2005.
524!>
525!> Transition dipole moment along direction 'd' is computed as following:
526!> \f[ t_d(spin) = Tr[evects^T dipole\_op\_mos\_occ(d, spin)] .\f]
527!> \endparblock
528! **************************************************************************************************
529 SUBROUTINE tddfpt_print_summary(log_unit, evects, evals, gs_mos, ostrength, mult, &
530 dipole_op_mos_occ, dipole_form)
531 INTEGER, INTENT(in) :: log_unit
532 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
533 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
534 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
535 POINTER :: gs_mos
536 REAL(kind=dp), DIMENSION(:), INTENT(inout) :: ostrength
537 INTEGER, INTENT(in) :: mult
538 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: dipole_op_mos_occ
539 INTEGER, INTENT(in) :: dipole_form
540
541 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_summary'
542
543 CHARACTER(len=1) :: lsd_str
544 CHARACTER(len=20) :: mult_str
545 INTEGER :: handle, i, ideriv, ispin, istate, j, &
546 nactive, nao, nocc, nspins, nstates
547 REAL(kind=dp) :: osc_strength
548 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: trans_dipoles
549 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
550 TYPE(cp_fm_type) :: dipact
551
552 CALL timeset(routinen, handle)
553
554 nspins = SIZE(evects, 1)
555 nstates = SIZE(evects, 2)
556
557 IF (nspins > 1) THEN
558 lsd_str = 'U'
559 ELSE
560 lsd_str = 'R'
561 END IF
562
563 ! *** summary header ***
564 IF (log_unit > 0) THEN
565 CALL integer_to_string(mult, mult_str)
566 WRITE (log_unit, '(/,1X,A1,A,1X,A)') lsd_str, "-TDDFPT states of multiplicity", trim(mult_str)
567 SELECT CASE (dipole_form)
569 WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using Berry operator formulation"
571 WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using length formulation"
573 WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using velocity formulation"
575 WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using old velocity formulation"
577 WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using SCF-MO moment formulation"
578 CASE DEFAULT
579 cpabort("Unimplemented form of the dipole operator")
580 END SELECT
581
582 WRITE (log_unit, '(T10,A,T19,A,T37,A,T69,A)') "State", "Excitation", &
583 "Transition dipole (a.u.)", "Oscillator"
584 WRITE (log_unit, '(T10,A,T19,A,T37,A,T49,A,T61,A,T67,A)') "number", "energy (eV)", &
585 "x", "y", "z", "strength (a.u.)"
586 WRITE (log_unit, '(T10,72("-"))')
587 END IF
588
589 ! transition dipole moment
590 ALLOCATE (trans_dipoles(nstates, nderivs, nspins))
591 trans_dipoles(:, :, :) = 0.0_dp
592
593 ! nspins == 1 .AND. mult == 3 : spin-flip transitions are forbidden due to symmetry reasons
594 IF (nspins > 1 .OR. mult == 1) THEN
595 DO ispin = 1, nspins
596 CALL cp_fm_get_info(dipole_op_mos_occ(1, ispin), nrow_global=nao, ncol_global=nocc)
597 CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive)
598 IF (nocc == nactive) THEN
599 DO ideriv = 1, nderivs
600 CALL cp_fm_trace(evects(ispin, :), dipole_op_mos_occ(ideriv, ispin), &
601 trans_dipoles(:, ideriv, ispin))
602 END DO
603 ELSE ! res
604 CALL cp_fm_get_info(evects(ispin, 1), matrix_struct=matrix_struct)
605 CALL cp_fm_create(dipact, matrix_struct)
606 DO ideriv = 1, nderivs
607 DO i = 1, nactive
608 j = gs_mos(ispin)%index_active(i)
609 CALL cp_fm_to_fm(dipole_op_mos_occ(ideriv, ispin), dipact, &
610 ncol=1, source_start=j, target_start=i)
611 END DO
612 CALL cp_fm_trace(evects(ispin, :), dipact, trans_dipoles(:, ideriv, ispin))
613 END DO
614 CALL cp_fm_release(dipact)
615 END IF
616 END DO
617
618 IF (nspins == 1) THEN
619 trans_dipoles(:, :, 1) = sqrt(2.0_dp)*trans_dipoles(:, :, 1)
620 ELSE
621 trans_dipoles(:, :, 1) = trans_dipoles(:, :, 1) + trans_dipoles(:, :, 2)
622 END IF
623 END IF
624
625 ! *** summary information ***
626 DO istate = 1, nstates
627
628 SELECT CASE (dipole_form)
630 osc_strength = 2.0_dp/3.0_dp*evals(istate)*sum(trans_dipoles(istate, :, 1)**2)
632 osc_strength = 2.0_dp/3.0_dp*evals(istate)*sum(trans_dipoles(istate, :, 1)**2)
634 osc_strength = 2.0_dp/3.0_dp/evals(istate)*sum(trans_dipoles(istate, :, 1)**2)
636 osc_strength = 2.0_dp/3.0_dp*evals(istate)*sum(trans_dipoles(istate, :, 1)**2)
637 CASE DEFAULT
638 cpabort("Unimplemented form of the dipole operator")
639 END SELECT
640
641 ostrength(istate) = osc_strength
642 IF (log_unit > 0) THEN
643 WRITE (log_unit, '(1X,A,T9,I7,T19,F11.5,T31,3(1X,ES11.4E2),T69,ES12.5E2)') &
644 "TDDFPT|", istate, evals(istate)*evolt, trans_dipoles(istate, 1:nderivs, 1), osc_strength
645 END IF
646 END DO
647
648 ! punch a checksum for the regs
649 IF (log_unit > 0) THEN
650 WRITE (log_unit, '(/,T2,A,E16.8)') 'TDDFPT : CheckSum E = ', sqrt(sum(evals**2))
651 WRITE (log_unit, '(/,T2,A,E16.8)') 'TDDFPT : CheckSum F = ', sqrt(sum(ostrength**2))
652 END IF
653
654 DEALLOCATE (trans_dipoles)
655
656 CALL timestop(handle)
657 END SUBROUTINE tddfpt_print_summary
658
659! **************************************************************************************************
660!> \brief Print excitation analysis.
661!> \param log_unit output unit
662!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
663!> SIZE(evects,2) -- number of excited states to print)
664!> \param evals TDDFPT eigenvalues
665!> \param gs_mos molecular orbitals optimised for the ground state
666!> \param matrix_s overlap matrix
667!> \param spinflip ...
668!> \param min_amplitude the smallest excitation amplitude to print
669!> \par History
670!> * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov]
671!> * 08.2018 splited of from 'tddfpt_print_summary' [Sergey Chulkov]
672! **************************************************************************************************
673 SUBROUTINE tddfpt_print_excitation_analysis(log_unit, evects, evals, gs_mos, matrix_s, spinflip, &
674 min_amplitude)
675 INTEGER, INTENT(in) :: log_unit
676 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
677 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
678 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
679 INTENT(in) :: gs_mos
680 TYPE(dbcsr_type), POINTER :: matrix_s
681 INTEGER :: spinflip
682 REAL(kind=dp), INTENT(in) :: min_amplitude
683
684 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_excitation_analysis'
685
686 CHARACTER(len=5) :: spin_label, spin_label2
687 INTEGER :: handle, icol, iproc, irow, ispin, &
688 istate, nao, ncols_local, nrows_local, &
689 nspins, nstates, spin2, state_spin, &
690 state_spin2
691 INTEGER(kind=int_8) :: iexc, imo_act, imo_occ, imo_virt, ind, &
692 nexcs, nexcs_local, nexcs_max_local, &
693 nmo_virt_occ, nmo_virt_occ_alpha
694 INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:) :: inds_local, inds_recv, nexcs_recv
695 INTEGER(kind=int_8), DIMENSION(1) :: nexcs_send
696 INTEGER(kind=int_8), DIMENSION(maxspins) :: nactive8, nmo_occ8, nmo_virt8
697 INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
698 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
699 INTEGER, DIMENSION(maxspins) :: nactive, nmo_occ, nmo_virt
700 LOGICAL :: do_exc_analysis
701 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights_local, weights_neg_abs_recv, &
702 weights_recv
703 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
704 POINTER :: local_data
705 TYPE(cp_blacs_env_type), POINTER :: blacs_env
706 TYPE(cp_fm_struct_type), POINTER :: fm_struct
707 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: s_mos_virt, weights_fm
708 TYPE(mp_para_env_type), POINTER :: para_env
709 TYPE(mp_request_type) :: send_handler, send_handler2
710 TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_handlers, recv_handlers2
711
712 CALL timeset(routinen, handle)
713
714 nspins = SIZE(gs_mos, 1)
715 nstates = SIZE(evects, 2)
716 do_exc_analysis = min_amplitude < 1.0_dp
717
718 CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env, para_env=para_env)
719 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
720
721 DO ispin = 1, nspins
722 nactive(ispin) = gs_mos(ispin)%nmo_active
723 nactive8(ispin) = int(nactive(ispin), kind=int_8)
724 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
725 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
726 nmo_occ8(ispin) = SIZE(gs_mos(ispin)%evals_occ, kind=int_8)
727 nmo_virt8(ispin) = SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
728 END DO
729
730 ! *** excitation analysis ***
731 IF (do_exc_analysis) THEN
732 cpassert(log_unit <= 0 .OR. para_env%is_source())
733 nmo_virt_occ_alpha = int(nmo_virt(1), int_8)*int(nmo_occ(1), int_8)
734
735 IF (log_unit > 0) THEN
736 WRITE (log_unit, "(1X,A)") "", &
737 "-------------------------------------------------------------------------------", &
738 "- Excitation analysis -", &
739 "-------------------------------------------------------------------------------"
740 WRITE (log_unit, '(8X,A,T27,A,T49,A,T69,A)') "State", "Occupied", "Virtual", "Excitation"
741 WRITE (log_unit, '(8X,A,T28,A,T49,A,T69,A)') "number", "orbital", "orbital", "amplitude"
742 WRITE (log_unit, '(1X,79("-"))')
743
744 IF (nspins == 1) THEN
745 state_spin = 1
746 state_spin2 = 2
747 spin_label = ' '
748 spin_label2 = ' '
749 ELSE IF (spinflip /= no_sf_tddfpt) THEN
750 state_spin = 1
751 state_spin2 = 2
752 spin_label = '(alp)'
753 spin_label2 = '(bet)'
754 END IF
755 END IF
756
757 ALLOCATE (s_mos_virt(SIZE(evects, 1)), weights_fm(SIZE(evects, 1)))
758 DO ispin = 1, SIZE(evects, 1)
759 IF (spinflip == no_sf_tddfpt) THEN
760 spin2 = ispin
761 ELSE
762 spin2 = 2
763 END IF
764 CALL cp_fm_get_info(gs_mos(spin2)%mos_virt, matrix_struct=fm_struct)
765 CALL cp_fm_create(s_mos_virt(ispin), fm_struct)
766 CALL cp_dbcsr_sm_fm_multiply(matrix_s, &
767 gs_mos(spin2)%mos_virt, &
768 s_mos_virt(ispin), &
769 ncol=nmo_virt(spin2), alpha=1.0_dp, beta=0.0_dp)
770
771 NULLIFY (fm_struct)
772 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(spin2), ncol_global=nactive(ispin), &
773 context=blacs_env)
774 CALL cp_fm_create(weights_fm(ispin), fm_struct)
775 CALL cp_fm_set_all(weights_fm(ispin), 0.0_dp)
776 CALL cp_fm_struct_release(fm_struct)
777 END DO
778
779 nexcs_max_local = 0
780 DO ispin = 1, SIZE(evects, 1)
781 CALL cp_fm_get_info(weights_fm(ispin), nrow_local=nrows_local, ncol_local=ncols_local)
782 nexcs_max_local = nexcs_max_local + int(nrows_local, int_8)*int(ncols_local, int_8)
783 END DO
784
785 ALLOCATE (weights_local(nexcs_max_local), inds_local(nexcs_max_local))
786
787 DO istate = 1, nstates
788 nexcs_local = 0
789 nmo_virt_occ = 0
790
791 ! analyse matrix elements locally and transfer only significant
792 ! excitations to the master node for subsequent ordering
793 DO ispin = 1, SIZE(evects, 1)
794 IF (spinflip == no_sf_tddfpt) THEN
795 spin2 = ispin
796 ELSE
797 spin2 = 2
798 END IF
799 ! compute excitation amplitudes
800 CALL parallel_gemm('T', 'N', nmo_virt(spin2), nactive(ispin), nao, 1.0_dp, s_mos_virt(ispin), &
801 evects(ispin, istate), 0.0_dp, weights_fm(ispin))
802
803 CALL cp_fm_get_info(weights_fm(ispin), nrow_local=nrows_local, ncol_local=ncols_local, &
804 row_indices=row_indices, col_indices=col_indices, local_data=local_data)
805
806 ! locate single excitations with significant amplitudes (>= min_amplitude)
807 DO icol = 1, ncols_local
808 DO irow = 1, nrows_local
809 IF (abs(local_data(irow, icol)) >= min_amplitude) THEN
810 ! number of non-negligible excitations
811 nexcs_local = nexcs_local + 1
812 ! excitation amplitude
813 weights_local(nexcs_local) = local_data(irow, icol)
814 ! index of single excitation (ivirt, iocc, ispin) in compressed form
815 inds_local(nexcs_local) = nmo_virt_occ + int(row_indices(irow), int_8) + &
816 int(col_indices(icol) - 1, int_8)*nmo_virt8(spin2)
817 END IF
818 END DO
819 END DO
820
821 nmo_virt_occ = nmo_virt_occ + nmo_virt8(spin2)*nmo_occ8(ispin)
822 END DO
823
824 IF (para_env%is_source()) THEN
825 ! master node
826 ALLOCATE (nexcs_recv(para_env%num_pe), recv_handlers(para_env%num_pe), recv_handlers2(para_env%num_pe))
827
828 ! collect number of non-negligible excitations from other nodes
829 DO iproc = 1, para_env%num_pe
830 IF (iproc - 1 /= para_env%mepos) THEN
831 CALL para_env%irecv(nexcs_recv(iproc:iproc), iproc - 1, recv_handlers(iproc), 0)
832 ELSE
833 nexcs_recv(iproc) = nexcs_local
834 END IF
835 END DO
836
837 DO iproc = 1, para_env%num_pe
838 IF (iproc - 1 /= para_env%mepos) &
839 CALL recv_handlers(iproc)%wait()
840 END DO
841
842 ! compute total number of non-negligible excitations
843 nexcs = 0
844 DO iproc = 1, para_env%num_pe
845 nexcs = nexcs + nexcs_recv(iproc)
846 END DO
847
848 ! receive indices and amplitudes of selected excitations
849 ALLOCATE (weights_recv(nexcs), weights_neg_abs_recv(nexcs))
850 ALLOCATE (inds_recv(nexcs), inds(nexcs))
851
852 nmo_virt_occ = 0
853 DO iproc = 1, para_env%num_pe
854 IF (nexcs_recv(iproc) > 0) THEN
855 IF (iproc - 1 /= para_env%mepos) THEN
856 ! excitation amplitudes
857 CALL para_env%irecv(weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), &
858 iproc - 1, recv_handlers(iproc), 1)
859 ! compressed indices
860 CALL para_env%irecv(inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), &
861 iproc - 1, recv_handlers2(iproc), 2)
862 ELSE
863 ! data on master node
864 weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = weights_local(1:nexcs_recv(iproc))
865 inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = inds_local(1:nexcs_recv(iproc))
866 END IF
867
868 nmo_virt_occ = nmo_virt_occ + nexcs_recv(iproc)
869 END IF
870 END DO
871
872 DO iproc = 1, para_env%num_pe
873 IF (iproc - 1 /= para_env%mepos .AND. nexcs_recv(iproc) > 0) THEN
874 CALL recv_handlers(iproc)%wait()
875 CALL recv_handlers2(iproc)%wait()
876 END IF
877 END DO
878
879 DEALLOCATE (nexcs_recv, recv_handlers, recv_handlers2)
880 ELSE
881 ! working node: send the number of selected excited states to the master node
882 nexcs_send(1) = nexcs_local
883 CALL para_env%isend(nexcs_send, para_env%source, send_handler, 0)
884 CALL send_handler%wait()
885
886 IF (nexcs_local > 0) THEN
887 ! send excitation amplitudes
888 CALL para_env%isend(weights_local(1:nexcs_local), para_env%source, send_handler, 1)
889 ! send compressed indices
890 CALL para_env%isend(inds_local(1:nexcs_local), para_env%source, send_handler2, 2)
891
892 CALL send_handler%wait()
893 CALL send_handler2%wait()
894 END IF
895 END IF
896
897 ! sort non-negligible excitations on the master node according to their amplitudes,
898 ! uncompress indices and print summary information
899 IF (para_env%is_source() .AND. log_unit > 0) THEN
900 weights_neg_abs_recv(:) = -abs(weights_recv)
901 CALL sort(weights_neg_abs_recv, int(nexcs), inds)
902
903 WRITE (log_unit, '(T7,I8,F10.5,A)') istate, evals(istate)*evolt, " eV"
904
905 ! This reinitialization is needed to prevent the intel fortran compiler from introduce
906 ! a bug when using optimization level 3 flag
907 state_spin = 1
908 state_spin2 = 1
909 IF (spinflip /= no_sf_tddfpt) THEN
910 state_spin = 1
911 state_spin2 = 2
912 END IF
913 DO iexc = 1, nexcs
914 ind = inds_recv(inds(iexc)) - 1
915 IF ((nspins > 1) .AND. (spinflip == no_sf_tddfpt)) THEN
916 IF (ind < nmo_virt_occ_alpha) THEN
917 state_spin = 1
918 state_spin2 = 1
919 spin_label = '(alp)'
920 spin_label2 = '(alp)'
921 ELSE
922 state_spin = 2
923 state_spin2 = 2
924 ind = ind - nmo_virt_occ_alpha
925 spin_label = '(bet)'
926 spin_label2 = '(bet)'
927 END IF
928 END IF
929 imo_act = ind/nmo_virt8(state_spin2) + 1
930 imo_occ = gs_mos(state_spin)%index_active(imo_act)
931 imo_virt = mod(ind, nmo_virt8(state_spin2)) + 1
932
933 WRITE (log_unit, '(T27,I8,1X,A5,T48,I8,1X,A5,T70,F9.6)') imo_occ, spin_label, &
934 nmo_occ8(state_spin2) + imo_virt, spin_label2, weights_recv(inds(iexc))
935 END DO
936 END IF
937
938 ! deallocate temporary arrays
939 IF (para_env%is_source()) &
940 DEALLOCATE (weights_recv, weights_neg_abs_recv, inds_recv, inds)
941 END DO
942
943 DEALLOCATE (weights_local, inds_local)
944 IF (log_unit > 0) THEN
945 WRITE (log_unit, "(1X,A)") &
946 "-------------------------------------------------------------------------------"
947 END IF
948 END IF
949
950 CALL cp_fm_release(weights_fm)
951 CALL cp_fm_release(s_mos_virt)
952
953 CALL timestop(handle)
954
956
957! **************************************************************************************************
958!> \brief Print natural transition orbital analysis.
959!> \param qs_env Information on Kinds and Particles
960!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
961!> SIZE(evects,2) -- number of excited states to print)
962!> \param evals TDDFPT eigenvalues
963!> \param ostrength ...
964!> \param gs_mos molecular orbitals optimised for the ground state
965!> \param matrix_s overlap matrix
966!> \param print_section ...
967!> \par History
968!> * 06.2019 created [JGH]
969! **************************************************************************************************
970 SUBROUTINE tddfpt_print_nto_analysis(qs_env, evects, evals, ostrength, gs_mos, matrix_s, print_section)
971 TYPE(qs_environment_type), POINTER :: qs_env
972 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
973 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals, ostrength
974 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
975 INTENT(in) :: gs_mos
976 TYPE(dbcsr_type), POINTER :: matrix_s
977 TYPE(section_vals_type), POINTER :: print_section
978
979 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_nto_analysis'
980 INTEGER, PARAMETER :: ntomax = 10
981
982 CHARACTER(LEN=20), DIMENSION(2) :: nto_name
983 INTEGER :: handle, i, ia, icg, iounit, ispin, &
984 istate, j, nao, nlist, nmax, nmo, &
985 nnto, nspins, nstates
986 INTEGER, DIMENSION(2) :: iv
987 INTEGER, DIMENSION(2, ntomax) :: ia_index
988 INTEGER, DIMENSION(:), POINTER :: slist, stride
989 LOGICAL :: append_cube, cube_file, explicit
990 REAL(kind=dp) :: os_threshold, sume, threshold
991 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvals
992 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvalues
993 REAL(kind=dp), DIMENSION(ntomax) :: ia_eval
994 TYPE(cell_type), POINTER :: cell
995 TYPE(cp_fm_struct_type), POINTER :: fm_mo_struct, fm_struct
996 TYPE(cp_fm_type) :: sev, smat, tmat, wmat, work, wvec
997 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: teig
998 TYPE(cp_logger_type), POINTER :: logger
999 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: nto_set
1000 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1001 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1002 TYPE(section_vals_type), POINTER :: molden_section, nto_section
1003
1004 CALL timeset(routinen, handle)
1005
1006 logger => cp_get_default_logger()
1007 iounit = cp_logger_get_default_io_unit(logger)
1008
1009 IF (btest(cp_print_key_should_output(logger%iter_info, print_section, &
1010 "NTO_ANALYSIS"), cp_p_file)) THEN
1011
1012 CALL cite_reference(martin2003)
1013
1014 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%THRESHOLD", r_val=threshold)
1015 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%INTENSITY_THRESHOLD", r_val=os_threshold)
1016 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", explicit=explicit)
1017
1018 IF (explicit) THEN
1019 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", i_vals=slist)
1020 nlist = SIZE(slist)
1021 ELSE
1022 nlist = 0
1023 END IF
1024
1025 IF (iounit > 0) THEN
1026 WRITE (iounit, "(1X,A)") "", &
1027 "-------------------------------------------------------------------------------", &
1028 "- Natural Orbital analysis -", &
1029 "-------------------------------------------------------------------------------"
1030 END IF
1031
1032 nspins = SIZE(evects, 1)
1033 nstates = SIZE(evects, 2)
1034 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1035
1036 DO istate = 1, nstates
1037 IF (os_threshold > ostrength(istate)) THEN
1038 IF (iounit > 0) THEN
1039 WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1040 END IF
1041 cycle
1042 END IF
1043 IF (nlist > 0) THEN
1044 IF (.NOT. any(slist == istate)) THEN
1045 IF (iounit > 0) THEN
1046 WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1047 END IF
1048 cycle
1049 END IF
1050 END IF
1051 IF (iounit > 0) THEN
1052 WRITE (iounit, "(1X,A,I6,T30,F10.5,A)") " STATE NR. ", istate, evals(istate)*evolt, " eV"
1053 END IF
1054 nmax = 0
1055 DO ispin = 1, nspins
1056 CALL cp_fm_get_info(evects(ispin, istate), matrix_struct=fm_struct, ncol_global=nmo)
1057 nmax = max(nmax, nmo)
1058 END DO
1059 ALLOCATE (eigenvalues(nmax, nspins))
1060 eigenvalues = 0.0_dp
1061 ! SET 1: Hole states
1062 ! SET 2: Particle states
1063 nto_name(1) = 'Hole_states'
1064 nto_name(2) = 'Particle_states'
1065 ALLOCATE (nto_set(2))
1066 DO i = 1, 2
1067 CALL allocate_mo_set(nto_set(i), nao, ntomax, 0, 0.0_dp, 1.0_dp, 0.0_dp)
1068 CALL cp_fm_get_info(evects(1, istate), matrix_struct=fm_struct)
1069 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1070 ncol_global=ntomax)
1071 CALL cp_fm_create(tmat, fm_mo_struct)
1072 CALL init_mo_set(nto_set(i), fm_ref=tmat, name=nto_name(i))
1073 CALL cp_fm_release(tmat)
1074 CALL cp_fm_struct_release(fm_mo_struct)
1075 END DO
1076 !
1077 ALLOCATE (teig(nspins))
1078 ! hole states
1079 ! Diagonalize X(T)*S*X
1080 DO ispin = 1, nspins
1081 associate(ev => evects(ispin, istate))
1082 CALL cp_fm_get_info(ev, matrix_struct=fm_struct, ncol_global=nmo)
1083 CALL cp_fm_create(sev, fm_struct)
1084 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1085 nrow_global=nmo, ncol_global=nmo)
1086 CALL cp_fm_create(tmat, fm_mo_struct)
1087 CALL cp_fm_create(teig(ispin), fm_mo_struct)
1088 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1089 CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, ev, sev, 0.0_dp, tmat)
1090 END associate
1091
1092 CALL choose_eigv_solver(tmat, teig(ispin), eigenvalues(1:nmo, ispin))
1093
1094 CALL cp_fm_struct_release(fm_mo_struct)
1095 CALL cp_fm_release(tmat)
1096 CALL cp_fm_release(sev)
1097 END DO
1098 ! find major determinants i->a
1099 ia_index = 0
1100 sume = 0.0_dp
1101 nnto = 0
1102 DO i = 1, ntomax
1103 iv = maxloc(eigenvalues)
1104 ia_eval(i) = eigenvalues(iv(1), iv(2))
1105 ia_index(1:2, i) = iv(1:2)
1106 sume = sume + ia_eval(i)
1107 eigenvalues(iv(1), iv(2)) = 0.0_dp
1108 nnto = nnto + 1
1109 IF (sume > threshold) EXIT
1110 END DO
1111 ! store hole states
1112 CALL set_mo_set(nto_set(1), nmo=nnto)
1113 DO i = 1, nnto
1114 ia = ia_index(1, i)
1115 ispin = ia_index(2, i)
1116 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=nmo)
1117 CALL cp_fm_get_info(teig(ispin), matrix_struct=fm_struct)
1118 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1119 nrow_global=nmo, ncol_global=1)
1120 CALL cp_fm_create(tmat, fm_mo_struct)
1121 CALL cp_fm_struct_release(fm_mo_struct)
1122 CALL cp_fm_get_info(gs_mos(1)%mos_occ, matrix_struct=fm_struct)
1123 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1124 ncol_global=1)
1125 CALL cp_fm_create(wvec, fm_mo_struct)
1126 CALL cp_fm_struct_release(fm_mo_struct)
1127 CALL cp_fm_to_fm(teig(ispin), tmat, 1, ia, 1)
1128 CALL parallel_gemm('N', 'N', nao, 1, nmo, 1.0_dp, gs_mos(ispin)%mos_occ, &
1129 tmat, 0.0_dp, wvec)
1130 CALL cp_fm_to_fm(wvec, nto_set(1)%mo_coeff, 1, 1, i)
1131 CALL cp_fm_release(wvec)
1132 CALL cp_fm_release(tmat)
1133 END DO
1134 ! particle states
1135 ! Solve generalized eigenvalue equation: (S*X)*(S*X)(T)*v = lambda*S*v
1136 CALL set_mo_set(nto_set(2), nmo=nnto)
1137 DO ispin = 1, nspins
1138 associate(ev => evects(ispin, istate))
1139 CALL cp_fm_get_info(ev, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
1140 ALLOCATE (eigvals(nao))
1141 eigvals = 0.0_dp
1142 CALL cp_fm_create(sev, fm_struct)
1143 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1144 END associate
1145 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1146 nrow_global=nao, ncol_global=nao)
1147 CALL cp_fm_create(tmat, fm_mo_struct)
1148 CALL cp_fm_create(smat, fm_mo_struct)
1149 CALL cp_fm_create(wmat, fm_mo_struct)
1150 CALL cp_fm_create(work, fm_mo_struct)
1151 CALL cp_fm_struct_release(fm_mo_struct)
1152 CALL copy_dbcsr_to_fm(matrix_s, smat)
1153 CALL parallel_gemm('N', 'T', nao, nao, nmo, 1.0_dp, sev, sev, 0.0_dp, tmat)
1154 CALL cp_fm_geeig(tmat, smat, wmat, eigvals, work)
1155 DO i = 1, nnto
1156 IF (ispin == ia_index(2, i)) THEN
1157 icg = 0
1158 DO j = 1, nao
1159 IF (abs(eigvals(j) - ia_eval(i)) < 1.e-6_dp) THEN
1160 icg = j
1161 EXIT
1162 END IF
1163 END DO
1164 IF (icg == 0) THEN
1165 CALL cp_warn(__location__, &
1166 "Could not locate particle state associated with hole state.")
1167 ELSE
1168 CALL cp_fm_to_fm(wmat, nto_set(2)%mo_coeff, 1, icg, i)
1169 END IF
1170 END IF
1171 END DO
1172 DEALLOCATE (eigvals)
1173 CALL cp_fm_release(sev)
1174 CALL cp_fm_release(tmat)
1175 CALL cp_fm_release(smat)
1176 CALL cp_fm_release(wmat)
1177 CALL cp_fm_release(work)
1178 END DO
1179 ! print
1180 IF (iounit > 0) THEN
1181 sume = 0.0_dp
1182 DO i = 1, nnto
1183 sume = sume + ia_eval(i)
1184 WRITE (iounit, "(T6,A,i2,T30,A,i1,T42,A,F8.5,T63,A,F8.5)") &
1185 "Particle-Hole state:", i, " Spin:", ia_index(2, i), &
1186 "Eigenvalue:", ia_eval(i), " Sum Eigv:", sume
1187 END DO
1188 END IF
1189 ! Cube and Molden files
1190 nto_section => section_vals_get_subs_vals(print_section, "NTO_ANALYSIS")
1191 CALL section_vals_val_get(nto_section, "CUBE_FILES", l_val=cube_file)
1192 CALL section_vals_val_get(nto_section, "STRIDE", i_vals=stride)
1193 CALL section_vals_val_get(nto_section, "APPEND", l_val=append_cube)
1194 IF (cube_file) THEN
1195 CALL print_nto_cubes(qs_env, nto_set, istate, stride, append_cube, nto_section)
1196 END IF
1197 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, cell=cell)
1198 molden_section => section_vals_get_subs_vals(print_section, "MOS_MOLDEN")
1199 CALL write_mos_molden(nto_set, qs_kind_set, particle_set, molden_section, cell=cell, qs_env=qs_env)
1200 !
1201 DEALLOCATE (eigenvalues)
1202 CALL cp_fm_release(teig)
1203 !
1204 DO i = 1, 2
1205 CALL deallocate_mo_set(nto_set(i))
1206 END DO
1207 DEALLOCATE (nto_set)
1208 END DO
1209
1210 IF (iounit > 0) THEN
1211 WRITE (iounit, "(1X,A)") &
1212 "-------------------------------------------------------------------------------"
1213 END IF
1214
1215 END IF
1216
1217 CALL timestop(handle)
1218
1219 END SUBROUTINE tddfpt_print_nto_analysis
1220
1221! **************************************************************************************************
1222!> \brief Print exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018)
1223!> \param log_unit output unit
1224!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
1225!> SIZE(evects,2) -- number of excited states to print)
1226!> \param gs_mos molecular orbitals optimised for the ground state
1227!> \param matrix_s overlap matrix
1228!> \param do_directional_exciton_descriptors flag for computing descriptors for each (cartesian) direction
1229!> \param qs_env Information on particles/geometry
1230!> \par History
1231!> * 12.2024 created as 'tddfpt_print_exciton_descriptors' [Maximilian Graml]
1232! **************************************************************************************************
1233 SUBROUTINE tddfpt_print_exciton_descriptors(log_unit, evects, gs_mos, matrix_s, &
1234 do_directional_exciton_descriptors, qs_env)
1235 INTEGER, INTENT(in) :: log_unit
1236 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
1237 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1238 INTENT(in) :: gs_mos
1239 TYPE(dbcsr_type), POINTER :: matrix_s
1240 LOGICAL, INTENT(IN) :: do_directional_exciton_descriptors
1241 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1242
1243 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_exciton_descriptors'
1244
1245 CHARACTER(LEN=4) :: prefix_output
1246 INTEGER :: handle, ispin, istate, n_moments_quad, &
1247 nactive, nao, nspins, nstates
1248 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt
1249 LOGICAL :: print_checkvalue
1250 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ref_point_multipole
1251 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1252 TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo_coeff, &
1253 fm_struct_s_mos_virt, fm_struct_x_ia_n
1254 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: eigvec_x_ia_n, fm_multipole_ab, &
1255 fm_multipole_ai, fm_multipole_ij, &
1256 s_mos_virt
1257 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_coeff
1258 TYPE(exciton_descr_type), ALLOCATABLE, &
1259 DIMENSION(:) :: exc_descr
1260
1261 CALL timeset(routinen, handle)
1262
1263 nspins = SIZE(evects, 1)
1264 nstates = SIZE(evects, 2)
1265
1266 cpassert(nspins == 1) ! Other spins are not yet implemented for exciton descriptors
1267
1268 CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env)
1269 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1270
1271 DO ispin = 1, nspins
1272 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
1273 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
1274 END DO
1275
1276 ! Prepare fm with all MO coefficents, i.e. nao x nao
1277 ALLOCATE (mo_coeff(nspins))
1278 CALL cp_fm_struct_create(fm_struct_mo_coeff, nrow_global=nao, ncol_global=nao, &
1279 context=blacs_env)
1280 DO ispin = 1, nspins
1281 CALL cp_fm_create(mo_coeff(ispin), fm_struct_mo_coeff)
1282 CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_occ, &
1283 mo_coeff(ispin), &
1284 nao, &
1285 nmo_occ(ispin), &
1286 1, &
1287 1, &
1288 1, &
1289 1, &
1290 blacs_env)
1291 CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_virt, &
1292 mo_coeff(ispin), &
1293 nao, &
1294 nmo_virt(ispin), &
1295 1, &
1296 1, &
1297 1, &
1298 nmo_occ(ispin) + 1, &
1299 blacs_env)
1300 END DO
1301 CALL cp_fm_struct_release(fm_struct_mo_coeff)
1302
1303 ! Compute multipole moments
1304 ! fm_multipole_XY have structure inherited by libint, i.e. x, y, z, xx, xy, xz, yy, yz, zz
1305 n_moments_quad = 9
1306 ALLOCATE (ref_point_multipole(3))
1307 ALLOCATE (fm_multipole_ij(n_moments_quad))
1308 ALLOCATE (fm_multipole_ab(n_moments_quad))
1309 ALLOCATE (fm_multipole_ai(n_moments_quad))
1310
1311 CALL get_multipoles_mo(fm_multipole_ai, fm_multipole_ij, fm_multipole_ab, &
1312 qs_env, mo_coeff, ref_point_multipole, 2, &
1313 nmo_occ(1), nmo_virt(1), blacs_env)
1314
1315 CALL cp_fm_release(mo_coeff)
1316
1317 ! Compute eigenvector X of the Casida equation from trial vectors
1318 ALLOCATE (s_mos_virt(nspins), eigvec_x_ia_n(nspins))
1319 DO ispin = 1, nspins
1320 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct_s_mos_virt)
1321 CALL cp_fm_create(s_mos_virt(ispin), fm_struct_s_mos_virt)
1322 NULLIFY (fm_struct_s_mos_virt)
1323 CALL cp_dbcsr_sm_fm_multiply(matrix_s, &
1324 gs_mos(ispin)%mos_virt, &
1325 s_mos_virt(ispin), &
1326 ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
1327
1328 CALL cp_fm_struct_create(fm_struct_x_ia_n, nrow_global=nmo_occ(ispin), ncol_global=nmo_virt(ispin), &
1329 context=blacs_env)
1330 CALL cp_fm_create(eigvec_x_ia_n(ispin), fm_struct_x_ia_n)
1331 CALL cp_fm_struct_release(fm_struct_x_ia_n)
1332 END DO
1333 ALLOCATE (exc_descr(nstates))
1334 DO istate = 1, nstates
1335 DO ispin = 1, nspins
1336 CALL cp_fm_set_all(eigvec_x_ia_n(ispin), 0.0_dp)
1337 ! compute eigenvectors X of the TDA equation
1338 ! Reshuffle multiplication from
1339 ! X_ai = S_ma ^T * C_mi
1340 ! to
1341 ! X_ia = C_mi ^T * S_ma
1342 ! for compatibility with the structure needed for get_exciton_descriptors of bse_properties.F
1343 CALL cp_fm_get_info(evects(ispin, istate), ncol_global=nactive)
1344 IF (nactive /= nmo_occ(ispin)) THEN
1345 CALL cp_abort(__location__, &
1346 "Reduced active space excitations not implemented")
1347 END IF
1348 CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_virt(ispin), nao, 1.0_dp, &
1349 evects(ispin, istate), s_mos_virt(ispin), 0.0_dp, eigvec_x_ia_n(ispin))
1350
1351 CALL get_exciton_descriptors(exc_descr, eigvec_x_ia_n(ispin), &
1352 fm_multipole_ij, fm_multipole_ab, &
1353 fm_multipole_ai, &
1354 istate, nmo_occ(ispin), nmo_virt(ispin))
1355 END DO
1356 END DO
1357 CALL cp_fm_release(eigvec_x_ia_n)
1358 CALL cp_fm_release(s_mos_virt)
1359 CALL cp_fm_release(fm_multipole_ai)
1360 CALL cp_fm_release(fm_multipole_ij)
1361 CALL cp_fm_release(fm_multipole_ab)
1362
1363 ! Actual printing
1364 print_checkvalue = .true.
1365 prefix_output = ' '
1366 CALL print_exciton_descriptors(exc_descr, ref_point_multipole, log_unit, &
1367 nstates, print_checkvalue, do_directional_exciton_descriptors, &
1368 prefix_output, qs_env)
1369
1370 DEALLOCATE (ref_point_multipole)
1371 DEALLOCATE (exc_descr)
1372
1373 CALL timestop(handle)
1374
1376
1377! **************************************************************************************************
1378!> \brief ...
1379!> \param vin ...
1380!> \param vout ...
1381!> \param mos_occ ...
1382!> \param matrix_s ...
1383! **************************************************************************************************
1384 SUBROUTINE project_vector(vin, vout, mos_occ, matrix_s)
1385 TYPE(dbcsr_type) :: vin, vout
1386 TYPE(cp_fm_type), INTENT(IN) :: mos_occ
1387 TYPE(dbcsr_type), POINTER :: matrix_s
1388
1389 CHARACTER(LEN=*), PARAMETER :: routinen = 'project_vector'
1390
1391 INTEGER :: handle, nao, nmo
1392 REAL(kind=dp) :: norm(1)
1393 TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_vec_struct
1394 TYPE(cp_fm_type) :: csvec, svec, vec
1395
1396 CALL timeset(routinen, handle)
1397
1398 CALL cp_fm_get_info(mos_occ, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
1399 CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
1400 nrow_global=nao, ncol_global=1)
1401 CALL cp_fm_create(vec, fm_vec_struct)
1402 CALL cp_fm_create(svec, fm_vec_struct)
1403 CALL cp_fm_struct_release(fm_vec_struct)
1404 CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
1405 nrow_global=nmo, ncol_global=1)
1406 CALL cp_fm_create(csvec, fm_vec_struct)
1407 CALL cp_fm_struct_release(fm_vec_struct)
1408
1409 CALL copy_dbcsr_to_fm(vin, vec)
1410 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vec, svec, ncol=1, alpha=1.0_dp, beta=0.0_dp)
1411 CALL parallel_gemm('T', 'N', nmo, 1, nao, 1.0_dp, mos_occ, svec, 0.0_dp, csvec)
1412 CALL parallel_gemm('N', 'N', nao, 1, nmo, -1.0_dp, mos_occ, csvec, 1.0_dp, vec)
1413 CALL cp_fm_vectorsnorm(vec, norm)
1414 cpassert(norm(1) > 1.e-14_dp)
1415 norm(1) = sqrt(1._dp/norm(1))
1416 CALL cp_fm_scale(norm(1), vec)
1417 CALL copy_fm_to_dbcsr(vec, vout, keep_sparsity=.false.)
1418
1419 CALL cp_fm_release(csvec)
1420 CALL cp_fm_release(svec)
1421 CALL cp_fm_release(vec)
1422
1423 CALL timestop(handle)
1424
1425 END SUBROUTINE project_vector
1426
1427! **************************************************************************************************
1428!> \brief ...
1429!> \param va ...
1430!> \param vb ...
1431!> \param res ...
1432! **************************************************************************************************
1433 SUBROUTINE vec_product(va, vb, res)
1434 TYPE(dbcsr_type) :: va, vb
1435 REAL(kind=dp), INTENT(OUT) :: res
1436
1437 CHARACTER(LEN=*), PARAMETER :: routinen = 'vec_product'
1438
1439 INTEGER :: handle, icol, irow
1440 LOGICAL :: found
1441 REAL(kind=dp), DIMENSION(:, :), POINTER :: vba, vbb
1442 TYPE(dbcsr_iterator_type) :: iter
1443 TYPE(mp_comm_type) :: group
1444
1445 CALL timeset(routinen, handle)
1446
1447 res = 0.0_dp
1448
1449 CALL dbcsr_get_info(va, group=group)
1450 CALL dbcsr_iterator_start(iter, va)
1451 DO WHILE (dbcsr_iterator_blocks_left(iter))
1452 CALL dbcsr_iterator_next_block(iter, irow, icol, vba)
1453 CALL dbcsr_get_block_p(vb, row=irow, col=icol, block=vbb, found=found)
1454 res = res + sum(vba*vbb)
1455 cpassert(found)
1456 END DO
1457 CALL dbcsr_iterator_stop(iter)
1458 CALL group%sum(res)
1459
1460 CALL timestop(handle)
1461
1462 END SUBROUTINE vec_product
1463
1464! **************************************************************************************************
1465!> \brief ...
1466!> \param qs_env ...
1467!> \param mos ...
1468!> \param istate ...
1469!> \param stride ...
1470!> \param append_cube ...
1471!> \param print_section ...
1472! **************************************************************************************************
1473 SUBROUTINE print_nto_cubes(qs_env, mos, istate, stride, append_cube, print_section)
1474
1475 TYPE(qs_environment_type), POINTER :: qs_env
1476 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1477 INTEGER, INTENT(IN) :: istate
1478 INTEGER, DIMENSION(:), POINTER :: stride
1479 LOGICAL, INTENT(IN) :: append_cube
1480 TYPE(section_vals_type), POINTER :: print_section
1481
1482 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
1483 INTEGER :: i, iset, nmo, unit_nr
1484 LOGICAL :: mpi_io
1485 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1486 TYPE(cell_type), POINTER :: cell
1487 TYPE(cp_fm_type), POINTER :: mo_coeff
1488 TYPE(cp_logger_type), POINTER :: logger
1489 TYPE(dft_control_type), POINTER :: dft_control
1490 TYPE(particle_list_type), POINTER :: particles
1491 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1492 TYPE(pw_c1d_gs_type) :: wf_g
1493 TYPE(pw_env_type), POINTER :: pw_env
1494 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1495 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1496 TYPE(pw_r3d_rs_type) :: wf_r
1497 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1498 TYPE(qs_subsys_type), POINTER :: subsys
1499
1500 logger => cp_get_default_logger()
1501
1502 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
1503 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1504 CALL auxbas_pw_pool%create_pw(wf_r)
1505 CALL auxbas_pw_pool%create_pw(wf_g)
1506
1507 CALL get_qs_env(qs_env, subsys=subsys)
1508 CALL qs_subsys_get(subsys, particles=particles)
1509
1510 my_pos_cube = "REWIND"
1511 IF (append_cube) THEN
1512 my_pos_cube = "APPEND"
1513 END IF
1514
1515 CALL get_qs_env(qs_env=qs_env, &
1516 atomic_kind_set=atomic_kind_set, &
1517 qs_kind_set=qs_kind_set, &
1518 cell=cell, &
1519 particle_set=particle_set)
1520
1521 DO iset = 1, 2
1522 CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo)
1523 DO i = 1, nmo
1524 CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1525 cell, dft_control, particle_set, pw_env)
1526 IF (iset == 1) THEN
1527 WRITE (filename, '(a4,I3.3,I2.2,a11)') "NTO_STATE", istate, i, "_Hole_State"
1528 ELSEIF (iset == 2) THEN
1529 WRITE (filename, '(a4,I3.3,I2.2,a15)') "NTO_STATE", istate, i, "_Particle_State"
1530 END IF
1531 mpi_io = .true.
1532 unit_nr = cp_print_key_unit_nr(logger, print_section, '', extension=".cube", &
1533 middle_name=trim(filename), file_position=my_pos_cube, &
1534 log_filename=.false., ignore_should_output=.true., mpi_io=mpi_io)
1535 IF (iset == 1) THEN
1536 WRITE (title, *) "Natural Transition Orbital Hole State", i
1537 ELSEIF (iset == 2) THEN
1538 WRITE (title, *) "Natural Transition Orbital Particle State", i
1539 END IF
1540 CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, stride=stride, mpi_io=mpi_io)
1541 CALL cp_print_key_finished_output(unit_nr, logger, print_section, '', &
1542 ignore_should_output=.true., mpi_io=mpi_io)
1543 END DO
1544 END DO
1545
1546 CALL auxbas_pw_pool%give_back_pw(wf_g)
1547 CALL auxbas_pw_pool%give_back_pw(wf_r)
1548
1549 END SUBROUTINE print_nto_cubes
1550
1551END MODULE qs_tddfpt2_properties
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 martin2003
Routines for printing information in context of the BSE calculation.
Definition bse_print.F:13
subroutine, public print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, num_print_exc_descr, print_checkvalue, print_directional_exc_descr, prefix_output, qs_env)
Prints exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018)
Definition bse_print.F:507
Routines for computing excitonic properties, e.g. exciton diameter, from the BSE.
subroutine, public get_exciton_descriptors(exc_descr, fm_x_ia, fm_multipole_ij_trunc, fm_multipole_ab_trunc, fm_multipole_ai_trunc, i_exc, homo, virtual, fm_y_ia)
...
Auxiliary routines for GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_util.F:13
subroutine, public get_multipoles_mo(fm_multipole_ai_trunc, fm_multipole_ij_trunc, fm_multipole_ab_trunc, qs_env, mo_coeff, rpoint, n_moments, homo_red, virtual_red, context_bse)
...
Definition bse_util.F:1715
Handles all functions related to the CELL.
Definition cell_types.F:15
methods related to the blacs parallel environment
Basic linear algebra operations for complex full matrices.
subroutine, public cp_cfm_solve(matrix_a, general_a, determinant)
Solve the system of linear equations A*b=A_general using LU decomposition. Pay attention that both ma...
Represents a complex full matrix distributed on many processors.
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_set(matrix, alpha)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
subroutine, public cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE. Use cuSOLVERMp directly when requested and large enough; otherwi...
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:245
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_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_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_vectorsnorm(matrix, norm_array)
find the inorm of each column norm_{j}= sqrt( \sum_{i} A_{ij}*A_{ij} )
subroutine, public cp_fm_to_fm_submat_general(source, destination, nrows, ncols, s_firstrow, s_firstcol, d_firstrow, d_firstcol, global_context)
General copy of a submatrix of fm matrix to a submatrix of another fm matrix. The two matrices can ha...
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
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,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io)
...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public tddfpt_dipole_berry
integer, parameter, public tddfpt_dipole_scf_moment
integer, parameter, public tddfpt_dipole_velocity
integer, parameter, public tddfpt_dipole_length
integer, parameter, public tddfpt_dipole_velocity_old
integer, parameter, public no_sf_tddfpt
objects that represent the structure of input sections and the data contained in an input section
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_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 int_8
Definition kinds.F:54
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_path_length
Definition kinds.F:58
Definition of mathematical constants and functions.
complex(kind=dp), parameter, public z_one
real(kind=dp), parameter, public twopi
complex(kind=dp), parameter, public z_zero
Interface to the message passing library MPI.
Functions handling the MOLDEN format. Split from mode_selective.
subroutine, public write_mos_molden(mos, qs_kind_set, particle_set, print_section, cell, unoccupied_orbs, unoccupied_evals, qs_env, calc_energies)
Write out the MOs in molden format for visualisation.
Calculates the moment integrals <a|r^m|b>
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
basic linear algebra operations for full matrixes
represent a simple array based list of the given type
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
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type)
maps a given wavefunction on the grid
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, mimic, 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, xcint_weights, 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.
Define the quickstep kind type and their sub types.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public set_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, uniform_occupation, kts, mu, flexible_electron_count)
Set the components of a MO set data structure.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name, counter)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition qs_moments.F:14
subroutine, public build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
...
Define the neighbor list data types and the corresponding functionality.
subroutine, public rrc_xyz_ao(op, qs_env, rc, order, minimum_image, soft)
Calculation of the components of the dipole operator in the length form by taking the relative positi...
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p, ext_kpoints)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:121
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
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.
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,...
All kind of helpful little routines.
Definition util.F:14
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...
Represent a complex full matrix.
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...
stores all the informations relevant to an mpi environment
contained for different pw related things
environment for the poisson solver
to create arrays of pools
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Ground state molecular orbitals.