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