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