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