(git:03cdb6f)
Loading...
Searching...
No Matches
qs_tddfpt2_utils.F
Go to the documentation of this file.
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
9 USE cell_types, ONLY: cell_type
14 USE cp_dbcsr_api, ONLY: dbcsr_add,&
31 USE cp_fm_types, ONLY: cp_fm_create,&
42 USE input_constants, ONLY: &
53 USE kinds, ONLY: dp,&
54 int_8
57 USE physcon, ONLY: evolt
61 USE qs_ks_types, ONLY: qs_ks_env_type,&
63 USE qs_mo_types, ONLY: allocate_mo_set,&
70 USE qs_scf_types, ONLY: ot_method_nr,&
73 USE util, ONLY: sort
74 USE xc_pot_saop, ONLY: add_saop_pot
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
92
93! **************************************************************************************************
94
95CONTAINS
96
97! **************************************************************************************************
98!> \brief Prepare MOs for TDDFPT Calculations
99!> \param qs_env Quickstep environment
100!> \param gs_mos ...
101!> \param iounit ...
102! **************************************************************************************************
103 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 REAL(kind=dp), DIMENSION(:), POINTER :: evals_virt_spin
116 TYPE(cell_type), POINTER :: cell
117 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 TARGET :: mos_virt
121 TYPE(cp_fm_type), POINTER :: mos_virt_spin
122 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 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 CALL timeset(routinen, handle)
132
133 CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, &
134 matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, scf_env=scf_env)
135 tddfpt_control => dft_control%tddfpt2_control
136 IF ((tddfpt_control%do_bse) .OR. (tddfpt_control%do_bse_w_only) .OR. &
137 (tddfpt_control%do_bse_gw_only)) THEN
138 NULLIFY (ks_env, ex_env)
139 CALL get_qs_env(qs_env, exstate_env=ex_env, ks_env=ks_env)
140 CALL dbcsr_copy(matrix_ks(1)%matrix, ex_env%matrix_ks(1)%matrix)
141 CALL set_ks_env(ks_env, matrix_ks=matrix_ks)
142 END IF
143
144 cpassert(.NOT. ASSOCIATED(gs_mos))
145 ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals
146 nspins = dft_control%nspins
147 ALLOCATE (gs_mos(nspins))
148
149 ! check if virtuals should be constructed for NAMD interface with NEWTONX
150 print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
151 CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals_newtonx)
152
153 ! when the number of unoccupied orbitals is limited and OT has been used
154 ! for the ground-state DFT calculation,
155 ! compute the missing unoccupied orbitals using OT as well.
156 NULLIFY (evals_virt, evals_virt_spin, mos_virt_spin)
157 IF (ASSOCIATED(scf_env)) THEN
158 IF ((scf_env%method == ot_method_nr .AND. tddfpt_control%nlumo > 0) .OR. &
159 (scf_env%method == ot_method_nr .AND. print_virtuals_newtonx)) THEN
160 ! As OT with ADDED_MOS/=0 is currently not implemented, the following block is equivalent to:
161 ! nmo_virt = tddfpt_control%nlumo
162 ! number of already computed unoccupied orbitals (added_mos) .
163 nmo_virt = huge(0)
164 DO ispin = 1, nspins
165 CALL get_mo_set(mos(ispin), nmo=nmo_avail, homo=nmo_occ)
166 nmo_virt = min(nmo_virt, nmo_avail - nmo_occ)
167 END DO
168 ! number of unoccupied orbitals to compute
169 nmo_virt = tddfpt_control%nlumo - nmo_virt
170 IF (.NOT. print_virtuals_newtonx) THEN
171 IF (nmo_virt > 0) THEN
172 ALLOCATE (evals_virt(nspins), mos_virt(nspins))
173 ! the number of actually computed unoccupied orbitals will be stored as 'nmo_avail'
174 CALL make_lumo_gpw(qs_env, scf_env, mos_virt, evals_virt, nmo_virt, nmo_avail)
175 END IF
176 END IF
177 END IF
178 END IF
179
180 DO ispin = 1, nspins
181 IF (ASSOCIATED(evals_virt)) THEN
182 evals_virt_spin => evals_virt(ispin)%array
183 ELSE
184 NULLIFY (evals_virt_spin)
185 END IF
186 IF (ALLOCATED(mos_virt)) THEN
187 mos_virt_spin => mos_virt(ispin)
188 ELSE
189 NULLIFY (mos_virt_spin)
190 END IF
191 CALL tddfpt_init_ground_state_mos(gs_mos=gs_mos(ispin), mo_set=mos(ispin), &
192 nlumo=tddfpt_control%nlumo, &
193 blacs_env=blacs_env, cholesky_method=cholesky_restore, &
194 matrix_ks=matrix_ks(ispin)%matrix, matrix_s=matrix_s(1)%matrix, &
195 mos_virt=mos_virt_spin, evals_virt=evals_virt_spin, &
196 qs_env=qs_env)
197 END DO
198
199 moc = 0
200 mvt = 0
201 DO ispin = 1, nspins
202 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, nrow_global=moc(1, ispin), ncol_global=moc(2, ispin))
203 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, nrow_global=mvt(1, ispin), ncol_global=mvt(2, ispin))
204 END DO
205 IF (iounit > 0) THEN
206 WRITE (iounit, "(T2,A,T36,A)") "TDDFPT| Molecular Orbitals:", &
207 " Spin AOs Occ Virt Total"
208 DO ispin = 1, nspins
209 WRITE (iounit, "(T2,A,T37,I4,4I10)") "TDDFPT| ", ispin, moc(1, ispin), moc(2, ispin), &
210 mvt(2, ispin), moc(2, ispin) + mvt(2, ispin)
211 END DO
212 END IF
213
214 IF (ASSOCIATED(evals_virt)) THEN
215 DO ispin = 1, SIZE(evals_virt)
216 IF (ASSOCIATED(evals_virt(ispin)%array)) DEALLOCATE (evals_virt(ispin)%array)
217 END DO
218 DEALLOCATE (evals_virt)
219 END IF
220
221 CALL cp_fm_release(mos_virt)
222
223 CALL timestop(handle)
224
225 END SUBROUTINE tddfpt_init_mos
226
227! **************************************************************************************************
228!> \brief Generate all virtual molecular orbitals for a given spin by diagonalising
229!> the corresponding Kohn-Sham matrix.
230!> \param gs_mos structure to store occupied and virtual molecular orbitals
231!> (allocated and initialised on exit)
232!> \param mo_set ground state molecular orbitals for a given spin
233!> \param nlumo number of unoccupied states to consider (-1 means all states)
234!> \param blacs_env BLACS parallel environment
235!> \param cholesky_method Cholesky method to compute the inverse overlap matrix
236!> \param matrix_ks Kohn-Sham matrix for a given spin
237!> \param matrix_s overlap matrix
238!> \param mos_virt precomputed (OT) expansion coefficients of virtual molecular orbitals
239!> (in addition to the ADDED_MOS, if present). NULL when no OT is in use.
240!> \param evals_virt orbital energies of precomputed (OT) virtual molecular orbitals.
241!> NULL when no OT is in use.
242!> \param qs_env ...
243!> \par History
244!> * 05.2016 created as tddfpt_lumos() [Sergey Chulkov]
245!> * 06.2016 renamed, altered prototype [Sergey Chulkov]
246!> * 04.2019 limit the number of unoccupied states, orbital energy correction [Sergey Chulkov]
247! **************************************************************************************************
248 SUBROUTINE tddfpt_init_ground_state_mos(gs_mos, mo_set, nlumo, blacs_env, cholesky_method, matrix_ks, matrix_s, &
249 mos_virt, evals_virt, qs_env)
250 TYPE(tddfpt_ground_state_mos) :: gs_mos
251 TYPE(mo_set_type), INTENT(IN) :: mo_set
252 INTEGER, INTENT(in) :: nlumo
253 TYPE(cp_blacs_env_type), POINTER :: blacs_env
254 INTEGER, INTENT(in) :: cholesky_method
255 TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_s
256 TYPE(cp_fm_type), INTENT(IN), POINTER :: mos_virt
257 REAL(kind=dp), DIMENSION(:), POINTER :: evals_virt
258 TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
259
260 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_init_ground_state_mos'
261 REAL(kind=dp), PARAMETER :: eps_dp = epsilon(0.0_dp)
262
263 INTEGER :: cholesky_method_inout, handle, icol_global, icol_local, imo, iounit, irow_global, &
264 irow_local, nao, ncol_local, nelectrons, nmo_occ, nmo_scf, nmo_virt, nrow_local, sign_int
265 INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, &
266 sum_sign_array
267 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
268 LOGICAL :: do_eigen, print_phases
269 REAL(kind=dp) :: element, maxocc
270 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
271 POINTER :: my_block
272 REAL(kind=dp), DIMENSION(:), POINTER :: mo_evals_extended, mo_occ_extended, &
273 mo_occ_scf
274 TYPE(cp_fm_struct_type), POINTER :: ao_ao_fm_struct, ao_mo_occ_fm_struct, &
275 ao_mo_virt_fm_struct, wfn_fm_struct
276 TYPE(cp_fm_type) :: matrix_ks_fm, ortho_fm, work_fm, &
277 work_fm_virt
278 TYPE(cp_fm_type), POINTER :: mo_coeff_extended
279 TYPE(cp_logger_type), POINTER :: logger
280 TYPE(mo_set_type), POINTER :: mos_extended
281 TYPE(mp_para_env_type), POINTER :: para_env
282 TYPE(section_vals_type), POINTER :: print_section
283
284 CALL timeset(routinen, handle)
285
286 NULLIFY (logger)
287 logger => cp_get_default_logger()
288 iounit = cp_logger_get_default_io_unit(logger)
289
290 CALL blacs_env%get(para_env=para_env)
291
292 CALL get_mo_set(mo_set, nao=nao, nmo=nmo_scf, homo=nmo_occ, maxocc=maxocc, &
293 nelectron=nelectrons, occupation_numbers=mo_occ_scf)
294
295 print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
296 CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases)
297
298 nmo_virt = nao - nmo_occ
299 IF (nlumo >= 0) &
300 nmo_virt = min(nmo_virt, nlumo)
301
302 IF (nmo_virt <= 0) &
303 CALL cp_abort(__location__, &
304 'At least one unoccupied molecular orbital is required to calculate excited states.')
305
306 do_eigen = .false.
307 ! diagonalise the Kohn-Sham matrix one more time if the number of available unoccupied states are too small
308 IF (ASSOCIATED(evals_virt)) THEN
309 cpassert(ASSOCIATED(mos_virt))
310 IF (nmo_virt > nmo_scf - nmo_occ + SIZE(evals_virt)) do_eigen = .true.
311 ELSE
312 IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .true.
313 END IF
314
315 ! ++ allocate storage space for gs_mos
316 NULLIFY (ao_mo_occ_fm_struct, ao_mo_virt_fm_struct)
317 ! Tiny fix (A.Sinyavskiy)
318 CALL cp_fm_struct_create(ao_mo_occ_fm_struct, template_fmstruct=mo_set%mo_coeff%matrix_struct, &
319 ncol_global=nmo_occ, context=blacs_env)
320 CALL cp_fm_struct_create(ao_mo_virt_fm_struct, template_fmstruct=mo_set%mo_coeff%matrix_struct, &
321 ncol_global=nmo_virt, context=blacs_env)
322
323 NULLIFY (gs_mos%mos_occ, gs_mos%mos_virt, gs_mos%evals_occ_matrix)
324 ALLOCATE (gs_mos%mos_occ, gs_mos%mos_virt)
325 CALL cp_fm_create(gs_mos%mos_occ, ao_mo_occ_fm_struct)
326 CALL cp_fm_create(gs_mos%mos_virt, ao_mo_virt_fm_struct)
327 gs_mos%nmo_occ = nmo_occ
328
329 ALLOCATE (gs_mos%evals_occ(nmo_occ))
330 ALLOCATE (gs_mos%evals_virt(nmo_virt))
331 ALLOCATE (gs_mos%phases_occ(nmo_occ))
332 ALLOCATE (gs_mos%phases_virt(nmo_virt))
333
334 ! ++ nullify pointers
335 NULLIFY (ao_ao_fm_struct, wfn_fm_struct)
336 NULLIFY (mos_extended, mo_coeff_extended, mo_evals_extended, mo_occ_extended)
337
338 IF (do_eigen) THEN
339 ! ++ set of molecular orbitals
340 CALL cp_fm_struct_create(ao_ao_fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
341 CALL cp_fm_struct_create(wfn_fm_struct, nrow_global=nao, ncol_global=nmo_occ + nmo_virt, context=blacs_env)
342 ALLOCATE (mos_extended)
343 CALL allocate_mo_set(mos_extended, nao, nmo_occ + nmo_virt, nelectrons, &
344 REAL(nelectrons, dp), maxocc, flexible_electron_count=0.0_dp)
345 CALL init_mo_set(mos_extended, fm_struct=wfn_fm_struct, name="mos-extended")
346 CALL cp_fm_struct_release(wfn_fm_struct)
347 CALL get_mo_set(mos_extended, mo_coeff=mo_coeff_extended, &
348 eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
349
350 ! use the explicit loop in order to avoid temporary arrays.
351 !
352 ! The assignment statement : mo_occ_extended(1:nmo_scf) = mo_occ_scf(1:nmo_scf)
353 ! implies temporary arrays as a compiler does not know in advance that the pointers
354 ! on both sides of the statement point to non-overlapped memory regions
355 DO imo = 1, nmo_scf
356 mo_occ_extended(imo) = mo_occ_scf(imo)
357 END DO
358 mo_occ_extended(nmo_scf + 1:) = 0.0_dp
359
360 ! ++ allocate temporary matrices
361 CALL cp_fm_create(matrix_ks_fm, ao_ao_fm_struct)
362 CALL cp_fm_create(ortho_fm, ao_ao_fm_struct)
363 CALL cp_fm_create(work_fm, ao_ao_fm_struct)
364 CALL cp_fm_struct_release(ao_ao_fm_struct)
365
366 ! some stuff from the subroutine general_eigenproblem()
367 CALL copy_dbcsr_to_fm(matrix_s, ortho_fm)
368 CALL copy_dbcsr_to_fm(matrix_ks, matrix_ks_fm)
369
370 IF (cholesky_method == cholesky_dbcsr) THEN
371 cpabort('CHOLESKY DBCSR_INVERSE is not implemented in TDDFT.')
372 ELSE IF (cholesky_method == cholesky_off) THEN
373 cpabort('CHOLESKY OFF is not implemented in TDDFT.')
374 ELSE
375 CALL cp_fm_cholesky_decompose(ortho_fm)
376 IF (cholesky_method == cholesky_inverse) THEN
377 CALL cp_fm_triangular_invert(ortho_fm)
378 END IF
379
380 ! need to store 'cholesky_method' in a temporary variable, as the subroutine eigensolver()
381 ! will update this variable
382 cholesky_method_inout = cholesky_method
383 CALL eigensolver(matrix_ks_fm=matrix_ks_fm, mo_set=mos_extended, ortho=ortho_fm, &
384 work=work_fm, cholesky_method=cholesky_method_inout, &
385 do_level_shift=.false., level_shift=0.0_dp, use_jacobi=.false.)
386 END IF
387
388 ! -- clean up needless matrices
389 CALL cp_fm_release(work_fm)
390 CALL cp_fm_release(ortho_fm)
391 CALL cp_fm_release(matrix_ks_fm)
392 ELSE
393 CALL get_mo_set(mo_set, mo_coeff=mo_coeff_extended, &
394 eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
395 END IF
396
397 ! compute the phase of molecular orbitals;
398 ! matrix work_fm holds occupied molecular orbital coefficients distributed among all the processors
399 !CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
400 CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
401 CALL cp_fm_struct_release(ao_mo_occ_fm_struct)
402
403 CALL cp_fm_to_fm(mo_coeff_extended, work_fm, ncol=nmo_occ, source_start=1, target_start=1)
404 CALL cp_fm_get_info(work_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
405 row_indices=row_indices, col_indices=col_indices, local_data=my_block)
406
407 ALLOCATE (minrow_neg_array(nmo_occ), minrow_pos_array(nmo_occ), sum_sign_array(nmo_occ))
408 minrow_neg_array(:) = nao
409 minrow_pos_array(:) = nao
410 sum_sign_array(:) = 0
411 DO icol_local = 1, ncol_local
412 icol_global = col_indices(icol_local)
413
414 DO irow_local = 1, nrow_local
415 element = my_block(irow_local, icol_local)
416
417 sign_int = 0
418 IF (element >= eps_dp) THEN
419 sign_int = 1
420 ELSE IF (element <= -eps_dp) THEN
421 sign_int = -1
422 END IF
423
424 sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
425
426 irow_global = row_indices(irow_local)
427 IF (sign_int > 0) THEN
428 IF (minrow_pos_array(icol_global) > irow_global) &
429 minrow_pos_array(icol_global) = irow_global
430 ELSE IF (sign_int < 0) THEN
431 IF (minrow_neg_array(icol_global) > irow_global) &
432 minrow_neg_array(icol_global) = irow_global
433 END IF
434 END DO
435 END DO
436
437 CALL para_env%sum(sum_sign_array)
438 CALL para_env%min(minrow_neg_array)
439 CALL para_env%min(minrow_pos_array)
440
441 DO icol_local = 1, nmo_occ
442 IF (sum_sign_array(icol_local) > 0) THEN
443 ! most of the expansion coefficients are positive => MO's phase = +1
444 gs_mos%phases_occ(icol_local) = 1.0_dp
445 ELSE IF (sum_sign_array(icol_local) < 0) THEN
446 ! most of the expansion coefficients are negative => MO's phase = -1
447 gs_mos%phases_occ(icol_local) = -1.0_dp
448 ELSE
449 ! equal number of positive and negative expansion coefficients
450 IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local)) THEN
451 ! the first positive expansion coefficient has a lower index then
452 ! the first negative expansion coefficient; MO's phase = +1
453 gs_mos%phases_occ(icol_local) = 1.0_dp
454 ELSE
455 ! MO's phase = -1
456 gs_mos%phases_occ(icol_local) = -1.0_dp
457 END IF
458 END IF
459 END DO
460
461 DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
462
463 ! return the requested occupied and virtual molecular orbitals and corresponding orbital energies
464 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_occ, ncol=nmo_occ, source_start=1, target_start=1)
465 gs_mos%evals_occ(1:nmo_occ) = mo_evals_extended(1:nmo_occ)
466
467 IF (ASSOCIATED(evals_virt) .AND. (.NOT. do_eigen) .AND. nmo_virt > nmo_scf - nmo_occ) THEN
468 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_scf - nmo_occ, &
469 source_start=nmo_occ + 1, target_start=1)
470 CALL cp_fm_to_fm(mos_virt, gs_mos%mos_virt, ncol=nmo_virt - (nmo_scf - nmo_occ), &
471 source_start=1, target_start=nmo_scf - nmo_occ + 1)
472 gs_mos%evals_virt(1:nmo_scf - nmo_occ) = evals_virt(nmo_occ + 1:nmo_occ + nmo_scf)
473 gs_mos%evals_virt(nmo_scf - nmo_occ + 1:nmo_virt) = evals_virt(1:nmo_virt - (nmo_scf - nmo_occ))
474 ELSE
475 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_virt, source_start=nmo_occ + 1, target_start=1)
476 gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt)
477 END IF
478
479 IF (print_phases) THEN
480 ! compute the phase of molecular orbitals;
481 ! matrix work_fm holds virtual molecular orbital coefficients distributed among all the processors
482 !CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
483 CALL cp_fm_create(work_fm_virt, ao_mo_virt_fm_struct)
484
485 CALL cp_fm_to_fm(gs_mos%mos_virt, work_fm_virt, ncol=nmo_virt, source_start=1, target_start=1)
486 CALL cp_fm_get_info(work_fm_virt, nrow_local=nrow_local, ncol_local=ncol_local, &
487 row_indices=row_indices, col_indices=col_indices, local_data=my_block)
488
489 ALLOCATE (minrow_neg_array(nmo_virt), minrow_pos_array(nmo_virt), sum_sign_array(nmo_virt))
490 minrow_neg_array(:) = nao
491 minrow_pos_array(:) = nao
492 sum_sign_array(:) = 0
493 DO icol_local = 1, ncol_local
494 icol_global = col_indices(icol_local)
495
496 DO irow_local = 1, nrow_local
497 element = my_block(irow_local, icol_local)
498
499 sign_int = 0
500 IF (element >= eps_dp) THEN
501 sign_int = 1
502 ELSE IF (element <= -eps_dp) THEN
503 sign_int = -1
504 END IF
505
506 sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
507
508 irow_global = row_indices(irow_local)
509 IF (sign_int > 0) THEN
510 IF (minrow_pos_array(icol_global) > irow_global) &
511 minrow_pos_array(icol_global) = irow_global
512 ELSE IF (sign_int < 0) THEN
513 IF (minrow_neg_array(icol_global) > irow_global) &
514 minrow_neg_array(icol_global) = irow_global
515 END IF
516 END DO
517 END DO
518
519 CALL para_env%sum(sum_sign_array)
520 CALL para_env%min(minrow_neg_array)
521 CALL para_env%min(minrow_pos_array)
522 DO icol_local = 1, nmo_virt
523 IF (sum_sign_array(icol_local) > 0) THEN
524 ! most of the expansion coefficients are positive => MO's phase = +1
525 gs_mos%phases_virt(icol_local) = 1.0_dp
526 ELSE IF (sum_sign_array(icol_local) < 0) THEN
527 ! most of the expansion coefficients are negative => MO's phase = -1
528 gs_mos%phases_virt(icol_local) = -1.0_dp
529 ELSE
530 ! equal number of positive and negative expansion coefficients
531 IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local)) THEN
532 ! the first positive expansion coefficient has a lower index then
533 ! the first negative expansion coefficient; MO's phase = +1
534 gs_mos%phases_virt(icol_local) = 1.0_dp
535 ELSE
536 ! MO's phase = -1
537 gs_mos%phases_virt(icol_local) = -1.0_dp
538 END IF
539 END IF
540 END DO
541 DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
542 CALL cp_fm_release(work_fm_virt)
543 END IF !print_phases
544 CALL cp_fm_struct_release(ao_mo_virt_fm_struct) ! here after print_phases
545
546 CALL cp_fm_release(work_fm)
547
548 IF (do_eigen) THEN
549 CALL deallocate_mo_set(mos_extended)
550 DEALLOCATE (mos_extended)
551 END IF
552
553 CALL timestop(handle)
554
555 END SUBROUTINE tddfpt_init_ground_state_mos
556
557! **************************************************************************************************
558!> \brief Release molecular orbitals.
559!> \param gs_mos structure that holds occupied and virtual molecular orbitals
560!> \par History
561!> * 06.2016 created [Sergey Chulkov]
562! **************************************************************************************************
564 TYPE(tddfpt_ground_state_mos), INTENT(inout) :: gs_mos
565
566 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_release_ground_state_mos'
567
568 INTEGER :: handle
569
570 CALL timeset(routinen, handle)
571
572 IF (ALLOCATED(gs_mos%phases_occ)) &
573 DEALLOCATE (gs_mos%phases_occ)
574
575 IF (ALLOCATED(gs_mos%evals_virt)) &
576 DEALLOCATE (gs_mos%evals_virt)
577
578 IF (ALLOCATED(gs_mos%evals_occ)) &
579 DEALLOCATE (gs_mos%evals_occ)
580
581 IF (ALLOCATED(gs_mos%phases_virt)) &
582 DEALLOCATE (gs_mos%phases_virt)
583
584 IF (ALLOCATED(gs_mos%index_active)) &
585 DEALLOCATE (gs_mos%index_active)
586
587 IF (ASSOCIATED(gs_mos%evals_occ_matrix)) THEN
588 CALL cp_fm_release(gs_mos%evals_occ_matrix)
589 DEALLOCATE (gs_mos%evals_occ_matrix)
590 END IF
591
592 IF (ASSOCIATED(gs_mos%mos_virt)) THEN
593 CALL cp_fm_release(gs_mos%mos_virt)
594 DEALLOCATE (gs_mos%mos_virt)
595 END IF
596
597 IF (ASSOCIATED(gs_mos%mos_occ)) THEN
598 CALL cp_fm_release(gs_mos%mos_occ)
599 DEALLOCATE (gs_mos%mos_occ)
600 END IF
601
602 IF (ASSOCIATED(gs_mos%mos_active)) THEN
603 CALL cp_fm_release(gs_mos%mos_active)
604 DEALLOCATE (gs_mos%mos_active)
605 END IF
606
607 CALL timestop(handle)
608
610
611! **************************************************************************************************
612!> \brief Callculate orbital corrected KS matrix for TDDFPT
613!> \param qs_env Quickstep environment
614!> \param gs_mos ...
615!> \param matrix_ks_oep ...
616! **************************************************************************************************
617 SUBROUTINE tddfpt_oecorr(qs_env, gs_mos, matrix_ks_oep)
618 TYPE(qs_environment_type), POINTER :: qs_env
619 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
620 POINTER :: gs_mos
621 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_oep
622
623 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_oecorr'
624
625 INTEGER :: handle, iounit, ispin, nao, nmo_occ, &
626 nspins
627 LOGICAL :: do_hfx
628 TYPE(cp_blacs_env_type), POINTER :: blacs_env
629 TYPE(cp_fm_struct_type), POINTER :: ao_mo_occ_fm_struct, &
630 mo_occ_mo_occ_fm_struct
631 TYPE(cp_fm_type) :: work_fm
632 TYPE(cp_logger_type), POINTER :: logger
633 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
634 TYPE(dft_control_type), POINTER :: dft_control
635 TYPE(section_vals_type), POINTER :: hfx_section, xc_fun_empty, &
636 xc_fun_original
637 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
638
639 CALL timeset(routinen, handle)
640
641 NULLIFY (logger)
642 logger => cp_get_default_logger()
643 iounit = cp_logger_get_default_io_unit(logger)
644
645 CALL get_qs_env(qs_env, blacs_env=blacs_env, dft_control=dft_control, matrix_ks=matrix_ks)
646 tddfpt_control => dft_control%tddfpt2_control
647
648 ! obtain corrected KS-matrix
649 ! We should 'save' the energy values?
650 nspins = SIZE(gs_mos)
651 NULLIFY (matrix_ks_oep)
652 IF (tddfpt_control%oe_corr /= oe_none) THEN
653 IF (iounit > 0) THEN
654 WRITE (iounit, "(1X,A)") "", &
655 "-------------------------------------------------------------------------------", &
656 "- Orbital Eigenvalue Correction Started -", &
657 "-------------------------------------------------------------------------------"
658 END IF
659
660 CALL cp_warn(__location__, &
661 "Orbital energy correction potential is an experimental feature. "// &
662 "Use it with extreme care")
663
664 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
665 CALL section_vals_get(hfx_section, explicit=do_hfx)
666 IF (do_hfx) THEN
667 CALL cp_abort(__location__, &
668 "Implementation of orbital energy correction XC-potentials is "// &
669 "currently incompatible with exact-exchange functionals")
670 END IF
671
672 CALL dbcsr_allocate_matrix_set(matrix_ks_oep, nspins)
673 DO ispin = 1, nspins
674 CALL dbcsr_init_p(matrix_ks_oep(ispin)%matrix)
675 CALL dbcsr_copy(matrix_ks_oep(ispin)%matrix, matrix_ks(ispin)%matrix)
676 END DO
677
678 ! KS-matrix without XC-terms
679 xc_fun_original => section_vals_get_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL")
680 CALL section_vals_retain(xc_fun_original)
681 NULLIFY (xc_fun_empty)
682 CALL section_vals_create(xc_fun_empty, xc_fun_original%section)
683 CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_empty)
684 CALL section_vals_release(xc_fun_empty)
685
686 IF (dft_control%qs_control%semi_empirical) THEN
687 cpabort("TDDFPT with SE not possible")
688 ELSEIF (dft_control%qs_control%dftb) THEN
689 cpabort("TDDFPT with DFTB not possible")
690 ELSEIF (dft_control%qs_control%xtb) THEN
691 CALL build_xtb_ks_matrix(qs_env, calculate_forces=.false., just_energy=.false., &
692 ext_ks_matrix=matrix_ks_oep)
693 ELSE
694 CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.false., just_energy=.false., &
695 ext_ks_matrix=matrix_ks_oep)
696 END IF
697
698 IF (tddfpt_control%oe_corr == oe_saop .OR. &
699 tddfpt_control%oe_corr == oe_lb .OR. &
700 tddfpt_control%oe_corr == oe_gllb) THEN
701 IF (iounit > 0) THEN
702 WRITE (iounit, "(T2,A)") " Orbital energy correction of SAOP type "
703 END IF
704 CALL add_saop_pot(matrix_ks_oep, qs_env, tddfpt_control%oe_corr)
705 ELSE IF (tddfpt_control%oe_corr == oe_shift) THEN
706 IF (iounit > 0) THEN
707 WRITE (iounit, "(T2,A,T71,F10.3)") &
708 " Virtual Orbital Eigenvalue Shift [eV] ", tddfpt_control%ev_shift*evolt
709 WRITE (iounit, "(T2,A,T71,F10.3)") &
710 " Open Shell Orbital Eigenvalue Shift [eV] ", tddfpt_control%eos_shift*evolt
711 END IF
712 CALL ev_shift_operator(qs_env, gs_mos, matrix_ks_oep, &
713 tddfpt_control%ev_shift, tddfpt_control%eos_shift)
714 ELSE
715 CALL cp_abort(__location__, &
716 "Unimplemented orbital energy correction potential")
717 END IF
718 CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_original)
719 CALL section_vals_release(xc_fun_original)
720
721 ! compute 'evals_occ_matrix'
722 CALL dbcsr_get_info(matrix_ks(1)%matrix, nfullrows_total=nao)
723 NULLIFY (mo_occ_mo_occ_fm_struct)
724 DO ispin = 1, nspins
725 nmo_occ = SIZE(gs_mos(ispin)%evals_occ)
726 CALL cp_fm_struct_create(mo_occ_mo_occ_fm_struct, nrow_global=nmo_occ, ncol_global=nmo_occ, &
727 context=blacs_env)
728 ALLOCATE (gs_mos(ispin)%evals_occ_matrix)
729 CALL cp_fm_create(gs_mos(ispin)%evals_occ_matrix, mo_occ_mo_occ_fm_struct)
730 CALL cp_fm_struct_release(mo_occ_mo_occ_fm_struct)
731 ! work_fm is a temporary [nao x nmo_occ] matrix
732 CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, &
733 context=blacs_env)
734 CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
735 CALL cp_fm_struct_release(ao_mo_occ_fm_struct)
736 CALL cp_dbcsr_sm_fm_multiply(matrix_ks_oep(ispin)%matrix, gs_mos(ispin)%mos_occ, &
737 work_fm, ncol=nmo_occ, alpha=1.0_dp, beta=0.0_dp)
738 CALL parallel_gemm('T', 'N', nmo_occ, nmo_occ, nao, 1.0_dp, gs_mos(ispin)%mos_occ, work_fm, &
739 0.0_dp, gs_mos(ispin)%evals_occ_matrix)
740 CALL cp_fm_release(work_fm)
741 END DO
742 IF (iounit > 0) THEN
743 WRITE (iounit, "(1X,A)") &
744 "-------------------------------------------------------------------------------"
745 END IF
746
747 END IF
748
749 CALL timestop(handle)
750
751 END SUBROUTINE tddfpt_oecorr
752
753! **************************************************************************************************
754!> \brief Compute the number of possible singly excited states (occ -> virt)
755!> \param tddfpt_control ...
756!> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state
757!> \return the number of possible single excitations
758!> \par History
759!> * 01.2017 created [Sergey Chulkov]
760! **************************************************************************************************
761 PURE FUNCTION tddfpt_total_number_of_states(tddfpt_control, gs_mos) RESULT(nstates_total)
762 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
763 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
764 INTENT(in) :: gs_mos
765 INTEGER(kind=int_8) :: nstates_total
766
767 INTEGER :: ispin, nspins
768
769 nstates_total = 0
770 nspins = SIZE(gs_mos)
771
772 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
773 ! Total number of possible excitations for spin-conserving TDDFT
774 DO ispin = 1, nspins
775 nstates_total = nstates_total + &
776 gs_mos(ispin)%nmo_active* &
777 SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
778 END DO
779 ELSE
780 ! Total number of possible excitations for spin-flip TDDFT
781 nstates_total = gs_mos(1)%nmo_active* &
782 SIZE(gs_mos(2)%evals_virt, kind=int_8)
783 END IF
785
786! **************************************************************************************************
787!> \brief Create a shift operator on virtual/open shell space
788!> Shift operator = Edelta*Q Q: projector on virtual space (1-PS)
789!> projector on open shell space PosS
790!> \param qs_env the qs_env that is perturbed by this p_env
791!> \param gs_mos ...
792!> \param matrix_ks ...
793!> \param ev_shift ...
794!> \param eos_shift ...
795!> \par History
796!> 02.04.2019 adapted for TDDFT use from p_env (JGH)
797!> \author JGH
798! **************************************************************************************************
799 SUBROUTINE ev_shift_operator(qs_env, gs_mos, matrix_ks, ev_shift, eos_shift)
800
801 TYPE(qs_environment_type), POINTER :: qs_env
802 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
803 POINTER :: gs_mos
804 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
805 REAL(kind=dp), INTENT(IN) :: ev_shift, eos_shift
806
807 CHARACTER(len=*), PARAMETER :: routinen = 'ev_shift_operator'
808
809 INTEGER :: handle, ispin, n_spins, na, nb, nhomo, &
810 nl, nos, nrow, nu, nvirt
811 TYPE(cp_fm_struct_type), POINTER :: fmstruct
812 TYPE(cp_fm_type) :: cmos, cvec
813 TYPE(cp_fm_type), POINTER :: coeff
814 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
815 TYPE(dbcsr_type), POINTER :: smat
816 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
817
818 CALL timeset(routinen, handle)
819
820 n_spins = SIZE(gs_mos)
821 cpassert(n_spins == SIZE(matrix_ks))
822
823 IF (eos_shift /= 0.0_dp .AND. n_spins > 1) THEN
824 cpabort("eos_shift not implemented")
825 CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
826 smat => matrix_s(1)%matrix
827 CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
828 CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
829 nl = min(na, nb)
830 nu = max(na, nb)
831 ! open shell orbital shift
832 DO ispin = 1, n_spins
833 coeff => gs_mos(ispin)%mos_occ
834 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
835 IF (nhomo == nu) THEN
836 ! downshift with -eos_shift using occupied orbitals
837 nos = nu - nl
838 CALL cp_fm_create(cmos, fmstruct)
839 CALL cp_fm_get_info(coeff, nrow_global=nrow)
840 CALL cp_fm_to_fm_submat(coeff, cmos, nrow, nos, 1, nl + 1, 1, 1)
841 CALL cp_fm_create(cvec, fmstruct)
842 CALL cp_dbcsr_sm_fm_multiply(smat, cmos, cvec, nos, 1.0_dp, 0.0_dp)
843 CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
844 alpha=-eos_shift, keep_sparsity=.true.)
845 CALL cp_fm_release(cmos)
846 CALL cp_fm_release(cvec)
847 ELSE
848 ! upshift with eos_shift using virtual orbitals
849 coeff => gs_mos(ispin)%mos_virt
850 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
851 nos = nu - nhomo
852 cpassert(nvirt >= nos)
853 CALL cp_fm_create(cvec, fmstruct)
854 CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nos, 1.0_dp, 0.0_dp)
855 CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
856 alpha=eos_shift, keep_sparsity=.true.)
857 CALL cp_fm_release(cvec)
858 END IF
859 END DO
860 ! virtual shift
861 IF (ev_shift /= 0.0_dp) THEN
862 DO ispin = 1, n_spins
863 CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
864 alpha_scalar=1.0_dp, beta_scalar=ev_shift)
865 coeff => gs_mos(ispin)%mos_occ
866 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
867 CALL cp_fm_create(cvec, fmstruct)
868 CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nhomo, 1.0_dp, 0.0_dp)
869 CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nhomo, &
870 alpha=-ev_shift, keep_sparsity=.true.)
871 CALL cp_fm_release(cvec)
872 IF (nhomo < nu) THEN
873 nos = nu - nhomo
874 coeff => gs_mos(ispin)%mos_virt
875 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
876 cpassert(nvirt >= nos)
877 CALL cp_fm_create(cvec, fmstruct)
878 CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nos, 1.0_dp, 0.0_dp)
879 CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
880 alpha=-ev_shift, keep_sparsity=.true.)
881 CALL cp_fm_release(cvec)
882 END IF
883 END DO
884 END IF
885 ELSE
886 ! virtual shift
887 IF (ev_shift /= 0.0_dp) THEN
888 CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
889 smat => matrix_s(1)%matrix
890 DO ispin = 1, n_spins
891 CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
892 alpha_scalar=1.0_dp, beta_scalar=ev_shift)
893 coeff => gs_mos(ispin)%mos_occ
894 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
895 CALL cp_fm_create(cvec, fmstruct)
896 CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nhomo, 1.0_dp, 0.0_dp)
897 CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nhomo, &
898 alpha=-ev_shift, keep_sparsity=.true.)
899 CALL cp_fm_release(cvec)
900 END DO
901 END IF
902 END IF
903 ! set eigenvalues
904 IF (eos_shift == 0.0_dp .OR. n_spins == 1) THEN
905 DO ispin = 1, n_spins
906 IF (ALLOCATED(gs_mos(ispin)%evals_virt)) THEN
907 gs_mos(ispin)%evals_virt(:) = gs_mos(ispin)%evals_virt(:) + ev_shift
908 END IF
909 END DO
910 ELSE
911 CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
912 CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
913 nl = min(na, nb)
914 nu = max(na, nb)
915 nos = nu - nl
916 IF (na == nu) THEN
917 IF (ALLOCATED(gs_mos(1)%evals_occ)) THEN
918 gs_mos(1)%evals_occ(nl + 1:nu) = gs_mos(1)%evals_occ(nl + 1:nu) - eos_shift
919 END IF
920 IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
921 gs_mos(1)%evals_virt(:) = gs_mos(1)%evals_virt(:) + ev_shift
922 END IF
923 IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
924 gs_mos(2)%evals_virt(1:nos) = gs_mos(2)%evals_virt(1:nos) + eos_shift
925 gs_mos(2)%evals_virt(nos + 1:) = gs_mos(2)%evals_virt(nos + 1:) + ev_shift
926 END IF
927 ELSE
928 IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
929 gs_mos(1)%evals_virt(1:nos) = gs_mos(1)%evals_virt(1:nos) + eos_shift
930 gs_mos(1)%evals_virt(nos + 1:) = gs_mos(1)%evals_virt(nos + 1:) + ev_shift
931 END IF
932 IF (ALLOCATED(gs_mos(2)%evals_occ)) THEN
933 gs_mos(2)%evals_occ(nl + 1:nu) = gs_mos(2)%evals_occ(nl + 1:nu) - eos_shift
934 END IF
935 IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
936 gs_mos(2)%evals_virt(:) = gs_mos(2)%evals_virt(:) + ev_shift
937 END IF
938 END IF
939 END IF
940
941 CALL timestop(handle)
942
943 END SUBROUTINE ev_shift_operator
944
945! **************************************************************************************************
946!> \brief Generate missed guess vectors.
947!> \param evects guess vectors distributed across all processors (initialised on exit)
948!> \param evals guessed transition energies (initialised on exit)
949!> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state
950!> \param log_unit output unit
951!> \param tddfpt_control ...
952!> \param fm_pool_ao_mo_active ...
953!> \param qs_env ...
954!> \param nspins ...
955!> \par History
956!> * 05.2016 created as tddfpt_guess() [Sergey Chulkov]
957!> * 06.2016 renamed, altered prototype, supports spin-polarised density [Sergey Chulkov]
958!> * 01.2017 simplified prototype, do not compute all possible singly-excited states
959!> [Sergey Chulkov]
960!> \note \parblock
961!> Based on the subroutine co_initial_guess() which was originally created by
962!> Thomas Chassaing on 06.2003.
963!>
964!> Only not associated guess vectors 'evects(spin, state)%matrix' are allocated and
965!> initialised; associated vectors assumed to be initialised elsewhere (e.g. using
966!> a restart file).
967!> \endparblock
968! **************************************************************************************************
969 SUBROUTINE tddfpt_guess_vectors(evects, evals, gs_mos, log_unit, tddfpt_control, &
970 fm_pool_ao_mo_active, qs_env, nspins)
971 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout) :: evects
972 REAL(kind=dp), DIMENSION(:), INTENT(inout) :: evals
973 INTEGER, INTENT(in) :: nspins
974 TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
975 TYPE(cp_fm_pool_p_type), DIMENSION(:), INTENT(in) :: fm_pool_ao_mo_active
976 TYPE(tddfpt2_control_type), INTENT(in), POINTER :: tddfpt_control
977 INTEGER, INTENT(in) :: log_unit
978 TYPE(tddfpt_ground_state_mos), DIMENSION(nspins), &
979 INTENT(in) :: gs_mos
980
981 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_guess_vectors'
982
983 CHARACTER(len=5) :: spin_label1, spin_label2
984 INTEGER :: handle, i, imo_occ, imo_virt, ind, ispin, istate, j, jspin, k, no, nstates, &
985 nstates_occ_virt_alpha, nstates_selected, nv, spin1, spin2
986 INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
987 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: reverse_index
988 INTEGER, DIMENSION(maxspins) :: nmo, nmo_occ_avail, nmo_occ_selected, &
989 nmo_virt_selected
990 REAL(kind=dp) :: e_occ
991 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: e_virt_minus_occ, ev_occ, ev_virt
992 TYPE(excited_energy_type), POINTER :: ex_env
993
994 CALL timeset(routinen, handle)
995
996 nstates = SIZE(evects, 2)
997
998 IF (debug_this_module) THEN
999 cpassert(nstates > 0)
1000 cpassert(nspins == 1 .OR. nspins == 2)
1001 END IF
1002
1003 NULLIFY (ex_env)
1004 CALL get_qs_env(qs_env, exstate_env=ex_env)
1005
1006 DO ispin = 1, nspins
1007 ! number of occupied orbitals for each spin component
1008 nmo_occ_avail(ispin) = gs_mos(ispin)%nmo_active
1009 nmo(ispin) = gs_mos(ispin)%nmo_occ
1010 ! number of occupied and virtual orbitals which can potentially
1011 ! contribute to the excited states in question.
1012 nmo_occ_selected(ispin) = min(nmo_occ_avail(ispin), nstates)
1013 nmo_virt_selected(ispin) = min(SIZE(gs_mos(ispin)%evals_virt), nstates)
1014 END DO
1015
1016 ! TO DO: the variable 'nstates_selected' should probably be declared as INTEGER(kind=int_8),
1017 ! however we need a special version of the subroutine sort() in order to do so
1018 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1019 nstates_selected = dot_product(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins))
1020 ELSE
1021 nstates_selected = nmo_occ_selected(1)*nmo_virt_selected(2)
1022 END IF
1023
1024 ALLOCATE (inds(nstates_selected))
1025 ALLOCATE (e_virt_minus_occ(nstates_selected))
1026
1027 istate = 0
1028 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1029 ! Guess for spin-conserving TDDFT
1030 DO ispin = 1, nspins
1031 no = nmo_occ_selected(ispin)
1032 nv = nmo_virt_selected(ispin)
1033 ALLOCATE (ev_virt(nv), ev_occ(no))
1034 ! if do_bse and do_gw, take gw zeroth order
1035 IF ((tddfpt_control%do_bse) .OR. (tddfpt_control%do_bse_w_only) .OR. &
1036 (tddfpt_control%do_bse_gw_only)) THEN
1037 ev_virt(1:nv) = ex_env%gw_eigen(nmo(ispin) + 1:nmo(ispin) + nv)
1038 DO i = 1, no
1039 j = nmo_occ_avail(ispin) - i + 1
1040 k = gs_mos(ispin)%index_active(j)
1041 ev_occ(i) = ex_env%gw_eigen(k)
1042 END DO
1043 ELSE
1044 ev_virt(1:nv) = gs_mos(ispin)%evals_virt(1:nv)
1045 DO i = 1, no
1046 j = nmo_occ_avail(ispin) - i + 1
1047 k = gs_mos(ispin)%index_active(j)
1048 ev_occ(i) = gs_mos(ispin)%evals_occ(k)
1049 END DO
1050 END IF
1051
1052 DO imo_occ = 1, nmo_occ_selected(ispin)
1053 ! Here imo_occ enumerate Occupied orbitals in inverse order (from the last to the first element)
1054 e_occ = ev_occ(imo_occ)
1055 !
1056 DO imo_virt = 1, nmo_virt_selected(ispin)
1057 istate = istate + 1
1058 e_virt_minus_occ(istate) = ev_virt(imo_virt) - e_occ
1059 END DO
1060 END DO
1061
1062 DEALLOCATE (ev_virt, ev_occ)
1063 END DO
1064 ELSE
1065 ! Guess for spin-flip TDDFT
1066 DO imo_occ = 1, nmo_occ_selected(1)
1067 ! Here imo_occ enumerate alpha Occupied orbitals in inverse order (from the last to the first element)
1068 i = gs_mos(1)%nmo_active - imo_occ + 1
1069 k = gs_mos(1)%index_active(i)
1070 e_occ = gs_mos(1)%evals_occ(k)
1071
1072 DO imo_virt = 1, nmo_virt_selected(2)
1073 istate = istate + 1
1074 e_virt_minus_occ(istate) = gs_mos(2)%evals_virt(imo_virt) - e_occ
1075 END DO
1076 END DO
1077 END IF
1078
1079 IF (debug_this_module) THEN
1080 cpassert(istate == nstates_selected)
1081 END IF
1082
1083 CALL sort(e_virt_minus_occ, nstates_selected, inds)
1084
1085 ! Labels and spin component for closed-shell
1086 IF (nspins == 1) THEN
1087 spin1 = 1
1088 spin2 = spin1
1089 spin_label1 = ' '
1090 spin_label2 = spin_label1
1091 ! Labels and spin component for spin-flip excitations
1092 ELSE IF (tddfpt_control%spinflip /= no_sf_tddfpt) THEN
1093 spin1 = 1
1094 spin2 = 2
1095 spin_label1 = '(alp)'
1096 spin_label2 = '(bet)'
1097 END IF
1098
1099 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1100 ! Calculate maximum number of alpha excitations
1101 nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1)
1102 ELSE
1103 ! Calculate maximum number of spin-flip excitations
1104 nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(2)
1105 END IF
1106 IF (log_unit > 0) THEN
1107 WRITE (log_unit, "(1X,A)") "", &
1108 "-------------------------------------------------------------------------------", &
1109 "- TDDFPT Initial Guess -", &
1110 "-------------------------------------------------------------------------------"
1111 WRITE (log_unit, '(T11,A)') "State Occupied -> Virtual Excitation"
1112 WRITE (log_unit, '(T11,A)') "number orbital orbital energy (eV)"
1113 WRITE (log_unit, '(1X,79("-"))')
1114 END IF
1115
1116 i = maxval(nmo(:))
1117 ALLOCATE (reverse_index(i, nspins))
1118 reverse_index = 0
1119 DO ispin = 1, nspins
1120 DO i = 1, SIZE(gs_mos(ispin)%index_active)
1121 j = gs_mos(ispin)%index_active(i)
1122 reverse_index(j, ispin) = i
1123 END DO
1124 END DO
1125
1126 DO istate = 1, nstates
1127 IF (ASSOCIATED(evects(1, istate)%matrix_struct)) THEN
1128 ! Initial guess vector read from restart file
1129 IF (log_unit > 0) &
1130 WRITE (log_unit, '(T7,I8,T28,A19,T60,F14.5)') &
1131 istate, "*** restarted ***", evals(istate)*evolt
1132 ELSE
1133 ! New initial guess vector
1134 !
1135 ! Index of excited state - 1
1136 ind = inds(istate) - 1
1137
1138 ! Labels and spin component for open-shell spin-conserving excitations
1139 IF ((nspins > 1) .AND. (tddfpt_control%spinflip == no_sf_tddfpt)) THEN
1140 IF (ind < nstates_occ_virt_alpha) THEN
1141 spin1 = 1
1142 spin2 = 1
1143 spin_label1 = '(alp)'
1144 spin_label2 = '(alp)'
1145 ELSE
1146 ind = ind - nstates_occ_virt_alpha
1147 spin1 = 2
1148 spin2 = 2
1149 spin_label1 = '(bet)'
1150 spin_label2 = '(bet)'
1151 END IF
1152 END IF
1153
1154 ! Recover index of occupied MO (imo_occ) and unoccupied MO (imo_virt)
1155 ! associated to the excited state index (ind+1)
1156 i = ind/nmo_virt_selected(spin2) + 1
1157 j = nmo_occ_avail(spin1) - i + 1
1158 imo_occ = gs_mos(spin1)%index_active(j)
1159 imo_virt = mod(ind, nmo_virt_selected(spin2)) + 1
1160 ! Assign initial guess for excitation energy
1161 evals(istate) = e_virt_minus_occ(istate)
1162
1163 IF (log_unit > 0) &
1164 WRITE (log_unit, '(T7,I8,T24,I8,T37,A5,T45,I8,T54,A5,T60,F14.5)') &
1165 istate, imo_occ, spin_label1, nmo(spin2) + imo_virt, spin_label2, e_virt_minus_occ(istate)*evolt
1166
1167 DO jspin = 1, SIZE(evects, 1)
1168 ! .NOT. ASSOCIATED(evects(jspin, istate)%matrix_struct))
1169 CALL fm_pool_create_fm(fm_pool_ao_mo_active(jspin)%pool, evects(jspin, istate))
1170 CALL cp_fm_set_all(evects(jspin, istate), 0.0_dp)
1171
1172 IF (jspin == spin1) THEN
1173 ! Half transform excitation vector to ao space:
1174 ! evects_mi = c_ma*X_ai
1175 i = reverse_index(imo_occ, spin1)
1176 CALL cp_fm_to_fm(gs_mos(spin2)%mos_virt, evects(spin1, istate), &
1177 ncol=1, source_start=imo_virt, target_start=i)
1178 END IF
1179 END DO
1180 END IF
1181 END DO
1182
1183 DEALLOCATE (reverse_index)
1184
1185 IF (log_unit > 0) THEN
1186 WRITE (log_unit, '(/,T7,A,T50,I24)') 'Number of active states:', &
1187 tddfpt_total_number_of_states(tddfpt_control, gs_mos)
1188 WRITE (log_unit, "(1X,A)") &
1189 "-------------------------------------------------------------------------------"
1190 END IF
1191
1192 DEALLOCATE (e_virt_minus_occ)
1193 DEALLOCATE (inds)
1194
1195 CALL timestop(handle)
1196
1197 END SUBROUTINE tddfpt_guess_vectors
1198
1199END MODULE qs_tddfpt2_utils
Handles all functions related to the CELL.
Definition cell_types.F:15
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
pool for for elements that are retained and released
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Types for excited states potential energies.
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public oe_saop
integer, parameter, public cholesky_restore
integer, parameter, public cholesky_dbcsr
integer, parameter, public cholesky_off
integer, parameter, public oe_none
integer, parameter, public oe_shift
integer, parameter, public cholesky_inverse
integer, parameter, public oe_lb
integer, parameter, public no_sf_tddfpt
integer, parameter, public oe_gllb
objects that represent the structure of input sections and the data contained in an input section
recursive subroutine, public section_vals_create(section_vals, section)
creates a object where to store the values of a section
subroutine, public section_vals_retain(section_vals)
retains the given section values (see doc/ReferenceCounting.html)
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_set_subs_vals(section_vals, subsection_name, new_section_vals, i_rep_section)
replaces of the requested subsection with the one given
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
recursive subroutine, public section_vals_release(section_vals)
releases the given object
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public int_8
Definition kinds.F:54
integer, parameter, public dp
Definition kinds.F:34
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public evolt
Definition physcon.F:183
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, print_active, ext_ks_matrix, ext_xc_section)
routine where the real calculations are made: the KS matrix is calculated
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
subroutine, public eigensolver(matrix_ks_fm, mo_set, ortho, work, cholesky_method, do_level_shift, level_shift, matrix_u_fm, use_jacobi)
Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- vectors and MO eigenvalues....
Does all kind of post scf calculations for GPW/GAPW.
subroutine, public make_lumo_gpw(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, nlumo, nlumos)
Gets the lumos, and eigenvalues for the lumos.
module that contains the definitions of the scf types
integer, parameter, public ot_method_nr
subroutine, public tddfpt_release_ground_state_mos(gs_mos)
Release molecular orbitals.
subroutine, public tddfpt_init_ground_state_mos(gs_mos, mo_set, nlumo, blacs_env, cholesky_method, matrix_ks, matrix_s, mos_virt, evals_virt, qs_env)
Generate all virtual molecular orbitals for a given spin by diagonalising the corresponding Kohn-Sham...
subroutine, public tddfpt_oecorr(qs_env, gs_mos, matrix_ks_oep)
Callculate orbital corrected KS matrix for TDDFPT.
subroutine, public tddfpt_guess_vectors(evects, evals, gs_mos, log_unit, tddfpt_control, fm_pool_ao_mo_active, qs_env, nspins)
Generate missed guess vectors.
subroutine, public tddfpt_init_mos(qs_env, gs_mos, iounit)
Prepare MOs for TDDFPT Calculations.
pure integer(kind=int_8) function, public tddfpt_total_number_of_states(tddfpt_control, gs_mos)
Compute the number of possible singly excited states (occ -> virt)
All kind of helpful little routines.
Definition util.F:14
Calculate the saop potential.
Definition xc_pot_saop.F:11
subroutine, public add_saop_pot(ks_matrix, qs_env, oe_corr)
...
Calculation of KS matrix in xTB Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov JCTC 1...
subroutine, public build_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
...
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
represent a pointer to a 1d array
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Contains information on the excited states energy.
stores all the informations relevant to an mpi environment
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Ground state molecular orbitals.