(git:5557c94)
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,&
68 USE kinds, ONLY: default_path_length,&
69 dp,&
70 int_8
71 USE mathconstants, ONLY: twopi,&
72 z_one,&
73 z_zero
74 USE message_passing, ONLY: mp_comm_type,&
82 USE physcon, ONLY: evolt
83 USE pw_env_types, ONLY: pw_env_get,&
86 USE pw_pool_types, ONLY: pw_pool_p_type,&
88 USE pw_types, ONLY: pw_c1d_gs_type,&
95 USE qs_mo_types, ONLY: allocate_mo_set,&
103 USE qs_operators_ao, ONLY: rrc_xyz_ao
105 USE qs_subsys_types, ONLY: qs_subsys_get,&
109 USE util, ONLY: sort
110#include "./base/base_uses.f90"
111
112 IMPLICIT NONE
113
114 PRIVATE
115
116 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_properties'
117
118 ! number of first derivative components (3: d/dx, d/dy, d/dz)
119 INTEGER, PARAMETER, PRIVATE :: nderivs = 3
120 INTEGER, PARAMETER, PRIVATE :: maxspins = 2
121
124
125! **************************************************************************************************
126
127CONTAINS
128
129! **************************************************************************************************
130!> \brief Compute the action of the dipole operator on the ground state wave function.
131!> \param dipole_op_mos_occ 2-D array [x,y,z ; spin] of matrices where to put the computed quantity
132!> (allocated and initialised on exit)
133!> \param tddfpt_control TDDFPT control parameters
134!> \param gs_mos molecular orbitals optimised for the ground state
135!> \param qs_env Quickstep environment
136!> \par History
137!> * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov]
138!> * 06.2018 dipole operator based on the Berry-phase formula [Sergey Chulkov]
139!> * 08.2018 splited of from 'tddfpt_print_summary' and merged with code from 'tddfpt'
140!> [Sergey Chulkov]
141!> \note \parblock
142!> Adapted version of the subroutine find_contributions() which was originally created
143!> by Thomas Chassaing on 02.2005.
144!>
145!> The relation between dipole integrals in velocity and length forms are the following:
146!> \f[<\psi_i|\nabla|\psi_a> = <\psi_i|\vec{r}|\hat{H}\psi_a> - <\hat{H}\psi_i|\vec{r}|\psi_a>
147!> = (\epsilon_a - \epsilon_i) <\psi_i|\vec{r}|\psi_a> .\f],
148!> due to the commutation identity:
149!> \f[\vec{r}\hat{H} - \hat{H}\vec{r} = [\vec{r},\hat{H}] = [\vec{r},-1/2 \nabla^2] = \nabla\f] .
150!> \endparblock
151! **************************************************************************************************
152 SUBROUTINE tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_env)
153 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :), &
154 INTENT(inout) :: dipole_op_mos_occ
155 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
156 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
157 INTENT(in) :: gs_mos
158 TYPE(qs_environment_type), POINTER :: qs_env
159
160 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_dipole_operator'
161
162 INTEGER :: handle, i_cos_sin, icol, ideriv, irow, &
163 ispin, jderiv, nao, ncols_local, &
164 ndim_periodic, nrows_local, nspins
165 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
166 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt
167 REAL(kind=dp) :: eval_occ
168 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
169 POINTER :: local_data_ediff, local_data_wfm
170 REAL(kind=dp), DIMENSION(3) :: kvec, reference_point
171 TYPE(cell_type), POINTER :: cell
172 TYPE(cp_blacs_env_type), POINTER :: blacs_env
173 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: gamma_00, gamma_inv_00
174 TYPE(cp_fm_struct_type), POINTER :: fm_struct
175 TYPE(cp_fm_type) :: ediff_inv, rrc_mos_occ, wfm_ao_ao, &
176 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_berry
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) is needed "// &
243 "for oscillator strengths based on the Berry phase formula")
244 END IF
245
246 NULLIFY (berry_cossin_xyz)
247 ! index: 1 = Re[exp(-i * G_t * t)],
248 ! 2 = Im[exp(-i * G_t * t)];
249 ! t = x,y,z
250 CALL dbcsr_allocate_matrix_set(berry_cossin_xyz, 2)
251
252 DO i_cos_sin = 1, 2
253 CALL dbcsr_init_p(berry_cossin_xyz(i_cos_sin)%matrix)
254 CALL dbcsr_copy(berry_cossin_xyz(i_cos_sin)%matrix, matrix_s(1)%matrix)
255 END DO
256
257 ! +++ allocate berry-phase-related work matrices
258 ALLOCATE (gamma_00(nspins), gamma_inv_00(nspins), gamma_real_imag(2, nspins), opvec(2, nspins))
259 ALLOCATE (dberry_mos_occ(nderivs, nspins))
260 DO ispin = 1, nspins
261 NULLIFY (fm_struct)
262 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_occ(ispin), &
263 ncol_global=nmo_occ(ispin), context=blacs_env)
264
265 CALL cp_cfm_create(gamma_00(ispin), fm_struct)
266 CALL cp_cfm_create(gamma_inv_00(ispin), fm_struct)
267
268 DO i_cos_sin = 1, 2
269 CALL cp_fm_create(gamma_real_imag(i_cos_sin, ispin), fm_struct)
270 END DO
271 CALL cp_fm_struct_release(fm_struct)
272
273 ! G_real C_0, G_imag C_0
274 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct)
275 DO i_cos_sin = 1, 2
276 CALL cp_fm_create(opvec(i_cos_sin, ispin), fm_struct)
277 END DO
278
279 ! dBerry * C_0
280 DO ideriv = 1, nderivs
281 CALL cp_fm_create(dberry_mos_occ(ideriv, ispin), fm_struct)
282 CALL cp_fm_set_all(dberry_mos_occ(ideriv, ispin), 0.0_dp)
283 END DO
284 END DO
285
286 DO ideriv = 1, nderivs
287 kvec(:) = twopi*cell%h_inv(ideriv, :)
288 DO i_cos_sin = 1, 2
289 CALL dbcsr_set(berry_cossin_xyz(i_cos_sin)%matrix, 0.0_dp)
290 END DO
291 CALL build_berry_moment_matrix(qs_env, berry_cossin_xyz(1)%matrix, &
292 berry_cossin_xyz(2)%matrix, kvec)
293
294 DO ispin = 1, nspins
295 ! i_cos_sin = 1: cos (real) component; opvec(1) = gamma_real C_0
296 ! i_cos_sin = 2: sin (imaginary) component; opvec(2) = gamma_imag C_0
297 DO i_cos_sin = 1, 2
298 CALL cp_dbcsr_sm_fm_multiply(berry_cossin_xyz(i_cos_sin)%matrix, &
299 gs_mos(ispin)%mos_occ, &
300 opvec(i_cos_sin, ispin), &
301 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
302 END DO
303
304 CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, &
305 1.0_dp, gs_mos(ispin)%mos_occ, opvec(1, ispin), &
306 0.0_dp, gamma_real_imag(1, ispin))
307
308 CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, &
309 -1.0_dp, gs_mos(ispin)%mos_occ, opvec(2, ispin), &
310 0.0_dp, gamma_real_imag(2, ispin))
311
312 CALL cp_fm_to_cfm(msourcer=gamma_real_imag(1, ispin), &
313 msourcei=gamma_real_imag(2, ispin), &
314 mtarget=gamma_00(ispin))
315
316 ! gamma_inv_00 = Q = [C_0^T (gamma_real - i gamma_imag) C_0] ^ {-1}
317 CALL cp_cfm_set_all(gamma_inv_00(ispin), z_zero, z_one)
318 CALL cp_cfm_solve(gamma_00(ispin), gamma_inv_00(ispin))
319
320 CALL cp_cfm_to_fm(msource=gamma_inv_00(ispin), &
321 mtargetr=gamma_real_imag(1, ispin), &
322 mtargeti=gamma_real_imag(2, ispin))
323
324 ! dBerry_mos_occ is identical to dBerry_psi0 from qs_linres_op % polar_operators()
325 CALL parallel_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), &
326 1.0_dp, opvec(1, ispin), gamma_real_imag(2, ispin), &
327 0.0_dp, dipole_op_mos_occ(1, ispin))
328 CALL parallel_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), &
329 -1.0_dp, opvec(2, ispin), gamma_real_imag(1, ispin), &
330 1.0_dp, dipole_op_mos_occ(1, ispin))
331
332 DO jderiv = 1, nderivs
333 CALL cp_fm_scale_and_add(1.0_dp, dberry_mos_occ(jderiv, ispin), &
334 cell%hmat(jderiv, ideriv), dipole_op_mos_occ(1, ispin))
335 END DO
336 END DO
337 END DO
338
339 ! --- release berry-phase-related work matrices
340 CALL cp_fm_release(opvec)
341 CALL cp_fm_release(gamma_real_imag)
342 DO ispin = nspins, 1, -1
343 CALL cp_cfm_release(gamma_inv_00(ispin))
344 CALL cp_cfm_release(gamma_00(ispin))
345 END DO
346 DEALLOCATE (gamma_00, gamma_inv_00)
347 CALL dbcsr_deallocate_matrix_set(berry_cossin_xyz)
348
349 NULLIFY (fm_struct)
350 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
351 CALL cp_fm_create(wfm_ao_ao, fm_struct)
352 CALL cp_fm_struct_release(fm_struct)
353
354 ! trans_dipole = 2|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00) +
355 ! 2|e|/|G_mu| * Tr Imag(C_0^T * (gamma_real - i gamma_imag) * evects * gamma_inv_00) ,
356 !
357 ! Taking into account the symmetry of the matrices 'gamma_real' and 'gamma_imag' and the fact
358 ! that the response wave-function is a real-valued function, the above expression can be simplified as
359 ! trans_dipole = 4|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00)
360 !
361 ! 1/|G_mu| = |lattice_vector_mu| / (2*pi) .
362 DO ispin = 1, nspins
363 ! wfm_ao_ao = S * mos_virt * mos_virt^T
364 CALL parallel_gemm('N', 'T', nao, nao, nmo_virt(ispin), &
365 1.0_dp/twopi, s_mos_virt(ispin), gs_mos(ispin)%mos_virt, &
366 0.0_dp, wfm_ao_ao)
367
368 DO ideriv = 1, nderivs
369 CALL parallel_gemm('N', 'N', nao, nmo_occ(ispin), nao, &
370 1.0_dp, wfm_ao_ao, dberry_mos_occ(ideriv, ispin), &
371 0.0_dp, dipole_op_mos_occ(ideriv, ispin))
372 END DO
373 END DO
374
375 CALL cp_fm_release(wfm_ao_ao)
376 CALL cp_fm_release(dberry_mos_occ)
377
379 IF (ndim_periodic /= 0) THEN
380 CALL cp_warn(__location__, &
381 "Non-periodic Poisson solver (PERIODIC none) is needed "// &
382 "for oscillator strengths based on the length operator")
383 END IF
384
385 ! compute components of the dipole operator in the length form
386 NULLIFY (rrc_xyz)
387 CALL dbcsr_allocate_matrix_set(rrc_xyz, nderivs)
388
389 DO ideriv = 1, nderivs
390 CALL dbcsr_init_p(rrc_xyz(ideriv)%matrix)
391 CALL dbcsr_copy(rrc_xyz(ideriv)%matrix, matrix_s(1)%matrix)
392 END DO
393
394 CALL get_reference_point(reference_point, qs_env=qs_env, &
395 reference=tddfpt_control%dipole_reference, &
396 ref_point=tddfpt_control%dipole_ref_point)
397
398 CALL rrc_xyz_ao(op=rrc_xyz, qs_env=qs_env, rc=reference_point, order=1, &
399 minimum_image=.false., soft=.false.)
400
401 NULLIFY (fm_struct)
402 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
403 CALL cp_fm_create(wfm_ao_ao, fm_struct)
404 CALL cp_fm_struct_release(fm_struct)
405
406 DO ispin = 1, nspins
407 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct)
408 CALL cp_fm_create(rrc_mos_occ, fm_struct)
409
410 ! wfm_ao_ao = S * mos_virt * mos_virt^T
411 CALL parallel_gemm('N', 'T', nao, nao, nmo_virt(ispin), &
412 1.0_dp, s_mos_virt(ispin), gs_mos(ispin)%mos_virt, &
413 0.0_dp, wfm_ao_ao)
414
415 DO ideriv = 1, nderivs
416 CALL cp_dbcsr_sm_fm_multiply(rrc_xyz(ideriv)%matrix, &
417 gs_mos(ispin)%mos_occ, &
418 rrc_mos_occ, &
419 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
420
421 CALL parallel_gemm('N', 'N', nao, nmo_occ(ispin), nao, &
422 1.0_dp, wfm_ao_ao, rrc_mos_occ, &
423 0.0_dp, dipole_op_mos_occ(ideriv, ispin))
424 END DO
425
426 CALL cp_fm_release(rrc_mos_occ)
427 END DO
428
429 CALL cp_fm_release(wfm_ao_ao)
430 CALL dbcsr_deallocate_matrix_set(rrc_xyz)
431
433 ! generate overlap derivatives
434 CALL get_qs_env(qs_env, ks_env=ks_env, sab_orb=sab_orb)
435 NULLIFY (scrm)
436 CALL build_overlap_matrix(ks_env, matrix_s=scrm, nderivative=1, &
437 basis_type_a="ORB", basis_type_b="ORB", &
438 sab_nl=sab_orb)
439
440 DO ispin = 1, nspins
441 NULLIFY (fm_struct)
442 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), &
443 ncol_global=nmo_occ(ispin), context=blacs_env)
444 CALL cp_fm_create(ediff_inv, fm_struct)
445 CALL cp_fm_create(wfm_mo_virt_mo_occ, fm_struct)
446 CALL cp_fm_struct_release(fm_struct)
447
448 CALL cp_fm_get_info(ediff_inv, nrow_local=nrows_local, ncol_local=ncols_local, &
449 row_indices=row_indices, col_indices=col_indices, local_data=local_data_ediff)
450 CALL cp_fm_get_info(wfm_mo_virt_mo_occ, local_data=local_data_wfm)
451
452!$OMP PARALLEL DO DEFAULT(NONE), &
453!$OMP PRIVATE(eval_occ, icol, irow), &
454!$OMP SHARED(col_indices, gs_mos, ispin, local_data_ediff, ncols_local, nrows_local, row_indices)
455 DO icol = 1, ncols_local
456 ! E_occ_i ; imo_occ = col_indices(icol)
457 eval_occ = gs_mos(ispin)%evals_occ(col_indices(icol))
458
459 DO irow = 1, nrows_local
460 ! ediff_inv_weights(a, i) = 1.0 / (E_virt_a - E_occ_i)
461 ! imo_virt = row_indices(irow)
462 local_data_ediff(irow, icol) = 1.0_dp/(gs_mos(ispin)%evals_virt(row_indices(irow)) - eval_occ)
463 END DO
464 END DO
465!$OMP END PARALLEL DO
466
467 DO ideriv = 1, nderivs
468 CALL cp_dbcsr_sm_fm_multiply(scrm(ideriv + 1)%matrix, &
469 gs_mos(ispin)%mos_occ, &
470 dipole_op_mos_occ(ideriv, ispin), &
471 ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
472
473 CALL parallel_gemm('T', 'N', nmo_virt(ispin), nmo_occ(ispin), nao, &
474 1.0_dp, gs_mos(ispin)%mos_virt, dipole_op_mos_occ(ideriv, ispin), &
475 0.0_dp, wfm_mo_virt_mo_occ)
476
477 ! in-place element-wise (Schur) product;
478 ! avoid allocation of a temporary [nmo_virt x nmo_occ] matrix which is needed
479 ! for cp_fm_schur_product() subroutine call
480
481!$OMP PARALLEL DO DEFAULT(NONE), &
482!$OMP PRIVATE(icol, irow), &
483!$OMP SHARED(ispin, local_data_ediff, local_data_wfm, ncols_local, nrows_local)
484 DO icol = 1, ncols_local
485 DO irow = 1, nrows_local
486 local_data_wfm(irow, icol) = local_data_wfm(irow, icol)*local_data_ediff(irow, icol)
487 END DO
488 END DO
489!$OMP END PARALLEL DO
490
491 CALL parallel_gemm('N', 'N', nao, nmo_occ(ispin), nmo_virt(ispin), &
492 1.0_dp, s_mos_virt(ispin), wfm_mo_virt_mo_occ, &
493 0.0_dp, dipole_op_mos_occ(ideriv, ispin))
494 END DO
495
496 CALL cp_fm_release(wfm_mo_virt_mo_occ)
497 CALL cp_fm_release(ediff_inv)
498 END DO
500
501 CASE DEFAULT
502 cpabort("Unimplemented form of the dipole operator")
503 END SELECT
504
505 ! --- release work matrices
506 CALL cp_fm_release(s_mos_virt)
507
508 CALL timestop(handle)
509 END SUBROUTINE tddfpt_dipole_operator
510
511! **************************************************************************************************
512!> \brief Print final TDDFPT excitation energies and oscillator strengths.
513!> \param log_unit output unit
514!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
515!> SIZE(evects,2) -- number of excited states to print)
516!> \param evals TDDFPT eigenvalues
517!> \param gs_mos ...
518!> \param ostrength TDDFPT oscillator strength
519!> \param mult multiplicity
520!> \param dipole_op_mos_occ action of the dipole operator on the ground state wave function
521!> [x,y,z ; spin]
522!> \param dipole_form ...
523!> \par History
524!> * 05.2016 created [Sergey Chulkov]
525!> * 06.2016 transition dipole moments and oscillator strengths [Sergey Chulkov]
526!> * 07.2016 spin-unpolarised electron density [Sergey Chulkov]
527!> * 08.2018 compute 'dipole_op_mos_occ' in a separate subroutine [Sergey Chulkov]
528!> \note \parblock
529!> Adapted version of the subroutine find_contributions() which was originally created
530!> by Thomas Chassaing on 02.2005.
531!>
532!> Transition dipole moment along direction 'd' is computed as following:
533!> \f[ t_d(spin) = Tr[evects^T dipole\_op\_mos\_occ(d, spin)] .\f]
534!> \endparblock
535! **************************************************************************************************
536 SUBROUTINE tddfpt_print_summary(log_unit, evects, evals, gs_mos, ostrength, mult, &
537 dipole_op_mos_occ, dipole_form)
538 INTEGER, INTENT(in) :: log_unit
539 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
540 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
541 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
542 POINTER :: gs_mos
543 REAL(kind=dp), DIMENSION(:), INTENT(inout) :: ostrength
544 INTEGER, INTENT(in) :: mult
545 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: dipole_op_mos_occ
546 INTEGER, INTENT(in) :: dipole_form
547
548 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_summary'
549
550 CHARACTER(len=1) :: lsd_str
551 CHARACTER(len=20) :: mult_str
552 INTEGER :: handle, i, ideriv, ispin, istate, j, &
553 nactive, nao, nocc, nspins, nstates
554 REAL(kind=dp) :: osc_strength
555 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: trans_dipoles
556 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
557 TYPE(cp_fm_type) :: dipact
558
559 CALL timeset(routinen, handle)
560
561 nspins = SIZE(evects, 1)
562 nstates = SIZE(evects, 2)
563
564 IF (nspins > 1) THEN
565 lsd_str = 'U'
566 ELSE
567 lsd_str = 'R'
568 END IF
569
570 ! *** summary header ***
571 IF (log_unit > 0) THEN
572 CALL integer_to_string(mult, mult_str)
573 WRITE (log_unit, '(/,1X,A1,A,1X,A)') lsd_str, "-TDDFPT states of multiplicity", trim(mult_str)
574 SELECT CASE (dipole_form)
576 WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using Berry operator formulation"
578 WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using length formulation"
580 WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using velocity formulation"
581 CASE DEFAULT
582 cpabort("Unimplemented form of the dipole operator")
583 END SELECT
584
585 WRITE (log_unit, '(T10,A,T19,A,T37,A,T69,A)') "State", "Excitation", &
586 "Transition dipole (a.u.)", "Oscillator"
587 WRITE (log_unit, '(T10,A,T19,A,T37,A,T49,A,T61,A,T67,A)') "number", "energy (eV)", &
588 "x", "y", "z", "strength (a.u.)"
589 WRITE (log_unit, '(T10,72("-"))')
590 END IF
591
592 ! transition dipole moment
593 ALLOCATE (trans_dipoles(nstates, nderivs, nspins))
594 trans_dipoles(:, :, :) = 0.0_dp
595
596 ! nspins == 1 .AND. mult == 3 : spin-flip transitions are forbidden due to symmetry reasons
597 IF (nspins > 1 .OR. mult == 1) THEN
598 DO ispin = 1, nspins
599 CALL cp_fm_get_info(dipole_op_mos_occ(1, ispin), nrow_global=nao, ncol_global=nocc)
600 CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive)
601 IF (nocc == nactive) THEN
602 DO ideriv = 1, nderivs
603 CALL cp_fm_trace(evects(ispin, :), dipole_op_mos_occ(ideriv, ispin), &
604 trans_dipoles(:, ideriv, ispin))
605 END DO
606 ELSE ! res
607 CALL cp_fm_get_info(evects(ispin, 1), matrix_struct=matrix_struct)
608 CALL cp_fm_create(dipact, matrix_struct)
609 DO ideriv = 1, nderivs
610 DO i = 1, nactive
611 j = gs_mos(ispin)%index_active(i)
612 CALL cp_fm_to_fm(dipole_op_mos_occ(ideriv, ispin), dipact, &
613 ncol=1, source_start=j, target_start=i)
614 END DO
615 CALL cp_fm_trace(evects(ispin, :), dipact, trans_dipoles(:, ideriv, ispin))
616 END DO
617 CALL cp_fm_release(dipact)
618 END IF
619 END DO
620
621 IF (nspins == 1) THEN
622 trans_dipoles(:, :, 1) = sqrt(2.0_dp)*trans_dipoles(:, :, 1)
623 ELSE
624 trans_dipoles(:, :, 1) = sqrt(trans_dipoles(:, :, 1)**2 + trans_dipoles(:, :, 2)**2)
625 END IF
626 END IF
627
628 ! *** summary information ***
629 DO istate = 1, nstates
630 osc_strength = 2.0_dp/3.0_dp*evals(istate)* &
631 accurate_dot_product(trans_dipoles(istate, :, 1), trans_dipoles(istate, :, 1))
632 ostrength(istate) = osc_strength
633 IF (log_unit > 0) THEN
634 WRITE (log_unit, '(1X,A,T9,I7,T19,F11.5,T31,3(1X,ES11.4E2),T69,ES12.5E2)') &
635 "TDDFPT|", istate, evals(istate)*evolt, trans_dipoles(istate, 1:nderivs, 1), osc_strength
636 END IF
637 END DO
638
639 ! punch a checksum for the regs
640 IF (log_unit > 0) THEN
641 WRITE (log_unit, '(/,T2,A,E14.6)') 'TDDFPT : CheckSum =', sqrt(sum(evals**2))
642 END IF
643
644 DEALLOCATE (trans_dipoles)
645
646 CALL timestop(handle)
647 END SUBROUTINE tddfpt_print_summary
648
649! **************************************************************************************************
650!> \brief Print excitation analysis.
651!> \param log_unit output unit
652!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
653!> SIZE(evects,2) -- number of excited states to print)
654!> \param evals TDDFPT eigenvalues
655!> \param gs_mos molecular orbitals optimised for the ground state
656!> \param matrix_s overlap matrix
657!> \param spinflip ...
658!> \param min_amplitude the smallest excitation amplitude to print
659!> \par History
660!> * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov]
661!> * 08.2018 splited of from 'tddfpt_print_summary' [Sergey Chulkov]
662! **************************************************************************************************
663 SUBROUTINE tddfpt_print_excitation_analysis(log_unit, evects, evals, gs_mos, matrix_s, spinflip, &
664 min_amplitude)
665 INTEGER, INTENT(in) :: log_unit
666 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
667 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
668 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
669 INTENT(in) :: gs_mos
670 TYPE(dbcsr_type), POINTER :: matrix_s
671 INTEGER :: spinflip
672 REAL(kind=dp), INTENT(in) :: min_amplitude
673
674 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_excitation_analysis'
675
676 CHARACTER(len=5) :: spin_label, spin_label2
677 INTEGER :: handle, icol, iproc, irow, ispin, &
678 istate, nao, ncols_local, nrows_local, &
679 nspins, nstates, spin2, state_spin, &
680 state_spin2
681 INTEGER(kind=int_8) :: iexc, imo_act, imo_occ, imo_virt, ind, &
682 nexcs, nexcs_local, nexcs_max_local, &
683 nmo_virt_occ, nmo_virt_occ_alpha
684 INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:) :: inds_local, inds_recv, nexcs_recv
685 INTEGER(kind=int_8), DIMENSION(1) :: nexcs_send
686 INTEGER(kind=int_8), DIMENSION(maxspins) :: nactive8, nmo_occ8, nmo_virt8
687 INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
688 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
689 INTEGER, DIMENSION(maxspins) :: nactive, nmo_occ, nmo_virt
690 LOGICAL :: do_exc_analysis
691 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights_local, weights_neg_abs_recv, &
692 weights_recv
693 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
694 POINTER :: local_data
695 TYPE(cp_blacs_env_type), POINTER :: blacs_env
696 TYPE(cp_fm_struct_type), POINTER :: fm_struct
697 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: s_mos_virt, weights_fm
698 TYPE(mp_para_env_type), POINTER :: para_env
699 TYPE(mp_request_type) :: send_handler, send_handler2
700 TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_handlers, recv_handlers2
701
702 CALL timeset(routinen, handle)
703
704 nspins = SIZE(gs_mos, 1)
705 nstates = SIZE(evects, 2)
706 do_exc_analysis = min_amplitude < 1.0_dp
707
708 CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env, para_env=para_env)
709 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
710
711 DO ispin = 1, nspins
712 nactive(ispin) = gs_mos(ispin)%nmo_active
713 nactive8(ispin) = int(nactive(ispin), kind=int_8)
714 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
715 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
716 nmo_occ8(ispin) = SIZE(gs_mos(ispin)%evals_occ, kind=int_8)
717 nmo_virt8(ispin) = SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
718 END DO
719
720 ! *** excitation analysis ***
721 IF (do_exc_analysis) THEN
722 cpassert(log_unit <= 0 .OR. para_env%is_source())
723 nmo_virt_occ_alpha = int(nmo_virt(1), int_8)*int(nmo_occ(1), int_8)
724
725 IF (log_unit > 0) THEN
726 WRITE (log_unit, "(1X,A)") "", &
727 "-------------------------------------------------------------------------------", &
728 "- Excitation analysis -", &
729 "-------------------------------------------------------------------------------"
730 WRITE (log_unit, '(8X,A,T27,A,T49,A,T69,A)') "State", "Occupied", "Virtual", "Excitation"
731 WRITE (log_unit, '(8X,A,T28,A,T49,A,T69,A)') "number", "orbital", "orbital", "amplitude"
732 WRITE (log_unit, '(1X,79("-"))')
733
734 IF (nspins == 1) THEN
735 state_spin = 1
736 state_spin2 = 2
737 spin_label = ' '
738 spin_label2 = ' '
739 ELSE IF (spinflip /= no_sf_tddfpt) THEN
740 state_spin = 1
741 state_spin2 = 2
742 spin_label = '(alp)'
743 spin_label2 = '(bet)'
744 END IF
745 END IF
746
747 ALLOCATE (s_mos_virt(SIZE(evects, 1)), weights_fm(SIZE(evects, 1)))
748 DO ispin = 1, SIZE(evects, 1)
749 IF (spinflip == no_sf_tddfpt) THEN
750 spin2 = ispin
751 ELSE
752 spin2 = 2
753 END IF
754 CALL cp_fm_get_info(gs_mos(spin2)%mos_virt, matrix_struct=fm_struct)
755 CALL cp_fm_create(s_mos_virt(ispin), fm_struct)
756 CALL cp_dbcsr_sm_fm_multiply(matrix_s, &
757 gs_mos(spin2)%mos_virt, &
758 s_mos_virt(ispin), &
759 ncol=nmo_virt(spin2), alpha=1.0_dp, beta=0.0_dp)
760
761 NULLIFY (fm_struct)
762 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(spin2), ncol_global=nactive(ispin), &
763 context=blacs_env)
764 CALL cp_fm_create(weights_fm(ispin), fm_struct)
765 CALL cp_fm_set_all(weights_fm(ispin), 0.0_dp)
766 CALL cp_fm_struct_release(fm_struct)
767 END DO
768
769 nexcs_max_local = 0
770 DO ispin = 1, SIZE(evects, 1)
771 CALL cp_fm_get_info(weights_fm(ispin), nrow_local=nrows_local, ncol_local=ncols_local)
772 nexcs_max_local = nexcs_max_local + int(nrows_local, int_8)*int(ncols_local, int_8)
773 END DO
774
775 ALLOCATE (weights_local(nexcs_max_local), inds_local(nexcs_max_local))
776
777 DO istate = 1, nstates
778 nexcs_local = 0
779 nmo_virt_occ = 0
780
781 ! analyse matrix elements locally and transfer only significant
782 ! excitations to the master node for subsequent ordering
783 DO ispin = 1, SIZE(evects, 1)
784 IF (spinflip == no_sf_tddfpt) THEN
785 spin2 = ispin
786 ELSE
787 spin2 = 2
788 END IF
789 ! compute excitation amplitudes
790 CALL parallel_gemm('T', 'N', nmo_virt(spin2), nactive(ispin), nao, 1.0_dp, s_mos_virt(ispin), &
791 evects(ispin, istate), 0.0_dp, weights_fm(ispin))
792
793 CALL cp_fm_get_info(weights_fm(ispin), nrow_local=nrows_local, ncol_local=ncols_local, &
794 row_indices=row_indices, col_indices=col_indices, local_data=local_data)
795
796 ! locate single excitations with significant amplitudes (>= min_amplitude)
797 DO icol = 1, ncols_local
798 DO irow = 1, nrows_local
799 IF (abs(local_data(irow, icol)) >= min_amplitude) THEN
800 ! number of non-negligible excitations
801 nexcs_local = nexcs_local + 1
802 ! excitation amplitude
803 weights_local(nexcs_local) = local_data(irow, icol)
804 ! index of single excitation (ivirt, iocc, ispin) in compressed form
805 inds_local(nexcs_local) = nmo_virt_occ + int(row_indices(irow), int_8) + &
806 int(col_indices(icol) - 1, int_8)*nmo_virt8(spin2)
807 END IF
808 END DO
809 END DO
810
811 nmo_virt_occ = nmo_virt_occ + nmo_virt8(spin2)*nmo_occ8(ispin)
812 END DO
813
814 IF (para_env%is_source()) THEN
815 ! master node
816 ALLOCATE (nexcs_recv(para_env%num_pe), recv_handlers(para_env%num_pe), recv_handlers2(para_env%num_pe))
817
818 ! collect number of non-negligible excitations from other nodes
819 DO iproc = 1, para_env%num_pe
820 IF (iproc - 1 /= para_env%mepos) THEN
821 CALL para_env%irecv(nexcs_recv(iproc:iproc), iproc - 1, recv_handlers(iproc), 0)
822 ELSE
823 nexcs_recv(iproc) = nexcs_local
824 END IF
825 END DO
826
827 DO iproc = 1, para_env%num_pe
828 IF (iproc - 1 /= para_env%mepos) &
829 CALL recv_handlers(iproc)%wait()
830 END DO
831
832 ! compute total number of non-negligible excitations
833 nexcs = 0
834 DO iproc = 1, para_env%num_pe
835 nexcs = nexcs + nexcs_recv(iproc)
836 END DO
837
838 ! receive indices and amplitudes of selected excitations
839 ALLOCATE (weights_recv(nexcs), weights_neg_abs_recv(nexcs))
840 ALLOCATE (inds_recv(nexcs), inds(nexcs))
841
842 nmo_virt_occ = 0
843 DO iproc = 1, para_env%num_pe
844 IF (nexcs_recv(iproc) > 0) THEN
845 IF (iproc - 1 /= para_env%mepos) THEN
846 ! excitation amplitudes
847 CALL para_env%irecv(weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), &
848 iproc - 1, recv_handlers(iproc), 1)
849 ! compressed indices
850 CALL para_env%irecv(inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), &
851 iproc - 1, recv_handlers2(iproc), 2)
852 ELSE
853 ! data on master node
854 weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = weights_local(1:nexcs_recv(iproc))
855 inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = inds_local(1:nexcs_recv(iproc))
856 END IF
857
858 nmo_virt_occ = nmo_virt_occ + nexcs_recv(iproc)
859 END IF
860 END DO
861
862 DO iproc = 1, para_env%num_pe
863 IF (iproc - 1 /= para_env%mepos .AND. nexcs_recv(iproc) > 0) THEN
864 CALL recv_handlers(iproc)%wait()
865 CALL recv_handlers2(iproc)%wait()
866 END IF
867 END DO
868
869 DEALLOCATE (nexcs_recv, recv_handlers, recv_handlers2)
870 ELSE
871 ! working node: send the number of selected excited states to the master node
872 nexcs_send(1) = nexcs_local
873 CALL para_env%isend(nexcs_send, para_env%source, send_handler, 0)
874 CALL send_handler%wait()
875
876 IF (nexcs_local > 0) THEN
877 ! send excitation amplitudes
878 CALL para_env%isend(weights_local(1:nexcs_local), para_env%source, send_handler, 1)
879 ! send compressed indices
880 CALL para_env%isend(inds_local(1:nexcs_local), para_env%source, send_handler2, 2)
881
882 CALL send_handler%wait()
883 CALL send_handler2%wait()
884 END IF
885 END IF
886
887 ! sort non-negligible excitations on the master node according to their amplitudes,
888 ! uncompress indices and print summary information
889 IF (para_env%is_source() .AND. log_unit > 0) THEN
890 weights_neg_abs_recv(:) = -abs(weights_recv)
891 CALL sort(weights_neg_abs_recv, int(nexcs), inds)
892
893 WRITE (log_unit, '(T7,I8,F10.5,A)') istate, evals(istate)*evolt, " eV"
894
895 ! This reinitialization is needed to prevent the intel fortran compiler from introduce
896 ! a bug when using optimization level 3 flag
897 state_spin = 1
898 state_spin2 = 1
899 IF (spinflip /= no_sf_tddfpt) THEN
900 state_spin = 1
901 state_spin2 = 2
902 END IF
903 DO iexc = 1, nexcs
904 ind = inds_recv(inds(iexc)) - 1
905 IF ((nspins > 1) .AND. (spinflip == no_sf_tddfpt)) THEN
906 IF (ind < nmo_virt_occ_alpha) THEN
907 spin_label = '(alp)'
908 spin_label2 = '(alp)'
909 ELSE
910 state_spin = 2
911 state_spin2 = 2
912 ind = ind - nmo_virt_occ_alpha
913 spin_label = '(bet)'
914 spin_label2 = '(bet)'
915 END IF
916 END IF
917 imo_act = ind/nmo_virt8(state_spin2) + 1
918 imo_occ = gs_mos(state_spin)%index_active(imo_act)
919 imo_virt = mod(ind, nmo_virt8(state_spin2)) + 1
920
921 WRITE (log_unit, '(T27,I8,1X,A5,T48,I8,1X,A5,T70,F9.6)') imo_occ, spin_label, &
922 nmo_occ8(state_spin2) + imo_virt, spin_label2, weights_recv(inds(iexc))
923 END DO
924 END IF
925
926 ! deallocate temporary arrays
927 IF (para_env%is_source()) &
928 DEALLOCATE (weights_recv, weights_neg_abs_recv, inds_recv, inds)
929 END DO
930
931 DEALLOCATE (weights_local, inds_local)
932 IF (log_unit > 0) THEN
933 WRITE (log_unit, "(1X,A)") &
934 "-------------------------------------------------------------------------------"
935 END IF
936 END IF
937
938 CALL cp_fm_release(weights_fm)
939 CALL cp_fm_release(s_mos_virt)
940
941 CALL timestop(handle)
942
944
945! **************************************************************************************************
946!> \brief Print natural transition orbital analysis.
947!> \param qs_env Information on Kinds and Particles
948!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
949!> SIZE(evects,2) -- number of excited states to print)
950!> \param evals TDDFPT eigenvalues
951!> \param ostrength ...
952!> \param gs_mos molecular orbitals optimised for the ground state
953!> \param matrix_s overlap matrix
954!> \param print_section ...
955!> \par History
956!> * 06.2019 created [JGH]
957! **************************************************************************************************
958 SUBROUTINE tddfpt_print_nto_analysis(qs_env, evects, evals, ostrength, gs_mos, matrix_s, print_section)
959 TYPE(qs_environment_type), POINTER :: qs_env
960 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
961 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals, ostrength
962 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
963 INTENT(in) :: gs_mos
964 TYPE(dbcsr_type), POINTER :: matrix_s
965 TYPE(section_vals_type), POINTER :: print_section
966
967 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_nto_analysis'
968 INTEGER, PARAMETER :: ntomax = 10
969
970 CHARACTER(LEN=20), DIMENSION(2) :: nto_name
971 INTEGER :: handle, i, ia, icg, iounit, ispin, &
972 istate, j, nao, nlist, nmax, nmo, &
973 nnto, nspins, nstates
974 INTEGER, DIMENSION(2) :: iv
975 INTEGER, DIMENSION(2, ntomax) :: ia_index
976 INTEGER, DIMENSION(:), POINTER :: slist, stride
977 LOGICAL :: append_cube, cube_file, explicit
978 REAL(kind=dp) :: os_threshold, sume, threshold
979 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvals
980 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvalues
981 REAL(kind=dp), DIMENSION(ntomax) :: ia_eval
982 TYPE(cell_type), POINTER :: cell
983 TYPE(cp_fm_struct_type), POINTER :: fm_mo_struct, fm_struct
984 TYPE(cp_fm_type) :: sev, smat, tmat, wmat, work, wvec
985 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: teig
986 TYPE(cp_logger_type), POINTER :: logger
987 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: nto_set
988 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
989 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
990 TYPE(section_vals_type), POINTER :: molden_section, nto_section
991
992 CALL timeset(routinen, handle)
993
994 logger => cp_get_default_logger()
995 iounit = cp_logger_get_default_io_unit(logger)
996
997 IF (btest(cp_print_key_should_output(logger%iter_info, print_section, &
998 "NTO_ANALYSIS"), cp_p_file)) THEN
999
1000 CALL cite_reference(martin2003)
1001
1002 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%THRESHOLD", r_val=threshold)
1003 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%INTENSITY_THRESHOLD", r_val=os_threshold)
1004 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", explicit=explicit)
1005
1006 IF (explicit) THEN
1007 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", i_vals=slist)
1008 nlist = SIZE(slist)
1009 ELSE
1010 nlist = 0
1011 END IF
1012
1013 IF (iounit > 0) THEN
1014 WRITE (iounit, "(1X,A)") "", &
1015 "-------------------------------------------------------------------------------", &
1016 "- Natural Orbital analysis -", &
1017 "-------------------------------------------------------------------------------"
1018 END IF
1019
1020 nspins = SIZE(evects, 1)
1021 nstates = SIZE(evects, 2)
1022 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1023
1024 DO istate = 1, nstates
1025 IF (os_threshold > ostrength(istate)) THEN
1026 IF (iounit > 0) THEN
1027 WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1028 END IF
1029 cycle
1030 END IF
1031 IF (nlist > 0) THEN
1032 IF (.NOT. any(slist == istate)) THEN
1033 IF (iounit > 0) THEN
1034 WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1035 END IF
1036 cycle
1037 END IF
1038 END IF
1039 IF (iounit > 0) THEN
1040 WRITE (iounit, "(1X,A,I6,T30,F10.5,A)") " STATE NR. ", istate, evals(istate)*evolt, " eV"
1041 END IF
1042 nmax = 0
1043 DO ispin = 1, nspins
1044 CALL cp_fm_get_info(evects(ispin, istate), matrix_struct=fm_struct, ncol_global=nmo)
1045 nmax = max(nmax, nmo)
1046 END DO
1047 ALLOCATE (eigenvalues(nmax, nspins))
1048 eigenvalues = 0.0_dp
1049 ! SET 1: Hole states
1050 ! SET 2: Particle states
1051 nto_name(1) = 'Hole_states'
1052 nto_name(2) = 'Particle_states'
1053 ALLOCATE (nto_set(2))
1054 DO i = 1, 2
1055 CALL allocate_mo_set(nto_set(i), nao, ntomax, 0, 0.0_dp, 1.0_dp, 0.0_dp)
1056 CALL cp_fm_get_info(evects(1, istate), matrix_struct=fm_struct)
1057 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1058 ncol_global=ntomax)
1059 CALL cp_fm_create(tmat, fm_mo_struct)
1060 CALL init_mo_set(nto_set(i), fm_ref=tmat, name=nto_name(i))
1061 CALL cp_fm_release(tmat)
1062 CALL cp_fm_struct_release(fm_mo_struct)
1063 END DO
1064 !
1065 ALLOCATE (teig(nspins))
1066 ! hole states
1067 ! Diagonalize X(T)*S*X
1068 DO ispin = 1, nspins
1069 associate(ev => evects(ispin, istate))
1070 CALL cp_fm_get_info(ev, matrix_struct=fm_struct, ncol_global=nmo)
1071 CALL cp_fm_create(sev, fm_struct)
1072 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1073 nrow_global=nmo, ncol_global=nmo)
1074 CALL cp_fm_create(tmat, fm_mo_struct)
1075 CALL cp_fm_create(teig(ispin), fm_mo_struct)
1076 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1077 CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, ev, sev, 0.0_dp, tmat)
1078 END associate
1079
1080 CALL choose_eigv_solver(tmat, teig(ispin), eigenvalues(1:nmo, ispin))
1081
1082 CALL cp_fm_struct_release(fm_mo_struct)
1083 CALL cp_fm_release(tmat)
1084 CALL cp_fm_release(sev)
1085 END DO
1086 ! find major determinants i->a
1087 ia_index = 0
1088 sume = 0.0_dp
1089 nnto = 0
1090 DO i = 1, ntomax
1091 iv = maxloc(eigenvalues)
1092 ia_eval(i) = eigenvalues(iv(1), iv(2))
1093 ia_index(1:2, i) = iv(1:2)
1094 sume = sume + ia_eval(i)
1095 eigenvalues(iv(1), iv(2)) = 0.0_dp
1096 nnto = nnto + 1
1097 IF (sume > threshold) EXIT
1098 END DO
1099 ! store hole states
1100 CALL set_mo_set(nto_set(1), nmo=nnto)
1101 DO i = 1, nnto
1102 ia = ia_index(1, i)
1103 ispin = ia_index(2, i)
1104 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=nmo)
1105 CALL cp_fm_get_info(teig(ispin), matrix_struct=fm_struct)
1106 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1107 nrow_global=nmo, ncol_global=1)
1108 CALL cp_fm_create(tmat, fm_mo_struct)
1109 CALL cp_fm_struct_release(fm_mo_struct)
1110 CALL cp_fm_get_info(gs_mos(1)%mos_occ, matrix_struct=fm_struct)
1111 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1112 ncol_global=1)
1113 CALL cp_fm_create(wvec, fm_mo_struct)
1114 CALL cp_fm_struct_release(fm_mo_struct)
1115 CALL cp_fm_to_fm(teig(ispin), tmat, 1, ia, 1)
1116 CALL parallel_gemm('N', 'N', nao, 1, nmo, 1.0_dp, gs_mos(ispin)%mos_occ, &
1117 tmat, 0.0_dp, wvec)
1118 CALL cp_fm_to_fm(wvec, nto_set(1)%mo_coeff, 1, 1, i)
1119 CALL cp_fm_release(wvec)
1120 CALL cp_fm_release(tmat)
1121 END DO
1122 ! particle states
1123 ! Solve generalized eigenvalue equation: (S*X)*(S*X)(T)*v = lambda*S*v
1124 CALL set_mo_set(nto_set(2), nmo=nnto)
1125 DO ispin = 1, nspins
1126 associate(ev => evects(ispin, istate))
1127 CALL cp_fm_get_info(ev, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
1128 ALLOCATE (eigvals(nao))
1129 eigvals = 0.0_dp
1130 CALL cp_fm_create(sev, fm_struct)
1131 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1132 END associate
1133 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1134 nrow_global=nao, ncol_global=nao)
1135 CALL cp_fm_create(tmat, fm_mo_struct)
1136 CALL cp_fm_create(smat, fm_mo_struct)
1137 CALL cp_fm_create(wmat, fm_mo_struct)
1138 CALL cp_fm_create(work, fm_mo_struct)
1139 CALL cp_fm_struct_release(fm_mo_struct)
1140 CALL copy_dbcsr_to_fm(matrix_s, smat)
1141 CALL parallel_gemm('N', 'T', nao, nao, nmo, 1.0_dp, sev, sev, 0.0_dp, tmat)
1142 CALL cp_fm_geeig(tmat, smat, wmat, eigvals, work)
1143 DO i = 1, nnto
1144 IF (ispin == ia_index(2, i)) THEN
1145 icg = 0
1146 DO j = 1, nao
1147 IF (abs(eigvals(j) - ia_eval(i)) < 1.e-6_dp) THEN
1148 icg = j
1149 EXIT
1150 END IF
1151 END DO
1152 IF (icg == 0) THEN
1153 CALL cp_warn(__location__, &
1154 "Could not locate particle state associated with hole state.")
1155 ELSE
1156 CALL cp_fm_to_fm(wmat, nto_set(2)%mo_coeff, 1, icg, i)
1157 END IF
1158 END IF
1159 END DO
1160 DEALLOCATE (eigvals)
1161 CALL cp_fm_release(sev)
1162 CALL cp_fm_release(tmat)
1163 CALL cp_fm_release(smat)
1164 CALL cp_fm_release(wmat)
1165 CALL cp_fm_release(work)
1166 END DO
1167 ! print
1168 IF (iounit > 0) THEN
1169 sume = 0.0_dp
1170 DO i = 1, nnto
1171 sume = sume + ia_eval(i)
1172 WRITE (iounit, "(T6,A,i2,T30,A,i1,T42,A,F8.5,T63,A,F8.5)") &
1173 "Particle-Hole state:", i, " Spin:", ia_index(2, i), &
1174 "Eigenvalue:", ia_eval(i), " Sum Eigv:", sume
1175 END DO
1176 END IF
1177 ! Cube and Molden files
1178 nto_section => section_vals_get_subs_vals(print_section, "NTO_ANALYSIS")
1179 CALL section_vals_val_get(nto_section, "CUBE_FILES", l_val=cube_file)
1180 CALL section_vals_val_get(nto_section, "STRIDE", i_vals=stride)
1181 CALL section_vals_val_get(nto_section, "APPEND", l_val=append_cube)
1182 IF (cube_file) THEN
1183 CALL print_nto_cubes(qs_env, nto_set, istate, stride, append_cube, nto_section)
1184 END IF
1185 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, cell=cell)
1186 molden_section => section_vals_get_subs_vals(print_section, "MOS_MOLDEN")
1187 CALL write_mos_molden(nto_set, qs_kind_set, particle_set, molden_section, cell=cell)
1188 !
1189 DEALLOCATE (eigenvalues)
1190 CALL cp_fm_release(teig)
1191 !
1192 DO i = 1, 2
1193 CALL deallocate_mo_set(nto_set(i))
1194 END DO
1195 DEALLOCATE (nto_set)
1196 END DO
1197
1198 IF (iounit > 0) THEN
1199 WRITE (iounit, "(1X,A)") &
1200 "-------------------------------------------------------------------------------"
1201 END IF
1202
1203 END IF
1204
1205 CALL timestop(handle)
1206
1207 END SUBROUTINE tddfpt_print_nto_analysis
1208
1209! **************************************************************************************************
1210!> \brief Print exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018)
1211!> \param log_unit output unit
1212!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
1213!> SIZE(evects,2) -- number of excited states to print)
1214!> \param gs_mos molecular orbitals optimised for the ground state
1215!> \param matrix_s overlap matrix
1216!> \param do_directional_exciton_descriptors flag for computing descriptors for each (cartesian) direction
1217!> \param qs_env Information on particles/geometry
1218!> \par History
1219!> * 12.2024 created as 'tddfpt_print_exciton_descriptors' [Maximilian Graml]
1220! **************************************************************************************************
1221 SUBROUTINE tddfpt_print_exciton_descriptors(log_unit, evects, gs_mos, matrix_s, &
1222 do_directional_exciton_descriptors, qs_env)
1223 INTEGER, INTENT(in) :: log_unit
1224 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
1225 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1226 INTENT(in) :: gs_mos
1227 TYPE(dbcsr_type), POINTER :: matrix_s
1228 LOGICAL, INTENT(IN) :: do_directional_exciton_descriptors
1229 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1230
1231 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_exciton_descriptors'
1232
1233 CHARACTER(LEN=4) :: prefix_output
1234 INTEGER :: handle, ispin, istate, n_moments_quad, &
1235 nactive, nao, nspins, nstates
1236 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt
1237 LOGICAL :: print_checkvalue
1238 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ref_point_multipole
1239 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1240 TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo_coeff, &
1241 fm_struct_s_mos_virt, fm_struct_x_ia_n
1242 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: eigvec_x_ia_n, fm_multipole_ab, &
1243 fm_multipole_ai, fm_multipole_ij, &
1244 s_mos_virt
1245 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_coeff
1246 TYPE(exciton_descr_type), ALLOCATABLE, &
1247 DIMENSION(:) :: exc_descr
1248
1249 CALL timeset(routinen, handle)
1250
1251 nspins = SIZE(evects, 1)
1252 nstates = SIZE(evects, 2)
1253
1254 cpassert(nspins == 1) ! Other spins are not yet implemented for exciton descriptors
1255
1256 CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env)
1257 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1258
1259 DO ispin = 1, nspins
1260 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
1261 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
1262 END DO
1263
1264 ! Prepare fm with all MO coefficents, i.e. nao x nao
1265 ALLOCATE (mo_coeff(nspins))
1266 CALL cp_fm_struct_create(fm_struct_mo_coeff, nrow_global=nao, ncol_global=nao, &
1267 context=blacs_env)
1268 DO ispin = 1, nspins
1269 CALL cp_fm_create(mo_coeff(ispin), fm_struct_mo_coeff)
1270 CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_occ, &
1271 mo_coeff(ispin), &
1272 nao, &
1273 nmo_occ(ispin), &
1274 1, &
1275 1, &
1276 1, &
1277 1, &
1278 blacs_env)
1279 CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_virt, &
1280 mo_coeff(ispin), &
1281 nao, &
1282 nmo_virt(ispin), &
1283 1, &
1284 1, &
1285 1, &
1286 nmo_occ(ispin) + 1, &
1287 blacs_env)
1288 END DO
1289 CALL cp_fm_struct_release(fm_struct_mo_coeff)
1290
1291 ! Compute multipole moments
1292 ! fm_multipole_XY have structure inherited by libint, i.e. x, y, z, xx, xy, xz, yy, yz, zz
1293 n_moments_quad = 9
1294 ALLOCATE (ref_point_multipole(3))
1295 ALLOCATE (fm_multipole_ij(n_moments_quad))
1296 ALLOCATE (fm_multipole_ab(n_moments_quad))
1297 ALLOCATE (fm_multipole_ai(n_moments_quad))
1298
1299 CALL get_multipoles_mo(fm_multipole_ai, fm_multipole_ij, fm_multipole_ab, &
1300 qs_env, mo_coeff, ref_point_multipole, 2, &
1301 nmo_occ(1), nmo_virt(1), blacs_env)
1302
1303 CALL cp_fm_release(mo_coeff)
1304
1305 ! Compute eigenvector X of the Casida equation from trial vectors
1306 ALLOCATE (s_mos_virt(nspins), eigvec_x_ia_n(nspins))
1307 DO ispin = 1, nspins
1308 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct_s_mos_virt)
1309 CALL cp_fm_create(s_mos_virt(ispin), fm_struct_s_mos_virt)
1310 NULLIFY (fm_struct_s_mos_virt)
1311 CALL cp_dbcsr_sm_fm_multiply(matrix_s, &
1312 gs_mos(ispin)%mos_virt, &
1313 s_mos_virt(ispin), &
1314 ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
1315
1316 CALL cp_fm_struct_create(fm_struct_x_ia_n, nrow_global=nmo_occ(ispin), ncol_global=nmo_virt(ispin), &
1317 context=blacs_env)
1318 CALL cp_fm_create(eigvec_x_ia_n(ispin), fm_struct_x_ia_n)
1319 CALL cp_fm_struct_release(fm_struct_x_ia_n)
1320 END DO
1321 ALLOCATE (exc_descr(nstates))
1322 DO istate = 1, nstates
1323 DO ispin = 1, nspins
1324 CALL cp_fm_set_all(eigvec_x_ia_n(ispin), 0.0_dp)
1325 ! compute eigenvectors X of the TDA equation
1326 ! Reshuffle multiplication from
1327 ! X_ai = S_ma ^T * C_mi
1328 ! to
1329 ! X_ia = C_mi ^T * S_ma
1330 ! for compatibility with the structure needed for get_exciton_descriptors of bse_properties.F
1331 CALL cp_fm_get_info(evects(ispin, istate), ncol_global=nactive)
1332 IF (nactive /= nmo_occ(ispin)) THEN
1333 CALL cp_abort(__location__, &
1334 "Reduced active space excitations not implemented")
1335 END IF
1336 CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_virt(ispin), nao, 1.0_dp, &
1337 evects(ispin, istate), s_mos_virt(ispin), 0.0_dp, eigvec_x_ia_n(ispin))
1338
1339 CALL get_exciton_descriptors(exc_descr, eigvec_x_ia_n(ispin), &
1340 fm_multipole_ij, fm_multipole_ab, &
1341 fm_multipole_ai, &
1342 istate, nmo_occ(ispin), nmo_virt(ispin))
1343 END DO
1344 END DO
1345 CALL cp_fm_release(eigvec_x_ia_n)
1346 CALL cp_fm_release(s_mos_virt)
1347 CALL cp_fm_release(fm_multipole_ai)
1348 CALL cp_fm_release(fm_multipole_ij)
1349 CALL cp_fm_release(fm_multipole_ab)
1350
1351 ! Actual printing
1352 print_checkvalue = .true.
1353 prefix_output = ' '
1354 CALL print_exciton_descriptors(exc_descr, ref_point_multipole, log_unit, &
1355 nstates, print_checkvalue, do_directional_exciton_descriptors, &
1356 prefix_output, qs_env)
1357
1358 DEALLOCATE (ref_point_multipole)
1359 DEALLOCATE (exc_descr)
1360
1361 CALL timestop(handle)
1362
1364
1365! **************************************************************************************************
1366!> \brief ...
1367!> \param vin ...
1368!> \param vout ...
1369!> \param mos_occ ...
1370!> \param matrix_s ...
1371! **************************************************************************************************
1372 SUBROUTINE project_vector(vin, vout, mos_occ, matrix_s)
1373 TYPE(dbcsr_type) :: vin, vout
1374 TYPE(cp_fm_type), INTENT(IN) :: mos_occ
1375 TYPE(dbcsr_type), POINTER :: matrix_s
1376
1377 CHARACTER(LEN=*), PARAMETER :: routinen = 'project_vector'
1378
1379 INTEGER :: handle, nao, nmo
1380 REAL(kind=dp) :: norm(1)
1381 TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_vec_struct
1382 TYPE(cp_fm_type) :: csvec, svec, vec
1383
1384 CALL timeset(routinen, handle)
1385
1386 CALL cp_fm_get_info(mos_occ, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
1387 CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
1388 nrow_global=nao, ncol_global=1)
1389 CALL cp_fm_create(vec, fm_vec_struct)
1390 CALL cp_fm_create(svec, fm_vec_struct)
1391 CALL cp_fm_struct_release(fm_vec_struct)
1392 CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
1393 nrow_global=nmo, ncol_global=1)
1394 CALL cp_fm_create(csvec, fm_vec_struct)
1395 CALL cp_fm_struct_release(fm_vec_struct)
1396
1397 CALL copy_dbcsr_to_fm(vin, vec)
1398 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vec, svec, ncol=1, alpha=1.0_dp, beta=0.0_dp)
1399 CALL parallel_gemm('T', 'N', nmo, 1, nao, 1.0_dp, mos_occ, svec, 0.0_dp, csvec)
1400 CALL parallel_gemm('N', 'N', nao, 1, nmo, -1.0_dp, mos_occ, csvec, 1.0_dp, vec)
1401 CALL cp_fm_vectorsnorm(vec, norm)
1402 cpassert(norm(1) > 1.e-14_dp)
1403 norm(1) = sqrt(1._dp/norm(1))
1404 CALL cp_fm_scale(norm(1), vec)
1405 CALL copy_fm_to_dbcsr(vec, vout, keep_sparsity=.false.)
1406
1407 CALL cp_fm_release(csvec)
1408 CALL cp_fm_release(svec)
1409 CALL cp_fm_release(vec)
1410
1411 CALL timestop(handle)
1412
1413 END SUBROUTINE project_vector
1414
1415! **************************************************************************************************
1416!> \brief ...
1417!> \param va ...
1418!> \param vb ...
1419!> \param res ...
1420! **************************************************************************************************
1421 SUBROUTINE vec_product(va, vb, res)
1422 TYPE(dbcsr_type) :: va, vb
1423 REAL(kind=dp), INTENT(OUT) :: res
1424
1425 CHARACTER(LEN=*), PARAMETER :: routinen = 'vec_product'
1426
1427 INTEGER :: handle, icol, irow
1428 LOGICAL :: found
1429 REAL(kind=dp), DIMENSION(:, :), POINTER :: vba, vbb
1430 TYPE(dbcsr_iterator_type) :: iter
1431 TYPE(mp_comm_type) :: group
1432
1433 CALL timeset(routinen, handle)
1434
1435 res = 0.0_dp
1436
1437 CALL dbcsr_get_info(va, group=group)
1438 CALL dbcsr_iterator_start(iter, va)
1439 DO WHILE (dbcsr_iterator_blocks_left(iter))
1440 CALL dbcsr_iterator_next_block(iter, irow, icol, vba)
1441 CALL dbcsr_get_block_p(vb, row=irow, col=icol, block=vbb, found=found)
1442 res = res + sum(vba*vbb)
1443 cpassert(found)
1444 END DO
1445 CALL dbcsr_iterator_stop(iter)
1446 CALL group%sum(res)
1447
1448 CALL timestop(handle)
1449
1450 END SUBROUTINE vec_product
1451
1452! **************************************************************************************************
1453!> \brief ...
1454!> \param qs_env ...
1455!> \param mos ...
1456!> \param istate ...
1457!> \param stride ...
1458!> \param append_cube ...
1459!> \param print_section ...
1460! **************************************************************************************************
1461 SUBROUTINE print_nto_cubes(qs_env, mos, istate, stride, append_cube, print_section)
1462
1463 TYPE(qs_environment_type), POINTER :: qs_env
1464 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1465 INTEGER, INTENT(IN) :: istate
1466 INTEGER, DIMENSION(:), POINTER :: stride
1467 LOGICAL, INTENT(IN) :: append_cube
1468 TYPE(section_vals_type), POINTER :: print_section
1469
1470 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
1471 INTEGER :: i, iset, nmo, unit_nr
1472 LOGICAL :: mpi_io
1473 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1474 TYPE(cell_type), POINTER :: cell
1475 TYPE(cp_fm_type), POINTER :: mo_coeff
1476 TYPE(cp_logger_type), POINTER :: logger
1477 TYPE(dft_control_type), POINTER :: dft_control
1478 TYPE(particle_list_type), POINTER :: particles
1479 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1480 TYPE(pw_c1d_gs_type) :: wf_g
1481 TYPE(pw_env_type), POINTER :: pw_env
1482 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1483 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1484 TYPE(pw_r3d_rs_type) :: wf_r
1485 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1486 TYPE(qs_subsys_type), POINTER :: subsys
1487
1488 logger => cp_get_default_logger()
1489
1490 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
1491 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1492 CALL auxbas_pw_pool%create_pw(wf_r)
1493 CALL auxbas_pw_pool%create_pw(wf_g)
1494
1495 CALL get_qs_env(qs_env, subsys=subsys)
1496 CALL qs_subsys_get(subsys, particles=particles)
1497
1498 my_pos_cube = "REWIND"
1499 IF (append_cube) THEN
1500 my_pos_cube = "APPEND"
1501 END IF
1502
1503 CALL get_qs_env(qs_env=qs_env, &
1504 atomic_kind_set=atomic_kind_set, &
1505 qs_kind_set=qs_kind_set, &
1506 cell=cell, &
1507 particle_set=particle_set)
1508
1509 DO iset = 1, 2
1510 CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo)
1511 DO i = 1, nmo
1512 CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1513 cell, dft_control, particle_set, pw_env)
1514 IF (iset == 1) THEN
1515 WRITE (filename, '(a4,I3.3,I2.2,a11)') "NTO_STATE", istate, i, "_Hole_State"
1516 ELSEIF (iset == 2) THEN
1517 WRITE (filename, '(a4,I3.3,I2.2,a15)') "NTO_STATE", istate, i, "_Particle_State"
1518 END IF
1519 mpi_io = .true.
1520 unit_nr = cp_print_key_unit_nr(logger, print_section, '', extension=".cube", &
1521 middle_name=trim(filename), file_position=my_pos_cube, &
1522 log_filename=.false., ignore_should_output=.true., mpi_io=mpi_io)
1523 IF (iset == 1) THEN
1524 WRITE (title, *) "Natural Transition Orbital Hole State", i
1525 ELSEIF (iset == 2) THEN
1526 WRITE (title, *) "Natural Transition Orbital Particle State", i
1527 END IF
1528 CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, stride=stride, mpi_io=mpi_io)
1529 CALL cp_print_key_finished_output(unit_nr, logger, print_section, '', &
1530 ignore_should_output=.true., mpi_io=mpi_io)
1531 END DO
1532 END DO
1533
1534 CALL auxbas_pw_pool%give_back_pw(wf_g)
1535 CALL auxbas_pw_pool%give_back_pw(wf_r)
1536
1537 END SUBROUTINE print_nto_cubes
1538
1539END 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...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
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_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 Single option version: Cholesky decomposition of B.
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:239
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_velocity
integer, parameter, public tddfpt_dipole_length
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
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
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)
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 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.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
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.