(git:2d33723)
Loading...
Searching...
No Matches
qs_mo_occupation.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 Set occupation of molecular orbitals
10!> \par History
11!> - set_mo_occupation subroutines moved from qs_mo_types (11.12.2014 MI)
12!> \author MI
13! **************************************************************************************************
14
15MODULE qs_mo_occupation
16
17 USE cp_control_types, ONLY: hairy_probes_type
18 USE cp_log_handling, ONLY: cp_to_string
19 USE hairy_probes, ONLY: probe_occupancy
20 USE input_constants, ONLY: smear_energy_window,&
21 smear_fermi_dirac,&
22 smear_gaussian,&
23 smear_list,&
24 smear_mp,&
25 smear_mv
26 USE kahan_sum, ONLY: accurate_sum
27 USE kinds, ONLY: dp
28 USE qs_mo_types, ONLY: get_mo_set,&
29 has_uniform_occupation,&
30 mo_set_type,&
31 set_mo_set
32 USE scf_control_types, ONLY: smear_type
33 USE smearing_utils, ONLY: smearfixed,&
34 smearfixedderivmv
35 USE util, ONLY: sort
36 USE xas_env_types, ONLY: get_xas_env,&
37 xas_environment_type
38#include "./base/base_uses.f90"
39
40 IMPLICIT NONE
41
42 PRIVATE
43
44 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mo_occupation'
45
46 PUBLIC :: set_mo_occupation
47
48 INTERFACE set_mo_occupation
49 MODULE PROCEDURE set_mo_occupation_1, set_mo_occupation_2
50 END INTERFACE
51
52CONTAINS
53
54! **************************************************************************************************
55!> \brief Occupation for smeared spin polarized electronic structures
56!> with relaxed multiplicity
57!>
58!> \param mo_array ...
59!> \param smear ...
60!> \date 10.03.2011 (MI)
61!> \author MI
62!> \version 1.0
63! **************************************************************************************************
64 SUBROUTINE set_mo_occupation_3(mo_array, smear)
65
66 TYPE(mo_set_type), DIMENSION(2), INTENT(INOUT) :: mo_array
67 TYPE(smear_type) :: smear
68
69 CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mo_occupation_3'
70
71 INTEGER :: all_nmo, handle, homo_a, homo_b, i, &
72 lfomo_a, lfomo_b, nmo_a, nmo_b, &
73 xas_estate
74 INTEGER, ALLOCATABLE, DIMENSION(:) :: all_index
75 LOGICAL :: is_large
76 REAL(KIND=dp) :: all_nelec, kts, mu, nelec_a, nelec_b, &
77 occ_estate
78 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: all_eigval, all_occ
79 REAL(KIND=dp), DIMENSION(:), POINTER :: eigval_a, eigval_b, occ_a, occ_b
80
81 CALL timeset(routinen, handle)
82
83 NULLIFY (eigval_a, eigval_b, occ_a, occ_b)
84 CALL get_mo_set(mo_set=mo_array(1), nmo=nmo_a, eigenvalues=eigval_a, &
85 occupation_numbers=occ_a)
86 CALL get_mo_set(mo_set=mo_array(2), nmo=nmo_b, eigenvalues=eigval_b, &
87 occupation_numbers=occ_b)
88 all_nmo = nmo_a + nmo_b
89 ALLOCATE (all_eigval(all_nmo))
90 ALLOCATE (all_occ(all_nmo))
91 ALLOCATE (all_index(all_nmo))
92
93 all_eigval(1:nmo_a) = eigval_a(1:nmo_a)
94 all_eigval(nmo_a + 1:all_nmo) = eigval_b(1:nmo_b)
95
96 CALL sort(all_eigval, all_nmo, all_index)
97
98 xas_estate = -1
99 occ_estate = 0.0_dp
100
101 nelec_a = 0.0_dp
102 nelec_b = 0.0_dp
103 all_nelec = 0.0_dp
104 nelec_a = accurate_sum(occ_a(:))
105 nelec_b = accurate_sum(occ_b(:))
106 all_nelec = nelec_a + nelec_b
107
108 DO i = 1, all_nmo
109 IF (all_index(i) <= nmo_a) THEN
110 all_occ(i) = occ_a(all_index(i))
111 ELSE
112 all_occ(i) = occ_b(all_index(i) - nmo_a)
113 END IF
114 END DO
115
116 CALL smearfixed(all_occ, mu, kts, all_eigval, all_nelec, &
117 smear%electronic_temperature, 1._dp, smear%method, xas_estate, occ_estate)
118
119 is_large = abs(maxval(all_occ) - 1.0_dp) > smear%eps_fermi_dirac
120 ! this is not a real problem, but the temperature might be a bit large
121 cpwarn_if(is_large, "Fermi-Dirac smearing includes the first MO")
122
123 is_large = abs(minval(all_occ)) > smear%eps_fermi_dirac
124 IF (is_large) &
125 CALL cp_warn(__location__, &
126 "Fermi-Dirac smearing includes the last MO => "// &
127 "Add more MOs for proper smearing.")
128
129 ! check that the total electron count is accurate
130 is_large = (abs(all_nelec - accurate_sum(all_occ(:))) > smear%eps_fermi_dirac*all_nelec)
131 cpwarn_if(is_large, "Total number of electrons is not accurate")
132
133 DO i = 1, all_nmo
134 IF (all_index(i) <= nmo_a) THEN
135 occ_a(all_index(i)) = all_occ(i)
136 eigval_a(all_index(i)) = all_eigval(i)
137 ELSE
138 occ_b(all_index(i) - nmo_a) = all_occ(i)
139 eigval_b(all_index(i) - nmo_a) = all_eigval(i)
140 END IF
141 END DO
142
143 nelec_a = accurate_sum(occ_a(:))
144 nelec_b = accurate_sum(occ_b(:))
145
146 DO i = 1, nmo_a
147 IF (occ_a(i) < 1.0_dp) THEN
148 lfomo_a = i
149 EXIT
150 END IF
151 END DO
152 DO i = 1, nmo_b
153 IF (occ_b(i) < 1.0_dp) THEN
154 lfomo_b = i
155 EXIT
156 END IF
157 END DO
158 homo_a = lfomo_a - 1
159 DO i = nmo_a, lfomo_a, -1
160 IF (occ_a(i) > smear%eps_fermi_dirac) THEN
161 homo_a = i
162 EXIT
163 END IF
164 END DO
165 homo_b = lfomo_b - 1
166 DO i = nmo_b, lfomo_b, -1
167 IF (occ_b(i) > smear%eps_fermi_dirac) THEN
168 homo_b = i
169 EXIT
170 END IF
171 END DO
172
173 CALL set_mo_set(mo_set=mo_array(1), kts=kts/2.0_dp, mu=mu, n_el_f=nelec_a, &
174 lfomo=lfomo_a, homo=homo_a, uniform_occupation=.false.)
175 CALL set_mo_set(mo_set=mo_array(2), kts=kts/2.0_dp, mu=mu, n_el_f=nelec_b, &
176 lfomo=lfomo_b, homo=homo_b, uniform_occupation=.false.)
177
178 CALL timestop(handle)
179
180 END SUBROUTINE set_mo_occupation_3
181
182! **************************************************************************************************
183!> \brief Prepare an occupation of alpha and beta MOs following an Aufbau
184!> principle, i.e. allowing a change in multiplicity.
185!> \param mo_array ...
186!> \param smear ...
187!> \param eval_deriv ...
188!> \param tot_zeff_corr ...
189!> \param probe ...
190!> \date 25.01.2010 (MK)
191!> \par History
192!> 10.2019 Added functionality to adjust mo occupation if the core
193!> charges are changed via CORE_CORRECTION during surface dipole
194!> calculation. Total number of electrons matches the total core
195!> charges if tot_zeff_corr is non-zero. Not yet implemented for
196!> OT type method. [Soumya Ghosh]
197!> \author Matthias Krack (MK)
198!> \version 1.0
199! **************************************************************************************************
200 SUBROUTINE set_mo_occupation_2(mo_array, smear, eval_deriv, tot_zeff_corr, probe)
201
202 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mo_array
203 TYPE(smear_type) :: smear
204 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eval_deriv
205 REAL(KIND=dp), OPTIONAL :: tot_zeff_corr
206 TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
207 POINTER :: probe
208
209 CHARACTER(LEN=*), PARAMETER :: routinen = 'set_mo_occupation_2'
210
211 INTEGER :: handle, i, lumo_a, lumo_b, &
212 multiplicity_new, multiplicity_old, &
213 nelec
214 REAL(kind=dp) :: nelec_f, threshold
215 REAL(kind=dp), DIMENSION(:), POINTER :: eigval_a, eigval_b
216
217 CALL timeset(routinen, handle)
218
219 ! Fall back for the case that we have only one MO set
220 IF (SIZE(mo_array) == 1) THEN
221 IF (PRESENT(probe)) THEN
222 CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe)
223 ELSE IF (PRESENT(eval_deriv)) THEN
224! Change of MO occupancy to account for CORE_CORRECTION is not yet implemented
225 CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv)
226 ELSE
227 IF (PRESENT(tot_zeff_corr)) THEN
228 CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
229 ELSE
230 CALL set_mo_occupation_1(mo_array(1), smear=smear)
231 END IF
232 END IF
233 CALL timestop(handle)
234 RETURN
235 END IF
236
237 IF (PRESENT(probe)) THEN
238 CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe)
239 CALL set_mo_occupation_1(mo_array(2), smear=smear, probe=probe)
240 END IF
241
242 IF (smear%do_smear) THEN
243 IF (smear%fixed_mag_mom < 0.0_dp) THEN
244 IF (PRESENT(tot_zeff_corr)) THEN
245 CALL cp_warn(__location__, &
246 "CORE_CORRECTION /= 0.0 might cause the cell to charge up "// &
247 "that will lead to application of different background "// &
248 "correction compared to the reference system. "// &
249 "Use FIXED_MAGNETIC_MOMENT >= 0.0 if using SMEAR keyword "// &
250 "to correct the electron density")
251 END IF
252 IF (smear%fixed_mag_mom /= -1.0_dp) THEN
253 cpassert(.NOT. (PRESENT(eval_deriv)))
254 CALL set_mo_occupation_3(mo_array, smear=smear)
255 CALL timestop(handle)
256 RETURN
257 END IF
258 ELSE
259 nelec_f = mo_array(1)%n_el_f + mo_array(2)%n_el_f
260 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
261 mo_array(1)%n_el_f = nelec_f/2.0_dp + smear%fixed_mag_mom/2.0_dp
262 mo_array(2)%n_el_f = nelec_f/2.0_dp - smear%fixed_mag_mom/2.0_dp
263 END IF
264 cpassert(.NOT. (PRESENT(eval_deriv)))
265 IF (PRESENT(tot_zeff_corr)) THEN
266 CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
267 CALL set_mo_occupation_1(mo_array(2), smear=smear, tot_zeff_corr=tot_zeff_corr)
268 ELSE
269 CALL set_mo_occupation_1(mo_array(1), smear=smear)
270 CALL set_mo_occupation_1(mo_array(2), smear=smear)
271 END IF
272 END IF
273 END IF
274
275 IF (.NOT. ((mo_array(1)%flexible_electron_count > 0.0_dp) .AND. &
276 (mo_array(2)%flexible_electron_count > 0.0_dp))) THEN
277 IF (PRESENT(probe)) THEN
278 CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe)
279 CALL set_mo_occupation_1(mo_array(2), smear=smear, probe=probe)
280 ELSE IF (PRESENT(eval_deriv)) THEN
281 CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv)
282 CALL set_mo_occupation_1(mo_array(2), smear=smear, eval_deriv=eval_deriv)
283 ELSE
284 IF (PRESENT(tot_zeff_corr)) THEN
285 CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
286 CALL set_mo_occupation_1(mo_array(2), smear=smear, tot_zeff_corr=tot_zeff_corr)
287 ELSE
288 CALL set_mo_occupation_1(mo_array(1), smear=smear)
289 CALL set_mo_occupation_1(mo_array(2), smear=smear)
290 END IF
291 END IF
292 CALL timestop(handle)
293 RETURN
294 END IF
295
296 nelec = mo_array(1)%nelectron + mo_array(2)%nelectron
297
298 multiplicity_old = mo_array(1)%nelectron - mo_array(2)%nelectron + 1
299
300 IF (mo_array(1)%nelectron >= mo_array(1)%nmo) &
301 CALL cp_warn(__location__, &
302 "All alpha MOs are occupied. Add more alpha MOs to "// &
303 "allow for a higher multiplicity")
304 IF ((mo_array(2)%nelectron >= mo_array(2)%nmo) .AND. (mo_array(2)%nelectron /= mo_array(1)%nelectron)) &
305 CALL cp_warn(__location__, "All beta MOs are occupied. Add more beta MOs to "// &
306 "allow for a lower multiplicity")
307
308 eigval_a => mo_array(1)%eigenvalues
309 eigval_b => mo_array(2)%eigenvalues
310
311 lumo_a = 1
312 lumo_b = 1
313
314 ! Apply Aufbau principle
315 DO i = 1, nelec
316 ! Threshold is needed to ensure a preference for alpha occupation in the case
317 ! of degeneracy
318 threshold = max(mo_array(1)%flexible_electron_count, mo_array(2)%flexible_electron_count)
319 IF ((eigval_a(lumo_a) - threshold) < eigval_b(lumo_b)) THEN
320 lumo_a = lumo_a + 1
321 ELSE
322 lumo_b = lumo_b + 1
323 END IF
324 IF (lumo_a > mo_array(1)%nmo) THEN
325 IF (i /= nelec) &
326 CALL cp_warn(__location__, &
327 "All alpha MOs are occupied. Add more alpha MOs to "// &
328 "allow for a higher multiplicity")
329 IF (i < nelec) THEN
330 lumo_a = lumo_a - 1
331 lumo_b = lumo_b + 1
332 END IF
333 END IF
334 IF (lumo_b > mo_array(2)%nmo) THEN
335 IF (lumo_b < lumo_a) &
336 CALL cp_warn(__location__, &
337 "All beta MOs are occupied. Add more beta MOs to "// &
338 "allow for a lower multiplicity")
339 IF (i < nelec) THEN
340 lumo_a = lumo_a + 1
341 lumo_b = lumo_b - 1
342 END IF
343 END IF
344 END DO
345
346 mo_array(1)%homo = lumo_a - 1
347 mo_array(2)%homo = lumo_b - 1
348
349 IF (mo_array(2)%homo > mo_array(1)%homo) THEN
350 CALL cp_warn(__location__, &
351 "More beta ("// &
352 trim(adjustl(cp_to_string(mo_array(2)%homo)))// &
353 ") than alpha ("// &
354 trim(adjustl(cp_to_string(mo_array(1)%homo)))// &
355 ") MOs are occupied. Resorting to low spin state")
356 mo_array(1)%homo = nelec/2 + modulo(nelec, 2)
357 mo_array(2)%homo = nelec/2
358 END IF
359
360 mo_array(1)%nelectron = mo_array(1)%homo
361 mo_array(2)%nelectron = mo_array(2)%homo
362 multiplicity_new = mo_array(1)%nelectron - mo_array(2)%nelectron + 1
363
364 IF (multiplicity_new /= multiplicity_old) &
365 CALL cp_warn(__location__, &
366 "Multiplicity changed from "// &
367 trim(adjustl(cp_to_string(multiplicity_old)))//" to "// &
368 trim(adjustl(cp_to_string(multiplicity_new))))
369
370 IF (PRESENT(probe)) THEN
371 CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe)
372 CALL set_mo_occupation_1(mo_array(2), smear=smear, probe=probe)
373 ELSE IF (PRESENT(eval_deriv)) THEN
374 CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv)
375 CALL set_mo_occupation_1(mo_array(2), smear=smear, eval_deriv=eval_deriv)
376 ELSE
377 IF (PRESENT(tot_zeff_corr)) THEN
378 CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
379 CALL set_mo_occupation_1(mo_array(2), smear=smear, tot_zeff_corr=tot_zeff_corr)
380 ELSE
381 CALL set_mo_occupation_1(mo_array(1), smear=smear)
382 CALL set_mo_occupation_1(mo_array(2), smear=smear)
383 END IF
384 END IF
385
386 CALL timestop(handle)
387
388 END SUBROUTINE set_mo_occupation_2
389
390! **************************************************************************************************
391!> \brief Smearing of the MO occupation with all kind of occupation numbers
392!> \param mo_set MO dataset structure
393!> \param smear optional smearing information
394!> \param eval_deriv on entry the derivative of the KS energy wrt to the occupation number
395!> on exit the derivative of the full free energy (i.e. KS and entropy) wrt to the eigenvalue
396!> \param xas_env ...
397!> \param tot_zeff_corr ...
398!> \param probe ...
399!> \date 17.04.2002 (v1.0), 26.08.2008 (v1.1)
400!> \par History
401!> 10.2019 Added functionality to adjust mo occupation if the core
402!> charges are changed via CORE_CORRECTION during surface dipole
403!> calculation. Total number of electrons matches the total core
404!> charges if tot_zeff_corr is non-zero. Not yet implemented for
405!> OT type method. [Soumya Ghosh]
406!> \author Matthias Krack
407!> \version 1.1
408! **************************************************************************************************
409 SUBROUTINE set_mo_occupation_1(mo_set, smear, eval_deriv, xas_env, tot_zeff_corr, probe)
410
411 TYPE(mo_set_type), INTENT(INOUT) :: mo_set
412 TYPE(smear_type), OPTIONAL :: smear
413 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eval_deriv
414 TYPE(xas_environment_type), OPTIONAL, POINTER :: xas_env
415 REAL(kind=dp), OPTIONAL :: tot_zeff_corr
416 TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
417 POINTER :: probe
418
419 CHARACTER(LEN=*), PARAMETER :: routinen = 'set_mo_occupation_1'
420
421 CHARACTER(LEN=20) :: method_label
422 INTEGER :: handle, i_first, imo, ir, irmo, nmo, &
423 nomo, xas_estate
424 LOGICAL :: equal_size, is_large
425 REAL(kind=dp) :: delectron, e1, e2, edelta, edist, &
426 el_count, nelec, occ_estate, &
427 total_zeff_corr, xas_nelectron
428 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tmp_v
429
430 CALL timeset(routinen, handle)
431
432 cpassert(ASSOCIATED(mo_set%eigenvalues))
433 cpassert(ASSOCIATED(mo_set%occupation_numbers))
434 mo_set%occupation_numbers(:) = 0.0_dp
435
436 ! Quick return, if no electrons are available
437 IF (mo_set%nelectron == 0) THEN
438 CALL timestop(handle)
439 RETURN
440 END IF
441
442 xas_estate = -1
443 occ_estate = 0.0_dp
444 IF (PRESENT(xas_env)) THEN
445 CALL get_xas_env(xas_env=xas_env, xas_nelectron=xas_nelectron, occ_estate=occ_estate, xas_estate=xas_estate)
446 nomo = ceiling(xas_nelectron + 1.0 - occ_estate - epsilon(0.0_dp))
447
448 mo_set%occupation_numbers(1:nomo) = mo_set%maxocc
449 IF (xas_estate > 0) mo_set%occupation_numbers(xas_estate) = occ_estate
450 el_count = sum(mo_set%occupation_numbers(1:nomo))
451 IF (el_count > xas_nelectron) &
452 mo_set%occupation_numbers(nomo) = mo_set%occupation_numbers(nomo) - (el_count - xas_nelectron)
453 el_count = sum(mo_set%occupation_numbers(1:nomo))
454 is_large = abs(el_count - xas_nelectron) > xas_nelectron*epsilon(el_count)
455 cpassert(.NOT. is_large)
456 ELSE
457 IF (modulo(mo_set%nelectron, int(mo_set%maxocc)) == 0) THEN
458 nomo = nint(mo_set%nelectron/mo_set%maxocc)
459 ! Initialize MO occupations
460 mo_set%occupation_numbers(1:nomo) = mo_set%maxocc
461 ELSE
462 nomo = int(mo_set%nelectron/mo_set%maxocc) + 1
463 ! Initialize MO occupations
464 mo_set%occupation_numbers(1:nomo - 1) = mo_set%maxocc
465 mo_set%occupation_numbers(nomo) = mo_set%nelectron - (nomo - 1)*mo_set%maxocc
466 END IF
467! introduce applied potential correction here
468! electron density is adjusted according to applied core correction
469! ref: SS, MT, MWF, JN PRL, 2018, 120, 246801
470! see whether both surface dipole correction and core correction is present in
471! the inputfile
472 IF (PRESENT(tot_zeff_corr)) THEN
473! find the additional core charges
474 total_zeff_corr = tot_zeff_corr
475 IF (int(mo_set%maxocc) == 1) total_zeff_corr = total_zeff_corr/2.0_dp
476 delectron = 0.0_dp
477 IF (total_zeff_corr < 0.0_dp) THEN
478! remove electron density from the mos
479 delectron = abs(total_zeff_corr) - real(mo_set%maxocc, kind=dp)
480 IF (delectron > 0.0_dp) THEN
481 mo_set%occupation_numbers(nomo) = 0.0_dp
482 irmo = ceiling(delectron/real(mo_set%maxocc, kind=dp))
483 DO ir = 1, irmo
484 delectron = delectron - real(mo_set%maxocc, kind=dp)
485 IF (delectron < 0.0_dp) THEN
486 mo_set%occupation_numbers(nomo - ir) = -delectron
487 ELSE
488 mo_set%occupation_numbers(nomo - ir) = 0.0_dp
489 END IF
490 END DO
491 nomo = nomo - irmo
492 IF (mo_set%occupation_numbers(nomo) == 0.0_dp) nomo = nomo - 1
493 ELSEIF (delectron < 0.0_dp) THEN
494 mo_set%occupation_numbers(nomo) = -delectron
495 ELSE
496 mo_set%occupation_numbers(nomo) = 0.0_dp
497 nomo = nomo - 1
498 END IF
499 ELSEIF (total_zeff_corr > 0.0_dp) THEN
500! add electron density to the mos
501 delectron = total_zeff_corr - real(mo_set%maxocc, kind=dp)
502 IF (delectron > 0.0_dp) THEN
503 mo_set%occupation_numbers(nomo + 1) = real(mo_set%maxocc, kind=dp)
504 nomo = nomo + 1
505 irmo = ceiling(delectron/real(mo_set%maxocc, kind=dp))
506 DO ir = 1, irmo
507 delectron = delectron - real(mo_set%maxocc, kind=dp)
508 IF (delectron < 0.0_dp) THEN
509 mo_set%occupation_numbers(nomo + ir) = delectron + real(mo_set%maxocc, kind=dp)
510 ELSE
511 mo_set%occupation_numbers(nomo + ir) = real(mo_set%maxocc, kind=dp)
512 END IF
513 END DO
514 nomo = nomo + irmo
515 ELSE
516 mo_set%occupation_numbers(nomo + 1) = total_zeff_corr
517 nomo = nomo + 1
518 END IF
519 END IF
520 END IF
521 END IF
522 nmo = SIZE(mo_set%eigenvalues)
523
524 cpassert(nmo >= nomo)
525 cpassert((SIZE(mo_set%occupation_numbers) == nmo))
526
527 mo_set%homo = nomo
528 mo_set%lfomo = nomo + 1
529 mo_set%mu = mo_set%eigenvalues(nomo)
530
531 ! Check consistency of the array lengths
532 IF (PRESENT(eval_deriv)) THEN
533 equal_size = (SIZE(mo_set%occupation_numbers, 1) == SIZE(eval_deriv, 1))
534 cpassert(equal_size)
535 END IF
536
537!calling of HP module HERE, before smear
538 IF (PRESENT(probe)) THEN
539 i_first = 1
540 IF (smear%fixed_mag_mom == -1.0_dp) THEN
541 nelec = real(mo_set%nelectron, dp)
542 ELSE
543 nelec = mo_set%n_el_f
544 END IF
545
546 mo_set%occupation_numbers(:) = 0.0_dp
547
548 CALL probe_occupancy(mo_set%occupation_numbers, mo_set%mu, mo_set%kTS, &
549 mo_set%eigenvalues, mo_set%mo_coeff, mo_set%maxocc, &
550 probe, n=nelec)
551 !NB: mu and T are taken from the hairy_probe type (defined in cp_control_types.F); these values are set in the input
552
553 ! Find the lowest fractional occupied MO (LFOMO)
554 DO imo = i_first, nmo
555 IF (mo_set%occupation_numbers(imo) < mo_set%maxocc) THEN
556 mo_set%lfomo = imo
557 EXIT
558 END IF
559 END DO
560 is_large = abs(maxval(mo_set%occupation_numbers) - mo_set%maxocc) > probe(1)%eps_hp
561 ! this is not a real problem, but the temperature might be a bit large
562 IF (is_large) &
563 cpwarn("Hair-probes occupancy distribution includes the first MO")
564
565 ! Find the highest (fractional) occupied MO which will be now the HOMO
566 DO imo = nmo, mo_set%lfomo, -1
567 IF (mo_set%occupation_numbers(imo) > probe(1)%eps_hp) THEN
568 mo_set%homo = imo
569 EXIT
570 END IF
571 END DO
572 is_large = abs(minval(mo_set%occupation_numbers)) > probe(1)%eps_hp
573 IF (is_large) &
574 CALL cp_warn(__location__, &
575 "Hair-probes occupancy distribution includes the last MO => "// &
576 "Add more MOs for proper smearing.")
577
578 ! check that the total electron count is accurate
579 is_large = (abs(nelec - accurate_sum(mo_set%occupation_numbers(:))) > probe(1)%eps_hp*nelec)
580 IF (is_large) &
581 cpwarn("Total number of electrons is not accurate")
582
583 END IF
584
585 ! Quick return, if no smearing information is supplied (TO BE FIXED, smear should become non-optional...)
586 IF (.NOT. PRESENT(smear)) THEN
587 ! there is no dependence of the energy on the eigenvalues
588 mo_set%uniform_occupation = .true.
589 IF (PRESENT(eval_deriv)) THEN
590 eval_deriv = 0.0_dp
591 END IF
592 CALL timestop(handle)
593 RETURN
594 END IF
595
596 ! Check if proper eigenvalues are already available
597 IF (smear%method /= smear_list) THEN
598 IF ((abs(mo_set%eigenvalues(1)) < 1.0e-12_dp) .AND. &
599 (abs(mo_set%eigenvalues(nmo)) < 1.0e-12_dp)) THEN
600 CALL timestop(handle)
601 RETURN
602 END IF
603 END IF
604
605 ! Perform smearing
606 IF (smear%do_smear) THEN
607 IF (PRESENT(xas_env)) THEN
608 i_first = xas_estate + 1
609 nelec = xas_nelectron
610 ELSE
611 i_first = 1
612 IF (smear%fixed_mag_mom == -1.0_dp) THEN
613 nelec = real(mo_set%nelectron, dp)
614 ELSE
615 nelec = mo_set%n_el_f
616 END IF
617 END IF
618 SELECT CASE (smear%method)
619 CASE (smear_fermi_dirac)
620 IF (.NOT. PRESENT(eval_deriv)) THEN
621 CALL smearfixed(mo_set%occupation_numbers(1:mo_set%nmo), mo_set%mu, mo_set%kTS, &
622 mo_set%eigenvalues(1:mo_set%nmo), nelec, &
623 smear%electronic_temperature, mo_set%maxocc, smear_fermi_dirac, &
624 xas_estate, occ_estate)
625 ELSE
626 IF (.NOT. ALLOCATED(tmp_v)) ALLOCATE (tmp_v(SIZE(eval_deriv)))
627 tmp_v(:) = eval_deriv - mo_set%eigenvalues + mo_set%mu
628 CALL smearfixedderivmv(eval_deriv, mo_set%occupation_numbers(1:mo_set%nmo), mo_set%mu, &
629 mo_set%kTS, mo_set%eigenvalues(1:mo_set%nmo), nelec, &
630 smear%electronic_temperature, mo_set%maxocc, smear_fermi_dirac, &
631 tmp_v, xas_estate, occ_estate)
632 END IF
633
634 ! Find the lowest fractional occupied MO (LFOMO)
635 DO imo = i_first, nmo
636 IF (mo_set%occupation_numbers(imo) < mo_set%maxocc) THEN
637 mo_set%lfomo = imo
638 EXIT
639 END IF
640 END DO
641 is_large = abs(maxval(mo_set%occupation_numbers) - mo_set%maxocc) > smear%eps_fermi_dirac
642 ! this is not a real problem, but the temperature might be a bit large
643 cpwarn_if(is_large, "Fermi-Dirac smearing includes the first MO")
644
645 ! Find the highest (fractional) occupied MO which will be now the HOMO
646 DO imo = nmo, mo_set%lfomo, -1
647 IF (mo_set%occupation_numbers(imo) > smear%eps_fermi_dirac) THEN
648 mo_set%homo = imo
649 EXIT
650 END IF
651 END DO
652 is_large = abs(minval(mo_set%occupation_numbers)) > smear%eps_fermi_dirac
653 IF (is_large) &
654 CALL cp_warn(__location__, &
655 "Fermi-Dirac smearing includes the last MO => "// &
656 "Add more MOs for proper smearing.")
657
658 ! check that the total electron count is accurate
659 is_large = (abs(nelec - accurate_sum(mo_set%occupation_numbers(:))) > smear%eps_fermi_dirac*nelec)
660 cpwarn_if(is_large, "Total number of electrons is not accurate")
661
662 CASE (smear_gaussian, smear_mp, smear_mv)
663 IF (.NOT. PRESENT(eval_deriv)) THEN
664 CALL smearfixed(mo_set%occupation_numbers(1:mo_set%nmo), mo_set%mu, mo_set%kTS, &
665 mo_set%eigenvalues(1:mo_set%nmo), nelec, &
666 smear%smearing_width, mo_set%maxocc, smear%method, &
667 xas_estate, occ_estate)
668 ELSE
669 IF (.NOT. ALLOCATED(tmp_v)) ALLOCATE (tmp_v(SIZE(eval_deriv)))
670 tmp_v(:) = eval_deriv - mo_set%eigenvalues + mo_set%mu
671 CALL smearfixedderivmv(eval_deriv, mo_set%occupation_numbers(1:mo_set%nmo), mo_set%mu, &
672 mo_set%kTS, mo_set%eigenvalues(1:mo_set%nmo), nelec, &
673 smear%smearing_width, mo_set%maxocc, smear%method, &
674 tmp_v, xas_estate, occ_estate)
675 END IF
676
677 ! Method label for warnings
678 SELECT CASE (smear%method)
679 CASE (smear_gaussian)
680 method_label = "Gaussian"
681 CASE (smear_mp)
682 method_label = "Methfessel-Paxton"
683 CASE (smear_mv)
684 method_label = "Marzari-Vanderbilt"
685 END SELECT
686
687 ! Find the lowest fractional occupied MO (LFOMO)
688 DO imo = i_first, nmo
689 IF (abs(mo_set%occupation_numbers(imo) - mo_set%maxocc) > smear%eps_fermi_dirac) THEN
690 mo_set%lfomo = imo
691 EXIT
692 END IF
693 END DO
694 is_large = abs(mo_set%occupation_numbers(1) - mo_set%maxocc) > smear%eps_fermi_dirac
695 cpwarn_if(is_large, trim(method_label)//" smearing includes the first MO")
696
697 ! Find the highest (fractional) occupied MO which will be now the HOMO
698 DO imo = nmo, mo_set%lfomo, -1
699 IF (abs(mo_set%occupation_numbers(imo)) > smear%eps_fermi_dirac) THEN
700 mo_set%homo = imo
701 EXIT
702 END IF
703 END DO
704 is_large = abs(mo_set%occupation_numbers(nmo)) > smear%eps_fermi_dirac
705 IF (is_large) &
706 CALL cp_warn(__location__, &
707 trim(method_label)//" smearing includes the last MO => "// &
708 "Add more MOs for proper smearing.")
709
710 ! Check that the total electron count is accurate
711 is_large = (abs(nelec - accurate_sum(mo_set%occupation_numbers(:))) > smear%eps_fermi_dirac*nelec)
712 cpwarn_if(is_large, "Total number of electrons is not accurate")
713
714 CASE (smear_energy_window)
715 ! not implemented
716 cpassert(.NOT. PRESENT(eval_deriv))
717
718 ! Define the energy window for the eigenvalues
719 e1 = mo_set%eigenvalues(mo_set%homo) - 0.5_dp*smear%window_size
720 IF (e1 <= mo_set%eigenvalues(1)) THEN
721 cpwarn("Energy window for smearing includes the first MO")
722 END IF
723
724 e2 = mo_set%eigenvalues(mo_set%homo) + 0.5_dp*smear%window_size
725 IF (e2 >= mo_set%eigenvalues(nmo)) &
726 CALL cp_warn(__location__, &
727 "Energy window for smearing includes the last MO => "// &
728 "Add more MOs for proper smearing.")
729
730 ! Find the lowest fractional occupied MO (LFOMO)
731 DO imo = i_first, nomo
732 IF (mo_set%eigenvalues(imo) > e1) THEN
733 mo_set%lfomo = imo
734 EXIT
735 END IF
736 END DO
737
738 ! Find the highest fractional occupied (non-zero) MO which will be the HOMO
739 DO imo = nmo, nomo, -1
740 IF (mo_set%eigenvalues(imo) < e2) THEN
741 mo_set%homo = imo
742 EXIT
743 END IF
744 END DO
745
746 ! Get the number of electrons to be smeared
747 edist = 0.0_dp
748 nelec = 0.0_dp
749
750 DO imo = mo_set%lfomo, mo_set%homo
751 nelec = nelec + mo_set%occupation_numbers(imo)
752 edist = edist + abs(e2 - mo_set%eigenvalues(imo))
753 END DO
754
755 ! Smear electrons inside the energy window
756 DO imo = mo_set%lfomo, mo_set%homo
757 edelta = abs(e2 - mo_set%eigenvalues(imo))
758 mo_set%occupation_numbers(imo) = min(mo_set%maxocc, nelec*edelta/edist)
759 nelec = nelec - mo_set%occupation_numbers(imo)
760 edist = edist - edelta
761 END DO
762
763 CASE (smear_list)
764 equal_size = SIZE(mo_set%occupation_numbers, 1) == SIZE(smear%list, 1)
765 cpassert(equal_size)
766 mo_set%occupation_numbers = smear%list
767 ! there is no dependence of the energy on the eigenvalues
768 IF (PRESENT(eval_deriv)) THEN
769 eval_deriv = 0.0_dp
770 END IF
771 ! most general case
772 mo_set%lfomo = 1
773 mo_set%homo = nmo
774 END SELECT
775
776 ! Check, if the smearing involves more than one MO
777 IF (mo_set%lfomo == mo_set%homo) THEN
778 mo_set%homo = nomo
779 mo_set%lfomo = nomo + 1
780 ELSE
781 mo_set%uniform_occupation = .false.
782 END IF
783
784 END IF ! do smear
785
786 ! zeros don't count as uniform
787 mo_set%uniform_occupation = has_uniform_occupation(mo_set=mo_set)
788
789 IF (ALLOCATED(tmp_v)) DEALLOCATE (tmp_v)
790 CALL timestop(handle)
791
792 END SUBROUTINE set_mo_occupation_1
793
794END MODULE qs_mo_occupation
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
kahan_sum::accurate_sum
Definition kahan_sum.F:41
qs_mo_occupation::set_mo_occupation
Definition qs_mo_occupation.F:48
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
hairy_probes
Definition hairy_probes.F:8
hairy_probes::probe_occupancy
subroutine, public probe_occupancy(occ, fermi, kts, energies, coeff, maxocc, probe, n)
subroutine to calculate occupation number and 'Fermi' level using the
Definition hairy_probes.F:52
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::smear_fermi_dirac
integer, parameter, public smear_fermi_dirac
Definition input_constants.F:910
input_constants::smear_energy_window
integer, parameter, public smear_energy_window
Definition input_constants.F:910
input_constants::smear_list
integer, parameter, public smear_list
Definition input_constants.F:910
input_constants::smear_gaussian
integer, parameter, public smear_gaussian
Definition input_constants.F:910
input_constants::smear_mv
integer, parameter, public smear_mv
Definition input_constants.F:910
input_constants::smear_mp
integer, parameter, public smear_mp
Definition input_constants.F:910
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
qs_mo_occupation
Set occupation of molecular orbitals.
Definition qs_mo_occupation.F:15
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::set_mo_set
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.
Definition qs_mo_types.F:454
qs_mo_types::has_uniform_occupation
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...
Definition qs_mo_types.F:505
qs_mo_types::get_mo_set
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.
Definition qs_mo_types.F:399
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
smearing_utils
Unified smearing module supporting four methods: smear_fermi_dirac — Fermi-Dirac distribution smear_g...
Definition smearing_utils.F:24
smearing_utils::smearfixed
subroutine, public smearfixed(f, mu, kts, e, n, sigma, maxocc, method, estate, festate)
Bisection search for the chemical potential mu such that the total electron count equals N,...
Definition smearing_utils.F:312
smearing_utils::smearfixedderivmv
subroutine, public smearfixedderivmv(result, f, mu, kts, e, n_el, sigma, maxocc, method, v, estate, festate)
Apply TRANSPOSE(df/de) to a vector WITHOUT forming the full N x N Jacobian. O(N) time and O(N) memory...
Definition smearing_utils.F:862
util
All kind of helpful little routines.
Definition util.F:14
xas_env_types
define create destroy get and put information in xas_env to calculate the x-ray absorption spectra
Definition xas_env_types.F:15
xas_env_types::get_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)
...
Definition xas_env_types.F:154
cp_control_types::hairy_probes_type
Definition cp_control_types.F:45
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
scf_control_types::smear_type
contains the parameters needed by a scf run
Definition scf_control_types.F:91
xas_env_types::xas_environment_type
Definition xas_env_types.F:69