Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : MODULE qs_tddfpt2_utils
9 : USE cell_types, ONLY: cell_type
10 : USE cp_array_utils, ONLY: cp_1d_r_p_type
11 : USE cp_blacs_env, ONLY: cp_blacs_env_type
12 : USE cp_control_types, ONLY: dft_control_type,&
13 : tddfpt2_control_type
14 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
15 : dbcsr_copy,&
16 : dbcsr_get_info,&
17 : dbcsr_init_p,&
18 : dbcsr_p_type,&
19 : dbcsr_type
20 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
21 : cp_dbcsr_plus_fm_fm_t,&
22 : cp_dbcsr_sm_fm_multiply,&
23 : dbcsr_allocate_matrix_set
24 : USE cp_fm_basic_linalg, ONLY: cp_fm_triangular_invert
25 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
26 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
27 : fm_pool_create_fm
28 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
29 : cp_fm_struct_release,&
30 : cp_fm_struct_type
31 : USE cp_fm_types, ONLY: cp_fm_create,&
32 : cp_fm_get_info,&
33 : cp_fm_release,&
34 : cp_fm_set_all,&
35 : cp_fm_to_fm,&
36 : cp_fm_to_fm_submat,&
37 : cp_fm_type
38 : USE cp_log_handling, ONLY: cp_get_default_logger,&
39 : cp_logger_get_default_io_unit,&
40 : cp_logger_type
41 : USE exstates_types, ONLY: excited_energy_type
42 : USE input_constants, ONLY: &
43 : cholesky_dbcsr, cholesky_inverse, cholesky_off, cholesky_restore, no_sf_tddfpt, oe_gllb, &
44 : oe_lb, oe_none, oe_saop, oe_shift
45 : USE input_section_types, ONLY: section_vals_create,&
46 : section_vals_get,&
47 : section_vals_get_subs_vals,&
48 : section_vals_release,&
49 : section_vals_retain,&
50 : section_vals_set_subs_vals,&
51 : section_vals_type,&
52 : section_vals_val_get
53 : USE kinds, ONLY: dp,&
54 : int_8
55 : USE message_passing, ONLY: mp_para_env_type
56 : USE parallel_gemm_api, ONLY: parallel_gemm
57 : USE physcon, ONLY: evolt
58 : USE qs_environment_types, ONLY: get_qs_env,&
59 : qs_environment_type
60 : USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
61 : USE qs_ks_types, ONLY: qs_ks_env_type,&
62 : set_ks_env
63 : USE qs_mo_types, ONLY: allocate_mo_set,&
64 : deallocate_mo_set,&
65 : get_mo_set,&
66 : init_mo_set,&
67 : mo_set_type
68 : USE qs_scf_methods, ONLY: eigensolver
69 : USE qs_scf_post_gpw, ONLY: make_lumo_gpw
70 : USE qs_scf_types, ONLY: ot_method_nr,&
71 : qs_scf_env_type
72 : USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos
73 : USE tblite_ks_matrix, ONLY: build_tblite_ks_matrix
74 : USE util, ONLY: sort
75 : USE xc_pot_saop, ONLY: add_saop_pot
76 : USE xtb_ks_matrix, ONLY: build_xtb_ks_matrix
77 : #include "./base/base_uses.f90"
78 :
79 : IMPLICIT NONE
80 :
81 : PRIVATE
82 :
83 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_utils'
84 :
85 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
86 : ! number of first derivative components (3: d/dx, d/dy, d/dz)
87 : INTEGER, PARAMETER, PRIVATE :: nderivs = 3
88 : INTEGER, PARAMETER, PRIVATE :: maxspins = 2
89 :
90 : PUBLIC :: tddfpt_init_ground_state_mos, tddfpt_release_ground_state_mos
91 : PUBLIC :: tddfpt_guess_vectors, tddfpt_init_mos, tddfpt_oecorr
92 : PUBLIC :: tddfpt_total_number_of_states
93 :
94 : ! **************************************************************************************************
95 :
96 : CONTAINS
97 :
98 : ! **************************************************************************************************
99 : !> \brief Prepare MOs for TDDFPT Calculations
100 : !> \param qs_env Quickstep environment
101 : !> \param gs_mos ...
102 : !> \param iounit ...
103 : ! **************************************************************************************************
104 1416 : SUBROUTINE tddfpt_init_mos(qs_env, gs_mos, iounit)
105 : TYPE(qs_environment_type), POINTER :: qs_env
106 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
107 : POINTER :: gs_mos
108 : INTEGER, INTENT(IN) :: iounit
109 :
110 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_init_mos'
111 :
112 : INTEGER :: handle, ispin, nmo_avail, nmo_occ, &
113 : nmo_virt, nspins
114 : INTEGER, DIMENSION(2, 2) :: moc, mvt
115 : LOGICAL :: print_virtuals_newtonx
116 1416 : REAL(kind=dp), DIMENSION(:), POINTER :: evals_virt_spin
117 : TYPE(cell_type), POINTER :: cell
118 1416 : TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: evals_virt
119 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
120 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
121 1416 : TARGET :: mos_virt
122 : TYPE(cp_fm_type), POINTER :: mos_virt_spin
123 1416 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
124 : TYPE(dft_control_type), POINTER :: dft_control
125 : TYPE(excited_energy_type), POINTER :: ex_env
126 1416 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
127 : TYPE(qs_ks_env_type), POINTER :: ks_env
128 : TYPE(qs_scf_env_type), POINTER :: scf_env
129 : TYPE(section_vals_type), POINTER :: print_section
130 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
131 :
132 1416 : CALL timeset(routineN, handle)
133 :
134 : CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, &
135 1416 : matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, scf_env=scf_env)
136 1416 : tddfpt_control => dft_control%tddfpt2_control
137 1416 : IF ((tddfpt_control%do_bse) .OR. (tddfpt_control%do_bse_w_only) .OR. &
138 : (tddfpt_control%do_bse_gw_only)) THEN
139 4 : NULLIFY (ks_env, ex_env)
140 4 : CALL get_qs_env(qs_env, exstate_env=ex_env, ks_env=ks_env)
141 4 : CALL dbcsr_copy(matrix_ks(1)%matrix, ex_env%matrix_ks(1)%matrix)
142 4 : CALL set_ks_env(ks_env, matrix_ks=matrix_ks)
143 : END IF
144 :
145 1416 : CPASSERT(.NOT. ASSOCIATED(gs_mos))
146 : ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals
147 1416 : nspins = dft_control%nspins
148 5858 : ALLOCATE (gs_mos(nspins))
149 :
150 : ! check if virtuals should be constructed for NAMD interface with NEWTONX
151 1416 : print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
152 1416 : CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals_newtonx)
153 :
154 : ! when the number of unoccupied orbitals is limited and OT has been used
155 : ! for the ground-state DFT calculation,
156 : ! compute the missing unoccupied orbitals using OT as well.
157 1416 : NULLIFY (evals_virt, evals_virt_spin, mos_virt_spin)
158 1416 : IF (ASSOCIATED(scf_env)) THEN
159 1416 : IF ((scf_env%method == ot_method_nr .AND. tddfpt_control%nlumo > 0) .OR. &
160 : (scf_env%method == ot_method_nr .AND. print_virtuals_newtonx)) THEN
161 : ! As OT with ADDED_MOS/=0 is currently not implemented, the following block is equivalent to:
162 : ! nmo_virt = tddfpt_control%nlumo
163 : ! number of already computed unoccupied orbitals (added_mos) .
164 2 : nmo_virt = HUGE(0)
165 4 : DO ispin = 1, nspins
166 2 : CALL get_mo_set(mos(ispin), nmo=nmo_avail, homo=nmo_occ)
167 4 : nmo_virt = MIN(nmo_virt, nmo_avail - nmo_occ)
168 : END DO
169 : ! number of unoccupied orbitals to compute
170 2 : nmo_virt = tddfpt_control%nlumo - nmo_virt
171 2 : IF (.NOT. print_virtuals_newtonx) THEN
172 0 : IF (nmo_virt > 0) THEN
173 0 : ALLOCATE (evals_virt(nspins), mos_virt(nspins))
174 : ! the number of actually computed unoccupied orbitals will be stored as 'nmo_avail'
175 0 : CALL make_lumo_gpw(qs_env, scf_env, mos_virt, evals_virt, nmo_virt, nmo_avail)
176 : END IF
177 : END IF
178 : END IF
179 : END IF
180 :
181 3026 : DO ispin = 1, nspins
182 1610 : IF (ASSOCIATED(evals_virt)) THEN
183 0 : evals_virt_spin => evals_virt(ispin)%array
184 : ELSE
185 1610 : NULLIFY (evals_virt_spin)
186 : END IF
187 1610 : IF (ALLOCATED(mos_virt)) THEN
188 0 : mos_virt_spin => mos_virt(ispin)
189 : ELSE
190 1610 : NULLIFY (mos_virt_spin)
191 : END IF
192 : CALL tddfpt_init_ground_state_mos(gs_mos=gs_mos(ispin), mo_set=mos(ispin), &
193 : nlumo=tddfpt_control%nlumo, &
194 : blacs_env=blacs_env, cholesky_method=cholesky_restore, &
195 : matrix_ks=matrix_ks(ispin)%matrix, matrix_s=matrix_s(1)%matrix, &
196 : mos_virt=mos_virt_spin, evals_virt=evals_virt_spin, &
197 3026 : qs_env=qs_env)
198 : END DO
199 :
200 1416 : moc = 0
201 1416 : mvt = 0
202 3026 : DO ispin = 1, nspins
203 1610 : CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, nrow_global=moc(1, ispin), ncol_global=moc(2, ispin))
204 3026 : CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, nrow_global=mvt(1, ispin), ncol_global=mvt(2, ispin))
205 : END DO
206 1416 : IF (iounit > 0) THEN
207 708 : WRITE (iounit, "(T2,A,T36,A)") "TDDFPT| Molecular Orbitals:", &
208 1416 : " Spin AOs Occ Virt Total"
209 1513 : DO ispin = 1, nspins
210 805 : WRITE (iounit, "(T2,A,T37,I4,4I10)") "TDDFPT| ", ispin, moc(1, ispin), moc(2, ispin), &
211 2318 : mvt(2, ispin), moc(2, ispin) + mvt(2, ispin)
212 : END DO
213 : END IF
214 :
215 1416 : IF (ASSOCIATED(evals_virt)) THEN
216 0 : DO ispin = 1, SIZE(evals_virt)
217 0 : IF (ASSOCIATED(evals_virt(ispin)%array)) DEALLOCATE (evals_virt(ispin)%array)
218 : END DO
219 0 : DEALLOCATE (evals_virt)
220 : END IF
221 :
222 1416 : CALL cp_fm_release(mos_virt)
223 :
224 1416 : CALL timestop(handle)
225 :
226 4248 : END SUBROUTINE tddfpt_init_mos
227 :
228 : ! **************************************************************************************************
229 : !> \brief Generate all virtual molecular orbitals for a given spin by diagonalising
230 : !> the corresponding Kohn-Sham matrix.
231 : !> \param gs_mos structure to store occupied and virtual molecular orbitals
232 : !> (allocated and initialised on exit)
233 : !> \param mo_set ground state molecular orbitals for a given spin
234 : !> \param nlumo number of unoccupied states to consider (-1 means all states)
235 : !> \param blacs_env BLACS parallel environment
236 : !> \param cholesky_method Cholesky method to compute the inverse overlap matrix
237 : !> \param matrix_ks Kohn-Sham matrix for a given spin
238 : !> \param matrix_s overlap matrix
239 : !> \param mos_virt precomputed (OT) expansion coefficients of virtual molecular orbitals
240 : !> (in addition to the ADDED_MOS, if present). NULL when no OT is in use.
241 : !> \param evals_virt orbital energies of precomputed (OT) virtual molecular orbitals.
242 : !> NULL when no OT is in use.
243 : !> \param qs_env ...
244 : !> \par History
245 : !> * 05.2016 created as tddfpt_lumos() [Sergey Chulkov]
246 : !> * 06.2016 renamed, altered prototype [Sergey Chulkov]
247 : !> * 04.2019 limit the number of unoccupied states, orbital energy correction [Sergey Chulkov]
248 : ! **************************************************************************************************
249 1610 : SUBROUTINE tddfpt_init_ground_state_mos(gs_mos, mo_set, nlumo, blacs_env, cholesky_method, matrix_ks, matrix_s, &
250 : mos_virt, evals_virt, qs_env)
251 : TYPE(tddfpt_ground_state_mos) :: gs_mos
252 : TYPE(mo_set_type), INTENT(IN) :: mo_set
253 : INTEGER, INTENT(in) :: nlumo
254 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
255 : INTEGER, INTENT(in) :: cholesky_method
256 : TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_s
257 : TYPE(cp_fm_type), INTENT(IN), POINTER :: mos_virt
258 : REAL(kind=dp), DIMENSION(:), POINTER :: evals_virt
259 : TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
260 :
261 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_init_ground_state_mos'
262 : REAL(kind=dp), PARAMETER :: eps_dp = EPSILON(0.0_dp)
263 :
264 : INTEGER :: cholesky_method_inout, handle, icol_global, icol_local, imo, iounit, irow_global, &
265 : irow_local, nao, ncol_local, nelectrons, nmo_occ, nmo_scf, nmo_virt, nrow_local, sign_int
266 1610 : INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, &
267 1610 : sum_sign_array
268 1610 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
269 : LOGICAL :: do_eigen, print_phases
270 : REAL(kind=dp) :: element, maxocc
271 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
272 1610 : POINTER :: my_block
273 1610 : REAL(kind=dp), DIMENSION(:), POINTER :: mo_evals_extended, mo_occ_extended, &
274 1610 : mo_occ_scf
275 : TYPE(cp_fm_struct_type), POINTER :: ao_ao_fm_struct, ao_mo_occ_fm_struct, &
276 : ao_mo_virt_fm_struct, wfn_fm_struct
277 : TYPE(cp_fm_type) :: matrix_ks_fm, ortho_fm, work_fm, &
278 : work_fm_virt
279 : TYPE(cp_fm_type), POINTER :: mo_coeff_extended
280 : TYPE(cp_logger_type), POINTER :: logger
281 : TYPE(mo_set_type), POINTER :: mos_extended
282 : TYPE(mp_para_env_type), POINTER :: para_env
283 : TYPE(section_vals_type), POINTER :: print_section
284 :
285 1610 : CALL timeset(routineN, handle)
286 :
287 1610 : NULLIFY (logger)
288 1610 : logger => cp_get_default_logger()
289 1610 : iounit = cp_logger_get_default_io_unit(logger)
290 :
291 1610 : CALL blacs_env%get(para_env=para_env)
292 :
293 : CALL get_mo_set(mo_set, nao=nao, nmo=nmo_scf, homo=nmo_occ, maxocc=maxocc, &
294 1610 : nelectron=nelectrons, occupation_numbers=mo_occ_scf)
295 :
296 1610 : print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
297 1610 : CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases)
298 :
299 1610 : nmo_virt = nao - nmo_occ
300 1610 : IF (nlumo >= 0) THEN
301 4 : nmo_virt = MIN(nmo_virt, nlumo)
302 : END IF
303 :
304 1610 : IF (nmo_virt <= 0) THEN
305 : CALL cp_abort(__LOCATION__, &
306 0 : 'At least one unoccupied molecular orbital is required to calculate excited states.')
307 : END IF
308 :
309 1610 : do_eigen = .FALSE.
310 : ! diagonalise the Kohn-Sham matrix one more time if the number of available unoccupied states are too small
311 1610 : IF (ASSOCIATED(evals_virt)) THEN
312 0 : CPASSERT(ASSOCIATED(mos_virt))
313 0 : IF (nmo_virt > nmo_scf - nmo_occ + SIZE(evals_virt)) do_eigen = .TRUE.
314 : ELSE
315 1610 : IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .TRUE.
316 : END IF
317 :
318 : ! ++ allocate storage space for gs_mos
319 1610 : NULLIFY (ao_mo_occ_fm_struct, ao_mo_virt_fm_struct)
320 : ! Tiny fix (A.Sinyavskiy)
321 : CALL cp_fm_struct_create(ao_mo_occ_fm_struct, template_fmstruct=mo_set%mo_coeff%matrix_struct, &
322 1610 : ncol_global=nmo_occ, context=blacs_env)
323 : CALL cp_fm_struct_create(ao_mo_virt_fm_struct, template_fmstruct=mo_set%mo_coeff%matrix_struct, &
324 1610 : ncol_global=nmo_virt, context=blacs_env)
325 :
326 1610 : NULLIFY (gs_mos%mos_occ, gs_mos%mos_virt, gs_mos%evals_occ_matrix)
327 1610 : ALLOCATE (gs_mos%mos_occ, gs_mos%mos_virt)
328 1610 : CALL cp_fm_create(gs_mos%mos_occ, ao_mo_occ_fm_struct)
329 1610 : CALL cp_fm_create(gs_mos%mos_virt, ao_mo_virt_fm_struct)
330 1610 : gs_mos%nmo_occ = nmo_occ
331 :
332 4830 : ALLOCATE (gs_mos%evals_occ(nmo_occ))
333 4830 : ALLOCATE (gs_mos%evals_virt(nmo_virt))
334 3220 : ALLOCATE (gs_mos%phases_occ(nmo_occ))
335 3220 : ALLOCATE (gs_mos%phases_virt(nmo_virt))
336 :
337 : ! ++ nullify pointers
338 1610 : NULLIFY (ao_ao_fm_struct, wfn_fm_struct)
339 1610 : NULLIFY (mos_extended, mo_coeff_extended, mo_evals_extended, mo_occ_extended)
340 :
341 1610 : IF (do_eigen) THEN
342 : ! ++ set of molecular orbitals
343 1600 : CALL cp_fm_struct_create(ao_ao_fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
344 1600 : CALL cp_fm_struct_create(wfn_fm_struct, nrow_global=nao, ncol_global=nmo_occ + nmo_virt, context=blacs_env)
345 1600 : ALLOCATE (mos_extended)
346 : CALL allocate_mo_set(mos_extended, nao, nmo_occ + nmo_virt, nelectrons, &
347 1600 : REAL(nelectrons, dp), maxocc, flexible_electron_count=0.0_dp)
348 1600 : CALL init_mo_set(mos_extended, fm_struct=wfn_fm_struct, name="mos-extended")
349 1600 : CALL cp_fm_struct_release(wfn_fm_struct)
350 : CALL get_mo_set(mos_extended, mo_coeff=mo_coeff_extended, &
351 1600 : eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
352 :
353 : ! use the explicit loop in order to avoid temporary arrays.
354 : !
355 : ! The assignment statement : mo_occ_extended(1:nmo_scf) = mo_occ_scf(1:nmo_scf)
356 : ! implies temporary arrays as a compiler does not know in advance that the pointers
357 : ! on both sides of the statement point to non-overlapped memory regions
358 9778 : DO imo = 1, nmo_scf
359 9778 : mo_occ_extended(imo) = mo_occ_scf(imo)
360 : END DO
361 28468 : mo_occ_extended(nmo_scf + 1:) = 0.0_dp
362 :
363 : ! ++ allocate temporary matrices
364 1600 : CALL cp_fm_create(matrix_ks_fm, ao_ao_fm_struct)
365 1600 : CALL cp_fm_create(ortho_fm, ao_ao_fm_struct)
366 1600 : CALL cp_fm_create(work_fm, ao_ao_fm_struct)
367 1600 : CALL cp_fm_struct_release(ao_ao_fm_struct)
368 :
369 : ! some stuff from the subroutine general_eigenproblem()
370 1600 : CALL copy_dbcsr_to_fm(matrix_s, ortho_fm)
371 1600 : CALL copy_dbcsr_to_fm(matrix_ks, matrix_ks_fm)
372 :
373 1600 : IF (cholesky_method == cholesky_dbcsr) THEN
374 0 : CPABORT('CHOLESKY DBCSR_INVERSE is not implemented in TDDFT.')
375 1600 : ELSE IF (cholesky_method == cholesky_off) THEN
376 0 : CPABORT('CHOLESKY OFF is not implemented in TDDFT.')
377 : ELSE
378 1600 : CALL cp_fm_cholesky_decompose(ortho_fm)
379 1600 : IF (cholesky_method == cholesky_inverse) THEN
380 0 : CALL cp_fm_triangular_invert(ortho_fm)
381 : END IF
382 :
383 : ! need to store 'cholesky_method' in a temporary variable, as the subroutine eigensolver()
384 : ! will update this variable
385 1600 : cholesky_method_inout = cholesky_method
386 : CALL eigensolver(matrix_ks_fm=matrix_ks_fm, mo_set=mos_extended, ortho=ortho_fm, &
387 : work=work_fm, cholesky_method=cholesky_method_inout, &
388 1600 : do_level_shift=.FALSE., level_shift=0.0_dp, use_jacobi=.FALSE.)
389 : END IF
390 :
391 : ! -- clean up needless matrices
392 1600 : CALL cp_fm_release(work_fm)
393 1600 : CALL cp_fm_release(ortho_fm)
394 1600 : CALL cp_fm_release(matrix_ks_fm)
395 : ELSE
396 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff_extended, &
397 10 : eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
398 : END IF
399 :
400 : ! compute the phase of molecular orbitals;
401 : ! matrix work_fm holds occupied molecular orbital coefficients distributed among all the processors
402 : !CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
403 1610 : CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
404 1610 : CALL cp_fm_struct_release(ao_mo_occ_fm_struct)
405 :
406 1610 : CALL cp_fm_to_fm(mo_coeff_extended, work_fm, ncol=nmo_occ, source_start=1, target_start=1)
407 : CALL cp_fm_get_info(work_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
408 1610 : row_indices=row_indices, col_indices=col_indices, local_data=my_block)
409 :
410 8050 : ALLOCATE (minrow_neg_array(nmo_occ), minrow_pos_array(nmo_occ), sum_sign_array(nmo_occ))
411 9812 : minrow_neg_array(:) = nao
412 9812 : minrow_pos_array(:) = nao
413 1610 : sum_sign_array(:) = 0
414 9812 : DO icol_local = 1, ncol_local
415 8202 : icol_global = col_indices(icol_local)
416 :
417 247683 : DO irow_local = 1, nrow_local
418 237871 : element = my_block(irow_local, icol_local)
419 :
420 237871 : sign_int = 0
421 237871 : IF (element >= eps_dp) THEN
422 : sign_int = 1
423 119191 : ELSE IF (element <= -eps_dp) THEN
424 118414 : sign_int = -1
425 : END IF
426 :
427 237871 : sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
428 :
429 237871 : irow_global = row_indices(irow_local)
430 246073 : IF (sign_int > 0) THEN
431 118680 : IF (minrow_pos_array(icol_global) > irow_global) THEN
432 8039 : minrow_pos_array(icol_global) = irow_global
433 : END IF
434 119191 : ELSE IF (sign_int < 0) THEN
435 118414 : IF (minrow_neg_array(icol_global) > irow_global) THEN
436 7885 : minrow_neg_array(icol_global) = irow_global
437 : END IF
438 : END IF
439 : END DO
440 : END DO
441 :
442 1610 : CALL para_env%sum(sum_sign_array)
443 1610 : CALL para_env%min(minrow_neg_array)
444 1610 : CALL para_env%min(minrow_pos_array)
445 :
446 9812 : DO icol_local = 1, nmo_occ
447 9812 : IF (sum_sign_array(icol_local) > 0) THEN
448 : ! most of the expansion coefficients are positive => MO's phase = +1
449 4098 : gs_mos%phases_occ(icol_local) = 1.0_dp
450 4104 : ELSE IF (sum_sign_array(icol_local) < 0) THEN
451 : ! most of the expansion coefficients are negative => MO's phase = -1
452 3742 : gs_mos%phases_occ(icol_local) = -1.0_dp
453 : ELSE
454 : ! equal number of positive and negative expansion coefficients
455 362 : IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local)) THEN
456 : ! the first positive expansion coefficient has a lower index then
457 : ! the first negative expansion coefficient; MO's phase = +1
458 204 : gs_mos%phases_occ(icol_local) = 1.0_dp
459 : ELSE
460 : ! MO's phase = -1
461 158 : gs_mos%phases_occ(icol_local) = -1.0_dp
462 : END IF
463 : END IF
464 : END DO
465 :
466 1610 : DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
467 :
468 : ! return the requested occupied and virtual molecular orbitals and corresponding orbital energies
469 1610 : CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_occ, ncol=nmo_occ, source_start=1, target_start=1)
470 9812 : gs_mos%evals_occ(1:nmo_occ) = mo_evals_extended(1:nmo_occ)
471 :
472 1610 : IF (ASSOCIATED(evals_virt) .AND. (.NOT. do_eigen) .AND. nmo_virt > nmo_scf - nmo_occ) THEN
473 : CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_scf - nmo_occ, &
474 0 : source_start=nmo_occ + 1, target_start=1)
475 : CALL cp_fm_to_fm(mos_virt, gs_mos%mos_virt, ncol=nmo_virt - (nmo_scf - nmo_occ), &
476 0 : source_start=1, target_start=nmo_scf - nmo_occ + 1)
477 0 : gs_mos%evals_virt(1:nmo_scf - nmo_occ) = evals_virt(nmo_occ + 1:nmo_occ + nmo_scf)
478 0 : gs_mos%evals_virt(nmo_scf - nmo_occ + 1:nmo_virt) = evals_virt(1:nmo_virt - (nmo_scf - nmo_occ))
479 : ELSE
480 1610 : CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_virt, source_start=nmo_occ + 1, target_start=1)
481 28650 : gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt)
482 : END IF
483 :
484 1610 : IF (print_phases) THEN
485 : ! compute the phase of molecular orbitals;
486 : ! matrix work_fm holds virtual molecular orbital coefficients distributed among all the processors
487 : !CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
488 0 : CALL cp_fm_create(work_fm_virt, ao_mo_virt_fm_struct)
489 :
490 0 : CALL cp_fm_to_fm(gs_mos%mos_virt, work_fm_virt, ncol=nmo_virt, source_start=1, target_start=1)
491 : CALL cp_fm_get_info(work_fm_virt, nrow_local=nrow_local, ncol_local=ncol_local, &
492 0 : row_indices=row_indices, col_indices=col_indices, local_data=my_block)
493 :
494 0 : ALLOCATE (minrow_neg_array(nmo_virt), minrow_pos_array(nmo_virt), sum_sign_array(nmo_virt))
495 0 : minrow_neg_array(:) = nao
496 0 : minrow_pos_array(:) = nao
497 0 : sum_sign_array(:) = 0
498 0 : DO icol_local = 1, ncol_local
499 0 : icol_global = col_indices(icol_local)
500 :
501 0 : DO irow_local = 1, nrow_local
502 0 : element = my_block(irow_local, icol_local)
503 :
504 0 : sign_int = 0
505 0 : IF (element >= eps_dp) THEN
506 : sign_int = 1
507 0 : ELSE IF (element <= -eps_dp) THEN
508 0 : sign_int = -1
509 : END IF
510 :
511 0 : sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
512 :
513 0 : irow_global = row_indices(irow_local)
514 0 : IF (sign_int > 0) THEN
515 0 : IF (minrow_pos_array(icol_global) > irow_global) THEN
516 0 : minrow_pos_array(icol_global) = irow_global
517 : END IF
518 0 : ELSE IF (sign_int < 0) THEN
519 0 : IF (minrow_neg_array(icol_global) > irow_global) THEN
520 0 : minrow_neg_array(icol_global) = irow_global
521 : END IF
522 : END IF
523 : END DO
524 : END DO
525 :
526 0 : CALL para_env%sum(sum_sign_array)
527 0 : CALL para_env%min(minrow_neg_array)
528 0 : CALL para_env%min(minrow_pos_array)
529 0 : DO icol_local = 1, nmo_virt
530 0 : IF (sum_sign_array(icol_local) > 0) THEN
531 : ! most of the expansion coefficients are positive => MO's phase = +1
532 0 : gs_mos%phases_virt(icol_local) = 1.0_dp
533 0 : ELSE IF (sum_sign_array(icol_local) < 0) THEN
534 : ! most of the expansion coefficients are negative => MO's phase = -1
535 0 : gs_mos%phases_virt(icol_local) = -1.0_dp
536 : ELSE
537 : ! equal number of positive and negative expansion coefficients
538 0 : IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local)) THEN
539 : ! the first positive expansion coefficient has a lower index then
540 : ! the first negative expansion coefficient; MO's phase = +1
541 0 : gs_mos%phases_virt(icol_local) = 1.0_dp
542 : ELSE
543 : ! MO's phase = -1
544 0 : gs_mos%phases_virt(icol_local) = -1.0_dp
545 : END IF
546 : END IF
547 : END DO
548 0 : DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
549 0 : CALL cp_fm_release(work_fm_virt)
550 : END IF !print_phases
551 1610 : CALL cp_fm_struct_release(ao_mo_virt_fm_struct) ! here after print_phases
552 :
553 1610 : CALL cp_fm_release(work_fm)
554 :
555 1610 : IF (do_eigen) THEN
556 1600 : CALL deallocate_mo_set(mos_extended)
557 1600 : DEALLOCATE (mos_extended)
558 : END IF
559 :
560 1610 : CALL timestop(handle)
561 :
562 9660 : END SUBROUTINE tddfpt_init_ground_state_mos
563 :
564 : ! **************************************************************************************************
565 : !> \brief Release molecular orbitals.
566 : !> \param gs_mos structure that holds occupied and virtual molecular orbitals
567 : !> \par History
568 : !> * 06.2016 created [Sergey Chulkov]
569 : ! **************************************************************************************************
570 1610 : SUBROUTINE tddfpt_release_ground_state_mos(gs_mos)
571 : TYPE(tddfpt_ground_state_mos), INTENT(inout) :: gs_mos
572 :
573 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_release_ground_state_mos'
574 :
575 : INTEGER :: handle
576 :
577 1610 : CALL timeset(routineN, handle)
578 :
579 1610 : IF (ALLOCATED(gs_mos%phases_occ)) THEN
580 1610 : DEALLOCATE (gs_mos%phases_occ)
581 : END IF
582 :
583 1610 : IF (ALLOCATED(gs_mos%evals_virt)) THEN
584 1610 : DEALLOCATE (gs_mos%evals_virt)
585 : END IF
586 :
587 1610 : IF (ALLOCATED(gs_mos%evals_occ)) THEN
588 1610 : DEALLOCATE (gs_mos%evals_occ)
589 : END IF
590 :
591 1610 : IF (ALLOCATED(gs_mos%phases_virt)) THEN
592 1610 : DEALLOCATE (gs_mos%phases_virt)
593 : END IF
594 :
595 1610 : IF (ALLOCATED(gs_mos%index_active)) THEN
596 1610 : DEALLOCATE (gs_mos%index_active)
597 : END IF
598 :
599 1610 : IF (ASSOCIATED(gs_mos%evals_occ_matrix)) THEN
600 38 : CALL cp_fm_release(gs_mos%evals_occ_matrix)
601 38 : DEALLOCATE (gs_mos%evals_occ_matrix)
602 : END IF
603 :
604 1610 : IF (ASSOCIATED(gs_mos%mos_virt)) THEN
605 1610 : CALL cp_fm_release(gs_mos%mos_virt)
606 1610 : DEALLOCATE (gs_mos%mos_virt)
607 : END IF
608 :
609 1610 : IF (ASSOCIATED(gs_mos%mos_occ)) THEN
610 1610 : CALL cp_fm_release(gs_mos%mos_occ)
611 1610 : DEALLOCATE (gs_mos%mos_occ)
612 : END IF
613 :
614 1610 : IF (ASSOCIATED(gs_mos%mos_active)) THEN
615 1610 : CALL cp_fm_release(gs_mos%mos_active)
616 1610 : DEALLOCATE (gs_mos%mos_active)
617 : END IF
618 :
619 1610 : CALL timestop(handle)
620 :
621 1610 : END SUBROUTINE tddfpt_release_ground_state_mos
622 :
623 : ! **************************************************************************************************
624 : !> \brief Callculate orbital corrected KS matrix for TDDFPT
625 : !> \param qs_env Quickstep environment
626 : !> \param gs_mos ...
627 : !> \param matrix_ks_oep ...
628 : ! **************************************************************************************************
629 1416 : SUBROUTINE tddfpt_oecorr(qs_env, gs_mos, matrix_ks_oep)
630 : TYPE(qs_environment_type), POINTER :: qs_env
631 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
632 : POINTER :: gs_mos
633 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_oep
634 :
635 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_oecorr'
636 :
637 : INTEGER :: handle, iounit, ispin, nao, nmo_occ, &
638 : nspins
639 : LOGICAL :: do_hfx
640 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
641 : TYPE(cp_fm_struct_type), POINTER :: ao_mo_occ_fm_struct, &
642 : mo_occ_mo_occ_fm_struct
643 : TYPE(cp_fm_type) :: work_fm
644 : TYPE(cp_logger_type), POINTER :: logger
645 1416 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
646 : TYPE(dft_control_type), POINTER :: dft_control
647 : TYPE(section_vals_type), POINTER :: hfx_section, xc_fun_empty, &
648 : xc_fun_original
649 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
650 :
651 1416 : CALL timeset(routineN, handle)
652 :
653 1416 : NULLIFY (logger)
654 1416 : logger => cp_get_default_logger()
655 1416 : iounit = cp_logger_get_default_io_unit(logger)
656 :
657 1416 : CALL get_qs_env(qs_env, blacs_env=blacs_env, dft_control=dft_control, matrix_ks=matrix_ks)
658 1416 : tddfpt_control => dft_control%tddfpt2_control
659 :
660 : ! obtain corrected KS-matrix
661 : ! We should 'save' the energy values?
662 1416 : nspins = SIZE(gs_mos)
663 1416 : NULLIFY (matrix_ks_oep)
664 1416 : IF (tddfpt_control%oe_corr /= oe_none) THEN
665 32 : IF (iounit > 0) THEN
666 16 : WRITE (iounit, "(1X,A)") "", &
667 16 : "-------------------------------------------------------------------------------", &
668 16 : "- Orbital Eigenvalue Correction Started -", &
669 32 : "-------------------------------------------------------------------------------"
670 : END IF
671 :
672 : CALL cp_warn(__LOCATION__, &
673 : "Orbital energy correction potential is an experimental feature. "// &
674 32 : "Use it with extreme care")
675 :
676 32 : hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
677 32 : CALL section_vals_get(hfx_section, explicit=do_hfx)
678 32 : IF (do_hfx) THEN
679 : CALL cp_abort(__LOCATION__, &
680 : "Implementation of orbital energy correction XC-potentials is "// &
681 0 : "currently incompatible with exact-exchange functionals")
682 : END IF
683 :
684 32 : CALL dbcsr_allocate_matrix_set(matrix_ks_oep, nspins)
685 70 : DO ispin = 1, nspins
686 38 : CALL dbcsr_init_p(matrix_ks_oep(ispin)%matrix)
687 70 : CALL dbcsr_copy(matrix_ks_oep(ispin)%matrix, matrix_ks(ispin)%matrix)
688 : END DO
689 :
690 : ! KS-matrix without XC-terms
691 32 : xc_fun_original => section_vals_get_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL")
692 32 : CALL section_vals_retain(xc_fun_original)
693 32 : NULLIFY (xc_fun_empty)
694 32 : CALL section_vals_create(xc_fun_empty, xc_fun_original%section)
695 32 : CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_empty)
696 32 : CALL section_vals_release(xc_fun_empty)
697 :
698 32 : IF (dft_control%qs_control%semi_empirical) THEN
699 0 : CPABORT("TDDFPT with SE not possible")
700 32 : ELSE IF (dft_control%qs_control%dftb) THEN
701 0 : CPABORT("TDDFPT with DFTB not possible")
702 32 : ELSE IF (dft_control%qs_control%xtb) THEN
703 18 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
704 : CALL build_tblite_ks_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
705 2 : ext_ks_matrix=matrix_ks_oep)
706 : ELSE
707 : CALL build_xtb_ks_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
708 16 : ext_ks_matrix=matrix_ks_oep)
709 : END IF
710 : ELSE
711 : CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
712 14 : ext_ks_matrix=matrix_ks_oep)
713 : END IF
714 :
715 : IF (tddfpt_control%oe_corr == oe_saop .OR. &
716 32 : tddfpt_control%oe_corr == oe_lb .OR. &
717 : tddfpt_control%oe_corr == oe_gllb) THEN
718 14 : IF (iounit > 0) THEN
719 7 : WRITE (iounit, "(T2,A)") " Orbital energy correction of SAOP type "
720 : END IF
721 14 : CALL add_saop_pot(matrix_ks_oep, qs_env, tddfpt_control%oe_corr)
722 18 : ELSE IF (tddfpt_control%oe_corr == oe_shift) THEN
723 18 : IF (iounit > 0) THEN
724 : WRITE (iounit, "(T2,A,T71,F10.3)") &
725 9 : " Virtual Orbital Eigenvalue Shift [eV] ", tddfpt_control%ev_shift*evolt
726 : WRITE (iounit, "(T2,A,T71,F10.3)") &
727 9 : " Open Shell Orbital Eigenvalue Shift [eV] ", tddfpt_control%eos_shift*evolt
728 : END IF
729 : CALL ev_shift_operator(qs_env, gs_mos, matrix_ks_oep, &
730 18 : tddfpt_control%ev_shift, tddfpt_control%eos_shift)
731 : ELSE
732 : CALL cp_abort(__LOCATION__, &
733 0 : "Unimplemented orbital energy correction potential")
734 : END IF
735 32 : CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_original)
736 32 : CALL section_vals_release(xc_fun_original)
737 :
738 : ! compute 'evals_occ_matrix'
739 32 : CALL dbcsr_get_info(matrix_ks(1)%matrix, nfullrows_total=nao)
740 32 : NULLIFY (mo_occ_mo_occ_fm_struct)
741 70 : DO ispin = 1, nspins
742 38 : nmo_occ = SIZE(gs_mos(ispin)%evals_occ)
743 : CALL cp_fm_struct_create(mo_occ_mo_occ_fm_struct, nrow_global=nmo_occ, ncol_global=nmo_occ, &
744 38 : context=blacs_env)
745 38 : ALLOCATE (gs_mos(ispin)%evals_occ_matrix)
746 38 : CALL cp_fm_create(gs_mos(ispin)%evals_occ_matrix, mo_occ_mo_occ_fm_struct)
747 38 : CALL cp_fm_struct_release(mo_occ_mo_occ_fm_struct)
748 : ! work_fm is a temporary [nao x nmo_occ] matrix
749 : CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, &
750 38 : context=blacs_env)
751 38 : CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
752 38 : CALL cp_fm_struct_release(ao_mo_occ_fm_struct)
753 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks_oep(ispin)%matrix, gs_mos(ispin)%mos_occ, &
754 38 : work_fm, ncol=nmo_occ, alpha=1.0_dp, beta=0.0_dp)
755 : CALL parallel_gemm('T', 'N', nmo_occ, nmo_occ, nao, 1.0_dp, gs_mos(ispin)%mos_occ, work_fm, &
756 38 : 0.0_dp, gs_mos(ispin)%evals_occ_matrix)
757 108 : CALL cp_fm_release(work_fm)
758 : END DO
759 32 : IF (iounit > 0) THEN
760 : WRITE (iounit, "(1X,A)") &
761 16 : "-------------------------------------------------------------------------------"
762 : END IF
763 :
764 : END IF
765 :
766 1416 : CALL timestop(handle)
767 :
768 1416 : END SUBROUTINE tddfpt_oecorr
769 :
770 : ! **************************************************************************************************
771 : !> \brief Compute the number of possible singly excited states (occ -> virt)
772 : !> \param tddfpt_control ...
773 : !> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state
774 : !> \return the number of possible single excitations
775 : !> \par History
776 : !> * 01.2017 created [Sergey Chulkov]
777 : ! **************************************************************************************************
778 2241 : PURE FUNCTION tddfpt_total_number_of_states(tddfpt_control, gs_mos) RESULT(nstates_total)
779 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
780 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
781 : INTENT(in) :: gs_mos
782 : INTEGER(kind=int_8) :: nstates_total
783 :
784 : INTEGER :: ispin, nspins
785 :
786 2241 : nstates_total = 0
787 2241 : nspins = SIZE(gs_mos)
788 :
789 2241 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
790 : ! Total number of possible excitations for spin-conserving TDDFT
791 4680 : DO ispin = 1, nspins
792 : nstates_total = nstates_total + &
793 : gs_mos(ispin)%nmo_active* &
794 4680 : SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
795 : END DO
796 : ELSE
797 : ! Total number of possible excitations for spin-flip TDDFT
798 : nstates_total = gs_mos(1)%nmo_active* &
799 37 : SIZE(gs_mos(2)%evals_virt, kind=int_8)
800 : END IF
801 2241 : END FUNCTION tddfpt_total_number_of_states
802 :
803 : ! **************************************************************************************************
804 : !> \brief Create a shift operator on virtual/open shell space
805 : !> Shift operator = Edelta*Q Q: projector on virtual space (1-PS)
806 : !> projector on open shell space PosS
807 : !> \param qs_env the qs_env that is perturbed by this p_env
808 : !> \param gs_mos ...
809 : !> \param matrix_ks ...
810 : !> \param ev_shift ...
811 : !> \param eos_shift ...
812 : !> \par History
813 : !> 02.04.2019 adapted for TDDFT use from p_env (JGH)
814 : !> \author JGH
815 : ! **************************************************************************************************
816 18 : SUBROUTINE ev_shift_operator(qs_env, gs_mos, matrix_ks, ev_shift, eos_shift)
817 :
818 : TYPE(qs_environment_type), POINTER :: qs_env
819 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
820 : POINTER :: gs_mos
821 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
822 : REAL(KIND=dp), INTENT(IN) :: ev_shift, eos_shift
823 :
824 : CHARACTER(len=*), PARAMETER :: routineN = 'ev_shift_operator'
825 :
826 : INTEGER :: handle, ispin, n_spins, na, nb, nhomo, &
827 : nl, nos, nrow, nu, nvirt
828 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
829 : TYPE(cp_fm_type) :: cmos, cvec
830 : TYPE(cp_fm_type), POINTER :: coeff
831 18 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
832 : TYPE(dbcsr_type), POINTER :: smat
833 18 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
834 :
835 18 : CALL timeset(routineN, handle)
836 :
837 18 : n_spins = SIZE(gs_mos)
838 18 : CPASSERT(n_spins == SIZE(matrix_ks))
839 :
840 18 : IF (eos_shift /= 0.0_dp .AND. n_spins > 1) THEN
841 0 : CPABORT("eos_shift not implemented")
842 0 : CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
843 0 : smat => matrix_s(1)%matrix
844 0 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
845 0 : CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
846 0 : nl = MIN(na, nb)
847 0 : nu = MAX(na, nb)
848 : ! open shell orbital shift
849 0 : DO ispin = 1, n_spins
850 0 : coeff => gs_mos(ispin)%mos_occ
851 0 : CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
852 0 : IF (nhomo == nu) THEN
853 : ! downshift with -eos_shift using occupied orbitals
854 0 : nos = nu - nl
855 0 : CALL cp_fm_create(cmos, fmstruct)
856 0 : CALL cp_fm_get_info(coeff, nrow_global=nrow)
857 0 : CALL cp_fm_to_fm_submat(coeff, cmos, nrow, nos, 1, nl + 1, 1, 1)
858 0 : CALL cp_fm_create(cvec, fmstruct)
859 0 : CALL cp_dbcsr_sm_fm_multiply(smat, cmos, cvec, nos, 1.0_dp, 0.0_dp)
860 : CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
861 0 : alpha=-eos_shift, keep_sparsity=.TRUE.)
862 0 : CALL cp_fm_release(cmos)
863 0 : CALL cp_fm_release(cvec)
864 : ELSE
865 : ! upshift with eos_shift using virtual orbitals
866 0 : coeff => gs_mos(ispin)%mos_virt
867 0 : CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
868 0 : nos = nu - nhomo
869 0 : CPASSERT(nvirt >= nos)
870 0 : CALL cp_fm_create(cvec, fmstruct)
871 0 : CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nos, 1.0_dp, 0.0_dp)
872 : CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
873 0 : alpha=eos_shift, keep_sparsity=.TRUE.)
874 0 : CALL cp_fm_release(cvec)
875 : END IF
876 : END DO
877 : ! virtual shift
878 0 : IF (ev_shift /= 0.0_dp) THEN
879 0 : DO ispin = 1, n_spins
880 : CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
881 0 : alpha_scalar=1.0_dp, beta_scalar=ev_shift)
882 0 : coeff => gs_mos(ispin)%mos_occ
883 0 : CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
884 0 : CALL cp_fm_create(cvec, fmstruct)
885 0 : CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nhomo, 1.0_dp, 0.0_dp)
886 : CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nhomo, &
887 0 : alpha=-ev_shift, keep_sparsity=.TRUE.)
888 0 : CALL cp_fm_release(cvec)
889 0 : IF (nhomo < nu) THEN
890 0 : nos = nu - nhomo
891 0 : coeff => gs_mos(ispin)%mos_virt
892 0 : CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
893 0 : CPASSERT(nvirt >= nos)
894 0 : CALL cp_fm_create(cvec, fmstruct)
895 0 : CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nos, 1.0_dp, 0.0_dp)
896 : CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
897 0 : alpha=-ev_shift, keep_sparsity=.TRUE.)
898 0 : CALL cp_fm_release(cvec)
899 : END IF
900 : END DO
901 : END IF
902 : ELSE
903 : ! virtual shift
904 18 : IF (ev_shift /= 0.0_dp) THEN
905 18 : CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
906 18 : smat => matrix_s(1)%matrix
907 36 : DO ispin = 1, n_spins
908 : CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
909 18 : alpha_scalar=1.0_dp, beta_scalar=ev_shift)
910 18 : coeff => gs_mos(ispin)%mos_occ
911 18 : CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
912 18 : CALL cp_fm_create(cvec, fmstruct)
913 18 : CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nhomo, 1.0_dp, 0.0_dp)
914 : CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nhomo, &
915 18 : alpha=-ev_shift, keep_sparsity=.TRUE.)
916 54 : CALL cp_fm_release(cvec)
917 : END DO
918 : END IF
919 : END IF
920 : ! set eigenvalues
921 18 : IF (eos_shift == 0.0_dp .OR. n_spins == 1) THEN
922 36 : DO ispin = 1, n_spins
923 36 : IF (ALLOCATED(gs_mos(ispin)%evals_virt)) THEN
924 1346 : gs_mos(ispin)%evals_virt(:) = gs_mos(ispin)%evals_virt(:) + ev_shift
925 : END IF
926 : END DO
927 : ELSE
928 0 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
929 0 : CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
930 0 : nl = MIN(na, nb)
931 0 : nu = MAX(na, nb)
932 0 : nos = nu - nl
933 0 : IF (na == nu) THEN
934 0 : IF (ALLOCATED(gs_mos(1)%evals_occ)) THEN
935 0 : gs_mos(1)%evals_occ(nl + 1:nu) = gs_mos(1)%evals_occ(nl + 1:nu) - eos_shift
936 : END IF
937 0 : IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
938 0 : gs_mos(1)%evals_virt(:) = gs_mos(1)%evals_virt(:) + ev_shift
939 : END IF
940 0 : IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
941 0 : gs_mos(2)%evals_virt(1:nos) = gs_mos(2)%evals_virt(1:nos) + eos_shift
942 0 : gs_mos(2)%evals_virt(nos + 1:) = gs_mos(2)%evals_virt(nos + 1:) + ev_shift
943 : END IF
944 : ELSE
945 0 : IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
946 0 : gs_mos(1)%evals_virt(1:nos) = gs_mos(1)%evals_virt(1:nos) + eos_shift
947 0 : gs_mos(1)%evals_virt(nos + 1:) = gs_mos(1)%evals_virt(nos + 1:) + ev_shift
948 : END IF
949 0 : IF (ALLOCATED(gs_mos(2)%evals_occ)) THEN
950 0 : gs_mos(2)%evals_occ(nl + 1:nu) = gs_mos(2)%evals_occ(nl + 1:nu) - eos_shift
951 : END IF
952 0 : IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
953 0 : gs_mos(2)%evals_virt(:) = gs_mos(2)%evals_virt(:) + ev_shift
954 : END IF
955 : END IF
956 : END IF
957 :
958 18 : CALL timestop(handle)
959 :
960 18 : END SUBROUTINE ev_shift_operator
961 :
962 : ! **************************************************************************************************
963 : !> \brief Generate missed guess vectors.
964 : !> \param evects guess vectors distributed across all processors (initialised on exit)
965 : !> \param evals guessed transition energies (initialised on exit)
966 : !> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state
967 : !> \param log_unit output unit
968 : !> \param tddfpt_control ...
969 : !> \param fm_pool_ao_mo_active ...
970 : !> \param qs_env ...
971 : !> \param nspins ...
972 : !> \par History
973 : !> * 05.2016 created as tddfpt_guess() [Sergey Chulkov]
974 : !> * 06.2016 renamed, altered prototype, supports spin-polarised density [Sergey Chulkov]
975 : !> * 01.2017 simplified prototype, do not compute all possible singly-excited states
976 : !> [Sergey Chulkov]
977 : !> \note \parblock
978 : !> Based on the subroutine co_initial_guess() which was originally created by
979 : !> Thomas Chassaing on 06.2003.
980 : !>
981 : !> Only not associated guess vectors 'evects(spin, state)%matrix' are allocated and
982 : !> initialised; associated vectors assumed to be initialised elsewhere (e.g. using
983 : !> a restart file).
984 : !> \endparblock
985 : ! **************************************************************************************************
986 1426 : SUBROUTINE tddfpt_guess_vectors(evects, evals, gs_mos, log_unit, tddfpt_control, &
987 1426 : fm_pool_ao_mo_active, qs_env, nspins)
988 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout) :: evects
989 : REAL(kind=dp), DIMENSION(:), INTENT(inout) :: evals
990 : INTEGER, INTENT(in) :: nspins
991 : TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
992 : TYPE(cp_fm_pool_p_type), DIMENSION(:), INTENT(in) :: fm_pool_ao_mo_active
993 : TYPE(tddfpt2_control_type), INTENT(in), POINTER :: tddfpt_control
994 : INTEGER, INTENT(in) :: log_unit
995 : TYPE(tddfpt_ground_state_mos), DIMENSION(nspins), &
996 : INTENT(in) :: gs_mos
997 :
998 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_guess_vectors'
999 :
1000 : CHARACTER(len=5) :: spin_label1, spin_label2
1001 : INTEGER :: handle, i, imo_occ, imo_virt, ind, ispin, istate, j, jspin, k, no, nstates, &
1002 : nstates_occ_virt_alpha, nstates_selected, nv, spin1, spin2
1003 1426 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
1004 1426 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: reverse_index
1005 : INTEGER, DIMENSION(maxspins) :: nmo, nmo_occ_avail, nmo_occ_selected, &
1006 : nmo_virt_selected
1007 : REAL(kind=dp) :: e_occ
1008 1426 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: e_virt_minus_occ, ev_occ, ev_virt
1009 : TYPE(excited_energy_type), POINTER :: ex_env
1010 :
1011 1426 : CALL timeset(routineN, handle)
1012 :
1013 1426 : nstates = SIZE(evects, 2)
1014 :
1015 : IF (debug_this_module) THEN
1016 : CPASSERT(nstates > 0)
1017 : CPASSERT(nspins == 1 .OR. nspins == 2)
1018 : END IF
1019 :
1020 1426 : NULLIFY (ex_env)
1021 1426 : CALL get_qs_env(qs_env, exstate_env=ex_env)
1022 :
1023 3046 : DO ispin = 1, nspins
1024 : ! number of occupied orbitals for each spin component
1025 1620 : nmo_occ_avail(ispin) = gs_mos(ispin)%nmo_active
1026 1620 : nmo(ispin) = gs_mos(ispin)%nmo_occ
1027 : ! number of occupied and virtual orbitals which can potentially
1028 : ! contribute to the excited states in question.
1029 1620 : nmo_occ_selected(ispin) = MIN(nmo_occ_avail(ispin), nstates)
1030 3046 : nmo_virt_selected(ispin) = MIN(SIZE(gs_mos(ispin)%evals_virt), nstates)
1031 : END DO
1032 :
1033 : ! TO DO: the variable 'nstates_selected' should probably be declared as INTEGER(kind=int_8),
1034 : ! however we need a special version of the subroutine sort() in order to do so
1035 1426 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1036 2980 : nstates_selected = DOT_PRODUCT(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins))
1037 : ELSE
1038 22 : nstates_selected = nmo_occ_selected(1)*nmo_virt_selected(2)
1039 : END IF
1040 :
1041 4278 : ALLOCATE (inds(nstates_selected))
1042 4278 : ALLOCATE (e_virt_minus_occ(nstates_selected))
1043 :
1044 1426 : istate = 0
1045 1426 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1046 : ! Guess for spin-conserving TDDFT
1047 2980 : DO ispin = 1, nspins
1048 1576 : no = nmo_occ_selected(ispin)
1049 1576 : nv = nmo_virt_selected(ispin)
1050 7880 : ALLOCATE (ev_virt(nv), ev_occ(no))
1051 : ! if do_bse and do_gw, take gw zeroth order
1052 1576 : IF ((tddfpt_control%do_bse) .OR. (tddfpt_control%do_bse_w_only) .OR. &
1053 : (tddfpt_control%do_bse_gw_only)) THEN
1054 44 : ev_virt(1:nv) = ex_env%gw_eigen(nmo(ispin) + 1:nmo(ispin) + nv)
1055 20 : DO i = 1, no
1056 16 : j = nmo_occ_avail(ispin) - i + 1
1057 16 : k = gs_mos(ispin)%index_active(j)
1058 20 : ev_occ(i) = ex_env%gw_eigen(k)
1059 : END DO
1060 : ELSE
1061 5496 : ev_virt(1:nv) = gs_mos(ispin)%evals_virt(1:nv)
1062 5158 : DO i = 1, no
1063 3586 : j = nmo_occ_avail(ispin) - i + 1
1064 3586 : k = gs_mos(ispin)%index_active(j)
1065 5158 : ev_occ(i) = gs_mos(ispin)%evals_occ(k)
1066 : END DO
1067 : END IF
1068 :
1069 5178 : DO imo_occ = 1, nmo_occ_selected(ispin)
1070 : ! Here imo_occ enumerate Occupied orbitals in inverse order (from the last to the first element)
1071 3602 : e_occ = ev_occ(imo_occ)
1072 : !
1073 17632 : DO imo_virt = 1, nmo_virt_selected(ispin)
1074 12454 : istate = istate + 1
1075 16056 : e_virt_minus_occ(istate) = ev_virt(imo_virt) - e_occ
1076 : END DO
1077 : END DO
1078 :
1079 2980 : DEALLOCATE (ev_virt, ev_occ)
1080 : END DO
1081 : ELSE
1082 : ! Guess for spin-flip TDDFT
1083 112 : DO imo_occ = 1, nmo_occ_selected(1)
1084 : ! Here imo_occ enumerate alpha Occupied orbitals in inverse order (from the last to the first element)
1085 90 : i = gs_mos(1)%nmo_active - imo_occ + 1
1086 90 : k = gs_mos(1)%index_active(i)
1087 90 : e_occ = gs_mos(1)%evals_occ(k)
1088 :
1089 502 : DO imo_virt = 1, nmo_virt_selected(2)
1090 390 : istate = istate + 1
1091 480 : e_virt_minus_occ(istate) = gs_mos(2)%evals_virt(imo_virt) - e_occ
1092 : END DO
1093 : END DO
1094 : END IF
1095 :
1096 : IF (debug_this_module) THEN
1097 : CPASSERT(istate == nstates_selected)
1098 : END IF
1099 :
1100 1426 : CALL sort(e_virt_minus_occ, nstates_selected, inds)
1101 :
1102 : ! Labels and spin component for closed-shell
1103 1426 : IF (nspins == 1) THEN
1104 1232 : spin1 = 1
1105 1232 : spin2 = spin1
1106 1232 : spin_label1 = ' '
1107 1232 : spin_label2 = spin_label1
1108 : ! Labels and spin component for spin-flip excitations
1109 194 : ELSE IF (tddfpt_control%spinflip /= no_sf_tddfpt) THEN
1110 22 : spin1 = 1
1111 22 : spin2 = 2
1112 22 : spin_label1 = '(alp)'
1113 22 : spin_label2 = '(bet)'
1114 : END IF
1115 :
1116 1426 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1117 : ! Calculate maximum number of alpha excitations
1118 1404 : nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1)
1119 : ELSE
1120 : ! Calculate maximum number of spin-flip excitations
1121 22 : nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(2)
1122 : END IF
1123 1426 : IF (log_unit > 0) THEN
1124 713 : WRITE (log_unit, "(1X,A)") "", &
1125 713 : "-------------------------------------------------------------------------------", &
1126 713 : "- TDDFPT Initial Guess -", &
1127 1426 : "-------------------------------------------------------------------------------"
1128 713 : WRITE (log_unit, '(T11,A)') "State Occupied -> Virtual Excitation"
1129 713 : WRITE (log_unit, '(T11,A)') "number orbital orbital energy (eV)"
1130 713 : WRITE (log_unit, '(1X,79("-"))')
1131 : END IF
1132 :
1133 4278 : i = MAXVAL(nmo(:))
1134 5704 : ALLOCATE (reverse_index(i, nspins))
1135 1426 : reverse_index = 0
1136 3046 : DO ispin = 1, nspins
1137 11250 : DO i = 1, SIZE(gs_mos(ispin)%index_active)
1138 8204 : j = gs_mos(ispin)%index_active(i)
1139 9824 : reverse_index(j, ispin) = i
1140 : END DO
1141 : END DO
1142 :
1143 5080 : DO istate = 1, nstates
1144 5080 : IF (ASSOCIATED(evects(1, istate)%matrix_struct)) THEN
1145 : ! Initial guess vector read from restart file
1146 6 : IF (log_unit > 0) THEN
1147 : WRITE (log_unit, '(T7,I8,T28,A19,T60,F14.5)') &
1148 3 : istate, "*** restarted ***", evals(istate)*evolt
1149 : END IF
1150 : ELSE
1151 : ! New initial guess vector
1152 : !
1153 : ! Index of excited state - 1
1154 3648 : ind = inds(istate) - 1
1155 :
1156 : ! Labels and spin component for open-shell spin-conserving excitations
1157 3648 : IF ((nspins > 1) .AND. (tddfpt_control%spinflip == no_sf_tddfpt)) THEN
1158 598 : IF (ind < nstates_occ_virt_alpha) THEN
1159 256 : spin1 = 1
1160 256 : spin2 = 1
1161 256 : spin_label1 = '(alp)'
1162 256 : spin_label2 = '(alp)'
1163 : ELSE
1164 342 : ind = ind - nstates_occ_virt_alpha
1165 342 : spin1 = 2
1166 342 : spin2 = 2
1167 342 : spin_label1 = '(bet)'
1168 342 : spin_label2 = '(bet)'
1169 : END IF
1170 : END IF
1171 :
1172 : ! Recover index of occupied MO (imo_occ) and unoccupied MO (imo_virt)
1173 : ! associated to the excited state index (ind+1)
1174 3648 : i = ind/nmo_virt_selected(spin2) + 1
1175 3648 : j = nmo_occ_avail(spin1) - i + 1
1176 3648 : imo_occ = gs_mos(spin1)%index_active(j)
1177 3648 : imo_virt = MOD(ind, nmo_virt_selected(spin2)) + 1
1178 : ! Assign initial guess for excitation energy
1179 3648 : evals(istate) = e_virt_minus_occ(istate)
1180 :
1181 3648 : IF (log_unit > 0) THEN
1182 : WRITE (log_unit, '(T7,I8,T24,I8,T37,A5,T45,I8,T54,A5,T60,F14.5)') &
1183 1824 : istate, imo_occ, spin_label1, nmo(spin2) + imo_virt, spin_label2, e_virt_minus_occ(istate)*evolt
1184 : END IF
1185 :
1186 7894 : DO jspin = 1, SIZE(evects, 1)
1187 : ! .NOT. ASSOCIATED(evects(jspin, istate)%matrix_struct))
1188 4246 : CALL fm_pool_create_fm(fm_pool_ao_mo_active(jspin)%pool, evects(jspin, istate))
1189 4246 : CALL cp_fm_set_all(evects(jspin, istate), 0.0_dp)
1190 :
1191 7894 : IF (jspin == spin1) THEN
1192 : ! Half transform excitation vector to ao space:
1193 : ! evects_mi = c_ma*X_ai
1194 3648 : i = reverse_index(imo_occ, spin1)
1195 : CALL cp_fm_to_fm(gs_mos(spin2)%mos_virt, evects(spin1, istate), &
1196 3648 : ncol=1, source_start=imo_virt, target_start=i)
1197 : END IF
1198 : END DO
1199 : END IF
1200 : END DO
1201 :
1202 1426 : DEALLOCATE (reverse_index)
1203 :
1204 1426 : IF (log_unit > 0) THEN
1205 713 : WRITE (log_unit, '(/,T7,A,T50,I24)') 'Number of active states:', &
1206 1426 : tddfpt_total_number_of_states(tddfpt_control, gs_mos)
1207 : WRITE (log_unit, "(1X,A)") &
1208 713 : "-------------------------------------------------------------------------------"
1209 : END IF
1210 :
1211 1426 : DEALLOCATE (e_virt_minus_occ)
1212 1426 : DEALLOCATE (inds)
1213 :
1214 1426 : CALL timestop(handle)
1215 :
1216 2852 : END SUBROUTINE tddfpt_guess_vectors
1217 :
1218 : END MODULE qs_tddfpt2_utils
|
