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