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