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