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