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