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