Line data Source code
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 :
8 : ! **************************************************************************************************
9 : !> \brief xas_scf for the tp method
10 : !> It is repeaated for every atom that have to be excited
11 : !> \par History
12 : !> created 05.2005
13 : !> \author MI (05.2005)
14 : ! **************************************************************************************************
15 : MODULE xas_tp_scf
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind
18 : USE cell_types, ONLY: cell_type,&
19 : pbc
20 : USE cp_array_utils, ONLY: cp_2d_r_p_type
21 : USE cp_blacs_env, ONLY: cp_blacs_env_type
22 : USE cp_control_types, ONLY: dft_control_type
23 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
24 : dbcsr_p_type
25 : USE cp_external_control, ONLY: external_control
26 : USE cp_fm_types, ONLY: cp_fm_get_submatrix,&
27 : cp_fm_init_random,&
28 : cp_fm_set_submatrix,&
29 : cp_fm_to_fm,&
30 : cp_fm_type
31 : USE cp_log_handling, ONLY: cp_get_default_logger,&
32 : cp_logger_get_default_io_unit,&
33 : cp_logger_type,&
34 : cp_to_string
35 : USE cp_output_handling, ONLY: cp_add_iter_level,&
36 : cp_iterate,&
37 : cp_p_file,&
38 : cp_print_key_finished_output,&
39 : cp_print_key_should_output,&
40 : cp_print_key_unit_nr,&
41 : cp_rm_iter_level
42 : USE input_constants, ONLY: ot_precond_full_kinetic,&
43 : ot_precond_solver_default,&
44 : xas_dscf
45 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
46 : section_vals_type
47 : USE kinds, ONLY: dp
48 : USE machine, ONLY: m_flush,&
49 : m_walltime
50 : USE message_passing, ONLY: mp_para_env_type
51 : USE particle_methods, ONLY: get_particle_set
52 : USE particle_types, ONLY: particle_type
53 : USE preconditioner, ONLY: make_preconditioner
54 : USE preconditioner_types, ONLY: destroy_preconditioner,&
55 : init_preconditioner,&
56 : preconditioner_type
57 : USE qs_charges_types, ONLY: qs_charges_type
58 : USE qs_density_matrices, ONLY: calculate_density_matrix
59 : USE qs_density_mixing_types, ONLY: broyden_mixing_nr,&
60 : direct_mixing_nr,&
61 : gspace_mixing_nr,&
62 : multisecant_mixing_nr,&
63 : no_mixing_nr,&
64 : pulay_mixing_nr
65 : USE qs_energy_types, ONLY: qs_energy_type
66 : USE qs_environment_types, ONLY: get_qs_env,&
67 : qs_environment_type
68 : USE qs_gspace_mixing, ONLY: gspace_mixing
69 : USE qs_kind_types, ONLY: qs_kind_type
70 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
71 : USE qs_ks_types, ONLY: qs_ks_did_change,&
72 : qs_ks_env_type
73 : USE qs_loc_main, ONLY: qs_loc_driver
74 : USE qs_loc_types, ONLY: get_qs_loc_env,&
75 : localized_wfn_control_type,&
76 : qs_loc_env_type
77 : USE qs_mixing_utils, ONLY: self_consistency_check
78 : USE qs_mo_io, ONLY: write_mo_set_to_restart
79 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
80 : USE qs_mo_occupation, ONLY: set_mo_occupation
81 : USE qs_mo_types, ONLY: get_mo_set,&
82 : mo_set_type
83 : USE qs_ot_eigensolver, ONLY: ot_eigensolver
84 : USE qs_rho_methods, ONLY: qs_rho_update_rho
85 : USE qs_rho_types, ONLY: qs_rho_get,&
86 : qs_rho_type
87 : USE qs_scf, ONLY: scf_env_cleanup,&
88 : scf_env_do_scf
89 : USE qs_scf_diagonalization, ONLY: general_eigenproblem
90 : USE qs_scf_initialization, ONLY: qs_scf_env_initialize
91 : USE qs_scf_methods, ONLY: scf_env_density_mixing
92 : USE qs_scf_output, ONLY: qs_scf_print_summary
93 : USE qs_scf_types, ONLY: general_diag_method_nr,&
94 : qs_scf_env_type
95 : USE scf_control_types, ONLY: scf_control_type
96 : USE xas_control, ONLY: xas_control_type
97 : USE xas_env_types, ONLY: get_xas_env,&
98 : set_xas_env,&
99 : xas_environment_type
100 : USE xas_restart, ONLY: find_excited_core_orbital,&
101 : xas_write_restart
102 : #include "./base/base_uses.f90"
103 :
104 : IMPLICIT NONE
105 : PRIVATE
106 :
107 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_tp_scf'
108 :
109 : ! *** Public subroutines ***
110 :
111 : PUBLIC :: xas_do_tp_scf, xes_scf_once
112 :
113 : CONTAINS
114 :
115 : ! **************************************************************************************************
116 : !> \brief perform an scf loop to calculate the xas spectrum
117 : !> given by the excitation of a inner state of a selected atom
118 : !> by using the transition potential method
119 : !> \param dft_control ...
120 : !> \param xas_env the environment for XAS calculations
121 : !> \param iatom ...
122 : !> \param istate keeps track of the state index for restart writing
123 : !> \param scf_env the scf_env where to perform the scf procedure
124 : !> \param qs_env the qs_env, the scf_env and xas_env live in
125 : !> \param xas_section ...
126 : !> \param scf_section ...
127 : !> \param converged ...
128 : !> \param should_stop ...
129 : !> \par History
130 : !> 05.2005 created [MI]
131 : !> \author MI
132 : ! **************************************************************************************************
133 160 : SUBROUTINE xas_do_tp_scf(dft_control, xas_env, iatom, istate, scf_env, qs_env, &
134 : xas_section, scf_section, converged, should_stop)
135 :
136 : TYPE(dft_control_type), POINTER :: dft_control
137 : TYPE(xas_environment_type), POINTER :: xas_env
138 : INTEGER, INTENT(IN) :: iatom, istate
139 : TYPE(qs_scf_env_type), POINTER :: scf_env
140 : TYPE(qs_environment_type), POINTER :: qs_env
141 : TYPE(section_vals_type), POINTER :: xas_section, scf_section
142 : LOGICAL, INTENT(OUT) :: converged, should_stop
143 :
144 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_do_tp_scf'
145 :
146 : INTEGER :: handle, handle2, hole_spin, ispin, &
147 : iter_count, my_spin, nspin, occ_spin, &
148 : output_unit
149 : LOGICAL :: diis_step, energy_only, exit_loop, gapw
150 : REAL(KIND=dp) :: t1, t2
151 80 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
152 : TYPE(cp_logger_type), POINTER :: logger
153 80 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
154 80 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
155 80 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
156 : TYPE(mp_para_env_type), POINTER :: para_env
157 80 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
158 : TYPE(qs_charges_type), POINTER :: qs_charges
159 : TYPE(qs_energy_type), POINTER :: energy
160 80 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
161 : TYPE(qs_ks_env_type), POINTER :: ks_env
162 : TYPE(qs_rho_type), POINTER :: rho
163 : TYPE(scf_control_type), POINTER :: scf_control
164 : TYPE(section_vals_type), POINTER :: dft_section
165 : TYPE(xas_control_type), POINTER :: xas_control
166 :
167 80 : CALL timeset(routineN, handle)
168 80 : NULLIFY (xas_control, matrix_s, matrix_ks, para_env, rho_ao_kp)
169 80 : NULLIFY (rho, energy, scf_control, logger, ks_env, mos, atomic_kind_set)
170 80 : NULLIFY (qs_charges, particle_set, qs_kind_set)
171 :
172 80 : logger => cp_get_default_logger()
173 80 : t1 = m_walltime()
174 80 : converged = .TRUE.
175 :
176 80 : CPASSERT(ASSOCIATED(xas_env))
177 80 : CPASSERT(ASSOCIATED(scf_env))
178 80 : CPASSERT(ASSOCIATED(qs_env))
179 :
180 : CALL get_qs_env(qs_env=qs_env, &
181 : atomic_kind_set=atomic_kind_set, &
182 : matrix_s=matrix_s, energy=energy, &
183 : qs_charges=qs_charges, &
184 80 : ks_env=ks_env, para_env=para_env)
185 :
186 80 : CALL get_xas_env(xas_env, scf_control=scf_control, spin_channel=my_spin)
187 :
188 80 : energy_only = .FALSE.
189 : output_unit = cp_print_key_unit_nr(logger, xas_section, "PRINT%PROGRAM_RUN_INFO", &
190 80 : extension=".xasLog")
191 80 : IF (output_unit > 0) THEN
192 40 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") "XAS_TP_SCF WAVEFUNCTION OPTIMIZATION"
193 : END IF
194 :
195 : ! GAPW method must be used
196 80 : gapw = dft_control%qs_control%gapw
197 80 : CPASSERT(gapw)
198 80 : xas_control => dft_control%xas_control
199 :
200 80 : CALL cp_add_iter_level(logger%iter_info, "XAS_SCF")
201 :
202 80 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks, rho=rho, mos=mos)
203 80 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
204 :
205 80 : iter_count = 0
206 80 : diis_step = .FALSE.
207 80 : nspin = dft_control%nspins
208 :
209 80 : IF (output_unit > 0) THEN
210 : WRITE (UNIT=output_unit, &
211 : FMT="(/,T3,A,T12,A,T31,A,T40,A,T60,A,T75,A/,T3,A)") &
212 40 : "Step", "Update method", "Time", "Convergence", "Total energy", "Change", &
213 80 : REPEAT("-", 77)
214 : END IF
215 :
216 : ! *** SCF loop ***
217 :
218 80 : energy%tot_old = 0.0_dp
219 1408 : scf_loop: DO
220 704 : CALL timeset(routineN//"_inner_loop", handle2)
221 :
222 704 : exit_loop = .FALSE.
223 704 : IF (output_unit > 0) CALL m_flush(output_unit)
224 :
225 704 : iter_count = iter_count + 1
226 704 : CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=iter_count)
227 :
228 : ! ** here qs_env%rho%rho_r and qs_env%rho%rho_g should be up to date
229 :
230 704 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=energy_only)
231 :
232 704 : SELECT CASE (scf_env%method)
233 : CASE DEFAULT
234 : CALL cp_abort(__LOCATION__, &
235 : "unknown scf method method for core level spectroscopy"// &
236 0 : cp_to_string(scf_env%method))
237 :
238 : CASE (general_diag_method_nr) ! diagonalisation (default)
239 :
240 704 : scf_env%iter_count = iter_count
241 : CALL general_eigenproblem(scf_env, mos, matrix_ks, &
242 : matrix_s, scf_control, scf_section, &
243 704 : diis_step)
244 :
245 704 : CALL find_excited_core_orbital(xas_env, mos, matrix_s)
246 :
247 : ! set occupation for the respective spin channels
248 704 : IF (my_spin == 1) THEN
249 : hole_spin = 1 ! hole is generated in channel 1
250 : occ_spin = 2
251 : ELSE
252 6 : hole_spin = 2
253 6 : occ_spin = 1
254 : END IF
255 : CALL set_mo_occupation(mo_set=mos(hole_spin), &
256 : smear=scf_control%smear, &
257 704 : xas_env=xas_env)
258 : CALL set_mo_occupation(mo_set=mos(occ_spin), &
259 704 : smear=scf_control%smear)
260 2112 : DO ispin = 1, nspin
261 : ! does not yet handle k-points
262 : CALL calculate_density_matrix(mos(ispin), &
263 2112 : scf_env%p_mix_new(ispin, 1)%matrix)
264 : END DO
265 704 : energy%kTS = 0.0_dp
266 704 : energy%efermi = 0.0_dp
267 2112 : DO ispin = 1, nspin
268 1408 : energy%kTS = energy%kTS + mos(ispin)%kTS
269 2112 : energy%efermi = energy%efermi + mos(ispin)%mu
270 : END DO
271 1408 : energy%efermi = energy%efermi/REAL(nspin, KIND=dp)
272 :
273 : END SELECT
274 :
275 1396 : SELECT CASE (scf_env%mixing_method)
276 : CASE (direct_mixing_nr)
277 : CALL scf_env_density_mixing(scf_env%p_mix_new, &
278 : scf_env%mixing_store, rho_ao_kp, para_env, scf_env%iter_delta, scf_env%iter_count, &
279 692 : diis=diis_step)
280 : CASE (gspace_mixing_nr, pulay_mixing_nr, broyden_mixing_nr, &
281 : multisecant_mixing_nr)
282 : ! Compute the difference p_out-p_in
283 : CALL self_consistency_check(rho_ao_kp, scf_env%p_delta, para_env, scf_env%p_mix_new, &
284 12 : delta=scf_env%iter_delta)
285 : CASE (no_mixing_nr)
286 : CASE DEFAULT
287 : CALL cp_abort(__LOCATION__, &
288 : "unknown scf mixing method: "// &
289 704 : cp_to_string(scf_env%mixing_method))
290 : END SELECT
291 :
292 704 : t2 = m_walltime()
293 :
294 704 : IF (output_unit > 0 .AND. scf_env%print_iter_line) THEN
295 : WRITE (UNIT=output_unit, &
296 : FMT="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,F14.8,1X,F20.10,1X,ES9.2)") &
297 352 : iter_count, TRIM(scf_env%iter_method), &
298 352 : scf_env%iter_param, t2 - t1, scf_env%iter_delta, energy%total, &
299 704 : energy%total - energy%tot_old
300 : END IF
301 704 : energy%tot_old = energy%total
302 :
303 : ! ** convergence check
304 : CALL external_control(should_stop, "XASSCF", target_time=qs_env%target_time, &
305 704 : start_time=qs_env%start_time)
306 704 : IF (scf_env%iter_delta < scf_control%eps_scf) THEN
307 46 : IF (output_unit > 0) THEN
308 : WRITE (UNIT=output_unit, FMT="(/,T3,A,I5,A/)") &
309 23 : "*** SCF run converged in ", iter_count, " steps ***"
310 : END IF
311 : exit_loop = .TRUE.
312 658 : ELSE IF (should_stop .OR. iter_count == scf_control%max_scf) THEN
313 34 : IF (output_unit > 0) THEN
314 : WRITE (UNIT=output_unit, FMT="(/,T3,A,/)") &
315 17 : "*** SCF run NOT converged ***"
316 : END IF
317 34 : converged = .FALSE.
318 : exit_loop = .TRUE.
319 : END IF
320 : ! *** Exit if we have finished with the SCF inner loop ***
321 : IF (exit_loop) THEN
322 : ! now, print out energies and charges corresponding to the obtained wfn
323 : ! (this actually is not 100% consistent at this point)!
324 80 : CALL qs_scf_print_summary(output_unit, qs_env)
325 80 : CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=iter_count)
326 : END IF
327 :
328 : ! ** Write restart file **
329 : CALL xas_write_restart(xas_env, xas_section, qs_env, xas_control%xas_method, &
330 704 : iatom, istate)
331 :
332 704 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
333 704 : CALL get_qs_env(qs_env=qs_env, mos=mos, particle_set=particle_set, qs_kind_set=qs_kind_set)
334 704 : CALL write_mo_set_to_restart(mos, particle_set, dft_section=dft_section, qs_kind_set=qs_kind_set)
335 :
336 704 : IF (exit_loop) THEN
337 80 : CALL timestop(handle2)
338 : EXIT scf_loop
339 : END IF
340 :
341 624 : IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, &
342 624 : xas_section, "PRINT%ITERATION_INFO/TIME_CUMUL"), cp_p_file)) t1 = m_walltime()
343 :
344 : ! *** mixing methods have the new density matrix in p_mix_new
345 624 : IF (scf_env%mixing_method > 0) THEN
346 1872 : DO ispin = 1, nspin
347 : ! does not yet handle k-points
348 1872 : CALL dbcsr_copy(rho_ao_kp(ispin, 1)%matrix, scf_env%p_mix_new(ispin, 1)%matrix)
349 : END DO
350 : END IF
351 :
352 : ! ** update qs_env%rho
353 624 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
354 624 : IF (scf_env%mixing_method >= gspace_mixing_nr) THEN
355 : CALL gspace_mixing(qs_env, scf_env%mixing_method, scf_env%mixing_store, rho, para_env, &
356 6 : scf_env%iter_count)
357 : END IF
358 :
359 624 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
360 624 : CALL timestop(handle2)
361 :
362 : END DO scf_loop
363 :
364 80 : IF (output_unit > 0) THEN
365 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T55,F25.14))") &
366 40 : "Ionization potential of the excited atom: ", xas_env%IP_energy
367 40 : CALL m_flush(output_unit)
368 : END IF
369 :
370 80 : CALL para_env%sync()
371 80 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
372 :
373 80 : CALL cls_prepare_states(xas_control, xas_env, qs_env, iatom, xas_section, output_unit)
374 :
375 80 : CALL para_env%sync()
376 :
377 : CALL cp_print_key_finished_output(output_unit, logger, xas_section, &
378 80 : "PRINT%PROGRAM_RUN_INFO")
379 80 : CALL cp_rm_iter_level(logger%iter_info, "XAS_SCF")
380 :
381 80 : CALL timestop(handle)
382 :
383 80 : END SUBROUTINE xas_do_tp_scf
384 : ! **************************************************************************************************
385 : !> \brief Post processing of the optimized wfn in XAS scheme, as preparation for
386 : !> the calculation of the spectrum
387 : !> \param xas_control ...
388 : !> \param xas_env the environment for XAS calculations
389 : !> \param qs_env the qs_env, the scf_env and xas_env live in
390 : !> \param iatom index of the excited atom
391 : !> \param xas_section ...
392 : !> \param output_unit ...
393 : !> \par History
394 : !> 05.2005 created [MI]
395 : !> \author MI
396 : ! **************************************************************************************************
397 80 : SUBROUTINE cls_prepare_states(xas_control, xas_env, qs_env, iatom, xas_section, output_unit)
398 :
399 : TYPE(xas_control_type) :: xas_control
400 : TYPE(xas_environment_type), POINTER :: xas_env
401 : TYPE(qs_environment_type), POINTER :: qs_env
402 : INTEGER, INTENT(IN) :: iatom
403 : TYPE(section_vals_type), POINTER :: xas_section
404 : INTEGER, INTENT(IN) :: output_unit
405 :
406 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cls_prepare_states'
407 :
408 : INTEGER :: handle, i, ikind, isgf, ispin, istate, j, my_kind, my_spin, my_state, nao, natom, &
409 : nexc_search, nmo, nvirtual2, uno_iter, xas_estate
410 80 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf
411 80 : INTEGER, DIMENSION(:), POINTER :: mykind_of_kind
412 80 : REAL(dp), DIMENSION(:, :), POINTER :: centers_wfn
413 : REAL(KIND=dp) :: component, dist, max_overlap, ra(3), &
414 : rac(3), rc(3), sto_state_overlap, &
415 : uno_eps
416 80 : REAL(KIND=dp), DIMENSION(:), POINTER :: all_evals, eigenvalues, uno_evals
417 80 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer, vecbuffer2
418 : TYPE(atomic_kind_type), POINTER :: atomic_kind
419 : TYPE(cell_type), POINTER :: cell
420 80 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: stogto_overlap
421 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
422 : TYPE(cp_fm_type), POINTER :: all_vectors, excvec_coeff, &
423 : excvec_overlap, mo_coeff, uno_orbs
424 80 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kinetic, matrix_ks, matrix_s
425 : TYPE(dft_control_type), POINTER :: dft_control
426 : TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
427 80 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
428 : TYPE(mp_para_env_type), POINTER :: para_env
429 80 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
430 : TYPE(preconditioner_type), POINTER :: local_preconditioner
431 80 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
432 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
433 : TYPE(scf_control_type), POINTER :: scf_control
434 : TYPE(section_vals_type), POINTER :: loc_section, print_loc_section
435 :
436 80 : CALL timeset(routineN, handle)
437 :
438 80 : NULLIFY (atomic_kind, dft_control, matrix_s, matrix_ks, qs_kind_set, kinetic)
439 80 : NULLIFY (cell, particle_set, local_preconditioner, vecbuffer, vecbuffer2)
440 80 : NULLIFY (dft_control, loc_section, mos, mo_coeff, eigenvalues)
441 80 : NULLIFY (centers_wfn, mykind_of_kind, qs_loc_env, localized_wfn_control, stogto_overlap)
442 80 : NULLIFY (all_evals, all_vectors, excvec_coeff, excvec_overlap, uno_evals, para_env, blacs_env)
443 :
444 80 : CPASSERT(ASSOCIATED(xas_env))
445 :
446 : CALL get_qs_env(qs_env, &
447 : cell=cell, &
448 : dft_control=dft_control, &
449 : matrix_ks=matrix_ks, &
450 : matrix_s=matrix_s, &
451 : kinetic=kinetic, &
452 : mos=mos, &
453 : particle_set=particle_set, &
454 : qs_kind_set=qs_kind_set, &
455 : para_env=para_env, &
456 80 : blacs_env=blacs_env)
457 :
458 : ! Some elements from the xas_env
459 : CALL get_xas_env(xas_env=xas_env, &
460 : all_vectors=all_vectors, all_evals=all_evals, &
461 : excvec_coeff=excvec_coeff, &
462 : nvirtual2=nvirtual2, xas_estate=xas_estate, &
463 : excvec_overlap=excvec_overlap, nexc_search=nexc_search, &
464 80 : spin_channel=my_spin, scf_control=scf_control)
465 80 : CPASSERT(ASSOCIATED(excvec_overlap))
466 :
467 : CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff, nao=nao, &
468 80 : eigenvalues=eigenvalues)
469 :
470 240 : ALLOCATE (vecbuffer(1, nao))
471 7160 : vecbuffer = 0.0_dp
472 160 : ALLOCATE (vecbuffer2(1, nao))
473 7160 : vecbuffer2 = 0.0_dp
474 80 : natom = SIZE(particle_set, 1)
475 240 : ALLOCATE (first_sgf(natom))
476 80 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf)
477 240 : ALLOCATE (centers_wfn(3, nexc_search))
478 1712 : centers_wfn = 0.0_dp
479 :
480 : ! Possible only for emission only
481 80 : IF (scf_control%use_ot) THEN
482 0 : IF (output_unit > 0) THEN
483 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") " Eigenstates are derived "// &
484 : "from the MOs optimized by OT. Follows localization of the core states"// &
485 0 : " to identify the excited orbital."
486 : END IF
487 :
488 : CALL get_xas_env(xas_env=xas_env, &
489 : mykind_of_kind=mykind_of_kind, qs_loc_env=qs_loc_env, &
490 0 : stogto_overlap=stogto_overlap)
491 : CALL get_qs_loc_env(qs_loc_env=qs_loc_env, &
492 0 : localized_wfn_control=localized_wfn_control)
493 0 : loc_section => section_vals_get_subs_vals(xas_section, "LOCALIZE")
494 0 : print_loc_section => section_vals_get_subs_vals(loc_section, "PRINT")
495 0 : CALL qs_loc_driver(qs_env, qs_loc_env, print_loc_section, myspin=1)
496 0 : ra(1:3) = particle_set(iatom)%r(1:3)
497 :
498 0 : NULLIFY (atomic_kind)
499 0 : atomic_kind => particle_set(iatom)%atomic_kind
500 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
501 0 : kind_number=ikind)
502 0 : my_kind = mykind_of_kind(ikind)
503 :
504 0 : my_state = 0
505 : CALL cp_fm_get_submatrix(mo_coeff, vecbuffer2, 1, my_state, &
506 0 : nao, 1, transpose=.TRUE.)
507 :
508 : ! Rotate the wfn to get the eigenstate of the KS hamiltonian
509 : ! Only ispin=1 should be needed
510 0 : DO ispin = 1, dft_control%nspins
511 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo, &
512 0 : eigenvalues=eigenvalues)
513 : CALL calculate_subspace_eigenvalues(mo_coeff, &
514 : matrix_ks(ispin)%matrix, eigenvalues, &
515 0 : do_rotation=.TRUE.)
516 : END DO ! ispin
517 :
518 : !Search for the core state to be excited
519 0 : max_overlap = 0.0_dp
520 0 : DO istate = 1, nexc_search
521 0 : centers_wfn(1, istate) = localized_wfn_control%centers_set(my_spin)%array(1, istate)
522 0 : centers_wfn(2, istate) = localized_wfn_control%centers_set(my_spin)%array(2, istate)
523 0 : centers_wfn(3, istate) = localized_wfn_control%centers_set(my_spin)%array(3, istate)
524 :
525 0 : rc(1:3) = centers_wfn(1:3, istate)
526 0 : rac = pbc(ra, rc, cell)
527 0 : dist = rac(1)*rac(1) + rac(2)*rac(2) + rac(3)*rac(3)
528 :
529 0 : IF (dist < 1.0_dp) THEN
530 : CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, istate, &
531 0 : nao, 1, transpose=.TRUE.)
532 0 : sto_state_overlap = 0.0_dp
533 0 : DO i = 1, SIZE(stogto_overlap(my_kind)%array, 1)
534 0 : component = 0.0_dp
535 0 : DO j = 1, SIZE(stogto_overlap(my_kind)%array, 2)
536 0 : isgf = first_sgf(iatom) + j - 1
537 0 : component = component + stogto_overlap(my_kind)%array(i, j)*vecbuffer(1, isgf)
538 : END DO ! j size
539 : sto_state_overlap = sto_state_overlap + &
540 0 : component*component
541 : END DO ! i size
542 :
543 0 : IF (sto_state_overlap .GT. max_overlap) THEN
544 0 : max_overlap = sto_state_overlap
545 0 : my_state = istate
546 : END IF
547 : END IF
548 0 : xas_estate = my_state
549 : END DO ! istate
550 :
551 0 : CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff)
552 : CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, xas_estate, &
553 0 : nao, 1, transpose=.TRUE.)
554 : CALL cp_fm_set_submatrix(excvec_coeff, vecbuffer2, 1, 1, &
555 0 : nao, 1, transpose=.TRUE.)
556 : !
557 : END IF
558 :
559 80 : CALL para_env%sync()
560 : !Calculate the virtual states from the KS matrix matrix_ks(1)
561 80 : IF (nvirtual2 .GT. 0) THEN
562 76 : NULLIFY (mo_coeff)
563 76 : CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff, nmo=nmo)
564 76 : IF (output_unit > 0) THEN
565 38 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A,I5,A,I6,A)") " Calculation of ", nvirtual2, &
566 : " additional virtual states of the subspace complementary to the"// &
567 76 : " lowest ", nmo, " states"
568 : END IF
569 :
570 76 : NULLIFY (uno_orbs, uno_evals, local_preconditioner)
571 76 : ALLOCATE (local_preconditioner)
572 : CALL init_preconditioner(local_preconditioner, para_env=para_env, &
573 76 : blacs_env=blacs_env)
574 :
575 : CALL make_preconditioner(local_preconditioner, &
576 : precon_type=ot_precond_full_kinetic, &
577 : solver_type=ot_precond_solver_default, &
578 : matrix_h=matrix_ks(my_spin)%matrix, &
579 : matrix_s=matrix_s(1)%matrix, &
580 : matrix_t=kinetic(1)%matrix, &
581 : convert_precond_to_dbcsr=.TRUE., &
582 76 : mo_set=mos(my_spin), energy_gap=0.2_dp)
583 :
584 : CALL get_xas_env(xas_env=xas_env, unoccupied_orbs=uno_orbs, &
585 76 : unoccupied_evals=uno_evals, unoccupied_eps=uno_eps, unoccupied_max_iter=uno_iter)
586 76 : CALL cp_fm_init_random(uno_orbs, nvirtual2)
587 :
588 : CALL ot_eigensolver(matrix_h=matrix_ks(my_spin)%matrix, matrix_s=matrix_s(1)%matrix, &
589 : matrix_c_fm=uno_orbs, matrix_orthogonal_space_fm=mo_coeff, &
590 : preconditioner=local_preconditioner, eps_gradient=uno_eps, &
591 76 : iter_max=uno_iter, size_ortho_space=nmo)
592 :
593 : CALL calculate_subspace_eigenvalues(uno_orbs, matrix_ks(my_spin)%matrix, &
594 76 : uno_evals, do_rotation=.TRUE.)
595 76 : CALL destroy_preconditioner(local_preconditioner)
596 :
597 76 : DEALLOCATE (local_preconditioner)
598 : END IF
599 :
600 80 : CALL para_env%sync()
601 : ! Prapare arrays for the calculation of the spectrum
602 80 : IF (.NOT. xas_control%xas_method == xas_dscf) THEN
603 : ! Copy the final vectors in the array
604 80 : NULLIFY (all_vectors, all_evals)
605 : CALL get_xas_env(xas_env=xas_env, all_vectors=all_vectors, &
606 80 : all_evals=all_evals)
607 : CALL get_mo_set(mos(my_spin), eigenvalues=eigenvalues, mo_coeff=mo_coeff, &
608 80 : nmo=nmo)
609 :
610 : CALL cp_fm_to_fm(mo_coeff, all_vectors, ncol=nmo, &
611 80 : source_start=1, target_start=1)
612 1230 : DO istate = 1, nmo
613 1230 : all_evals(istate) = eigenvalues(istate)
614 : END DO
615 80 : IF (nvirtual2 .GT. 0) THEN
616 : CALL cp_fm_to_fm(uno_orbs, all_vectors, ncol=nvirtual2, &
617 76 : source_start=1, target_start=1 + nmo)
618 2128 : DO istate = 1, nvirtual2
619 2128 : all_evals(istate + nmo) = uno_evals(istate)
620 : END DO
621 : END IF
622 : END IF
623 :
624 80 : DEALLOCATE (vecbuffer)
625 80 : DEALLOCATE (vecbuffer2)
626 80 : DEALLOCATE (centers_wfn, first_sgf)
627 :
628 80 : CALL timestop(handle)
629 :
630 240 : END SUBROUTINE cls_prepare_states
631 :
632 : ! **************************************************************************************************
633 : !> \brief SCF for emission spectra calculations: vacancy in valence
634 : !> \param qs_env the qs_env, the scf_env and xas_env live in
635 : !> \param xas_env the environment for XAS calculations
636 : !> \param converged ...
637 : !> \param should_stop ...
638 : !> \par History
639 : !> 05.2005 created [MI]
640 : !> \author MI
641 : ! **************************************************************************************************
642 20 : SUBROUTINE xes_scf_once(qs_env, xas_env, converged, should_stop)
643 :
644 : TYPE(qs_environment_type), POINTER :: qs_env
645 : TYPE(xas_environment_type), POINTER :: xas_env
646 : LOGICAL, INTENT(OUT) :: converged, should_stop
647 :
648 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xes_scf_once'
649 :
650 : INTEGER :: handle, ispin, istate, my_spin, nmo, &
651 : nvirtual, nvirtual2, output_unit, &
652 : tsteps
653 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: all_evals, eigenvalues
654 : TYPE(cp_fm_type), POINTER :: all_vectors, mo_coeff
655 : TYPE(cp_logger_type), POINTER :: logger
656 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
657 : TYPE(dft_control_type), POINTER :: dft_control
658 4 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
659 : TYPE(mp_para_env_type), POINTER :: para_env
660 : TYPE(qs_scf_env_type), POINTER :: scf_env
661 : TYPE(scf_control_type), POINTER :: scf_control
662 : TYPE(section_vals_type), POINTER :: dft_section, scf_section, xas_section
663 : TYPE(xas_control_type), POINTER :: xas_control
664 :
665 4 : NULLIFY (dft_control, scf_control, scf_env, matrix_ks, mos, para_env, xas_control)
666 4 : NULLIFY (dft_section, xas_section, scf_section, all_vectors, mo_coeff, all_evals, eigenvalues)
667 4 : NULLIFY (logger)
668 8 : logger => cp_get_default_logger()
669 4 : output_unit = cp_logger_get_default_io_unit(logger)
670 :
671 4 : CALL timeset(routineN, handle)
672 :
673 4 : CPASSERT(ASSOCIATED(xas_env))
674 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
675 4 : matrix_ks=matrix_ks, mos=mos, para_env=para_env)
676 :
677 4 : xas_control => dft_control%xas_control
678 4 : CALL get_xas_env(xas_env, scf_control=scf_control, spin_channel=my_spin)
679 :
680 4 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
681 4 : xas_section => section_vals_get_subs_vals(dft_section, "XAS")
682 4 : scf_section => section_vals_get_subs_vals(xas_section, "SCF")
683 :
684 4 : IF (xas_env%homo_occ == 0) THEN
685 :
686 2 : NULLIFY (scf_env)
687 2 : CALL get_xas_env(xas_env, scf_env=scf_env)
688 2 : IF (.NOT. ASSOCIATED(scf_env)) THEN
689 2 : CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
690 2 : CALL set_xas_env(xas_env, scf_env=scf_env)
691 : ELSE
692 0 : CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
693 : END IF
694 :
695 : CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
696 2 : converged=converged, should_stop=should_stop, total_scf_steps=tsteps)
697 2 : CALL scf_env_cleanup(scf_env)
698 :
699 : END IF
700 :
701 : ! The eigenstate of the KS Hamiltonian are nedeed
702 4 : NULLIFY (mo_coeff, eigenvalues)
703 4 : IF (scf_control%use_ot) THEN
704 0 : IF (output_unit > 0) THEN
705 : WRITE (UNIT=output_unit, FMT="(/,T10,A,/)") &
706 0 : "Get eigenstates and eigenvalues from ground state MOs"
707 : END IF
708 0 : DO ispin = 1, SIZE(mos)
709 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
710 0 : eigenvalues=eigenvalues)
711 : CALL calculate_subspace_eigenvalues(mo_coeff, &
712 : matrix_ks(ispin)%matrix, eigenvalues, &
713 0 : do_rotation=.TRUE.)
714 : END DO
715 : END IF
716 : CALL get_xas_env(xas_env=xas_env, all_vectors=all_vectors, &
717 4 : all_evals=all_evals, nvirtual2=nvirtual2)
718 4 : CALL get_mo_set(mos(my_spin), eigenvalues=eigenvalues, mo_coeff=mo_coeff, nmo=nmo)
719 :
720 : CALL cp_fm_to_fm(mo_coeff, all_vectors, ncol=nmo, &
721 4 : source_start=1, target_start=1)
722 36 : DO istate = 1, nmo
723 36 : all_evals(istate) = eigenvalues(istate)
724 : END DO
725 :
726 4 : IF (nvirtual2 /= 0) THEN
727 4 : IF (output_unit > 0) THEN
728 : WRITE (UNIT=output_unit, FMT="(/,T10,A,/)") &
729 2 : "WARNING: for this XES calculation additional unoccupied MOs are not needed"
730 : END IF
731 4 : nvirtual2 = 0
732 4 : nvirtual = nmo
733 4 : CALL set_xas_env(xas_env=xas_env, nvirtual=nvirtual, nvirtual2=nvirtual2)
734 : END IF
735 :
736 4 : CALL timestop(handle)
737 :
738 4 : END SUBROUTINE xes_scf_once
739 :
740 : END MODULE xas_tp_scf
|
