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