Line data Source code
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_api, ONLY: dbcsr_add,&
15 : dbcsr_copy,&
16 : dbcsr_get_info,&
17 : dbcsr_init_p,&
18 : dbcsr_p_type,&
19 : dbcsr_type
20 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
21 : cp_dbcsr_plus_fm_fm_t,&
22 : cp_dbcsr_sm_fm_multiply,&
23 : dbcsr_allocate_matrix_set
24 : USE cp_fm_basic_linalg, ONLY: cp_fm_triangular_invert
25 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
26 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
27 : cp_fm_struct_p_type,&
28 : cp_fm_struct_release,&
29 : cp_fm_struct_type
30 : USE cp_fm_types, ONLY: cp_fm_create,&
31 : cp_fm_get_info,&
32 : cp_fm_release,&
33 : cp_fm_set_all,&
34 : cp_fm_to_fm,&
35 : cp_fm_to_fm_submat,&
36 : cp_fm_type
37 : USE cp_log_handling, ONLY: cp_get_default_logger,&
38 : cp_logger_get_default_io_unit,&
39 : cp_logger_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_ks_matrix, 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 1054 : 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 1054 : REAL(kind=dp), DIMENSION(:), POINTER :: evals_virt_spin
112 : TYPE(cell_type), POINTER :: cell
113 1054 : 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 1054 : TARGET :: mos_virt
117 : TYPE(cp_fm_type), POINTER :: mos_virt_spin
118 1054 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
119 : TYPE(dft_control_type), POINTER :: dft_control
120 1054 : 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 1054 : CALL timeset(routineN, handle)
126 :
127 : CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, &
128 1054 : matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, scf_env=scf_env)
129 1054 : tddfpt_control => dft_control%tddfpt2_control
130 :
131 1054 : CPASSERT(.NOT. ASSOCIATED(gs_mos))
132 : ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals
133 1054 : nspins = dft_control%nspins
134 4340 : ALLOCATE (gs_mos(nspins))
135 :
136 : ! check if virtuals should be constructed for NAMD interface with NEWTONX
137 1054 : print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
138 1054 : 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 1054 : NULLIFY (evals_virt, evals_virt_spin, mos_virt_spin)
144 1054 : IF (ASSOCIATED(scf_env)) THEN
145 1054 : 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 2 : nmo_virt = HUGE(0)
151 4 : DO ispin = 1, nspins
152 2 : CALL get_mo_set(mos(ispin), nmo=nmo_avail, homo=nmo_occ)
153 4 : nmo_virt = MIN(nmo_virt, nmo_avail - nmo_occ)
154 : END DO
155 : ! number of unoccupied orbitals to compute
156 2 : nmo_virt = tddfpt_control%nlumo - nmo_virt
157 2 : IF (.NOT. print_virtuals_newtonx) THEN
158 0 : IF (nmo_virt > 0) THEN
159 0 : ALLOCATE (evals_virt(nspins), mos_virt(nspins))
160 : ! the number of actually computed unoccupied orbitals will be stored as 'nmo_avail'
161 0 : 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 2232 : DO ispin = 1, nspins
168 1178 : IF (ASSOCIATED(evals_virt)) THEN
169 0 : evals_virt_spin => evals_virt(ispin)%array
170 : ELSE
171 1178 : NULLIFY (evals_virt_spin)
172 : END IF
173 1178 : IF (ALLOCATED(mos_virt)) THEN
174 0 : mos_virt_spin => mos_virt(ispin)
175 : ELSE
176 1178 : 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 2232 : qs_env=qs_env)
184 : END DO
185 :
186 1054 : moc = 0
187 1054 : mvt = 0
188 2232 : DO ispin = 1, nspins
189 1178 : CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, nrow_global=moc(1, ispin), ncol_global=moc(2, ispin))
190 2232 : CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, nrow_global=mvt(1, ispin), ncol_global=mvt(2, ispin))
191 : END DO
192 1054 : IF (iounit > 0) THEN
193 527 : WRITE (iounit, "(T2,A,T36,A)") "TDDFPT| Molecular Orbitals:", &
194 1054 : " Spin AOs Occ Virt Total"
195 1116 : DO ispin = 1, nspins
196 589 : WRITE (iounit, "(T2,A,T37,I4,4I10)") "TDDFPT| ", ispin, moc(1, ispin), moc(2, ispin), &
197 1705 : mvt(2, ispin), moc(2, ispin) + mvt(2, ispin)
198 : END DO
199 : END IF
200 :
201 1054 : IF (ASSOCIATED(evals_virt)) THEN
202 0 : DO ispin = 1, SIZE(evals_virt)
203 0 : IF (ASSOCIATED(evals_virt(ispin)%array)) DEALLOCATE (evals_virt(ispin)%array)
204 : END DO
205 0 : DEALLOCATE (evals_virt)
206 : END IF
207 :
208 1054 : CALL cp_fm_release(mos_virt)
209 :
210 1054 : CALL timestop(handle)
211 :
212 3162 : 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 1178 : 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 1178 : INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, &
253 1178 : sum_sign_array
254 1178 : 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 1178 : POINTER :: my_block
259 1178 : REAL(kind=dp), DIMENSION(:), POINTER :: mo_evals_extended, mo_occ_extended, &
260 1178 : 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 1178 : CALL timeset(routineN, handle)
272 :
273 1178 : NULLIFY (logger)
274 1178 : logger => cp_get_default_logger()
275 1178 : iounit = cp_logger_get_default_io_unit(logger)
276 :
277 1178 : 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 1178 : nelectron=nelectrons, occupation_numbers=mo_occ_scf)
281 :
282 1178 : print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
283 1178 : CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases)
284 :
285 1178 : nmo_virt = nao - nmo_occ
286 1178 : IF (nlumo >= 0) &
287 0 : nmo_virt = MIN(nmo_virt, nlumo)
288 :
289 1178 : IF (nmo_virt <= 0) &
290 : CALL cp_abort(__LOCATION__, &
291 0 : 'At least one unoccupied molecular orbital is required to calculate excited states.')
292 :
293 1178 : do_eigen = .FALSE.
294 : ! diagonalise the Kohn-Sham matrix one more time if the number of available unoccupied states are too small
295 1178 : IF (ASSOCIATED(evals_virt)) THEN
296 0 : CPASSERT(ASSOCIATED(mos_virt))
297 0 : IF (nmo_virt > nmo_scf - nmo_occ + SIZE(evals_virt)) do_eigen = .TRUE.
298 : ELSE
299 1178 : IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .TRUE.
300 : END IF
301 :
302 : ! ++ allocate storage space for gs_mos
303 1178 : NULLIFY (ao_mo_occ_fm_struct, ao_mo_virt_fm_struct)
304 1178 : CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
305 1178 : CALL cp_fm_struct_create(ao_mo_virt_fm_struct, nrow_global=nao, ncol_global=nmo_virt, context=blacs_env)
306 :
307 1178 : NULLIFY (gs_mos%mos_occ, gs_mos%mos_virt, gs_mos%evals_occ_matrix)
308 1178 : ALLOCATE (gs_mos%mos_occ, gs_mos%mos_virt)
309 1178 : CALL cp_fm_create(gs_mos%mos_occ, ao_mo_occ_fm_struct)
310 1178 : CALL cp_fm_create(gs_mos%mos_virt, ao_mo_virt_fm_struct)
311 :
312 3534 : ALLOCATE (gs_mos%evals_occ(nmo_occ))
313 3534 : ALLOCATE (gs_mos%evals_virt(nmo_virt))
314 2356 : ALLOCATE (gs_mos%phases_occ(nmo_occ))
315 2356 : ALLOCATE (gs_mos%phases_virt(nmo_virt))
316 :
317 : ! ++ nullify pointers
318 1178 : NULLIFY (ao_ao_fm_struct, wfn_fm_struct)
319 1178 : NULLIFY (mos_extended, mo_coeff_extended, mo_evals_extended, mo_occ_extended)
320 1178 : CALL cp_fm_struct_create(ao_ao_fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
321 :
322 1178 : IF (do_eigen) THEN
323 : ! ++ set of molecular orbitals
324 1178 : CALL cp_fm_struct_create(wfn_fm_struct, nrow_global=nao, ncol_global=nmo_occ + nmo_virt, context=blacs_env)
325 1178 : ALLOCATE (mos_extended)
326 : CALL allocate_mo_set(mos_extended, nao, nmo_occ + nmo_virt, nelectrons, &
327 1178 : REAL(nelectrons, dp), maxocc, flexible_electron_count=0.0_dp)
328 1178 : CALL init_mo_set(mos_extended, fm_struct=wfn_fm_struct, name="mos-extended")
329 1178 : CALL cp_fm_struct_release(wfn_fm_struct)
330 : CALL get_mo_set(mos_extended, mo_coeff=mo_coeff_extended, &
331 1178 : 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 7448 : DO imo = 1, nmo_scf
339 7448 : mo_occ_extended(imo) = mo_occ_scf(imo)
340 : END DO
341 20430 : mo_occ_extended(nmo_scf + 1:) = 0.0_dp
342 :
343 : ! ++ allocate temporary matrices
344 1178 : CALL cp_fm_create(matrix_ks_fm, ao_ao_fm_struct)
345 1178 : CALL cp_fm_create(ortho_fm, ao_ao_fm_struct)
346 1178 : CALL cp_fm_create(work_fm, ao_ao_fm_struct)
347 1178 : CALL cp_fm_struct_release(ao_ao_fm_struct)
348 :
349 : ! some stuff from the subroutine general_eigenproblem()
350 1178 : CALL copy_dbcsr_to_fm(matrix_s, ortho_fm)
351 1178 : CALL copy_dbcsr_to_fm(matrix_ks, matrix_ks_fm)
352 :
353 1178 : IF (cholesky_method == cholesky_dbcsr) THEN
354 0 : CPABORT('CHOLESKY DBCSR_INVERSE is not implemented in TDDFT.')
355 1178 : ELSE IF (cholesky_method == cholesky_off) THEN
356 0 : CPABORT('CHOLESKY OFF is not implemented in TDDFT.')
357 : ELSE
358 1178 : CALL cp_fm_cholesky_decompose(ortho_fm)
359 1178 : IF (cholesky_method == cholesky_inverse) THEN
360 0 : 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 1178 : 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 1178 : do_level_shift=.FALSE., level_shift=0.0_dp, use_jacobi=.FALSE.)
369 : END IF
370 :
371 : ! -- clean up needless matrices
372 1178 : CALL cp_fm_release(work_fm)
373 1178 : CALL cp_fm_release(ortho_fm)
374 1178 : CALL cp_fm_release(matrix_ks_fm)
375 : ELSE
376 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff_extended, &
377 0 : 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 1178 : CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
384 1178 : CALL cp_fm_struct_release(ao_mo_occ_fm_struct)
385 :
386 1178 : 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 1178 : row_indices=row_indices, col_indices=col_indices, local_data=my_block)
389 :
390 5890 : ALLOCATE (minrow_neg_array(nmo_occ), minrow_pos_array(nmo_occ), sum_sign_array(nmo_occ))
391 7448 : minrow_neg_array(:) = nao
392 7448 : minrow_pos_array(:) = nao
393 7448 : sum_sign_array(:) = 0
394 7448 : DO icol_local = 1, ncol_local
395 6270 : icol_global = col_indices(icol_local)
396 :
397 222770 : DO irow_local = 1, nrow_local
398 215322 : element = my_block(irow_local, icol_local)
399 :
400 215322 : sign_int = 0
401 215322 : IF (element >= eps_dp) THEN
402 : sign_int = 1
403 108838 : ELSE IF (element <= -eps_dp) THEN
404 108319 : sign_int = -1
405 : END IF
406 :
407 215322 : sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
408 :
409 215322 : irow_global = row_indices(irow_local)
410 221592 : IF (sign_int > 0) THEN
411 106484 : IF (minrow_pos_array(icol_global) > irow_global) &
412 6147 : minrow_pos_array(icol_global) = irow_global
413 108838 : ELSE IF (sign_int < 0) THEN
414 108319 : IF (minrow_neg_array(icol_global) > irow_global) &
415 6034 : minrow_neg_array(icol_global) = irow_global
416 : END IF
417 : END DO
418 : END DO
419 :
420 1178 : CALL para_env%sum(sum_sign_array)
421 1178 : CALL para_env%min(minrow_neg_array)
422 1178 : CALL para_env%min(minrow_pos_array)
423 :
424 7448 : DO icol_local = 1, nmo_occ
425 7448 : IF (sum_sign_array(icol_local) > 0) THEN
426 : ! most of the expansion coefficients are positive => MO's phase = +1
427 2766 : gs_mos%phases_occ(icol_local) = 1.0_dp
428 3504 : ELSE IF (sum_sign_array(icol_local) < 0) THEN
429 : ! most of the expansion coefficients are negative => MO's phase = -1
430 3258 : gs_mos%phases_occ(icol_local) = -1.0_dp
431 : ELSE
432 : ! equal number of positive and negative expansion coefficients
433 246 : 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 128 : gs_mos%phases_occ(icol_local) = 1.0_dp
437 : ELSE
438 : ! MO's phase = -1
439 118 : gs_mos%phases_occ(icol_local) = -1.0_dp
440 : END IF
441 : END IF
442 : END DO
443 :
444 1178 : 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 1178 : CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_occ, ncol=nmo_occ, source_start=1, target_start=1)
448 7448 : gs_mos%evals_occ(1:nmo_occ) = mo_evals_extended(1:nmo_occ)
449 :
450 1178 : 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 0 : 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 0 : source_start=1, target_start=nmo_scf - nmo_occ + 1)
455 0 : gs_mos%evals_virt(1:nmo_scf - nmo_occ) = evals_virt(nmo_occ + 1:nmo_occ + nmo_scf)
456 0 : gs_mos%evals_virt(nmo_scf - nmo_occ + 1:nmo_virt) = evals_virt(1:nmo_virt - (nmo_scf - nmo_occ))
457 : ELSE
458 1178 : CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_virt, source_start=nmo_occ + 1, target_start=1)
459 20430 : gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt)
460 : END IF
461 :
462 1178 : 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 0 : CALL cp_fm_create(work_fm_virt, ao_mo_virt_fm_struct)
467 :
468 0 : 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 0 : row_indices=row_indices, col_indices=col_indices, local_data=my_block)
471 :
472 0 : ALLOCATE (minrow_neg_array(nmo_virt), minrow_pos_array(nmo_virt), sum_sign_array(nmo_virt))
473 0 : minrow_neg_array(:) = nao
474 0 : minrow_pos_array(:) = nao
475 0 : sum_sign_array(:) = 0
476 0 : DO icol_local = 1, ncol_local
477 0 : icol_global = col_indices(icol_local)
478 :
479 0 : DO irow_local = 1, nrow_local
480 0 : element = my_block(irow_local, icol_local)
481 :
482 0 : sign_int = 0
483 0 : IF (element >= eps_dp) THEN
484 : sign_int = 1
485 0 : ELSE IF (element <= -eps_dp) THEN
486 0 : sign_int = -1
487 : END IF
488 :
489 0 : sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
490 :
491 0 : irow_global = row_indices(irow_local)
492 0 : IF (sign_int > 0) THEN
493 0 : IF (minrow_pos_array(icol_global) > irow_global) &
494 0 : minrow_pos_array(icol_global) = irow_global
495 0 : ELSE IF (sign_int < 0) THEN
496 0 : IF (minrow_neg_array(icol_global) > irow_global) &
497 0 : minrow_neg_array(icol_global) = irow_global
498 : END IF
499 : END DO
500 : END DO
501 :
502 0 : CALL para_env%sum(sum_sign_array)
503 0 : CALL para_env%min(minrow_neg_array)
504 0 : CALL para_env%min(minrow_pos_array)
505 0 : DO icol_local = 1, nmo_virt
506 0 : IF (sum_sign_array(icol_local) > 0) THEN
507 : ! most of the expansion coefficients are positive => MO's phase = +1
508 0 : gs_mos%phases_virt(icol_local) = 1.0_dp
509 0 : ELSE IF (sum_sign_array(icol_local) < 0) THEN
510 : ! most of the expansion coefficients are negative => MO's phase = -1
511 0 : gs_mos%phases_virt(icol_local) = -1.0_dp
512 : ELSE
513 : ! equal number of positive and negative expansion coefficients
514 0 : 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 0 : gs_mos%phases_virt(icol_local) = 1.0_dp
518 : ELSE
519 : ! MO's phase = -1
520 0 : gs_mos%phases_virt(icol_local) = -1.0_dp
521 : END IF
522 : END IF
523 : END DO
524 0 : DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
525 0 : CALL cp_fm_release(work_fm_virt)
526 : END IF !print_phases
527 1178 : CALL cp_fm_struct_release(ao_mo_virt_fm_struct) ! here after print_phases
528 :
529 1178 : CALL cp_fm_release(work_fm)
530 :
531 1178 : IF (do_eigen) THEN
532 1178 : CALL deallocate_mo_set(mos_extended)
533 1178 : DEALLOCATE (mos_extended)
534 : END IF
535 :
536 1178 : CALL timestop(handle)
537 :
538 7068 : 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 1178 : 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 1178 : CALL timeset(routineN, handle)
554 :
555 1178 : IF (ALLOCATED(gs_mos%phases_occ)) &
556 1178 : DEALLOCATE (gs_mos%phases_occ)
557 :
558 1178 : IF (ALLOCATED(gs_mos%evals_virt)) &
559 1178 : DEALLOCATE (gs_mos%evals_virt)
560 :
561 1178 : IF (ALLOCATED(gs_mos%evals_occ)) &
562 1178 : DEALLOCATE (gs_mos%evals_occ)
563 :
564 1178 : IF (ALLOCATED(gs_mos%phases_virt)) &
565 1178 : DEALLOCATE (gs_mos%phases_virt)
566 :
567 1178 : IF (ASSOCIATED(gs_mos%evals_occ_matrix)) THEN
568 36 : CALL cp_fm_release(gs_mos%evals_occ_matrix)
569 36 : DEALLOCATE (gs_mos%evals_occ_matrix)
570 : END IF
571 :
572 1178 : IF (ASSOCIATED(gs_mos%mos_virt)) THEN
573 1178 : CALL cp_fm_release(gs_mos%mos_virt)
574 1178 : DEALLOCATE (gs_mos%mos_virt)
575 : END IF
576 :
577 1178 : IF (ASSOCIATED(gs_mos%mos_occ)) THEN
578 1178 : CALL cp_fm_release(gs_mos%mos_occ)
579 1178 : DEALLOCATE (gs_mos%mos_occ)
580 : END IF
581 :
582 1178 : CALL timestop(handle)
583 1178 : 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 1054 : 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 1054 : 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 1054 : CALL timeset(routineN, handle)
614 :
615 1054 : NULLIFY (logger)
616 1054 : logger => cp_get_default_logger()
617 1054 : iounit = cp_logger_get_default_io_unit(logger)
618 :
619 1054 : CALL get_qs_env(qs_env, blacs_env=blacs_env, dft_control=dft_control, matrix_ks=matrix_ks)
620 1054 : tddfpt_control => dft_control%tddfpt2_control
621 :
622 : ! obtain corrected KS-matrix
623 : ! We should 'save' the energy values?
624 1054 : nspins = SIZE(gs_mos)
625 1054 : NULLIFY (matrix_ks_oep)
626 1054 : IF (tddfpt_control%oe_corr /= oe_none) THEN
627 30 : IF (iounit > 0) THEN
628 15 : WRITE (iounit, "(1X,A)") "", &
629 15 : "-------------------------------------------------------------------------------", &
630 15 : "- Orbital Eigenvalue Correction Started -", &
631 30 : "-------------------------------------------------------------------------------"
632 : END IF
633 :
634 : CALL cp_warn(__LOCATION__, &
635 : "Orbital energy correction potential is an experimental feature. "// &
636 30 : "Use it with extreme care")
637 :
638 30 : hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
639 30 : CALL section_vals_get(hfx_section, explicit=do_hfx)
640 30 : IF (do_hfx) THEN
641 : CALL cp_abort(__LOCATION__, &
642 : "Implementation of orbital energy correction XC-potentials is "// &
643 0 : "currently incompatible with exact-exchange functionals")
644 : END IF
645 :
646 30 : CALL dbcsr_allocate_matrix_set(matrix_ks_oep, nspins)
647 66 : DO ispin = 1, nspins
648 36 : CALL dbcsr_init_p(matrix_ks_oep(ispin)%matrix)
649 66 : CALL dbcsr_copy(matrix_ks_oep(ispin)%matrix, matrix_ks(ispin)%matrix)
650 : END DO
651 :
652 : ! KS-matrix without XC-terms
653 30 : xc_fun_original => section_vals_get_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL")
654 30 : CALL section_vals_retain(xc_fun_original)
655 30 : NULLIFY (xc_fun_empty)
656 30 : CALL section_vals_create(xc_fun_empty, xc_fun_original%section)
657 30 : CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_empty)
658 30 : CALL section_vals_release(xc_fun_empty)
659 :
660 30 : IF (dft_control%qs_control%semi_empirical) THEN
661 0 : CPABORT("TDDFPT with SE not possible")
662 30 : ELSEIF (dft_control%qs_control%dftb) THEN
663 0 : CPABORT("TDDFPT with DFTB not possible")
664 30 : ELSEIF (dft_control%qs_control%xtb) THEN
665 : CALL build_xtb_ks_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
666 16 : 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 14 : ext_ks_matrix=matrix_ks_oep)
670 : END IF
671 :
672 : IF (tddfpt_control%oe_corr == oe_saop .OR. &
673 30 : tddfpt_control%oe_corr == oe_lb .OR. &
674 : tddfpt_control%oe_corr == oe_gllb) THEN
675 14 : IF (iounit > 0) THEN
676 7 : WRITE (iounit, "(T2,A)") " Orbital energy correction of SAOP type "
677 : END IF
678 14 : CALL add_saop_pot(matrix_ks_oep, qs_env, tddfpt_control%oe_corr)
679 16 : ELSE IF (tddfpt_control%oe_corr == oe_shift) THEN
680 16 : IF (iounit > 0) THEN
681 : WRITE (iounit, "(T2,A,T71,F10.3)") &
682 8 : " Virtual Orbital Eigenvalue Shift [eV] ", tddfpt_control%ev_shift*evolt
683 : WRITE (iounit, "(T2,A,T71,F10.3)") &
684 8 : " 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 16 : tddfpt_control%ev_shift, tddfpt_control%eos_shift)
688 : ELSE
689 : CALL cp_abort(__LOCATION__, &
690 0 : "Unimplemented orbital energy correction potential")
691 : END IF
692 30 : CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_original)
693 30 : CALL section_vals_release(xc_fun_original)
694 :
695 : ! compute 'evals_occ_matrix'
696 30 : CALL dbcsr_get_info(matrix_ks(1)%matrix, nfullrows_total=nao)
697 30 : NULLIFY (mo_occ_mo_occ_fm_struct)
698 66 : DO ispin = 1, nspins
699 36 : 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 36 : context=blacs_env)
702 36 : ALLOCATE (gs_mos(ispin)%evals_occ_matrix)
703 36 : CALL cp_fm_create(gs_mos(ispin)%evals_occ_matrix, mo_occ_mo_occ_fm_struct)
704 36 : 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 36 : context=blacs_env)
708 36 : CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
709 36 : 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 36 : 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 36 : 0.0_dp, gs_mos(ispin)%evals_occ_matrix)
714 102 : CALL cp_fm_release(work_fm)
715 : END DO
716 30 : IF (iounit > 0) THEN
717 : WRITE (iounit, "(1X,A)") &
718 15 : "-------------------------------------------------------------------------------"
719 : END IF
720 :
721 : END IF
722 :
723 1054 : CALL timestop(handle)
724 :
725 1054 : 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 2653 : 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 2653 : nstates_total = 0
742 2653 : nspins = SIZE(gs_mos)
743 :
744 5614 : DO ispin = 1, nspins
745 : nstates_total = nstates_total + &
746 : SIZE(gs_mos(ispin)%evals_occ, kind=int_8)* &
747 5614 : SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
748 : END DO
749 2653 : 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 16 : 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 16 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
780 : TYPE(dbcsr_type), POINTER :: smat
781 16 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
782 :
783 16 : CALL timeset(routineN, handle)
784 :
785 16 : n_spins = SIZE(gs_mos)
786 16 : CPASSERT(n_spins == SIZE(matrix_ks))
787 :
788 16 : IF (eos_shift /= 0.0_dp .AND. n_spins > 1) THEN
789 0 : CPABORT("eos_shift not implemented")
790 0 : CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
791 0 : smat => matrix_s(1)%matrix
792 0 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
793 0 : CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
794 0 : nl = MIN(na, nb)
795 0 : nu = MAX(na, nb)
796 : ! open shell orbital shift
797 0 : DO ispin = 1, n_spins
798 0 : coeff => gs_mos(ispin)%mos_occ
799 0 : CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
800 0 : IF (nhomo == nu) THEN
801 : ! downshift with -eos_shift using occupied orbitals
802 0 : nos = nu - nl
803 0 : CALL cp_fm_create(cmos, fmstruct)
804 0 : CALL cp_fm_get_info(coeff, nrow_global=nrow)
805 0 : CALL cp_fm_to_fm_submat(coeff, cmos, nrow, nos, 1, nl + 1, 1, 1)
806 0 : CALL cp_fm_create(cvec, fmstruct)
807 0 : 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 0 : alpha=-eos_shift, keep_sparsity=.TRUE.)
810 0 : CALL cp_fm_release(cmos)
811 0 : CALL cp_fm_release(cvec)
812 : ELSE
813 : ! upshift with eos_shift using virtual orbitals
814 0 : coeff => gs_mos(ispin)%mos_virt
815 0 : CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
816 0 : nos = nu - nhomo
817 0 : CPASSERT(nvirt >= nos)
818 0 : CALL cp_fm_create(cvec, fmstruct)
819 0 : 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 0 : alpha=eos_shift, keep_sparsity=.TRUE.)
822 0 : CALL cp_fm_release(cvec)
823 : END IF
824 : END DO
825 : ! virtual shift
826 0 : IF (ev_shift /= 0.0_dp) THEN
827 0 : DO ispin = 1, n_spins
828 : CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
829 0 : alpha_scalar=1.0_dp, beta_scalar=ev_shift)
830 0 : coeff => gs_mos(ispin)%mos_occ
831 0 : CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
832 0 : CALL cp_fm_create(cvec, fmstruct)
833 0 : 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 0 : alpha=-ev_shift, keep_sparsity=.TRUE.)
836 0 : CALL cp_fm_release(cvec)
837 0 : IF (nhomo < nu) THEN
838 0 : nos = nu - nhomo
839 0 : coeff => gs_mos(ispin)%mos_virt
840 0 : CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
841 0 : CPASSERT(nvirt >= nos)
842 0 : CALL cp_fm_create(cvec, fmstruct)
843 0 : 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 0 : alpha=-ev_shift, keep_sparsity=.TRUE.)
846 0 : CALL cp_fm_release(cvec)
847 : END IF
848 : END DO
849 : END IF
850 : ELSE
851 : ! virtual shift
852 16 : IF (ev_shift /= 0.0_dp) THEN
853 16 : CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
854 16 : smat => matrix_s(1)%matrix
855 32 : DO ispin = 1, n_spins
856 : CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
857 16 : alpha_scalar=1.0_dp, beta_scalar=ev_shift)
858 16 : coeff => gs_mos(ispin)%mos_occ
859 16 : CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
860 16 : CALL cp_fm_create(cvec, fmstruct)
861 16 : 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 16 : alpha=-ev_shift, keep_sparsity=.TRUE.)
864 48 : CALL cp_fm_release(cvec)
865 : END DO
866 : END IF
867 : END IF
868 : ! set eigenvalues
869 16 : IF (eos_shift == 0.0_dp .OR. n_spins == 1) THEN
870 32 : DO ispin = 1, n_spins
871 32 : IF (ALLOCATED(gs_mos(ispin)%evals_virt)) THEN
872 1332 : gs_mos(ispin)%evals_virt(:) = gs_mos(ispin)%evals_virt(:) + ev_shift
873 : END IF
874 : END DO
875 : ELSE
876 0 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
877 0 : CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
878 0 : nl = MIN(na, nb)
879 0 : nu = MAX(na, nb)
880 0 : nos = nu - nl
881 0 : IF (na == nu) THEN
882 0 : IF (ALLOCATED(gs_mos(1)%evals_occ)) THEN
883 0 : gs_mos(1)%evals_occ(nl + 1:nu) = gs_mos(1)%evals_occ(nl + 1:nu) - eos_shift
884 : END IF
885 0 : IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
886 0 : gs_mos(1)%evals_virt(:) = gs_mos(1)%evals_virt(:) + ev_shift
887 : END IF
888 0 : IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
889 0 : gs_mos(2)%evals_virt(1:nos) = gs_mos(2)%evals_virt(1:nos) + eos_shift
890 0 : gs_mos(2)%evals_virt(nos + 1:) = gs_mos(2)%evals_virt(nos + 1:) + ev_shift
891 : END IF
892 : ELSE
893 0 : IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
894 0 : gs_mos(1)%evals_virt(1:nos) = gs_mos(1)%evals_virt(1:nos) + eos_shift
895 0 : gs_mos(1)%evals_virt(nos + 1:) = gs_mos(1)%evals_virt(nos + 1:) + ev_shift
896 : END IF
897 0 : IF (ALLOCATED(gs_mos(2)%evals_occ)) THEN
898 0 : gs_mos(2)%evals_occ(nl + 1:nu) = gs_mos(2)%evals_occ(nl + 1:nu) - eos_shift
899 : END IF
900 0 : IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
901 0 : gs_mos(2)%evals_virt(:) = gs_mos(2)%evals_virt(:) + ev_shift
902 : END IF
903 : END IF
904 : END IF
905 :
906 16 : CALL timestop(handle)
907 :
908 16 : 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 1062 : 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 1062 : 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 1062 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: e_virt_minus_occ
949 3186 : TYPE(cp_fm_struct_p_type), DIMENSION(maxspins) :: fm_struct_evects
950 :
951 1062 : CALL timeset(routineN, handle)
952 :
953 1062 : nspins = SIZE(evects, 1)
954 1062 : 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 2248 : DO ispin = 1, nspins
962 : ! number of occupied orbitals for each spin component
963 1186 : 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 1186 : nmo_occ_selected(ispin) = MIN(nmo_occ_avail(ispin), nstates)
967 1186 : nmo_virt_selected(ispin) = MIN(SIZE(gs_mos(ispin)%evals_virt), nstates)
968 :
969 2248 : 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 2248 : nstates_selected = DOT_PRODUCT(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins))
975 :
976 3186 : ALLOCATE (inds(nstates_selected))
977 3186 : ALLOCATE (e_virt_minus_occ(nstates_selected))
978 :
979 1062 : istate = 0
980 2248 : DO ispin = 1, nspins
981 5066 : 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 2818 : e_occ = gs_mos(ispin)%evals_occ(nmo_occ_avail(ispin) - imo_occ + 1)
984 :
985 13930 : DO imo_virt = 1, nmo_virt_selected(ispin)
986 9926 : istate = istate + 1
987 12744 : 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 1062 : CALL sort(e_virt_minus_occ, nstates_selected, inds)
997 :
998 1062 : IF (nspins == 1) THEN
999 938 : ispin = 1
1000 938 : spin_label = ' '
1001 : END IF
1002 :
1003 1062 : nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1)
1004 1062 : IF (log_unit > 0) THEN
1005 531 : WRITE (log_unit, "(1X,A)") "", &
1006 531 : "-------------------------------------------------------------------------------", &
1007 531 : "- TDDFPT Initial Guess -", &
1008 1062 : "-------------------------------------------------------------------------------"
1009 531 : WRITE (log_unit, '(T11,A)') "State Occupied -> Virtual Excitation"
1010 531 : WRITE (log_unit, '(T11,A)') "number orbital orbital energy (eV)"
1011 531 : WRITE (log_unit, '(1X,79("-"))')
1012 : END IF
1013 :
1014 3936 : DO istate = 1, nstates
1015 3936 : IF (ASSOCIATED(evects(1, istate)%matrix_struct)) THEN
1016 6 : IF (log_unit > 0) &
1017 : WRITE (log_unit, '(T7,I8,T28,A19,T60,F14.5)') &
1018 3 : istate, "*** restarted ***", evals(istate)*evolt
1019 : ELSE
1020 2868 : ind = inds(istate) - 1
1021 2868 : IF (nspins > 1) THEN
1022 372 : IF (ind < nstates_occ_virt_alpha) THEN
1023 156 : ispin = 1
1024 156 : spin_label = '(alp)'
1025 : ELSE
1026 216 : ispin = 2
1027 216 : ind = ind - nstates_occ_virt_alpha
1028 216 : spin_label = '(bet)'
1029 : END IF
1030 : END IF
1031 :
1032 2868 : imo_occ = nmo_occ_avail(ispin) - ind/nmo_virt_selected(ispin)
1033 2868 : imo_virt = MOD(ind, nmo_virt_selected(ispin)) + 1
1034 2868 : evals(istate) = e_virt_minus_occ(istate)
1035 :
1036 2868 : IF (log_unit > 0) &
1037 : WRITE (log_unit, '(T7,I8,T24,I8,T37,A5,T45,I8,T54,A5,T60,F14.5)') &
1038 1434 : istate, imo_occ, spin_label, nmo_occ_avail(ispin) + imo_virt, spin_label, e_virt_minus_occ(istate)*evolt
1039 :
1040 6108 : DO jspin = 1, nspins
1041 : ! .NOT. ASSOCIATED(evects(jspin, istate)%matrix_struct))
1042 3240 : CALL cp_fm_create(evects(jspin, istate), fm_struct_evects(jspin)%struct)
1043 3240 : CALL cp_fm_set_all(evects(jspin, istate), 0.0_dp)
1044 :
1045 3240 : IF (jspin == ispin) &
1046 : CALL cp_fm_to_fm(gs_mos(ispin)%mos_virt, evects(ispin, istate), &
1047 5736 : ncol=1, source_start=imo_virt, target_start=imo_occ)
1048 : END DO
1049 : END IF
1050 : END DO
1051 :
1052 1062 : IF (log_unit > 0) THEN
1053 531 : 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 531 : "-------------------------------------------------------------------------------"
1056 : END IF
1057 :
1058 1062 : DEALLOCATE (e_virt_minus_occ)
1059 1062 : DEALLOCATE (inds)
1060 :
1061 1062 : CALL timestop(handle)
1062 :
1063 1062 : END SUBROUTINE tddfpt_guess_vectors
1064 :
1065 : END MODULE qs_tddfpt2_utils
|
