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