(git:4e96ea1)
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
607 CALL cp_fm_get_info(evects(ispin, 1), matrix_struct=matrix_struct)
608 CALL cp_fm_create(dipact, matrix_struct)
609 DO ideriv = 1, nderivs
610 DO i = 1, nactive
611 j = gs_mos(ispin)%index_active(i)
612 CALL cp_fm_to_fm(dipole_op_mos_occ(ideriv, ispin), dipact, &
613 ncol=1, source_start=j, target_start=i)
614 END DO
615 CALL cp_fm_trace(evects(ispin, :), dipact, trans_dipoles(:, ideriv, ispin))
616 END DO
617 CALL cp_fm_release(dipact)
618 END IF
619 END DO
620
621 IF (nspins == 1) THEN
622 trans_dipoles(:, :, 1) = sqrt(2.0_dp)*trans_dipoles(:, :, 1)
623 ELSE
624 trans_dipoles(:, :, 1) = sqrt(trans_dipoles(:, :, 1)**2 + trans_dipoles(:, :, 2)**2)
625 END IF
626 END IF
627
628 ! *** summary information ***
629 DO istate = 1, nstates
630 osc_strength = 2.0_dp/3.0_dp*evals(istate)* &
631 accurate_dot_product(trans_dipoles(istate, :, 1), trans_dipoles(istate, :, 1))
632 ostrength(istate) = osc_strength
633 IF (log_unit > 0) THEN
634 WRITE (log_unit, '(1X,A,T9,I7,T19,F11.5,T31,3(1X,ES11.4E2),T69,ES12.5E2)') &
635 "TDDFPT|", istate, evals(istate)*evolt, trans_dipoles(istate, 1:nderivs, 1), osc_strength
636 END IF
637 END DO
638
639 ! punch a checksum for the regs
640 IF (log_unit > 0) THEN
641 WRITE (log_unit, '(/,T2,A,E14.6)') 'TDDFPT : CheckSum =', sqrt(sum(evals**2))
642 END IF
643
644 DEALLOCATE (trans_dipoles)
645
646 CALL timestop(handle)
647 END SUBROUTINE tddfpt_print_summary
648
649! **************************************************************************************************
650!> \brief Print excitation analysis.
651!> \param log_unit output unit
652!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
653!> SIZE(evects,2) -- number of excited states to print)
654!> \param evals TDDFPT eigenvalues
655!> \param gs_mos molecular orbitals optimised for the ground state
656!> \param matrix_s overlap matrix
657!> \param spinflip ...
658!> \param min_amplitude the smallest excitation amplitude to print
659!> \par History
660!> * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov]
661!> * 08.2018 splited of from 'tddfpt_print_summary' [Sergey Chulkov]
662! **************************************************************************************************
663 SUBROUTINE tddfpt_print_excitation_analysis(log_unit, evects, evals, gs_mos, matrix_s, spinflip, &
664 min_amplitude)
665 INTEGER, INTENT(in) :: log_unit
666 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
667 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
668 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
669 INTENT(in) :: gs_mos
670 TYPE(dbcsr_type), POINTER :: matrix_s
671 INTEGER :: spinflip
672 REAL(kind=dp), INTENT(in) :: min_amplitude
673
674 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_excitation_analysis'
675
676 CHARACTER(len=5) :: spin_label, spin_label2
677 INTEGER :: handle, icol, iproc, irow, ispin, &
678 istate, nao, ncols_local, nrows_local, &
679 nspins, nstates, spin2, state_spin, &
680 state_spin2
681 INTEGER(kind=int_8) :: iexc, imo_act, imo_occ, imo_virt, ind, &
682 nexcs, nexcs_local, nexcs_max_local, &
683 nmo_virt_occ, nmo_virt_occ_alpha
684 INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:) :: inds_local, inds_recv, nexcs_recv
685 INTEGER(kind=int_8), DIMENSION(1) :: nexcs_send
686 INTEGER(kind=int_8), DIMENSION(maxspins) :: nactive8, nmo_occ8, nmo_virt8
687 INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
688 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
689 INTEGER, DIMENSION(maxspins) :: nactive, nmo_occ, nmo_virt
690 LOGICAL :: do_exc_analysis
691 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights_local, weights_neg_abs_recv, &
692 weights_recv
693 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
694 POINTER :: local_data
695 TYPE(cp_blacs_env_type), POINTER :: blacs_env
696 TYPE(cp_fm_struct_type), POINTER :: fm_struct
697 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: s_mos_virt, weights_fm
698 TYPE(mp_para_env_type), POINTER :: para_env
699 TYPE(mp_request_type) :: send_handler, send_handler2
700 TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_handlers, recv_handlers2
701
702 CALL timeset(routinen, handle)
703
704 nspins = SIZE(gs_mos, 1)
705 nstates = SIZE(evects, 2)
706 do_exc_analysis = min_amplitude < 1.0_dp
707
708 CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env, para_env=para_env)
709 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
710
711 DO ispin = 1, nspins
712 nactive(ispin) = gs_mos(ispin)%nmo_active
713 nactive8(ispin) = int(nactive(ispin), kind=int_8)
714 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
715 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
716 nmo_occ8(ispin) = SIZE(gs_mos(ispin)%evals_occ, kind=int_8)
717 nmo_virt8(ispin) = SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
718 END DO
719
720 ! *** excitation analysis ***
721 IF (do_exc_analysis) THEN
722 cpassert(log_unit <= 0 .OR. para_env%is_source())
723 nmo_virt_occ_alpha = int(nmo_virt(1), int_8)*int(nmo_occ(1), int_8)
724
725 IF (log_unit > 0) THEN
726 WRITE (log_unit, "(1X,A)") "", &
727 "-------------------------------------------------------------------------------", &
728 "- Excitation analysis -", &
729 "-------------------------------------------------------------------------------"
730 WRITE (log_unit, '(8X,A,T27,A,T49,A,T69,A)') "State", "Occupied", "Virtual", "Excitation"
731 WRITE (log_unit, '(8X,A,T28,A,T49,A,T69,A)') "number", "orbital", "orbital", "amplitude"
732 WRITE (log_unit, '(1X,79("-"))')
733
734 IF (nspins == 1) THEN
735 state_spin = 1
736 state_spin2 = 2
737 spin_label = ' '
738 spin_label2 = ' '
739 ELSE IF (spinflip /= no_sf_tddfpt) THEN
740 state_spin = 1
741 state_spin2 = 2
742 spin_label = '(alp)'
743 spin_label2 = '(bet)'
744 END IF
745 END IF
746
747 ALLOCATE (s_mos_virt(SIZE(evects, 1)), weights_fm(SIZE(evects, 1)))
748 DO ispin = 1, SIZE(evects, 1)
749 IF (spinflip == no_sf_tddfpt) THEN
750 spin2 = ispin
751 ELSE
752 spin2 = 2
753 END IF
754 CALL cp_fm_get_info(gs_mos(spin2)%mos_virt, matrix_struct=fm_struct)
755 CALL cp_fm_create(s_mos_virt(ispin), fm_struct)
756 CALL cp_dbcsr_sm_fm_multiply(matrix_s, &
757 gs_mos(spin2)%mos_virt, &
758 s_mos_virt(ispin), &
759 ncol=nmo_virt(spin2), alpha=1.0_dp, beta=0.0_dp)
760
761 NULLIFY (fm_struct)
762 CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(spin2), ncol_global=nactive(ispin), &
763 context=blacs_env)
764 CALL cp_fm_create(weights_fm(ispin), fm_struct)
765 CALL cp_fm_set_all(weights_fm(ispin), 0.0_dp)
766 CALL cp_fm_struct_release(fm_struct)
767 END DO
768
769 nexcs_max_local = 0
770 DO ispin = 1, SIZE(evects, 1)
771 CALL cp_fm_get_info(weights_fm(ispin), nrow_local=nrows_local, ncol_local=ncols_local)
772 nexcs_max_local = nexcs_max_local + int(nrows_local, int_8)*int(ncols_local, int_8)
773 END DO
774
775 ALLOCATE (weights_local(nexcs_max_local), inds_local(nexcs_max_local))
776
777 DO istate = 1, nstates
778 nexcs_local = 0
779 nmo_virt_occ = 0
780
781 ! analyse matrix elements locally and transfer only significant
782 ! excitations to the master node for subsequent ordering
783 DO ispin = 1, SIZE(evects, 1)
784 IF (spinflip == no_sf_tddfpt) THEN
785 spin2 = ispin
786 ELSE
787 spin2 = 2
788 END IF
789 ! compute excitation amplitudes
790 CALL parallel_gemm('T', 'N', nmo_virt(spin2), nactive(ispin), nao, 1.0_dp, s_mos_virt(ispin), &
791 evects(ispin, istate), 0.0_dp, weights_fm(ispin))
792
793 CALL cp_fm_get_info(weights_fm(ispin), nrow_local=nrows_local, ncol_local=ncols_local, &
794 row_indices=row_indices, col_indices=col_indices, local_data=local_data)
795
796 ! locate single excitations with significant amplitudes (>= min_amplitude)
797 DO icol = 1, ncols_local
798 DO irow = 1, nrows_local
799 IF (abs(local_data(irow, icol)) >= min_amplitude) THEN
800 ! number of non-negligible excitations
801 nexcs_local = nexcs_local + 1
802 ! excitation amplitude
803 weights_local(nexcs_local) = local_data(irow, icol)
804 ! index of single excitation (ivirt, iocc, ispin) in compressed form
805 inds_local(nexcs_local) = nmo_virt_occ + int(row_indices(irow), int_8) + &
806 int(col_indices(icol) - 1, int_8)*nmo_virt8(spin2)
807 END IF
808 END DO
809 END DO
810
811 nmo_virt_occ = nmo_virt_occ + nmo_virt8(spin2)*nmo_occ8(ispin)
812 END DO
813
814 IF (para_env%is_source()) THEN
815 ! master node
816 ALLOCATE (nexcs_recv(para_env%num_pe), recv_handlers(para_env%num_pe), recv_handlers2(para_env%num_pe))
817
818 ! collect number of non-negligible excitations from other nodes
819 DO iproc = 1, para_env%num_pe
820 IF (iproc - 1 /= para_env%mepos) THEN
821 CALL para_env%irecv(nexcs_recv(iproc:iproc), iproc - 1, recv_handlers(iproc), 0)
822 ELSE
823 nexcs_recv(iproc) = nexcs_local
824 END IF
825 END DO
826
827 DO iproc = 1, para_env%num_pe
828 IF (iproc - 1 /= para_env%mepos) &
829 CALL recv_handlers(iproc)%wait()
830 END DO
831
832 ! compute total number of non-negligible excitations
833 nexcs = 0
834 DO iproc = 1, para_env%num_pe
835 nexcs = nexcs + nexcs_recv(iproc)
836 END DO
837
838 ! receive indices and amplitudes of selected excitations
839 ALLOCATE (weights_recv(nexcs), weights_neg_abs_recv(nexcs))
840 ALLOCATE (inds_recv(nexcs), inds(nexcs))
841
842 nmo_virt_occ = 0
843 DO iproc = 1, para_env%num_pe
844 IF (nexcs_recv(iproc) > 0) THEN
845 IF (iproc - 1 /= para_env%mepos) THEN
846 ! excitation amplitudes
847 CALL para_env%irecv(weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), &
848 iproc - 1, recv_handlers(iproc), 1)
849 ! compressed indices
850 CALL para_env%irecv(inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), &
851 iproc - 1, recv_handlers2(iproc), 2)
852 ELSE
853 ! data on master node
854 weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = weights_local(1:nexcs_recv(iproc))
855 inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = inds_local(1:nexcs_recv(iproc))
856 END IF
857
858 nmo_virt_occ = nmo_virt_occ + nexcs_recv(iproc)
859 END IF
860 END DO
861
862 DO iproc = 1, para_env%num_pe
863 IF (iproc - 1 /= para_env%mepos .AND. nexcs_recv(iproc) > 0) THEN
864 CALL recv_handlers(iproc)%wait()
865 CALL recv_handlers2(iproc)%wait()
866 END IF
867 END DO
868
869 DEALLOCATE (nexcs_recv, recv_handlers, recv_handlers2)
870 ELSE
871 ! working node: send the number of selected excited states to the master node
872 nexcs_send(1) = nexcs_local
873 CALL para_env%isend(nexcs_send, para_env%source, send_handler, 0)
874 CALL send_handler%wait()
875
876 IF (nexcs_local > 0) THEN
877 ! send excitation amplitudes
878 CALL para_env%isend(weights_local(1:nexcs_local), para_env%source, send_handler, 1)
879 ! send compressed indices
880 CALL para_env%isend(inds_local(1:nexcs_local), para_env%source, send_handler2, 2)
881
882 CALL send_handler%wait()
883 CALL send_handler2%wait()
884 END IF
885 END IF
886
887 ! sort non-negligible excitations on the master node according to their amplitudes,
888 ! uncompress indices and print summary information
889 IF (para_env%is_source() .AND. log_unit > 0) THEN
890 weights_neg_abs_recv(:) = -abs(weights_recv)
891 CALL sort(weights_neg_abs_recv, int(nexcs), inds)
892
893 WRITE (log_unit, '(T7,I8,F10.5,A)') istate, evals(istate)*evolt, " eV"
894
895 ! This reinitialization is needed to prevent the intel fortran compiler from introduce
896 ! a bug when using optimization level 3 flag
897 state_spin = 1
898 state_spin2 = 1
899 IF (spinflip /= no_sf_tddfpt) THEN
900 state_spin = 1
901 state_spin2 = 2
902 END IF
903 DO iexc = 1, nexcs
904 ind = inds_recv(inds(iexc)) - 1
905 IF ((nspins > 1) .AND. (spinflip == no_sf_tddfpt)) THEN
906 IF (ind < nmo_virt_occ_alpha) THEN
907 spin_label = '(alp)'
908 spin_label2 = '(alp)'
909 ELSE
910 state_spin = 2
911 state_spin2 = 2
912 ind = ind - nmo_virt_occ_alpha
913 spin_label = '(bet)'
914 spin_label2 = '(bet)'
915 END IF
916 END IF
917 imo_act = ind/nmo_virt8(state_spin2) + 1
918 imo_occ = gs_mos(state_spin)%index_active(imo_act)
919 imo_virt = mod(ind, nmo_virt8(state_spin2)) + 1
920
921 WRITE (log_unit, '(T27,I8,1X,A5,T48,I8,1X,A5,T70,F9.6)') imo_occ, spin_label, &
922 nmo_occ8(state_spin2) + imo_virt, spin_label2, weights_recv(inds(iexc))
923 END DO
924 END IF
925
926 ! deallocate temporary arrays
927 IF (para_env%is_source()) &
928 DEALLOCATE (weights_recv, weights_neg_abs_recv, inds_recv, inds)
929 END DO
930
931 DEALLOCATE (weights_local, inds_local)
932 IF (log_unit > 0) THEN
933 WRITE (log_unit, "(1X,A)") &
934 "-------------------------------------------------------------------------------"
935 END IF
936 END IF
937
938 CALL cp_fm_release(weights_fm)
939 CALL cp_fm_release(s_mos_virt)
940
941 CALL timestop(handle)
942
944
945! **************************************************************************************************
946!> \brief Print natural transition orbital analysis.
947!> \param qs_env Information on Kinds and Particles
948!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
949!> SIZE(evects,2) -- number of excited states to print)
950!> \param evals TDDFPT eigenvalues
951!> \param ostrength ...
952!> \param gs_mos molecular orbitals optimised for the ground state
953!> \param matrix_s overlap matrix
954!> \param print_section ...
955!> \par History
956!> * 06.2019 created [JGH]
957! **************************************************************************************************
958 SUBROUTINE tddfpt_print_nto_analysis(qs_env, evects, evals, ostrength, gs_mos, matrix_s, print_section)
959 TYPE(qs_environment_type), POINTER :: qs_env
960 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
961 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals, ostrength
962 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
963 INTENT(in) :: gs_mos
964 TYPE(dbcsr_type), POINTER :: matrix_s
965 TYPE(section_vals_type), POINTER :: print_section
966
967 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_nto_analysis'
968 INTEGER, PARAMETER :: ntomax = 10
969
970 CHARACTER(LEN=20), DIMENSION(2) :: nto_name
971 INTEGER :: handle, i, ia, icg, iounit, ispin, &
972 istate, j, nao, nlist, nmax, nmo, &
973 nnto, nspins, nstates
974 INTEGER, DIMENSION(2) :: iv
975 INTEGER, DIMENSION(2, ntomax) :: ia_index
976 INTEGER, DIMENSION(:), POINTER :: slist, stride
977 LOGICAL :: append_cube, cube_file, explicit
978 REAL(kind=dp) :: os_threshold, sume, threshold
979 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvals
980 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvalues
981 REAL(kind=dp), DIMENSION(ntomax) :: ia_eval
982 TYPE(cp_fm_struct_type), POINTER :: fm_mo_struct, fm_struct
983 TYPE(cp_fm_type) :: sev, smat, tmat, wmat, work, wvec
984 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: teig
985 TYPE(cp_logger_type), POINTER :: logger
986 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: nto_set
987 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
988 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
989 TYPE(section_vals_type), POINTER :: molden_section, nto_section
990
991 CALL timeset(routinen, handle)
992
993 logger => cp_get_default_logger()
994 iounit = cp_logger_get_default_io_unit(logger)
995
996 IF (btest(cp_print_key_should_output(logger%iter_info, print_section, &
997 "NTO_ANALYSIS"), cp_p_file)) THEN
998
999 CALL cite_reference(martin2003)
1000
1001 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%THRESHOLD", r_val=threshold)
1002 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%INTENSITY_THRESHOLD", r_val=os_threshold)
1003 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", explicit=explicit)
1004
1005 IF (explicit) THEN
1006 CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", i_vals=slist)
1007 nlist = SIZE(slist)
1008 ELSE
1009 nlist = 0
1010 END IF
1011
1012 IF (iounit > 0) THEN
1013 WRITE (iounit, "(1X,A)") "", &
1014 "-------------------------------------------------------------------------------", &
1015 "- Natural Orbital analysis -", &
1016 "-------------------------------------------------------------------------------"
1017 END IF
1018
1019 nspins = SIZE(evects, 1)
1020 nstates = SIZE(evects, 2)
1021 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1022
1023 DO istate = 1, nstates
1024 IF (os_threshold > ostrength(istate)) THEN
1025 IF (iounit > 0) THEN
1026 WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1027 END IF
1028 cycle
1029 END IF
1030 IF (nlist > 0) THEN
1031 IF (.NOT. any(slist == istate)) THEN
1032 IF (iounit > 0) THEN
1033 WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1034 END IF
1035 cycle
1036 END IF
1037 END IF
1038 IF (iounit > 0) THEN
1039 WRITE (iounit, "(1X,A,I6,T30,F10.5,A)") " STATE NR. ", istate, evals(istate)*evolt, " eV"
1040 END IF
1041 nmax = 0
1042 DO ispin = 1, nspins
1043 CALL cp_fm_get_info(evects(ispin, istate), matrix_struct=fm_struct, ncol_global=nmo)
1044 nmax = max(nmax, nmo)
1045 END DO
1046 ALLOCATE (eigenvalues(nmax, nspins))
1047 eigenvalues = 0.0_dp
1048 ! SET 1: Hole states
1049 ! SET 2: Particle states
1050 nto_name(1) = 'Hole_states'
1051 nto_name(2) = 'Particle_states'
1052 ALLOCATE (nto_set(2))
1053 DO i = 1, 2
1054 CALL allocate_mo_set(nto_set(i), nao, ntomax, 0, 0.0_dp, 1.0_dp, 0.0_dp)
1055 CALL cp_fm_get_info(evects(1, istate), matrix_struct=fm_struct)
1056 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1057 ncol_global=ntomax)
1058 CALL cp_fm_create(tmat, fm_mo_struct)
1059 CALL init_mo_set(nto_set(i), fm_ref=tmat, name=nto_name(i))
1060 CALL cp_fm_release(tmat)
1061 CALL cp_fm_struct_release(fm_mo_struct)
1062 END DO
1063 !
1064 ALLOCATE (teig(nspins))
1065 ! hole states
1066 ! Diagonalize X(T)*S*X
1067 DO ispin = 1, nspins
1068 associate(ev => evects(ispin, istate))
1069 CALL cp_fm_get_info(ev, matrix_struct=fm_struct, ncol_global=nmo)
1070 CALL cp_fm_create(sev, fm_struct)
1071 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1072 nrow_global=nmo, ncol_global=nmo)
1073 CALL cp_fm_create(tmat, fm_mo_struct)
1074 CALL cp_fm_create(teig(ispin), fm_mo_struct)
1075 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1076 CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, ev, sev, 0.0_dp, tmat)
1077 END associate
1078
1079 CALL choose_eigv_solver(tmat, teig(ispin), eigenvalues(1:nmo, ispin))
1080
1081 CALL cp_fm_struct_release(fm_mo_struct)
1082 CALL cp_fm_release(tmat)
1083 CALL cp_fm_release(sev)
1084 END DO
1085 ! find major determinants i->a
1086 ia_index = 0
1087 sume = 0.0_dp
1088 nnto = 0
1089 DO i = 1, ntomax
1090 iv = maxloc(eigenvalues)
1091 ia_eval(i) = eigenvalues(iv(1), iv(2))
1092 ia_index(1:2, i) = iv(1:2)
1093 sume = sume + ia_eval(i)
1094 eigenvalues(iv(1), iv(2)) = 0.0_dp
1095 nnto = nnto + 1
1096 IF (sume > threshold) EXIT
1097 END DO
1098 ! store hole states
1099 CALL set_mo_set(nto_set(1), nmo=nnto)
1100 DO i = 1, nnto
1101 ia = ia_index(1, i)
1102 ispin = ia_index(2, i)
1103 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=nmo)
1104 CALL cp_fm_get_info(teig(ispin), matrix_struct=fm_struct)
1105 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1106 nrow_global=nmo, ncol_global=1)
1107 CALL cp_fm_create(tmat, fm_mo_struct)
1108 CALL cp_fm_struct_release(fm_mo_struct)
1109 CALL cp_fm_get_info(gs_mos(1)%mos_occ, matrix_struct=fm_struct)
1110 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1111 ncol_global=1)
1112 CALL cp_fm_create(wvec, fm_mo_struct)
1113 CALL cp_fm_struct_release(fm_mo_struct)
1114 CALL cp_fm_to_fm(teig(ispin), tmat, 1, ia, 1)
1115 CALL parallel_gemm('N', 'N', nao, 1, nmo, 1.0_dp, gs_mos(ispin)%mos_occ, &
1116 tmat, 0.0_dp, wvec)
1117 CALL cp_fm_to_fm(wvec, nto_set(1)%mo_coeff, 1, 1, i)
1118 CALL cp_fm_release(wvec)
1119 CALL cp_fm_release(tmat)
1120 END DO
1121 ! particle states
1122 ! Solve generalized eigenvalue equation: (S*X)*(S*X)(T)*v = lambda*S*v
1123 CALL set_mo_set(nto_set(2), nmo=nnto)
1124 DO ispin = 1, nspins
1125 associate(ev => evects(ispin, istate))
1126 CALL cp_fm_get_info(ev, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
1127 ALLOCATE (eigvals(nao))
1128 eigvals = 0.0_dp
1129 CALL cp_fm_create(sev, fm_struct)
1130 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1131 END associate
1132 CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1133 nrow_global=nao, ncol_global=nao)
1134 CALL cp_fm_create(tmat, fm_mo_struct)
1135 CALL cp_fm_create(smat, fm_mo_struct)
1136 CALL cp_fm_create(wmat, fm_mo_struct)
1137 CALL cp_fm_create(work, fm_mo_struct)
1138 CALL cp_fm_struct_release(fm_mo_struct)
1139 CALL copy_dbcsr_to_fm(matrix_s, smat)
1140 CALL parallel_gemm('N', 'T', nao, nao, nmo, 1.0_dp, sev, sev, 0.0_dp, tmat)
1141 CALL cp_fm_geeig(tmat, smat, wmat, eigvals, work)
1142 DO i = 1, nnto
1143 IF (ispin == ia_index(2, i)) THEN
1144 icg = 0
1145 DO j = 1, nao
1146 IF (abs(eigvals(j) - ia_eval(i)) < 1.e-6_dp) THEN
1147 icg = j
1148 EXIT
1149 END IF
1150 END DO
1151 IF (icg == 0) THEN
1152 CALL cp_warn(__location__, &
1153 "Could not locate particle state associated with hole state.")
1154 ELSE
1155 CALL cp_fm_to_fm(wmat, nto_set(2)%mo_coeff, 1, icg, i)
1156 END IF
1157 END IF
1158 END DO
1159 DEALLOCATE (eigvals)
1160 CALL cp_fm_release(sev)
1161 CALL cp_fm_release(tmat)
1162 CALL cp_fm_release(smat)
1163 CALL cp_fm_release(wmat)
1164 CALL cp_fm_release(work)
1165 END DO
1166 ! print
1167 IF (iounit > 0) THEN
1168 sume = 0.0_dp
1169 DO i = 1, nnto
1170 sume = sume + ia_eval(i)
1171 WRITE (iounit, "(T6,A,i2,T30,A,i1,T42,A,F8.5,T63,A,F8.5)") &
1172 "Particle-Hole state:", i, " Spin:", ia_index(2, i), &
1173 "Eigenvalue:", ia_eval(i), " Sum Eigv:", sume
1174 END DO
1175 END IF
1176 ! Cube and Molden files
1177 nto_section => section_vals_get_subs_vals(print_section, "NTO_ANALYSIS")
1178 CALL section_vals_val_get(nto_section, "CUBE_FILES", l_val=cube_file)
1179 CALL section_vals_val_get(nto_section, "STRIDE", i_vals=stride)
1180 CALL section_vals_val_get(nto_section, "APPEND", l_val=append_cube)
1181 IF (cube_file) THEN
1182 CALL print_nto_cubes(qs_env, nto_set, istate, stride, append_cube, nto_section)
1183 END IF
1184 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set)
1185 molden_section => section_vals_get_subs_vals(print_section, "MOS_MOLDEN")
1186 CALL write_mos_molden(nto_set, qs_kind_set, particle_set, molden_section)
1187 !
1188 DEALLOCATE (eigenvalues)
1189 CALL cp_fm_release(teig)
1190 !
1191 DO i = 1, 2
1192 CALL deallocate_mo_set(nto_set(i))
1193 END DO
1194 DEALLOCATE (nto_set)
1195 END DO
1196
1197 IF (iounit > 0) THEN
1198 WRITE (iounit, "(1X,A)") &
1199 "-------------------------------------------------------------------------------"
1200 END IF
1201
1202 END IF
1203
1204 CALL timestop(handle)
1205
1206 END SUBROUTINE tddfpt_print_nto_analysis
1207
1208! **************************************************************************************************
1209!> \brief Print exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018)
1210!> \param log_unit output unit
1211!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
1212!> SIZE(evects,2) -- number of excited states to print)
1213!> \param gs_mos molecular orbitals optimised for the ground state
1214!> \param matrix_s overlap matrix
1215!> \param do_directional_exciton_descriptors flag for computing descriptors for each (cartesian) direction
1216!> \param qs_env Information on particles/geometry
1217!> \par History
1218!> * 12.2024 created as 'tddfpt_print_exciton_descriptors' [Maximilian Graml]
1219! **************************************************************************************************
1220 SUBROUTINE tddfpt_print_exciton_descriptors(log_unit, evects, gs_mos, matrix_s, &
1221 do_directional_exciton_descriptors, qs_env)
1222 INTEGER, INTENT(in) :: log_unit
1223 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
1224 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1225 INTENT(in) :: gs_mos
1226 TYPE(dbcsr_type), POINTER :: matrix_s
1227 LOGICAL, INTENT(IN) :: do_directional_exciton_descriptors
1228 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1229
1230 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_print_exciton_descriptors'
1231
1232 CHARACTER(LEN=4) :: prefix_output
1233 INTEGER :: handle, ispin, istate, n_moments_quad, &
1234 nactive, nao, nspins, nstates
1235 INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt
1236 LOGICAL :: print_checkvalue
1237 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ref_point_multipole
1238 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1239 TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo_coeff, &
1240 fm_struct_s_mos_virt, fm_struct_x_ia_n
1241 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: eigvec_x_ia_n, fm_multipole_ab, &
1242 fm_multipole_ai, fm_multipole_ij, &
1243 s_mos_virt
1244 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_coeff
1245 TYPE(exciton_descr_type), ALLOCATABLE, &
1246 DIMENSION(:) :: exc_descr
1247
1248 CALL timeset(routinen, handle)
1249
1250 nspins = SIZE(evects, 1)
1251 nstates = SIZE(evects, 2)
1252
1253 cpassert(nspins == 1) ! Other spins are not yet implemented for exciton descriptors
1254
1255 CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env)
1256 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1257
1258 DO ispin = 1, nspins
1259 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
1260 nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
1261 END DO
1262
1263 ! Prepare fm with all MO coefficents, i.e. nao x nao
1264 ALLOCATE (mo_coeff(nspins))
1265 CALL cp_fm_struct_create(fm_struct_mo_coeff, nrow_global=nao, ncol_global=nao, &
1266 context=blacs_env)
1267 DO ispin = 1, nspins
1268 CALL cp_fm_create(mo_coeff(ispin), fm_struct_mo_coeff)
1269 CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_occ, &
1270 mo_coeff(ispin), &
1271 nao, &
1272 nmo_occ(ispin), &
1273 1, &
1274 1, &
1275 1, &
1276 1, &
1277 blacs_env)
1278 CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_virt, &
1279 mo_coeff(ispin), &
1280 nao, &
1281 nmo_virt(ispin), &
1282 1, &
1283 1, &
1284 1, &
1285 nmo_occ(ispin) + 1, &
1286 blacs_env)
1287 END DO
1288 CALL cp_fm_struct_release(fm_struct_mo_coeff)
1289
1290 ! Compute multipole moments
1291 ! fm_multipole_XY have structure inherited by libint, i.e. x, y, z, xx, xy, xz, yy, yz, zz
1292 n_moments_quad = 9
1293 ALLOCATE (ref_point_multipole(3))
1294 ALLOCATE (fm_multipole_ij(n_moments_quad))
1295 ALLOCATE (fm_multipole_ab(n_moments_quad))
1296 ALLOCATE (fm_multipole_ai(n_moments_quad))
1297
1298 CALL get_multipoles_mo(fm_multipole_ai, fm_multipole_ij, fm_multipole_ab, &
1299 qs_env, mo_coeff, ref_point_multipole, 2, &
1300 nmo_occ(1), nmo_virt(1), blacs_env)
1301
1302 CALL cp_fm_release(mo_coeff)
1303
1304 ! Compute eigenvector X of the Casida equation from trial vectors
1305 ALLOCATE (s_mos_virt(nspins), eigvec_x_ia_n(nspins))
1306 DO ispin = 1, nspins
1307 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct_s_mos_virt)
1308 CALL cp_fm_create(s_mos_virt(ispin), fm_struct_s_mos_virt)
1309 NULLIFY (fm_struct_s_mos_virt)
1310 CALL cp_dbcsr_sm_fm_multiply(matrix_s, &
1311 gs_mos(ispin)%mos_virt, &
1312 s_mos_virt(ispin), &
1313 ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
1314
1315 CALL cp_fm_struct_create(fm_struct_x_ia_n, nrow_global=nmo_occ(ispin), ncol_global=nmo_virt(ispin), &
1316 context=blacs_env)
1317 CALL cp_fm_create(eigvec_x_ia_n(ispin), fm_struct_x_ia_n)
1318 CALL cp_fm_struct_release(fm_struct_x_ia_n)
1319 END DO
1320 ALLOCATE (exc_descr(nstates))
1321 DO istate = 1, nstates
1322 DO ispin = 1, nspins
1323 CALL cp_fm_set_all(eigvec_x_ia_n(ispin), 0.0_dp)
1324 ! compute eigenvectors X of the TDA equation
1325 ! Reshuffle multiplication from
1326 ! X_ai = S_ma ^T * C_mi
1327 ! to
1328 ! X_ia = C_mi ^T * S_ma
1329 ! for compatibility with the structure needed for get_exciton_descriptors of bse_properties.F
1330 CALL cp_fm_get_info(evects(ispin, istate), ncol_global=nactive)
1331 IF (nactive /= nmo_occ(ispin)) THEN
1332 CALL cp_abort(__location__, &
1333 "Reduced active space excitations not implemented")
1334 END IF
1335 CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_virt(ispin), nao, 1.0_dp, &
1336 evects(ispin, istate), s_mos_virt(ispin), 0.0_dp, eigvec_x_ia_n(ispin))
1337
1338 CALL get_exciton_descriptors(exc_descr, eigvec_x_ia_n(ispin), &
1339 fm_multipole_ij, fm_multipole_ab, &
1340 fm_multipole_ai, &
1341 istate, nmo_occ(ispin), nmo_virt(ispin))
1342 END DO
1343 END DO
1344 CALL cp_fm_release(eigvec_x_ia_n)
1345 CALL cp_fm_release(s_mos_virt)
1346 CALL cp_fm_release(fm_multipole_ai)
1347 CALL cp_fm_release(fm_multipole_ij)
1348 CALL cp_fm_release(fm_multipole_ab)
1349
1350 ! Actual printing
1351 print_checkvalue = .true.
1352 prefix_output = ' '
1353 CALL print_exciton_descriptors(exc_descr, ref_point_multipole, log_unit, &
1354 nstates, print_checkvalue, do_directional_exciton_descriptors, &
1355 prefix_output, qs_env)
1356
1357 DEALLOCATE (ref_point_multipole)
1358 DEALLOCATE (exc_descr)
1359
1360 CALL timestop(handle)
1361
1363
1364! **************************************************************************************************
1365!> \brief ...
1366!> \param vin ...
1367!> \param vout ...
1368!> \param mos_occ ...
1369!> \param matrix_s ...
1370! **************************************************************************************************
1371 SUBROUTINE project_vector(vin, vout, mos_occ, matrix_s)
1372 TYPE(dbcsr_type) :: vin, vout
1373 TYPE(cp_fm_type), INTENT(IN) :: mos_occ
1374 TYPE(dbcsr_type), POINTER :: matrix_s
1375
1376 CHARACTER(LEN=*), PARAMETER :: routinen = 'project_vector'
1377
1378 INTEGER :: handle, nao, nmo
1379 REAL(kind=dp) :: norm(1)
1380 TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_vec_struct
1381 TYPE(cp_fm_type) :: csvec, svec, vec
1382
1383 CALL timeset(routinen, handle)
1384
1385 CALL cp_fm_get_info(mos_occ, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
1386 CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
1387 nrow_global=nao, ncol_global=1)
1388 CALL cp_fm_create(vec, fm_vec_struct)
1389 CALL cp_fm_create(svec, fm_vec_struct)
1390 CALL cp_fm_struct_release(fm_vec_struct)
1391 CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
1392 nrow_global=nmo, ncol_global=1)
1393 CALL cp_fm_create(csvec, fm_vec_struct)
1394 CALL cp_fm_struct_release(fm_vec_struct)
1395
1396 CALL copy_dbcsr_to_fm(vin, vec)
1397 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vec, svec, ncol=1, alpha=1.0_dp, beta=0.0_dp)
1398 CALL parallel_gemm('T', 'N', nmo, 1, nao, 1.0_dp, mos_occ, svec, 0.0_dp, csvec)
1399 CALL parallel_gemm('N', 'N', nao, 1, nmo, -1.0_dp, mos_occ, csvec, 1.0_dp, vec)
1400 CALL cp_fm_vectorsnorm(vec, norm)
1401 cpassert(norm(1) > 1.e-14_dp)
1402 norm(1) = sqrt(1._dp/norm(1))
1403 CALL cp_fm_scale(norm(1), vec)
1404 CALL copy_fm_to_dbcsr(vec, vout, keep_sparsity=.false.)
1405
1406 CALL cp_fm_release(csvec)
1407 CALL cp_fm_release(svec)
1408 CALL cp_fm_release(vec)
1409
1410 CALL timestop(handle)
1411
1412 END SUBROUTINE project_vector
1413
1414! **************************************************************************************************
1415!> \brief ...
1416!> \param va ...
1417!> \param vb ...
1418!> \param res ...
1419! **************************************************************************************************
1420 SUBROUTINE vec_product(va, vb, res)
1421 TYPE(dbcsr_type) :: va, vb
1422 REAL(kind=dp), INTENT(OUT) :: res
1423
1424 CHARACTER(LEN=*), PARAMETER :: routinen = 'vec_product'
1425
1426 INTEGER :: handle, icol, irow
1427 LOGICAL :: found
1428 REAL(kind=dp), DIMENSION(:, :), POINTER :: vba, vbb
1429 TYPE(dbcsr_iterator_type) :: iter
1430 TYPE(mp_comm_type) :: group
1431
1432 CALL timeset(routinen, handle)
1433
1434 res = 0.0_dp
1435
1436 CALL dbcsr_get_info(va, group=group)
1437 CALL dbcsr_iterator_start(iter, va)
1438 DO WHILE (dbcsr_iterator_blocks_left(iter))
1439 CALL dbcsr_iterator_next_block(iter, irow, icol, vba)
1440 CALL dbcsr_get_block_p(vb, row=irow, col=icol, block=vbb, found=found)
1441 res = res + sum(vba*vbb)
1442 cpassert(found)
1443 END DO
1444 CALL dbcsr_iterator_stop(iter)
1445 CALL group%sum(res)
1446
1447 CALL timestop(handle)
1448
1449 END SUBROUTINE vec_product
1450
1451! **************************************************************************************************
1452!> \brief ...
1453!> \param qs_env ...
1454!> \param mos ...
1455!> \param istate ...
1456!> \param stride ...
1457!> \param append_cube ...
1458!> \param print_section ...
1459! **************************************************************************************************
1460 SUBROUTINE print_nto_cubes(qs_env, mos, istate, stride, append_cube, print_section)
1461
1462 TYPE(qs_environment_type), POINTER :: qs_env
1463 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1464 INTEGER, INTENT(IN) :: istate
1465 INTEGER, DIMENSION(:), POINTER :: stride
1466 LOGICAL, INTENT(IN) :: append_cube
1467 TYPE(section_vals_type), POINTER :: print_section
1468
1469 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
1470 INTEGER :: i, iset, nmo, unit_nr
1471 LOGICAL :: mpi_io
1472 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1473 TYPE(cell_type), POINTER :: cell
1474 TYPE(cp_fm_type), POINTER :: mo_coeff
1475 TYPE(cp_logger_type), POINTER :: logger
1476 TYPE(dft_control_type), POINTER :: dft_control
1477 TYPE(particle_list_type), POINTER :: particles
1478 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1479 TYPE(pw_c1d_gs_type) :: wf_g
1480 TYPE(pw_env_type), POINTER :: pw_env
1481 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1482 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1483 TYPE(pw_r3d_rs_type) :: wf_r
1484 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1485 TYPE(qs_subsys_type), POINTER :: subsys
1486
1487 logger => cp_get_default_logger()
1488
1489 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
1490 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1491 CALL auxbas_pw_pool%create_pw(wf_r)
1492 CALL auxbas_pw_pool%create_pw(wf_g)
1493
1494 CALL get_qs_env(qs_env, subsys=subsys)
1495 CALL qs_subsys_get(subsys, particles=particles)
1496
1497 my_pos_cube = "REWIND"
1498 IF (append_cube) THEN
1499 my_pos_cube = "APPEND"
1500 END IF
1501
1502 CALL get_qs_env(qs_env=qs_env, &
1503 atomic_kind_set=atomic_kind_set, &
1504 qs_kind_set=qs_kind_set, &
1505 cell=cell, &
1506 particle_set=particle_set)
1507
1508 DO iset = 1, 2
1509 CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo)
1510 DO i = 1, nmo
1511 CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1512 cell, dft_control, particle_set, pw_env)
1513 IF (iset == 1) THEN
1514 WRITE (filename, '(a4,I3.3,I2.2,a11)') "NTO_STATE", istate, i, "_Hole_State"
1515 ELSEIF (iset == 2) THEN
1516 WRITE (filename, '(a4,I3.3,I2.2,a15)') "NTO_STATE", istate, i, "_Particle_State"
1517 END IF
1518 mpi_io = .true.
1519 unit_nr = cp_print_key_unit_nr(logger, print_section, '', extension=".cube", &
1520 middle_name=trim(filename), file_position=my_pos_cube, &
1521 log_filename=.false., ignore_should_output=.true., mpi_io=mpi_io)
1522 IF (iset == 1) THEN
1523 WRITE (title, *) "Natural Transition Orbital Hole State", i
1524 ELSEIF (iset == 2) THEN
1525 WRITE (title, *) "Natural Transition Orbital Particle State", i
1526 END IF
1527 CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, stride=stride, mpi_io=mpi_io)
1528 CALL cp_print_key_finished_output(unit_nr, logger, print_section, '', &
1529 ignore_should_output=.true., mpi_io=mpi_io)
1530 END DO
1531 END DO
1532
1533 CALL auxbas_pw_pool%give_back_pw(wf_g)
1534 CALL auxbas_pw_pool%give_back_pw(wf_r)
1535
1536 END SUBROUTINE print_nto_cubes
1537
1538END 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:234
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)
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, 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)
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_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.