Line data Source code
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 :
15 : MODULE 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 :
54 : CONTAINS
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 2552 : 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 2552 : 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 2552 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigval_a, eigval_b, occ_a, occ_b
85 :
86 2552 : CALL timeset(routineN, handle)
87 :
88 2552 : 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 2552 : occupation_numbers=occ_a)
91 : CALL get_mo_set(mo_set=mo_array(2), nmo=nmo_b, eigenvalues=eigval_b, &
92 2552 : occupation_numbers=occ_b)
93 2552 : all_nmo = nmo_a + nmo_b
94 7656 : ALLOCATE (all_eigval(all_nmo))
95 5104 : ALLOCATE (all_occ(all_nmo))
96 7656 : ALLOCATE (all_index(all_nmo))
97 :
98 35148 : all_eigval(1:nmo_a) = eigval_a(1:nmo_a)
99 33216 : all_eigval(nmo_a + 1:all_nmo) = eigval_b(1:nmo_b)
100 :
101 2552 : CALL sort(all_eigval, all_nmo, all_index)
102 :
103 2552 : IF (PRESENT(gce)) THEN
104 64 : do_gce = gce%do_gce
105 : ELSE
106 : do_gce = .FALSE.
107 : END IF
108 :
109 4940 : SELECT CASE (smear%method)
110 : CASE (smear_fermi_dirac)
111 2388 : smear_width = smear%electronic_temperature
112 2388 : method_label = "Fermi-Dirac"
113 : CASE (smear_gaussian)
114 164 : smear_width = smear%smearing_width
115 164 : method_label = "Gaussian"
116 : CASE (smear_mp)
117 0 : smear_width = smear%smearing_width
118 0 : method_label = "Methfessel-Paxton"
119 : CASE (smear_mv)
120 0 : smear_width = smear%smearing_width
121 0 : method_label = "Marzari-Vanderbilt"
122 : CASE DEFAULT
123 2552 : CPABORT("set_mo_occupation_3: unsupported smearing method")
124 : END SELECT
125 :
126 2552 : IF (.NOT. do_gce) THEN
127 2488 : xas_estate = -1
128 2488 : 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 2488 : nelec_a = accurate_sum(occ_a(:))
134 2488 : nelec_b = accurate_sum(occ_b(:))
135 2488 : all_nelec = nelec_a + nelec_b
136 :
137 63316 : DO i = 1, all_nmo
138 63316 : IF (all_index(i) <= nmo_a) THEN
139 31380 : all_occ(i) = occ_a(all_index(i))
140 : ELSE
141 29448 : 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 2488 : smear_width, 1._dp, smear%method, xas_estate, occ_estate)
147 : ELSE
148 64 : gce%prev_workfunction = gce%ref_esp - mo_array(1)%mu
149 64 : mu = gce%ref_esp - ((1.0_dp - gce%mixing_coef)*gce%prev_workfunction + gce%mixing_coef*gce%target_workfunction)
150 64 : CALL SmearOcc(all_occ, all_nelec, kTS, all_eigval, mu, smear_width, 1._dp, smear%method)
151 : END IF
152 :
153 2552 : 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 2552 : CPWARN_IF(is_large, TRIM(method_label)//" smearing includes the first MO")
156 :
157 2552 : is_large = ABS(all_occ(all_nmo)) > smear%eps_fermi_dirac
158 2552 : IF (is_large) THEN
159 : CALL cp_warn(__LOCATION__, &
160 : TRIM(method_label)//" smearing includes the last MO => "// &
161 20 : "Add more MOs for proper smearing.")
162 : END IF
163 2552 : IF (.NOT. do_gce) THEN
164 : ! check that the total electron count is accurate
165 2488 : is_large = (ABS(all_nelec - accurate_sum(all_occ(:))) > smear%eps_fermi_dirac*all_nelec)
166 2488 : CPWARN_IF(is_large, "Total number of electrons is not accurate")
167 : END IF
168 :
169 65812 : DO i = 1, all_nmo
170 65812 : IF (all_index(i) <= nmo_a) THEN
171 32596 : occ_a(all_index(i)) = all_occ(i)
172 32596 : eigval_a(all_index(i)) = all_eigval(i)
173 : ELSE
174 30664 : occ_b(all_index(i) - nmo_a) = all_occ(i)
175 30664 : eigval_b(all_index(i) - nmo_a) = all_eigval(i)
176 : END IF
177 : END DO
178 :
179 2552 : nelec_a = accurate_sum(occ_a(:))
180 2552 : nelec_b = accurate_sum(occ_b(:))
181 :
182 2552 : lfomo_a = nmo_a + 1
183 2552 : lfomo_b = nmo_b + 1
184 : homo_a = 0
185 : homo_b = 0
186 :
187 21236 : DO i = 1, nmo_a
188 21236 : IF (occ_a(i) < 1.0_dp) THEN
189 2552 : lfomo_a = i
190 2552 : EXIT
191 : END IF
192 : END DO
193 20864 : DO i = 1, nmo_b
194 20864 : IF (occ_b(i) < 1.0_dp) THEN
195 2552 : lfomo_b = i
196 2552 : EXIT
197 : END IF
198 : END DO
199 2552 : homo_a = lfomo_a - 1
200 13032 : DO i = nmo_a, lfomo_a, -1
201 13032 : IF (occ_a(i) > smear%eps_fermi_dirac) THEN
202 1876 : homo_a = i
203 1876 : EXIT
204 : END IF
205 : END DO
206 2552 : homo_b = lfomo_b - 1
207 12944 : DO i = nmo_b, lfomo_b, -1
208 12944 : IF (occ_b(i) > smear%eps_fermi_dirac) THEN
209 1134 : homo_b = i
210 1134 : 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 2552 : 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 2552 : lfomo=lfomo_b, homo=homo_b, uniform_occupation=.FALSE.)
218 :
219 2552 : CALL timestop(handle)
220 :
221 5104 : 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 117695 : 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 117695 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigval_a, eigval_b
259 :
260 117695 : CALL timeset(routineN, handle)
261 :
262 : ! Fall back for the case that we have only one MO set
263 117695 : IF (SIZE(mo_array) == 1) THEN
264 100173 : IF (PRESENT(probe)) THEN
265 14 : CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe)
266 100159 : ELSE IF (PRESENT(eval_deriv)) THEN
267 : ! Change of MO occupancy to account for CORE_CORRECTION is not yet implemented
268 0 : IF (PRESENT(gce)) THEN
269 0 : CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv, gce=gce)
270 : ELSE
271 0 : CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv)
272 : END IF
273 : ELSE
274 100159 : IF (PRESENT(tot_zeff_corr)) THEN
275 20 : CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
276 : ELSE
277 100139 : IF (PRESENT(gce)) THEN
278 0 : CALL set_mo_occupation_1(mo_array(1), smear=smear, gce=gce)
279 : ELSE
280 100139 : CALL set_mo_occupation_1(mo_array(1), smear=smear)
281 : END IF
282 : END IF
283 : END IF
284 100173 : CALL timestop(handle)
285 : RETURN
286 : END IF
287 :
288 17522 : IF (PRESENT(probe)) THEN
289 0 : CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe)
290 0 : CALL set_mo_occupation_1(mo_array(2), smear=smear, probe=probe)
291 : END IF
292 :
293 17522 : IF (smear%do_smear) THEN
294 4258 : IF (smear%fixed_mag_mom < 0.0_dp) THEN
295 2906 : 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 0 : "to correct the electron density")
302 : END IF
303 2906 : IF (smear%fixed_mag_mom /= -1.0_dp) THEN
304 2552 : CPASSERT(.NOT. (PRESENT(eval_deriv)))
305 2552 : CALL set_mo_occupation_3(mo_array, smear=smear, gce=gce)
306 2552 : CALL timestop(handle)
307 2552 : RETURN
308 : END IF
309 : ELSE
310 1352 : nelec_f = mo_array(1)%n_el_f + mo_array(2)%n_el_f
311 1352 : 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 2 : mo_array(1)%n_el_f = nelec_f/2.0_dp + smear%fixed_mag_mom/2.0_dp
313 2 : mo_array(2)%n_el_f = nelec_f/2.0_dp - smear%fixed_mag_mom/2.0_dp
314 : END IF
315 1352 : CPASSERT(.NOT. (PRESENT(eval_deriv)))
316 1352 : IF (PRESENT(tot_zeff_corr)) THEN
317 20 : CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
318 20 : CALL set_mo_occupation_1(mo_array(2), smear=smear, tot_zeff_corr=tot_zeff_corr)
319 : ELSE
320 1332 : CALL set_mo_occupation_1(mo_array(1), smear=smear)
321 1332 : CALL set_mo_occupation_1(mo_array(2), smear=smear)
322 : END IF
323 : END IF
324 : END IF
325 :
326 14970 : IF (.NOT. ((mo_array(1)%flexible_electron_count > 0.0_dp) .AND. &
327 : (mo_array(2)%flexible_electron_count > 0.0_dp))) THEN
328 14784 : IF (PRESENT(probe)) THEN
329 0 : CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe)
330 0 : CALL set_mo_occupation_1(mo_array(2), smear=smear, probe=probe)
331 14784 : ELSE IF (PRESENT(eval_deriv)) THEN
332 0 : CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv)
333 0 : CALL set_mo_occupation_1(mo_array(2), smear=smear, eval_deriv=eval_deriv)
334 : ELSE
335 14784 : IF (PRESENT(tot_zeff_corr)) THEN
336 20 : CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
337 20 : CALL set_mo_occupation_1(mo_array(2), smear=smear, tot_zeff_corr=tot_zeff_corr)
338 : ELSE
339 14764 : CALL set_mo_occupation_1(mo_array(1), smear=smear)
340 14764 : CALL set_mo_occupation_1(mo_array(2), smear=smear)
341 : END IF
342 : END IF
343 14784 : CALL timestop(handle)
344 14784 : RETURN
345 : END IF
346 :
347 186 : nelec = mo_array(1)%nelectron + mo_array(2)%nelectron
348 :
349 186 : multiplicity_old = mo_array(1)%nelectron - mo_array(2)%nelectron + 1
350 :
351 186 : 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 0 : "allow for a higher multiplicity")
355 : END IF
356 186 : 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 0 : "allow for a lower multiplicity")
359 : END IF
360 :
361 186 : eigval_a => mo_array(1)%eigenvalues
362 186 : eigval_b => mo_array(2)%eigenvalues
363 :
364 186 : lumo_a = 1
365 186 : lumo_b = 1
366 :
367 : ! Apply Aufbau principle
368 2306 : DO i = 1, nelec
369 : ! Threshold is needed to ensure a preference for alpha occupation in the case
370 : ! of degeneracy
371 2120 : threshold = MAX(mo_array(1)%flexible_electron_count, mo_array(2)%flexible_electron_count)
372 2120 : IF ((eigval_a(lumo_a) - threshold) < eigval_b(lumo_b)) THEN
373 1212 : lumo_a = lumo_a + 1
374 : ELSE
375 908 : lumo_b = lumo_b + 1
376 : END IF
377 2120 : IF (lumo_a > mo_array(1)%nmo) THEN
378 0 : IF (i /= nelec) THEN
379 : CALL cp_warn(__LOCATION__, &
380 : "All alpha MOs are occupied. Add more alpha MOs to "// &
381 0 : "allow for a higher multiplicity")
382 : END IF
383 0 : IF (i < nelec) THEN
384 0 : lumo_a = lumo_a - 1
385 0 : lumo_b = lumo_b + 1
386 : END IF
387 : END IF
388 2306 : IF (lumo_b > mo_array(2)%nmo) THEN
389 34 : IF (lumo_b < lumo_a) THEN
390 : CALL cp_warn(__LOCATION__, &
391 : "All beta MOs are occupied. Add more beta MOs to "// &
392 0 : "allow for a lower multiplicity")
393 : END IF
394 34 : IF (i < nelec) THEN
395 6 : lumo_a = lumo_a + 1
396 6 : lumo_b = lumo_b - 1
397 : END IF
398 : END IF
399 : END DO
400 :
401 186 : mo_array(1)%homo = lumo_a - 1
402 186 : mo_array(2)%homo = lumo_b - 1
403 :
404 186 : 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 0 : ") MOs are occupied. Resorting to low spin state")
411 0 : mo_array(1)%homo = nelec/2 + MODULO(nelec, 2)
412 0 : mo_array(2)%homo = nelec/2
413 : END IF
414 :
415 186 : mo_array(1)%nelectron = mo_array(1)%homo
416 186 : mo_array(2)%nelectron = mo_array(2)%homo
417 186 : multiplicity_new = mo_array(1)%nelectron - mo_array(2)%nelectron + 1
418 :
419 186 : 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 8 : TRIM(ADJUSTL(cp_to_string(multiplicity_new))))
424 : END IF
425 :
426 186 : IF (PRESENT(probe)) THEN
427 0 : CALL set_mo_occupation_1(mo_array(1), smear=smear, probe=probe)
428 0 : CALL set_mo_occupation_1(mo_array(2), smear=smear, probe=probe)
429 186 : ELSE IF (PRESENT(eval_deriv)) THEN
430 0 : CALL set_mo_occupation_1(mo_array(1), smear=smear, eval_deriv=eval_deriv)
431 0 : CALL set_mo_occupation_1(mo_array(2), smear=smear, eval_deriv=eval_deriv)
432 : ELSE
433 186 : IF (PRESENT(tot_zeff_corr)) THEN
434 0 : CALL set_mo_occupation_1(mo_array(1), smear=smear, tot_zeff_corr=tot_zeff_corr)
435 0 : CALL set_mo_occupation_1(mo_array(2), smear=smear, tot_zeff_corr=tot_zeff_corr)
436 : ELSE
437 186 : CALL set_mo_occupation_1(mo_array(1), smear=smear)
438 186 : CALL set_mo_occupation_1(mo_array(2), smear=smear)
439 : END IF
440 : END IF
441 :
442 186 : CALL timestop(handle)
443 :
444 117695 : 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 263344 : 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 263344 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp_v
492 :
493 263344 : CALL timeset(routineN, handle)
494 :
495 263344 : my_emit_warnings = .TRUE.
496 263344 : IF (PRESENT(emit_warnings)) my_emit_warnings = emit_warnings
497 :
498 263344 : CPASSERT(ASSOCIATED(mo_set%eigenvalues))
499 263344 : CPASSERT(ASSOCIATED(mo_set%occupation_numbers))
500 3023805 : mo_set%occupation_numbers(:) = 0.0_dp
501 :
502 : ! Quick return, if no electrons are available
503 263344 : IF (mo_set%nelectron == 0) THEN
504 2528 : CALL timestop(handle)
505 2528 : RETURN
506 : END IF
507 :
508 260816 : xas_estate = -1
509 260816 : occ_estate = 0.0_dp
510 260816 : IF (PRESENT(xas_env)) THEN
511 798 : CALL get_xas_env(xas_env=xas_env, xas_nelectron=xas_nelectron, occ_estate=occ_estate, xas_estate=xas_estate)
512 798 : nomo = CEILING(xas_nelectron + 1.0 - occ_estate - EPSILON(0.0_dp))
513 :
514 8102 : mo_set%occupation_numbers(1:nomo) = mo_set%maxocc
515 798 : IF (xas_estate > 0) mo_set%occupation_numbers(xas_estate) = occ_estate
516 8102 : el_count = SUM(mo_set%occupation_numbers(1:nomo))
517 798 : IF (el_count > xas_nelectron) THEN
518 98 : mo_set%occupation_numbers(nomo) = mo_set%occupation_numbers(nomo) - (el_count - xas_nelectron)
519 : END IF
520 8102 : el_count = SUM(mo_set%occupation_numbers(1:nomo))
521 798 : is_large = ABS(el_count - xas_nelectron) > xas_nelectron*EPSILON(el_count)
522 798 : CPASSERT(.NOT. is_large)
523 : ELSE
524 260018 : IF (PRESENT(gce)) THEN
525 4 : do_gce = gce%do_gce
526 : ELSE
527 : do_gce = .FALSE.
528 : END IF
529 : ! GCE workfunction (Fermi energy) mixing
530 4 : IF (do_gce) THEN
531 4 : IF (smear%method /= smear_fermi_dirac) THEN
532 0 : CPABORT("Grand canonical ensemble DFT SCF now only support Fermi Dirac smearing.")
533 : END IF
534 4 : IF (gce%prev_workfunction < -1000.0_dp) THEN
535 2 : 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 2 : smear%electronic_temperature, mo_set%maxocc, smear_fermi_dirac)
538 2 : gce%prev_workfunction = -501.0_dp
539 2 : ELSE IF (gce%prev_workfunction < -500.0_dp) THEN
540 2 : 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 2 : smear%electronic_temperature, mo_set%maxocc, smear_fermi_dirac)
543 2 : 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 0 : + gce%mixing_coef*gce%target_workfunction)
547 : CALL SmearOcc(mo_set%occupation_numbers, my_nelec, mo_set%kTS, mo_set%eigenvalues, gce_mu, &
548 0 : smear%electronic_temperature, mo_set%maxocc, smear_fermi_dirac)
549 0 : mo_set%mu = gce_mu
550 0 : gce%prev_workfunction = gce%ref_esp - mo_set%mu
551 0 : is_large = ABS(MAXVAL(mo_set%occupation_numbers) - mo_set%maxocc) > smear%eps_fermi_dirac
552 0 : CPWARN_IF(is_large .AND. my_emit_warnings, "Fermi-Dirac smearing includes the first MO")
553 : END IF
554 4 : DO i = 1, SIZE(mo_set%occupation_numbers)
555 4 : IF (mo_set%occupation_numbers(i) < mo_set%maxocc) THEN
556 4 : mo_set%lfomo = i
557 4 : EXIT
558 : END IF
559 : END DO
560 4 : DO i = SIZE(mo_set%occupation_numbers), 1, -1
561 4 : IF (mo_set%occupation_numbers(i) > smear%eps_fermi_dirac) THEN
562 4 : mo_set%homo = i
563 4 : EXIT
564 : END IF
565 : END DO
566 4 : mo_set%uniform_occupation = .FALSE.
567 4 : mo_set%n_el_f = my_nelec
568 4 : CALL timestop(handle)
569 4 : RETURN
570 : END IF
571 :
572 260014 : IF (MODULO(mo_set%nelectron, INT(mo_set%maxocc)) == 0) THEN
573 258940 : nomo = NINT(mo_set%nelectron/mo_set%maxocc)
574 : ! Initialize MO occupations
575 2399149 : mo_set%occupation_numbers(1:nomo) = mo_set%maxocc
576 : ELSE
577 1074 : nomo = INT(mo_set%nelectron/mo_set%maxocc) + 1
578 : ! Initialize MO occupations
579 6608 : mo_set%occupation_numbers(1:nomo - 1) = mo_set%maxocc
580 1074 : 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 260014 : IF (PRESENT(tot_zeff_corr)) THEN
588 : ! find the additional core charges
589 106 : total_zeff_corr = tot_zeff_corr
590 106 : IF (INT(mo_set%maxocc) == 1) total_zeff_corr = total_zeff_corr/2.0_dp
591 106 : delectron = 0.0_dp
592 106 : IF (total_zeff_corr < 0.0_dp) THEN
593 : ! remove electron density from the mos
594 106 : delectron = ABS(total_zeff_corr) - REAL(mo_set%maxocc, KIND=dp)
595 106 : IF (delectron > 0.0_dp) THEN
596 0 : mo_set%occupation_numbers(nomo) = 0.0_dp
597 0 : irmo = CEILING(delectron/REAL(mo_set%maxocc, KIND=dp))
598 0 : DO ir = 1, irmo
599 0 : delectron = delectron - REAL(mo_set%maxocc, KIND=dp)
600 0 : IF (delectron < 0.0_dp) THEN
601 0 : mo_set%occupation_numbers(nomo - ir) = -delectron
602 : ELSE
603 0 : mo_set%occupation_numbers(nomo - ir) = 0.0_dp
604 : END IF
605 : END DO
606 0 : nomo = nomo - irmo
607 0 : IF (mo_set%occupation_numbers(nomo) == 0.0_dp) nomo = nomo - 1
608 106 : ELSE IF (delectron < 0.0_dp) THEN
609 106 : mo_set%occupation_numbers(nomo) = -delectron
610 : ELSE
611 0 : mo_set%occupation_numbers(nomo) = 0.0_dp
612 0 : nomo = nomo - 1
613 : END IF
614 0 : ELSE IF (total_zeff_corr > 0.0_dp) THEN
615 : ! add electron density to the mos
616 0 : delectron = total_zeff_corr - REAL(mo_set%maxocc, KIND=dp)
617 0 : IF (delectron > 0.0_dp) THEN
618 0 : mo_set%occupation_numbers(nomo + 1) = REAL(mo_set%maxocc, KIND=dp)
619 0 : nomo = nomo + 1
620 0 : irmo = CEILING(delectron/REAL(mo_set%maxocc, KIND=dp))
621 0 : DO ir = 1, irmo
622 0 : delectron = delectron - REAL(mo_set%maxocc, KIND=dp)
623 0 : IF (delectron < 0.0_dp) THEN
624 0 : mo_set%occupation_numbers(nomo + ir) = delectron + REAL(mo_set%maxocc, KIND=dp)
625 : ELSE
626 0 : mo_set%occupation_numbers(nomo + ir) = REAL(mo_set%maxocc, KIND=dp)
627 : END IF
628 : END DO
629 0 : nomo = nomo + irmo
630 : ELSE
631 0 : mo_set%occupation_numbers(nomo + 1) = total_zeff_corr
632 0 : nomo = nomo + 1
633 : END IF
634 : END IF
635 : END IF
636 : END IF
637 260812 : nmo = SIZE(mo_set%eigenvalues)
638 :
639 260812 : CPASSERT(nmo >= nomo)
640 260812 : CPASSERT((SIZE(mo_set%occupation_numbers) == nmo))
641 :
642 260812 : mo_set%homo = nomo
643 260812 : mo_set%lfomo = nomo + 1
644 260812 : mo_set%mu = mo_set%eigenvalues(nomo)
645 :
646 : ! Check consistency of the array lengths
647 260812 : IF (PRESENT(eval_deriv)) THEN
648 0 : equal_size = (SIZE(mo_set%occupation_numbers, 1) == SIZE(eval_deriv, 1))
649 0 : CPASSERT(equal_size)
650 : END IF
651 :
652 : !calling of HP module HERE, before smear
653 260812 : IF (PRESENT(probe)) THEN
654 14 : i_first = 1
655 14 : IF (smear%fixed_mag_mom == -1.0_dp) THEN
656 0 : nelec = REAL(mo_set%nelectron, dp)
657 : ELSE
658 14 : nelec = mo_set%n_el_f
659 : END IF
660 :
661 294 : 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 14 : 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 198 : DO imo = i_first, nmo
670 198 : IF (mo_set%occupation_numbers(imo) < mo_set%maxocc) THEN
671 14 : mo_set%lfomo = imo
672 14 : EXIT
673 : END IF
674 : END DO
675 294 : 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 14 : IF (is_large .AND. my_emit_warnings) THEN
678 0 : 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 22 : DO imo = nmo, mo_set%lfomo, -1
683 22 : IF (mo_set%occupation_numbers(imo) > probe(1)%eps_hp) THEN
684 14 : mo_set%homo = imo
685 14 : EXIT
686 : END IF
687 : END DO
688 294 : is_large = ABS(MINVAL(mo_set%occupation_numbers)) > probe(1)%eps_hp
689 14 : IF (is_large .AND. my_emit_warnings) THEN
690 : CALL cp_warn(__LOCATION__, &
691 : "Hair-probes occupancy distribution includes the last MO => "// &
692 6 : "Add more MOs for proper smearing.")
693 : END IF
694 :
695 : ! check that the total electron count is accurate
696 14 : is_large = (ABS(nelec - accurate_sum(mo_set%occupation_numbers(:))) > probe(1)%eps_hp*nelec)
697 14 : IF (is_large .AND. my_emit_warnings) THEN
698 0 : 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 260812 : IF (.NOT. PRESENT(smear)) THEN
705 : ! there is no dependence of the energy on the eigenvalues
706 230 : mo_set%uniform_occupation = .TRUE.
707 230 : IF (PRESENT(eval_deriv)) THEN
708 0 : eval_deriv = 0.0_dp
709 : END IF
710 230 : CALL timestop(handle)
711 230 : RETURN
712 : END IF
713 :
714 : ! Check if proper eigenvalues are already available
715 260582 : IF (smear%method /= smear_list) THEN
716 260558 : IF ((ABS(mo_set%eigenvalues(1)) < 1.0E-12_dp) .AND. &
717 : (ABS(mo_set%eigenvalues(nmo)) < 1.0E-12_dp)) THEN
718 57736 : CALL timestop(handle)
719 57736 : RETURN
720 : END IF
721 : END IF
722 :
723 : ! Perform smearing
724 202846 : IF (smear%do_smear) THEN
725 23832 : IF (PRESENT(xas_env)) THEN
726 30 : i_first = xas_estate + 1
727 30 : nelec = xas_nelectron
728 : ELSE
729 23802 : i_first = 1
730 23802 : IF (smear%fixed_mag_mom == -1.0_dp) THEN
731 782 : nelec = REAL(mo_set%nelectron, dp)
732 : ELSE
733 23020 : nelec = mo_set%n_el_f
734 : END IF
735 : END IF
736 22290 : SELECT CASE (smear%method)
737 : CASE (smear_fermi_dirac)
738 22290 : 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 22290 : xas_estate, occ_estate)
743 : ELSE
744 0 : IF (.NOT. ALLOCATED(tmp_v)) ALLOCATE (tmp_v(SIZE(eval_deriv)))
745 0 : 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 0 : tmp_v, xas_estate, occ_estate)
750 : END IF
751 :
752 : ! Find the lowest fractional occupied MO (LFOMO)
753 216536 : DO imo = i_first, nmo
754 216536 : IF (mo_set%occupation_numbers(imo) < mo_set%maxocc) THEN
755 22290 : mo_set%lfomo = imo
756 22290 : EXIT
757 : END IF
758 : END DO
759 22290 : IF (i_first <= nmo) THEN
760 22290 : 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 22290 : 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 455272 : DO imo = nmo, mo_set%lfomo, -1
769 455272 : IF (mo_set%occupation_numbers(imo) > smear%eps_fermi_dirac) THEN
770 14474 : mo_set%homo = imo
771 14474 : EXIT
772 : END IF
773 : END DO
774 710890 : is_large = ABS(MINVAL(mo_set%occupation_numbers)) > smear%eps_fermi_dirac
775 22290 : IF (is_large .AND. my_emit_warnings) THEN
776 : CALL cp_warn(__LOCATION__, &
777 : "Fermi-Dirac smearing includes the last MO => "// &
778 684 : "Add more MOs for proper smearing.")
779 : END IF
780 :
781 : ! check that the total electron count is accurate
782 22290 : is_large = (ABS(nelec - accurate_sum(mo_set%occupation_numbers(:))) > smear%eps_fermi_dirac*nelec)
783 22290 : 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 1360 : 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 1360 : xas_estate, occ_estate)
791 : ELSE
792 0 : IF (.NOT. ALLOCATED(tmp_v)) ALLOCATE (tmp_v(SIZE(eval_deriv)))
793 0 : 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 0 : tmp_v, xas_estate, occ_estate)
798 : END IF
799 :
800 : ! Method label for warnings
801 2720 : SELECT CASE (smear%method)
802 : CASE (smear_gaussian)
803 1360 : method_label = "Gaussian"
804 : CASE (smear_mp)
805 0 : method_label = "Methfessel-Paxton"
806 : CASE (smear_mv)
807 1360 : method_label = "Marzari-Vanderbilt"
808 : END SELECT
809 :
810 : ! Find the lowest fractional occupied MO (LFOMO)
811 17500 : DO imo = i_first, nmo
812 17500 : IF (ABS(mo_set%occupation_numbers(imo) - mo_set%maxocc) > smear%eps_fermi_dirac) THEN
813 1360 : mo_set%lfomo = imo
814 1360 : EXIT
815 : END IF
816 : END DO
817 1360 : IF (i_first <= nmo) THEN
818 1360 : 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 1360 : IF (is_large .AND. my_emit_warnings) THEN
823 0 : 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 16670 : DO imo = nmo, mo_set%lfomo, -1
828 16670 : IF (ABS(mo_set%occupation_numbers(imo)) > smear%eps_fermi_dirac) THEN
829 928 : mo_set%homo = imo
830 928 : EXIT
831 : END IF
832 : END DO
833 1360 : is_large = ABS(mo_set%occupation_numbers(nmo)) > smear%eps_fermi_dirac
834 1360 : IF (is_large .AND. my_emit_warnings) THEN
835 : CALL cp_warn(__LOCATION__, &
836 : TRIM(method_label)//" smearing includes the last MO => "// &
837 0 : "Add more MOs for proper smearing.")
838 : END IF
839 :
840 : ! Check that the total electron count is accurate
841 1360 : is_large = (ABS(nelec - accurate_sum(mo_set%occupation_numbers(:))) > smear%eps_fermi_dirac*nelec)
842 1360 : 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 158 : CPASSERT(.NOT. PRESENT(eval_deriv))
847 :
848 : ! Define the energy window for the eigenvalues
849 158 : e1 = mo_set%eigenvalues(mo_set%homo) - 0.5_dp*smear%window_size
850 158 : IF (e1 <= mo_set%eigenvalues(1) .AND. my_emit_warnings) THEN
851 0 : CPWARN("Energy window for smearing includes the first MO")
852 : END IF
853 :
854 158 : e2 = mo_set%eigenvalues(mo_set%homo) + 0.5_dp*smear%window_size
855 158 : 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 0 : "Add more MOs for proper smearing.")
859 : END IF
860 :
861 : ! Find the lowest fractional occupied MO (LFOMO)
862 2636 : DO imo = i_first, nomo
863 2636 : IF (mo_set%eigenvalues(imo) > e1) THEN
864 158 : mo_set%lfomo = imo
865 158 : EXIT
866 : END IF
867 : END DO
868 :
869 : ! Find the highest fractional occupied (non-zero) MO which will be the HOMO
870 1344 : DO imo = nmo, nomo, -1
871 1344 : IF (mo_set%eigenvalues(imo) < e2) THEN
872 158 : mo_set%homo = imo
873 158 : EXIT
874 : END IF
875 : END DO
876 :
877 : ! Get the number of electrons to be smeared
878 158 : edist = 0.0_dp
879 158 : nelec = 0.0_dp
880 :
881 390 : DO imo = mo_set%lfomo, mo_set%homo
882 232 : nelec = nelec + mo_set%occupation_numbers(imo)
883 390 : edist = edist + ABS(e2 - mo_set%eigenvalues(imo))
884 : END DO
885 :
886 : ! Smear electrons inside the energy window
887 390 : DO imo = mo_set%lfomo, mo_set%homo
888 232 : edelta = ABS(e2 - mo_set%eigenvalues(imo))
889 232 : mo_set%occupation_numbers(imo) = MIN(mo_set%maxocc, nelec*edelta/edist)
890 232 : nelec = nelec - mo_set%occupation_numbers(imo)
891 390 : edist = edist - edelta
892 : END DO
893 :
894 : CASE (smear_list)
895 24 : equal_size = SIZE(mo_set%occupation_numbers, 1) == SIZE(smear%list, 1)
896 24 : CPASSERT(equal_size)
897 168 : mo_set%occupation_numbers = smear%list
898 : ! there is no dependence of the energy on the eigenvalues
899 24 : IF (PRESENT(eval_deriv)) THEN
900 0 : eval_deriv = 0.0_dp
901 : END IF
902 : ! most general case
903 24 : mo_set%lfomo = 1
904 23856 : mo_set%homo = nmo
905 : END SELECT
906 :
907 : ! Check, if the smearing involves more than one MO
908 23832 : IF (mo_set%lfomo == mo_set%homo) THEN
909 1578 : mo_set%homo = nomo
910 1578 : mo_set%lfomo = nomo + 1
911 : ELSE
912 22254 : mo_set%uniform_occupation = .FALSE.
913 : END IF
914 :
915 : END IF ! do smear
916 :
917 : ! zeros don't count as uniform
918 202846 : mo_set%uniform_occupation = has_uniform_occupation(mo_set=mo_set)
919 :
920 202846 : IF (ALLOCATED(tmp_v)) DEALLOCATE (tmp_v)
921 202846 : CALL timestop(handle)
922 :
923 202846 : END SUBROUTINE set_mo_occupation_1
924 :
925 : END MODULE qs_mo_occupation
|
