(git:98357aa)
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) THEN
839 CALL recv_handlers(iproc)%wait()
840 END IF
841 END DO
842
843 ! compute total number of non-negligible excitations
844 nexcs = 0
845 DO iproc = 1, para_env%num_pe
846 nexcs = nexcs + nexcs_recv(iproc)
847 END DO
848
849 ! receive indices and amplitudes of selected excitations
850 ALLOCATE (weights_recv(nexcs), weights_neg_abs_recv(nexcs))
851 ALLOCATE (inds_recv(nexcs), inds(nexcs))
852
853 nmo_virt_occ = 0
854 DO iproc = 1, para_env%num_pe
855 IF (nexcs_recv(iproc) > 0) THEN
856 IF (iproc - 1 /= para_env%mepos) THEN
857 ! excitation amplitudes
858 CALL para_env%irecv(weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), &
859 iproc - 1, recv_handlers(iproc), 1)
860 ! compressed indices
861 CALL para_env%irecv(inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), &
862 iproc - 1, recv_handlers2(iproc), 2)
863 ELSE
864 ! data on master node
865 weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = weights_local(1:nexcs_recv(iproc))
866 inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = inds_local(1:nexcs_recv(iproc))
867 END IF
868
869 nmo_virt_occ = nmo_virt_occ + nexcs_recv(iproc)
870 END IF
871 END DO
872
873 DO iproc = 1, para_env%num_pe
874 IF (iproc - 1 /= para_env%mepos .AND. nexcs_recv(iproc) > 0) THEN
875 CALL recv_handlers(iproc)%wait()
876 CALL recv_handlers2(iproc)%wait()
877 END IF
878 END DO
879
880 DEALLOCATE (nexcs_recv, recv_handlers, recv_handlers2)
881 ELSE
882 ! working node: send the number of selected excited states to the master node
883 nexcs_send(1) = nexcs_local
884 CALL para_env%isend(nexcs_send, para_env%source, send_handler, 0)
885 CALL send_handler%wait()
886
887 IF (nexcs_local > 0) THEN
888 ! send excitation amplitudes
889 CALL para_env%isend(weights_local(1:nexcs_local), para_env%source, send_handler, 1)
890 ! send compressed indices
891 CALL para_env%isend(inds_local(1:nexcs_local), para_env%source, send_handler2, 2)
892
893 CALL send_handler%wait()
894 CALL send_handler2%wait()
895 END IF
896 END IF
897
898 ! sort non-negligible excitations on the master node according to their amplitudes,
899 ! uncompress indices and print summary information
900 IF (para_env%is_source() .AND. log_unit > 0) THEN
901 weights_neg_abs_recv(:) = -abs(weights_recv)
902 CALL sort(weights_neg_abs_recv, int(nexcs), inds)
903
904 WRITE (log_unit, '(T7,I8,F10.5,A)') istate, evals(istate)*evolt, " eV"
905
906 ! This reinitialization is needed to prevent the intel fortran compiler from introduce
907 ! a bug when using optimization level 3 flag
908 state_spin = 1
909 state_spin2 = 1
910 IF (spinflip /= no_sf_tddfpt) THEN
911 state_spin = 1
912 state_spin2 = 2
913 END IF
914 DO iexc = 1, nexcs
915 ind = inds_recv(inds(iexc)) - 1
916 IF ((nspins > 1) .AND. (spinflip == no_sf_tddfpt)) THEN
917 IF (ind < nmo_virt_occ_alpha) THEN
918 state_spin = 1
919 state_spin2 = 1
920 spin_label = '(alp)'
921 spin_label2 = '(alp)'
922 ELSE
923 state_spin = 2
924 state_spin2 = 2
925 ind = ind - nmo_virt_occ_alpha
926 spin_label = '(bet)'
927 spin_label2 = '(bet)'
928 END IF
929 END IF
930 imo_act = ind/nmo_virt8(state_spin2) + 1
931 imo_occ = gs_mos(state_spin)%index_active(imo_act)
932 imo_virt = mod(ind, nmo_virt8(state_spin2)) + 1
933
934 WRITE (log_unit, '(T27,I8,1X,A5,T48,I8,1X,A5,T70,F9.6)') imo_occ, spin_label, &
935 nmo_occ8(state_spin2) + imo_virt, spin_label2, weights_recv(inds(iexc))
936 END DO
937 END IF
938
939 ! deallocate temporary arrays
940 IF (para_env%is_source()) THEN
941 DEALLOCATE (weights_recv, weights_neg_abs_recv, inds_recv, inds)
942 END IF
943 END DO
944
945 DEALLOCATE (weights_local, inds_local)
946 IF (log_unit > 0) THEN
947 WRITE (log_unit, "(1X,A)") &
948 "-------------------------------------------------------------------------------"
949 END IF
950 END IF
951
952 CALL cp_fm_release(weights_fm)
953 CALL cp_fm_release(s_mos_virt)
954
955 CALL timestop(handle)
956
958
959! **************************************************************************************************
960!> \brief Print natural transition orbital analysis.
961!> \param qs_env Information on Kinds and Particles
962!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
963!> SIZE(evects,2) -- number of excited states to print)
964!> \param evals TDDFPT eigenvalues
965!> \param ostrength ...
966!> \param gs_mos molecular orbitals optimised for the ground state
967!> \param matrix_s overlap matrix
968!> \param print_section ...
969!> \par History
970!> * 06.2019 created [JGH]
971! **************************************************************************************************
972 SUBROUTINE tddfpt_print_nto_analysis(qs_env, evects, evals, ostrength, gs_mos, matrix_s, print_section)
973 TYPE(qs_environment_type), POINTER :: qs_env
974 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
975 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals, ostrength
976 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
977 INTENT(in) :: gs_mos
978 TYPE(dbcsr_type), POINTER :: matrix_s
979 TYPE(section_vals_type), POINTER :: print_section
980
981 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_nto_analysis'
982 INTEGER, PARAMETER :: ntomax = 10
983
984 CHARACTER(LEN=20), DIMENSION(2) :: nto_name
985 INTEGER :: handle, i, ia, icg, iounit, ispin, &
986 istate, j, nao, nlist, nmax, nmo, &
987 nnto, nspins, nstates
988 INTEGER, DIMENSION(2) :: iv
989 INTEGER, DIMENSION(2, ntomax) :: ia_index
990 INTEGER, DIMENSION(:), POINTER :: slist, stride
991 LOGICAL :: append_cube, cube_file, explicit
992 REAL(kind=dp) :: os_threshold, sume, threshold
993 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvals
994 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvalues
995 REAL(kind=dp), DIMENSION(ntomax) :: ia_eval
996 TYPE(cell_type), POINTER :: cell
997 TYPE(cp_fm_struct_type), POINTER :: fm_mo_struct, fm_struct
998 TYPE(cp_fm_type) :: sev, smat, tmat, wmat, work, wvec
999 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: teig
1000 TYPE(cp_logger_type), POINTER :: logger
1001 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: nto_set
1002 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1003 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1004 TYPE(section_vals_type), POINTER :: molden_section, nto_section
1005
1006 CALL timeset(routinen, handle)
1007
1008 logger => cp_get_default_logger()
1009 iounit = cp_logger_get_default_io_unit(logger)
1010
1011 IF (btest(cp_print_key_should_output(logger%iter_info, print_section, &
1012 "NTO_ANALYSIS"), cp_p_file)) THEN
1013
1014 CALL cite_reference(martin2003)
1015
1016 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%THRESHOLD", r_val=threshold)
1017 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%INTENSITY_THRESHOLD", r_val=os_threshold)
1018 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", explicit=explicit)
1019
1020 IF (explicit) THEN
1021 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", i_vals=slist)
1022 nlist = SIZE(slist)
1023 ELSE
1024 nlist = 0
1025 END IF
1026
1027 IF (iounit > 0) THEN
1028 WRITE (iounit, "(1X,A)") "", &
1029 "-------------------------------------------------------------------------------", &
1030 "- Natural Orbital analysis -", &
1031 "-------------------------------------------------------------------------------"
1032 END IF
1033
1034 nspins = SIZE(evects, 1)
1035 nstates = SIZE(evects, 2)
1036 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1037
1038 DO istate = 1, nstates
1039 IF (os_threshold > ostrength(istate)) THEN
1040 IF (iounit > 0) THEN
1041 WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1042 END IF
1043 cycle
1044 END IF
1045 IF (nlist > 0) THEN
1046 IF (.NOT. any(slist == istate)) THEN
1047 IF (iounit > 0) THEN
1048 WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1049 END IF
1050 cycle
1051 END IF
1052 END IF
1053 IF (iounit > 0) THEN
1054 WRITE (iounit, "(1X,A,I6,T30,F10.5,A)") " STATE NR. ", istate, evals(istate)*evolt, " eV"
1055 END IF
1056 nmax = 0
1057 DO ispin = 1, nspins
1058 CALL cp_fm_get_info(evects(ispin, istate), matrix_struct=fm_struct, ncol_global=nmo)
1059 nmax = max(nmax, nmo)
1060 END DO
1061 ALLOCATE (eigenvalues(nmax, nspins))
1062 eigenvalues = 0.0_dp
1063 ! SET 1: Hole states
1064 ! SET 2: Particle states
1065 nto_name(1) = 'Hole_states'
1066 nto_name(2) = 'Particle_states'
1067 ALLOCATE (nto_set(2))
1068 DO i = 1, 2
1069 CALL allocate_mo_set(nto_set(i), nao, ntomax, 0, 0.0_dp, 1.0_dp, 0.0_dp)
1070 CALL cp_fm_get_info(evects(1, istate), matrix_struct=fm_struct)
1071 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1072 ncol_global=ntomax)
1073 CALL cp_fm_create(tmat, fm_mo_struct)
1074 CALL init_mo_set(nto_set(i), fm_ref=tmat, name=nto_name(i))
1075 CALL cp_fm_release(tmat)
1076 CALL cp_fm_struct_release(fm_mo_struct)
1077 END DO
1078 !
1079 ALLOCATE (teig(nspins))
1080 ! hole states
1081 ! Diagonalize X(T)*S*X
1082 DO ispin = 1, nspins
1083 associate(ev => evects(ispin, istate))
1084 CALL cp_fm_get_info(ev, matrix_struct=fm_struct, ncol_global=nmo)
1085 CALL cp_fm_create(sev, fm_struct)
1086 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1087 nrow_global=nmo, ncol_global=nmo)
1088 CALL cp_fm_create(tmat, fm_mo_struct)
1089 CALL cp_fm_create(teig(ispin), fm_mo_struct)
1090 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1091 CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, ev, sev, 0.0_dp, tmat)
1092 END associate
1093
1094 CALL choose_eigv_solver(tmat, teig(ispin), eigenvalues(1:nmo, ispin))
1095
1096 CALL cp_fm_struct_release(fm_mo_struct)
1097 CALL cp_fm_release(tmat)
1098 CALL cp_fm_release(sev)
1099 END DO
1100 ! find major determinants i->a
1101 ia_index = 0
1102 sume = 0.0_dp
1103 nnto = 0
1104 DO i = 1, ntomax
1105 iv = maxloc(eigenvalues)
1106 ia_eval(i) = eigenvalues(iv(1), iv(2))
1107 ia_index(1:2, i) = iv(1:2)
1108 sume = sume + ia_eval(i)
1109 eigenvalues(iv(1), iv(2)) = 0.0_dp
1110 nnto = nnto + 1
1111 IF (sume > threshold) EXIT
1112 END DO
1113 ! store hole states
1114 CALL set_mo_set(nto_set(1), nmo=nnto)
1115 DO i = 1, nnto
1116 ia = ia_index(1, i)
1117 ispin = ia_index(2, i)
1118 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=nmo)
1119 CALL cp_fm_get_info(teig(ispin), matrix_struct=fm_struct)
1120 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1121 nrow_global=nmo, ncol_global=1)
1122 CALL cp_fm_create(tmat, fm_mo_struct)
1123 CALL cp_fm_struct_release(fm_mo_struct)
1124 CALL cp_fm_get_info(gs_mos(1)%mos_occ, matrix_struct=fm_struct)
1125 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1126 ncol_global=1)
1127 CALL cp_fm_create(wvec, fm_mo_struct)
1128 CALL cp_fm_struct_release(fm_mo_struct)
1129 CALL cp_fm_to_fm(teig(ispin), tmat, 1, ia, 1)
1130 CALL parallel_gemm('N', 'N', nao, 1, nmo, 1.0_dp, gs_mos(ispin)%mos_occ, &
1131 tmat, 0.0_dp, wvec)
1132 CALL cp_fm_to_fm(wvec, nto_set(1)%mo_coeff, 1, 1, i)
1133 CALL cp_fm_release(wvec)
1134 CALL cp_fm_release(tmat)
1135 END DO
1136 ! particle states
1137 ! Solve generalized eigenvalue equation: (S*X)*(S*X)(T)*v = lambda*S*v
1138 CALL set_mo_set(nto_set(2), nmo=nnto)
1139 DO ispin = 1, nspins
1140 associate(ev => evects(ispin, istate))
1141 CALL cp_fm_get_info(ev, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
1142 ALLOCATE (eigvals(nao))
1143 eigvals = 0.0_dp
1144 CALL cp_fm_create(sev, fm_struct)
1145 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1146 END associate
1147 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1148 nrow_global=nao, ncol_global=nao)
1149 CALL cp_fm_create(tmat, fm_mo_struct)
1150 CALL cp_fm_create(smat, fm_mo_struct)
1151 CALL cp_fm_create(wmat, fm_mo_struct)
1152 CALL cp_fm_create(work, fm_mo_struct)
1153 CALL cp_fm_struct_release(fm_mo_struct)
1154 CALL copy_dbcsr_to_fm(matrix_s, smat)
1155 CALL parallel_gemm('N', 'T', nao, nao, nmo, 1.0_dp, sev, sev, 0.0_dp, tmat)
1156 CALL cp_fm_geeig(tmat, smat, wmat, eigvals, work)
1157 DO i = 1, nnto
1158 IF (ispin == ia_index(2, i)) THEN
1159 icg = 0
1160 DO j = 1, nao
1161 IF (abs(eigvals(j) - ia_eval(i)) < 1.e-6_dp) THEN
1162 icg = j
1163 EXIT
1164 END IF
1165 END DO
1166 IF (icg == 0) THEN
1167 CALL cp_warn(__location__, &
1168 "Could not locate particle state associated with hole state.")
1169 ELSE
1170 CALL cp_fm_to_fm(wmat, nto_set(2)%mo_coeff, 1, icg, i)
1171 END IF
1172 END IF
1173 END DO
1174 DEALLOCATE (eigvals)
1175 CALL cp_fm_release(sev)
1176 CALL cp_fm_release(tmat)
1177 CALL cp_fm_release(smat)
1178 CALL cp_fm_release(wmat)
1179 CALL cp_fm_release(work)
1180 END DO
1181 ! print
1182 IF (iounit > 0) THEN
1183 sume = 0.0_dp
1184 DO i = 1, nnto
1185 sume = sume + ia_eval(i)
1186 WRITE (iounit, "(T6,A,i2,T30,A,i1,T42,A,F8.5,T63,A,F8.5)") &
1187 "Particle-Hole state:", i, " Spin:", ia_index(2, i), &
1188 "Eigenvalue:", ia_eval(i), " Sum Eigv:", sume
1189 END DO
1190 END IF
1191 ! Cube and Molden files
1192 nto_section => section_vals_get_subs_vals(print_section, "NTO_ANALYSIS")
1193 CALL section_vals_val_get(nto_section, "CUBE_FILES", l_val=cube_file)
1194 CALL section_vals_val_get(nto_section, "STRIDE", i_vals=stride)
1195 CALL section_vals_val_get(nto_section, "APPEND", l_val=append_cube)
1196 IF (cube_file) THEN
1197 CALL print_nto_cubes(qs_env, nto_set, istate, stride, append_cube, nto_section)
1198 END IF
1199 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, cell=cell)
1200 molden_section => section_vals_get_subs_vals(print_section, "MOS_MOLDEN")
1201 CALL write_mos_molden(nto_set, qs_kind_set, particle_set, molden_section, cell=cell, qs_env=qs_env)
1202 !
1203 DEALLOCATE (eigenvalues)
1204 CALL cp_fm_release(teig)
1205 !
1206 DO i = 1, 2
1207 CALL deallocate_mo_set(nto_set(i))
1208 END DO
1209 DEALLOCATE (nto_set)
1210 END DO
1211
1212 IF (iounit > 0) THEN
1213 WRITE (iounit, "(1X,A)") &
1214 "-------------------------------------------------------------------------------"
1215 END IF
1216
1217 END IF
1218
1219 CALL timestop(handle)
1220
1221 END SUBROUTINE tddfpt_print_nto_analysis
1222
1223! **************************************************************************************************
1224!> \brief Print exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018)
1225!> \param log_unit output unit
1226!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
1227!> SIZE(evects,2) -- number of excited states to print)
1228!> \param gs_mos molecular orbitals optimised for the ground state
1229!> \param matrix_s overlap matrix
1230!> \param do_directional_exciton_descriptors flag for computing descriptors for each (cartesian) direction
1231!> \param qs_env Information on particles/geometry
1232!> \par History
1233!> * 12.2024 created as 'tddfpt_print_exciton_descriptors' [Maximilian Graml]
1234! **************************************************************************************************
1235 SUBROUTINE tddfpt_print_exciton_descriptors(log_unit, evects, gs_mos, matrix_s, &
1236 do_directional_exciton_descriptors, qs_env)
1237 INTEGER, INTENT(in) :: log_unit
1238 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
1239 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1240 INTENT(in) :: gs_mos
1241 TYPE(dbcsr_type), POINTER :: matrix_s
1242 LOGICAL, INTENT(IN) :: do_directional_exciton_descriptors
1243 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1244
1245 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_exciton_descriptors'
1246
1247 CHARACTER(LEN=4) :: prefix_output
1248 INTEGER :: handle, ispin, istate, n_moments_quad, &
1249 nactive, nao, nspins, nstates
1250 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt
1251 LOGICAL :: print_checkvalue
1252 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ref_point_multipole
1253 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1254 TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo_coeff, &
1255 fm_struct_s_mos_virt, fm_struct_x_ia_n
1256 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: eigvec_x_ia_n, fm_multipole_ab, &
1257 fm_multipole_ai, fm_multipole_ij, &
1258 s_mos_virt
1259 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_coeff
1260 TYPE(exciton_descr_type), ALLOCATABLE, &
1261 DIMENSION(:) :: exc_descr
1262
1263 CALL timeset(routinen, handle)
1264
1265 nspins = SIZE(evects, 1)
1266 nstates = SIZE(evects, 2)
1267
1268 cpassert(nspins == 1) ! Other spins are not yet implemented for exciton descriptors
1269
1270 CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env)
1271 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1272
1273 DO ispin = 1, nspins
1274 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
1275 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
1276 END DO
1277
1278 ! Prepare fm with all MO coefficents, i.e. nao x nao
1279 ALLOCATE (mo_coeff(nspins))
1280 CALL cp_fm_struct_create(fm_struct_mo_coeff, nrow_global=nao, ncol_global=nao, &
1281 context=blacs_env)
1282 DO ispin = 1, nspins
1283 CALL cp_fm_create(mo_coeff(ispin), fm_struct_mo_coeff)
1284 CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_occ, &
1285 mo_coeff(ispin), &
1286 nao, &
1287 nmo_occ(ispin), &
1288 1, &
1289 1, &
1290 1, &
1291 1, &
1292 blacs_env)
1293 CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_virt, &
1294 mo_coeff(ispin), &
1295 nao, &
1296 nmo_virt(ispin), &
1297 1, &
1298 1, &
1299 1, &
1300 nmo_occ(ispin) + 1, &
1301 blacs_env)
1302 END DO
1303 CALL cp_fm_struct_release(fm_struct_mo_coeff)
1304
1305 ! Compute multipole moments
1306 ! fm_multipole_XY have structure inherited by libint, i.e. x, y, z, xx, xy, xz, yy, yz, zz
1307 n_moments_quad = 9
1308 ALLOCATE (ref_point_multipole(3))
1309 ALLOCATE (fm_multipole_ij(n_moments_quad))
1310 ALLOCATE (fm_multipole_ab(n_moments_quad))
1311 ALLOCATE (fm_multipole_ai(n_moments_quad))
1312
1313 CALL get_multipoles_mo(fm_multipole_ai, fm_multipole_ij, fm_multipole_ab, &
1314 qs_env, mo_coeff, ref_point_multipole, 2, &
1315 nmo_occ(1), nmo_virt(1), blacs_env)
1316
1317 CALL cp_fm_release(mo_coeff)
1318
1319 ! Compute eigenvector X of the Casida equation from trial vectors
1320 ALLOCATE (s_mos_virt(nspins), eigvec_x_ia_n(nspins))
1321 DO ispin = 1, nspins
1322 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct_s_mos_virt)
1323 CALL cp_fm_create(s_mos_virt(ispin), fm_struct_s_mos_virt)
1324 NULLIFY (fm_struct_s_mos_virt)
1325 CALL cp_dbcsr_sm_fm_multiply(matrix_s, &
1326 gs_mos(ispin)%mos_virt, &
1327 s_mos_virt(ispin), &
1328 ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
1329
1330 CALL cp_fm_struct_create(fm_struct_x_ia_n, nrow_global=nmo_occ(ispin), ncol_global=nmo_virt(ispin), &
1331 context=blacs_env)
1332 CALL cp_fm_create(eigvec_x_ia_n(ispin), fm_struct_x_ia_n)
1333 CALL cp_fm_struct_release(fm_struct_x_ia_n)
1334 END DO
1335 ALLOCATE (exc_descr(nstates))
1336 DO istate = 1, nstates
1337 DO ispin = 1, nspins
1338 CALL cp_fm_set_all(eigvec_x_ia_n(ispin), 0.0_dp)
1339 ! compute eigenvectors X of the TDA equation
1340 ! Reshuffle multiplication from
1341 ! X_ai = S_ma ^T * C_mi
1342 ! to
1343 ! X_ia = C_mi ^T * S_ma
1344 ! for compatibility with the structure needed for get_exciton_descriptors of bse_properties.F
1345 CALL cp_fm_get_info(evects(ispin, istate), ncol_global=nactive)
1346 IF (nactive /= nmo_occ(ispin)) THEN
1347 CALL cp_abort(__location__, &
1348 "Reduced active space excitations not implemented")
1349 END IF
1350 CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_virt(ispin), nao, 1.0_dp, &
1351 evects(ispin, istate), s_mos_virt(ispin), 0.0_dp, eigvec_x_ia_n(ispin))
1352
1353 CALL get_exciton_descriptors(exc_descr, eigvec_x_ia_n(ispin), &
1354 fm_multipole_ij, fm_multipole_ab, &
1355 fm_multipole_ai, &
1356 istate, nmo_occ(ispin), nmo_virt(ispin))
1357 END DO
1358 END DO
1359 CALL cp_fm_release(eigvec_x_ia_n)
1360 CALL cp_fm_release(s_mos_virt)
1361 CALL cp_fm_release(fm_multipole_ai)
1362 CALL cp_fm_release(fm_multipole_ij)
1363 CALL cp_fm_release(fm_multipole_ab)
1364
1365 ! Actual printing
1366 print_checkvalue = .true.
1367 prefix_output = ' '
1368 CALL print_exciton_descriptors(exc_descr, ref_point_multipole, log_unit, &
1369 nstates, print_checkvalue, do_directional_exciton_descriptors, &
1370 prefix_output, qs_env)
1371
1372 DEALLOCATE (ref_point_multipole)
1373 DEALLOCATE (exc_descr)
1374
1375 CALL timestop(handle)
1376
1378
1379! **************************************************************************************************
1380!> \brief ...
1381!> \param vin ...
1382!> \param vout ...
1383!> \param mos_occ ...
1384!> \param matrix_s ...
1385! **************************************************************************************************
1386 SUBROUTINE project_vector(vin, vout, mos_occ, matrix_s)
1387 TYPE(dbcsr_type) :: vin, vout
1388 TYPE(cp_fm_type), INTENT(IN) :: mos_occ
1389 TYPE(dbcsr_type), POINTER :: matrix_s
1390
1391 CHARACTER(LEN=*), PARAMETER :: routinen = 'project_vector'
1392
1393 INTEGER :: handle, nao, nmo
1394 REAL(kind=dp) :: norm(1)
1395 TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_vec_struct
1396 TYPE(cp_fm_type) :: csvec, svec, vec
1397
1398 CALL timeset(routinen, handle)
1399
1400 CALL cp_fm_get_info(mos_occ, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
1401 CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
1402 nrow_global=nao, ncol_global=1)
1403 CALL cp_fm_create(vec, fm_vec_struct)
1404 CALL cp_fm_create(svec, fm_vec_struct)
1405 CALL cp_fm_struct_release(fm_vec_struct)
1406 CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
1407 nrow_global=nmo, ncol_global=1)
1408 CALL cp_fm_create(csvec, fm_vec_struct)
1409 CALL cp_fm_struct_release(fm_vec_struct)
1410
1411 CALL copy_dbcsr_to_fm(vin, vec)
1412 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vec, svec, ncol=1, alpha=1.0_dp, beta=0.0_dp)
1413 CALL parallel_gemm('T', 'N', nmo, 1, nao, 1.0_dp, mos_occ, svec, 0.0_dp, csvec)
1414 CALL parallel_gemm('N', 'N', nao, 1, nmo, -1.0_dp, mos_occ, csvec, 1.0_dp, vec)
1415 CALL cp_fm_vectorsnorm(vec, norm)
1416 cpassert(norm(1) > 1.e-14_dp)
1417 norm(1) = sqrt(1._dp/norm(1))
1418 CALL cp_fm_scale(norm(1), vec)
1419 CALL copy_fm_to_dbcsr(vec, vout, keep_sparsity=.false.)
1420
1421 CALL cp_fm_release(csvec)
1422 CALL cp_fm_release(svec)
1423 CALL cp_fm_release(vec)
1424
1425 CALL timestop(handle)
1426
1427 END SUBROUTINE project_vector
1428
1429! **************************************************************************************************
1430!> \brief ...
1431!> \param va ...
1432!> \param vb ...
1433!> \param res ...
1434! **************************************************************************************************
1435 SUBROUTINE vec_product(va, vb, res)
1436 TYPE(dbcsr_type) :: va, vb
1437 REAL(kind=dp), INTENT(OUT) :: res
1438
1439 CHARACTER(LEN=*), PARAMETER :: routinen = 'vec_product'
1440
1441 INTEGER :: handle, icol, irow
1442 LOGICAL :: found
1443 REAL(kind=dp), DIMENSION(:, :), POINTER :: vba, vbb
1444 TYPE(dbcsr_iterator_type) :: iter
1445 TYPE(mp_comm_type) :: group
1446
1447 CALL timeset(routinen, handle)
1448
1449 res = 0.0_dp
1450
1451 CALL dbcsr_get_info(va, group=group)
1452 CALL dbcsr_iterator_start(iter, va)
1453 DO WHILE (dbcsr_iterator_blocks_left(iter))
1454 CALL dbcsr_iterator_next_block(iter, irow, icol, vba)
1455 CALL dbcsr_get_block_p(vb, row=irow, col=icol, block=vbb, found=found)
1456 res = res + sum(vba*vbb)
1457 cpassert(found)
1458 END DO
1459 CALL dbcsr_iterator_stop(iter)
1460 CALL group%sum(res)
1461
1462 CALL timestop(handle)
1463
1464 END SUBROUTINE vec_product
1465
1466! **************************************************************************************************
1467!> \brief ...
1468!> \param qs_env ...
1469!> \param mos ...
1470!> \param istate ...
1471!> \param stride ...
1472!> \param append_cube ...
1473!> \param print_section ...
1474! **************************************************************************************************
1475 SUBROUTINE print_nto_cubes(qs_env, mos, istate, stride, append_cube, print_section)
1476
1477 TYPE(qs_environment_type), POINTER :: qs_env
1478 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1479 INTEGER, INTENT(IN) :: istate
1480 INTEGER, DIMENSION(:), POINTER :: stride
1481 LOGICAL, INTENT(IN) :: append_cube
1482 TYPE(section_vals_type), POINTER :: print_section
1483
1484 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
1485 INTEGER :: i, iset, nmo, unit_nr
1486 LOGICAL :: mpi_io
1487 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1488 TYPE(cell_type), POINTER :: cell
1489 TYPE(cp_fm_type), POINTER :: mo_coeff
1490 TYPE(cp_logger_type), POINTER :: logger
1491 TYPE(dft_control_type), POINTER :: dft_control
1492 TYPE(particle_list_type), POINTER :: particles
1493 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1494 TYPE(pw_c1d_gs_type) :: wf_g
1495 TYPE(pw_env_type), POINTER :: pw_env
1496 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1497 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1498 TYPE(pw_r3d_rs_type) :: wf_r
1499 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1500 TYPE(qs_subsys_type), POINTER :: subsys
1501
1502 logger => cp_get_default_logger()
1503
1504 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
1505 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1506 CALL auxbas_pw_pool%create_pw(wf_r)
1507 CALL auxbas_pw_pool%create_pw(wf_g)
1508
1509 CALL get_qs_env(qs_env, subsys=subsys)
1510 CALL qs_subsys_get(subsys, particles=particles)
1511
1512 my_pos_cube = "REWIND"
1513 IF (append_cube) THEN
1514 my_pos_cube = "APPEND"
1515 END IF
1516
1517 CALL get_qs_env(qs_env=qs_env, &
1518 atomic_kind_set=atomic_kind_set, &
1519 qs_kind_set=qs_kind_set, &
1520 cell=cell, &
1521 particle_set=particle_set)
1522
1523 DO iset = 1, 2
1524 CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo)
1525 DO i = 1, nmo
1526 CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1527 cell, dft_control, particle_set, pw_env)
1528 IF (iset == 1) THEN
1529 WRITE (filename, '(a4,I3.3,I2.2,a11)') "NTO_STATE", istate, i, "_Hole_State"
1530 ELSE IF (iset == 2) THEN
1531 WRITE (filename, '(a4,I3.3,I2.2,a15)') "NTO_STATE", istate, i, "_Particle_State"
1532 END IF
1533 mpi_io = .true.
1534 unit_nr = cp_print_key_unit_nr(logger, print_section, '', extension=".cube", &
1535 middle_name=trim(filename), file_position=my_pos_cube, &
1536 log_filename=.false., ignore_should_output=.true., mpi_io=mpi_io)
1537 IF (iset == 1) THEN
1538 WRITE (title, *) "Natural Transition Orbital Hole State", i
1539 ELSE IF (iset == 2) THEN
1540 WRITE (title, *) "Natural Transition Orbital Particle State", i
1541 END IF
1542 CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, stride=stride, mpi_io=mpi_io)
1543 CALL cp_print_key_finished_output(unit_nr, logger, print_section, '', &
1544 ignore_should_output=.true., mpi_io=mpi_io)
1545 END DO
1546 END DO
1547
1548 CALL auxbas_pw_pool%give_back_pw(wf_g)
1549 CALL auxbas_pw_pool%give_back_pw(wf_r)
1550
1551 END SUBROUTINE print_nto_cubes
1552
1553END 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:573
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, ispin)
...
Definition bse_util.F:1886
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_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
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
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, matrix_vhxc, 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.