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 : MODULE hairy_probes
9 :
10 : USE atomic_kind_types, ONLY: atomic_kind_type,&
11 : get_atomic_kind,&
12 : get_atomic_kind_set
13 : USE basis_set_types, ONLY: get_gto_basis_set,&
14 : gto_basis_set_type
15 : USE cp_control_types, ONLY: hairy_probes_type
16 : USE cp_fm_types, ONLY: cp_fm_get_info,&
17 : cp_fm_get_submatrix,&
18 : cp_fm_type
19 : USE kahan_sum, ONLY: accurate_sum
20 : USE kinds, ONLY: dp
21 : USE orbital_pointers, ONLY: nso
22 : USE particle_types, ONLY: particle_type
23 : USE qs_kind_types, ONLY: get_qs_kind,&
24 : get_qs_kind_set,&
25 : qs_kind_type
26 : #include "./base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 : PRIVATE
30 :
31 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hairy_probes'
32 :
33 : INTEGER, PARAMETER, PRIVATE :: BISECT_MAX_ITER = 400
34 : PUBLIC :: probe_occupancy, probe_occupancy_kp, AO_boundaries
35 :
36 : CONTAINS
37 :
38 : !**************************************************************************************
39 : !> \brief subroutine to calculate occupation number and 'Fermi' level using the
40 : !> \brief HAIR PROBE approach; gamma point calculation.
41 : !> \param occ occupation numbers
42 : !> \param fermi fermi level
43 : !> \param kTS entropic energy contribution
44 : !> \param energies MOs eigenvalues
45 : !> \param coeff MOs coefficient
46 : !> \param maxocc maximum allowed occupation number of an MO (1 or 2)
47 : !> \param probe hairy probe
48 : !> \param N number of electrons
49 : ! **************************************************************************************************
50 :
51 14 : SUBROUTINE probe_occupancy(occ, fermi, kTS, energies, coeff, maxocc, probe, N)
52 :
53 : !i/o variables and arrays
54 : REAL(KIND=dp), INTENT(out) :: occ(:), fermi, kTS
55 : REAL(KIND=dp), INTENT(IN) :: energies(:)
56 : TYPE(cp_fm_type), INTENT(IN), POINTER :: coeff
57 : REAL(KIND=dp), INTENT(IN) :: maxocc
58 : TYPE(hairy_probes_type), INTENT(INOUT) :: probe(:)
59 : REAL(KIND=dp), INTENT(IN) :: N
60 :
61 : REAL(KIND=dp), PARAMETER :: epsocc = 1.0e-12_dp
62 :
63 : INTEGER :: iter, ncol_global, nrow_global
64 : REAL(KIND=dp) :: de, delta_fermi, fermi_fit, fermi_half, &
65 : fermi_max, fermi_min, h, N_fit, &
66 : N_half, N_max, N_min, N_now, y0, y1, y2
67 : REAL(KIND=dp), ALLOCATABLE :: smatrix_squared(:, :)
68 : REAL(KIND=dp), POINTER :: smatrix(:, :)
69 :
70 : !subroutine variables and arrays
71 : !smatrix: Spherical MOs matrix
72 : !squared Spherical MOs matrix
73 :
74 : CALL cp_fm_get_info(coeff, &
75 : nrow_global=nrow_global, &
76 14 : ncol_global=ncol_global)
77 56 : ALLOCATE (smatrix(nrow_global, ncol_global))
78 14 : CALL cp_fm_get_submatrix(coeff, smatrix)
79 :
80 42 : ALLOCATE (smatrix_squared(nrow_global, ncol_global))
81 5894 : smatrix_squared(:, :) = smatrix(:, :)**2.0d0
82 :
83 : !**************************************************************************************
84 : !1)calculate Fermi energy using the hairy probe formula
85 : !**************************************************************************************
86 14 : de = probe(1)%T*LOG((1.0_dp - epsocc)/epsocc)
87 14 : de = MAX(de, 0.5_dp)
88 56 : fermi_max = MAXVAL(probe%mu) + de
89 56 : fermi_min = MINVAL(probe%mu) - de
90 :
91 : CALL HP_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, &
92 14 : maxocc=maxocc, fermi=fermi_max, occupancy=occ, kTS=kTS)
93 14 : N_max = accurate_sum(occ)
94 :
95 : CALL HP_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, &
96 14 : maxocc=maxocc, fermi=fermi_min, occupancy=occ, kTS=kTS)
97 14 : N_min = accurate_sum(occ)
98 :
99 14 : iter = 0
100 124 : DO WHILE (ABS(N_max - N_min) > N*epsocc)
101 124 : iter = iter + 1
102 :
103 124 : fermi_half = (fermi_max + fermi_min)/2.0_dp
104 : CALL HP_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, &
105 124 : maxocc=maxocc, fermi=fermi_half, occupancy=occ, kTS=kTS)
106 124 : N_half = accurate_sum(occ)
107 :
108 124 : h = fermi_half - fermi_min
109 124 : IF (h > N*epsocc*100) THEN
110 124 : y0 = N_min - N
111 124 : y1 = N_half - N
112 124 : y2 = N_max - N
113 :
114 124 : CALL three_point_zero(y0, y1, y2, h, delta_fermi)
115 124 : fermi_fit = fermi_min + delta_fermi
116 :
117 : CALL HP_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, &
118 124 : maxocc=maxocc, fermi=fermi_fit, occupancy=occ, kTS=kTS)
119 124 : N_fit = accurate_sum(occ)
120 : END IF
121 :
122 : !define 1st bracked using fermi_half
123 124 : IF (N_half < N) THEN
124 52 : fermi_min = fermi_half
125 52 : N_min = N_half
126 72 : ELSE IF (N_half > N) THEN
127 72 : fermi_max = fermi_half
128 72 : N_max = N_half
129 : ELSE
130 0 : fermi_min = fermi_half
131 0 : N_min = N_half
132 0 : fermi_max = fermi_half
133 0 : N_max = N_half
134 0 : h = 0.0d0
135 : END IF
136 :
137 : !define 2nd bracker using fermi_fit
138 124 : IF (h > N*epsocc*100) THEN
139 124 : IF (fermi_fit >= fermi_min .AND. fermi_fit <= fermi_max) THEN
140 124 : IF (N_fit < N) THEN
141 68 : fermi_min = fermi_fit
142 68 : N_min = N_fit
143 56 : ELSE IF (N_fit > N) THEN
144 56 : fermi_max = fermi_fit
145 56 : N_max = N_fit
146 : END IF
147 : END IF
148 : END IF
149 :
150 124 : IF (ABS(N_max - N) < N*epsocc) THEN
151 6 : fermi = fermi_max
152 6 : EXIT
153 118 : ELSE IF (ABS(N_min - N) < N*epsocc) THEN
154 8 : fermi = fermi_min
155 8 : EXIT
156 : END IF
157 :
158 110 : IF (iter > BISECT_MAX_ITER) THEN
159 0 : CPWARN("Maximum number of iterations reached while finding the Fermi energy")
160 0 : EXIT
161 : END IF
162 : END DO
163 :
164 : !**************************************************************************************
165 : !2)calculate occupation numbers according to hairy probe formula
166 : !**************************************************************************************
167 294 : occ(:) = 0.0_dp
168 14 : N_now = 0.0_dp
169 : CALL HP_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, &
170 14 : maxocc=maxocc, fermi=fermi, occupancy=occ, kTS=kTS)
171 14 : N_now = accurate_sum(occ)
172 :
173 14 : IF (ABS(N_now - N) > N*epsocc) CPWARN("Total number of electrons is not accurate - HP")
174 :
175 14 : DEALLOCATE (smatrix, smatrix_squared)
176 :
177 14 : END SUBROUTINE probe_occupancy
178 :
179 : !**************************************************************************************
180 : !> \brief subroutine to calculate occupation number and 'Fermi' level using the
181 : !> \brief HAIR PROBE approach; kpoints calculation.
182 : !> \param occ occupation numbers
183 : !> \param fermi fermi level
184 : !> \param kTS entropic energy contribution
185 : !> \param energies eigenvalues
186 : !> \param rcoeff ...
187 : !> \param icoeff ...
188 : !> \param maxocc maximum allowed occupation number of an MO (1 or 2)
189 : !> \param probe hairy probe
190 : !> \param N number of electrons
191 : !> \param wk weight of kpoints
192 : ! **************************************************************************************************
193 0 : SUBROUTINE probe_occupancy_kp(occ, fermi, kTS, energies, rcoeff, icoeff, maxocc, probe, N, wk)
194 :
195 : REAL(KIND=dp), INTENT(OUT) :: occ(:, :, :), fermi, kTS
196 : REAL(KIND=dp), INTENT(IN) :: energies(:, :, :), rcoeff(:, :, :, :), &
197 : icoeff(:, :, :, :), maxocc
198 : TYPE(hairy_probes_type), INTENT(IN) :: probe(:)
199 : REAL(KIND=dp), INTENT(IN) :: N, wk(:)
200 :
201 : CHARACTER(LEN=*), PARAMETER :: routineN = 'probe_occupancy_kp'
202 : REAL(KIND=dp), PARAMETER :: epsocc = 1.0e-12_dp
203 :
204 : INTEGER :: handle, ikp, ispin, iter, nAO, nkp, nMO, &
205 : nspin
206 : REAL(KIND=dp) :: de, delta_fermi, fermi_fit, fermi_half, &
207 : fermi_max, fermi_min, h, kTS_kp, &
208 : N_fit, N_half, N_max, N_min, N_now, &
209 : y0, y1, y2
210 0 : REAL(KIND=dp), ALLOCATABLE :: coeff_squared(:, :, :, :)
211 :
212 0 : CALL timeset(routineN, handle)
213 :
214 : !**************************************************************************************
215 : !1)calculate Fermi energy using the hairy probe formula
216 : !**************************************************************************************
217 0 : nAO = SIZE(rcoeff, 1)
218 0 : nMO = SIZE(rcoeff, 2)
219 0 : nkp = SIZE(rcoeff, 3)
220 0 : nspin = SIZE(rcoeff, 4)
221 :
222 0 : ALLOCATE (coeff_squared(nAO, nMO, nkp, nspin))
223 0 : coeff_squared(:, :, :, :) = rcoeff(:, :, :, :)**2.0d0 + icoeff(:, :, :, :)**2.0d0
224 :
225 0 : occ(:, :, :) = 0.0_dp
226 :
227 : !define initial brackets
228 0 : de = probe(1)%T*LOG((1.0_dp - epsocc)/epsocc)
229 0 : de = MAX(de, 0.5_dp)
230 0 : fermi_max = MAXVAL(probe%mu) + de
231 0 : fermi_min = MINVAL(probe%mu) - de
232 :
233 0 : N_max = 0.0_dp
234 0 : kTS = 0.0_dp
235 : !***HP loop
236 0 : DO ispin = 1, nspin
237 0 : DO ikp = 1, nkp
238 : CALL HP_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), &
239 0 : maxocc, fermi_max, occ(:, ikp, ispin), kTS_kp)
240 :
241 0 : kTS = kTS + kTS_kp*wk(ikp) !entropic contribution
242 0 : N_max = N_max + accurate_sum(occ(1:nMo, ikp, ispin))*wk(ikp)
243 : END DO
244 : END DO
245 : !***HP loop
246 :
247 0 : N_min = 0.0_dp
248 0 : kTS = 0.0_dp
249 : !***HP loop
250 0 : DO ispin = 1, nspin
251 0 : DO ikp = 1, nkp
252 : CALL HP_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), &
253 0 : maxocc, fermi_min, occ(:, ikp, ispin), kTS_kp)
254 :
255 0 : kTS = kTS + kTS_kp*wk(ikp) !entropic contribution
256 0 : N_min = N_min + accurate_sum(occ(1:nMo, ikp, ispin))*wk(ikp)
257 : END DO
258 : END DO
259 : !***HP loop
260 :
261 : iter = 0
262 0 : DO WHILE (ABS(N_max - N_min) > N*epsocc)
263 0 : iter = iter + 1
264 0 : fermi_half = (fermi_max + fermi_min)/2.0_dp
265 0 : N_half = 0.0_dp
266 0 : kTS = 0.0_dp
267 : !***HP loop
268 0 : DO ispin = 1, nspin
269 0 : DO ikp = 1, nkp
270 : CALL HP_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), &
271 0 : maxocc, fermi_half, occ(:, ikp, ispin), kTS_kp)
272 :
273 0 : kTS = kTS + kTS_kp*wk(ikp) !entropic contribution
274 0 : N_half = N_half + accurate_sum(occ(1:nMo, ikp, ispin))*wk(ikp)
275 : END DO
276 : END DO
277 : !***HP loop
278 :
279 0 : h = fermi_half - fermi_min
280 0 : IF (h > N*epsocc*100) THEN
281 0 : y0 = N_min - N
282 0 : y1 = N_half - N
283 0 : y2 = N_max - N
284 :
285 0 : CALL three_point_zero(y0, y1, y2, h, delta_fermi)
286 0 : fermi_fit = fermi_min + delta_fermi
287 0 : N_fit = 0.0_dp
288 0 : kTS = 0.0_dp
289 :
290 : !***HP loop
291 0 : DO ispin = 1, nspin
292 0 : DO ikp = 1, nkp
293 : CALL HP_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), &
294 0 : maxocc, fermi_fit, occ(:, ikp, ispin), kTS_kp)
295 :
296 0 : kTS = kTS + kTS_kp*wk(ikp) !entropic contribution
297 0 : N_fit = N_fit + accurate_sum(occ(1:nMo, ikp, ispin))*wk(ikp)
298 : END DO
299 : END DO
300 : !***HP loop
301 :
302 : END IF
303 :
304 : !define 1st bracked using fermi_half
305 0 : IF (N_half < N) THEN
306 0 : fermi_min = fermi_half
307 0 : N_min = N_half
308 0 : ELSE IF (N_half > N) THEN
309 0 : fermi_max = fermi_half
310 0 : N_max = N_half
311 : ELSE
312 0 : fermi_min = fermi_half
313 0 : N_min = N_half
314 0 : fermi_max = fermi_half
315 0 : N_max = N_half
316 0 : h = 0.0d0
317 : END IF
318 :
319 : !define 2nd bracker using fermi_fit
320 0 : IF (h > N*epsocc*100) THEN
321 0 : IF (fermi_fit >= fermi_min .AND. fermi_fit <= fermi_max) THEN
322 0 : IF (N_fit < N) THEN
323 0 : fermi_min = fermi_fit
324 0 : N_min = N_fit
325 0 : ELSE IF (N_fit > N) THEN
326 0 : fermi_max = fermi_fit
327 0 : N_max = N_fit
328 : END IF
329 : END IF
330 : END IF
331 :
332 0 : IF (ABS(N_max - N) < N*epsocc) THEN
333 0 : fermi = fermi_max
334 0 : EXIT
335 0 : ELSE IF (ABS(N_min - N) < N*epsocc) THEN
336 0 : fermi = fermi_min
337 0 : EXIT
338 : END IF
339 :
340 0 : IF (iter > BISECT_MAX_ITER) THEN
341 0 : CPWARN("Maximum number of iterations reached while finding the Fermi energy")
342 0 : EXIT
343 : END IF
344 : END DO
345 :
346 : !**************************************************************************************
347 : !3)calculate occupation numbers using the hairy probe formula
348 : !**************************************************************************************
349 0 : N_now = 0.0_dp
350 0 : kTS = 0.0_dp
351 :
352 : !***HP loop
353 0 : DO ispin = 1, nspin
354 0 : DO ikp = 1, nkp
355 : CALL HP_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), &
356 0 : maxocc, fermi, occ(:, ikp, ispin), kTS_kp)
357 :
358 0 : kTS = kTS + kTS_kp*wk(ikp) !entropic contribution
359 0 : N_now = N_now + accurate_sum(occ(1:nMo, ikp, ispin))*wk(ikp)
360 : END DO
361 : END DO
362 : !***HP loop
363 :
364 0 : DEALLOCATE (coeff_squared)
365 :
366 0 : CALL timestop(handle)
367 0 : END SUBROUTINE probe_occupancy_kp
368 :
369 : !**************************************************************************************
370 : !
371 : !
372 : !
373 : !**************************************************************************************
374 : ! **************************************************************************************************
375 : !> \brief ...
376 : !> \param probe ...
377 : !> \param matrix ...
378 : !> \param energies ...
379 : !> \param maxocc ...
380 : !> \param fermi ...
381 : !> \param occupancy ...
382 : !> \param kTS ...
383 : ! **************************************************************************************************
384 290 : SUBROUTINE HP_occupancy(probe, matrix, energies, maxocc, fermi, occupancy, kTS)
385 :
386 : TYPE(hairy_probes_type), INTENT(IN) :: probe(:)
387 : REAL(KIND=dp), INTENT(IN) :: matrix(:, :), energies(:), maxocc, fermi
388 : REAL(KIND=dp), INTENT(OUT) :: occupancy(:), kTS
389 :
390 : INTEGER :: iMO, ip, nMO, np
391 : REAL(KIND=dp) :: alpha, C, f, fermi_fun, fermi_fun_sol, &
392 : mu, s, sum_coeff, sum_coeff_sol, &
393 : sum_fermi_fun, sum_fermi_fun_sol
394 :
395 : !squared coefficient matrix
396 :
397 290 : nMO = SIZE(matrix, 2)
398 290 : np = SIZE(probe)
399 290 : kTS = 0.0_dp
400 290 : alpha = 1.0_dp
401 : ! kTS is the entropic contribution to the electronic energy
402 :
403 6090 : mos2: DO iMO = 1, nMO
404 :
405 5800 : sum_fermi_fun = 0.0_dp
406 5800 : sum_coeff = 0.0_dp
407 :
408 5800 : sum_fermi_fun_sol = 0.0_dp
409 5800 : sum_coeff_sol = 0.0_dp
410 5800 : fermi_fun_sol = 0.0_dp
411 :
412 17400 : probes2: DO ip = 1, np
413 17400 : IF (probe(ip)%alpha < 1.0_dp) THEN
414 0 : alpha = probe(ip)%alpha
415 : !fermi distribution, solution probes
416 0 : CALL fermi_distribution(energies(iMO), fermi, probe(ip)%T, fermi_fun_sol)
417 : !sum of coefficients, solution probes
418 0 : sum_coeff_sol = sum_coeff_sol + SUM(matrix(probe(ip)%first_ao:probe(ip)%last_ao, iMO))
419 : ELSE
420 127600 : C = SUM(matrix(probe(ip)%first_ao:probe(ip)%last_ao, iMO))
421 : !bias probes
422 11600 : mu = fermi - probe(ip)%mu
423 : !fermi distribution, main probes
424 11600 : CALL fermi_distribution(energies(iMO), mu, probe(ip)%T, fermi_fun)
425 : !sum fermi distribution * coefficients
426 11600 : sum_fermi_fun = sum_fermi_fun + (fermi_fun*C)
427 : !sum cofficients
428 11600 : sum_coeff = sum_coeff + C
429 : END IF
430 : END DO probes2
431 :
432 5800 : sum_fermi_fun_sol = alpha*fermi_fun_sol*sum_coeff_sol
433 5800 : sum_coeff_sol = alpha*sum_coeff_sol
434 5800 : f = (sum_fermi_fun_sol + sum_fermi_fun)/(sum_coeff_sol + sum_coeff)
435 5800 : occupancy(iMO) = f*maxocc
436 :
437 : !entropy kTS= kT*[f ln f + (1-f) ln (1-f)]
438 5800 : IF (f == 0.0d0 .OR. f == 1.0d0) THEN
439 : s = 0.0d0
440 : ELSE
441 1732 : s = f*LOG(f) + (1.0d0 - f)*LOG(1.0d0 - f)
442 : END IF
443 6090 : kTS = kTS + probe(np)%T*maxocc*s
444 : END DO MOs2
445 :
446 290 : END SUBROUTINE HP_occupancy
447 :
448 : !**************************************************************************************
449 : !
450 : !
451 : !
452 : !**************************************************************************************
453 : ! **************************************************************************************************
454 : !> \brief ...
455 : !> \param probe ...
456 : !> \param atomic_kind_set ...
457 : !> \param qs_kind_set ...
458 : !> \param particle_set ...
459 : !> \param nAO ...
460 : ! **************************************************************************************************
461 16 : SUBROUTINE AO_boundaries(probe, atomic_kind_set, qs_kind_set, particle_set, nAO)
462 :
463 : TYPE(hairy_probes_type), INTENT(INOUT) :: probe
464 : TYPE(atomic_kind_type), INTENT(IN), POINTER :: atomic_kind_set(:)
465 : TYPE(qs_kind_type), INTENT(IN), POINTER :: qs_kind_set(:)
466 : TYPE(particle_type), INTENT(IN), POINTER :: particle_set(:)
467 : INTEGER, INTENT(IN) :: nAO
468 :
469 : INTEGER :: iatom, ii, ikind, iset, isgf, ishell, &
470 : lshell, natom, nset, nsgf, p_atom
471 8 : INTEGER, DIMENSION(:), POINTER :: nshell
472 8 : INTEGER, DIMENSION(:, :), POINTER :: l
473 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
474 :
475 : !from get_atomic_kind_set
476 : !from get_atomic_kind
477 : !from get_qs_kind_set
478 : !subroutine variables and arrays
479 : !from get_qs_kind
480 : !from get_gto_basis_set
481 : !from get_gto_basis_set
482 : !from get_gto_basis_set
483 :
484 : !get sets from different modules:
485 8 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
486 8 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, nsgf=nsgf) !indexes for orbital symbols
487 :
488 : !get iAO boundaries:
489 8 : p_atom = SIZE(probe%atom_ids)
490 8 : probe%first_ao = nsgf !nAO
491 8 : probe%last_ao = 1
492 8 : isgf = 0
493 :
494 24 : atoms: DO iatom = 1, natom
495 16 : NULLIFY (orb_basis_set)
496 : !1) iatom is used to find the correct atom kind and its index (ikind) is the particle set
497 16 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
498 :
499 : !2) ikind is used to find the basis set associate to that atomic kind
500 16 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
501 :
502 : !3) orb_basis_set is used to get the gto basis set variables
503 16 : IF (ASSOCIATED(orb_basis_set)) THEN
504 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
505 16 : nset=nset, nshell=nshell, l=l) !? ,cgf_symbol=bcgf_symbol )
506 : END IF
507 :
508 : !4) get iAO boundaries
509 56 : sets: DO iset = 1, nset
510 :
511 96 : shells: DO ishell = 1, nshell(iset)
512 64 : lshell = l(ishell, iset)
513 :
514 64 : isgf = isgf + nso(lshell)
515 :
516 144 : boundaries: DO ii = 1, p_atom
517 128 : IF (iatom /= probe%atom_ids(ii)) THEN
518 : CYCLE boundaries
519 : ELSE
520 : !defines iAO boundaries***********
521 32 : probe%first_ao = MIN(probe%first_ao, isgf)
522 64 : probe%last_ao = MAX(probe%last_ao, isgf)
523 : END IF
524 : END DO boundaries
525 :
526 : END DO shells
527 : END DO sets
528 : END DO atoms
529 :
530 8 : IF (isgf /= nAO) CPWARN("row count does not correspond to nAO, number of rows in mo_coeff")
531 :
532 8 : END SUBROUTINE AO_boundaries
533 :
534 : !*******************************************************************************************************
535 : !*******************************************************************************************************
536 :
537 : ! **************************************************************************************************
538 : !> \brief ...
539 : !> \param E ...
540 : !> \param pot ...
541 : !> \param temp ...
542 : !> \param f_p ...
543 : ! **************************************************************************************************
544 11600 : SUBROUTINE fermi_distribution(E, pot, temp, f_p)
545 :
546 : REAL(kind=dp), INTENT(IN) :: E, pot, temp
547 : REAL(kind=dp), INTENT(OUT) :: f_p
548 :
549 : REAL(KIND=dp) :: arg, exponential, exponential_plus_1, f, &
550 : one_minus_f
551 :
552 : ! have the result of exp go to zero instead of overflowing
553 11600 : IF (E > pot) THEN
554 1542 : arg = -(E - pot)/temp
555 : ! exponential is smaller than 1
556 1542 : exponential = EXP(arg)
557 1542 : exponential_plus_1 = exponential + 1.0_dp
558 :
559 1542 : one_minus_f = exponential/exponential_plus_1
560 1542 : f = 1.0_dp/exponential_plus_1
561 1542 : f_p = one_minus_f
562 : ELSE
563 10058 : arg = (E - pot)/temp
564 : ! exponential is smaller than 1
565 10058 : exponential = EXP(arg)
566 10058 : exponential_plus_1 = exponential + 1.0_dp
567 :
568 10058 : f = 1.0_dp/exponential_plus_1
569 10058 : one_minus_f = exponential/exponential_plus_1
570 10058 : f_p = f
571 : END IF
572 :
573 11600 : END SUBROUTINE fermi_distribution
574 :
575 : !*******************************************************************************************************
576 : !*******************************************************************************************************
577 : ! **************************************************************************************************
578 : !> \brief ...
579 : !> \param y0 ...
580 : !> \param y1 ...
581 : !> \param y2 ...
582 : !> \param h ...
583 : !> \param delta ...
584 : ! **************************************************************************************************
585 124 : SUBROUTINE three_point_zero(y0, y1, y2, h, delta)
586 :
587 : REAL(kind=dp), INTENT(IN) :: y0, y1, y2, h
588 : REAL(kind=dp), INTENT(OUT) :: delta
589 :
590 : REAL(kind=dp) :: a, b, c, d, u0
591 :
592 124 : a = (y2 - 2.0_dp*y1 + y0)/(2.0_dp*h*h)
593 124 : b = (4.0_dp*y1 - 3.0_dp*y0 - y2)/(2.0_dp*h)
594 124 : c = y0
595 :
596 124 : IF (ABS(a) < 1.0e-15_dp) THEN
597 0 : delta = 0.0_dp
598 0 : RETURN
599 : END IF
600 :
601 124 : d = SQRT(b*b - 4.0_dp*a*c)
602 :
603 124 : u0 = (-b + d)/(2.0_dp*a)
604 :
605 124 : IF (u0 >= 0.0_dp .AND. u0 <= 2.0_dp*h) THEN
606 124 : delta = u0
607 : !RETURN
608 : ELSE
609 0 : u0 = (-b - d)/(2.0_dp*a)
610 0 : IF (u0 >= 0.0_dp .AND. u0 <= 2.0_dp*h) THEN
611 0 : delta = u0
612 : !RETURN
613 : ELSE
614 0 : IF (y1 < 0.0_dp) delta = 2.0_dp*h
615 0 : IF (y1 >= 0.0_dp) delta = 0.0_dp
616 : END IF
617 : END IF
618 :
619 : END SUBROUTINE three_point_zero
620 :
621 : END MODULE hairy_probes
|
