(git:34ef472)
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-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 MODULE qs_tddfpt2_utils
9  USE cell_types, ONLY: cell_type
10  USE cp_array_utils, ONLY: cp_1d_r_p_type
11  USE cp_blacs_env, ONLY: cp_blacs_env_type
12  USE cp_control_types, ONLY: dft_control_type,&
13  tddfpt2_control_type
14  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
15  cp_dbcsr_plus_fm_fm_t,&
16  cp_dbcsr_sm_fm_multiply,&
17  dbcsr_allocate_matrix_set
18  USE cp_fm_basic_linalg, ONLY: cp_fm_triangular_invert
19  USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
20  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
21  cp_fm_struct_p_type,&
22  cp_fm_struct_release,&
23  cp_fm_struct_type
24  USE cp_fm_types, ONLY: cp_fm_create,&
25  cp_fm_get_info,&
26  cp_fm_release,&
27  cp_fm_set_all,&
28  cp_fm_to_fm,&
29  cp_fm_to_fm_submat,&
30  cp_fm_type
31  USE cp_log_handling, ONLY: cp_get_default_logger,&
32  cp_logger_get_default_io_unit,&
33  cp_logger_type
34  USE dbcsr_api, ONLY: dbcsr_add,&
35  dbcsr_copy,&
36  dbcsr_get_info,&
37  dbcsr_init_p,&
38  dbcsr_p_type,&
39  dbcsr_type
40  USE input_constants, ONLY: &
41  cholesky_dbcsr, cholesky_inverse, cholesky_off, cholesky_restore, oe_gllb, oe_lb, oe_none, &
42  oe_saop, oe_shift
43  USE input_section_types, ONLY: section_vals_create,&
44  section_vals_get,&
45  section_vals_get_subs_vals,&
46  section_vals_release,&
47  section_vals_retain,&
48  section_vals_set_subs_vals,&
49  section_vals_type,&
50  section_vals_val_get
51  USE kinds, ONLY: dp,&
52  int_8
53  USE message_passing, ONLY: mp_para_env_type
54  USE parallel_gemm_api, ONLY: parallel_gemm
55  USE physcon, ONLY: evolt
56  USE qs_environment_types, ONLY: get_qs_env,&
57  qs_environment_type
58  USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
59  USE qs_mo_types, ONLY: allocate_mo_set,&
60  deallocate_mo_set,&
61  get_mo_set,&
62  init_mo_set,&
63  mo_set_type
64  USE qs_scf_methods, ONLY: eigensolver
65  USE qs_scf_post_gpw, ONLY: make_lumo_gpw
66  USE qs_scf_types, ONLY: ot_method_nr,&
67  qs_scf_env_type
68  USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos
69  USE util, ONLY: sort
70  USE xc_pot_saop, ONLY: add_saop_pot
71  USE xtb_matrices, ONLY: build_xtb_ks_matrix
72 #include "./base/base_uses.f90"
73 
74  IMPLICIT NONE
75 
76  PRIVATE
77 
78  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_utils'
79 
80  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
81  ! number of first derivative components (3: d/dx, d/dy, d/dz)
82  INTEGER, PARAMETER, PRIVATE :: nderivs = 3
83  INTEGER, PARAMETER, PRIVATE :: maxspins = 2
84 
85  PUBLIC :: tddfpt_init_ground_state_mos, tddfpt_release_ground_state_mos
86  PUBLIC :: tddfpt_guess_vectors, tddfpt_init_mos, tddfpt_oecorr
87  PUBLIC :: tddfpt_total_number_of_states
88 
89 ! **************************************************************************************************
90 
91 CONTAINS
92 
93 ! **************************************************************************************************
94 !> \brief Prepare MOs for TDDFPT Calculations
95 !> \param qs_env Quickstep environment
96 !> \param gs_mos ...
97 !> \param iounit ...
98 ! **************************************************************************************************
99  SUBROUTINE tddfpt_init_mos(qs_env, gs_mos, iounit)
100  TYPE(qs_environment_type), POINTER :: qs_env
101  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
102  POINTER :: gs_mos
103  INTEGER, INTENT(IN) :: iounit
104 
105  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_init_mos'
106 
107  INTEGER :: handle, ispin, nmo_avail, nmo_occ, &
108  nmo_virt, nspins
109  INTEGER, DIMENSION(2, 2) :: moc, mvt
110  LOGICAL :: print_virtuals_newtonx
111  REAL(kind=dp), DIMENSION(:), POINTER :: evals_virt_spin
112  TYPE(cell_type), POINTER :: cell
113  TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: evals_virt
114  TYPE(cp_blacs_env_type), POINTER :: blacs_env
115  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
116  TARGET :: mos_virt
117  TYPE(cp_fm_type), POINTER :: mos_virt_spin
118  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
119  TYPE(dft_control_type), POINTER :: dft_control
120  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
121  TYPE(qs_scf_env_type), POINTER :: scf_env
122  TYPE(section_vals_type), POINTER :: print_section
123  TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
124 
125  CALL timeset(routinen, handle)
126 
127  CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, &
128  matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, scf_env=scf_env)
129  tddfpt_control => dft_control%tddfpt2_control
130 
131  cpassert(.NOT. ASSOCIATED(gs_mos))
132  ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals
133  nspins = dft_control%nspins
134  ALLOCATE (gs_mos(nspins))
135 
136  ! check if virtuals should be constructed for NAMD interface with NEWTONX
137  print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
138  CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals_newtonx)
139 
140  ! when the number of unoccupied orbitals is limited and OT has been used
141  ! for the ground-state DFT calculation,
142  ! compute the missing unoccupied orbitals using OT as well.
143  NULLIFY (evals_virt, evals_virt_spin, mos_virt_spin)
144  IF (ASSOCIATED(scf_env)) THEN
145  IF ((scf_env%method == ot_method_nr .AND. tddfpt_control%nlumo > 0) .OR. &
146  (scf_env%method == ot_method_nr .AND. print_virtuals_newtonx)) THEN
147  ! As OT with ADDED_MOS/=0 is currently not implemented, the following block is equivalent to:
148  ! nmo_virt = tddfpt_control%nlumo
149  ! number of already computed unoccupied orbitals (added_mos) .
150  nmo_virt = huge(0)
151  DO ispin = 1, nspins
152  CALL get_mo_set(mos(ispin), nmo=nmo_avail, homo=nmo_occ)
153  nmo_virt = min(nmo_virt, nmo_avail - nmo_occ)
154  END DO
155  ! number of unoccupied orbitals to compute
156  nmo_virt = tddfpt_control%nlumo - nmo_virt
157  IF (.NOT. print_virtuals_newtonx) THEN
158  IF (nmo_virt > 0) THEN
159  ALLOCATE (evals_virt(nspins), mos_virt(nspins))
160  ! the number of actually computed unoccupied orbitals will be stored as 'nmo_avail'
161  CALL make_lumo_gpw(qs_env, scf_env, mos_virt, evals_virt, nmo_virt, nmo_avail)
162  END IF
163  END IF
164  END IF
165  END IF
166 
167  DO ispin = 1, nspins
168  IF (ASSOCIATED(evals_virt)) THEN
169  evals_virt_spin => evals_virt(ispin)%array
170  ELSE
171  NULLIFY (evals_virt_spin)
172  END IF
173  IF (ALLOCATED(mos_virt)) THEN
174  mos_virt_spin => mos_virt(ispin)
175  ELSE
176  NULLIFY (mos_virt_spin)
177  END IF
178  CALL tddfpt_init_ground_state_mos(gs_mos=gs_mos(ispin), mo_set=mos(ispin), &
179  nlumo=tddfpt_control%nlumo, &
180  blacs_env=blacs_env, cholesky_method=cholesky_restore, &
181  matrix_ks=matrix_ks(ispin)%matrix, matrix_s=matrix_s(1)%matrix, &
182  mos_virt=mos_virt_spin, evals_virt=evals_virt_spin, &
183  qs_env=qs_env)
184  END DO
185 
186  moc = 0
187  mvt = 0
188  DO ispin = 1, nspins
189  CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, nrow_global=moc(1, ispin), ncol_global=moc(2, ispin))
190  CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, nrow_global=mvt(1, ispin), ncol_global=mvt(2, ispin))
191  END DO
192  IF (iounit > 0) THEN
193  WRITE (iounit, "(T2,A,T36,A)") "TDDFPT| Molecular Orbitals:", &
194  " Spin AOs Occ Virt Total"
195  DO ispin = 1, nspins
196  WRITE (iounit, "(T2,A,T37,I4,4I10)") "TDDFPT| ", ispin, moc(1, ispin), moc(2, ispin), &
197  mvt(2, ispin), moc(2, ispin) + mvt(2, ispin)
198  END DO
199  END IF
200 
201  IF (ASSOCIATED(evals_virt)) THEN
202  DO ispin = 1, SIZE(evals_virt)
203  IF (ASSOCIATED(evals_virt(ispin)%array)) DEALLOCATE (evals_virt(ispin)%array)
204  END DO
205  DEALLOCATE (evals_virt)
206  END IF
207 
208  CALL cp_fm_release(mos_virt)
209 
210  CALL timestop(handle)
211 
212  END SUBROUTINE tddfpt_init_mos
213 
214 ! **************************************************************************************************
215 !> \brief Generate all virtual molecular orbitals for a given spin by diagonalising
216 !> the corresponding Kohn-Sham matrix.
217 !> \param gs_mos structure to store occupied and virtual molecular orbitals
218 !> (allocated and initialised on exit)
219 !> \param mo_set ground state molecular orbitals for a given spin
220 !> \param nlumo number of unoccupied states to consider (-1 means all states)
221 !> \param blacs_env BLACS parallel environment
222 !> \param cholesky_method Cholesky method to compute the inverse overlap matrix
223 !> \param matrix_ks Kohn-Sham matrix for a given spin
224 !> \param matrix_s overlap matrix
225 !> \param mos_virt precomputed (OT) expansion coefficients of virtual molecular orbitals
226 !> (in addition to the ADDED_MOS, if present). NULL when no OT is in use.
227 !> \param evals_virt orbital energies of precomputed (OT) virtual molecular orbitals.
228 !> NULL when no OT is in use.
229 !> \param qs_env ...
230 !> \par History
231 !> * 05.2016 created as tddfpt_lumos() [Sergey Chulkov]
232 !> * 06.2016 renamed, altered prototype [Sergey Chulkov]
233 !> * 04.2019 limit the number of unoccupied states, orbital energy correction [Sergey Chulkov]
234 ! **************************************************************************************************
235  SUBROUTINE tddfpt_init_ground_state_mos(gs_mos, mo_set, nlumo, blacs_env, cholesky_method, matrix_ks, matrix_s, &
236  mos_virt, evals_virt, qs_env)
237  TYPE(tddfpt_ground_state_mos) :: gs_mos
238  TYPE(mo_set_type), INTENT(IN) :: mo_set
239  INTEGER, INTENT(in) :: nlumo
240  TYPE(cp_blacs_env_type), POINTER :: blacs_env
241  INTEGER, INTENT(in) :: cholesky_method
242  TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_s
243  TYPE(cp_fm_type), INTENT(IN), POINTER :: mos_virt
244  REAL(kind=dp), DIMENSION(:), POINTER :: evals_virt
245  TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
246 
247  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_init_ground_state_mos'
248  REAL(kind=dp), PARAMETER :: eps_dp = epsilon(0.0_dp)
249 
250  INTEGER :: cholesky_method_inout, handle, icol_global, icol_local, imo, iounit, irow_global, &
251  irow_local, nao, ncol_local, nelectrons, nmo_occ, nmo_scf, nmo_virt, nrow_local, sign_int
252  INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, &
253  sum_sign_array
254  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
255  LOGICAL :: do_eigen, print_phases
256  REAL(kind=dp) :: element, maxocc
257  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
258  POINTER :: my_block
259  REAL(kind=dp), DIMENSION(:), POINTER :: mo_evals_extended, mo_occ_extended, &
260  mo_occ_scf
261  TYPE(cp_fm_struct_type), POINTER :: ao_ao_fm_struct, ao_mo_occ_fm_struct, &
262  ao_mo_virt_fm_struct, wfn_fm_struct
263  TYPE(cp_fm_type) :: matrix_ks_fm, ortho_fm, work_fm, &
264  work_fm_virt
265  TYPE(cp_fm_type), POINTER :: mo_coeff_extended
266  TYPE(cp_logger_type), POINTER :: logger
267  TYPE(mo_set_type), POINTER :: mos_extended
268  TYPE(mp_para_env_type), POINTER :: para_env
269  TYPE(section_vals_type), POINTER :: print_section
270 
271  CALL timeset(routinen, handle)
272 
273  NULLIFY (logger)
274  logger => cp_get_default_logger()
275  iounit = cp_logger_get_default_io_unit(logger)
276 
277  CALL blacs_env%get(para_env=para_env)
278 
279  CALL get_mo_set(mo_set, nao=nao, nmo=nmo_scf, homo=nmo_occ, maxocc=maxocc, &
280  nelectron=nelectrons, occupation_numbers=mo_occ_scf)
281 
282  print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
283  CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases)
284 
285  nmo_virt = nao - nmo_occ
286  IF (nlumo >= 0) &
287  nmo_virt = min(nmo_virt, nlumo)
288 
289  IF (nmo_virt <= 0) &
290  CALL cp_abort(__location__, &
291  'At least one unoccupied molecular orbital is required to calculate excited states.')
292 
293  do_eigen = .false.
294  ! diagonalise the Kohn-Sham matrix one more time if the number of available unoccupied states are too small
295  IF (ASSOCIATED(evals_virt)) THEN
296  cpassert(ASSOCIATED(mos_virt))
297  IF (nmo_virt > nmo_scf - nmo_occ + SIZE(evals_virt)) do_eigen = .true.
298  ELSE
299  IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .true.
300  END IF
301 
302  ! ++ allocate storage space for gs_mos
303  NULLIFY (ao_mo_occ_fm_struct, ao_mo_virt_fm_struct)
304  CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
305  CALL cp_fm_struct_create(ao_mo_virt_fm_struct, nrow_global=nao, ncol_global=nmo_virt, context=blacs_env)
306 
307  NULLIFY (gs_mos%mos_occ, gs_mos%mos_virt, gs_mos%evals_occ_matrix)
308  ALLOCATE (gs_mos%mos_occ, gs_mos%mos_virt)
309  CALL cp_fm_create(gs_mos%mos_occ, ao_mo_occ_fm_struct)
310  CALL cp_fm_create(gs_mos%mos_virt, ao_mo_virt_fm_struct)
311 
312  ALLOCATE (gs_mos%evals_occ(nmo_occ))
313  ALLOCATE (gs_mos%evals_virt(nmo_virt))
314  ALLOCATE (gs_mos%phases_occ(nmo_occ))
315  ALLOCATE (gs_mos%phases_virt(nmo_virt))
316 
317  ! ++ nullify pointers
318  NULLIFY (ao_ao_fm_struct, wfn_fm_struct)
319  NULLIFY (mos_extended, mo_coeff_extended, mo_evals_extended, mo_occ_extended)
320  CALL cp_fm_struct_create(ao_ao_fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
321 
322  IF (do_eigen) THEN
323  ! ++ set of molecular orbitals
324  CALL cp_fm_struct_create(wfn_fm_struct, nrow_global=nao, ncol_global=nmo_occ + nmo_virt, context=blacs_env)
325  ALLOCATE (mos_extended)
326  CALL allocate_mo_set(mos_extended, nao, nmo_occ + nmo_virt, nelectrons, &
327  REAL(nelectrons, dp), maxocc, flexible_electron_count=0.0_dp)
328  CALL init_mo_set(mos_extended, fm_struct=wfn_fm_struct, name="mos-extended")
329  CALL cp_fm_struct_release(wfn_fm_struct)
330  CALL get_mo_set(mos_extended, mo_coeff=mo_coeff_extended, &
331  eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
332 
333  ! use the explicit loop in order to avoid temporary arrays.
334  !
335  ! The assignment statement : mo_occ_extended(1:nmo_scf) = mo_occ_scf(1:nmo_scf)
336  ! implies temporary arrays as a compiler does not know in advance that the pointers
337  ! on both sides of the statement point to non-overlapped memory regions
338  DO imo = 1, nmo_scf
339  mo_occ_extended(imo) = mo_occ_scf(imo)
340  END DO
341  mo_occ_extended(nmo_scf + 1:) = 0.0_dp
342 
343  ! ++ allocate temporary matrices
344  CALL cp_fm_create(matrix_ks_fm, ao_ao_fm_struct)
345  CALL cp_fm_create(ortho_fm, ao_ao_fm_struct)
346  CALL cp_fm_create(work_fm, ao_ao_fm_struct)
347  CALL cp_fm_struct_release(ao_ao_fm_struct)
348 
349  ! some stuff from the subroutine general_eigenproblem()
350  CALL copy_dbcsr_to_fm(matrix_s, ortho_fm)
351  CALL copy_dbcsr_to_fm(matrix_ks, matrix_ks_fm)
352 
353  IF (cholesky_method == cholesky_dbcsr) THEN
354  cpabort('CHOLESKY DBCSR_INVERSE is not implemented in TDDFT.')
355  ELSE IF (cholesky_method == cholesky_off) THEN
356  cpabort('CHOLESKY OFF is not implemented in TDDFT.')
357  ELSE
358  CALL cp_fm_cholesky_decompose(ortho_fm)
359  IF (cholesky_method == cholesky_inverse) THEN
360  CALL cp_fm_triangular_invert(ortho_fm)
361  END IF
362 
363  ! need to store 'cholesky_method' in a temporary variable, as the subroutine eigensolver()
364  ! will update this variable
365  cholesky_method_inout = cholesky_method
366  CALL eigensolver(matrix_ks_fm=matrix_ks_fm, mo_set=mos_extended, ortho=ortho_fm, &
367  work=work_fm, cholesky_method=cholesky_method_inout, &
368  do_level_shift=.false., level_shift=0.0_dp, use_jacobi=.false.)
369  END IF
370 
371  ! -- clean up needless matrices
372  CALL cp_fm_release(work_fm)
373  CALL cp_fm_release(ortho_fm)
374  CALL cp_fm_release(matrix_ks_fm)
375  ELSE
376  CALL get_mo_set(mo_set, mo_coeff=mo_coeff_extended, &
377  eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
378  END IF
379 
380  ! compute the phase of molecular orbitals;
381  ! matrix work_fm holds occupied molecular orbital coefficients distributed among all the processors
382  !CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
383  CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
384  CALL cp_fm_struct_release(ao_mo_occ_fm_struct)
385 
386  CALL cp_fm_to_fm(mo_coeff_extended, work_fm, ncol=nmo_occ, source_start=1, target_start=1)
387  CALL cp_fm_get_info(work_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
388  row_indices=row_indices, col_indices=col_indices, local_data=my_block)
389 
390  ALLOCATE (minrow_neg_array(nmo_occ), minrow_pos_array(nmo_occ), sum_sign_array(nmo_occ))
391  minrow_neg_array(:) = nao
392  minrow_pos_array(:) = nao
393  sum_sign_array(:) = 0
394  DO icol_local = 1, ncol_local
395  icol_global = col_indices(icol_local)
396 
397  DO irow_local = 1, nrow_local
398  element = my_block(irow_local, icol_local)
399 
400  sign_int = 0
401  IF (element >= eps_dp) THEN
402  sign_int = 1
403  ELSE IF (element <= -eps_dp) THEN
404  sign_int = -1
405  END IF
406 
407  sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
408 
409  irow_global = row_indices(irow_local)
410  IF (sign_int > 0) THEN
411  IF (minrow_pos_array(icol_global) > irow_global) &
412  minrow_pos_array(icol_global) = irow_global
413  ELSE IF (sign_int < 0) THEN
414  IF (minrow_neg_array(icol_global) > irow_global) &
415  minrow_neg_array(icol_global) = irow_global
416  END IF
417  END DO
418  END DO
419 
420  CALL para_env%sum(sum_sign_array)
421  CALL para_env%min(minrow_neg_array)
422  CALL para_env%min(minrow_pos_array)
423 
424  DO icol_local = 1, nmo_occ
425  IF (sum_sign_array(icol_local) > 0) THEN
426  ! most of the expansion coefficients are positive => MO's phase = +1
427  gs_mos%phases_occ(icol_local) = 1.0_dp
428  ELSE IF (sum_sign_array(icol_local) < 0) THEN
429  ! most of the expansion coefficients are negative => MO's phase = -1
430  gs_mos%phases_occ(icol_local) = -1.0_dp
431  ELSE
432  ! equal number of positive and negative expansion coefficients
433  IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local)) THEN
434  ! the first positive expansion coefficient has a lower index then
435  ! the first negative expansion coefficient; MO's phase = +1
436  gs_mos%phases_occ(icol_local) = 1.0_dp
437  ELSE
438  ! MO's phase = -1
439  gs_mos%phases_occ(icol_local) = -1.0_dp
440  END IF
441  END IF
442  END DO
443 
444  DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
445 
446  ! return the requested occupied and virtual molecular orbitals and corresponding orbital energies
447  CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_occ, ncol=nmo_occ, source_start=1, target_start=1)
448  gs_mos%evals_occ(1:nmo_occ) = mo_evals_extended(1:nmo_occ)
449 
450  IF (ASSOCIATED(evals_virt) .AND. (.NOT. do_eigen) .AND. nmo_virt > nmo_scf - nmo_occ) THEN
451  CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_scf - nmo_occ, &
452  source_start=nmo_occ + 1, target_start=1)
453  CALL cp_fm_to_fm(mos_virt, gs_mos%mos_virt, ncol=nmo_virt - (nmo_scf - nmo_occ), &
454  source_start=1, target_start=nmo_scf - nmo_occ + 1)
455  gs_mos%evals_virt(1:nmo_scf - nmo_occ) = evals_virt(nmo_occ + 1:nmo_occ + nmo_scf)
456  gs_mos%evals_virt(nmo_scf - nmo_occ + 1:nmo_virt) = evals_virt(1:nmo_virt - (nmo_scf - nmo_occ))
457  ELSE
458  CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_virt, source_start=nmo_occ + 1, target_start=1)
459  gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt)
460  END IF
461 
462  IF (print_phases) THEN
463  ! compute the phase of molecular orbitals;
464  ! matrix work_fm holds virtual molecular orbital coefficients distributed among all the processors
465  !CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
466  CALL cp_fm_create(work_fm_virt, ao_mo_virt_fm_struct)
467 
468  CALL cp_fm_to_fm(gs_mos%mos_virt, work_fm_virt, ncol=nmo_virt, source_start=1, target_start=1)
469  CALL cp_fm_get_info(work_fm_virt, nrow_local=nrow_local, ncol_local=ncol_local, &
470  row_indices=row_indices, col_indices=col_indices, local_data=my_block)
471 
472  ALLOCATE (minrow_neg_array(nmo_virt), minrow_pos_array(nmo_virt), sum_sign_array(nmo_virt))
473  minrow_neg_array(:) = nao
474  minrow_pos_array(:) = nao
475  sum_sign_array(:) = 0
476  DO icol_local = 1, ncol_local
477  icol_global = col_indices(icol_local)
478 
479  DO irow_local = 1, nrow_local
480  element = my_block(irow_local, icol_local)
481 
482  sign_int = 0
483  IF (element >= eps_dp) THEN
484  sign_int = 1
485  ELSE IF (element <= -eps_dp) THEN
486  sign_int = -1
487  END IF
488 
489  sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
490 
491  irow_global = row_indices(irow_local)
492  IF (sign_int > 0) THEN
493  IF (minrow_pos_array(icol_global) > irow_global) &
494  minrow_pos_array(icol_global) = irow_global
495  ELSE IF (sign_int < 0) THEN
496  IF (minrow_neg_array(icol_global) > irow_global) &
497  minrow_neg_array(icol_global) = irow_global
498  END IF
499  END DO
500  END DO
501 
502  CALL para_env%sum(sum_sign_array)
503  CALL para_env%min(minrow_neg_array)
504  CALL para_env%min(minrow_pos_array)
505  DO icol_local = 1, nmo_virt
506  IF (sum_sign_array(icol_local) > 0) THEN
507  ! most of the expansion coefficients are positive => MO's phase = +1
508  gs_mos%phases_virt(icol_local) = 1.0_dp
509  ELSE IF (sum_sign_array(icol_local) < 0) THEN
510  ! most of the expansion coefficients are negative => MO's phase = -1
511  gs_mos%phases_virt(icol_local) = -1.0_dp
512  ELSE
513  ! equal number of positive and negative expansion coefficients
514  IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local)) THEN
515  ! the first positive expansion coefficient has a lower index then
516  ! the first negative expansion coefficient; MO's phase = +1
517  gs_mos%phases_virt(icol_local) = 1.0_dp
518  ELSE
519  ! MO's phase = -1
520  gs_mos%phases_virt(icol_local) = -1.0_dp
521  END IF
522  END IF
523  END DO
524  DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
525  CALL cp_fm_release(work_fm_virt)
526  END IF !print_phases
527  CALL cp_fm_struct_release(ao_mo_virt_fm_struct) ! here after print_phases
528 
529  CALL cp_fm_release(work_fm)
530 
531  IF (do_eigen) THEN
532  CALL deallocate_mo_set(mos_extended)
533  DEALLOCATE (mos_extended)
534  END IF
535 
536  CALL timestop(handle)
537 
538  END SUBROUTINE tddfpt_init_ground_state_mos
539 
540 ! **************************************************************************************************
541 !> \brief Release molecular orbitals.
542 !> \param gs_mos structure that holds occupied and virtual molecular orbitals
543 !> \par History
544 !> * 06.2016 created [Sergey Chulkov]
545 ! **************************************************************************************************
546  SUBROUTINE tddfpt_release_ground_state_mos(gs_mos)
547  TYPE(tddfpt_ground_state_mos), INTENT(inout) :: gs_mos
548 
549  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_release_ground_state_mos'
550 
551  INTEGER :: handle
552 
553  CALL timeset(routinen, handle)
554 
555  IF (ALLOCATED(gs_mos%phases_occ)) &
556  DEALLOCATE (gs_mos%phases_occ)
557 
558  IF (ALLOCATED(gs_mos%evals_virt)) &
559  DEALLOCATE (gs_mos%evals_virt)
560 
561  IF (ALLOCATED(gs_mos%evals_occ)) &
562  DEALLOCATE (gs_mos%evals_occ)
563 
564  IF (ALLOCATED(gs_mos%phases_virt)) &
565  DEALLOCATE (gs_mos%phases_virt)
566 
567  IF (ASSOCIATED(gs_mos%evals_occ_matrix)) THEN
568  CALL cp_fm_release(gs_mos%evals_occ_matrix)
569  DEALLOCATE (gs_mos%evals_occ_matrix)
570  END IF
571 
572  IF (ASSOCIATED(gs_mos%mos_virt)) THEN
573  CALL cp_fm_release(gs_mos%mos_virt)
574  DEALLOCATE (gs_mos%mos_virt)
575  END IF
576 
577  IF (ASSOCIATED(gs_mos%mos_occ)) THEN
578  CALL cp_fm_release(gs_mos%mos_occ)
579  DEALLOCATE (gs_mos%mos_occ)
580  END IF
581 
582  CALL timestop(handle)
583  END SUBROUTINE tddfpt_release_ground_state_mos
584 
585 ! **************************************************************************************************
586 !> \brief Callculate orbital corrected KS matrix for TDDFPT
587 !> \param qs_env Quickstep environment
588 !> \param gs_mos ...
589 !> \param matrix_ks_oep ...
590 ! **************************************************************************************************
591  SUBROUTINE tddfpt_oecorr(qs_env, gs_mos, matrix_ks_oep)
592  TYPE(qs_environment_type), POINTER :: qs_env
593  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
594  POINTER :: gs_mos
595  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_oep
596 
597  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_oecorr'
598 
599  INTEGER :: handle, iounit, ispin, nao, nmo_occ, &
600  nspins
601  LOGICAL :: do_hfx
602  TYPE(cp_blacs_env_type), POINTER :: blacs_env
603  TYPE(cp_fm_struct_type), POINTER :: ao_mo_occ_fm_struct, &
604  mo_occ_mo_occ_fm_struct
605  TYPE(cp_fm_type) :: work_fm
606  TYPE(cp_logger_type), POINTER :: logger
607  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
608  TYPE(dft_control_type), POINTER :: dft_control
609  TYPE(section_vals_type), POINTER :: hfx_section, xc_fun_empty, &
610  xc_fun_original
611  TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
612 
613  CALL timeset(routinen, handle)
614 
615  NULLIFY (logger)
616  logger => cp_get_default_logger()
617  iounit = cp_logger_get_default_io_unit(logger)
618 
619  CALL get_qs_env(qs_env, blacs_env=blacs_env, dft_control=dft_control, matrix_ks=matrix_ks)
620  tddfpt_control => dft_control%tddfpt2_control
621 
622  ! obtain corrected KS-matrix
623  ! We should 'save' the energy values?
624  nspins = SIZE(gs_mos)
625  NULLIFY (matrix_ks_oep)
626  IF (tddfpt_control%oe_corr /= oe_none) THEN
627  IF (iounit > 0) THEN
628  WRITE (iounit, "(1X,A)") "", &
629  "-------------------------------------------------------------------------------", &
630  "- Orbital Eigenvalue Correction Started -", &
631  "-------------------------------------------------------------------------------"
632  END IF
633 
634  CALL cp_warn(__location__, &
635  "Orbital energy correction potential is an experimental feature. "// &
636  "Use it with extreme care")
637 
638  hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
639  CALL section_vals_get(hfx_section, explicit=do_hfx)
640  IF (do_hfx) THEN
641  CALL cp_abort(__location__, &
642  "Implementation of orbital energy correction XC-potentials is "// &
643  "currently incompatible with exact-exchange functionals")
644  END IF
645 
646  CALL dbcsr_allocate_matrix_set(matrix_ks_oep, nspins)
647  DO ispin = 1, nspins
648  CALL dbcsr_init_p(matrix_ks_oep(ispin)%matrix)
649  CALL dbcsr_copy(matrix_ks_oep(ispin)%matrix, matrix_ks(ispin)%matrix)
650  END DO
651 
652  ! KS-matrix without XC-terms
653  xc_fun_original => section_vals_get_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL")
654  CALL section_vals_retain(xc_fun_original)
655  NULLIFY (xc_fun_empty)
656  CALL section_vals_create(xc_fun_empty, xc_fun_original%section)
657  CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_empty)
658  CALL section_vals_release(xc_fun_empty)
659 
660  IF (dft_control%qs_control%semi_empirical) THEN
661  cpabort("TDDFPT with SE not possible")
662  ELSEIF (dft_control%qs_control%dftb) THEN
663  cpabort("TDDFPT with DFTB not possible")
664  ELSEIF (dft_control%qs_control%xtb) THEN
665  CALL build_xtb_ks_matrix(qs_env, calculate_forces=.false., just_energy=.false., &
666  ext_ks_matrix=matrix_ks_oep)
667  ELSE
668  CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.false., just_energy=.false., &
669  ext_ks_matrix=matrix_ks_oep)
670  END IF
671 
672  IF (tddfpt_control%oe_corr == oe_saop .OR. &
673  tddfpt_control%oe_corr == oe_lb .OR. &
674  tddfpt_control%oe_corr == oe_gllb) THEN
675  IF (iounit > 0) THEN
676  WRITE (iounit, "(T2,A)") " Orbital energy correction of SAOP type "
677  END IF
678  CALL add_saop_pot(matrix_ks_oep, qs_env, tddfpt_control%oe_corr)
679  ELSE IF (tddfpt_control%oe_corr == oe_shift) THEN
680  IF (iounit > 0) THEN
681  WRITE (iounit, "(T2,A,T71,F10.3)") &
682  " Virtual Orbital Eigenvalue Shift [eV] ", tddfpt_control%ev_shift*evolt
683  WRITE (iounit, "(T2,A,T71,F10.3)") &
684  " Open Shell Orbital Eigenvalue Shift [eV] ", tddfpt_control%eos_shift*evolt
685  END IF
686  CALL ev_shift_operator(qs_env, gs_mos, matrix_ks_oep, &
687  tddfpt_control%ev_shift, tddfpt_control%eos_shift)
688  ELSE
689  CALL cp_abort(__location__, &
690  "Unimplemented orbital energy correction potential")
691  END IF
692  CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_original)
693  CALL section_vals_release(xc_fun_original)
694 
695  ! compute 'evals_occ_matrix'
696  CALL dbcsr_get_info(matrix_ks(1)%matrix, nfullrows_total=nao)
697  NULLIFY (mo_occ_mo_occ_fm_struct)
698  DO ispin = 1, nspins
699  nmo_occ = SIZE(gs_mos(ispin)%evals_occ)
700  CALL cp_fm_struct_create(mo_occ_mo_occ_fm_struct, nrow_global=nmo_occ, ncol_global=nmo_occ, &
701  context=blacs_env)
702  ALLOCATE (gs_mos(ispin)%evals_occ_matrix)
703  CALL cp_fm_create(gs_mos(ispin)%evals_occ_matrix, mo_occ_mo_occ_fm_struct)
704  CALL cp_fm_struct_release(mo_occ_mo_occ_fm_struct)
705  ! work_fm is a temporary [nao x nmo_occ] matrix
706  CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, &
707  context=blacs_env)
708  CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
709  CALL cp_fm_struct_release(ao_mo_occ_fm_struct)
710  CALL cp_dbcsr_sm_fm_multiply(matrix_ks_oep(ispin)%matrix, gs_mos(ispin)%mos_occ, &
711  work_fm, ncol=nmo_occ, alpha=1.0_dp, beta=0.0_dp)
712  CALL parallel_gemm('T', 'N', nmo_occ, nmo_occ, nao, 1.0_dp, gs_mos(ispin)%mos_occ, work_fm, &
713  0.0_dp, gs_mos(ispin)%evals_occ_matrix)
714  CALL cp_fm_release(work_fm)
715  END DO
716  IF (iounit > 0) THEN
717  WRITE (iounit, "(1X,A)") &
718  "-------------------------------------------------------------------------------"
719  END IF
720 
721  END IF
722 
723  CALL timestop(handle)
724 
725  END SUBROUTINE tddfpt_oecorr
726 
727 ! **************************************************************************************************
728 !> \brief Compute the number of possible singly excited states (occ -> virt)
729 !> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state
730 !> \return the number of possible single excitations
731 !> \par History
732 !> * 01.2017 created [Sergey Chulkov]
733 ! **************************************************************************************************
734  PURE FUNCTION tddfpt_total_number_of_states(gs_mos) RESULT(nstates_total)
735  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
736  INTENT(in) :: gs_mos
737  INTEGER(kind=int_8) :: nstates_total
738 
739  INTEGER :: ispin, nspins
740 
741  nstates_total = 0
742  nspins = SIZE(gs_mos)
743 
744  DO ispin = 1, nspins
745  nstates_total = nstates_total + &
746  SIZE(gs_mos(ispin)%evals_occ, kind=int_8)* &
747  SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
748  END DO
749  END FUNCTION tddfpt_total_number_of_states
750 
751 ! **************************************************************************************************
752 !> \brief Create a shift operator on virtual/open shell space
753 !> Shift operator = Edelta*Q Q: projector on virtual space (1-PS)
754 !> projector on open shell space PosS
755 !> \param qs_env the qs_env that is perturbed by this p_env
756 !> \param gs_mos ...
757 !> \param matrix_ks ...
758 !> \param ev_shift ...
759 !> \param eos_shift ...
760 !> \par History
761 !> 02.04.2019 adapted for TDDFT use from p_env (JGH)
762 !> \author JGH
763 ! **************************************************************************************************
764  SUBROUTINE ev_shift_operator(qs_env, gs_mos, matrix_ks, ev_shift, eos_shift)
765 
766  TYPE(qs_environment_type), POINTER :: qs_env
767  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
768  POINTER :: gs_mos
769  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
770  REAL(kind=dp), INTENT(IN) :: ev_shift, eos_shift
771 
772  CHARACTER(len=*), PARAMETER :: routinen = 'ev_shift_operator'
773 
774  INTEGER :: handle, ispin, n_spins, na, nb, nhomo, &
775  nl, nos, nrow, nu, nvirt
776  TYPE(cp_fm_struct_type), POINTER :: fmstruct
777  TYPE(cp_fm_type) :: cmos, cvec
778  TYPE(cp_fm_type), POINTER :: coeff
779  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
780  TYPE(dbcsr_type), POINTER :: smat
781  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
782 
783  CALL timeset(routinen, handle)
784 
785  n_spins = SIZE(gs_mos)
786  cpassert(n_spins == SIZE(matrix_ks))
787 
788  IF (eos_shift /= 0.0_dp .AND. n_spins > 1) THEN
789  cpabort("eos_shift not implemented")
790  CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
791  smat => matrix_s(1)%matrix
792  CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
793  CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
794  nl = min(na, nb)
795  nu = max(na, nb)
796  ! open shell orbital shift
797  DO ispin = 1, n_spins
798  coeff => gs_mos(ispin)%mos_occ
799  CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
800  IF (nhomo == nu) THEN
801  ! downshift with -eos_shift using occupied orbitals
802  nos = nu - nl
803  CALL cp_fm_create(cmos, fmstruct)
804  CALL cp_fm_get_info(coeff, nrow_global=nrow)
805  CALL cp_fm_to_fm_submat(coeff, cmos, nrow, nos, 1, nl + 1, 1, 1)
806  CALL cp_fm_create(cvec, fmstruct)
807  CALL cp_dbcsr_sm_fm_multiply(smat, cmos, cvec, nos, 1.0_dp, 0.0_dp)
808  CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
809  alpha=-eos_shift, keep_sparsity=.true.)
810  CALL cp_fm_release(cmos)
811  CALL cp_fm_release(cvec)
812  ELSE
813  ! upshift with eos_shift using virtual orbitals
814  coeff => gs_mos(ispin)%mos_virt
815  CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
816  nos = nu - nhomo
817  cpassert(nvirt >= nos)
818  CALL cp_fm_create(cvec, fmstruct)
819  CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nos, 1.0_dp, 0.0_dp)
820  CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
821  alpha=eos_shift, keep_sparsity=.true.)
822  CALL cp_fm_release(cvec)
823  END IF
824  END DO
825  ! virtual shift
826  IF (ev_shift /= 0.0_dp) THEN
827  DO ispin = 1, n_spins
828  CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
829  alpha_scalar=1.0_dp, beta_scalar=ev_shift)
830  coeff => gs_mos(ispin)%mos_occ
831  CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
832  CALL cp_fm_create(cvec, fmstruct)
833  CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nhomo, 1.0_dp, 0.0_dp)
834  CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nhomo, &
835  alpha=-ev_shift, keep_sparsity=.true.)
836  CALL cp_fm_release(cvec)
837  IF (nhomo < nu) THEN
838  nos = nu - nhomo
839  coeff => gs_mos(ispin)%mos_virt
840  CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
841  cpassert(nvirt >= nos)
842  CALL cp_fm_create(cvec, fmstruct)
843  CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nos, 1.0_dp, 0.0_dp)
844  CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
845  alpha=-ev_shift, keep_sparsity=.true.)
846  CALL cp_fm_release(cvec)
847  END IF
848  END DO
849  END IF
850  ELSE
851  ! virtual shift
852  IF (ev_shift /= 0.0_dp) THEN
853  CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
854  smat => matrix_s(1)%matrix
855  DO ispin = 1, n_spins
856  CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
857  alpha_scalar=1.0_dp, beta_scalar=ev_shift)
858  coeff => gs_mos(ispin)%mos_occ
859  CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
860  CALL cp_fm_create(cvec, fmstruct)
861  CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nhomo, 1.0_dp, 0.0_dp)
862  CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nhomo, &
863  alpha=-ev_shift, keep_sparsity=.true.)
864  CALL cp_fm_release(cvec)
865  END DO
866  END IF
867  END IF
868  ! set eigenvalues
869  IF (eos_shift == 0.0_dp .OR. n_spins == 1) THEN
870  DO ispin = 1, n_spins
871  IF (ALLOCATED(gs_mos(ispin)%evals_virt)) THEN
872  gs_mos(ispin)%evals_virt(:) = gs_mos(ispin)%evals_virt(:) + ev_shift
873  END IF
874  END DO
875  ELSE
876  CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
877  CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
878  nl = min(na, nb)
879  nu = max(na, nb)
880  nos = nu - nl
881  IF (na == nu) THEN
882  IF (ALLOCATED(gs_mos(1)%evals_occ)) THEN
883  gs_mos(1)%evals_occ(nl + 1:nu) = gs_mos(1)%evals_occ(nl + 1:nu) - eos_shift
884  END IF
885  IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
886  gs_mos(1)%evals_virt(:) = gs_mos(1)%evals_virt(:) + ev_shift
887  END IF
888  IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
889  gs_mos(2)%evals_virt(1:nos) = gs_mos(2)%evals_virt(1:nos) + eos_shift
890  gs_mos(2)%evals_virt(nos + 1:) = gs_mos(2)%evals_virt(nos + 1:) + ev_shift
891  END IF
892  ELSE
893  IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
894  gs_mos(1)%evals_virt(1:nos) = gs_mos(1)%evals_virt(1:nos) + eos_shift
895  gs_mos(1)%evals_virt(nos + 1:) = gs_mos(1)%evals_virt(nos + 1:) + ev_shift
896  END IF
897  IF (ALLOCATED(gs_mos(2)%evals_occ)) THEN
898  gs_mos(2)%evals_occ(nl + 1:nu) = gs_mos(2)%evals_occ(nl + 1:nu) - eos_shift
899  END IF
900  IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
901  gs_mos(2)%evals_virt(:) = gs_mos(2)%evals_virt(:) + ev_shift
902  END IF
903  END IF
904  END IF
905 
906  CALL timestop(handle)
907 
908  END SUBROUTINE ev_shift_operator
909 
910 ! **************************************************************************************************
911 !> \brief Generate missed guess vectors.
912 !> \param evects guess vectors distributed across all processors (initialised on exit)
913 !> \param evals guessed transition energies (initialised on exit)
914 !> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state
915 !> \param log_unit output unit
916 !> \par History
917 !> * 05.2016 created as tddfpt_guess() [Sergey Chulkov]
918 !> * 06.2016 renamed, altered prototype, supports spin-polarised density [Sergey Chulkov]
919 !> * 01.2017 simplified prototype, do not compute all possible singly-excited states
920 !> [Sergey Chulkov]
921 !> \note \parblock
922 !> Based on the subroutine co_initial_guess() which was originally created by
923 !> Thomas Chassaing on 06.2003.
924 !>
925 !> Only not associated guess vectors 'evects(spin, state)%matrix' are allocated and
926 !> initialised; associated vectors assumed to be initialised elsewhere (e.g. using
927 !> a restart file).
928 !> \endparblock
929 ! **************************************************************************************************
930  SUBROUTINE tddfpt_guess_vectors(evects, evals, gs_mos, log_unit)
931  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout) :: evects
932  REAL(kind=dp), DIMENSION(:), INTENT(inout) :: evals
933  TYPE(tddfpt_ground_state_mos), &
934  DIMENSION(SIZE(evects, 1)), INTENT(in) :: gs_mos
935  INTEGER, INTENT(in) :: log_unit
936 
937  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_guess_vectors'
938 
939  CHARACTER(len=5) :: spin_label
940  INTEGER :: handle, imo_occ, imo_virt, ind, ispin, &
941  istate, jspin, nspins, nstates, &
942  nstates_occ_virt_alpha, &
943  nstates_selected
944  INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
945  INTEGER, DIMENSION(maxspins) :: nmo_occ_avail, nmo_occ_selected, &
946  nmo_virt_selected
947  REAL(kind=dp) :: e_occ
948  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: e_virt_minus_occ
949  TYPE(cp_fm_struct_p_type), DIMENSION(maxspins) :: fm_struct_evects
950 
951  CALL timeset(routinen, handle)
952 
953  nspins = SIZE(evects, 1)
954  nstates = SIZE(evects, 2)
955 
956  IF (debug_this_module) THEN
957  cpassert(nstates > 0)
958  cpassert(nspins == 1 .OR. nspins == 2)
959  END IF
960 
961  DO ispin = 1, nspins
962  ! number of occupied orbitals for each spin component
963  nmo_occ_avail(ispin) = SIZE(gs_mos(ispin)%evals_occ)
964  ! number of occupied and virtual orbitals which can potentially
965  ! contribute to the excited states in question.
966  nmo_occ_selected(ispin) = min(nmo_occ_avail(ispin), nstates)
967  nmo_virt_selected(ispin) = min(SIZE(gs_mos(ispin)%evals_virt), nstates)
968 
969  CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct_evects(ispin)%struct)
970  END DO
971 
972  ! TO DO: the variable 'nstates_selected' should probably be declared as INTEGER(kind=int_8),
973  ! however we need a special version of the subroutine sort() in order to do so
974  nstates_selected = dot_product(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins))
975 
976  ALLOCATE (inds(nstates_selected))
977  ALLOCATE (e_virt_minus_occ(nstates_selected))
978 
979  istate = 0
980  DO ispin = 1, nspins
981  DO imo_occ = 1, nmo_occ_selected(ispin)
982  ! Here imo_occ enumerate Occupied orbitals in inverse order (from the last to the first element)
983  e_occ = gs_mos(ispin)%evals_occ(nmo_occ_avail(ispin) - imo_occ + 1)
984 
985  DO imo_virt = 1, nmo_virt_selected(ispin)
986  istate = istate + 1
987  e_virt_minus_occ(istate) = gs_mos(ispin)%evals_virt(imo_virt) - e_occ
988  END DO
989  END DO
990  END DO
991 
992  IF (debug_this_module) THEN
993  cpassert(istate == nstates_selected)
994  END IF
995 
996  CALL sort(e_virt_minus_occ, nstates_selected, inds)
997 
998  IF (nspins == 1) THEN
999  ispin = 1
1000  spin_label = ' '
1001  END IF
1002 
1003  nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1)
1004  IF (log_unit > 0) THEN
1005  WRITE (log_unit, "(1X,A)") "", &
1006  "-------------------------------------------------------------------------------", &
1007  "- TDDFPT Initial Guess -", &
1008  "-------------------------------------------------------------------------------"
1009  WRITE (log_unit, '(T11,A)') "State Occupied -> Virtual Excitation"
1010  WRITE (log_unit, '(T11,A)') "number orbital orbital energy (eV)"
1011  WRITE (log_unit, '(1X,79("-"))')
1012  END IF
1013 
1014  DO istate = 1, nstates
1015  IF (ASSOCIATED(evects(1, istate)%matrix_struct)) THEN
1016  IF (log_unit > 0) &
1017  WRITE (log_unit, '(T7,I8,T28,A19,T60,F14.5)') &
1018  istate, "*** restarted ***", evals(istate)*evolt
1019  ELSE
1020  ind = inds(istate) - 1
1021  IF (nspins > 1) THEN
1022  IF (ind < nstates_occ_virt_alpha) THEN
1023  ispin = 1
1024  spin_label = '(alp)'
1025  ELSE
1026  ispin = 2
1027  ind = ind - nstates_occ_virt_alpha
1028  spin_label = '(bet)'
1029  END IF
1030  END IF
1031 
1032  imo_occ = nmo_occ_avail(ispin) - ind/nmo_virt_selected(ispin)
1033  imo_virt = mod(ind, nmo_virt_selected(ispin)) + 1
1034  evals(istate) = e_virt_minus_occ(istate)
1035 
1036  IF (log_unit > 0) &
1037  WRITE (log_unit, '(T7,I8,T24,I8,T37,A5,T45,I8,T54,A5,T60,F14.5)') &
1038  istate, imo_occ, spin_label, nmo_occ_avail(ispin) + imo_virt, spin_label, e_virt_minus_occ(istate)*evolt
1039 
1040  DO jspin = 1, nspins
1041  ! .NOT. ASSOCIATED(evects(jspin, istate)%matrix_struct))
1042  CALL cp_fm_create(evects(jspin, istate), fm_struct_evects(jspin)%struct)
1043  CALL cp_fm_set_all(evects(jspin, istate), 0.0_dp)
1044 
1045  IF (jspin == ispin) &
1046  CALL cp_fm_to_fm(gs_mos(ispin)%mos_virt, evects(ispin, istate), &
1047  ncol=1, source_start=imo_virt, target_start=imo_occ)
1048  END DO
1049  END IF
1050  END DO
1051 
1052  IF (log_unit > 0) THEN
1053  WRITE (log_unit, '(/,T7,A,T50,I24)') 'Number of active states:', tddfpt_total_number_of_states(gs_mos)
1054  WRITE (log_unit, "(1X,A)") &
1055  "-------------------------------------------------------------------------------"
1056  END IF
1057 
1058  DEALLOCATE (e_virt_minus_occ)
1059  DEALLOCATE (inds)
1060 
1061  CALL timestop(handle)
1062 
1063  END SUBROUTINE tddfpt_guess_vectors
1064 
1065 END MODULE qs_tddfpt2_utils
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition: cp_array_utils.F:18
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
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,...
Definition: cp_dbcsr_operations.F:887
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_triangular_invert
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
Definition: cp_fm_basic_linalg.F:2193
cp_fm_cholesky
various cholesky decomposition related routines
Definition: cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_decompose
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,...
Definition: cp_fm_cholesky.F:50
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
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
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
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
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_to_fm_submat
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
Definition: cp_fm_types.F:1473
cp_fm_types::cp_fm_set_all
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
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
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...
Definition: cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::oe_saop
integer, parameter, public oe_saop
Definition: input_constants.F:584
input_constants::cholesky_restore
integer, parameter, public cholesky_restore
Definition: input_constants.F:195
input_constants::cholesky_dbcsr
integer, parameter, public cholesky_dbcsr
Definition: input_constants.F:195
input_constants::cholesky_off
integer, parameter, public cholesky_off
Definition: input_constants.F:195
input_constants::oe_none
integer, parameter, public oe_none
Definition: input_constants.F:584
input_constants::oe_shift
integer, parameter, public oe_shift
Definition: input_constants.F:584
input_constants::cholesky_inverse
integer, parameter, public cholesky_inverse
Definition: input_constants.F:195
input_constants::oe_lb
integer, parameter, public oe_lb
Definition: input_constants.F:584
input_constants::oe_gllb
integer, parameter, public oe_gllb
Definition: input_constants.F:584
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_create
recursive subroutine, public section_vals_create(section_vals, section)
creates a object where to store the values of a section
Definition: input_section_types.F:601
input_section_types::section_vals_retain
subroutine, public section_vals_retain(section_vals)
retains the given section values (see doc/ReferenceCounting.html)
Definition: input_section_types.F:631
input_section_types::section_vals_get_subs_vals
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
Definition: input_section_types.F:724
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition: input_section_types.F:697
input_section_types::section_vals_set_subs_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
Definition: input_section_types.F:1679
input_section_types::section_vals_val_get
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
Definition: input_section_types.F:1040
input_section_types::section_vals_release
recursive subroutine, public section_vals_release(section_vals)
releases the given object
Definition: input_section_types.F:646
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition: kinds.F:54
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition: physcon.F:183
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_ks_methods
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_methods.F:22
qs_ks_methods::qs_ks_build_kohn_sham_matrix
subroutine, public qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, print_active, ext_ks_matrix)
routine where the real calculations are made: the KS matrix is calculated
Definition: qs_ks_methods.F:173
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::allocate_mo_set
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...
Definition: qs_mo_types.F:206
qs_mo_types::deallocate_mo_set
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition: qs_mo_types.F:352
qs_mo_types::get_mo_set
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.
Definition: qs_mo_types.F:397
qs_mo_types::init_mo_set
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 ...
Definition: qs_mo_types.F:245
qs_scf_methods
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
Definition: qs_scf_methods.F:19
qs_scf_methods::eigensolver
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....
Definition: qs_scf_methods.F:168
qs_scf_post_gpw
Does all kind of post scf calculations for GPW/GAPW.
Definition: qs_scf_post_gpw.F:15
qs_scf_post_gpw::make_lumo_gpw
subroutine, public make_lumo_gpw(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, nlumo, nlumos)
Gets the lumos, and eigenvalues for the lumos.
Definition: qs_scf_post_gpw.F:853
qs_scf_types
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14
qs_scf_types::ot_method_nr
integer, parameter, public ot_method_nr
Definition: qs_scf_types.F:51
qs_tddfpt2_types
Definition: qs_tddfpt2_types.F:8
qs_tddfpt2_utils
Definition: qs_tddfpt2_utils.F:8
qs_tddfpt2_utils::tddfpt_release_ground_state_mos
subroutine, public tddfpt_release_ground_state_mos(gs_mos)
Release molecular orbitals.
Definition: qs_tddfpt2_utils.F:547
qs_tddfpt2_utils::tddfpt_total_number_of_states
pure integer(kind=int_8) function, public tddfpt_total_number_of_states(gs_mos)
Compute the number of possible singly excited states (occ -> virt)
Definition: qs_tddfpt2_utils.F:735
qs_tddfpt2_utils::tddfpt_init_ground_state_mos
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...
Definition: qs_tddfpt2_utils.F:237
qs_tddfpt2_utils::tddfpt_guess_vectors
subroutine, public tddfpt_guess_vectors(evects, evals, gs_mos, log_unit)
Generate missed guess vectors.
Definition: qs_tddfpt2_utils.F:931
qs_tddfpt2_utils::tddfpt_oecorr
subroutine, public tddfpt_oecorr(qs_env, gs_mos, matrix_ks_oep)
Callculate orbital corrected KS matrix for TDDFPT.
Definition: qs_tddfpt2_utils.F:592
qs_tddfpt2_utils::tddfpt_init_mos
subroutine, public tddfpt_init_mos(qs_env, gs_mos, iounit)
Prepare MOs for TDDFPT Calculations.
Definition: qs_tddfpt2_utils.F:100
util
All kind of helpful little routines.
Definition: util.F:14
xc_pot_saop
Calculate the saop potential.
Definition: xc_pot_saop.F:11
xc_pot_saop::add_saop_pot
subroutine, public add_saop_pot(ks_matrix, qs_env, oe_corr)
...
Definition: xc_pot_saop.F:126
xtb_matrices
Calculation of Overlap and Hamiltonian matrices in xTB Reference: Stefan Grimme, Christoph Bannwarth,...
Definition: xtb_matrices.F:15
xtb_matrices::build_xtb_ks_matrix
subroutine, public build_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
...
Definition: xtb_matrices.F:1521