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 : ! **************************************************************************************************
9 : !> \brief Initialize the XAS orbitals for specific core excitations
10 : !> Either the GS orbitals are used as initial guess, or the
11 : !> xas mos are read from a previous calculation.
12 : !> In the latter case, the core-hole potetial should be the same.
13 : !> \note
14 : !> The restart with the same core-hole potential should be checked
15 : !> and a wrong restart should stop the program
16 : !> \par History
17 : !> created 09.2006
18 : !> \author MI (09.2006)
19 : ! **************************************************************************************************
20 : MODULE xas_restart
21 :
22 : USE cp_control_types, ONLY: dft_control_type
23 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
24 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
25 : USE cp_files, ONLY: close_file,&
26 : open_file
27 : USE cp_fm_types, ONLY: cp_fm_create,&
28 : cp_fm_get_info,&
29 : cp_fm_get_submatrix,&
30 : cp_fm_release,&
31 : cp_fm_set_all,&
32 : cp_fm_set_submatrix,&
33 : cp_fm_type,&
34 : cp_fm_write_unformatted
35 : USE cp_log_handling, ONLY: cp_get_default_logger,&
36 : cp_logger_type,&
37 : cp_to_string
38 : USE cp_output_handling, ONLY: cp_p_file,&
39 : cp_print_key_finished_output,&
40 : cp_print_key_generate_filename,&
41 : cp_print_key_should_output,&
42 : cp_print_key_unit_nr
43 : USE input_section_types, ONLY: section_vals_type
44 : USE kinds, ONLY: default_path_length,&
45 : default_string_length,&
46 : dp
47 : USE message_passing, ONLY: mp_para_env_type
48 : USE parallel_gemm_api, ONLY: parallel_gemm
49 : USE particle_types, ONLY: particle_type
50 : USE qs_density_matrices, ONLY: calculate_density_matrix
51 : USE qs_environment_types, ONLY: get_qs_env,&
52 : qs_environment_type
53 : USE qs_kind_types, ONLY: qs_kind_type
54 : USE qs_ks_types, ONLY: qs_ks_did_change
55 : USE qs_mixing_utils, ONLY: mixing_init
56 : USE qs_mo_io, ONLY: wfn_restart_file_name,&
57 : write_mo_set_low
58 : USE qs_mo_occupation, ONLY: set_mo_occupation
59 : USE qs_mo_types, ONLY: get_mo_set,&
60 : mo_set_type,&
61 : set_mo_set
62 : USE qs_rho_atom_types, ONLY: rho_atom_type
63 : USE qs_rho_methods, ONLY: qs_rho_update_rho
64 : USE qs_rho_types, ONLY: qs_rho_get,&
65 : qs_rho_type
66 : USE qs_scf_types, ONLY: qs_scf_env_type
67 : USE scf_control_types, ONLY: scf_control_type
68 : USE string_utilities, ONLY: xstring
69 : USE xas_env_types, ONLY: get_xas_env,&
70 : set_xas_env,&
71 : xas_environment_type
72 : #include "./base/base_uses.f90"
73 :
74 : IMPLICIT NONE
75 : PRIVATE
76 :
77 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_restart'
78 :
79 : ! *** Public subroutines ***
80 :
81 : PUBLIC :: xas_read_restart, xas_write_restart, xas_initialize_rho, find_excited_core_orbital
82 :
83 : CONTAINS
84 :
85 : ! **************************************************************************************************
86 : !> \brief Set up for reading the restart
87 : !> corresponding to the excitation of iatom
88 : !> If the corresponding restart file does not exist
89 : !> the GS orbitals are used as initial guess
90 : !> \param xas_env ...
91 : !> \param xas_section input section for XAS calculations
92 : !> qs_env:
93 : !> \param qs_env ...
94 : !> \param xas_method ...
95 : !> \param iatom index of the absorbing atom
96 : !> \param estate index of the core-hole orbital
97 : !> \param istate counter of excited states per atom
98 : !> error:
99 : !> \par History
100 : !> 09.2006 created [MI]
101 : !> \author MI
102 : ! **************************************************************************************************
103 12 : SUBROUTINE xas_read_restart(xas_env, xas_section, qs_env, xas_method, iatom, estate, istate)
104 :
105 : TYPE(xas_environment_type), POINTER :: xas_env
106 : TYPE(section_vals_type), POINTER :: xas_section
107 : TYPE(qs_environment_type), POINTER :: qs_env
108 : INTEGER, INTENT(IN) :: xas_method, iatom
109 : INTEGER, INTENT(OUT) :: estate
110 : INTEGER, INTENT(IN) :: istate
111 :
112 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_read_restart'
113 :
114 : CHARACTER(LEN=default_path_length) :: filename
115 : INTEGER :: handle, i, ia, ie, ispin, my_spin, nao, nao_read, nelectron, nexc_atoms, &
116 : nexc_atoms_read, nexc_search, nexc_search_read, nmo, nmo_read, output_unit, rst_unit, &
117 : xas_estate, xas_estate_read, xas_method_read
118 : LOGICAL :: file_exists
119 : REAL(dp) :: occ_estate, occ_estate_read, &
120 : xas_nelectron, xas_nelectron_read
121 12 : REAL(dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
122 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig_read, occ_read
123 12 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
124 : TYPE(cp_fm_type), POINTER :: mo_coeff
125 : TYPE(cp_logger_type), POINTER :: logger
126 12 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
127 12 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
128 : TYPE(mp_para_env_type), POINTER :: para_env
129 :
130 12 : CALL timeset(routineN, handle)
131 :
132 12 : file_exists = .FALSE.
133 12 : rst_unit = -1
134 :
135 12 : NULLIFY (eigenvalues, matrix_s, mos, occupation_numbers, vecbuffer)
136 12 : NULLIFY (logger)
137 12 : logger => cp_get_default_logger()
138 :
139 : output_unit = cp_print_key_unit_nr(logger, xas_section, &
140 12 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
141 :
142 12 : CALL get_qs_env(qs_env=qs_env, para_env=para_env)
143 :
144 12 : IF (para_env%is_source()) THEN
145 : CALL wfn_restart_file_name(filename, file_exists, xas_section, logger, &
146 6 : xas=.TRUE.)
147 :
148 6 : CALL xstring(filename, ia, ie)
149 : filename = filename(ia:ie)//'-at'//TRIM(ADJUSTL(cp_to_string(iatom)))// &
150 6 : '_st'//TRIM(ADJUSTL(cp_to_string(istate)))//'.rst'
151 :
152 6 : INQUIRE (FILE=filename, EXIST=file_exists)
153 : ! open file
154 6 : IF (file_exists) THEN
155 :
156 : CALL open_file(file_name=TRIM(filename), &
157 : file_action="READ", &
158 : file_form="UNFORMATTED", &
159 : file_position="REWIND", &
160 : file_status="OLD", &
161 6 : unit_number=rst_unit)
162 :
163 6 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,T20,A,I5,/)") &
164 6 : "Read restart file for atom ", iatom
165 :
166 : ELSE IF (.NOT. file_exists) THEN
167 0 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,T10,A,I5,A,/)") &
168 0 : "Restart file for atom ", iatom, &
169 0 : " not available. Initialization done with GS orbitals"
170 : END IF
171 : END IF
172 12 : CALL para_env%bcast(file_exists)
173 :
174 : CALL get_xas_env(xas_env=xas_env, occ_estate=occ_estate, xas_estate=xas_estate, &
175 : xas_nelectron=xas_nelectron, nexc_search=nexc_search, &
176 12 : nexc_atoms=nexc_atoms, spin_channel=my_spin)
177 :
178 12 : IF (file_exists) THEN
179 12 : CALL get_qs_env(qs_env=qs_env, mos=mos, matrix_s=matrix_s)
180 :
181 12 : IF (rst_unit > 0) THEN
182 6 : READ (rst_unit) xas_method_read
183 6 : READ (rst_unit) nexc_search_read, nexc_atoms_read, occ_estate_read, xas_nelectron_read
184 6 : READ (rst_unit) xas_estate_read
185 :
186 6 : IF (xas_method_read /= xas_method) &
187 0 : CPABORT("READ XAS RESTART: restart with different XAS method is not possible.")
188 6 : IF (nexc_atoms_read /= nexc_atoms) &
189 : CALL cp_abort(__LOCATION__, &
190 : "READ XAS RESTART: restart with different excited atoms "// &
191 0 : "is not possible. Start instead a new XAS run with the new set of atoms.")
192 : END IF
193 :
194 12 : CALL para_env%bcast(xas_estate_read)
195 12 : CALL set_xas_env(xas_env=xas_env, xas_estate=xas_estate_read)
196 12 : estate = xas_estate_read
197 :
198 12 : CALL get_mo_set(mo_set=mos(my_spin), nao=nao)
199 36 : ALLOCATE (vecbuffer(1, nao))
200 :
201 36 : DO ispin = 1, SIZE(mos)
202 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, eigenvalues=eigenvalues, &
203 24 : occupation_numbers=occupation_numbers, mo_coeff=mo_coeff, nelectron=nelectron)
204 232 : eigenvalues = 0.0_dp
205 232 : occupation_numbers = 0.0_dp
206 24 : CALL cp_fm_set_all(mo_coeff, 0.0_dp)
207 24 : IF (para_env%is_source()) THEN
208 12 : READ (rst_unit) nao_read, nmo_read
209 12 : IF (nao /= nao_read) &
210 0 : CPABORT("To change basis is not possible. ")
211 48 : ALLOCATE (eig_read(nmo_read), occ_read(nmo_read))
212 116 : eig_read = 0.0_dp
213 116 : occ_read = 0.0_dp
214 12 : nmo = MIN(nmo, nmo_read)
215 12 : READ (rst_unit) eig_read(1:nmo_read), occ_read(1:nmo_read)
216 116 : eigenvalues(1:nmo) = eig_read(1:nmo)
217 116 : occupation_numbers(1:nmo) = occ_read(1:nmo)
218 12 : IF (nmo_read > nmo) THEN
219 0 : IF (occupation_numbers(nmo) >= EPSILON(0.0_dp)) &
220 : CALL cp_warn(__LOCATION__, &
221 : "The number of occupied MOs on the restart unit is larger than "// &
222 0 : "the allocated MOs.")
223 :
224 : END IF
225 12 : DEALLOCATE (eig_read, occ_read)
226 : END IF
227 440 : CALL para_env%bcast(eigenvalues)
228 440 : CALL para_env%bcast(occupation_numbers)
229 :
230 232 : DO i = 1, nmo
231 208 : IF (para_env%is_source()) THEN
232 5928 : READ (rst_unit) vecbuffer
233 : ELSE
234 3016 : vecbuffer(1, :) = 0.0_dp
235 : END IF
236 23504 : CALL para_env%bcast(vecbuffer)
237 : CALL cp_fm_set_submatrix(mo_coeff, &
238 232 : vecbuffer, 1, i, nao, 1, transpose=.TRUE.)
239 : END DO
240 : ! Skip extra MOs if there any
241 60 : IF (para_env%is_source()) THEN
242 12 : DO i = nmo + 1, nmo_read
243 12 : READ (rst_unit) vecbuffer
244 : END DO
245 : END IF
246 :
247 : END DO ! ispin
248 :
249 24 : DEALLOCATE (vecbuffer)
250 :
251 : ! nspin = SIZE(mos,1)
252 : ! DO ispin = 1,nspin
253 : ! ! ortho so that one can restart for different positions (basis sets?)
254 : ! NULLIFY(mo_coeff)
255 : ! CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff,homo=homo)
256 : ! CALL make_basis_sm(mo_coeff,homo,matrix_s(1)%matrix)
257 : ! END DO
258 : END IF !file_exist
259 :
260 12 : IF (para_env%is_source()) THEN
261 6 : IF (file_exists) CALL close_file(unit_number=rst_unit)
262 : END IF
263 :
264 12 : CALL timestop(handle)
265 :
266 24 : END SUBROUTINE xas_read_restart
267 :
268 : ! **************************************************************************************************
269 : !> \brief ...
270 : !> \param xas_env ...
271 : !> \param xas_section ...
272 : !> \param qs_env ...
273 : !> \param xas_method ...
274 : !> \param iatom ...
275 : !> \param istate ...
276 : ! **************************************************************************************************
277 1408 : SUBROUTINE xas_write_restart(xas_env, xas_section, qs_env, xas_method, iatom, istate)
278 :
279 : TYPE(xas_environment_type), POINTER :: xas_env
280 : TYPE(section_vals_type), POINTER :: xas_section
281 : TYPE(qs_environment_type), POINTER :: qs_env
282 : INTEGER, INTENT(IN) :: xas_method, iatom, istate
283 :
284 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_write_restart'
285 :
286 : CHARACTER(LEN=default_path_length) :: filename
287 : CHARACTER(LEN=default_string_length) :: my_middle
288 : INTEGER :: handle, ispin, nao, nexc_atoms, &
289 : nexc_search, nmo, output_unit, &
290 : rst_unit, xas_estate
291 : REAL(dp) :: occ_estate, xas_nelectron
292 704 : REAL(dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
293 : TYPE(cp_fm_type), POINTER :: mo_coeff
294 : TYPE(cp_logger_type), POINTER :: logger
295 704 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
296 704 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
297 704 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
298 : TYPE(section_vals_type), POINTER :: print_key
299 :
300 704 : CALL timeset(routineN, handle)
301 704 : NULLIFY (mos, logger, print_key, particle_set, qs_kind_set)
302 704 : logger => cp_get_default_logger()
303 :
304 : CALL get_xas_env(xas_env=xas_env, occ_estate=occ_estate, xas_estate=xas_estate, &
305 704 : xas_nelectron=xas_nelectron, nexc_search=nexc_search, nexc_atoms=nexc_atoms)
306 :
307 704 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
308 : xas_section, "PRINT%RESTART", used_print_key=print_key), &
309 : cp_p_file)) THEN
310 :
311 : output_unit = cp_print_key_unit_nr(logger, xas_section, &
312 604 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
313 :
314 604 : CALL get_qs_env(qs_env=qs_env, mos=mos)
315 :
316 : ! Open file
317 604 : rst_unit = -1
318 604 : my_middle = 'at'//TRIM(ADJUSTL(cp_to_string(iatom)))//'_st'//TRIM(ADJUSTL(cp_to_string(istate)))
319 : rst_unit = cp_print_key_unit_nr(logger, xas_section, "PRINT%RESTART", &
320 : extension=".rst", file_status="REPLACE", file_action="WRITE", &
321 604 : file_form="UNFORMATTED", middle_name=TRIM(my_middle))
322 :
323 : filename = cp_print_key_generate_filename(logger, print_key, &
324 : middle_name=TRIM(my_middle), extension=".rst", &
325 604 : my_local=.FALSE.)
326 :
327 604 : IF (output_unit > 0) THEN
328 : WRITE (UNIT=output_unit, FMT="(/,T10,A,I5,A,A,/)") &
329 302 : "Xas orbitals for the absorbing atom ", iatom, &
330 604 : " are written in ", TRIM(filename)
331 :
332 : END IF
333 :
334 : ! Write mos
335 604 : IF (rst_unit > 0) THEN
336 302 : WRITE (rst_unit) xas_method
337 302 : WRITE (rst_unit) nexc_search, nexc_atoms, occ_estate, xas_nelectron
338 302 : WRITE (rst_unit) xas_estate
339 : END IF
340 1812 : DO ispin = 1, SIZE(mos)
341 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo, &
342 1208 : eigenvalues=eigenvalues, occupation_numbers=occupation_numbers)
343 1208 : IF ((rst_unit > 0)) THEN
344 604 : WRITE (rst_unit) nao, nmo
345 604 : WRITE (rst_unit) eigenvalues(1:nmo), &
346 1208 : occupation_numbers(1:nmo)
347 : END IF
348 3020 : CALL cp_fm_write_unformatted(mo_coeff, rst_unit)
349 : END DO
350 :
351 : ! Close file
352 : CALL cp_print_key_finished_output(rst_unit, logger, xas_section, &
353 604 : "PRINT%RESTART")
354 : END IF
355 :
356 704 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
357 : xas_section, "PRINT%FULL_RESTART", used_print_key=print_key), &
358 : cp_p_file)) THEN
359 : rst_unit = cp_print_key_unit_nr(logger, xas_section, "PRINT%FULL_RESTART", &
360 : extension="_full.rst", file_status="REPLACE", file_action="WRITE", &
361 6 : file_form="UNFORMATTED", middle_name=TRIM(my_middle))
362 :
363 6 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set)
364 : CALL write_mo_set_low(mos, particle_set=particle_set, &
365 6 : qs_kind_set=qs_kind_set, ires=rst_unit)
366 6 : CALL cp_print_key_finished_output(rst_unit, logger, xas_section, "PRINT%FULL_RESTART")
367 :
368 : END IF
369 :
370 704 : CALL timestop(handle)
371 :
372 704 : END SUBROUTINE xas_write_restart
373 :
374 : !****f* xas_restart/xas_initialize_rho [1.0] *
375 :
376 : ! **************************************************************************************************
377 : !> \brief Once the mos and the occupation numbers are initialized
378 : !> the electronic density of the excited state can be calclated
379 : !> \param qs_env ...
380 : !> \param scf_env ...
381 : !> \param scf_control ...
382 : !> \par History
383 : !> 09-2006 MI created
384 : !> \author MI
385 : ! **************************************************************************************************
386 82 : SUBROUTINE xas_initialize_rho(qs_env, scf_env, scf_control)
387 :
388 : TYPE(qs_environment_type), POINTER :: qs_env
389 : TYPE(qs_scf_env_type), POINTER :: scf_env
390 : TYPE(scf_control_type), POINTER :: scf_control
391 :
392 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_initialize_rho'
393 :
394 : INTEGER :: handle, ispin, my_spin, nelectron
395 82 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
396 : TYPE(dft_control_type), POINTER :: dft_control
397 82 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
398 : TYPE(mp_para_env_type), POINTER :: para_env
399 : TYPE(qs_rho_type), POINTER :: rho
400 82 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
401 : TYPE(xas_environment_type), POINTER :: xas_env
402 :
403 82 : CALL timeset(routineN, handle)
404 :
405 82 : NULLIFY (mos, rho, xas_env, para_env, rho_ao)
406 :
407 : CALL get_qs_env(qs_env, &
408 : mos=mos, &
409 : rho=rho, &
410 : xas_env=xas_env, &
411 82 : para_env=para_env)
412 :
413 82 : my_spin = xas_env%spin_channel
414 82 : CALL qs_rho_get(rho, rho_ao=rho_ao)
415 246 : DO ispin = 1, SIZE(mos)
416 164 : IF (ispin == my_spin) THEN
417 82 : IF (xas_env%homo_occ == 0) THEN
418 2 : CALL get_mo_set(mos(ispin), nelectron=nelectron)
419 2 : nelectron = nelectron - 1
420 2 : CALL set_mo_set(mos(ispin), nelectron=nelectron)
421 : END IF
422 : CALL set_mo_occupation(mo_set=qs_env%mos(ispin), smear=scf_control%smear, &
423 82 : xas_env=xas_env)
424 : ELSE
425 82 : CALL set_mo_occupation(mo_set=qs_env%mos(ispin), smear=scf_control%smear)
426 : END IF
427 : CALL calculate_density_matrix(mo_set=mos(ispin), &
428 246 : density_matrix=rho_ao(ispin)%matrix)
429 : END DO
430 :
431 82 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
432 82 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
433 :
434 82 : IF (scf_env%mixing_method > 1) THEN
435 6 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
436 6 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
437 0 : CPABORT('TB Code not available')
438 6 : ELSE IF (dft_control%qs_control%semi_empirical) THEN
439 0 : CPABORT('SE Code not possible')
440 : ELSE
441 6 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
442 : CALL mixing_init(scf_env%mixing_method, rho, scf_env%mixing_store, &
443 6 : para_env, rho_atom=rho_atom)
444 : END IF
445 : END IF
446 :
447 82 : CALL timestop(handle)
448 :
449 82 : END SUBROUTINE xas_initialize_rho
450 :
451 : ! **************************************************************************************************
452 : !> \brief Find the index of the core orbital that has been excited by XAS
453 : !> \param xas_env ...
454 : !> \param mos ...
455 : !> \param matrix_s ...
456 : !> \par History
457 : !> 03-2010 MI created
458 : !> \author MI
459 : ! **************************************************************************************************
460 :
461 704 : SUBROUTINE find_excited_core_orbital(xas_env, mos, matrix_s)
462 :
463 : TYPE(xas_environment_type), POINTER :: xas_env
464 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
465 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
466 :
467 : INTEGER :: i, ic_max, ir_max, m, my_spin, n, nao, &
468 : nexc_search, nmo, xas_estate
469 704 : INTEGER, DIMENSION(:), POINTER :: col_indices
470 : REAL(dp) :: a_max, b_max, ip_energy, occ_estate
471 704 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
472 704 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer, vecbuffer2
473 : TYPE(cp_fm_type) :: fm_work
474 : TYPE(cp_fm_type), POINTER :: excvec_coeff, excvec_overlap, mo_coeff
475 :
476 704 : NULLIFY (excvec_coeff, excvec_overlap, mo_coeff)
477 : ! Some elements from the xas_env
478 : CALL get_xas_env(xas_env=xas_env, excvec_coeff=excvec_coeff, &
479 : excvec_overlap=excvec_overlap, nexc_search=nexc_search, &
480 704 : xas_estate=xas_estate, occ_estate=occ_estate, spin_channel=my_spin)
481 704 : CPASSERT(ASSOCIATED(excvec_overlap))
482 :
483 : CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff, nao=nao, nmo=nmo, &
484 704 : eigenvalues=eigenvalues, occupation_numbers=occupation_numbers)
485 2112 : ALLOCATE (vecbuffer(1, nao))
486 47928 : vecbuffer = 0.0_dp
487 2112 : ALLOCATE (vecbuffer2(1, nexc_search))
488 7448 : vecbuffer2 = 0.0_dp
489 :
490 : ! ** use the maximum overlap criterion to find the index of the excited orbital
491 704 : CALL cp_fm_create(fm_work, mo_coeff%matrix_struct)
492 704 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, fm_work, ncol=nmo)
493 : CALL parallel_gemm("T", "N", 1, xas_env%nexc_search, nao, 1.0_dp, excvec_coeff, &
494 704 : fm_work, 0.0_dp, excvec_overlap, b_first_col=1)
495 : CALL cp_fm_get_info(matrix=excvec_overlap, col_indices=col_indices, &
496 704 : nrow_global=m, ncol_global=n)
497 : CALL cp_fm_get_submatrix(excvec_overlap, vecbuffer2, 1, 1, &
498 704 : 1, nexc_search, transpose=.FALSE.)
499 704 : CALL cp_fm_release(fm_work)
500 :
501 704 : b_max = 0.0_dp
502 704 : ic_max = xas_estate
503 4076 : DO i = 1, nexc_search
504 3372 : a_max = ABS(vecbuffer2(1, i))
505 4076 : IF (a_max > b_max) THEN
506 1282 : ic_max = i
507 :
508 1282 : b_max = a_max
509 : END IF
510 : END DO
511 :
512 704 : IF (ic_max /= xas_estate) THEN
513 30 : ir_max = xas_estate
514 30 : xas_estate = ic_max
515 30 : occupation_numbers(xas_estate) = occ_estate
516 30 : occupation_numbers(ir_max) = 1.0_dp
517 : END IF
518 :
519 : ! Ionization Potential
520 704 : iP_energy = eigenvalues(xas_estate)
521 704 : CALL set_xas_env(xas_env=xas_env, xas_estate=xas_estate, ip_energy=ip_energy)
522 :
523 : CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, xas_estate, &
524 704 : nao, 1, transpose=.TRUE.)
525 : CALL cp_fm_set_submatrix(excvec_coeff, vecbuffer, 1, 1, &
526 704 : nao, 1, transpose=.TRUE.)
527 :
528 704 : DEALLOCATE (vecbuffer, vecbuffer2)
529 :
530 2816 : END SUBROUTINE find_excited_core_orbital
531 :
532 : END MODULE xas_restart
|
