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