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