33 #include "./base/base_uses.f90"
39 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_mo_occupation'
41 PUBLIC :: set_mo_occupation
43 INTERFACE set_mo_occupation
44 MODULE PROCEDURE set_mo_occupation_1, set_mo_occupation_2
59 SUBROUTINE set_mo_occupation_3(mo_array, smear)
61 TYPE(mo_set_type),
DIMENSION(2),
INTENT(INOUT) :: mo_array
62 TYPE(smear_type),
POINTER :: smear
64 CHARACTER(LEN=*),
PARAMETER :: routineN =
'set_mo_occupation_3'
66 INTEGER :: all_nmo, handle, homo_a, homo_b, i, &
67 lfomo_a, lfomo_b, nmo_a, nmo_b, &
69 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: all_index
71 REAL(KIND=
dp) :: all_nelec, kts, mu, nelec_a, nelec_b, &
73 REAL(KIND=
dp),
ALLOCATABLE,
DIMENSION(:) :: all_eigval, all_occ
74 REAL(KIND=
dp),
DIMENSION(:),
POINTER :: eigval_a, eigval_b, occ_a, occ_b
76 CALL timeset(routinen, handle)
78 NULLIFY (eigval_a, eigval_b, occ_a, occ_b)
79 CALL get_mo_set(mo_set=mo_array(1), nmo=nmo_a, eigenvalues=eigval_a, &
80 occupation_numbers=occ_a)
81 CALL get_mo_set(mo_set=mo_array(2), nmo=nmo_b, eigenvalues=eigval_b, &
82 occupation_numbers=occ_b)
83 all_nmo = nmo_a + nmo_b
84 ALLOCATE (all_eigval(all_nmo))
85 ALLOCATE (all_occ(all_nmo))
86 ALLOCATE (all_index(all_nmo))
88 all_eigval(1:nmo_a) = eigval_a(1:nmo_a)
89 all_eigval(nmo_a + 1:all_nmo) = eigval_b(1:nmo_b)
91 CALL sort(all_eigval, all_nmo, all_index)
99 nelec_a = accurate_sum(occ_a(:))
100 nelec_b = accurate_sum(occ_b(:))
101 all_nelec = nelec_a + nelec_b
104 IF (all_index(i) <= nmo_a)
THEN
105 all_occ(i) = occ_a(all_index(i))
107 all_occ(i) = occ_b(all_index(i) - nmo_a)
111 CALL fermifixed(all_occ, mu, kts, all_eigval, all_nelec, &
112 smear%electronic_temperature, 1._dp, xas_estate, occ_estate)
114 is_large = abs(maxval(all_occ) - 1.0_dp) > smear%eps_fermi_dirac
117 cpwarn(
"Fermi-Dirac smearing includes the first MO")
119 is_large = abs(minval(all_occ)) > smear%eps_fermi_dirac
121 CALL cp_warn(__location__, &
122 "Fermi-Dirac smearing includes the last MO => "// &
123 "Add more MOs for proper smearing.")
126 is_large = (abs(all_nelec - accurate_sum(all_occ(:))) > smear%eps_fermi_dirac*all_nelec)
128 cpwarn(
"Total number of electrons is not accurate")
131 IF (all_index(i) <= nmo_a)
THEN
132 occ_a(all_index(i)) = all_occ(i)
133 eigval_a(all_index(i)) = all_eigval(i)
135 occ_b(all_index(i) - nmo_a) = all_occ(i)
136 eigval_b(all_index(i) - nmo_a) = all_eigval(i)
140 nelec_a = accurate_sum(occ_a(:))
141 nelec_b = accurate_sum(occ_b(:))
144 IF (occ_a(i) < 1.0_dp)
THEN
150 IF (occ_b(i) < 1.0_dp)
THEN
156 DO i = nmo_a, lfomo_a, -1
157 IF (occ_a(i) > smear%eps_fermi_dirac)
THEN
163 DO i = nmo_b, lfomo_b, -1
164 IF (occ_b(i) > smear%eps_fermi_dirac)
THEN
170 CALL set_mo_set(mo_set=mo_array(1), kts=kts/2.0_dp, mu=mu, n_el_f=nelec_a, &
171 lfomo=lfomo_a, homo=homo_a, uniform_occupation=.false.)
172 CALL set_mo_set(mo_set=mo_array(2), kts=kts/2.0_dp, mu=mu, n_el_f=nelec_b, &
173 lfomo=lfomo_b, homo=homo_b, uniform_occupation=.false.)
175 CALL timestop(handle)
177 END SUBROUTINE set_mo_occupation_3
196 SUBROUTINE set_mo_occupation_2(mo_array, smear, eval_deriv, tot_zeff_corr)
198 TYPE(mo_set_type),
DIMENSION(:),
INTENT(INOUT) :: mo_array
199 TYPE(smear_type),
POINTER :: smear
200 REAL(KIND=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: eval_deriv
201 REAL(KIND=
dp),
OPTIONAL :: tot_zeff_corr
203 CHARACTER(LEN=*),
PARAMETER :: routineN =
'set_mo_occupation_2'
205 INTEGER :: handle, i, lumo_a, lumo_b, &
206 multiplicity_new, multiplicity_old, &
208 REAL(KIND=
dp) :: nelec_f, threshold
209 REAL(KIND=
dp),
DIMENSION(:),
POINTER :: eigval_a, eigval_b
211 CALL timeset(routinen, handle)
214 IF (
SIZE(mo_array) == 1)
THEN
215 IF (
PRESENT(eval_deriv))
THEN
217 CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv)
219 IF (
PRESENT(tot_zeff_corr))
THEN
220 CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
222 CALL set_mo_occupation_1(mo_array(1), smear=smear)
225 CALL timestop(handle)
229 IF (smear%do_smear)
THEN
230 IF (smear%fixed_mag_mom < 0.0_dp)
THEN
231 IF (
PRESENT(tot_zeff_corr))
THEN
232 CALL cp_warn(__location__, &
233 "CORE_CORRECTION /= 0.0 might cause the cell to charge up "// &
234 "that will lead to application of different background "// &
235 "correction compared to the reference system. "// &
236 "Use FIXED_MAGNETIC_MOMENT >= 0.0 if using SMEAR keyword "// &
237 "to correct the electron density")
239 IF (smear%fixed_mag_mom /= -1.0_dp)
THEN
240 cpassert(.NOT. (
PRESENT(eval_deriv)))
241 CALL set_mo_occupation_3(mo_array, smear=smear)
242 CALL timestop(handle)
246 nelec_f = mo_array(1)%n_el_f + mo_array(2)%n_el_f
247 IF (abs((mo_array(1)%n_el_f - mo_array(2)%n_el_f) - smear%fixed_mag_mom) > smear%eps_fermi_dirac*nelec_f)
THEN
248 mo_array(1)%n_el_f = nelec_f/2.0_dp + smear%fixed_mag_mom/2.0_dp
249 mo_array(2)%n_el_f = nelec_f/2.0_dp - smear%fixed_mag_mom/2.0_dp
251 cpassert(.NOT. (
PRESENT(eval_deriv)))
252 IF (
PRESENT(tot_zeff_corr))
THEN
253 CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
254 CALL set_mo_occupation_1(mo_array(2), smear=smear, tot_zeff_corr=tot_zeff_corr)
256 CALL set_mo_occupation_1(mo_array(1), smear=smear)
257 CALL set_mo_occupation_1(mo_array(2), smear=smear)
262 IF (.NOT. ((mo_array(1)%flexible_electron_count > 0.0_dp) .AND. &
263 (mo_array(2)%flexible_electron_count > 0.0_dp)))
THEN
264 IF (
PRESENT(eval_deriv))
THEN
265 CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv)
266 CALL set_mo_occupation_1(mo_array(2), smear=smear, eval_deriv=eval_deriv)
268 IF (
PRESENT(tot_zeff_corr))
THEN
269 CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
270 CALL set_mo_occupation_1(mo_array(2), smear=smear, tot_zeff_corr=tot_zeff_corr)
272 CALL set_mo_occupation_1(mo_array(1), smear=smear)
273 CALL set_mo_occupation_1(mo_array(2), smear=smear)
276 CALL timestop(handle)
280 nelec = mo_array(1)%nelectron + mo_array(2)%nelectron
282 multiplicity_old = mo_array(1)%nelectron - mo_array(2)%nelectron + 1
284 IF (mo_array(1)%nelectron >= mo_array(1)%nmo) &
285 CALL cp_warn(__location__, &
286 "All alpha MOs are occupied. Add more alpha MOs to "// &
287 "allow for a higher multiplicity")
288 IF ((mo_array(2)%nelectron >= mo_array(2)%nmo) .AND. (mo_array(2)%nelectron /= mo_array(1)%nelectron)) &
289 CALL cp_warn(__location__,
"All beta MOs are occupied. Add more beta MOs to "// &
290 "allow for a lower multiplicity")
292 eigval_a => mo_array(1)%eigenvalues
293 eigval_b => mo_array(2)%eigenvalues
302 threshold = max(mo_array(1)%flexible_electron_count, mo_array(2)%flexible_electron_count)
303 IF ((eigval_a(lumo_a) - threshold) < eigval_b(lumo_b))
THEN
308 IF (lumo_a > mo_array(1)%nmo)
THEN
310 CALL cp_warn(__location__, &
311 "All alpha MOs are occupied. Add more alpha MOs to "// &
312 "allow for a higher multiplicity")
318 IF (lumo_b > mo_array(2)%nmo)
THEN
319 IF (lumo_b < lumo_a) &
320 CALL cp_warn(__location__, &
321 "All beta MOs are occupied. Add more beta MOs to "// &
322 "allow for a lower multiplicity")
330 mo_array(1)%homo = lumo_a - 1
331 mo_array(2)%homo = lumo_b - 1
333 IF (mo_array(2)%homo > mo_array(1)%homo)
THEN
334 CALL cp_warn(__location__, &
336 trim(adjustl(cp_to_string(mo_array(2)%homo)))// &
338 trim(adjustl(cp_to_string(mo_array(1)%homo)))// &
339 ") MOs are occupied. Resorting to low spin state")
340 mo_array(1)%homo = nelec/2 +
modulo(nelec, 2)
341 mo_array(2)%homo = nelec/2
344 mo_array(1)%nelectron = mo_array(1)%homo
345 mo_array(2)%nelectron = mo_array(2)%homo
346 multiplicity_new = mo_array(1)%nelectron - mo_array(2)%nelectron + 1
348 IF (multiplicity_new /= multiplicity_old) &
349 CALL cp_warn(__location__, &
350 "Multiplicity changed from "// &
351 trim(adjustl(cp_to_string(multiplicity_old)))//
" to "// &
352 trim(adjustl(cp_to_string(multiplicity_new))))
354 IF (
PRESENT(eval_deriv))
THEN
355 CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv)
356 CALL set_mo_occupation_1(mo_array(2), smear=smear, eval_deriv=eval_deriv)
358 IF (
PRESENT(tot_zeff_corr))
THEN
359 CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
360 CALL set_mo_occupation_1(mo_array(2), smear=smear, tot_zeff_corr=tot_zeff_corr)
362 CALL set_mo_occupation_1(mo_array(1), smear=smear)
363 CALL set_mo_occupation_1(mo_array(2), smear=smear)
367 CALL timestop(handle)
369 END SUBROUTINE set_mo_occupation_2
389 SUBROUTINE set_mo_occupation_1(mo_set, smear, eval_deriv, xas_env, tot_zeff_corr)
391 TYPE(mo_set_type),
INTENT(INOUT) :: mo_set
392 TYPE(smear_type),
OPTIONAL,
POINTER :: smear
393 REAL(KIND=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: eval_deriv
394 TYPE(xas_environment_type),
OPTIONAL,
POINTER :: xas_env
395 REAL(KIND=
dp),
OPTIONAL :: tot_zeff_corr
397 CHARACTER(LEN=*),
PARAMETER :: routineN =
'set_mo_occupation_1'
399 INTEGER :: handle, i_first, imo, ir, irmo, nmo, &
401 LOGICAL :: equal_size, is_large
402 REAL(KIND=
dp) :: delectron, e1, e2, edelta, edist, &
403 el_count, lengthscale, nelec, &
404 occ_estate, total_zeff_corr, &
406 REAL(KIND=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dfde
408 CALL timeset(routinen, handle)
410 cpassert(
ASSOCIATED(mo_set%eigenvalues))
411 cpassert(
ASSOCIATED(mo_set%occupation_numbers))
412 mo_set%occupation_numbers(:) = 0.0_dp
415 IF (mo_set%nelectron == 0)
THEN
416 CALL timestop(handle)
422 IF (
PRESENT(xas_env))
THEN
423 CALL get_xas_env(xas_env=xas_env, xas_nelectron=xas_nelectron, occ_estate=occ_estate, xas_estate=xas_estate)
424 nomo = ceiling(xas_nelectron + 1.0 - occ_estate - epsilon(0.0_dp))
426 mo_set%occupation_numbers(1:nomo) = mo_set%maxocc
427 IF (xas_estate > 0) mo_set%occupation_numbers(xas_estate) = occ_estate
428 el_count = sum(mo_set%occupation_numbers(1:nomo))
429 IF (el_count > xas_nelectron) &
430 mo_set%occupation_numbers(nomo) = mo_set%occupation_numbers(nomo) - (el_count - xas_nelectron)
431 el_count = sum(mo_set%occupation_numbers(1:nomo))
432 is_large = abs(el_count - xas_nelectron) > xas_nelectron*epsilon(el_count)
433 cpassert(.NOT. is_large)
435 IF (
modulo(mo_set%nelectron, int(mo_set%maxocc)) == 0)
THEN
436 nomo = nint(mo_set%nelectron/mo_set%maxocc)
438 mo_set%occupation_numbers(1:nomo) = mo_set%maxocc
440 nomo = int(mo_set%nelectron/mo_set%maxocc) + 1
442 mo_set%occupation_numbers(1:nomo - 1) = mo_set%maxocc
443 mo_set%occupation_numbers(nomo) = mo_set%nelectron - (nomo - 1)*mo_set%maxocc
450 IF (
PRESENT(tot_zeff_corr))
THEN
452 total_zeff_corr = tot_zeff_corr
453 IF (int(mo_set%maxocc) == 1) total_zeff_corr = total_zeff_corr/2.0_dp
455 IF (total_zeff_corr < 0.0_dp)
THEN
457 delectron = abs(total_zeff_corr) - real(mo_set%maxocc, kind=
dp)
458 IF (delectron > 0.0_dp)
THEN
459 mo_set%occupation_numbers(nomo) = 0.0_dp
460 irmo = ceiling(delectron/real(mo_set%maxocc, kind=
dp))
462 delectron = delectron - real(mo_set%maxocc, kind=
dp)
463 IF (delectron < 0.0_dp)
THEN
464 mo_set%occupation_numbers(nomo - ir) = -delectron
466 mo_set%occupation_numbers(nomo - ir) = 0.0_dp
470 IF (mo_set%occupation_numbers(nomo) == 0.0_dp) nomo = nomo - 1
471 ELSEIF (delectron < 0.0_dp)
THEN
472 mo_set%occupation_numbers(nomo) = -delectron
474 mo_set%occupation_numbers(nomo) = 0.0_dp
477 ELSEIF (total_zeff_corr > 0.0_dp)
THEN
479 delectron = total_zeff_corr - real(mo_set%maxocc, kind=
dp)
480 IF (delectron > 0.0_dp)
THEN
481 mo_set%occupation_numbers(nomo + 1) = real(mo_set%maxocc, kind=
dp)
483 irmo = ceiling(delectron/real(mo_set%maxocc, kind=
dp))
485 delectron = delectron - real(mo_set%maxocc, kind=
dp)
486 IF (delectron < 0.0_dp)
THEN
487 mo_set%occupation_numbers(nomo + ir) = delectron + real(mo_set%maxocc, kind=
dp)
489 mo_set%occupation_numbers(nomo + ir) = real(mo_set%maxocc, kind=
dp)
494 mo_set%occupation_numbers(nomo + 1) = total_zeff_corr
500 nmo =
SIZE(mo_set%eigenvalues)
502 cpassert(nmo >= nomo)
503 cpassert((
SIZE(mo_set%occupation_numbers) == nmo))
506 mo_set%lfomo = nomo + 1
507 mo_set%mu = mo_set%eigenvalues(nomo)
510 IF (
PRESENT(eval_deriv))
THEN
511 equal_size = (
SIZE(mo_set%occupation_numbers, 1) ==
SIZE(eval_deriv, 1))
516 IF (.NOT.
PRESENT(smear))
THEN
518 mo_set%uniform_occupation = .true.
519 IF (
PRESENT(eval_deriv))
THEN
522 CALL timestop(handle)
528 IF ((abs(mo_set%eigenvalues(1)) < 1.0e-12_dp) .AND. &
529 (abs(mo_set%eigenvalues(nmo)) < 1.0e-12_dp))
THEN
530 CALL timestop(handle)
536 IF (smear%do_smear)
THEN
537 IF (
PRESENT(xas_env))
THEN
538 i_first = xas_estate + 1
539 nelec = xas_nelectron
542 IF (smear%fixed_mag_mom == -1.0_dp)
THEN
543 nelec = real(mo_set%nelectron,
dp)
545 nelec = mo_set%n_el_f
548 SELECT CASE (smear%method)
550 IF (.NOT.
PRESENT(eval_deriv))
THEN
551 CALL fermifixed(mo_set%occupation_numbers, mo_set%mu, mo_set%kTS, mo_set%eigenvalues, nelec, &
552 smear%electronic_temperature, mo_set%maxocc, xas_estate, occ_estate)
556 ALLOCATE (dfde(nmo, nmo))
558 lengthscale = 10*smear%electronic_temperature
560 CALL fermifixedderiv(dfde, mo_set%occupation_numbers, mo_set%mu, mo_set%kTS, mo_set%eigenvalues, nelec, &
561 smear%electronic_temperature, mo_set%maxocc, lengthscale, xas_estate, occ_estate)
564 eval_deriv = eval_deriv - mo_set%eigenvalues + mo_set%mu
566 eval_deriv = matmul(transpose(dfde), eval_deriv)
572 DO imo = i_first, nmo
573 IF (mo_set%occupation_numbers(imo) < mo_set%maxocc)
THEN
578 is_large = abs(maxval(mo_set%occupation_numbers) - mo_set%maxocc) > smear%eps_fermi_dirac
581 cpwarn(
"Fermi-Dirac smearing includes the first MO")
584 DO imo = nmo, mo_set%lfomo, -1
585 IF (mo_set%occupation_numbers(imo) > smear%eps_fermi_dirac)
THEN
590 is_large = abs(minval(mo_set%occupation_numbers)) > smear%eps_fermi_dirac
592 CALL cp_warn(__location__, &
593 "Fermi-Dirac smearing includes the last MO => "// &
594 "Add more MOs for proper smearing.")
597 is_large = (abs(nelec - accurate_sum(mo_set%occupation_numbers(:))) > smear%eps_fermi_dirac*nelec)
599 cpwarn(
"Total number of electrons is not accurate")
603 cpassert(.NOT.
PRESENT(eval_deriv))
606 e1 = mo_set%eigenvalues(mo_set%homo) - 0.5_dp*smear%window_size
607 IF (e1 <= mo_set%eigenvalues(1)) &
608 cpwarn(
"Energy window for smearing includes the first MO")
610 e2 = mo_set%eigenvalues(mo_set%homo) + 0.5_dp*smear%window_size
611 IF (e2 >= mo_set%eigenvalues(nmo)) &
612 CALL cp_warn(__location__, &
613 "Energy window for smearing includes the last MO => "// &
614 "Add more MOs for proper smearing.")
617 DO imo = i_first, nomo
618 IF (mo_set%eigenvalues(imo) > e1)
THEN
625 DO imo = nmo, nomo, -1
626 IF (mo_set%eigenvalues(imo) < e2)
THEN
636 DO imo = mo_set%lfomo, mo_set%homo
637 nelec = nelec + mo_set%occupation_numbers(imo)
638 edist = edist + abs(e2 - mo_set%eigenvalues(imo))
642 DO imo = mo_set%lfomo, mo_set%homo
643 edelta = abs(e2 - mo_set%eigenvalues(imo))
644 mo_set%occupation_numbers(imo) = min(mo_set%maxocc, nelec*edelta/edist)
645 nelec = nelec - mo_set%occupation_numbers(imo)
646 edist = edist - edelta
650 equal_size =
SIZE(mo_set%occupation_numbers, 1) ==
SIZE(smear%list, 1)
652 mo_set%occupation_numbers = smear%list
654 IF (
PRESENT(eval_deriv))
THEN
663 IF (mo_set%lfomo == mo_set%homo)
THEN
665 mo_set%lfomo = nomo + 1
667 mo_set%uniform_occupation = .false.
675 CALL timestop(handle)
677 END SUBROUTINE set_mo_occupation_1
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
various routines to log and control the output. The idea is that decisions about where to log should ...
deal with the Fermi distribution, compute it, fix mu, get derivs
subroutine, public fermifixedderiv(dfde, f, mu, kTS, e, N, T, maxocc, l, estate, festate)
returns f and dfde for a given set of energies and number of electrons it is a numerical derivative,...
subroutine, public fermifixed(f, mu, kTS, e, N, T, maxocc, estate, festate)
returns occupations according to Fermi-Dirac statistics for a given set of energies and number of ele...
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Defines the basic variable types.
integer, parameter, public dp
Set occupation of molecular orbitals.
Definition and initialisation of the mo data type.
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.
logical function, public has_uniform_occupation(mo_set, first_mo, last_mo, occupation, tolerance)
Check if the set of MOs in mo_set specifed by the MO index range [first_mo,last_mo] an integer occupa...
parameters that control an scf iteration
All kind of helpful little routines.
define create destroy get and put information in xas_env to calculate the x-ray absorption spectra
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)
...