Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
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_set, 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_velocity
64 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
65 : section_vals_type,&
66 : section_vals_val_get
67 : USE kahan_sum, ONLY: accurate_dot_product
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,&
75 : mp_para_env_type,&
76 : mp_request_type
77 : USE molden_utils, ONLY: write_mos_molden
78 : USE moments_utils, ONLY: get_reference_point
79 : USE parallel_gemm_api, ONLY: parallel_gemm
80 : USE particle_list_types, ONLY: particle_list_type
81 : USE particle_types, ONLY: particle_type
82 : USE physcon, ONLY: evolt
83 : USE pw_env_types, ONLY: pw_env_get,&
84 : pw_env_type
85 : USE pw_poisson_types, ONLY: pw_poisson_type
86 : USE pw_pool_types, ONLY: pw_pool_p_type,&
87 : pw_pool_type
88 : USE pw_types, ONLY: pw_c1d_gs_type,&
89 : pw_r3d_rs_type
90 : USE qs_collocate_density, ONLY: calculate_wavefunction
91 : USE qs_environment_types, ONLY: get_qs_env,&
92 : qs_environment_type
93 : USE qs_kind_types, ONLY: qs_kind_type
94 : USE qs_ks_types, ONLY: qs_ks_env_type
95 : USE qs_mo_types, ONLY: allocate_mo_set,&
96 : deallocate_mo_set,&
97 : get_mo_set,&
98 : init_mo_set,&
99 : mo_set_type,&
100 : set_mo_set
101 : USE qs_moments, ONLY: build_berry_moment_matrix
102 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
103 : USE qs_operators_ao, ONLY: rRc_xyz_ao
104 : USE qs_overlap, ONLY: build_overlap_matrix
105 : USE qs_subsys_types, ONLY: qs_subsys_get,&
106 : qs_subsys_type
107 : USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos
108 : USE string_utilities, ONLY: integer_to_string
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 :
122 : PUBLIC :: tddfpt_dipole_operator, tddfpt_print_summary, tddfpt_print_excitation_analysis, &
123 : tddfpt_print_nto_analysis, tddfpt_print_exciton_descriptors
124 :
125 : ! **************************************************************************************************
126 :
127 : CONTAINS
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 1116 : 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 1116 : 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 1116 : 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 1116 : 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 1116 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: S_mos_virt
178 1116 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: dBerry_mos_occ, gamma_real_imag, opvec
179 1116 : 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 1116 : 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 1116 : CALL timeset(routineN, handle)
188 :
189 1116 : NULLIFY (blacs_env, cell, matrix_s, pw_env)
190 1116 : CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, matrix_s=matrix_s, pw_env=pw_env)
191 :
192 1116 : nspins = SIZE(gs_mos)
193 1116 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
194 2384 : DO ispin = 1, nspins
195 1268 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
196 2384 : nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
197 : END DO
198 :
199 : ! +++ allocate dipole operator matrices (must be deallocated elsewhere)
200 8420 : ALLOCATE (dipole_op_mos_occ(nderivs, nspins))
201 2384 : DO ispin = 1, nspins
202 1268 : CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct)
203 :
204 6188 : DO ideriv = 1, nderivs
205 5072 : CALL cp_fm_create(dipole_op_mos_occ(ideriv, ispin), fm_struct)
206 : END DO
207 : END DO
208 :
209 : ! +++ allocate work matrices
210 4616 : ALLOCATE (S_mos_virt(nspins))
211 2384 : DO ispin = 1, nspins
212 1268 : CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct)
213 1268 : 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 2384 : 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 1116 : CALL pw_env_get(pw_env, poisson_env=poisson_env)
222 4464 : ndim_periodic = COUNT(poisson_env%parameters%periodic == 1)
223 :
224 : ! select default for dipole form
225 1116 : IF (tddfpt_control%dipole_form == 0) THEN
226 524 : CALL get_qs_env(qs_env, dft_control=dft_control)
227 524 : IF (dft_control%qs_control%xtb) THEN
228 32 : IF (ndim_periodic == 0) THEN
229 0 : tddfpt_control%dipole_form = tddfpt_dipole_length
230 : ELSE
231 32 : tddfpt_control%dipole_form = tddfpt_dipole_berry
232 : END IF
233 : ELSE
234 492 : tddfpt_control%dipole_form = tddfpt_dipole_velocity
235 : END IF
236 : END IF
237 :
238 1316 : SELECT CASE (tddfpt_control%dipole_form)
239 : CASE (tddfpt_dipole_berry)
240 200 : IF (ndim_periodic /= 3) THEN
241 : CALL cp_warn(__LOCATION__, &
242 : "Fully periodic Poisson solver (PERIODIC xyz) is needed "// &
243 0 : "for oscillator strengths based on the Berry phase formula")
244 : END IF
245 :
246 200 : 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 200 : CALL dbcsr_allocate_matrix_set(berry_cossin_xyz, 2)
251 :
252 600 : DO i_cos_sin = 1, 2
253 400 : CALL dbcsr_init_p(berry_cossin_xyz(i_cos_sin)%matrix)
254 600 : 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 3000 : ALLOCATE (gamma_00(nspins), gamma_inv_00(nspins), gamma_real_imag(2, nspins), opvec(2, nspins))
259 1200 : ALLOCATE (dBerry_mos_occ(nderivs, nspins))
260 400 : DO ispin = 1, nspins
261 200 : NULLIFY (fm_struct)
262 : CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_occ(ispin), &
263 200 : ncol_global=nmo_occ(ispin), context=blacs_env)
264 :
265 200 : CALL cp_cfm_create(gamma_00(ispin), fm_struct)
266 200 : CALL cp_cfm_create(gamma_inv_00(ispin), fm_struct)
267 :
268 600 : DO i_cos_sin = 1, 2
269 600 : CALL cp_fm_create(gamma_real_imag(i_cos_sin, ispin), fm_struct)
270 : END DO
271 200 : CALL cp_fm_struct_release(fm_struct)
272 :
273 : ! G_real C_0, G_imag C_0
274 200 : CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct)
275 600 : DO i_cos_sin = 1, 2
276 600 : CALL cp_fm_create(opvec(i_cos_sin, ispin), fm_struct)
277 : END DO
278 :
279 : ! dBerry * C_0
280 1000 : DO ideriv = 1, nderivs
281 600 : CALL cp_fm_create(dBerry_mos_occ(ideriv, ispin), fm_struct)
282 800 : CALL cp_fm_set_all(dBerry_mos_occ(ideriv, ispin), 0.0_dp)
283 : END DO
284 : END DO
285 :
286 800 : DO ideriv = 1, nderivs
287 2400 : kvec(:) = twopi*cell%h_inv(ideriv, :)
288 1800 : DO i_cos_sin = 1, 2
289 1800 : 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 600 : berry_cossin_xyz(2)%matrix, kvec)
293 :
294 1400 : 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 1800 : 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 1800 : 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 600 : 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 600 : 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 600 : mtarget=gamma_00(ispin))
315 :
316 : ! gamma_inv_00 = Q = [C_0^T (gamma_real - i gamma_imag) C_0] ^ {-1}
317 600 : CALL cp_cfm_set_all(gamma_inv_00(ispin), z_zero, z_one)
318 600 : 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 600 : 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 600 : 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 600 : 1.0_dp, dipole_op_mos_occ(1, ispin))
331 :
332 3000 : DO jderiv = 1, nderivs
333 : CALL cp_fm_scale_and_add(1.0_dp, dBerry_mos_occ(jderiv, ispin), &
334 2400 : 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 200 : CALL cp_fm_release(opvec)
341 200 : CALL cp_fm_release(gamma_real_imag)
342 400 : DO ispin = nspins, 1, -1
343 200 : CALL cp_cfm_release(gamma_inv_00(ispin))
344 400 : CALL cp_cfm_release(gamma_00(ispin))
345 : END DO
346 200 : DEALLOCATE (gamma_00, gamma_inv_00)
347 200 : CALL dbcsr_deallocate_matrix_set(berry_cossin_xyz)
348 :
349 200 : NULLIFY (fm_struct)
350 200 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
351 200 : CALL cp_fm_create(wfm_ao_ao, fm_struct)
352 200 : 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 400 : 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 200 : 0.0_dp, wfm_ao_ao)
367 :
368 1000 : 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 800 : 0.0_dp, dipole_op_mos_occ(ideriv, ispin))
372 : END DO
373 : END DO
374 :
375 200 : CALL cp_fm_release(wfm_ao_ao)
376 400 : CALL cp_fm_release(dBerry_mos_occ)
377 :
378 : CASE (tddfpt_dipole_length)
379 8 : IF (ndim_periodic /= 0) THEN
380 : CALL cp_warn(__LOCATION__, &
381 : "Non-periodic Poisson solver (PERIODIC none) is needed "// &
382 0 : "for oscillator strengths based on the length operator")
383 : END IF
384 :
385 : ! compute components of the dipole operator in the length form
386 8 : NULLIFY (rRc_xyz)
387 8 : CALL dbcsr_allocate_matrix_set(rRc_xyz, nderivs)
388 :
389 32 : DO ideriv = 1, nderivs
390 24 : CALL dbcsr_init_p(rRc_xyz(ideriv)%matrix)
391 32 : 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 8 : 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 8 : minimum_image=.FALSE., soft=.FALSE.)
400 :
401 8 : NULLIFY (fm_struct)
402 8 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
403 8 : CALL cp_fm_create(wfm_ao_ao, fm_struct)
404 8 : CALL cp_fm_struct_release(fm_struct)
405 :
406 16 : DO ispin = 1, nspins
407 8 : CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct)
408 8 : 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 8 : 0.0_dp, wfm_ao_ao)
414 :
415 32 : 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 24 : 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 32 : 0.0_dp, dipole_op_mos_occ(ideriv, ispin))
424 : END DO
425 :
426 24 : CALL cp_fm_release(rRc_mos_occ)
427 : END DO
428 :
429 8 : CALL cp_fm_release(wfm_ao_ao)
430 8 : CALL dbcsr_deallocate_matrix_set(rRc_xyz)
431 :
432 : CASE (tddfpt_dipole_velocity)
433 : ! generate overlap derivatives
434 908 : CALL get_qs_env(qs_env, ks_env=ks_env, sab_orb=sab_orb)
435 908 : NULLIFY (scrm)
436 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, nderivative=1, &
437 : basis_type_a="ORB", basis_type_b="ORB", &
438 908 : sab_nl=sab_orb)
439 :
440 1968 : DO ispin = 1, nspins
441 1060 : NULLIFY (fm_struct)
442 : CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), &
443 1060 : ncol_global=nmo_occ(ispin), context=blacs_env)
444 1060 : CALL cp_fm_create(ediff_inv, fm_struct)
445 1060 : CALL cp_fm_create(wfm_mo_virt_mo_occ, fm_struct)
446 1060 : 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 1060 : row_indices=row_indices, col_indices=col_indices, local_data=local_data_ediff)
450 1060 : 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 1060 : !$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 4240 : 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 3180 : 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 3180 : 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 3180 : !$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 4240 : 0.0_dp, dipole_op_mos_occ(ideriv, ispin))
494 : END DO
495 :
496 1060 : CALL cp_fm_release(wfm_mo_virt_mo_occ)
497 4088 : CALL cp_fm_release(ediff_inv)
498 : END DO
499 908 : CALL dbcsr_deallocate_matrix_set(scrm)
500 :
501 : CASE DEFAULT
502 1124 : CPABORT("Unimplemented form of the dipole operator")
503 : END SELECT
504 :
505 : ! --- release work matrices
506 1116 : CALL cp_fm_release(S_mos_virt)
507 :
508 1116 : CALL timestop(handle)
509 3348 : 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 2248 : SUBROUTINE tddfpt_print_summary(log_unit, evects, evals, gs_mos, ostrength, mult, &
537 1124 : 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 1124 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: trans_dipoles
556 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
557 : TYPE(cp_fm_type) :: dipact
558 :
559 1124 : CALL timeset(routineN, handle)
560 :
561 1124 : nspins = SIZE(evects, 1)
562 1124 : nstates = SIZE(evects, 2)
563 :
564 1124 : IF (nspins > 1) THEN
565 130 : lsd_str = 'U'
566 : ELSE
567 994 : lsd_str = 'R'
568 : END IF
569 :
570 : ! *** summary header ***
571 1124 : IF (log_unit > 0) THEN
572 562 : CALL integer_to_string(mult, mult_str)
573 562 : WRITE (log_unit, '(/,1X,A1,A,1X,A)') lsd_str, "-TDDFPT states of multiplicity", TRIM(mult_str)
574 662 : SELECT CASE (dipole_form)
575 : CASE (tddfpt_dipole_berry)
576 100 : WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using Berry operator formulation"
577 : CASE (tddfpt_dipole_length)
578 7 : WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using length formulation"
579 : CASE (tddfpt_dipole_velocity)
580 455 : WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using velocity formulation"
581 : CASE DEFAULT
582 562 : CPABORT("Unimplemented form of the dipole operator")
583 : END SELECT
584 :
585 562 : WRITE (log_unit, '(T10,A,T19,A,T37,A,T69,A)') "State", "Excitation", &
586 1124 : "Transition dipole (a.u.)", "Oscillator"
587 562 : WRITE (log_unit, '(T10,A,T19,A,T37,A,T49,A,T61,A,T67,A)') "number", "energy (eV)", &
588 1124 : "x", "y", "z", "strength (a.u.)"
589 562 : WRITE (log_unit, '(T10,72("-"))')
590 : END IF
591 :
592 : ! transition dipole moment
593 4496 : ALLOCATE (trans_dipoles(nstates, nderivs, nspins))
594 16574 : trans_dipoles(:, :, :) = 0.0_dp
595 :
596 : ! nspins == 1 .AND. mult == 3 : spin-flip transitions are forbidden due to symmetry reasons
597 1124 : IF (nspins > 1 .OR. mult == 1) THEN
598 1982 : DO ispin = 1, nspins
599 1056 : CALL cp_fm_get_info(dipole_op_mos_occ(1, ispin), nrow_global=nao, ncol_global=nocc)
600 1056 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive)
601 1982 : IF (nocc == nactive) THEN
602 4176 : DO ideriv = 1, nderivs
603 : CALL cp_fm_trace(evects(ispin, :), dipole_op_mos_occ(ideriv, ispin), &
604 4176 : trans_dipoles(:, ideriv, ispin))
605 : END DO
606 : ELSE
607 12 : CALL cp_fm_get_info(evects(ispin, 1), matrix_struct=matrix_struct)
608 12 : CALL cp_fm_create(dipact, matrix_struct)
609 48 : DO ideriv = 1, nderivs
610 90 : DO i = 1, nactive
611 54 : j = gs_mos(ispin)%index_active(i)
612 : CALL cp_fm_to_fm(dipole_op_mos_occ(ideriv, ispin), dipact, &
613 90 : ncol=1, source_start=j, target_start=i)
614 : END DO
615 48 : CALL cp_fm_trace(evects(ispin, :), dipact, trans_dipoles(:, ideriv, ispin))
616 : END DO
617 12 : CALL cp_fm_release(dipact)
618 : END IF
619 : END DO
620 :
621 926 : IF (nspins == 1) THEN
622 9094 : trans_dipoles(:, :, 1) = SQRT(2.0_dp)*trans_dipoles(:, :, 1)
623 : ELSE
624 1732 : trans_dipoles(:, :, 1) = SQRT(trans_dipoles(:, :, 1)**2 + trans_dipoles(:, :, 2)**2)
625 : END IF
626 : END IF
627 :
628 : ! *** summary information ***
629 4198 : DO istate = 1, nstates
630 : osc_strength = 2.0_dp/3.0_dp*evals(istate)* &
631 3074 : accurate_dot_product(trans_dipoles(istate, :, 1), trans_dipoles(istate, :, 1))
632 3074 : ostrength(istate) = osc_strength
633 4198 : 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 1537 : "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 1124 : IF (log_unit > 0) THEN
641 2099 : WRITE (log_unit, '(/,T2,A,E14.6)') 'TDDFPT : CheckSum =', SQRT(SUM(evals**2))
642 : END IF
643 :
644 1124 : DEALLOCATE (trans_dipoles)
645 :
646 1124 : CALL timestop(handle)
647 1124 : 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 1124 : 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 1124 : 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 1124 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
688 1124 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
689 : INTEGER, DIMENSION(maxspins) :: nactive, nmo_occ, nmo_virt
690 : LOGICAL :: do_exc_analysis
691 1124 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights_local, weights_neg_abs_recv, &
692 1124 : weights_recv
693 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
694 1124 : POINTER :: local_data
695 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
696 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
697 1124 : 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 1124 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_handlers, recv_handlers2
701 :
702 1124 : CALL timeset(routineN, handle)
703 :
704 1124 : nspins = SIZE(gs_mos, 1)
705 1124 : nstates = SIZE(evects, 2)
706 1124 : do_exc_analysis = min_amplitude < 1.0_dp
707 :
708 1124 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env, para_env=para_env)
709 1124 : CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
710 :
711 2400 : DO ispin = 1, nspins
712 1276 : nactive(ispin) = gs_mos(ispin)%nmo_active
713 : nactive8(ispin) = INT(nactive(ispin), kind=int_8)
714 1276 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
715 1276 : nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
716 1276 : nmo_occ8(ispin) = SIZE(gs_mos(ispin)%evals_occ, kind=int_8)
717 2400 : nmo_virt8(ispin) = SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
718 : END DO
719 :
720 : ! *** excitation analysis ***
721 1124 : IF (do_exc_analysis) THEN
722 1124 : CPASSERT(log_unit <= 0 .OR. para_env%is_source())
723 1124 : nmo_virt_occ_alpha = INT(nmo_virt(1), int_8)*INT(nmo_occ(1), int_8)
724 :
725 1124 : IF (log_unit > 0) THEN
726 562 : WRITE (log_unit, "(1X,A)") "", &
727 562 : "-------------------------------------------------------------------------------", &
728 562 : "- Excitation analysis -", &
729 1124 : "-------------------------------------------------------------------------------"
730 562 : WRITE (log_unit, '(8X,A,T27,A,T49,A,T69,A)') "State", "Occupied", "Virtual", "Excitation"
731 562 : WRITE (log_unit, '(8X,A,T28,A,T49,A,T69,A)') "number", "orbital", "orbital", "amplitude"
732 562 : WRITE (log_unit, '(1X,79("-"))')
733 :
734 562 : IF (nspins == 1) THEN
735 486 : state_spin = 1
736 486 : state_spin2 = 2
737 486 : spin_label = ' '
738 486 : spin_label2 = ' '
739 76 : ELSE IF (spinflip /= no_sf_tddfpt) THEN
740 11 : state_spin = 1
741 11 : state_spin2 = 2
742 11 : spin_label = '(alp)'
743 11 : spin_label2 = '(bet)'
744 : END IF
745 : END IF
746 :
747 7004 : ALLOCATE (S_mos_virt(SIZE(evects, 1)), weights_fm(SIZE(evects, 1)))
748 2378 : DO ispin = 1, SIZE(evects, 1)
749 1254 : IF (spinflip == no_sf_tddfpt) THEN
750 : spin2 = ispin
751 : ELSE
752 22 : spin2 = 2
753 : END IF
754 1254 : CALL cp_fm_get_info(gs_mos(spin2)%mos_virt, matrix_struct=fm_struct)
755 1254 : 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 1254 : ncol=nmo_virt(spin2), alpha=1.0_dp, beta=0.0_dp)
760 :
761 1254 : NULLIFY (fm_struct)
762 : CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(spin2), ncol_global=nactive(ispin), &
763 1254 : context=blacs_env)
764 1254 : CALL cp_fm_create(weights_fm(ispin), fm_struct)
765 1254 : CALL cp_fm_set_all(weights_fm(ispin), 0.0_dp)
766 2378 : CALL cp_fm_struct_release(fm_struct)
767 : END DO
768 :
769 2378 : nexcs_max_local = 0
770 2378 : DO ispin = 1, SIZE(evects, 1)
771 1254 : CALL cp_fm_get_info(weights_fm(ispin), nrow_local=nrows_local, ncol_local=ncols_local)
772 2378 : nexcs_max_local = nexcs_max_local + INT(nrows_local, int_8)*INT(ncols_local, int_8)
773 : END DO
774 :
775 4496 : ALLOCATE (weights_local(nexcs_max_local), inds_local(nexcs_max_local))
776 :
777 4198 : DO istate = 1, nstates
778 6552 : nexcs_local = 0
779 6552 : nmo_virt_occ = 0
780 :
781 : ! analyse matrix elements locally and transfer only significant
782 : ! excitations to the master node for subsequent ordering
783 6552 : DO ispin = 1, SIZE(evects, 1)
784 3478 : IF (spinflip == no_sf_tddfpt) THEN
785 : spin2 = ispin
786 : ELSE
787 90 : 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 3478 : 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 3478 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
795 :
796 : ! locate single excitations with significant amplitudes (>= min_amplitude)
797 19250 : DO icol = 1, ncols_local
798 204900 : DO irow = 1, nrows_local
799 201422 : IF (ABS(local_data(irow, icol)) >= min_amplitude) THEN
800 : ! number of non-negligible excitations
801 2441 : nexcs_local = nexcs_local + 1
802 : ! excitation amplitude
803 2441 : 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 2441 : INT(col_indices(icol) - 1, int_8)*nmo_virt8(spin2)
807 : END IF
808 : END DO
809 : END DO
810 :
811 10030 : nmo_virt_occ = nmo_virt_occ + nmo_virt8(spin2)*nmo_occ8(ispin)
812 : END DO
813 :
814 3074 : IF (para_env%is_source()) THEN
815 : ! master node
816 15370 : 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 4611 : DO iproc = 1, para_env%num_pe
820 4611 : IF (iproc - 1 /= para_env%mepos) THEN
821 1537 : CALL para_env%irecv(nexcs_recv(iproc:iproc), iproc - 1, recv_handlers(iproc), 0)
822 : ELSE
823 1537 : nexcs_recv(iproc) = nexcs_local
824 : END IF
825 : END DO
826 :
827 4611 : DO iproc = 1, para_env%num_pe
828 3074 : IF (iproc - 1 /= para_env%mepos) &
829 3074 : CALL recv_handlers(iproc)%wait()
830 : END DO
831 :
832 : ! compute total number of non-negligible excitations
833 1537 : nexcs = 0
834 4611 : DO iproc = 1, para_env%num_pe
835 4611 : nexcs = nexcs + nexcs_recv(iproc)
836 : END DO
837 :
838 : ! receive indices and amplitudes of selected excitations
839 6148 : ALLOCATE (weights_recv(nexcs), weights_neg_abs_recv(nexcs))
840 6148 : ALLOCATE (inds_recv(nexcs), inds(nexcs))
841 :
842 4611 : nmo_virt_occ = 0
843 4611 : DO iproc = 1, para_env%num_pe
844 4611 : IF (nexcs_recv(iproc) > 0) THEN
845 1617 : 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 236 : 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 236 : iproc - 1, recv_handlers2(iproc), 2)
852 : ELSE
853 : ! data on master node
854 3549 : weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = weights_local(1:nexcs_recv(iproc))
855 3549 : inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = inds_local(1:nexcs_recv(iproc))
856 : END IF
857 :
858 1617 : nmo_virt_occ = nmo_virt_occ + nexcs_recv(iproc)
859 : END IF
860 : END DO
861 :
862 4611 : DO iproc = 1, para_env%num_pe
863 4611 : IF (iproc - 1 /= para_env%mepos .AND. nexcs_recv(iproc) > 0) THEN
864 236 : CALL recv_handlers(iproc)%wait()
865 236 : CALL recv_handlers2(iproc)%wait()
866 : END IF
867 : END DO
868 :
869 1537 : DEALLOCATE (nexcs_recv, recv_handlers, recv_handlers2)
870 : ELSE
871 : ! working node: send the number of selected excited states to the master node
872 1537 : nexcs_send(1) = nexcs_local
873 1537 : CALL para_env%isend(nexcs_send, para_env%source, send_handler, 0)
874 1537 : CALL send_handler%wait()
875 :
876 1537 : IF (nexcs_local > 0) THEN
877 : ! send excitation amplitudes
878 236 : CALL para_env%isend(weights_local(1:nexcs_local), para_env%source, send_handler, 1)
879 : ! send compressed indices
880 236 : CALL para_env%isend(inds_local(1:nexcs_local), para_env%source, send_handler2, 2)
881 :
882 236 : CALL send_handler%wait()
883 236 : 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 3074 : IF (para_env%is_source() .AND. log_unit > 0) THEN
890 3978 : weights_neg_abs_recv(:) = -ABS(weights_recv)
891 1537 : CALL sort(weights_neg_abs_recv, INT(nexcs), inds)
892 :
893 1537 : 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 1537 : state_spin = 1
898 1537 : state_spin2 = 1
899 1537 : IF (spinflip /= no_sf_tddfpt) THEN
900 45 : state_spin = 1
901 45 : state_spin2 = 2
902 : END IF
903 3978 : DO iexc = 1, nexcs
904 2441 : ind = inds_recv(inds(iexc)) - 1
905 2441 : IF ((nspins > 1) .AND. (spinflip == no_sf_tddfpt)) THEN
906 443 : IF (ind < nmo_virt_occ_alpha) THEN
907 200 : spin_label = '(alp)'
908 200 : spin_label2 = '(alp)'
909 : ELSE
910 243 : state_spin = 2
911 243 : state_spin2 = 2
912 243 : ind = ind - nmo_virt_occ_alpha
913 243 : spin_label = '(bet)'
914 243 : spin_label2 = '(bet)'
915 : END IF
916 : END IF
917 2441 : imo_act = ind/nmo_virt8(state_spin2) + 1
918 2441 : imo_occ = gs_mos(state_spin)%index_active(imo_act)
919 2441 : imo_virt = MOD(ind, nmo_virt8(state_spin2)) + 1
920 :
921 2441 : WRITE (log_unit, '(T27,I8,1X,A5,T48,I8,1X,A5,T70,F9.6)') imo_occ, spin_label, &
922 6419 : nmo_occ8(state_spin2) + imo_virt, spin_label2, weights_recv(inds(iexc))
923 : END DO
924 : END IF
925 :
926 : ! deallocate temporary arrays
927 3074 : IF (para_env%is_source()) &
928 2661 : DEALLOCATE (weights_recv, weights_neg_abs_recv, inds_recv, inds)
929 : END DO
930 :
931 1124 : DEALLOCATE (weights_local, inds_local)
932 1124 : IF (log_unit > 0) THEN
933 : WRITE (log_unit, "(1X,A)") &
934 562 : "-------------------------------------------------------------------------------"
935 : END IF
936 : END IF
937 :
938 1124 : CALL cp_fm_release(weights_fm)
939 1124 : CALL cp_fm_release(S_mos_virt)
940 :
941 1124 : CALL timestop(handle)
942 :
943 2248 : END SUBROUTINE tddfpt_print_excitation_analysis
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 1124 : 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 1124 : INTEGER, DIMENSION(:), POINTER :: slist, stride
977 : LOGICAL :: append_cube, cube_file, explicit
978 : REAL(KIND=dp) :: os_threshold, sume, threshold
979 1124 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigvals
980 1124 : 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 1124 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: teig
985 : TYPE(cp_logger_type), POINTER :: logger
986 1124 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: nto_set
987 1124 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
988 1124 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
989 : TYPE(section_vals_type), POINTER :: molden_section, nto_section
990 :
991 1124 : CALL timeset(routineN, handle)
992 :
993 1124 : logger => cp_get_default_logger()
994 1124 : iounit = cp_logger_get_default_io_unit(logger)
995 :
996 1124 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_section, &
997 : "NTO_ANALYSIS"), cp_p_file)) THEN
998 :
999 144 : CALL cite_reference(Martin2003)
1000 :
1001 144 : CALL section_vals_val_get(print_section, "NTO_ANALYSIS%THRESHOLD", r_val=threshold)
1002 144 : CALL section_vals_val_get(print_section, "NTO_ANALYSIS%INTENSITY_THRESHOLD", r_val=os_threshold)
1003 144 : CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", EXPLICIT=explicit)
1004 :
1005 144 : IF (explicit) THEN
1006 4 : CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", i_vals=slist)
1007 4 : nlist = SIZE(slist)
1008 : ELSE
1009 : nlist = 0
1010 : END IF
1011 :
1012 144 : IF (iounit > 0) THEN
1013 72 : WRITE (iounit, "(1X,A)") "", &
1014 72 : "-------------------------------------------------------------------------------", &
1015 72 : "- Natural Orbital analysis -", &
1016 144 : "-------------------------------------------------------------------------------"
1017 : END IF
1018 :
1019 144 : nspins = SIZE(evects, 1)
1020 144 : nstates = SIZE(evects, 2)
1021 144 : CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1022 :
1023 344 : DO istate = 1, nstates
1024 200 : IF (os_threshold > ostrength(istate)) THEN
1025 44 : IF (iounit > 0) THEN
1026 22 : WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1027 : END IF
1028 : CYCLE
1029 : END IF
1030 156 : IF (nlist > 0) THEN
1031 12 : IF (.NOT. ANY(slist == istate)) THEN
1032 0 : IF (iounit > 0) THEN
1033 0 : WRITE (iounit, "(1X,A,I6)") " Skipping state ", istate
1034 : END IF
1035 : CYCLE
1036 : END IF
1037 : END IF
1038 156 : IF (iounit > 0) THEN
1039 78 : WRITE (iounit, "(1X,A,I6,T30,F10.5,A)") " STATE NR. ", istate, evals(istate)*evolt, " eV"
1040 : END IF
1041 : nmax = 0
1042 326 : DO ispin = 1, nspins
1043 170 : CALL cp_fm_get_info(evects(ispin, istate), matrix_struct=fm_struct, ncol_global=nmo)
1044 326 : nmax = MAX(nmax, nmo)
1045 : END DO
1046 624 : ALLOCATE (eigenvalues(nmax, nspins))
1047 1434 : eigenvalues = 0.0_dp
1048 : ! SET 1: Hole states
1049 : ! SET 2: Particle states
1050 156 : nto_name(1) = 'Hole_states'
1051 156 : nto_name(2) = 'Particle_states'
1052 468 : ALLOCATE (nto_set(2))
1053 468 : DO i = 1, 2
1054 312 : CALL allocate_mo_set(nto_set(i), nao, ntomax, 0, 0.0_dp, 1.0_dp, 0.0_dp)
1055 312 : 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 312 : ncol_global=ntomax)
1058 312 : CALL cp_fm_create(tmat, fm_mo_struct)
1059 312 : CALL init_mo_set(nto_set(i), fm_ref=tmat, name=nto_name(i))
1060 312 : CALL cp_fm_release(tmat)
1061 780 : CALL cp_fm_struct_release(fm_mo_struct)
1062 : END DO
1063 : !
1064 638 : ALLOCATE (teig(nspins))
1065 : ! hole states
1066 : ! Diagonalize X(T)*S*X
1067 326 : DO ispin = 1, nspins
1068 : ASSOCIATE (ev => evects(ispin, istate))
1069 170 : CALL cp_fm_get_info(ev, matrix_struct=fm_struct, ncol_global=nmo)
1070 170 : CALL cp_fm_create(Sev, fm_struct)
1071 : CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
1072 170 : nrow_global=nmo, ncol_global=nmo)
1073 170 : CALL cp_fm_create(tmat, fm_mo_struct)
1074 170 : CALL cp_fm_create(teig(ispin), fm_mo_struct)
1075 170 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, Sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1076 170 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, ev, Sev, 0.0_dp, tmat)
1077 : END ASSOCIATE
1078 :
1079 170 : CALL choose_eigv_solver(tmat, teig(ispin), eigenvalues(1:nmo, ispin))
1080 :
1081 170 : CALL cp_fm_struct_release(fm_mo_struct)
1082 170 : CALL cp_fm_release(tmat)
1083 666 : CALL cp_fm_release(Sev)
1084 : END DO
1085 : ! find major determinants i->a
1086 156 : ia_index = 0
1087 156 : sume = 0.0_dp
1088 156 : nnto = 0
1089 186 : DO i = 1, ntomax
1090 186 : iv = MAXLOC(eigenvalues)
1091 186 : ia_eval(i) = eigenvalues(iv(1), iv(2))
1092 558 : ia_index(1:2, i) = iv(1:2)
1093 186 : sume = sume + ia_eval(i)
1094 186 : eigenvalues(iv(1), iv(2)) = 0.0_dp
1095 186 : nnto = nnto + 1
1096 186 : IF (sume > threshold) EXIT
1097 : END DO
1098 : ! store hole states
1099 156 : CALL set_mo_set(nto_set(1), nmo=nnto)
1100 342 : DO i = 1, nnto
1101 186 : ia = ia_index(1, i)
1102 186 : ispin = ia_index(2, i)
1103 186 : CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=nmo)
1104 186 : 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 186 : nrow_global=nmo, ncol_global=1)
1107 186 : CALL cp_fm_create(tmat, fm_mo_struct)
1108 186 : CALL cp_fm_struct_release(fm_mo_struct)
1109 186 : 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 186 : ncol_global=1)
1112 186 : CALL cp_fm_create(wvec, fm_mo_struct)
1113 186 : CALL cp_fm_struct_release(fm_mo_struct)
1114 186 : 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 186 : tmat, 0.0_dp, wvec)
1117 186 : CALL cp_fm_to_fm(wvec, nto_set(1)%mo_coeff, 1, 1, i)
1118 186 : CALL cp_fm_release(wvec)
1119 900 : 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 156 : CALL set_mo_set(nto_set(2), nmo=nnto)
1124 326 : DO ispin = 1, nspins
1125 : ASSOCIATE (ev => evects(ispin, istate))
1126 170 : CALL cp_fm_get_info(ev, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
1127 510 : ALLOCATE (eigvals(nao))
1128 4534 : eigvals = 0.0_dp
1129 170 : CALL cp_fm_create(Sev, fm_struct)
1130 340 : 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 170 : nrow_global=nao, ncol_global=nao)
1134 170 : CALL cp_fm_create(tmat, fm_mo_struct)
1135 170 : CALL cp_fm_create(smat, fm_mo_struct)
1136 170 : CALL cp_fm_create(wmat, fm_mo_struct)
1137 170 : CALL cp_fm_create(work, fm_mo_struct)
1138 170 : CALL cp_fm_struct_release(fm_mo_struct)
1139 170 : CALL copy_dbcsr_to_fm(matrix_s, smat)
1140 170 : CALL parallel_gemm('N', 'T', nao, nao, nmo, 1.0_dp, Sev, Sev, 0.0_dp, tmat)
1141 170 : CALL cp_fm_geeig(tmat, smat, wmat, eigvals, work)
1142 376 : DO i = 1, nnto
1143 376 : IF (ispin == ia_index(2, i)) THEN
1144 186 : icg = 0
1145 4586 : DO j = 1, nao
1146 4586 : IF (ABS(eigvals(j) - ia_eval(i)) < 1.E-6_dp) THEN
1147 186 : icg = j
1148 186 : EXIT
1149 : END IF
1150 : END DO
1151 186 : IF (icg == 0) THEN
1152 : CALL cp_warn(__LOCATION__, &
1153 0 : "Could not locate particle state associated with hole state.")
1154 : ELSE
1155 186 : CALL cp_fm_to_fm(wmat, nto_set(2)%mo_coeff, 1, icg, i)
1156 : END IF
1157 : END IF
1158 : END DO
1159 170 : DEALLOCATE (eigvals)
1160 170 : CALL cp_fm_release(Sev)
1161 170 : CALL cp_fm_release(tmat)
1162 170 : CALL cp_fm_release(smat)
1163 170 : CALL cp_fm_release(wmat)
1164 496 : CALL cp_fm_release(work)
1165 : END DO
1166 : ! print
1167 156 : IF (iounit > 0) THEN
1168 78 : sume = 0.0_dp
1169 171 : DO i = 1, nnto
1170 93 : sume = sume + ia_eval(i)
1171 : WRITE (iounit, "(T6,A,i2,T30,A,i1,T42,A,F8.5,T63,A,F8.5)") &
1172 93 : "Particle-Hole state:", i, " Spin:", ia_index(2, i), &
1173 264 : "Eigenvalue:", ia_eval(i), " Sum Eigv:", sume
1174 : END DO
1175 : END IF
1176 : ! Cube and Molden files
1177 156 : nto_section => section_vals_get_subs_vals(print_section, "NTO_ANALYSIS")
1178 156 : CALL section_vals_val_get(nto_section, "CUBE_FILES", l_val=cube_file)
1179 156 : CALL section_vals_val_get(nto_section, "STRIDE", i_vals=stride)
1180 156 : CALL section_vals_val_get(nto_section, "APPEND", l_val=append_cube)
1181 156 : IF (cube_file) THEN
1182 16 : CALL print_nto_cubes(qs_env, nto_set, istate, stride, append_cube, nto_section)
1183 : END IF
1184 156 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set)
1185 156 : molden_section => section_vals_get_subs_vals(print_section, "MOS_MOLDEN")
1186 156 : CALL write_mos_molden(nto_set, qs_kind_set, particle_set, molden_section)
1187 : !
1188 156 : DEALLOCATE (eigenvalues)
1189 156 : CALL cp_fm_release(teig)
1190 : !
1191 468 : DO i = 1, 2
1192 468 : CALL deallocate_mo_set(nto_set(i))
1193 : END DO
1194 612 : DEALLOCATE (nto_set)
1195 : END DO
1196 :
1197 144 : IF (iounit > 0) THEN
1198 : WRITE (iounit, "(1X,A)") &
1199 72 : "-------------------------------------------------------------------------------"
1200 : END IF
1201 :
1202 : END IF
1203 :
1204 1124 : CALL timestop(handle)
1205 :
1206 2248 : 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 2 : 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 2 : 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 2 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: eigvec_X_ia_n, fm_multipole_ab, &
1242 2 : fm_multipole_ai, fm_multipole_ij, &
1243 2 : S_mos_virt
1244 2 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_coeff
1245 : TYPE(exciton_descr_type), ALLOCATABLE, &
1246 2 : DIMENSION(:) :: exc_descr
1247 :
1248 2 : CALL timeset(routineN, handle)
1249 :
1250 2 : nspins = SIZE(evects, 1)
1251 2 : nstates = SIZE(evects, 2)
1252 :
1253 2 : CPASSERT(nspins == 1) ! Other spins are not yet implemented for exciton descriptors
1254 :
1255 2 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env)
1256 2 : CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
1257 :
1258 4 : DO ispin = 1, nspins
1259 2 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
1260 4 : 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 8 : ALLOCATE (mo_coeff(nspins))
1265 : CALL cp_fm_struct_create(fm_struct_mo_coeff, nrow_global=nao, ncol_global=nao, &
1266 2 : context=blacs_env)
1267 4 : DO ispin = 1, nspins
1268 2 : 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 2 : 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 4 : blacs_env)
1287 : END DO
1288 2 : 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 2 : n_moments_quad = 9
1293 2 : ALLOCATE (ref_point_multipole(3))
1294 20 : ALLOCATE (fm_multipole_ij(n_moments_quad))
1295 20 : ALLOCATE (fm_multipole_ab(n_moments_quad))
1296 20 : 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 2 : nmo_occ(1), nmo_virt(1), blacs_env)
1301 :
1302 2 : CALL cp_fm_release(mo_coeff)
1303 :
1304 : ! Compute eigenvector X of the Casida equation from trial vectors
1305 14 : ALLOCATE (S_mos_virt(nspins), eigvec_X_ia_n(nspins))
1306 4 : DO ispin = 1, nspins
1307 2 : CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct_S_mos_virt)
1308 2 : CALL cp_fm_create(S_mos_virt(ispin), fm_struct_S_mos_virt)
1309 2 : 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 2 : 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 2 : context=blacs_env)
1317 2 : CALL cp_fm_create(eigvec_X_ia_n(ispin), fm_struct_X_ia_n)
1318 4 : CALL cp_fm_struct_release(fm_struct_X_ia_n)
1319 : END DO
1320 172 : ALLOCATE (exc_descr(nstates))
1321 12 : DO istate = 1, nstates
1322 22 : DO ispin = 1, nspins
1323 10 : 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 10 : CALL cp_fm_get_info(evects(ispin, istate), ncol_global=nactive)
1331 10 : IF (nactive /= nmo_occ(ispin)) THEN
1332 : CALL cp_abort(__LOCATION__, &
1333 0 : "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 10 : 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 30 : istate, nmo_occ(ispin), nmo_virt(ispin))
1342 : END DO
1343 : END DO
1344 2 : CALL cp_fm_release(eigvec_X_ia_n)
1345 2 : CALL cp_fm_release(S_mos_virt)
1346 2 : CALL cp_fm_release(fm_multipole_ai)
1347 2 : CALL cp_fm_release(fm_multipole_ij)
1348 2 : CALL cp_fm_release(fm_multipole_ab)
1349 :
1350 : ! Actual printing
1351 2 : print_checkvalue = .TRUE.
1352 2 : prefix_output = ' '
1353 : CALL print_exciton_descriptors(exc_descr, ref_point_multipole, log_unit, &
1354 : nstates, print_checkvalue, do_directional_exciton_descriptors, &
1355 2 : prefix_output, qs_env)
1356 :
1357 2 : DEALLOCATE (ref_point_multipole)
1358 2 : DEALLOCATE (exc_descr)
1359 :
1360 2 : CALL timestop(handle)
1361 :
1362 6 : END SUBROUTINE tddfpt_print_exciton_descriptors
1363 :
1364 : ! **************************************************************************************************
1365 : !> \brief ...
1366 : !> \param vin ...
1367 : !> \param vout ...
1368 : !> \param mos_occ ...
1369 : !> \param matrix_s ...
1370 : ! **************************************************************************************************
1371 0 : 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 0 : CALL timeset(routineN, handle)
1384 :
1385 0 : 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 0 : nrow_global=nao, ncol_global=1)
1388 0 : CALL cp_fm_create(vec, fm_vec_struct)
1389 0 : CALL cp_fm_create(svec, fm_vec_struct)
1390 0 : CALL cp_fm_struct_release(fm_vec_struct)
1391 : CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
1392 0 : nrow_global=nmo, ncol_global=1)
1393 0 : CALL cp_fm_create(csvec, fm_vec_struct)
1394 0 : CALL cp_fm_struct_release(fm_vec_struct)
1395 :
1396 0 : CALL copy_dbcsr_to_fm(vin, vec)
1397 0 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, vec, svec, ncol=1, alpha=1.0_dp, beta=0.0_dp)
1398 0 : CALL parallel_gemm('T', 'N', nmo, 1, nao, 1.0_dp, mos_occ, svec, 0.0_dp, csvec)
1399 0 : CALL parallel_gemm('N', 'N', nao, 1, nmo, -1.0_dp, mos_occ, csvec, 1.0_dp, vec)
1400 0 : CALL cp_fm_vectorsnorm(vec, norm)
1401 0 : CPASSERT(norm(1) > 1.e-14_dp)
1402 0 : norm(1) = SQRT(1._dp/norm(1))
1403 0 : CALL cp_fm_scale(norm(1), vec)
1404 0 : CALL copy_fm_to_dbcsr(vec, vout, keep_sparsity=.FALSE.)
1405 :
1406 0 : CALL cp_fm_release(csvec)
1407 0 : CALL cp_fm_release(svec)
1408 0 : CALL cp_fm_release(vec)
1409 :
1410 0 : CALL timestop(handle)
1411 :
1412 0 : END SUBROUTINE project_vector
1413 :
1414 : ! **************************************************************************************************
1415 : !> \brief ...
1416 : !> \param va ...
1417 : !> \param vb ...
1418 : !> \param res ...
1419 : ! **************************************************************************************************
1420 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vba, vbb
1429 : TYPE(dbcsr_iterator_type) :: iter
1430 : TYPE(mp_comm_type) :: group
1431 :
1432 0 : CALL timeset(routineN, handle)
1433 :
1434 0 : res = 0.0_dp
1435 :
1436 0 : CALL dbcsr_get_info(va, group=group)
1437 0 : CALL dbcsr_iterator_start(iter, va)
1438 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1439 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, vba)
1440 0 : CALL dbcsr_get_block_p(vb, row=irow, col=icol, block=vbb, found=found)
1441 0 : res = res + SUM(vba*vbb)
1442 0 : CPASSERT(found)
1443 : END DO
1444 0 : CALL dbcsr_iterator_stop(iter)
1445 0 : CALL group%sum(res)
1446 :
1447 0 : CALL timestop(handle)
1448 :
1449 0 : 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 16 : 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 16 : 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 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1479 : TYPE(pw_c1d_gs_type) :: wf_g
1480 : TYPE(pw_env_type), POINTER :: pw_env
1481 16 : 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 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1485 : TYPE(qs_subsys_type), POINTER :: subsys
1486 :
1487 32 : logger => cp_get_default_logger()
1488 :
1489 16 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
1490 16 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1491 16 : CALL auxbas_pw_pool%create_pw(wf_r)
1492 16 : CALL auxbas_pw_pool%create_pw(wf_g)
1493 :
1494 16 : CALL get_qs_env(qs_env, subsys=subsys)
1495 16 : CALL qs_subsys_get(subsys, particles=particles)
1496 :
1497 16 : my_pos_cube = "REWIND"
1498 16 : IF (append_cube) THEN
1499 0 : 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 16 : particle_set=particle_set)
1507 :
1508 48 : DO iset = 1, 2
1509 32 : CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo)
1510 92 : DO i = 1, nmo
1511 : CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1512 44 : cell, dft_control, particle_set, pw_env)
1513 44 : IF (iset == 1) THEN
1514 22 : WRITE (filename, '(a4,I3.3,I2.2,a11)') "NTO_STATE", istate, i, "_Hole_State"
1515 22 : ELSEIF (iset == 2) THEN
1516 22 : WRITE (filename, '(a4,I3.3,I2.2,a15)') "NTO_STATE", istate, i, "_Particle_State"
1517 : END IF
1518 44 : 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 44 : log_filename=.FALSE., ignore_should_output=.TRUE., mpi_io=mpi_io)
1522 44 : IF (iset == 1) THEN
1523 22 : WRITE (title, *) "Natural Transition Orbital Hole State", i
1524 22 : ELSEIF (iset == 2) THEN
1525 22 : WRITE (title, *) "Natural Transition Orbital Particle State", i
1526 : END IF
1527 44 : 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 76 : ignore_should_output=.TRUE., mpi_io=mpi_io)
1530 : END DO
1531 : END DO
1532 :
1533 16 : CALL auxbas_pw_pool%give_back_pw(wf_g)
1534 16 : CALL auxbas_pw_pool%give_back_pw(wf_r)
1535 :
1536 16 : END SUBROUTINE print_nto_cubes
1537 :
1538 : END MODULE qs_tddfpt2_properties
|
