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