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