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