(git:477b1f1)
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) THEN
137 NULLIFY (ks_env, ex_env)
138 CALL get_qs_env(qs_env, exstate_env=ex_env, ks_env=ks_env)
139 CALL dbcsr_copy(matrix_ks(1)%matrix, ex_env%matrix_ks(1)%matrix)
140 CALL set_ks_env(ks_env, matrix_ks=matrix_ks)
141 END IF
142
143 cpassert(.NOT. ASSOCIATED(gs_mos))
144 ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals
145 nspins = dft_control%nspins
146 ALLOCATE (gs_mos(nspins))
147
148 ! check if virtuals should be constructed for NAMD interface with NEWTONX
149 print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
150 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 NULLIFY (evals_virt, evals_virt_spin, mos_virt_spin)
156 IF (ASSOCIATED(scf_env)) THEN
157 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 nmo_virt = huge(0)
163 DO ispin = 1, nspins
164 CALL get_mo_set(mos(ispin), nmo=nmo_avail, homo=nmo_occ)
165 nmo_virt = min(nmo_virt, nmo_avail - nmo_occ)
166 END DO
167 ! number of unoccupied orbitals to compute
168 nmo_virt = tddfpt_control%nlumo - nmo_virt
169 IF (.NOT. print_virtuals_newtonx) THEN
170 IF (nmo_virt > 0) THEN
171 ALLOCATE (evals_virt(nspins), mos_virt(nspins))
172 ! the number of actually computed unoccupied orbitals will be stored as 'nmo_avail'
173 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 DO ispin = 1, nspins
180 IF (ASSOCIATED(evals_virt)) THEN
181 evals_virt_spin => evals_virt(ispin)%array
182 ELSE
183 NULLIFY (evals_virt_spin)
184 END IF
185 IF (ALLOCATED(mos_virt)) THEN
186 mos_virt_spin => mos_virt(ispin)
187 ELSE
188 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 qs_env=qs_env)
196 END DO
197
198 moc = 0
199 mvt = 0
200 DO ispin = 1, nspins
201 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, nrow_global=moc(1, ispin), ncol_global=moc(2, ispin))
202 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, nrow_global=mvt(1, ispin), ncol_global=mvt(2, ispin))
203 END DO
204 IF (iounit > 0) THEN
205 WRITE (iounit, "(T2,A,T36,A)") "TDDFPT| Molecular Orbitals:", &
206 " Spin AOs Occ Virt Total"
207 DO ispin = 1, nspins
208 WRITE (iounit, "(T2,A,T37,I4,4I10)") "TDDFPT| ", ispin, moc(1, ispin), moc(2, ispin), &
209 mvt(2, ispin), moc(2, ispin) + mvt(2, ispin)
210 END DO
211 END IF
212
213 IF (ASSOCIATED(evals_virt)) THEN
214 DO ispin = 1, SIZE(evals_virt)
215 IF (ASSOCIATED(evals_virt(ispin)%array)) DEALLOCATE (evals_virt(ispin)%array)
216 END DO
217 DEALLOCATE (evals_virt)
218 END IF
219
220 CALL cp_fm_release(mos_virt)
221
222 CALL timestop(handle)
223
224 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, &
265 sum_sign_array
266 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 POINTER :: my_block
271 REAL(kind=dp), DIMENSION(:), POINTER :: mo_evals_extended, mo_occ_extended, &
272 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 CALL timeset(routinen, handle)
284
285 NULLIFY (logger)
286 logger => cp_get_default_logger()
287 iounit = cp_logger_get_default_io_unit(logger)
288
289 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 nelectron=nelectrons, occupation_numbers=mo_occ_scf)
293
294 print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
295 CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases)
296
297 nmo_virt = nao - nmo_occ
298 IF (nlumo >= 0) &
299 nmo_virt = min(nmo_virt, nlumo)
300
301 IF (nmo_virt <= 0) &
302 CALL cp_abort(__location__, &
303 'At least one unoccupied molecular orbital is required to calculate excited states.')
304
305 do_eigen = .false.
306 ! diagonalise the Kohn-Sham matrix one more time if the number of available unoccupied states are too small
307 IF (ASSOCIATED(evals_virt)) THEN
308 cpassert(ASSOCIATED(mos_virt))
309 IF (nmo_virt > nmo_scf - nmo_occ + SIZE(evals_virt)) do_eigen = .true.
310 ELSE
311 IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .true.
312 END IF
313
314 ! ++ allocate storage space for gs_mos
315 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 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 ncol_global=nmo_virt, context=blacs_env)
321
322 NULLIFY (gs_mos%mos_occ, gs_mos%mos_virt, gs_mos%evals_occ_matrix)
323 ALLOCATE (gs_mos%mos_occ, gs_mos%mos_virt)
324 CALL cp_fm_create(gs_mos%mos_occ, ao_mo_occ_fm_struct)
325 CALL cp_fm_create(gs_mos%mos_virt, ao_mo_virt_fm_struct)
326 gs_mos%nmo_occ = nmo_occ
327
328 ALLOCATE (gs_mos%evals_occ(nmo_occ))
329 ALLOCATE (gs_mos%evals_virt(nmo_virt))
330 ALLOCATE (gs_mos%phases_occ(nmo_occ))
331 ALLOCATE (gs_mos%phases_virt(nmo_virt))
332
333 ! ++ nullify pointers
334 NULLIFY (ao_ao_fm_struct, wfn_fm_struct)
335 NULLIFY (mos_extended, mo_coeff_extended, mo_evals_extended, mo_occ_extended)
336
337 IF (do_eigen) THEN
338 ! ++ set of molecular orbitals
339 CALL cp_fm_struct_create(ao_ao_fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
340 CALL cp_fm_struct_create(wfn_fm_struct, nrow_global=nao, ncol_global=nmo_occ + nmo_virt, context=blacs_env)
341 ALLOCATE (mos_extended)
342 CALL allocate_mo_set(mos_extended, nao, nmo_occ + nmo_virt, nelectrons, &
343 REAL(nelectrons, dp), maxocc, flexible_electron_count=0.0_dp)
344 CALL init_mo_set(mos_extended, fm_struct=wfn_fm_struct, name="mos-extended")
345 CALL cp_fm_struct_release(wfn_fm_struct)
346 CALL get_mo_set(mos_extended, mo_coeff=mo_coeff_extended, &
347 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 DO imo = 1, nmo_scf
355 mo_occ_extended(imo) = mo_occ_scf(imo)
356 END DO
357 mo_occ_extended(nmo_scf + 1:) = 0.0_dp
358
359 ! ++ allocate temporary matrices
360 CALL cp_fm_create(matrix_ks_fm, ao_ao_fm_struct)
361 CALL cp_fm_create(ortho_fm, ao_ao_fm_struct)
362 CALL cp_fm_create(work_fm, ao_ao_fm_struct)
363 CALL cp_fm_struct_release(ao_ao_fm_struct)
364
365 ! some stuff from the subroutine general_eigenproblem()
366 CALL copy_dbcsr_to_fm(matrix_s, ortho_fm)
367 CALL copy_dbcsr_to_fm(matrix_ks, matrix_ks_fm)
368
369 IF (cholesky_method == cholesky_dbcsr) THEN
370 cpabort('CHOLESKY DBCSR_INVERSE is not implemented in TDDFT.')
371 ELSE IF (cholesky_method == cholesky_off) THEN
372 cpabort('CHOLESKY OFF is not implemented in TDDFT.')
373 ELSE
374 CALL cp_fm_cholesky_decompose(ortho_fm)
375 IF (cholesky_method == cholesky_inverse) THEN
376 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 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 do_level_shift=.false., level_shift=0.0_dp, use_jacobi=.false.)
385 END IF
386
387 ! -- clean up needless matrices
388 CALL cp_fm_release(work_fm)
389 CALL cp_fm_release(ortho_fm)
390 CALL cp_fm_release(matrix_ks_fm)
391 ELSE
392 CALL get_mo_set(mo_set, mo_coeff=mo_coeff_extended, &
393 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 CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
400 CALL cp_fm_struct_release(ao_mo_occ_fm_struct)
401
402 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 row_indices=row_indices, col_indices=col_indices, local_data=my_block)
405
406 ALLOCATE (minrow_neg_array(nmo_occ), minrow_pos_array(nmo_occ), sum_sign_array(nmo_occ))
407 minrow_neg_array(:) = nao
408 minrow_pos_array(:) = nao
409 sum_sign_array(:) = 0
410 DO icol_local = 1, ncol_local
411 icol_global = col_indices(icol_local)
412
413 DO irow_local = 1, nrow_local
414 element = my_block(irow_local, icol_local)
415
416 sign_int = 0
417 IF (element >= eps_dp) THEN
418 sign_int = 1
419 ELSE IF (element <= -eps_dp) THEN
420 sign_int = -1
421 END IF
422
423 sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
424
425 irow_global = row_indices(irow_local)
426 IF (sign_int > 0) THEN
427 IF (minrow_pos_array(icol_global) > irow_global) &
428 minrow_pos_array(icol_global) = irow_global
429 ELSE IF (sign_int < 0) THEN
430 IF (minrow_neg_array(icol_global) > irow_global) &
431 minrow_neg_array(icol_global) = irow_global
432 END IF
433 END DO
434 END DO
435
436 CALL para_env%sum(sum_sign_array)
437 CALL para_env%min(minrow_neg_array)
438 CALL para_env%min(minrow_pos_array)
439
440 DO icol_local = 1, nmo_occ
441 IF (sum_sign_array(icol_local) > 0) THEN
442 ! most of the expansion coefficients are positive => MO's phase = +1
443 gs_mos%phases_occ(icol_local) = 1.0_dp
444 ELSE IF (sum_sign_array(icol_local) < 0) THEN
445 ! most of the expansion coefficients are negative => MO's phase = -1
446 gs_mos%phases_occ(icol_local) = -1.0_dp
447 ELSE
448 ! equal number of positive and negative expansion coefficients
449 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 gs_mos%phases_occ(icol_local) = 1.0_dp
453 ELSE
454 ! MO's phase = -1
455 gs_mos%phases_occ(icol_local) = -1.0_dp
456 END IF
457 END IF
458 END DO
459
460 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 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_occ, ncol=nmo_occ, source_start=1, target_start=1)
464 gs_mos%evals_occ(1:nmo_occ) = mo_evals_extended(1:nmo_occ)
465
466 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 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 source_start=1, target_start=nmo_scf - nmo_occ + 1)
471 gs_mos%evals_virt(1:nmo_scf - nmo_occ) = evals_virt(nmo_occ + 1:nmo_occ + nmo_scf)
472 gs_mos%evals_virt(nmo_scf - nmo_occ + 1:nmo_virt) = evals_virt(1:nmo_virt - (nmo_scf - nmo_occ))
473 ELSE
474 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_virt, source_start=nmo_occ + 1, target_start=1)
475 gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt)
476 END IF
477
478 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 CALL cp_fm_create(work_fm_virt, ao_mo_virt_fm_struct)
483
484 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 row_indices=row_indices, col_indices=col_indices, local_data=my_block)
487
488 ALLOCATE (minrow_neg_array(nmo_virt), minrow_pos_array(nmo_virt), sum_sign_array(nmo_virt))
489 minrow_neg_array(:) = nao
490 minrow_pos_array(:) = nao
491 sum_sign_array(:) = 0
492 DO icol_local = 1, ncol_local
493 icol_global = col_indices(icol_local)
494
495 DO irow_local = 1, nrow_local
496 element = my_block(irow_local, icol_local)
497
498 sign_int = 0
499 IF (element >= eps_dp) THEN
500 sign_int = 1
501 ELSE IF (element <= -eps_dp) THEN
502 sign_int = -1
503 END IF
504
505 sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
506
507 irow_global = row_indices(irow_local)
508 IF (sign_int > 0) THEN
509 IF (minrow_pos_array(icol_global) > irow_global) &
510 minrow_pos_array(icol_global) = irow_global
511 ELSE IF (sign_int < 0) THEN
512 IF (minrow_neg_array(icol_global) > irow_global) &
513 minrow_neg_array(icol_global) = irow_global
514 END IF
515 END DO
516 END DO
517
518 CALL para_env%sum(sum_sign_array)
519 CALL para_env%min(minrow_neg_array)
520 CALL para_env%min(minrow_pos_array)
521 DO icol_local = 1, nmo_virt
522 IF (sum_sign_array(icol_local) > 0) THEN
523 ! most of the expansion coefficients are positive => MO's phase = +1
524 gs_mos%phases_virt(icol_local) = 1.0_dp
525 ELSE IF (sum_sign_array(icol_local) < 0) THEN
526 ! most of the expansion coefficients are negative => MO's phase = -1
527 gs_mos%phases_virt(icol_local) = -1.0_dp
528 ELSE
529 ! equal number of positive and negative expansion coefficients
530 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 gs_mos%phases_virt(icol_local) = 1.0_dp
534 ELSE
535 ! MO's phase = -1
536 gs_mos%phases_virt(icol_local) = -1.0_dp
537 END IF
538 END IF
539 END DO
540 DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
541 CALL cp_fm_release(work_fm_virt)
542 END IF !print_phases
543 CALL cp_fm_struct_release(ao_mo_virt_fm_struct) ! here after print_phases
544
545 CALL cp_fm_release(work_fm)
546
547 IF (do_eigen) THEN
548 CALL deallocate_mo_set(mos_extended)
549 DEALLOCATE (mos_extended)
550 END IF
551
552 CALL timestop(handle)
553
554 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! **************************************************************************************************
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 CALL timeset(routinen, handle)
570
571 IF (ALLOCATED(gs_mos%phases_occ)) &
572 DEALLOCATE (gs_mos%phases_occ)
573
574 IF (ALLOCATED(gs_mos%evals_virt)) &
575 DEALLOCATE (gs_mos%evals_virt)
576
577 IF (ALLOCATED(gs_mos%evals_occ)) &
578 DEALLOCATE (gs_mos%evals_occ)
579
580 IF (ALLOCATED(gs_mos%phases_virt)) &
581 DEALLOCATE (gs_mos%phases_virt)
582
583 IF (ALLOCATED(gs_mos%index_active)) &
584 DEALLOCATE (gs_mos%index_active)
585
586 IF (ASSOCIATED(gs_mos%evals_occ_matrix)) THEN
587 CALL cp_fm_release(gs_mos%evals_occ_matrix)
588 DEALLOCATE (gs_mos%evals_occ_matrix)
589 END IF
590
591 IF (ASSOCIATED(gs_mos%mos_virt)) THEN
592 CALL cp_fm_release(gs_mos%mos_virt)
593 DEALLOCATE (gs_mos%mos_virt)
594 END IF
595
596 IF (ASSOCIATED(gs_mos%mos_occ)) THEN
597 CALL cp_fm_release(gs_mos%mos_occ)
598 DEALLOCATE (gs_mos%mos_occ)
599 END IF
600
601 IF (ASSOCIATED(gs_mos%mos_active)) THEN
602 CALL cp_fm_release(gs_mos%mos_active)
603 DEALLOCATE (gs_mos%mos_active)
604 END IF
605
606 CALL timestop(handle)
607
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 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 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 CALL timeset(routinen, handle)
639
640 NULLIFY (logger)
641 logger => cp_get_default_logger()
642 iounit = cp_logger_get_default_io_unit(logger)
643
644 CALL get_qs_env(qs_env, blacs_env=blacs_env, dft_control=dft_control, matrix_ks=matrix_ks)
645 tddfpt_control => dft_control%tddfpt2_control
646
647 ! obtain corrected KS-matrix
648 ! We should 'save' the energy values?
649 nspins = SIZE(gs_mos)
650 NULLIFY (matrix_ks_oep)
651 IF (tddfpt_control%oe_corr /= oe_none) THEN
652 IF (iounit > 0) THEN
653 WRITE (iounit, "(1X,A)") "", &
654 "-------------------------------------------------------------------------------", &
655 "- Orbital Eigenvalue Correction Started -", &
656 "-------------------------------------------------------------------------------"
657 END IF
658
659 CALL cp_warn(__location__, &
660 "Orbital energy correction potential is an experimental feature. "// &
661 "Use it with extreme care")
662
663 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
664 CALL section_vals_get(hfx_section, explicit=do_hfx)
665 IF (do_hfx) THEN
666 CALL cp_abort(__location__, &
667 "Implementation of orbital energy correction XC-potentials is "// &
668 "currently incompatible with exact-exchange functionals")
669 END IF
670
671 CALL dbcsr_allocate_matrix_set(matrix_ks_oep, nspins)
672 DO ispin = 1, nspins
673 CALL dbcsr_init_p(matrix_ks_oep(ispin)%matrix)
674 CALL dbcsr_copy(matrix_ks_oep(ispin)%matrix, matrix_ks(ispin)%matrix)
675 END DO
676
677 ! KS-matrix without XC-terms
678 xc_fun_original => section_vals_get_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL")
679 CALL section_vals_retain(xc_fun_original)
680 NULLIFY (xc_fun_empty)
681 CALL section_vals_create(xc_fun_empty, xc_fun_original%section)
682 CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_empty)
683 CALL section_vals_release(xc_fun_empty)
684
685 IF (dft_control%qs_control%semi_empirical) THEN
686 cpabort("TDDFPT with SE not possible")
687 ELSEIF (dft_control%qs_control%dftb) THEN
688 cpabort("TDDFPT with DFTB not possible")
689 ELSEIF (dft_control%qs_control%xtb) THEN
690 CALL build_xtb_ks_matrix(qs_env, calculate_forces=.false., just_energy=.false., &
691 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 ext_ks_matrix=matrix_ks_oep)
695 END IF
696
697 IF (tddfpt_control%oe_corr == oe_saop .OR. &
698 tddfpt_control%oe_corr == oe_lb .OR. &
699 tddfpt_control%oe_corr == oe_gllb) THEN
700 IF (iounit > 0) THEN
701 WRITE (iounit, "(T2,A)") " Orbital energy correction of SAOP type "
702 END IF
703 CALL add_saop_pot(matrix_ks_oep, qs_env, tddfpt_control%oe_corr)
704 ELSE IF (tddfpt_control%oe_corr == oe_shift) THEN
705 IF (iounit > 0) THEN
706 WRITE (iounit, "(T2,A,T71,F10.3)") &
707 " Virtual Orbital Eigenvalue Shift [eV] ", tddfpt_control%ev_shift*evolt
708 WRITE (iounit, "(T2,A,T71,F10.3)") &
709 " 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 tddfpt_control%ev_shift, tddfpt_control%eos_shift)
713 ELSE
714 CALL cp_abort(__location__, &
715 "Unimplemented orbital energy correction potential")
716 END IF
717 CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_original)
718 CALL section_vals_release(xc_fun_original)
719
720 ! compute 'evals_occ_matrix'
721 CALL dbcsr_get_info(matrix_ks(1)%matrix, nfullrows_total=nao)
722 NULLIFY (mo_occ_mo_occ_fm_struct)
723 DO ispin = 1, nspins
724 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 context=blacs_env)
727 ALLOCATE (gs_mos(ispin)%evals_occ_matrix)
728 CALL cp_fm_create(gs_mos(ispin)%evals_occ_matrix, mo_occ_mo_occ_fm_struct)
729 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 context=blacs_env)
733 CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
734 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 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 0.0_dp, gs_mos(ispin)%evals_occ_matrix)
739 CALL cp_fm_release(work_fm)
740 END DO
741 IF (iounit > 0) THEN
742 WRITE (iounit, "(1X,A)") &
743 "-------------------------------------------------------------------------------"
744 END IF
745
746 END IF
747
748 CALL timestop(handle)
749
750 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 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 nstates_total = 0
769 nspins = SIZE(gs_mos)
770
771 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
772 ! Total number of possible excitations for spin-conserving TDDFT
773 DO ispin = 1, nspins
774 nstates_total = nstates_total + &
775 gs_mos(ispin)%nmo_active* &
776 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 SIZE(gs_mos(2)%evals_virt, kind=int_8)
782 END IF
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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
814 TYPE(dbcsr_type), POINTER :: smat
815 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
816
817 CALL timeset(routinen, handle)
818
819 n_spins = SIZE(gs_mos)
820 cpassert(n_spins == SIZE(matrix_ks))
821
822 IF (eos_shift /= 0.0_dp .AND. n_spins > 1) THEN
823 cpabort("eos_shift not implemented")
824 CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
825 smat => matrix_s(1)%matrix
826 CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
827 CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
828 nl = min(na, nb)
829 nu = max(na, nb)
830 ! open shell orbital shift
831 DO ispin = 1, n_spins
832 coeff => gs_mos(ispin)%mos_occ
833 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
834 IF (nhomo == nu) THEN
835 ! downshift with -eos_shift using occupied orbitals
836 nos = nu - nl
837 CALL cp_fm_create(cmos, fmstruct)
838 CALL cp_fm_get_info(coeff, nrow_global=nrow)
839 CALL cp_fm_to_fm_submat(coeff, cmos, nrow, nos, 1, nl + 1, 1, 1)
840 CALL cp_fm_create(cvec, fmstruct)
841 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 alpha=-eos_shift, keep_sparsity=.true.)
844 CALL cp_fm_release(cmos)
845 CALL cp_fm_release(cvec)
846 ELSE
847 ! upshift with eos_shift using virtual orbitals
848 coeff => gs_mos(ispin)%mos_virt
849 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
850 nos = nu - nhomo
851 cpassert(nvirt >= nos)
852 CALL cp_fm_create(cvec, fmstruct)
853 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 alpha=eos_shift, keep_sparsity=.true.)
856 CALL cp_fm_release(cvec)
857 END IF
858 END DO
859 ! virtual shift
860 IF (ev_shift /= 0.0_dp) THEN
861 DO ispin = 1, n_spins
862 CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
863 alpha_scalar=1.0_dp, beta_scalar=ev_shift)
864 coeff => gs_mos(ispin)%mos_occ
865 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
866 CALL cp_fm_create(cvec, fmstruct)
867 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 alpha=-ev_shift, keep_sparsity=.true.)
870 CALL cp_fm_release(cvec)
871 IF (nhomo < nu) THEN
872 nos = nu - nhomo
873 coeff => gs_mos(ispin)%mos_virt
874 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
875 cpassert(nvirt >= nos)
876 CALL cp_fm_create(cvec, fmstruct)
877 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 alpha=-ev_shift, keep_sparsity=.true.)
880 CALL cp_fm_release(cvec)
881 END IF
882 END DO
883 END IF
884 ELSE
885 ! virtual shift
886 IF (ev_shift /= 0.0_dp) THEN
887 CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
888 smat => matrix_s(1)%matrix
889 DO ispin = 1, n_spins
890 CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
891 alpha_scalar=1.0_dp, beta_scalar=ev_shift)
892 coeff => gs_mos(ispin)%mos_occ
893 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
894 CALL cp_fm_create(cvec, fmstruct)
895 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 alpha=-ev_shift, keep_sparsity=.true.)
898 CALL cp_fm_release(cvec)
899 END DO
900 END IF
901 END IF
902 ! set eigenvalues
903 IF (eos_shift == 0.0_dp .OR. n_spins == 1) THEN
904 DO ispin = 1, n_spins
905 IF (ALLOCATED(gs_mos(ispin)%evals_virt)) THEN
906 gs_mos(ispin)%evals_virt(:) = gs_mos(ispin)%evals_virt(:) + ev_shift
907 END IF
908 END DO
909 ELSE
910 CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
911 CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
912 nl = min(na, nb)
913 nu = max(na, nb)
914 nos = nu - nl
915 IF (na == nu) THEN
916 IF (ALLOCATED(gs_mos(1)%evals_occ)) THEN
917 gs_mos(1)%evals_occ(nl + 1:nu) = gs_mos(1)%evals_occ(nl + 1:nu) - eos_shift
918 END IF
919 IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
920 gs_mos(1)%evals_virt(:) = gs_mos(1)%evals_virt(:) + ev_shift
921 END IF
922 IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
923 gs_mos(2)%evals_virt(1:nos) = gs_mos(2)%evals_virt(1:nos) + eos_shift
924 gs_mos(2)%evals_virt(nos + 1:) = gs_mos(2)%evals_virt(nos + 1:) + ev_shift
925 END IF
926 ELSE
927 IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
928 gs_mos(1)%evals_virt(1:nos) = gs_mos(1)%evals_virt(1:nos) + eos_shift
929 gs_mos(1)%evals_virt(nos + 1:) = gs_mos(1)%evals_virt(nos + 1:) + ev_shift
930 END IF
931 IF (ALLOCATED(gs_mos(2)%evals_occ)) THEN
932 gs_mos(2)%evals_occ(nl + 1:nu) = gs_mos(2)%evals_occ(nl + 1:nu) - eos_shift
933 END IF
934 IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
935 gs_mos(2)%evals_virt(:) = gs_mos(2)%evals_virt(:) + ev_shift
936 END IF
937 END IF
938 END IF
939
940 CALL timestop(handle)
941
942 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 SUBROUTINE tddfpt_guess_vectors(evects, evals, gs_mos, log_unit, tddfpt_control, &
969 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
986 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: e_virt_minus_occ, ev_occ, ev_virt
991 TYPE(excited_energy_type), POINTER :: ex_env
992
993 CALL timeset(routinen, handle)
994
995 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 NULLIFY (ex_env)
1003 CALL get_qs_env(qs_env, exstate_env=ex_env)
1004
1005 DO ispin = 1, nspins
1006 ! number of occupied orbitals for each spin component
1007 nmo_occ_avail(ispin) = gs_mos(ispin)%nmo_active
1008 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 nmo_occ_selected(ispin) = min(nmo_occ_avail(ispin), nstates)
1012 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 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1018 nstates_selected = dot_product(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins))
1019 ELSE
1020 nstates_selected = nmo_occ_selected(1)*nmo_virt_selected(2)
1021 END IF
1022
1023 ALLOCATE (inds(nstates_selected))
1024 ALLOCATE (e_virt_minus_occ(nstates_selected))
1025
1026 istate = 0
1027 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1028 ! Guess for spin-conserving TDDFT
1029 DO ispin = 1, nspins
1030 no = nmo_occ_selected(ispin)
1031 nv = nmo_virt_selected(ispin)
1032 ALLOCATE (ev_virt(nv), ev_occ(no))
1033 ! if do_bse and do_gw, take gw zeroth order
1034 IF (tddfpt_control%do_bse) THEN
1035 ev_virt(1:nv) = ex_env%gw_eigen(nmo(ispin) + 1:nmo(ispin) + nv)
1036 DO i = 1, no
1037 j = nmo_occ_avail(ispin) - i + 1
1038 k = gs_mos(ispin)%index_active(j)
1039 ev_occ(i) = ex_env%gw_eigen(k)
1040 END DO
1041 ELSE
1042 ev_virt(1:nv) = gs_mos(ispin)%evals_virt(1:nv)
1043 DO i = 1, no
1044 j = nmo_occ_avail(ispin) - i + 1
1045 k = gs_mos(ispin)%index_active(j)
1046 ev_occ(i) = gs_mos(ispin)%evals_occ(k)
1047 END DO
1048 END IF
1049
1050 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 e_occ = ev_occ(imo_occ)
1053 !
1054 DO imo_virt = 1, nmo_virt_selected(ispin)
1055 istate = istate + 1
1056 e_virt_minus_occ(istate) = ev_virt(imo_virt) - e_occ
1057 END DO
1058 END DO
1059
1060 DEALLOCATE (ev_virt, ev_occ)
1061 END DO
1062 ELSE
1063 ! Guess for spin-flip TDDFT
1064 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 i = gs_mos(1)%nmo_active - imo_occ + 1
1067 k = gs_mos(1)%index_active(i)
1068 e_occ = gs_mos(1)%evals_occ(k)
1069
1070 DO imo_virt = 1, nmo_virt_selected(2)
1071 istate = istate + 1
1072 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 CALL sort(e_virt_minus_occ, nstates_selected, inds)
1082
1083 ! Labels and spin component for closed-shell
1084 IF (nspins == 1) THEN
1085 spin1 = 1
1086 spin2 = spin1
1087 spin_label1 = ' '
1088 spin_label2 = spin_label1
1089 ! Labels and spin component for spin-flip excitations
1090 ELSE IF (tddfpt_control%spinflip /= no_sf_tddfpt) THEN
1091 spin1 = 1
1092 spin2 = 2
1093 spin_label1 = '(alp)'
1094 spin_label2 = '(bet)'
1095 END IF
1096
1097 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
1098 ! Calculate maximum number of alpha excitations
1099 nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1)
1100 ELSE
1101 ! Calculate maximum number of spin-flip excitations
1102 nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(2)
1103 END IF
1104 IF (log_unit > 0) THEN
1105 WRITE (log_unit, "(1X,A)") "", &
1106 "-------------------------------------------------------------------------------", &
1107 "- TDDFPT Initial Guess -", &
1108 "-------------------------------------------------------------------------------"
1109 WRITE (log_unit, '(T11,A)') "State Occupied -> Virtual Excitation"
1110 WRITE (log_unit, '(T11,A)') "number orbital orbital energy (eV)"
1111 WRITE (log_unit, '(1X,79("-"))')
1112 END IF
1113
1114 i = maxval(nmo(:))
1115 ALLOCATE (reverse_index(i, nspins))
1116 reverse_index = 0
1117 DO ispin = 1, nspins
1118 DO i = 1, SIZE(gs_mos(ispin)%index_active)
1119 j = gs_mos(ispin)%index_active(i)
1120 reverse_index(j, ispin) = i
1121 END DO
1122 END DO
1123
1124 DO istate = 1, nstates
1125 IF (ASSOCIATED(evects(1, istate)%matrix_struct)) THEN
1126 ! Initial guess vector read from restart file
1127 IF (log_unit > 0) &
1128 WRITE (log_unit, '(T7,I8,T28,A19,T60,F14.5)') &
1129 istate, "*** restarted ***", evals(istate)*evolt
1130 ELSE
1131 ! New initial guess vector
1132 !
1133 ! Index of excited state - 1
1134 ind = inds(istate) - 1
1135
1136 ! Labels and spin component for open-shell spin-conserving excitations
1137 IF ((nspins > 1) .AND. (tddfpt_control%spinflip == no_sf_tddfpt)) THEN
1138 IF (ind < nstates_occ_virt_alpha) THEN
1139 spin1 = 1
1140 spin2 = 1
1141 spin_label1 = '(alp)'
1142 spin_label2 = '(alp)'
1143 ELSE
1144 ind = ind - nstates_occ_virt_alpha
1145 spin1 = 2
1146 spin2 = 2
1147 spin_label1 = '(bet)'
1148 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 i = ind/nmo_virt_selected(spin2) + 1
1155 j = nmo_occ_avail(spin1) - i + 1
1156 imo_occ = gs_mos(spin1)%index_active(j)
1157 imo_virt = mod(ind, nmo_virt_selected(spin2)) + 1
1158 ! Assign initial guess for excitation energy
1159 evals(istate) = e_virt_minus_occ(istate)
1160
1161 IF (log_unit > 0) &
1162 WRITE (log_unit, '(T7,I8,T24,I8,T37,A5,T45,I8,T54,A5,T60,F14.5)') &
1163 istate, imo_occ, spin_label1, nmo(spin2) + imo_virt, spin_label2, e_virt_minus_occ(istate)*evolt
1164
1165 DO jspin = 1, SIZE(evects, 1)
1166 ! .NOT. ASSOCIATED(evects(jspin, istate)%matrix_struct))
1167 CALL fm_pool_create_fm(fm_pool_ao_mo_active(jspin)%pool, evects(jspin, istate))
1168 CALL cp_fm_set_all(evects(jspin, istate), 0.0_dp)
1169
1170 IF (jspin == spin1) THEN
1171 ! Half transform excitation vector to ao space:
1172 ! evects_mi = c_ma*X_ai
1173 i = reverse_index(imo_occ, spin1)
1174 CALL cp_fm_to_fm(gs_mos(spin2)%mos_virt, evects(spin1, istate), &
1175 ncol=1, source_start=imo_virt, target_start=i)
1176 END IF
1177 END DO
1178 END IF
1179 END DO
1180
1181 DEALLOCATE (reverse_index)
1182
1183 IF (log_unit > 0) THEN
1184 WRITE (log_unit, '(/,T7,A,T50,I24)') 'Number of active states:', &
1185 tddfpt_total_number_of_states(tddfpt_control, gs_mos)
1186 WRITE (log_unit, "(1X,A)") &
1187 "-------------------------------------------------------------------------------"
1188 END IF
1189
1190 DEALLOCATE (e_virt_minus_occ)
1191 DEALLOCATE (inds)
1192
1193 CALL timestop(handle)
1194
1195 END SUBROUTINE tddfpt_guess_vectors
1196
1197END 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, 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, 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, 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.