(git:9cbee8b)
Loading...
Searching...
No Matches
xas_methods.F
Go to the documentation of this file.
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 driver for the xas calculation and xas_scf for the tp method
10!> \par History
11!> created 05.2005
12!> replace overlap integral routine [07.2014,JGH]
13!> \author MI (05.2005)
14! **************************************************************************************************
16
17 USE ai_contraction, ONLY: block_add,&
19 USE ai_overlap, ONLY: overlap_ab
22 USE basis_set_types, ONLY: &
26 USE cell_types, ONLY: cell_type,&
27 pbc
30 USE cp_dbcsr_api, ONLY: dbcsr_convert_offsets_to_sizes,&
35 dbcsr_set,&
37 dbcsr_type_antisymmetric
47 USE cp_fm_types, ONLY: cp_fm_create,&
59 USE cp_output_handling, ONLY: cp_p_file,&
63 USE input_constants, ONLY: &
72 USE kinds, ONLY: default_string_length,&
73 dp
76 USE orbital_pointers, ONLY: ncoset
80 USE periodic_table, ONLY: ptable
81 USE physcon, ONLY: evolt
82 USE qs_diis, ONLY: qs_diis_b_clear,&
87 USE qs_kind_types, ONLY: get_qs_kind,&
89 USE qs_loc_main, ONLY: qs_loc_driver
98 USE qs_matrix_pools, ONLY: mpools_get,&
102 USE qs_mo_types, ONLY: get_mo_set,&
106 USE qs_operators_ao, ONLY: p_xyz_ao,&
109 USE qs_scf, ONLY: scf_env_cleanup
111 USE qs_scf_types, ONLY: qs_scf_env_type,&
116 USE xas_control, ONLY: read_xas_control,&
120 USE xas_env_types, ONLY: get_xas_env,&
126 USE xas_tp_scf, ONLY: xas_do_tp_scf,&
128#include "./base/base_uses.f90"
129
130 IMPLICIT NONE
131 PRIVATE
132
133 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_methods'
134
135! *** Public subroutines ***
136
137 PUBLIC :: xas, calc_stogto_overlap
138
139CONTAINS
140
141! **************************************************************************************************
142!> \brief Driver for xas calculations
143!> The initial mos are prepared
144!> A loop on the atoms to be excited is started
145!> For each atom the state to be excited is identified
146!> An scf optimization using the TP scheme or TD-DFT is used
147!> to evaluate the spectral energies and oscillator strengths
148!> \param qs_env the qs_env, the xas_env lives in
149!> \param dft_control ...
150!> \par History
151!> 05.2005 created [MI]
152!> \author MI
153!> \note
154!> the iteration counter is not finalized yet
155!> only the transition potential approach is active
156!> the localization can be switched off, otherwise
157!> it uses by default the berry phase approach
158!> The number of states to be localized is xas_control%nexc_search
159!> In general only the core states are needed
160! **************************************************************************************************
161 SUBROUTINE xas(qs_env, dft_control)
162
163 TYPE(qs_environment_type), POINTER :: qs_env
164 TYPE(dft_control_type), POINTER :: dft_control
165
166 CHARACTER(LEN=*), PARAMETER :: routinen = 'xas'
167
168 INTEGER :: handle, homo, i, iat, iatom, ispin, istate, my_homo(2), my_nelectron(2), my_spin, &
169 nao, nexc_atoms, nexc_search, nmo, nspins, output_unit, state_to_be_excited
170 INTEGER, DIMENSION(2) :: added_mos
171 INTEGER, DIMENSION(:), POINTER :: nexc_states
172 INTEGER, DIMENSION(:, :), POINTER :: state_of_atom
173 LOGICAL :: ch_method_flags, converged, my_uocc(2), &
174 should_stop, skip_scf, &
175 transition_potential
176 REAL(dp) :: maxocc, occ_estate, tmp, xas_nelectron
177 REAL(dp), DIMENSION(:), POINTER :: eigenvalues
178 REAL(dp), DIMENSION(:, :), POINTER :: vecbuffer
179 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
180 TYPE(cell_type), POINTER :: cell
181 TYPE(cp_fm_type), DIMENSION(:), POINTER :: groundstate_coeff
182 TYPE(cp_fm_type), POINTER :: all_vectors, excvec_coeff, mo_coeff
183 TYPE(cp_logger_type), POINTER :: logger
184 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, op_sm, ostrength_sm
185 TYPE(dbcsr_type), POINTER :: mo_coeff_b
186 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
187 TYPE(mp_para_env_type), POINTER :: para_env
188 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
189 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
190 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
191 TYPE(qs_scf_env_type), POINTER :: scf_env
192 TYPE(scf_control_type), POINTER :: scf_control
193 TYPE(section_vals_type), POINTER :: dft_section, loc_section, &
194 print_loc_section, scf_section, &
195 xas_section
196 TYPE(xas_control_type), POINTER :: xas_control
197 TYPE(xas_environment_type), POINTER :: xas_env
198
199 CALL timeset(routinen, handle)
200
201 transition_potential = .false.
202 skip_scf = .false.
203 converged = .true.
204 should_stop = .false.
205 ch_method_flags = .false.
206
207 NULLIFY (logger)
208 logger => cp_get_default_logger()
209 output_unit = cp_logger_get_default_io_unit(logger)
210
211 NULLIFY (xas_env, groundstate_coeff, ostrength_sm, op_sm)
212 NULLIFY (excvec_coeff, qs_loc_env, cell, scf_env)
213 NULLIFY (matrix_ks)
214 NULLIFY (all_vectors, state_of_atom, nexc_states, xas_control)
215 NULLIFY (vecbuffer, op_sm, mo_coeff_b)
216 NULLIFY (dft_section, xas_section, scf_section, loc_section, print_loc_section)
217 dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
218 xas_section => section_vals_get_subs_vals(dft_section, "XAS")
219 scf_section => section_vals_get_subs_vals(xas_section, "SCF")
220 loc_section => section_vals_get_subs_vals(xas_section, "LOCALIZE")
221 print_loc_section => section_vals_get_subs_vals(loc_section, "PRINT")
222
223 output_unit = cp_print_key_unit_nr(logger, xas_section, "PRINT%PROGRAM_RUN_INFO", &
224 extension=".Log")
225 IF (output_unit > 0) THEN
226 WRITE (unit=output_unit, fmt="(/,T3,A,/,T25,A,/,T3,A,/)") &
227 repeat("=", 77), &
228 "START CORE LEVEL SPECTROSCOPY CALCULATION", &
229 repeat("=", 77)
230 END IF
231
232! Create the xas environment
233 CALL get_qs_env(qs_env, xas_env=xas_env)
234 IF (.NOT. ASSOCIATED(xas_env)) THEN
235 IF (output_unit > 0) THEN
236 WRITE (unit=output_unit, fmt="(/,T5,A)") &
237 "Create and initialize the xas environment"
238 END IF
239 ALLOCATE (xas_env)
240 CALL xas_env_create(xas_env)
241 CALL xas_env_init(xas_env, qs_env, dft_section, logger)
242 xas_control => dft_control%xas_control
243 CALL set_qs_env(qs_env, xas_env=xas_env)
244 END IF
245
246! Initialize the type of calculation
247 NULLIFY (atomic_kind_set, qs_kind_set, scf_control, mos, para_env, particle_set)
248 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
249 cell=cell, scf_control=scf_control, &
250 matrix_ks=matrix_ks, mos=mos, para_env=para_env, &
251 particle_set=particle_set)
252
253! The eigenstate of the KS Hamiltonian are nedeed
254 NULLIFY (mo_coeff, eigenvalues)
255 IF (scf_control%use_ot) THEN
256 IF (output_unit > 0) THEN
257 WRITE (unit=output_unit, fmt="(/,T10,A,/)") &
258 "Get eigenstates and eigenvalues from ground state MOs"
259 END IF
260 DO ispin = 1, dft_control%nspins
261 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo, &
262 eigenvalues=eigenvalues, homo=homo)
263 CALL calculate_subspace_eigenvalues(mo_coeff, &
264 matrix_ks(ispin)%matrix, eigenvalues, &
265 do_rotation=.true.)
266 END DO
267 END IF
268! In xas SCF we need to use the same number of MOS as for GS
269 added_mos = scf_control%added_mos
270 NULLIFY (scf_control)
271! Consider to use get function for this
272 CALL get_xas_env(xas_env, scf_control=scf_control)
273 scf_control%added_mos = added_mos
274
275! Set initial occupation numbers, and store the original ones
276 my_homo = 0
277 my_nelectron = 0
278 DO ispin = 1, dft_control%nspins
279 CALL get_mo_set(mos(ispin), nelectron=my_nelectron(ispin), maxocc=maxocc, &
280 homo=my_homo(ispin), uniform_occupation=my_uocc(ispin))
281 END DO
282
283 nspins = dft_control%nspins
284! at the moment the only implemented method for XAS and XES calculations
285 transition_potential = .true. !(xas_control%xas_method==xas_tp_hh).OR.&
286 ! (xas_control%xas_method==xas_tp_fh).OR.&
287 ! (xas_control%xas_method==xas_tp_xhh).OR.&
288 ! (xas_control%xas_method==xas_tp_xfh).OR.&
289 ! (xas_control%xas_method==xas_dscf)
290 IF (nspins == 1 .AND. transition_potential) THEN
291 cpabort("XAS with TP method requires LSD calculations")
292 END IF
293
294 CALL get_xas_env(xas_env=xas_env, &
295 all_vectors=all_vectors, &
296 groundstate_coeff=groundstate_coeff, excvec_coeff=excvec_coeff, &
297 nexc_atoms=nexc_atoms, &
298 spin_channel=my_spin)
299
300! Set of states among which there is the state to be excited
301 CALL get_mo_set(mos(my_spin), nao=nao, homo=homo)
302 IF (xas_control%nexc_search < 0) xas_control%nexc_search = homo
303 nexc_search = xas_control%nexc_search
304
305 CALL set_xas_env(xas_env=xas_env, nexc_search=nexc_search)
306
307 !Define the qs_loc_env : to find centers, spread and possibly localize them
308 CALL get_xas_env(xas_env=xas_env, qs_loc_env=qs_loc_env)
309 IF (qs_loc_env%do_localize) THEN
310 IF (output_unit > 0) THEN
311 WRITE (unit=output_unit, fmt="(/,T2,A34,I3,A36/)") &
312 "Localize a sub-set of MOs of spin ", my_spin, ","// &
313 " to better identify the core states"
314 IF ( &
315 qs_loc_env%localized_wfn_control%set_of_states == state_loc_range) THEN
316 WRITE (unit=output_unit, fmt="( A , I7, A, I7)") " The sub-set contains states from ", &
317 qs_loc_env%localized_wfn_control%lu_bound_states(1, my_spin), " to ", &
318 qs_loc_env%localized_wfn_control%lu_bound_states(2, my_spin)
319 ELSE IF (qs_loc_env%localized_wfn_control%set_of_states == state_loc_list) THEN
320 WRITE (unit=output_unit, fmt="( A )") " The sub-set contains states given in the input list"
321 END IF
322
323 END IF
324 CALL qs_loc_driver(qs_env, qs_loc_env, print_loc_section, myspin=my_spin)
325 END IF
326
327 cpassert(ASSOCIATED(groundstate_coeff))
328 DO ispin = 1, nspins
329 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, nmo=nmo)
330 CALL cp_fm_to_fm(mo_coeff, groundstate_coeff(ispin), nmo, 1, 1)
331 IF (ASSOCIATED(mo_coeff_b)) THEN
332
333 END IF
334 END DO
335
336! SCF for only XES using occupied core and empty homo (only one SCF)
337! Probably better not to do the localization in this case, but only single out the
338! core orbital for the specific atom for which the spectrum is computed
339 IF (xas_control%xas_method == xes_tp_val .AND. &
340 xas_control%xes_core_occupation == 1.0_dp) THEN
341 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(/,/,T10,A)') &
342 "START Core Level Spectroscopy Calculation for the Emission Spectrum"
343 IF (xas_control%xes_homo_occupation == 1) THEN
344 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(T10,A,/,A)') &
345 "The core state is fully occupied and XES from ground state calculation.", &
346 " No SCF is needed, MOS already available"
347 ELSE IF (xas_control%xes_homo_occupation == 0) THEN
348 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(T10,A,/,A)') &
349 "The core state is fully occupied and the homo is empty", &
350 " (final state of the core hole decay). Only one SCF is needed (not one per atom)"
351 END IF
352 skip_scf = .true.
353
354 CALL set_xas_env(xas_env=xas_env, xas_estate=-1, homo_occ=xas_control%xes_homo_occupation)
355 CALL xes_scf_once(qs_env, xas_env, converged, should_stop)
356
357 IF (converged .AND. .NOT. should_stop .AND. xas_control%xes_homo_occupation == 0) THEN
358 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(/,T10,A,I6)') &
359 "SCF with empty homo converged "
360 ELSE IF (.NOT. converged .OR. should_stop) THEN
361 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(/,T10,A,I6)') &
362 "SCF with empty homo NOT converged"
363 ! Release what has to be released
364 IF (ASSOCIATED(vecbuffer)) THEN
365 DEALLOCATE (vecbuffer)
366 DEALLOCATE (op_sm)
367 END IF
368
369 DO ispin = 1, dft_control%nspins
370 CALL set_mo_set(mos(ispin), homo=my_homo(ispin), &
371 uniform_occupation=my_uocc(ispin), nelectron=my_nelectron(ispin))
372 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
373 CALL cp_fm_to_fm(groundstate_coeff(ispin), mos(ispin)%mo_coeff, nmo, 1, 1)
374 END DO
375
376 IF (output_unit > 0) THEN
377 WRITE (unit=output_unit, fmt="(/,T3,A,/,T25,A,/,T3,A,/)") &
378 repeat("=", 77), &
379 "END CORE LEVEL SPECTROSCOPY CALCULATION", &
380 repeat("=", 77)
381 END IF
382
383 CALL xas_env_release(qs_env%xas_env)
384 DEALLOCATE (qs_env%xas_env)
385 NULLIFY (qs_env%xas_env)
386
387 CALL cp_print_key_finished_output(output_unit, logger, xas_section, &
388 "PRINT%PROGRAM_RUN_INFO")
389 CALL timestop(handle)
390 RETURN
391 END IF
392 END IF
393
394 ! Assign the character of the selected core states
395 ! through the overlap with atomic-like states
396 CALL cls_assign_core_states(xas_control, xas_env, qs_loc_env%localized_wfn_control, &
397 qs_env)
398 CALL get_xas_env(xas_env=xas_env, &
399 state_of_atom=state_of_atom, nexc_states=nexc_states)
400
401 IF (skip_scf) THEN
402 CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff)
403 CALL cp_fm_to_fm(mo_coeff, all_vectors, ncol=nexc_search, &
404 source_start=1, target_start=1)
405 END IF
406
407 ALLOCATE (vecbuffer(1, nao))
408 ALLOCATE (op_sm(3))
409
410 ! copy the coefficients of the mos in a temporary fm with the right structure
411 IF (transition_potential) THEN
412 ! Calculate the operator
413 CALL get_xas_env(xas_env=xas_env, ostrength_sm=ostrength_sm)
414 DO i = 1, 3
415 NULLIFY (op_sm(i)%matrix)
416 op_sm(i)%matrix => ostrength_sm(i)%matrix
417 END DO
418 IF (xas_control%dipole_form == xas_dip_vel) THEN
419 CALL p_xyz_ao(op_sm, qs_env)
420 END IF
421 END IF
422
423 ! DO SCF if required
424 DO iat = 1, nexc_atoms
425 iatom = xas_env%exc_atoms(iat)
426 DO istate = 1, nexc_states(iat)
427 ! determine which state has to be excited in the global list
428 state_to_be_excited = state_of_atom(iat, istate)
429
430 ! Take the state_to_be_excited vector from the full set and copy into excvec_coeff
431 CALL get_mo_set(mos(my_spin), nmo=nmo)
432 CALL get_xas_env(xas_env, occ_estate=occ_estate, xas_nelectron=xas_nelectron)
433 tmp = xas_nelectron + 1.0_dp - occ_estate
434 IF (nmo < tmp) THEN
435 cpabort("CLS: the required method needs added_mos to the ground state")
436 END IF
437 ! If the restart file for this atom exists, the mos and the
438 ! occupation numbers are overwritten
439 ! It is necessary that the restart is for the same xas method
440 ! otherwise the number of electrons and the occupation numbers
441 ! may not be consistent
442 IF (xas_control%xas_restart) THEN
443 CALL xas_read_restart(xas_env, xas_section, qs_env, xas_control%xas_method, iatom, &
444 state_to_be_excited, istate)
445 END IF
446 CALL set_xas_env(xas_env=xas_env, xas_estate=state_to_be_excited)
447 CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff)
448 cpassert(ASSOCIATED(excvec_coeff))
449 CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, state_to_be_excited, &
450 nao, 1, transpose=.true.)
451 CALL cp_fm_set_submatrix(excvec_coeff, vecbuffer, 1, 1, &
452 nao, 1, transpose=.true.)
453
454 IF (transition_potential) THEN
455
456 IF (.NOT. skip_scf) THEN
457 IF (output_unit > 0) THEN
458 WRITE (unit=output_unit, fmt='(/,T5,A)') repeat("-", 75)
459 IF (xas_control%xas_method == xas_dscf) THEN
460 WRITE (unit=output_unit, fmt='(/,/,T10,A,I6)') &
461 "START DeltaSCF for the first excited state from the core state of ATOM ", iatom
462 ELSE
463 WRITE (unit=output_unit, fmt='(/,T10,A,I6)') &
464 "Start Core Level Spectroscopy Calculation with TP approach for ATOM ", iatom
465 WRITE (unit=output_unit, fmt='(/,T10,A,I6,T34,A,T54,I6)') &
466 "Excited state", istate, "out of", nexc_states(iat)
467 WRITE (unit=output_unit, fmt='(T10,A,T50,f10.4)') "Occupation of the core orbital", &
468 occ_estate
469 WRITE (unit=output_unit, fmt='(T10,A28,I3, T50,F10.4)') "Number of electrons in Spin ", &
470 my_spin, xas_nelectron
471 END IF
472 END IF
473
474 CALL get_xas_env(xas_env=xas_env, scf_env=scf_env)
475 IF (.NOT. ASSOCIATED(scf_env)) THEN
476 CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
477 ! Moved here from qs_scf_env_initialize to be able to have more scf_env
478 CALL set_xas_env(xas_env, scf_env=scf_env)
479 ELSE
480 CALL qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
481 END IF
482
483 DO ispin = 1, SIZE(mos)
484 IF (ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN !fm->dbcsr
485 CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
486 mos(ispin)%mo_coeff_b) !fm->dbcsr
487 END IF !fm->dbcsr
488 END DO !fm->dbcsr
489
490 IF (.NOT. scf_env%skip_diis) THEN
491 IF (.NOT. ASSOCIATED(scf_env%scf_diis_buffer)) THEN
492 ALLOCATE (scf_env%scf_diis_buffer)
493 CALL qs_diis_b_create(scf_env%scf_diis_buffer, nbuffer=scf_control%max_diis)
494 END IF
495 CALL qs_diis_b_clear(scf_env%scf_diis_buffer)
496 END IF
497
498 CALL xas_do_tp_scf(dft_control, xas_env, iatom, istate, scf_env, qs_env, &
499 xas_section, scf_section, converged, should_stop)
500
501 CALL external_control(should_stop, "CLS", target_time=qs_env%target_time, &
502 start_time=qs_env%start_time)
503 IF (should_stop) THEN
504 CALL scf_env_cleanup(scf_env)
505 EXIT
506 END IF
507
508 END IF
509 ! SCF DONE
510
511 ! Write last wavefunction to screen
512 IF (SIZE(mos) > 1) THEN
513 CALL write_mo_set_to_output_unit(mos(1), qs_kind_set, particle_set, dft_section, &
514 4, 0, final_mos=.false., spin="XAS ALPHA")
515 CALL write_mo_set_to_output_unit(mos(2), qs_kind_set, particle_set, dft_section, &
516 4, 0, final_mos=.false., spin="XAS BETA")
517 ELSE
518 CALL write_mo_set_to_output_unit(mos(1), qs_kind_set, particle_set, dft_section, &
519 4, 0, final_mos=.false., spin="XAS")
520 END IF
521
522 ELSE
523 ! Core level spectroscopy by TDDFT is not yet implemented
524 ! the states defined by the rotation are the ground state orbitals
525 ! the initial state from which I excite should be localized
526 ! I take the excitations from lumo to nmo
527 END IF
528
529 IF (converged) THEN
530 CALL cls_calculate_spectrum(xas_control, xas_env, qs_env, xas_section, &
531 iatom, istate)
532 ELSE
533 IF (output_unit > 0) WRITE (unit=output_unit, fmt='(/,/,T10,A,I6)') &
534 "SCF with core hole NOT converged for ATOM ", iatom
535 END IF
536
537 IF (.NOT. skip_scf) THEN
538 ! Reset the initial core orbitals.
539 ! The valence orbitals are taken from the last SCF,
540 ! it should be a better initial guess
541 CALL get_qs_env(qs_env, mos=mos)
542 DO ispin = 1, nspins
543 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
544 CALL cp_fm_to_fm(groundstate_coeff(ispin), mos(ispin)%mo_coeff, nmo, 1, 1)
545 END DO
546 IF (iat == nexc_atoms) THEN
547 CALL scf_env_cleanup(scf_env)
548 CALL scf_env_release(xas_env%scf_env)
549 DEALLOCATE (xas_env%scf_env)
550 END IF
551 END IF
552
553 END DO ! istate
554 END DO ! iat = 1,nexc_atoms
555
556 ! END of Calculation
557
558 ! Release what has to be released
559 IF (ASSOCIATED(vecbuffer)) THEN
560 DEALLOCATE (vecbuffer)
561 DEALLOCATE (op_sm)
562 END IF
563
564 DO ispin = 1, dft_control%nspins
565 CALL set_mo_set(mos(ispin), homo=my_homo(ispin), &
566 uniform_occupation=my_uocc(ispin), nelectron=my_nelectron(ispin))
567 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
568 CALL cp_fm_to_fm(groundstate_coeff(ispin), mos(ispin)%mo_coeff, nmo, 1, 1)
569 END DO
570
571 IF (output_unit > 0) THEN
572 WRITE (unit=output_unit, fmt="(/,T3,A,/,T25,A,/,T3,A,/)") &
573 repeat("=", 77), &
574 "END CORE LEVEL SPECTROSCOPY CALCULATION", &
575 repeat("=", 77)
576 END IF
577
578 CALL xas_env_release(qs_env%xas_env)
579 DEALLOCATE (qs_env%xas_env)
580 NULLIFY (qs_env%xas_env)
581
582 CALL cp_print_key_finished_output(output_unit, logger, xas_section, &
583 "PRINT%PROGRAM_RUN_INFO")
584 CALL timestop(handle)
585
586 END SUBROUTINE xas
587
588! **************************************************************************************************
589!> \brief allocate and initialize the structure needed for the xas calculation
590!> \param xas_env the environment for XAS calculations
591!> \param qs_env the qs_env, the xas_env lives in
592!> \param dft_section ...
593!> \param logger ...
594!> \par History
595!> 05.2005 created [MI]
596!> \author MI
597! **************************************************************************************************
598 SUBROUTINE xas_env_init(xas_env, qs_env, dft_section, logger)
599
600 TYPE(xas_environment_type), POINTER :: xas_env
601 TYPE(qs_environment_type), POINTER :: qs_env
602 TYPE(section_vals_type), POINTER :: dft_section
603 TYPE(cp_logger_type), POINTER :: logger
604
605 CHARACTER(LEN=default_string_length) :: name_sto
606 INTEGER :: homo, i, iat, iatom, ik, ikind, ispin, j, l, lfomo, my_spin, n_mo(2), n_rep, nao, &
607 natom, ncubes, nelectron, nexc_atoms, nexc_search, nj, nk, nkind, nmo, nmoloc(2), &
608 nsgf_gto, nsgf_sto, nspins, nvirtual, nvirtual2
609 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_type_tmp, kind_z_tmp, &
610 last_sgf
611 INTEGER, DIMENSION(4, 7) :: ne
612 INTEGER, DIMENSION(:), POINTER :: bounds, list, lq, nq, row_blk_sizes
613 LOGICAL :: ihavethis
614 REAL(dp) :: nele, occ_estate, occ_homo, &
615 occ_homo_plus, zatom
616 REAL(dp), DIMENSION(:), POINTER :: sto_zet
617 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
618 TYPE(atomic_kind_type), POINTER :: atomic_kind
619 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
620 TYPE(cp_fm_type), POINTER :: mo_coeff
621 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
622 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
623 TYPE(dft_control_type), POINTER :: dft_control
624 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
625 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
626 TYPE(mp_para_env_type), POINTER :: para_env
627 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
628 POINTER :: sab_orb
629 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
630 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
631 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
632 TYPE(qs_matrix_pools_type), POINTER :: mpools
633 TYPE(scf_control_type), POINTER :: scf_control
634 TYPE(section_vals_type), POINTER :: loc_section, xas_section
635 TYPE(sto_basis_set_type), POINTER :: sto_basis_set
636 TYPE(xas_control_type), POINTER :: xas_control
637
638 n_mo(1:2) = 0
639 cpassert(ASSOCIATED(xas_env))
640
641 NULLIFY (atomic_kind_set, qs_kind_set, dft_control, scf_control, matrix_s, mos, mpools)
642 NULLIFY (para_env, particle_set, xas_control)
643 NULLIFY (qs_loc_env)
644 NULLIFY (sab_orb)
645 CALL get_qs_env(qs_env=qs_env, &
646 atomic_kind_set=atomic_kind_set, &
647 qs_kind_set=qs_kind_set, &
648 dft_control=dft_control, &
649 mpools=mpools, &
650 matrix_s=matrix_s, mos=mos, &
651 para_env=para_env, particle_set=particle_set, &
652 sab_orb=sab_orb, &
653 dbcsr_dist=dbcsr_dist)
654
655 xas_section => section_vals_get_subs_vals(dft_section, "XAS")
656 ALLOCATE (dft_control%xas_control)
657 CALL xas_control_create(dft_control%xas_control)
658 CALL read_xas_control(dft_control%xas_control, xas_section)
659 CALL write_xas_control(dft_control%xas_control, dft_section)
660 xas_control => dft_control%xas_control
661 ALLOCATE (scf_control)
662 CALL scf_c_create(scf_control)
663 CALL scf_c_read_parameters(scf_control, xas_section)
664 CALL set_xas_env(xas_env, scf_control=scf_control)
665
666 my_spin = xas_control%spin_channel
667 nexc_search = xas_control%nexc_search
668 IF (nexc_search < 0) THEN
669 ! ground state occupation
670 CALL get_mo_set(mos(my_spin), nmo=nmo, lfomo=lfomo)
671 nexc_search = lfomo - 1
672 END IF
673 nexc_atoms = xas_control%nexc_atoms
674 ALLOCATE (xas_env%exc_atoms(nexc_atoms))
675 xas_env%exc_atoms = xas_control%exc_atoms
676 CALL set_xas_env(xas_env=xas_env, nexc_search=nexc_search, &
677 nexc_atoms=nexc_atoms, spin_channel=my_spin)
678
679 CALL mpools_get(mpools, ao_mo_fm_pools=xas_env%ao_mo_fm_pools)
680
681 NULLIFY (mo_coeff)
682 CALL get_mo_set(mos(my_spin), nao=nao, homo=homo, nmo=nmo, mo_coeff=mo_coeff, nelectron=nelectron)
683
684 nvirtual2 = 0
685 IF (xas_control%added_mos > 0) THEN
686 nvirtual2 = min(xas_control%added_mos, nao - nmo)
687 xas_env%unoccupied_eps = xas_control%eps_added
688 xas_env%unoccupied_max_iter = xas_control%max_iter_added
689 END IF
690 nvirtual = nmo + nvirtual2
691
692 n_mo(1:2) = nmo
693
694 ALLOCATE (xas_env%centers_wfn(3, nexc_search))
695 ALLOCATE (xas_env%atom_of_state(nexc_search))
696 ALLOCATE (xas_env%type_of_state(nexc_search))
697 ALLOCATE (xas_env%state_of_atom(nexc_atoms, nexc_search))
698 ALLOCATE (xas_env%nexc_states(nexc_atoms))
699 ALLOCATE (xas_env%mykind_of_atom(nexc_atoms))
700 nkind = SIZE(atomic_kind_set, 1)
701 ALLOCATE (xas_env%mykind_of_kind(nkind))
702 xas_env%mykind_of_kind = 0
703
704 ! create a new matrix structure nao x 1
705 NULLIFY (tmp_fm_struct)
706 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
707 ncol_global=1, para_env=para_env, context=mo_coeff%matrix_struct%context)
708 ALLOCATE (xas_env%excvec_coeff)
709 CALL cp_fm_create(xas_env%excvec_coeff, tmp_fm_struct)
710 CALL cp_fm_struct_release(tmp_fm_struct)
711
712 NULLIFY (tmp_fm_struct)
713 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=1, &
714 ncol_global=nexc_search, para_env=para_env, &
715 context=mo_coeff%matrix_struct%context)
716 ALLOCATE (xas_env%excvec_overlap)
717 CALL cp_fm_create(xas_env%excvec_overlap, tmp_fm_struct)
718 CALL cp_fm_struct_release(tmp_fm_struct)
719
720 nspins = SIZE(mos, 1)
721
722 ! initialize operators for the calculation of the oscillator strengths
723 IF (xas_control%xas_method == xas_tp_hh) THEN
724 occ_estate = 0.5_dp
725 nele = real(nelectron, dp) - 0.5_dp
726 occ_homo = 1.0_dp
727 occ_homo_plus = 0._dp
728 ELSE IF (xas_control%xas_method == xas_tp_xhh) THEN
729 occ_estate = 0.5_dp
730 nele = real(nelectron, dp)
731 occ_homo = 1.0_dp
732 occ_homo_plus = 0.5_dp
733 ELSE IF (xas_control%xas_method == xas_tp_fh) THEN
734 occ_estate = 0.0_dp
735 nele = real(nelectron, dp) - 1.0_dp
736 occ_homo = 1.0_dp
737 occ_homo_plus = 0._dp
738 ELSE IF (xas_control%xas_method == xas_tp_xfh) THEN
739 occ_estate = 0.0_dp
740 nele = real(nelectron, dp)
741 occ_homo = 1.0_dp
742 occ_homo_plus = 1._dp
743 ELSE IF (xas_control%xas_method == xes_tp_val) THEN
744 occ_estate = xas_control%xes_core_occupation
745 nele = real(nelectron, dp) - xas_control%xes_core_occupation
746 occ_homo = xas_control%xes_homo_occupation
747 ELSE IF (xas_control%xas_method == xas_dscf) THEN
748 occ_estate = 0.0_dp
749 nele = real(nelectron, dp)
750 occ_homo = 1.0_dp
751 occ_homo_plus = 1._dp
752 ELSE IF (xas_control%xas_method == xas_tp_flex) THEN
753 nele = real(xas_control%nel_tot, dp)
754 occ_estate = real(xas_control%xas_core_occupation, dp)
755 IF (nele < 0.0_dp) nele = real(nelectron, dp) - (1.0_dp - occ_estate)
756 occ_homo = 1.0_dp
757 END IF
758 CALL set_xas_env(xas_env=xas_env, occ_estate=occ_estate, xas_nelectron=nele, &
759 nvirtual2=nvirtual2, nvirtual=nvirtual, homo_occ=occ_homo)
760
761 ! Initialize the list of orbitals for cube files printing
762 IF (btest(cp_print_key_should_output(logger%iter_info, xas_section, &
763 "PRINT%CLS_FUNCTION_CUBES"), cp_p_file)) THEN
764 NULLIFY (bounds, list)
765 CALL section_vals_val_get(xas_section, &
766 "PRINT%CLS_FUNCTION_CUBES%CUBES_LU_BOUNDS", &
767 i_vals=bounds)
768 ncubes = bounds(2) - bounds(1) + 1
769 IF (ncubes > 0) THEN
770 ALLOCATE (xas_control%list_cubes(ncubes))
771
772 DO ik = 1, ncubes
773 xas_control%list_cubes(ik) = bounds(1) + (ik - 1)
774 END DO
775 END IF
776
777 IF (.NOT. ASSOCIATED(xas_control%list_cubes)) THEN
778 CALL section_vals_val_get(xas_section, &
779 "PRINT%CLS_FUNCTION_CUBES%CUBES_LIST", &
780 n_rep_val=n_rep)
781 ncubes = 0
782 DO ik = 1, n_rep
783 NULLIFY (list)
784 CALL section_vals_val_get(xas_section, &
785 "PRINT%CLS_FUNCTION_CUBES%CUBES_LIST", &
786 i_rep_val=ik, i_vals=list)
787 IF (ASSOCIATED(list)) THEN
788 CALL reallocate(xas_control%list_cubes, 1, ncubes + SIZE(list))
789 DO i = 1, SIZE(list)
790 xas_control%list_cubes(i + ncubes) = list(i)
791 END DO
792 ncubes = ncubes + SIZE(list)
793 END IF
794 END DO ! ik
795 END IF
796
797 IF (.NOT. ASSOCIATED(xas_control%list_cubes)) THEN
798 ncubes = max(10, xas_control%added_mos/10)
799 ncubes = min(ncubes, xas_control%added_mos)
800 ALLOCATE (xas_control%list_cubes(ncubes))
801 DO ik = 1, ncubes
802 xas_control%list_cubes(ik) = homo + ik
803 END DO
804 END IF
805 ELSE
806 NULLIFY (xas_control%list_cubes)
807 END IF
808
809 NULLIFY (tmp_fm_struct)
810 ALLOCATE (xas_env%groundstate_coeff(nspins))
811 DO ispin = 1, nspins
812 CALL get_mo_set(mos(ispin), nao=nao, nmo=nmo)
813 CALL fm_pool_create_fm(xas_env%ao_mo_fm_pools(ispin)%pool, &
814 xas_env%groundstate_coeff(ispin), &
815 name="xas_env%mo0"//trim(adjustl(cp_to_string(ispin))))
816 END DO ! ispin
817
818 NULLIFY (tmp_fm_struct)
819 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=1, &
820 ncol_global=nvirtual, para_env=para_env, &
821 context=mo_coeff%matrix_struct%context)
822 ALLOCATE (xas_env%dip_fm_set(2, 3))
823 DO i = 1, 3
824 DO j = 1, 2
825 CALL cp_fm_create(xas_env%dip_fm_set(j, i), tmp_fm_struct)
826 END DO
827 END DO
828 CALL cp_fm_struct_release(tmp_fm_struct)
829
830 !Array to store all the eigenstates: occupied and the required not occupied
831 IF (nvirtual2 > 0) THEN
832 ALLOCATE (xas_env%unoccupied_evals(nvirtual2))
833 NULLIFY (tmp_fm_struct)
834 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
835 ncol_global=nvirtual2, &
836 para_env=para_env, context=mo_coeff%matrix_struct%context)
837 ALLOCATE (xas_env%unoccupied_orbs)
838 CALL cp_fm_create(xas_env%unoccupied_orbs, tmp_fm_struct)
839 CALL cp_fm_struct_release(tmp_fm_struct)
840 END IF
841
842 NULLIFY (tmp_fm_struct)
843 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
844 ncol_global=nvirtual, &
845 para_env=para_env, context=mo_coeff%matrix_struct%context)
846 ALLOCATE (xas_env%all_vectors)
847 CALL cp_fm_create(xas_env%all_vectors, tmp_fm_struct)
848 CALL cp_fm_struct_release(tmp_fm_struct)
849
850 ! Array to store all the energies needed for the spectrum
851 ALLOCATE (xas_env%all_evals(nvirtual))
852
853 IF (xas_control%dipole_form == xas_dip_len) THEN
854 CALL dbcsr_allocate_matrix_set(xas_env%ostrength_sm, 3)
855 DO i = 1, 3
856 ALLOCATE (xas_env%ostrength_sm(i)%matrix)
857 CALL dbcsr_copy(xas_env%ostrength_sm(i)%matrix, matrix_s(1)%matrix, &
858 "xas_env%ostrength_sm-"//trim(adjustl(cp_to_string(i))))
859 CALL dbcsr_set(xas_env%ostrength_sm(i)%matrix, 0.0_dp)
860 END DO
861 ELSE IF (xas_control%dipole_form == xas_dip_vel) THEN
862 !
863 ! prepare for allocation
864 natom = SIZE(particle_set, 1)
865 ALLOCATE (first_sgf(natom))
866 ALLOCATE (last_sgf(natom))
867 CALL get_particle_set(particle_set, qs_kind_set, &
868 first_sgf=first_sgf, &
869 last_sgf=last_sgf)
870 ALLOCATE (row_blk_sizes(natom))
871 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
872 DEALLOCATE (first_sgf)
873 DEALLOCATE (last_sgf)
874 !
875 !
876 CALL dbcsr_allocate_matrix_set(xas_env%ostrength_sm, 3)
877 ALLOCATE (xas_env%ostrength_sm(1)%matrix)
878 CALL dbcsr_create(matrix=xas_env%ostrength_sm(1)%matrix, &
879 name="xas_env%ostrength_sm", &
880 dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
881 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
882 mutable_work=.true.)
883 CALL cp_dbcsr_alloc_block_from_nbl(xas_env%ostrength_sm(1)%matrix, sab_orb)
884 CALL dbcsr_set(xas_env%ostrength_sm(1)%matrix, 0.0_dp)
885 DO i = 2, 3
886 ALLOCATE (xas_env%ostrength_sm(i)%matrix)
887 CALL dbcsr_copy(xas_env%ostrength_sm(i)%matrix, xas_env%ostrength_sm(1)%matrix, &
888 "xas_env%ostrength_sm-"//trim(adjustl(cp_to_string(i))))
889 CALL dbcsr_set(xas_env%ostrength_sm(i)%matrix, 0.0_dp)
890 END DO
891
892 DEALLOCATE (row_blk_sizes)
893 END IF
894
895 ! Define the qs_loc_env : to find centers, spread and possibly localize them
896 IF (.NOT. (ASSOCIATED(xas_env%qs_loc_env))) THEN
897 ALLOCATE (qs_loc_env)
898 CALL qs_loc_env_create(qs_loc_env)
899 CALL set_xas_env(xas_env=xas_env, qs_loc_env=qs_loc_env)
900 loc_section => section_vals_get_subs_vals(xas_section, "LOCALIZE")
901
902 CALL qs_loc_control_init(qs_loc_env, loc_section, do_homo=.true., &
903 do_xas=.true., nloc_xas=nexc_search, spin_xas=my_spin)
904
905 IF (.NOT. qs_loc_env%do_localize) THEN
906 qs_loc_env%localized_wfn_control%localization_method = do_loc_none
907
908 ELSE
909 nmoloc = qs_loc_env%localized_wfn_control%nloc_states
910 CALL set_loc_wfn_lists(qs_loc_env%localized_wfn_control, nmoloc, n_mo, nspins, my_spin)
911 CALL set_loc_centers(qs_loc_env%localized_wfn_control, nmoloc, nspins)
912 CALL qs_loc_env_init(qs_loc_env, qs_loc_env%localized_wfn_control, &
913 qs_env, myspin=my_spin, do_localize=qs_loc_env%do_localize)
914 END IF
915 END IF
916
917 !Type of state
918 ALLOCATE (nq(1), lq(1), sto_zet(1))
919 IF (xas_control%state_type == xas_1s_type) THEN
920 nq(1) = 1
921 lq(1) = 0
922 ELSE IF (xas_control%state_type == xas_2s_type) THEN
923 nq(1) = 2
924 lq(1) = 0
925 ELSE IF (xas_control%state_type == xas_2p_type) THEN
926 nq(1) = 2
927 lq(1) = 1
928 ELSE IF (xas_control%state_type == xas_3s_type) THEN
929 nq(1) = 3
930 lq(1) = 0
931 ELSE IF (xas_control%state_type == xas_3p_type) THEN
932 nq(1) = 3
933 lq(1) = 1
934 ELSE IF (xas_control%state_type == xas_3d_type) THEN
935 nq(1) = 3
936 lq(1) = 2
937 ELSE IF (xas_control%state_type == xas_4s_type) THEN
938 nq(1) = 4
939 lq(1) = 0
940 ELSE IF (xas_control%state_type == xas_4p_type) THEN
941 nq(1) = 4
942 lq(1) = 1
943 ELSE IF (xas_control%state_type == xas_4d_type) THEN
944 nq(1) = 4
945 lq(1) = 2
946 ELSE IF (xas_control%state_type == xas_4f_type) THEN
947 nq(1) = 4
948 lq(1) = 3
949 ELSE
950 cpabort("XAS type of state not implemented")
951 END IF
952
953! Find core orbitals of right angular momentum
954 ALLOCATE (kind_type_tmp(nkind))
955 ALLOCATE (kind_z_tmp(nkind))
956 kind_type_tmp = 0
957 kind_z_tmp = 0
958 nk = 0
959 DO iat = 1, nexc_atoms
960 iatom = xas_env%exc_atoms(iat)
961 NULLIFY (atomic_kind)
962 atomic_kind => particle_set(iatom)%atomic_kind
963 CALL get_atomic_kind(atomic_kind=atomic_kind, kind_number=ikind)
964 CALL get_qs_kind(qs_kind_set(ikind), zeff=zatom)
965 ihavethis = .false.
966 DO ik = 1, nk
967 IF (ikind == kind_type_tmp(ik)) THEN
968 ihavethis = .true.
969 xas_env%mykind_of_atom(iat) = ik
970 EXIT
971 END IF
972 END DO
973 IF (.NOT. ihavethis) THEN
974 nk = nk + 1
975 kind_type_tmp(nk) = ikind
976 kind_z_tmp(nk) = int(zatom)
977 xas_env%mykind_of_atom(iat) = nk
978 xas_env%mykind_of_kind(ikind) = nk
979 END IF
980 END DO ! iat
981
982 ALLOCATE (xas_env%my_gto_basis(nk))
983 ALLOCATE (xas_env%stogto_overlap(nk))
984 DO ik = 1, nk
985 NULLIFY (xas_env%my_gto_basis(ik)%gto_basis_set, sto_basis_set)
986 ne = 0
987 DO l = 1, lq(1) + 1
988 nj = 2*(l - 1) + 1
989 DO i = l, nq(1)
990 ne(l, i) = ptable(kind_z_tmp(ik))%e_conv(l - 1) - 2*nj*(i - l)
991 ne(l, i) = max(ne(l, i), 0)
992 ne(l, i) = min(ne(l, i), 2*nj)
993 END DO
994 END DO
995
996 sto_zet(1) = srules(kind_z_tmp(ik), ne, nq(1), lq(1))
997 CALL allocate_sto_basis_set(sto_basis_set)
998 name_sto = 'xas_tmp_sto'
999 CALL set_sto_basis_set(sto_basis_set, nshell=1, nq=nq, &
1000 lq=lq, zet=sto_zet, name=name_sto)
1001 CALL create_gto_from_sto_basis(sto_basis_set, &
1002 xas_env%my_gto_basis(ik)%gto_basis_set, xas_control%ngauss)
1003 CALL deallocate_sto_basis_set(sto_basis_set)
1004 xas_env%my_gto_basis(ik)%gto_basis_set%norm_type = 2
1005 CALL init_orb_basis_set(xas_env%my_gto_basis(ik)%gto_basis_set)
1006
1007 ikind = kind_type_tmp(ik)
1008 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1009
1010 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nsgf_gto)
1011 CALL get_gto_basis_set(gto_basis_set=xas_env%my_gto_basis(ik)%gto_basis_set, nsgf=nsgf_sto)
1012 ALLOCATE (xas_env%stogto_overlap(ik)%array(nsgf_sto, nsgf_gto))
1013
1014 CALL calc_stogto_overlap(xas_env%my_gto_basis(ik)%gto_basis_set, orb_basis_set, &
1015 xas_env%stogto_overlap(ik)%array)
1016 END DO
1017
1018 DEALLOCATE (nq, lq, sto_zet)
1019 DEALLOCATE (kind_type_tmp, kind_z_tmp)
1020
1021 END SUBROUTINE xas_env_init
1022
1023! **************************************************************************************************
1024!> \brief Calculate and write the spectrum relative to the core level excitation
1025!> of a specific atom. It works for TP approach, because of the definition
1026!> of the oscillator strengths as matrix elements of the dipole operator
1027!> \param xas_control ...
1028!> \param xas_env ...
1029!> \param qs_env ...
1030!> \param xas_section ...
1031!> \param iatom index of the excited atom
1032!> \param istate ...
1033!> \par History
1034!> 03.2006 created [MI]
1035!> \author MI
1036!> \note
1037!> for the tddft calculation should be re-thought
1038! **************************************************************************************************
1039 SUBROUTINE cls_calculate_spectrum(xas_control, xas_env, qs_env, xas_section, &
1040 iatom, istate)
1041
1043 TYPE(xas_environment_type), POINTER :: xas_env
1044 TYPE(qs_environment_type), POINTER :: qs_env
1045 TYPE(section_vals_type), POINTER :: xas_section
1046 INTEGER, INTENT(IN) :: iatom, istate
1047
1048 INTEGER :: homo, i, lfomo, my_spin, nabs, nmo, &
1049 nvirtual, output_unit, xas_estate
1050 LOGICAL :: append_cube, length
1051 REAL(dp) :: rc(3)
1052 REAL(dp), DIMENSION(:), POINTER :: all_evals
1053 REAL(dp), DIMENSION(:, :), POINTER :: sp_ab, sp_em
1054 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dip_fm_set
1055 TYPE(cp_fm_type), POINTER :: all_vectors, excvec_coeff
1056 TYPE(cp_logger_type), POINTER :: logger
1057 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_sm, ostrength_sm
1058 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1059 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1060
1061 NULLIFY (logger)
1062 logger => cp_get_default_logger()
1063 output_unit = cp_logger_get_default_io_unit(logger)
1064
1065 NULLIFY (ostrength_sm, op_sm, dip_fm_set)
1066 NULLIFY (all_evals, all_vectors, excvec_coeff)
1067 NULLIFY (mos, particle_set, sp_em, sp_ab)
1068 ALLOCATE (op_sm(3))
1069
1070 CALL get_qs_env(qs_env=qs_env, &
1071 mos=mos, particle_set=particle_set)
1072
1073 CALL get_xas_env(xas_env=xas_env, all_vectors=all_vectors, xas_estate=xas_estate, &
1074 all_evals=all_evals, dip_fm_set=dip_fm_set, excvec_coeff=excvec_coeff, &
1075 ostrength_sm=ostrength_sm, nvirtual=nvirtual, spin_channel=my_spin)
1076 CALL get_mo_set(mos(my_spin), homo=homo, lfomo=lfomo, nmo=nmo)
1077
1078 nabs = nvirtual - lfomo + 1
1079 ALLOCATE (sp_em(6, homo))
1080 ALLOCATE (sp_ab(6, nabs))
1081 cpassert(ASSOCIATED(excvec_coeff))
1082
1083 IF (.NOT. xas_control%xas_method == xas_dscf) THEN
1084 ! Calculate the spectrum
1085 IF (xas_control%dipole_form == xas_dip_len) THEN
1086 rc(1:3) = particle_set(iatom)%r(1:3)
1087 DO i = 1, 3
1088 NULLIFY (op_sm(i)%matrix)
1089 op_sm(i)%matrix => ostrength_sm(i)%matrix
1090 END DO
1091 CALL rrc_xyz_ao(op_sm, qs_env, rc, order=1, minimum_image=.true.)
1092 CALL spectrum_dip_vel(dip_fm_set, op_sm, mos, excvec_coeff, &
1093 all_vectors, all_evals, &
1094 sp_em, sp_ab, xas_estate, nvirtual, my_spin)
1095 DO i = 1, SIZE(ostrength_sm, 1)
1096 CALL dbcsr_set(ostrength_sm(i)%matrix, 0.0_dp)
1097 END DO
1098 ELSE
1099 DO i = 1, 3
1100 NULLIFY (op_sm(i)%matrix)
1101 op_sm(i)%matrix => ostrength_sm(i)%matrix
1102 END DO
1103 CALL spectrum_dip_vel(dip_fm_set, op_sm, mos, excvec_coeff, &
1104 all_vectors, all_evals, &
1105 sp_em, sp_ab, xas_estate, nvirtual, my_spin)
1106 END IF
1107 END IF
1108
1109 CALL get_mo_set(mos(my_spin), lfomo=lfomo)
1110 ! write the spectrum, if the file exists it is appended
1111 IF (.NOT. xas_control%xas_method == xas_dscf) THEN
1112 length = (.NOT. xas_control%dipole_form == xas_dip_vel)
1113 CALL xas_write(sp_em, sp_ab, xas_estate, &
1114 xas_section, iatom, istate, lfomo, length=length)
1115 END IF
1116
1117 DEALLOCATE (sp_em)
1118 DEALLOCATE (sp_ab)
1119
1120 IF (btest(cp_print_key_should_output(logger%iter_info, xas_section, &
1121 "PRINT%CLS_FUNCTION_CUBES"), cp_p_file)) THEN
1122 append_cube = section_get_lval(xas_section, "PRINT%CLS_FUNCTION_CUBES%APPEND")
1123 CALL xas_print_cubes(xas_control, qs_env, xas_section, mos, all_vectors, &
1124 iatom, append_cube)
1125 END IF
1126
1127 IF (btest(cp_print_key_should_output(logger%iter_info, xas_section, &
1128 "PRINT%PDOS"), cp_p_file)) THEN
1129 CALL xas_pdos(qs_env, xas_section, mos, iatom)
1130 END IF
1131
1132 DEALLOCATE (op_sm)
1133
1134 END SUBROUTINE cls_calculate_spectrum
1135
1136! **************************************************************************************************
1137!> \brief write the spectrum for each atom in a different output file
1138!> \param sp_em ...
1139!> \param sp_ab ...
1140!> \param estate ...
1141!> \param xas_section ...
1142!> \param iatom index of the excited atom
1143!> \param state_to_be_excited ...
1144!> \param lfomo ...
1145!> \param length ...
1146!> \par History
1147!> 05.2005 created [MI]
1148!> \author MI
1149!> \note
1150!> the iteration counter is not finilized yet
1151! **************************************************************************************************
1152 SUBROUTINE xas_write(sp_em, sp_ab, estate, xas_section, iatom, state_to_be_excited, &
1153 lfomo, length)
1154
1155 REAL(dp), DIMENSION(:, :), POINTER :: sp_em, sp_ab
1156 INTEGER, INTENT(IN) :: estate
1157 TYPE(section_vals_type), POINTER :: xas_section
1158 INTEGER, INTENT(IN) :: iatom, state_to_be_excited, lfomo
1159 LOGICAL, INTENT(IN) :: length
1160
1161 CHARACTER(LEN=default_string_length) :: mittle_ab, mittle_em, my_act, my_pos
1162 INTEGER :: i, istate, out_sp_ab, out_sp_em
1163 REAL(dp) :: ene2
1164 TYPE(cp_logger_type), POINTER :: logger
1165
1166 NULLIFY (logger)
1167 logger => cp_get_default_logger()
1168
1169 my_pos = "APPEND"
1170 my_act = "WRITE"
1171
1172 mittle_em = "xes_at"//trim(adjustl(cp_to_string(iatom)))//"_st"//trim(adjustl(cp_to_string(state_to_be_excited)))
1173
1174 out_sp_em = cp_print_key_unit_nr(logger, xas_section, "PRINT%XES_SPECTRUM", &
1175 extension=".spectrum", file_position=my_pos, file_action=my_act, &
1176 file_form="FORMATTED", middle_name=trim(mittle_em))
1177
1178 IF (out_sp_em > 0) THEN
1179 WRITE (out_sp_em, '(A,I6,A,I6,A,I6)') " Emission spectrum for atom ", iatom, &
1180 ", index of excited core MO is", estate, ", # of lines ", SIZE(sp_em, 2)
1181 ene2 = 1.0_dp
1182 DO istate = estate, SIZE(sp_em, 2)
1183 IF (length) ene2 = sp_em(1, istate)*sp_em(1, istate)
1184 WRITE (out_sp_em, '(I6,5F16.8,F10.5)') istate, sp_em(1, istate)*evolt, &
1185 sp_em(2, istate)*ene2, sp_em(3, istate)*ene2, &
1186 sp_em(4, istate)*ene2, sp_em(5, istate)*ene2, sp_em(6, istate)
1187 END DO
1188 END IF
1189 CALL cp_print_key_finished_output(out_sp_em, logger, xas_section, &
1190 "PRINT%XES_SPECTRUM")
1191
1192 mittle_ab = "xas_at"//trim(adjustl(cp_to_string(iatom)))//"_st"//trim(adjustl(cp_to_string(state_to_be_excited)))
1193 out_sp_ab = cp_print_key_unit_nr(logger, xas_section, "PRINT%XAS_SPECTRUM", &
1194 extension=".spectrum", file_position=my_pos, file_action=my_act, &
1195 file_form="FORMATTED", middle_name=trim(mittle_ab))
1196
1197 IF (out_sp_ab > 0) THEN
1198 WRITE (out_sp_ab, '(A,I6,A,I6,A,I6)') " Absorption spectrum for atom ", iatom, &
1199 ", index of excited core MO is", estate, ", # of lines ", SIZE(sp_ab, 2)
1200 ene2 = 1.0_dp
1201 DO i = 1, SIZE(sp_ab, 2)
1202 istate = lfomo - 1 + i
1203 IF (length) ene2 = sp_ab(1, i)*sp_ab(1, i)
1204 WRITE (out_sp_ab, '(I6,5F16.8,F10.5)') istate, sp_ab(1, i)*evolt, &
1205 sp_ab(2, i)*ene2, sp_ab(3, i)*ene2, &
1206 sp_ab(4, i)*ene2, sp_ab(5, i)*ene2, sp_ab(6, i)
1207 END DO
1208 END IF
1209
1210 CALL cp_print_key_finished_output(out_sp_ab, logger, xas_section, &
1211 "PRINT%XAS_SPECTRUM")
1212
1213 END SUBROUTINE xas_write
1214
1215! **************************************************************************************************
1216!> \brief write the cube files for a set of selected states
1217!> \param xas_control provide number ant indexes of the states to be printed
1218!> \param qs_env ...
1219!> \param xas_section ...
1220!> \param mos mos from which the states to be printed are extracted
1221!> \param all_vectors ...
1222!> \param iatom index of the atom that has been excited
1223!> \param append_cube ...
1224!> \par History
1225!> 08.2005 created [MI]
1226!> \author MI
1227! **************************************************************************************************
1228 SUBROUTINE xas_print_cubes(xas_control, qs_env, xas_section, &
1229 mos, all_vectors, iatom, append_cube)
1230
1232 TYPE(qs_environment_type), POINTER :: qs_env
1233 TYPE(section_vals_type), POINTER :: xas_section
1234 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1235 TYPE(cp_fm_type), INTENT(IN) :: all_vectors
1236 INTEGER, INTENT(IN) :: iatom
1237 LOGICAL, INTENT(IN) :: append_cube
1238
1239 CHARACTER(LEN=default_string_length) :: my_mittle, my_pos
1240 INTEGER :: homo, istate0, my_spin, nspins, nstates
1241 REAL(dp), DIMENSION(:, :), POINTER :: centers
1242 TYPE(section_vals_type), POINTER :: print_key
1243
1244 nspins = SIZE(mos)
1245
1246 print_key => section_vals_get_subs_vals(xas_section, "PRINT%CLS_FUNCTION_CUBES")
1247 my_mittle = 'at'//trim(adjustl(cp_to_string(iatom)))
1248 nstates = SIZE(xas_control%list_cubes, 1)
1249
1250 IF (xas_control%do_centers) THEN
1251 ! one might like to calculate the centers of the xas orbital (without localizing them)
1252 ELSE
1253 ALLOCATE (centers(6, nstates))
1254 centers = 0.0_dp
1255 END IF
1256 my_spin = xas_control%spin_channel
1257
1258 CALL get_mo_set(mos(my_spin), homo=homo)
1259 istate0 = 0
1260
1261 my_pos = "REWIND"
1262 IF (append_cube) THEN
1263 my_pos = "APPEND"
1264 END IF
1265
1266 CALL qs_print_cubes(qs_env, all_vectors, nstates, xas_control%list_cubes, &
1267 centers, print_key, my_mittle, state0=istate0, file_position=my_pos)
1268
1269 DEALLOCATE (centers)
1270
1271 END SUBROUTINE xas_print_cubes
1272
1273! **************************************************************************************************
1274!> \brief write the PDOS after the XAS SCF, i.e., with one excited core
1275!> \param qs_env ...
1276!> \param xas_section ...
1277!> \param mos mos from which the eigenvalues and expansion coeffiecients are obtained
1278!> \param iatom index of the atom that has been excited
1279!> \par History
1280!> 03.2016 created [MI]
1281!> \author MI
1282! **************************************************************************************************
1283
1284 SUBROUTINE xas_pdos(qs_env, xas_section, mos, iatom)
1285
1286 TYPE(qs_environment_type), POINTER :: qs_env
1287 TYPE(section_vals_type), POINTER :: xas_section
1288 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1289 INTEGER, INTENT(IN) :: iatom
1290
1291 CHARACTER(LEN=default_string_length) :: xas_mittle
1292 INTEGER :: ispin
1293 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1294 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1295 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1296
1297 NULLIFY (atomic_kind_set, particle_set, qs_kind_set)
1298 xas_mittle = 'xasat'//trim(adjustl(cp_to_string(iatom)))//'_'
1299
1300 CALL get_qs_env(qs_env, &
1301 atomic_kind_set=atomic_kind_set, &
1302 particle_set=particle_set, &
1303 qs_kind_set=qs_kind_set)
1304
1305 DO ispin = 1, 2
1306 CALL calculate_projected_dos(mos(ispin), atomic_kind_set, qs_kind_set, particle_set, qs_env, &
1307 xas_section, ispin, xas_mittle)
1308 END DO
1309
1310 END SUBROUTINE xas_pdos
1311! **************************************************************************************************
1312!> \brief Calculation of the spectrum when the dipole approximation
1313!> in the velocity form is used.
1314!> \param fm_set components of the position operator in a full matrix form
1315!> already multiplied by the coefficiets
1316!> only the terms <C_i Op C_f> are calculated where
1317!> C_i are the coefficients of the excited state
1318!> \param op_sm components of the position operator for the dipole
1319!> in a sparse matrix form (cos and sin)
1320!> calculated for the basis functions
1321!> \param mos wavefunctions coefficients
1322!> \param excvec coefficients of the excited orbital
1323!> \param all_vectors ...
1324!> \param all_evals ...
1325!> \param sp_em ...
1326!> \param sp_ab ...
1327!> \param estate index of the excited state
1328!> \param nstate ...
1329!> \param my_spin ...
1330!> \par History
1331!> 06.2005 created [MI]
1332!> \author MI
1333! **************************************************************************************************
1334 SUBROUTINE spectrum_dip_vel(fm_set, op_sm, mos, excvec, &
1335 all_vectors, all_evals, sp_em, sp_ab, estate, nstate, my_spin)
1336
1337 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: fm_set
1338 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_sm
1339 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1340 TYPE(cp_fm_type), INTENT(IN) :: excvec, all_vectors
1341 REAL(dp), DIMENSION(:), POINTER :: all_evals
1342 REAL(dp), DIMENSION(:, :), POINTER :: sp_em, sp_ab
1343 INTEGER, INTENT(IN) :: estate, nstate, my_spin
1344
1345 INTEGER :: homo, i, i_abs, istate, lfomo, nao, nmo
1346 REAL(dp) :: dip(3), ene_f, ene_i
1347 REAL(dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
1348 TYPE(cp_fm_type) :: fm_work
1349
1350 cpassert(ASSOCIATED(fm_set))
1351 NULLIFY (eigenvalues, occupation_numbers)
1352
1353 CALL get_mo_set(mos(my_spin), eigenvalues=eigenvalues, occupation_numbers=occupation_numbers, &
1354 nao=nao, nmo=nmo, homo=homo, lfomo=lfomo)
1355
1356 CALL cp_fm_create(fm_work, all_vectors%matrix_struct)
1357 DO i = 1, SIZE(fm_set, 2)
1358 CALL cp_fm_set_all(fm_set(my_spin, i), 0.0_dp)
1359 CALL cp_fm_set_all(fm_work, 0.0_dp)
1360 CALL cp_dbcsr_sm_fm_multiply(op_sm(i)%matrix, all_vectors, fm_work, ncol=nstate)
1361 CALL parallel_gemm("T", "N", 1, nstate, nao, 1.0_dp, excvec, &
1362 fm_work, 0.0_dp, fm_set(my_spin, i), b_first_col=1)
1363 END DO
1364 CALL cp_fm_release(fm_work)
1365
1366 sp_em = 0.0_dp
1367 sp_ab = 0.0_dp
1368 ene_i = eigenvalues(estate)
1369 DO istate = 1, nstate
1370 ene_f = all_evals(istate)
1371 DO i = 1, 3
1372 CALL cp_fm_get_element(fm_set(my_spin, i), 1, istate, dip(i))
1373 END DO
1374 IF (istate <= homo) THEN
1375 sp_em(1, istate) = ene_f - ene_i
1376 sp_em(2, istate) = dip(1)
1377 sp_em(3, istate) = dip(2)
1378 sp_em(4, istate) = dip(3)
1379 sp_em(5, istate) = dip(1)*dip(1) + dip(2)*dip(2) + dip(3)*dip(3)
1380 sp_em(6, istate) = occupation_numbers(istate)
1381 END IF
1382 IF (istate >= lfomo) THEN
1383 i_abs = istate - lfomo + 1
1384 sp_ab(1, i_abs) = ene_f - ene_i
1385 sp_ab(2, i_abs) = dip(1)
1386 sp_ab(3, i_abs) = dip(2)
1387 sp_ab(4, i_abs) = dip(3)
1388 sp_ab(5, i_abs) = dip(1)*dip(1) + dip(2)*dip(2) + dip(3)*dip(3)
1389 IF (istate <= nmo) sp_ab(6, i_abs) = occupation_numbers(istate)
1390 END IF
1391
1392 END DO
1393
1394 END SUBROUTINE spectrum_dip_vel
1395
1396! **************************************************************************************************
1397!> \brief ...
1398!> \param base_a ...
1399!> \param base_b ...
1400!> \param matrix ...
1401! **************************************************************************************************
1402 SUBROUTINE calc_stogto_overlap(base_a, base_b, matrix)
1403
1404 TYPE(gto_basis_set_type), POINTER :: base_a, base_b
1405 REAL(dp), DIMENSION(:, :), POINTER :: matrix
1406
1407 INTEGER :: iset, jset, ldsab, maxcoa, maxcob, maxl, &
1408 maxla, maxlb, na, nb, nseta, nsetb, &
1409 nsgfa, nsgfb, sgfa, sgfb
1410 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1411 npgfb, nsgfa_set, nsgfb_set
1412 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1413 REAL(dp) :: rab(3)
1414 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
1415 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, sphi_a, &
1416 sphi_b, zeta, zetb
1417
1418 NULLIFY (la_max, la_min, lb_max, lb_min)
1419 NULLIFY (npgfa, npgfb, nsgfa_set, nsgfb_set)
1420 NULLIFY (first_sgfa, first_sgfb)
1421 NULLIFY (rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
1422
1423 CALL get_gto_basis_set(gto_basis_set=base_a, nsgf=nsgfa, nsgf_set=nsgfa_set, lmax=la_max, &
1424 lmin=la_min, npgf=npgfa, pgf_radius=rpgfa, &
1425 sphi=sphi_a, scon=scon_a, zet=zeta, first_sgf=first_sgfa, &
1426 maxco=maxcoa, nset=nseta, maxl=maxla)
1427
1428 CALL get_gto_basis_set(gto_basis_set=base_b, nsgf=nsgfb, nsgf_set=nsgfb_set, lmax=lb_max, &
1429 lmin=lb_min, npgf=npgfb, pgf_radius=rpgfb, &
1430 sphi=sphi_b, scon=scon_b, zet=zetb, first_sgf=first_sgfb, &
1431 maxco=maxcob, nset=nsetb, maxl=maxlb)
1432 ! Initialize and allocate
1433 rab = 0.0_dp
1434 matrix = 0.0_dp
1435
1436 ldsab = max(maxcoa, maxcob, nsgfa, nsgfb)
1437 maxl = max(maxla, maxlb)
1438
1439 ALLOCATE (sab(ldsab, ldsab))
1440 ALLOCATE (work(ldsab, ldsab))
1441
1442 DO iset = 1, nseta
1443
1444 na = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1445 sgfa = first_sgfa(1, iset)
1446
1447 DO jset = 1, nsetb
1448 nb = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1449 sgfb = first_sgfb(1, jset)
1450
1451 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1452 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1453 rab, sab)
1454 CALL contraction(sab, work, ca=scon_a(:, sgfa:), na=na, ma=nsgfa_set(iset), &
1455 cb=scon_b(:, sgfb:), nb=nb, mb=nsgfb_set(jset))
1456 CALL block_add("IN", work, nsgfa_set(iset), nsgfb_set(jset), matrix, sgfa, sgfb)
1457
1458 END DO ! jset
1459 END DO ! iset
1460 DEALLOCATE (sab, work)
1461
1462 END SUBROUTINE calc_stogto_overlap
1463
1464! **************************************************************************************************
1465!> \brief Starting from a set of mos, determine on which atom are centered
1466!> and if they are of the right type (1s,2s ...)
1467!> to be used in the specific core level spectrum calculation
1468!> The set of states need to be from the core, otherwise the
1469!> characterization of the type is not valid, since it assumes that
1470!> the orbital is localizad on a specific atom
1471!> It is probably reccomandable to run a localization cycle before
1472!> proceeding to the assignment of the type
1473!> The type is determined by computing the overalp with a
1474!> type specific, minimal, STO bais set
1475!> \param xas_control ...
1476!> \param xas_env ...
1477!> \param localized_wfn_control ...
1478!> \param qs_env ...
1479!> \par History
1480!> 03.2006 created [MI]
1481!> \author MI
1482! **************************************************************************************************
1483 SUBROUTINE cls_assign_core_states(xas_control, xas_env, localized_wfn_control, qs_env)
1484
1486 TYPE(xas_environment_type), POINTER :: xas_env
1487 TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
1488 TYPE(qs_environment_type), POINTER :: qs_env
1489
1490 INTEGER :: chosen_state, homo, i, iat, iatom, &
1491 ikind, isgf, istate, j, my_kind, &
1492 my_spin, nao, natom, nexc_atoms, &
1493 nexc_search, output_unit
1494 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf
1495 INTEGER, DIMENSION(3) :: perd0
1496 INTEGER, DIMENSION(:), POINTER :: atom_of_state, mykind_of_kind, &
1497 nexc_states, state_of_mytype, &
1498 type_of_state
1499 INTEGER, DIMENSION(:, :), POINTER :: state_of_atom
1500 REAL(dp) :: component, dist, distmin, maxocc, ra(3), &
1501 rac(3), rc(3)
1502 REAL(dp), DIMENSION(:), POINTER :: max_overlap, sto_state_overlap
1503 REAL(dp), DIMENSION(:, :), POINTER :: centers_wfn
1504 REAL(kind=dp), DIMENSION(:, :), POINTER :: vecbuffer
1505 TYPE(atomic_kind_type), POINTER :: atomic_kind
1506 TYPE(cell_type), POINTER :: cell
1507 TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: stogto_overlap
1508 TYPE(cp_fm_type), POINTER :: mo_coeff
1509 TYPE(cp_logger_type), POINTER :: logger
1510 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1511 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1512 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1513
1514 NULLIFY (cell, mos, particle_set)
1515 NULLIFY (atom_of_state, centers_wfn, mykind_of_kind, state_of_atom, nexc_states)
1516 NULLIFY (stogto_overlap, type_of_state, max_overlap, qs_kind_set)
1517 NULLIFY (state_of_mytype, type_of_state, sto_state_overlap)
1518
1519 NULLIFY (logger)
1520 logger => cp_get_default_logger()
1521 output_unit = cp_logger_get_default_io_unit(logger)
1522
1523 CALL get_qs_env(qs_env=qs_env, cell=cell, mos=mos, particle_set=particle_set, &
1524 qs_kind_set=qs_kind_set)
1525
1526 ! The Berry operator can be used only for periodic systems
1527 ! If an isolated system is used the periodicity is overimposed
1528 perd0(1:3) = cell%perd(1:3)
1529 cell%perd(1:3) = 1
1530
1531 CALL get_xas_env(xas_env=xas_env, &
1532 centers_wfn=centers_wfn, atom_of_state=atom_of_state, &
1533 mykind_of_kind=mykind_of_kind, &
1534 type_of_state=type_of_state, state_of_atom=state_of_atom, &
1535 stogto_overlap=stogto_overlap, nexc_atoms=nexc_atoms, &
1536 spin_channel=my_spin, nexc_search=nexc_search, nexc_states=nexc_states)
1537
1538 CALL get_mo_set(mos(my_spin), mo_coeff=mo_coeff, maxocc=maxocc, nao=nao, homo=homo)
1539
1540 ! scratch array for the state
1541 ALLOCATE (vecbuffer(1, nao))
1542 natom = SIZE(particle_set)
1543
1544 ALLOCATE (first_sgf(natom))
1545 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf)
1546 ALLOCATE (sto_state_overlap(nexc_search))
1547 ALLOCATE (max_overlap(natom))
1548 max_overlap = 0.0_dp
1549 ALLOCATE (state_of_mytype(natom))
1550 state_of_mytype = 0
1551 atom_of_state = 0
1552 nexc_states = 1
1553 state_of_atom = 0
1554
1555 IF (xas_control%orbital_list(1) < 0) THEN !Checks for manually selected orbitals from the localized set
1556
1557 DO istate = 1, nexc_search
1558 centers_wfn(1, istate) = localized_wfn_control%centers_set(my_spin)%array(1, istate)
1559 centers_wfn(2, istate) = localized_wfn_control%centers_set(my_spin)%array(2, istate)
1560 centers_wfn(3, istate) = localized_wfn_control%centers_set(my_spin)%array(3, istate)
1561
1562 ! Assign the state to the closest atom
1563 distmin = 100.0_dp
1564 DO iat = 1, nexc_atoms
1565 iatom = xas_control%exc_atoms(iat)
1566 ra(1:3) = particle_set(iatom)%r(1:3)
1567 rc(1:3) = centers_wfn(1:3, istate)
1568 rac = pbc(ra, rc, cell)
1569 dist = rac(1)*rac(1) + rac(2)*rac(2) + rac(3)*rac(3)
1570
1571 IF (dist < distmin) THEN
1572
1573 atom_of_state(istate) = iatom
1574 distmin = dist
1575 END IF
1576 END DO
1577 IF (atom_of_state(istate) /= 0) THEN
1578 !Character of the state
1579 CALL cp_fm_get_submatrix(mo_coeff, vecbuffer, 1, istate, &
1580 nao, 1, transpose=.true.)
1581
1582 iatom = atom_of_state(istate)
1583
1584 NULLIFY (atomic_kind)
1585 atomic_kind => particle_set(iatom)%atomic_kind
1586 CALL get_atomic_kind(atomic_kind=atomic_kind, &
1587 kind_number=ikind)
1588
1589 my_kind = mykind_of_kind(ikind)
1590
1591 sto_state_overlap(istate) = 0.0_dp
1592 DO i = 1, SIZE(stogto_overlap(my_kind)%array, 1)
1593 component = 0.0_dp
1594 DO j = 1, SIZE(stogto_overlap(my_kind)%array, 2)
1595 isgf = first_sgf(iatom) + j - 1
1596 component = component + stogto_overlap(my_kind)%array(i, j)*vecbuffer(1, isgf)
1597 END DO
1598 sto_state_overlap(istate) = sto_state_overlap(istate) + &
1599 component*component
1600 END DO
1601
1602 IF (sto_state_overlap(istate) > max_overlap(iatom)) THEN
1603 state_of_mytype(iatom) = istate
1604 max_overlap(iatom) = sto_state_overlap(istate)
1605 END IF
1606 END IF
1607 END DO ! istate
1608
1609 ! Includes all states within the chosen threshold relative to the maximum overlap
1610 IF (xas_control%overlap_threshold < 1) THEN
1611 DO iat = 1, nexc_atoms
1612 iatom = xas_control%exc_atoms(iat)
1613 DO istate = 1, nexc_search
1614 IF (atom_of_state(istate) == iatom) THEN
1615 IF (sto_state_overlap(istate) > max_overlap(iatom)*xas_control%overlap_threshold &
1616 .AND. istate /= state_of_mytype(iat)) THEN
1617 nexc_states(iat) = nexc_states(iat) + 1
1618 state_of_atom(iat, nexc_states(iat)) = istate
1619 END IF
1620 END IF
1621 END DO
1622 END DO
1623 END IF
1624
1625 ! In the set of states, assign the index of the state to be excited for iatom
1626 IF (output_unit > 0) THEN
1627 WRITE (unit=output_unit, fmt="(/,T10,A,/)") &
1628 "List the atoms to be excited and the relative of MOs index "
1629 END IF
1630
1631 DO iat = 1, nexc_atoms
1632 iatom = xas_env%exc_atoms(iat)
1633 state_of_atom(iat, 1) = state_of_mytype(iatom) ! Place the state with maximum overlap first in the list
1634 IF (output_unit > 0) THEN
1635 WRITE (unit=output_unit, fmt="(T10,A,I3,T26,A)", advance='NO') &
1636 'Atom: ', iatom, "MO index:"
1637 END IF
1638 DO istate = 1, nexc_states(iat)
1639 IF (istate < nexc_states(iat)) THEN
1640 IF (output_unit > 0) WRITE (unit=output_unit, fmt="(I4)", advance='NO') state_of_atom(iat, istate)
1641 ELSE
1642 IF (output_unit > 0) WRITE (unit=output_unit, fmt="(I4)") state_of_atom(iat, istate)
1643 END IF
1644 END DO
1645 IF (state_of_atom(iat, 1) == 0 .OR. state_of_atom(iat, 1) > homo) THEN
1646 cpabort("A wrong state has been selected for excitation, check the Wannier centers")
1647 END IF
1648 END DO
1649
1650 IF (xas_control%overlap_threshold < 1) THEN
1651 DO iat = 1, nexc_atoms
1652 IF (output_unit > 0) THEN
1653 WRITE (unit=output_unit, fmt="(/,T10,A,I6)") &
1654 'Overlap integrals for Atom: ', iat
1655 DO istate = 1, nexc_states(iat)
1656 WRITE (unit=output_unit, fmt="(T10,A,I3,T26,A,T38,f10.8)") &
1657 'State: ', state_of_atom(iat, istate), "Overlap:", sto_state_overlap(state_of_atom(iat, istate))
1658 END DO
1659 END IF
1660 END DO
1661 END IF
1662
1663 ! Scales down the 2d-array to the minimal size
1664 CALL reallocate(xas_env%state_of_atom, 1, nexc_atoms, 1, maxval(nexc_states))
1665
1666 ELSE ! Manually selected orbital indices
1667
1668 ! Reallocate nexc_states and state_of_atom to include any atom
1669 CALL reallocate(xas_env%nexc_states, 1, natom)
1670 CALL reallocate(xas_env%state_of_atom, 1, natom, 1, SIZE(xas_control%orbital_list))
1671 CALL get_xas_env(xas_env, nexc_states=nexc_states, state_of_atom=state_of_atom)
1672
1673 nexc_states = 0
1674 state_of_atom = 0
1675 nexc_atoms = natom !To include all possible atoms in the spectrum calculation
1676
1677 DO istate = 1, SIZE(xas_control%orbital_list)
1678
1679 chosen_state = xas_control%orbital_list(istate)
1680 nexc_atoms = 1
1681 centers_wfn(1, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(1, chosen_state)
1682 centers_wfn(2, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(2, chosen_state)
1683 centers_wfn(3, chosen_state) = localized_wfn_control%centers_set(my_spin)%array(3, chosen_state)
1684
1685 distmin = 100.0_dp
1686 DO iat = 1, natom
1687 ra(1:3) = particle_set(iat)%r(1:3)
1688 rc(1:3) = centers_wfn(1:3, chosen_state)
1689 rac = pbc(ra, rc, cell)
1690 dist = rac(1)*rac(1) + rac(2)*rac(2) + rac(3)*rac(3)
1691 IF (dist < distmin) THEN
1692 atom_of_state(chosen_state) = iat !?
1693 distmin = dist
1694 END IF
1695 END DO ! iat
1696
1697 nexc_states(atom_of_state(chosen_state)) = nexc_states(atom_of_state(chosen_state)) + 1
1698 state_of_atom(atom_of_state(chosen_state), nexc_states(atom_of_state(chosen_state))) = chosen_state
1699
1700 END DO !istate
1701
1702 ! In the set of states, assign the index of the state to be excited for iatom
1703 IF (output_unit > 0) THEN
1704 WRITE (unit=output_unit, fmt="(/,T10,A,/)") &
1705 "List the atoms to be excited and the relative of MOs index "
1706 END IF
1707
1708 DO iat = 1, natom
1709 IF (output_unit > 0 .AND. state_of_atom(iat, 1) /= 0) THEN
1710 WRITE (unit=output_unit, fmt="(T10,A,I3,T26,A)", advance='NO') &
1711 'Atom: ', iat, "MO index:"
1712 DO i = 1, nexc_states(iat)
1713 IF (i < nexc_states(iat)) THEN
1714 WRITE (unit=output_unit, fmt="(I4)", advance='NO') state_of_atom(iat, i)
1715 ELSE
1716 WRITE (unit=output_unit, fmt="(I4)") state_of_atom(iat, i)
1717 END IF
1718 END DO
1719 END IF
1720 IF (state_of_atom(iat, 1) > homo) THEN
1721 cpabort("A wrong state has been selected for excitation, check the Wannier centers")
1722 END IF
1723 END DO
1724
1725 CALL reallocate(xas_env%state_of_atom, 1, natom, 1, maxval(nexc_states)) ! Scales down the 2d-array to the minimal size
1726
1727 END IF !Checks for manually selected orbitals from the localized set
1728
1729 ! Set back the correct periodicity
1730 cell%perd(1:3) = perd0(1:3)
1731
1732 DEALLOCATE (vecbuffer)
1733 DEALLOCATE (first_sgf)
1734 DEALLOCATE (sto_state_overlap)
1735 DEALLOCATE (max_overlap)
1736 DEALLOCATE (state_of_mytype)
1737
1738 END SUBROUTINE cls_assign_core_states
1739
1740END MODULE xas_methods
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:681
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
pure real(dp) function, public srules(z, ne, n, l)
...
subroutine, public deallocate_sto_basis_set(sto_basis_set)
...
subroutine, public allocate_sto_basis_set(sto_basis_set)
...
subroutine, public create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss, ortho)
...
subroutine, public set_sto_basis_set(sto_basis_set, name, nshell, symbol, nq, lq, zet)
...
subroutine, public init_orb_basis_set(gto_basis_set)
Initialise a Gaussian-type orbital (GTO) basis set data set.
Handles all functions related to the CELL.
Definition cell_types.F:15
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_set(matrix, alpha)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Routines to handle the external control of CP2K.
subroutine, public external_control(should_stop, flag, globenv, target_time, start_time, force_check)
External manipulations during a run : when the <PROJECT_NAME>.EXIT_$runtype command is sent the progr...
pool for for elements that are retained and released
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public xas_1s_type
integer, parameter, public xas_4p_type
integer, parameter, public xas_3s_type
integer, parameter, public do_loc_none
integer, parameter, public xas_3d_type
integer, parameter, public xas_dscf
integer, parameter, public xas_4s_type
integer, parameter, public xas_tp_xhh
integer, parameter, public xas_2p_type
integer, parameter, public xas_dip_len
integer, parameter, public xas_dip_vel
integer, parameter, public xas_2s_type
integer, parameter, public xas_tp_xfh
integer, parameter, public xas_3p_type
integer, parameter, public xas_tp_fh
integer, parameter, public xas_tp_flex
integer, parameter, public xas_4f_type
integer, parameter, public state_loc_list
integer, parameter, public state_loc_range
integer, parameter, public xas_4d_type
integer, parameter, public xas_tp_hh
integer, parameter, public xes_tp_val
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
logical function, public section_get_lval(section_vals, keyword_name)
...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
Utility routines for the memory handling.
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
basic linear algebra operations for full matrixes
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Define the data structure for the particle information.
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public evolt
Definition physcon.F:183
Apply the direct inversion in the iterative subspace (DIIS) of Pulay in the framework of an SCF itera...
Definition qs_diis.F:21
pure subroutine, public qs_diis_b_clear(diis_buffer)
clears the buffer
Definition qs_diis.F:519
subroutine, public qs_diis_b_create(diis_buffer, nbuffer)
Allocates an SCF DIIS buffer.
Definition qs_diis.F:103
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, mimic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Driver for the localization that should be general for all the methods available and all the definiti...
Definition qs_loc_main.F:22
subroutine, public qs_loc_driver(qs_env, qs_loc_env, print_loc_section, myspin, ext_mo_coeff)
set up the calculation of localized orbitals
Definition qs_loc_main.F:96
Driver for the localization that should be general for all the methods available and all the definiti...
subroutine, public qs_print_cubes(qs_env, mo_coeff, nstates, state_list, centers, print_key, root, ispin, idir, state0, file_position)
write the cube files for a set of selected states
New version of the module for the localization of the molecular orbitals This should be able to use d...
subroutine, public qs_loc_env_create(qs_loc_env)
...
Some utilities for the construction of the localization environment.
subroutine, public set_loc_wfn_lists(localized_wfn_control, nmoloc, nmo, nspins, my_spin)
create the lists of mos that are taken into account
subroutine, public qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, myspin, do_localize, loc_coeff, mo_loc_history)
allocates the data, and initializes the operators
subroutine, public set_loc_centers(localized_wfn_control, nmoloc, nspins)
create the center and spread array and the file names for the output
subroutine, public qs_loc_control_init(qs_loc_env, loc_section, do_homo, do_mixed, do_xas, nloc_xas, spin_xas)
initializes everything needed for localization of the HOMOs
wrapper for the pools of matrixes
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Definition and initialisation of the mo data type.
Definition qs_mo_io.F:21
subroutine, public write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, dft_section, before, kpoint, final_mos, spin, solver_method, rtp, cpart, sim_step, umo_set, qs_env)
Write MO information to output file (eigenvalues, occupation numbers, coefficients)
Definition qs_mo_io.F:1023
collects routines that perform operations directly related to MOs
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public set_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, uniform_occupation, kts, mu, flexible_electron_count)
Set the components of a MO set data structure.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Define the neighbor list data types and the corresponding functionality.
subroutine, public p_xyz_ao(op, qs_env, minimum_image)
Calculation of the components of the dipole operator in the velocity form The elements of the sparse ...
subroutine, public rrc_xyz_ao(op, qs_env, rc, order, minimum_image, soft)
Calculation of the components of the dipole operator in the length form by taking the relative positi...
Calculation and writing of projected density of states The DOS is computed per angular momentum and p...
Definition qs_pdos.F:15
subroutine, public calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, qs_env, dft_section, ispin, xas_mittle, external_matrix_shalf, unoccupied_orbs, unoccupied_evals, pdos_print_key, write_pdos, write_pdos_curve)
Compute and write projected density of states.
Definition qs_pdos.F:155
Utility routines for qs_scf.
subroutine, public qs_scf_env_initialize(qs_env, scf_env, scf_control, scf_section)
initializes input parameters if needed or restores values from previous runs to fill scf_env with the...
module that contains the definitions of the scf types
subroutine, public scf_env_release(scf_env)
releases an scf_env (see doc/ReferenceCounting.html)
Routines for the Quickstep SCF run.
Definition qs_scf.F:47
subroutine, public scf_env_cleanup(scf_env)
perform cleanup operations (like releasing temporary storage) at the end of the scf
Definition qs_scf.F:1034
parameters that control an scf iteration
subroutine, public scf_c_read_parameters(scf_control, inp_section)
reads the parameters of the scf section into the given scf_control
subroutine, public scf_c_create(scf_control)
allocates and initializes an scf control object with the default values
Defines control structures, which contain the parameters and the settings for the calculations.
Definition xas_control.F:12
subroutine, public read_xas_control(xas_control, xas_section)
read from input the instructions for a xes/xas calculation
Definition xas_control.F:85
subroutine, public xas_control_create(xas_control)
create retain release the xas_control_type
subroutine, public write_xas_control(xas_control, dft_section)
write on the instructions for a xes/xas calculation
define create destroy get and put information in xas_env to calculate the x-ray absorption spectra
subroutine, public xas_env_release(xas_env)
...
subroutine, public set_xas_env(xas_env, nexc_search, spin_channel, nexc_atoms, nvirtual, nvirtual2, ip_energy, occ_estate, qs_loc_env, xas_estate, xas_nelectron, homo_occ, scf_env, scf_control)
...
subroutine, public xas_env_create(xas_env)
...
subroutine, public get_xas_env(xas_env, exc_state, nao, nvirtual, nvirtual2, centers_wfn, atom_of_state, exc_atoms, nexc_states, type_of_state, mykind_of_atom, mykind_of_kind, state_of_atom, spectrum, groundstate_coeff, ostrength_sm, dip_fm_set, excvec_coeff, excvec_overlap, unoccupied_orbs, unoccupied_evals, unoccupied_max_iter, unoccupied_eps, all_vectors, all_evals, my_gto_basis, qs_loc_env, stogto_overlap, occ_estate, xas_nelectron, xas_estate, nexc_atoms, nexc_search, spin_channel, scf_env, scf_control)
...
driver for the xas calculation and xas_scf for the tp method
Definition xas_methods.F:15
subroutine, public calc_stogto_overlap(base_a, base_b, matrix)
...
subroutine, public xas(qs_env, dft_control)
Driver for xas calculations The initial mos are prepared A loop on the atoms to be excited is started...
Initialize the XAS orbitals for specific core excitations Either the GS orbitals are used as initial ...
Definition xas_restart.F:20
subroutine, public xas_read_restart(xas_env, xas_section, qs_env, xas_method, iatom, estate, istate)
Set up for reading the restart corresponding to the excitation of iatom If the corresponding restart ...
xas_scf for the tp method It is repeaated for every atom that have to be excited
Definition xas_tp_scf.F:15
subroutine, public xas_do_tp_scf(dft_control, xas_env, iatom, istate, scf_env, qs_env, xas_section, scf_section, converged, should_stop)
perform an scf loop to calculate the xas spectrum given by the excitation of a inner state of a selec...
Definition xas_tp_scf.F:137
subroutine, public xes_scf_once(qs_env, xas_env, converged, should_stop)
SCF for emission spectra calculations: vacancy in valence.
Definition xas_tp_scf.F:643
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
represent a pointer to a 2d array
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.
A type that holds controlling information for the calculation of the spread of wfn and the optimizati...
contains all the info needed by quickstep to calculate the spread of a selected set of orbitals and i...
container for the pools of matrixes used by qs
A type that holds controlling information for a xas calculation.
Definition xas_control.F:40