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