Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
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) THEN
187 0 : CPABORT("READ XAS RESTART: restart with different XAS method is not possible.")
188 : END IF
189 6 : IF (nexc_atoms_read /= nexc_atoms) THEN
190 : CALL cp_abort(__LOCATION__, &
191 : "READ XAS RESTART: restart with different excited atoms "// &
192 0 : "is not possible. Start instead a new XAS run with the new set of atoms.")
193 : END IF
194 : END IF
195 :
196 12 : CALL para_env%bcast(xas_estate_read)
197 12 : CALL set_xas_env(xas_env=xas_env, xas_estate=xas_estate_read)
198 12 : estate = xas_estate_read
199 :
200 12 : CALL get_mo_set(mo_set=mos(my_spin), nao=nao)
201 36 : ALLOCATE (vecbuffer(1, nao))
202 :
203 36 : DO ispin = 1, SIZE(mos)
204 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, eigenvalues=eigenvalues, &
205 24 : occupation_numbers=occupation_numbers, mo_coeff=mo_coeff, nelectron=nelectron)
206 232 : eigenvalues = 0.0_dp
207 232 : occupation_numbers = 0.0_dp
208 24 : CALL cp_fm_set_all(mo_coeff, 0.0_dp)
209 24 : IF (para_env%is_source()) THEN
210 12 : READ (rst_unit) nao_read, nmo_read
211 12 : IF (nao /= nao_read) THEN
212 0 : CPABORT("To change basis is not possible. ")
213 : END IF
214 48 : ALLOCATE (eig_read(nmo_read), occ_read(nmo_read))
215 12 : eig_read = 0.0_dp
216 12 : occ_read = 0.0_dp
217 12 : nmo = MIN(nmo, nmo_read)
218 12 : READ (rst_unit) eig_read(1:nmo_read), occ_read(1:nmo_read)
219 116 : eigenvalues(1:nmo) = eig_read(1:nmo)
220 116 : occupation_numbers(1:nmo) = occ_read(1:nmo)
221 12 : IF (nmo_read > nmo) THEN
222 0 : IF (occupation_numbers(nmo) >= EPSILON(0.0_dp)) THEN
223 : CALL cp_warn(__LOCATION__, &
224 : "The number of occupied MOs on the restart unit is larger than "// &
225 0 : "the allocated MOs.")
226 : END IF
227 :
228 : END IF
229 12 : DEALLOCATE (eig_read, occ_read)
230 : END IF
231 440 : CALL para_env%bcast(eigenvalues)
232 440 : CALL para_env%bcast(occupation_numbers)
233 :
234 232 : DO i = 1, nmo
235 208 : IF (para_env%is_source()) THEN
236 5928 : READ (rst_unit) vecbuffer
237 : ELSE
238 3016 : vecbuffer(1, :) = 0.0_dp
239 : END IF
240 23504 : CALL para_env%bcast(vecbuffer)
241 : CALL cp_fm_set_submatrix(mo_coeff, &
242 232 : vecbuffer, 1, i, nao, 1, transpose=.TRUE.)
243 : END DO
244 : ! Skip extra MOs if there any
245 60 : IF (para_env%is_source()) THEN
246 12 : DO i = nmo + 1, nmo_read
247 12 : READ (rst_unit) vecbuffer
248 : END DO
249 : END IF
250 :
251 : END DO ! ispin
252 :
253 24 : DEALLOCATE (vecbuffer)
254 :
255 : ! nspin = SIZE(mos,1)
256 : ! DO ispin = 1,nspin
257 : ! ! ortho so that one can restart for different positions (basis sets?)
258 : ! NULLIFY(mo_coeff)
259 : ! CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff,homo=homo)
260 : ! CALL make_basis_sm(mo_coeff,homo,matrix_s(1)%matrix)
261 : ! END DO
262 : END IF !file_exist
263 :
264 12 : IF (para_env%is_source()) THEN
265 6 : IF (file_exists) CALL close_file(unit_number=rst_unit)
266 : END IF
267 :
268 12 : CALL timestop(handle)
269 :
270 24 : END SUBROUTINE xas_read_restart
271 :
272 : ! **************************************************************************************************
273 : !> \brief ...
274 : !> \param xas_env ...
275 : !> \param xas_section ...
276 : !> \param qs_env ...
277 : !> \param xas_method ...
278 : !> \param iatom ...
279 : !> \param istate ...
280 : ! **************************************************************************************************
281 1432 : SUBROUTINE xas_write_restart(xas_env, xas_section, qs_env, xas_method, iatom, istate)
282 :
283 : TYPE(xas_environment_type), POINTER :: xas_env
284 : TYPE(section_vals_type), POINTER :: xas_section
285 : TYPE(qs_environment_type), POINTER :: qs_env
286 : INTEGER, INTENT(IN) :: xas_method, iatom, istate
287 :
288 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_write_restart'
289 :
290 : CHARACTER(LEN=default_path_length) :: filename
291 : CHARACTER(LEN=default_string_length) :: my_middle
292 : INTEGER :: handle, ispin, nao, nexc_atoms, &
293 : nexc_search, nmo, output_unit, &
294 : rst_unit, xas_estate
295 : REAL(dp) :: occ_estate, xas_nelectron
296 716 : REAL(dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
297 : TYPE(cp_fm_type), POINTER :: mo_coeff
298 : TYPE(cp_logger_type), POINTER :: logger
299 716 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
300 716 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
301 716 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
302 : TYPE(section_vals_type), POINTER :: print_key
303 :
304 716 : CALL timeset(routineN, handle)
305 716 : NULLIFY (mos, logger, print_key, particle_set, qs_kind_set)
306 716 : logger => cp_get_default_logger()
307 :
308 : CALL get_xas_env(xas_env=xas_env, occ_estate=occ_estate, xas_estate=xas_estate, &
309 716 : xas_nelectron=xas_nelectron, nexc_search=nexc_search, nexc_atoms=nexc_atoms)
310 :
311 716 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
312 : xas_section, "PRINT%RESTART", used_print_key=print_key), &
313 : cp_p_file)) THEN
314 :
315 : output_unit = cp_print_key_unit_nr(logger, xas_section, &
316 604 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
317 :
318 604 : CALL get_qs_env(qs_env=qs_env, mos=mos)
319 :
320 : ! Open file
321 604 : rst_unit = -1
322 604 : my_middle = 'at'//TRIM(ADJUSTL(cp_to_string(iatom)))//'_st'//TRIM(ADJUSTL(cp_to_string(istate)))
323 : rst_unit = cp_print_key_unit_nr(logger, xas_section, "PRINT%RESTART", &
324 : extension=".rst", file_status="REPLACE", file_action="WRITE", &
325 604 : file_form="UNFORMATTED", middle_name=TRIM(my_middle))
326 :
327 : filename = cp_print_key_generate_filename(logger, print_key, &
328 : middle_name=TRIM(my_middle), extension=".rst", &
329 604 : my_local=.FALSE.)
330 :
331 604 : IF (output_unit > 0) THEN
332 : WRITE (UNIT=output_unit, FMT="(/,T10,A,I5,A,A,/)") &
333 302 : "Xas orbitals for the absorbing atom ", iatom, &
334 604 : " are written in ", TRIM(filename)
335 :
336 : END IF
337 :
338 : ! Write mos
339 604 : IF (rst_unit > 0) THEN
340 302 : WRITE (rst_unit) xas_method
341 302 : WRITE (rst_unit) nexc_search, nexc_atoms, occ_estate, xas_nelectron
342 302 : WRITE (rst_unit) xas_estate
343 : END IF
344 1812 : DO ispin = 1, SIZE(mos)
345 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo, &
346 1208 : eigenvalues=eigenvalues, occupation_numbers=occupation_numbers)
347 1208 : IF ((rst_unit > 0)) THEN
348 604 : WRITE (rst_unit) nao, nmo
349 6618 : WRITE (rst_unit) eigenvalues(1:nmo), &
350 7222 : occupation_numbers(1:nmo)
351 : END IF
352 3020 : CALL cp_fm_write_unformatted(mo_coeff, rst_unit)
353 : END DO
354 :
355 : ! Close file
356 : CALL cp_print_key_finished_output(rst_unit, logger, xas_section, &
357 604 : "PRINT%RESTART")
358 : END IF
359 :
360 716 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
361 : xas_section, "PRINT%FULL_RESTART", used_print_key=print_key), &
362 : cp_p_file)) THEN
363 : rst_unit = cp_print_key_unit_nr(logger, xas_section, "PRINT%FULL_RESTART", &
364 : extension="_full.rst", file_status="REPLACE", file_action="WRITE", &
365 6 : file_form="UNFORMATTED", middle_name=TRIM(my_middle))
366 :
367 6 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set)
368 : CALL write_mo_set_low(mos, particle_set=particle_set, &
369 6 : qs_kind_set=qs_kind_set, ires=rst_unit)
370 6 : CALL cp_print_key_finished_output(rst_unit, logger, xas_section, "PRINT%FULL_RESTART")
371 :
372 : END IF
373 :
374 716 : CALL timestop(handle)
375 :
376 716 : END SUBROUTINE xas_write_restart
377 :
378 : !****f* xas_restart/xas_initialize_rho [1.0] *
379 :
380 : ! **************************************************************************************************
381 : !> \brief Once the mos and the occupation numbers are initialized
382 : !> the electronic density of the excited state can be calclated
383 : !> \param qs_env ...
384 : !> \param scf_env ...
385 : !> \param scf_control ...
386 : !> \par History
387 : !> 09-2006 MI created
388 : !> \author MI
389 : ! **************************************************************************************************
390 82 : SUBROUTINE xas_initialize_rho(qs_env, scf_env, scf_control)
391 :
392 : TYPE(qs_environment_type), POINTER :: qs_env
393 : TYPE(qs_scf_env_type), POINTER :: scf_env
394 : TYPE(scf_control_type), POINTER :: scf_control
395 :
396 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_initialize_rho'
397 :
398 : INTEGER :: handle, ispin, my_spin, nelectron
399 82 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
400 : TYPE(dft_control_type), POINTER :: dft_control
401 82 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
402 : TYPE(mp_para_env_type), POINTER :: para_env
403 : TYPE(qs_rho_type), POINTER :: rho
404 82 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
405 : TYPE(xas_environment_type), POINTER :: xas_env
406 :
407 82 : CALL timeset(routineN, handle)
408 :
409 82 : NULLIFY (mos, rho, xas_env, para_env, rho_ao)
410 :
411 : CALL get_qs_env(qs_env, &
412 : mos=mos, &
413 : rho=rho, &
414 : xas_env=xas_env, &
415 82 : para_env=para_env)
416 :
417 82 : my_spin = xas_env%spin_channel
418 82 : CALL qs_rho_get(rho, rho_ao=rho_ao)
419 246 : DO ispin = 1, SIZE(mos)
420 164 : IF (ispin == my_spin) THEN
421 82 : IF (xas_env%homo_occ == 0) THEN
422 2 : CALL get_mo_set(mos(ispin), nelectron=nelectron)
423 2 : nelectron = nelectron - 1
424 2 : CALL set_mo_set(mos(ispin), nelectron=nelectron)
425 : END IF
426 : CALL set_mo_occupation(mo_set=qs_env%mos(ispin), smear=scf_control%smear, &
427 82 : xas_env=xas_env)
428 : ELSE
429 82 : CALL set_mo_occupation(mo_set=qs_env%mos(ispin), smear=scf_control%smear)
430 : END IF
431 : CALL calculate_density_matrix(mo_set=mos(ispin), &
432 246 : density_matrix=rho_ao(ispin)%matrix)
433 : END DO
434 :
435 82 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
436 82 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
437 :
438 82 : IF (scf_env%mixing_method > 1) THEN
439 6 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
440 6 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
441 0 : CPABORT('TB Code not available')
442 6 : ELSE IF (dft_control%qs_control%semi_empirical) THEN
443 0 : CPABORT('SE Code not possible')
444 : ELSE
445 6 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
446 : CALL mixing_init(scf_env%mixing_method, rho, scf_env%mixing_store, &
447 6 : para_env, rho_atom=rho_atom)
448 : END IF
449 : END IF
450 :
451 82 : CALL timestop(handle)
452 :
453 82 : END SUBROUTINE xas_initialize_rho
454 :
455 : ! **************************************************************************************************
456 : !> \brief Find the index of the core orbital that has been excited by XAS
457 : !> \param xas_env ...
458 : !> \param mos ...
459 : !> \param matrix_s ...
460 : !> \par History
461 : !> 03-2010 MI created
462 : !> \author MI
463 : ! **************************************************************************************************
464 :
465 716 : SUBROUTINE find_excited_core_orbital(xas_env, mos, matrix_s)
466 :
467 : TYPE(xas_environment_type), POINTER :: xas_env
468 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
469 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
470 :
471 : INTEGER :: i, ic_max, ir_max, m, my_spin, n, nao, &
472 : nexc_search, nmo, xas_estate
473 716 : INTEGER, DIMENSION(:), POINTER :: col_indices
474 : REAL(dp) :: a_max, b_max, ip_energy, occ_estate
475 716 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
476 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer, vecbuffer2
477 : TYPE(cp_fm_type) :: fm_work
478 : TYPE(cp_fm_type), POINTER :: excvec_coeff, excvec_overlap, mo_coeff
479 :
480 716 : NULLIFY (excvec_coeff, excvec_overlap, mo_coeff)
481 : ! Some elements from the xas_env
482 : CALL get_xas_env(xas_env=xas_env, excvec_coeff=excvec_coeff, &
483 : excvec_overlap=excvec_overlap, nexc_search=nexc_search, &
484 716 : xas_estate=xas_estate, occ_estate=occ_estate, spin_channel=my_spin)
485 716 : CPASSERT(ASSOCIATED(excvec_overlap))
486 :
487 : CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff, nao=nao, nmo=nmo, &
488 716 : eigenvalues=eigenvalues, occupation_numbers=occupation_numbers)
489 2148 : ALLOCATE (vecbuffer(1, nao))
490 50628 : vecbuffer = 0.0_dp
491 2148 : ALLOCATE (vecbuffer2(1, nexc_search))
492 7844 : vecbuffer2 = 0.0_dp
493 :
494 : ! ** use the maximum overlap criterion to find the index of the excited orbital
495 716 : CALL cp_fm_create(fm_work, mo_coeff%matrix_struct)
496 716 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, fm_work, ncol=nmo)
497 : CALL parallel_gemm("T", "N", 1, xas_env%nexc_search, nao, 1.0_dp, excvec_coeff, &
498 716 : fm_work, 0.0_dp, excvec_overlap, b_first_col=1)
499 : CALL cp_fm_get_info(matrix=excvec_overlap, col_indices=col_indices, &
500 716 : nrow_global=m, ncol_global=n)
501 : CALL cp_fm_get_submatrix(excvec_overlap, vecbuffer2, 1, 1, &
502 716 : 1, nexc_search, transpose=.FALSE.)
503 716 : CALL cp_fm_release(fm_work)
504 :
505 716 : b_max = 0.0_dp
506 716 : ic_max = xas_estate
507 4280 : DO i = 1, nexc_search
508 3564 : a_max = ABS(vecbuffer2(1, i))
509 4280 : IF (a_max > b_max) THEN
510 1294 : ic_max = i
511 :
512 1294 : b_max = a_max
513 : END IF
514 : END DO
515 :
516 716 : IF (ic_max /= xas_estate) THEN
517 30 : ir_max = xas_estate
518 30 : xas_estate = ic_max
519 30 : occupation_numbers(xas_estate) = occ_estate
520 30 : occupation_numbers(ir_max) = 1.0_dp
521 : END IF
522 :
523 : ! Ionization Potential
524 716 : iP_energy = eigenvalues(xas_estate)
525 716 : CALL set_xas_env(xas_env=xas_env, xas_estate=xas_estate, ip_energy=ip_energy)
526 :
527 : CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, xas_estate, &
528 716 : nao, 1, transpose=.TRUE.)
529 : CALL cp_fm_set_submatrix(excvec_coeff, vecbuffer, 1, 1, &
530 716 : nao, 1, transpose=.TRUE.)
531 :
532 716 : DEALLOCATE (vecbuffer, vecbuffer2)
533 :
534 2864 : END SUBROUTINE find_excited_core_orbital
535 :
536 : END MODULE xas_restart
|
