(git:03cdb6f)
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 state_spin = 1
908 state_spin2 = 1
909 spin_label = '(alp)'
910 spin_label2 = '(alp)'
911 ELSE
912 state_spin = 2
913 state_spin2 = 2
914 ind = ind - nmo_virt_occ_alpha
915 spin_label = '(bet)'
916 spin_label2 = '(bet)'
917 END IF
918 END IF
919 imo_act = ind/nmo_virt8(state_spin2) + 1
920 imo_occ = gs_mos(state_spin)%index_active(imo_act)
921 imo_virt = mod(ind, nmo_virt8(state_spin2)) + 1
922
923 WRITE (log_unit, '(T27,I8,1X,A5,T48,I8,1X,A5,T70,F9.6)') imo_occ, spin_label, &
924 nmo_occ8(state_spin2) + imo_virt, spin_label2, weights_recv(inds(iexc))
925 END DO
926 END IF
927
928 ! deallocate temporary arrays
929 IF (para_env%is_source()) &
930 DEALLOCATE (weights_recv, weights_neg_abs_recv, inds_recv, inds)
931 END DO
932
933 DEALLOCATE (weights_local, inds_local)
934 IF (log_unit > 0) THEN
935 WRITE (log_unit, "(1X,A)") &
936 "-------------------------------------------------------------------------------"
937 END IF
938 END IF
939
940 CALL cp_fm_release(weights_fm)
941 CALL cp_fm_release(s_mos_virt)
942
943 CALL timestop(handle)
944
946
947! **************************************************************************************************
948!> \brief Print natural transition orbital analysis.
949!> \param qs_env Information on Kinds and Particles
950!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
951!> SIZE(evects,2) -- number of excited states to print)
952!> \param evals TDDFPT eigenvalues
953!> \param ostrength ...
954!> \param gs_mos molecular orbitals optimised for the ground state
955!> \param matrix_s overlap matrix
956!> \param print_section ...
957!> \par History
958!> * 06.2019 created [JGH]
959! **************************************************************************************************
960 SUBROUTINE tddfpt_print_nto_analysis(qs_env, evects, evals, ostrength, gs_mos, matrix_s, print_section)
961 TYPE(qs_environment_type), POINTER :: qs_env
962 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
963 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals, ostrength
964 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
965 INTENT(in) :: gs_mos
966 TYPE(dbcsr_type), POINTER :: matrix_s
967 TYPE(section_vals_type), POINTER :: print_section
968
969 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_nto_analysis'
970 INTEGER, PARAMETER :: ntomax = 10
971
972 CHARACTER(LEN=20), DIMENSION(2) :: nto_name
973 INTEGER :: handle, i, ia, icg, iounit, ispin, &
974 istate, j, nao, nlist, nmax, nmo, &
975 nnto, nspins, nstates
976 INTEGER, DIMENSION(2) :: iv
977 INTEGER, DIMENSION(2, ntomax) :: ia_index
978 INTEGER, DIMENSION(:), POINTER :: slist, stride
979 LOGICAL :: append_cube, cube_file, explicit
980 REAL(kind=dp) :: os_threshold, sume, threshold
981 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvals
982 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvalues
983 REAL(kind=dp), DIMENSION(ntomax) :: ia_eval
984 TYPE(cell_type), POINTER :: cell
985 TYPE(cp_fm_struct_type), POINTER :: fm_mo_struct, fm_struct
986 TYPE(cp_fm_type) :: sev, smat, tmat, wmat, work, wvec
987 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: teig
988 TYPE(cp_logger_type), POINTER :: logger
989 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: nto_set
990 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
991 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
992 TYPE(section_vals_type), POINTER :: molden_section, nto_section
993
994 CALL timeset(routinen, handle)
995
996 logger => cp_get_default_logger()
997 iounit = cp_logger_get_default_io_unit(logger)
998
999 IF (btest(cp_print_key_should_output(logger%iter_info, print_section, &
1000 "NTO_ANALYSIS"), cp_p_file)) THEN
1001
1002 CALL cite_reference(martin2003)
1003
1004 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%THRESHOLD", r_val=threshold)
1005 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%INTENSITY_THRESHOLD", r_val=os_threshold)
1006 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", explicit=explicit)
1007
1008 IF (explicit) THEN
1009 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", i_vals=slist)
1010 nlist = SIZE(slist)
1011 ELSE
1012 nlist = 0
1013 END IF
1014
1015 IF (iounit > 0) THEN
1016 WRITE (iounit, "(1X,A)") "", &
1017 "-------------------------------------------------------------------------------", &
1018 "- Natural Orbital analysis -", &
1019 "-------------------------------------------------------------------------------"
1020 END IF
1021
1022 nspins = SIZE(evects, 1)
1023 nstates = SIZE(evects, 2)
1024 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1025
1026 DO istate = 1, nstates
1027 IF (os_threshold > ostrength(istate)) THEN
1028 IF (iounit > 0) THEN
1029 WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1030 END IF
1031 cycle
1032 END IF
1033 IF (nlist > 0) THEN
1034 IF (.NOT. any(slist == istate)) THEN
1035 IF (iounit > 0) THEN
1036 WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1037 END IF
1038 cycle
1039 END IF
1040 END IF
1041 IF (iounit > 0) THEN
1042 WRITE (iounit, "(1X,A,I6,T30,F10.5,A)") " STATE NR. ", istate, evals(istate)*evolt, " eV"
1043 END IF
1044 nmax = 0
1045 DO ispin = 1, nspins
1046 CALL cp_fm_get_info(evects(ispin, istate), matrix_struct=fm_struct, ncol_global=nmo)
1047 nmax = max(nmax, nmo)
1048 END DO
1049 ALLOCATE (eigenvalues(nmax, nspins))
1050 eigenvalues = 0.0_dp
1051 ! SET 1: Hole states
1052 ! SET 2: Particle states
1053 nto_name(1) = 'Hole_states'
1054 nto_name(2) = 'Particle_states'
1055 ALLOCATE (nto_set(2))
1056 DO i = 1, 2
1057 CALL allocate_mo_set(nto_set(i), nao, ntomax, 0, 0.0_dp, 1.0_dp, 0.0_dp)
1058 CALL cp_fm_get_info(evects(1, istate), matrix_struct=fm_struct)
1059 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1060 ncol_global=ntomax)
1061 CALL cp_fm_create(tmat, fm_mo_struct)
1062 CALL init_mo_set(nto_set(i), fm_ref=tmat, name=nto_name(i))
1063 CALL cp_fm_release(tmat)
1064 CALL cp_fm_struct_release(fm_mo_struct)
1065 END DO
1066 !
1067 ALLOCATE (teig(nspins))
1068 ! hole states
1069 ! Diagonalize X(T)*S*X
1070 DO ispin = 1, nspins
1071 associate(ev => evects(ispin, istate))
1072 CALL cp_fm_get_info(ev, matrix_struct=fm_struct, ncol_global=nmo)
1073 CALL cp_fm_create(sev, fm_struct)
1074 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1075 nrow_global=nmo, ncol_global=nmo)
1076 CALL cp_fm_create(tmat, fm_mo_struct)
1077 CALL cp_fm_create(teig(ispin), fm_mo_struct)
1078 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1079 CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, ev, sev, 0.0_dp, tmat)
1080 END associate
1081
1082 CALL choose_eigv_solver(tmat, teig(ispin), eigenvalues(1:nmo, ispin))
1083
1084 CALL cp_fm_struct_release(fm_mo_struct)
1085 CALL cp_fm_release(tmat)
1086 CALL cp_fm_release(sev)
1087 END DO
1088 ! find major determinants i->a
1089 ia_index = 0
1090 sume = 0.0_dp
1091 nnto = 0
1092 DO i = 1, ntomax
1093 iv = maxloc(eigenvalues)
1094 ia_eval(i) = eigenvalues(iv(1), iv(2))
1095 ia_index(1:2, i) = iv(1:2)
1096 sume = sume + ia_eval(i)
1097 eigenvalues(iv(1), iv(2)) = 0.0_dp
1098 nnto = nnto + 1
1099 IF (sume > threshold) EXIT
1100 END DO
1101 ! store hole states
1102 CALL set_mo_set(nto_set(1), nmo=nnto)
1103 DO i = 1, nnto
1104 ia = ia_index(1, i)
1105 ispin = ia_index(2, i)
1106 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=nmo)
1107 CALL cp_fm_get_info(teig(ispin), matrix_struct=fm_struct)
1108 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1109 nrow_global=nmo, ncol_global=1)
1110 CALL cp_fm_create(tmat, fm_mo_struct)
1111 CALL cp_fm_struct_release(fm_mo_struct)
1112 CALL cp_fm_get_info(gs_mos(1)%mos_occ, matrix_struct=fm_struct)
1113 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1114 ncol_global=1)
1115 CALL cp_fm_create(wvec, fm_mo_struct)
1116 CALL cp_fm_struct_release(fm_mo_struct)
1117 CALL cp_fm_to_fm(teig(ispin), tmat, 1, ia, 1)
1118 CALL parallel_gemm('N', 'N', nao, 1, nmo, 1.0_dp, gs_mos(ispin)%mos_occ, &
1119 tmat, 0.0_dp, wvec)
1120 CALL cp_fm_to_fm(wvec, nto_set(1)%mo_coeff, 1, 1, i)
1121 CALL cp_fm_release(wvec)
1122 CALL cp_fm_release(tmat)
1123 END DO
1124 ! particle states
1125 ! Solve generalized eigenvalue equation: (S*X)*(S*X)(T)*v = lambda*S*v
1126 CALL set_mo_set(nto_set(2), nmo=nnto)
1127 DO ispin = 1, nspins
1128 associate(ev => evects(ispin, istate))
1129 CALL cp_fm_get_info(ev, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
1130 ALLOCATE (eigvals(nao))
1131 eigvals = 0.0_dp
1132 CALL cp_fm_create(sev, fm_struct)
1133 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1134 END associate
1135 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1136 nrow_global=nao, ncol_global=nao)
1137 CALL cp_fm_create(tmat, fm_mo_struct)
1138 CALL cp_fm_create(smat, fm_mo_struct)
1139 CALL cp_fm_create(wmat, fm_mo_struct)
1140 CALL cp_fm_create(work, fm_mo_struct)
1141 CALL cp_fm_struct_release(fm_mo_struct)
1142 CALL copy_dbcsr_to_fm(matrix_s, smat)
1143 CALL parallel_gemm('N', 'T', nao, nao, nmo, 1.0_dp, sev, sev, 0.0_dp, tmat)
1144 CALL cp_fm_geeig(tmat, smat, wmat, eigvals, work)
1145 DO i = 1, nnto
1146 IF (ispin == ia_index(2, i)) THEN
1147 icg = 0
1148 DO j = 1, nao
1149 IF (abs(eigvals(j) - ia_eval(i)) < 1.e-6_dp) THEN
1150 icg = j
1151 EXIT
1152 END IF
1153 END DO
1154 IF (icg == 0) THEN
1155 CALL cp_warn(__location__, &
1156 "Could not locate particle state associated with hole state.")
1157 ELSE
1158 CALL cp_fm_to_fm(wmat, nto_set(2)%mo_coeff, 1, icg, i)
1159 END IF
1160 END IF
1161 END DO
1162 DEALLOCATE (eigvals)
1163 CALL cp_fm_release(sev)
1164 CALL cp_fm_release(tmat)
1165 CALL cp_fm_release(smat)
1166 CALL cp_fm_release(wmat)
1167 CALL cp_fm_release(work)
1168 END DO
1169 ! print
1170 IF (iounit > 0) THEN
1171 sume = 0.0_dp
1172 DO i = 1, nnto
1173 sume = sume + ia_eval(i)
1174 WRITE (iounit, "(T6,A,i2,T30,A,i1,T42,A,F8.5,T63,A,F8.5)") &
1175 "Particle-Hole state:", i, " Spin:", ia_index(2, i), &
1176 "Eigenvalue:", ia_eval(i), " Sum Eigv:", sume
1177 END DO
1178 END IF
1179 ! Cube and Molden files
1180 nto_section => section_vals_get_subs_vals(print_section, "NTO_ANALYSIS")
1181 CALL section_vals_val_get(nto_section, "CUBE_FILES", l_val=cube_file)
1182 CALL section_vals_val_get(nto_section, "STRIDE", i_vals=stride)
1183 CALL section_vals_val_get(nto_section, "APPEND", l_val=append_cube)
1184 IF (cube_file) THEN
1185 CALL print_nto_cubes(qs_env, nto_set, istate, stride, append_cube, nto_section)
1186 END IF
1187 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, cell=cell)
1188 molden_section => section_vals_get_subs_vals(print_section, "MOS_MOLDEN")
1189 CALL write_mos_molden(nto_set, qs_kind_set, particle_set, molden_section, cell=cell, qs_env=qs_env)
1190 !
1191 DEALLOCATE (eigenvalues)
1192 CALL cp_fm_release(teig)
1193 !
1194 DO i = 1, 2
1195 CALL deallocate_mo_set(nto_set(i))
1196 END DO
1197 DEALLOCATE (nto_set)
1198 END DO
1199
1200 IF (iounit > 0) THEN
1201 WRITE (iounit, "(1X,A)") &
1202 "-------------------------------------------------------------------------------"
1203 END IF
1204
1205 END IF
1206
1207 CALL timestop(handle)
1208
1209 END SUBROUTINE tddfpt_print_nto_analysis
1210
1211! **************************************************************************************************
1212!> \brief Print exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018)
1213!> \param log_unit output unit
1214!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
1215!> SIZE(evects,2) -- number of excited states to print)
1216!> \param gs_mos molecular orbitals optimised for the ground state
1217!> \param matrix_s overlap matrix
1218!> \param do_directional_exciton_descriptors flag for computing descriptors for each (cartesian) direction
1219!> \param qs_env Information on particles/geometry
1220!> \par History
1221!> * 12.2024 created as 'tddfpt_print_exciton_descriptors' [Maximilian Graml]
1222! **************************************************************************************************
1223 SUBROUTINE tddfpt_print_exciton_descriptors(log_unit, evects, gs_mos, matrix_s, &
1224 do_directional_exciton_descriptors, qs_env)
1225 INTEGER, INTENT(in) :: log_unit
1226 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
1227 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1228 INTENT(in) :: gs_mos
1229 TYPE(dbcsr_type), POINTER :: matrix_s
1230 LOGICAL, INTENT(IN) :: do_directional_exciton_descriptors
1231 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1232
1233 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_exciton_descriptors'
1234
1235 CHARACTER(LEN=4) :: prefix_output
1236 INTEGER :: handle, ispin, istate, n_moments_quad, &
1237 nactive, nao, nspins, nstates
1238 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt
1239 LOGICAL :: print_checkvalue
1240 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ref_point_multipole
1241 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1242 TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo_coeff, &
1243 fm_struct_s_mos_virt, fm_struct_x_ia_n
1244 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: eigvec_x_ia_n, fm_multipole_ab, &
1245 fm_multipole_ai, fm_multipole_ij, &
1246 s_mos_virt
1247 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_coeff
1248 TYPE(exciton_descr_type), ALLOCATABLE, &
1249 DIMENSION(:) :: exc_descr
1250
1251 CALL timeset(routinen, handle)
1252
1253 nspins = SIZE(evects, 1)
1254 nstates = SIZE(evects, 2)
1255
1256 cpassert(nspins == 1) ! Other spins are not yet implemented for exciton descriptors
1257
1258 CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env)
1259 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1260
1261 DO ispin = 1, nspins
1262 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
1263 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
1264 END DO
1265
1266 ! Prepare fm with all MO coefficents, i.e. nao x nao
1267 ALLOCATE (mo_coeff(nspins))
1268 CALL cp_fm_struct_create(fm_struct_mo_coeff, nrow_global=nao, ncol_global=nao, &
1269 context=blacs_env)
1270 DO ispin = 1, nspins
1271 CALL cp_fm_create(mo_coeff(ispin), fm_struct_mo_coeff)
1272 CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_occ, &
1273 mo_coeff(ispin), &
1274 nao, &
1275 nmo_occ(ispin), &
1276 1, &
1277 1, &
1278 1, &
1279 1, &
1280 blacs_env)
1281 CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_virt, &
1282 mo_coeff(ispin), &
1283 nao, &
1284 nmo_virt(ispin), &
1285 1, &
1286 1, &
1287 1, &
1288 nmo_occ(ispin) + 1, &
1289 blacs_env)
1290 END DO
1291 CALL cp_fm_struct_release(fm_struct_mo_coeff)
1292
1293 ! Compute multipole moments
1294 ! fm_multipole_XY have structure inherited by libint, i.e. x, y, z, xx, xy, xz, yy, yz, zz
1295 n_moments_quad = 9
1296 ALLOCATE (ref_point_multipole(3))
1297 ALLOCATE (fm_multipole_ij(n_moments_quad))
1298 ALLOCATE (fm_multipole_ab(n_moments_quad))
1299 ALLOCATE (fm_multipole_ai(n_moments_quad))
1300
1301 CALL get_multipoles_mo(fm_multipole_ai, fm_multipole_ij, fm_multipole_ab, &
1302 qs_env, mo_coeff, ref_point_multipole, 2, &
1303 nmo_occ(1), nmo_virt(1), blacs_env)
1304
1305 CALL cp_fm_release(mo_coeff)
1306
1307 ! Compute eigenvector X of the Casida equation from trial vectors
1308 ALLOCATE (s_mos_virt(nspins), eigvec_x_ia_n(nspins))
1309 DO ispin = 1, nspins
1310 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct_s_mos_virt)
1311 CALL cp_fm_create(s_mos_virt(ispin), fm_struct_s_mos_virt)
1312 NULLIFY (fm_struct_s_mos_virt)
1313 CALL cp_dbcsr_sm_fm_multiply(matrix_s, &
1314 gs_mos(ispin)%mos_virt, &
1315 s_mos_virt(ispin), &
1316 ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
1317
1318 CALL cp_fm_struct_create(fm_struct_x_ia_n, nrow_global=nmo_occ(ispin), ncol_global=nmo_virt(ispin), &
1319 context=blacs_env)
1320 CALL cp_fm_create(eigvec_x_ia_n(ispin), fm_struct_x_ia_n)
1321 CALL cp_fm_struct_release(fm_struct_x_ia_n)
1322 END DO
1323 ALLOCATE (exc_descr(nstates))
1324 DO istate = 1, nstates
1325 DO ispin = 1, nspins
1326 CALL cp_fm_set_all(eigvec_x_ia_n(ispin), 0.0_dp)
1327 ! compute eigenvectors X of the TDA equation
1328 ! Reshuffle multiplication from
1329 ! X_ai = S_ma ^T * C_mi
1330 ! to
1331 ! X_ia = C_mi ^T * S_ma
1332 ! for compatibility with the structure needed for get_exciton_descriptors of bse_properties.F
1333 CALL cp_fm_get_info(evects(ispin, istate), ncol_global=nactive)
1334 IF (nactive /= nmo_occ(ispin)) THEN
1335 CALL cp_abort(__location__, &
1336 "Reduced active space excitations not implemented")
1337 END IF
1338 CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_virt(ispin), nao, 1.0_dp, &
1339 evects(ispin, istate), s_mos_virt(ispin), 0.0_dp, eigvec_x_ia_n(ispin))
1340
1341 CALL get_exciton_descriptors(exc_descr, eigvec_x_ia_n(ispin), &
1342 fm_multipole_ij, fm_multipole_ab, &
1343 fm_multipole_ai, &
1344 istate, nmo_occ(ispin), nmo_virt(ispin))
1345 END DO
1346 END DO
1347 CALL cp_fm_release(eigvec_x_ia_n)
1348 CALL cp_fm_release(s_mos_virt)
1349 CALL cp_fm_release(fm_multipole_ai)
1350 CALL cp_fm_release(fm_multipole_ij)
1351 CALL cp_fm_release(fm_multipole_ab)
1352
1353 ! Actual printing
1354 print_checkvalue = .true.
1355 prefix_output = ' '
1356 CALL print_exciton_descriptors(exc_descr, ref_point_multipole, log_unit, &
1357 nstates, print_checkvalue, do_directional_exciton_descriptors, &
1358 prefix_output, qs_env)
1359
1360 DEALLOCATE (ref_point_multipole)
1361 DEALLOCATE (exc_descr)
1362
1363 CALL timestop(handle)
1364
1366
1367! **************************************************************************************************
1368!> \brief ...
1369!> \param vin ...
1370!> \param vout ...
1371!> \param mos_occ ...
1372!> \param matrix_s ...
1373! **************************************************************************************************
1374 SUBROUTINE project_vector(vin, vout, mos_occ, matrix_s)
1375 TYPE(dbcsr_type) :: vin, vout
1376 TYPE(cp_fm_type), INTENT(IN) :: mos_occ
1377 TYPE(dbcsr_type), POINTER :: matrix_s
1378
1379 CHARACTER(LEN=*), PARAMETER :: routinen = 'project_vector'
1380
1381 INTEGER :: handle, nao, nmo
1382 REAL(kind=dp) :: norm(1)
1383 TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_vec_struct
1384 TYPE(cp_fm_type) :: csvec, svec, vec
1385
1386 CALL timeset(routinen, handle)
1387
1388 CALL cp_fm_get_info(mos_occ, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
1389 CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
1390 nrow_global=nao, ncol_global=1)
1391 CALL cp_fm_create(vec, fm_vec_struct)
1392 CALL cp_fm_create(svec, fm_vec_struct)
1393 CALL cp_fm_struct_release(fm_vec_struct)
1394 CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
1395 nrow_global=nmo, ncol_global=1)
1396 CALL cp_fm_create(csvec, fm_vec_struct)
1397 CALL cp_fm_struct_release(fm_vec_struct)
1398
1399 CALL copy_dbcsr_to_fm(vin, vec)
1400 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vec, svec, ncol=1, alpha=1.0_dp, beta=0.0_dp)
1401 CALL parallel_gemm('T', 'N', nmo, 1, nao, 1.0_dp, mos_occ, svec, 0.0_dp, csvec)
1402 CALL parallel_gemm('N', 'N', nao, 1, nmo, -1.0_dp, mos_occ, csvec, 1.0_dp, vec)
1403 CALL cp_fm_vectorsnorm(vec, norm)
1404 cpassert(norm(1) > 1.e-14_dp)
1405 norm(1) = sqrt(1._dp/norm(1))
1406 CALL cp_fm_scale(norm(1), vec)
1407 CALL copy_fm_to_dbcsr(vec, vout, keep_sparsity=.false.)
1408
1409 CALL cp_fm_release(csvec)
1410 CALL cp_fm_release(svec)
1411 CALL cp_fm_release(vec)
1412
1413 CALL timestop(handle)
1414
1415 END SUBROUTINE project_vector
1416
1417! **************************************************************************************************
1418!> \brief ...
1419!> \param va ...
1420!> \param vb ...
1421!> \param res ...
1422! **************************************************************************************************
1423 SUBROUTINE vec_product(va, vb, res)
1424 TYPE(dbcsr_type) :: va, vb
1425 REAL(kind=dp), INTENT(OUT) :: res
1426
1427 CHARACTER(LEN=*), PARAMETER :: routinen = 'vec_product'
1428
1429 INTEGER :: handle, icol, irow
1430 LOGICAL :: found
1431 REAL(kind=dp), DIMENSION(:, :), POINTER :: vba, vbb
1432 TYPE(dbcsr_iterator_type) :: iter
1433 TYPE(mp_comm_type) :: group
1434
1435 CALL timeset(routinen, handle)
1436
1437 res = 0.0_dp
1438
1439 CALL dbcsr_get_info(va, group=group)
1440 CALL dbcsr_iterator_start(iter, va)
1441 DO WHILE (dbcsr_iterator_blocks_left(iter))
1442 CALL dbcsr_iterator_next_block(iter, irow, icol, vba)
1443 CALL dbcsr_get_block_p(vb, row=irow, col=icol, block=vbb, found=found)
1444 res = res + sum(vba*vbb)
1445 cpassert(found)
1446 END DO
1447 CALL dbcsr_iterator_stop(iter)
1448 CALL group%sum(res)
1449
1450 CALL timestop(handle)
1451
1452 END SUBROUTINE vec_product
1453
1454! **************************************************************************************************
1455!> \brief ...
1456!> \param qs_env ...
1457!> \param mos ...
1458!> \param istate ...
1459!> \param stride ...
1460!> \param append_cube ...
1461!> \param print_section ...
1462! **************************************************************************************************
1463 SUBROUTINE print_nto_cubes(qs_env, mos, istate, stride, append_cube, print_section)
1464
1465 TYPE(qs_environment_type), POINTER :: qs_env
1466 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1467 INTEGER, INTENT(IN) :: istate
1468 INTEGER, DIMENSION(:), POINTER :: stride
1469 LOGICAL, INTENT(IN) :: append_cube
1470 TYPE(section_vals_type), POINTER :: print_section
1471
1472 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
1473 INTEGER :: i, iset, nmo, unit_nr
1474 LOGICAL :: mpi_io
1475 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1476 TYPE(cell_type), POINTER :: cell
1477 TYPE(cp_fm_type), POINTER :: mo_coeff
1478 TYPE(cp_logger_type), POINTER :: logger
1479 TYPE(dft_control_type), POINTER :: dft_control
1480 TYPE(particle_list_type), POINTER :: particles
1481 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1482 TYPE(pw_c1d_gs_type) :: wf_g
1483 TYPE(pw_env_type), POINTER :: pw_env
1484 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1485 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1486 TYPE(pw_r3d_rs_type) :: wf_r
1487 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1488 TYPE(qs_subsys_type), POINTER :: subsys
1489
1490 logger => cp_get_default_logger()
1491
1492 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
1493 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1494 CALL auxbas_pw_pool%create_pw(wf_r)
1495 CALL auxbas_pw_pool%create_pw(wf_g)
1496
1497 CALL get_qs_env(qs_env, subsys=subsys)
1498 CALL qs_subsys_get(subsys, particles=particles)
1499
1500 my_pos_cube = "REWIND"
1501 IF (append_cube) THEN
1502 my_pos_cube = "APPEND"
1503 END IF
1504
1505 CALL get_qs_env(qs_env=qs_env, &
1506 atomic_kind_set=atomic_kind_set, &
1507 qs_kind_set=qs_kind_set, &
1508 cell=cell, &
1509 particle_set=particle_set)
1510
1511 DO iset = 1, 2
1512 CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo)
1513 DO i = 1, nmo
1514 CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1515 cell, dft_control, particle_set, pw_env)
1516 IF (iset == 1) THEN
1517 WRITE (filename, '(a4,I3.3,I2.2,a11)') "NTO_STATE", istate, i, "_Hole_State"
1518 ELSEIF (iset == 2) THEN
1519 WRITE (filename, '(a4,I3.3,I2.2,a15)') "NTO_STATE", istate, i, "_Particle_State"
1520 END IF
1521 mpi_io = .true.
1522 unit_nr = cp_print_key_unit_nr(logger, print_section, '', extension=".cube", &
1523 middle_name=trim(filename), file_position=my_pos_cube, &
1524 log_filename=.false., ignore_should_output=.true., mpi_io=mpi_io)
1525 IF (iset == 1) THEN
1526 WRITE (title, *) "Natural Transition Orbital Hole State", i
1527 ELSEIF (iset == 2) THEN
1528 WRITE (title, *) "Natural Transition Orbital Particle State", i
1529 END IF
1530 CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, stride=stride, mpi_io=mpi_io)
1531 CALL cp_print_key_finished_output(unit_nr, logger, print_section, '', &
1532 ignore_should_output=.true., mpi_io=mpi_io)
1533 END DO
1534 END DO
1535
1536 CALL auxbas_pw_pool%give_back_pw(wf_g)
1537 CALL auxbas_pw_pool%give_back_pw(wf_r)
1538
1539 END SUBROUTINE print_nto_cubes
1540
1541END 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, qs_env, calc_energies)
Write out the MOs in molden format for visualisation.
Calculates the moment integrals <a|r^m|b>
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
basic linear algebra operations for full matrixes
represent a simple array based list of the given type
Define the data structure for the particle information.
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public evolt
Definition physcon.F:183
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type)
maps a given wavefunction on the grid
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public set_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, uniform_occupation, kts, mu, flexible_electron_count)
Set the components of a MO set data structure.
subroutine, public 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.