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