(git:e5b1968)
Loading...
Searching...
No Matches
hairy_probes.F
Go to the documentation of this file.
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
8MODULE 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
36CONTAINS
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 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 ncol_global=ncol_global)
77 ALLOCATE (smatrix(nrow_global, ncol_global))
78 CALL cp_fm_get_submatrix(coeff, smatrix)
79
80 ALLOCATE (smatrix_squared(nrow_global, ncol_global))
81 smatrix_squared(:, :) = smatrix(:, :)**2.0d0
82
83!**************************************************************************************
84!1)calculate Fermi energy using the hairy probe formula
85!**************************************************************************************
86 de = probe(1)%T*log((1.0_dp - epsocc)/epsocc)
87 de = max(de, 0.5_dp)
88 fermi_max = maxval(probe%mu) + de
89 fermi_min = minval(probe%mu) - de
90
91 CALL hp_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, &
92 maxocc=maxocc, fermi=fermi_max, occupancy=occ, kts=kts)
93 n_max = accurate_sum(occ)
94
95 CALL hp_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, &
96 maxocc=maxocc, fermi=fermi_min, occupancy=occ, kts=kts)
97 n_min = accurate_sum(occ)
98
99 iter = 0
100 DO WHILE (abs(n_max - n_min) > n*epsocc)
101 iter = iter + 1
102
103 fermi_half = (fermi_max + fermi_min)/2.0_dp
104 CALL hp_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, &
105 maxocc=maxocc, fermi=fermi_half, occupancy=occ, kts=kts)
106 n_half = accurate_sum(occ)
107
108 h = fermi_half - fermi_min
109 IF (h .GT. n*epsocc*100) THEN
110 y0 = n_min - n
111 y1 = n_half - n
112 y2 = n_max - n
113
114 CALL three_point_zero(y0, y1, y2, h, delta_fermi)
115 fermi_fit = fermi_min + delta_fermi
116
117 CALL hp_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, &
118 maxocc=maxocc, fermi=fermi_fit, occupancy=occ, kts=kts)
119 n_fit = accurate_sum(occ)
120 END IF
121
122 !define 1st bracked using fermi_half
123 IF (n_half < n) THEN
124 fermi_min = fermi_half
125 n_min = n_half
126 ELSE IF (n_half > n) THEN
127 fermi_max = fermi_half
128 n_max = n_half
129 ELSE
130 fermi_min = fermi_half
131 n_min = n_half
132 fermi_max = fermi_half
133 n_max = n_half
134 h = 0.0d0
135 END IF
136
137 !define 2nd bracker using fermi_fit
138 IF (h .GT. n*epsocc*100) THEN
139 IF (fermi_fit .GE. fermi_min .AND. fermi_fit .LE. fermi_max) THEN
140 IF (n_fit < n) THEN
141 fermi_min = fermi_fit
142 n_min = n_fit
143 ELSE IF (n_fit > n) THEN
144 fermi_max = fermi_fit
145 n_max = n_fit
146 END IF
147 END IF
148 END IF
149
150 IF (abs(n_max - n) < n*epsocc) THEN
151 fermi = fermi_max
152 EXIT
153 ELSE IF (abs(n_min - n) < n*epsocc) THEN
154 fermi = fermi_min
155 EXIT
156 END IF
157
158 IF (iter > bisect_max_iter) THEN
159 cpwarn("Maximum number of iterations reached while finding the Fermi energy")
160 EXIT
161 END IF
162 END DO
163
164!**************************************************************************************
165!2)calculate occupation numbers according to hairy probe formula
166!**************************************************************************************
167 occ(:) = 0.0_dp
168 n_now = 0.0_dp
169 CALL hp_occupancy(probe=probe, matrix=smatrix_squared, energies=energies, &
170 maxocc=maxocc, fermi=fermi, occupancy=occ, kts=kts)
171 n_now = accurate_sum(occ)
172
173 IF (abs(n_now - n) > n*epsocc) cpwarn("Total number of electrons is not accurate - HP")
174
175 DEALLOCATE (smatrix, smatrix_squared)
176
177 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 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 REAL(kind=dp), ALLOCATABLE :: coeff_squared(:, :, :, :)
211
212 CALL timeset(routinen, handle)
213
214!**************************************************************************************
215!1)calculate Fermi energy using the hairy probe formula
216!**************************************************************************************
217 nao = SIZE(rcoeff, 1)
218 nmo = SIZE(rcoeff, 2)
219 nkp = SIZE(rcoeff, 3)
220 nspin = SIZE(rcoeff, 4)
221
222 ALLOCATE (coeff_squared(nao, nmo, nkp, nspin))
223 coeff_squared(:, :, :, :) = rcoeff(:, :, :, :)**2.0d0 + icoeff(:, :, :, :)**2.0d0
224
225 occ(:, :, :) = 0.0_dp
226
227 !define initial brackets
228 de = probe(1)%T*log((1.0_dp - epsocc)/epsocc)
229 de = max(de, 0.5_dp)
230 fermi_max = maxval(probe%mu) + de
231 fermi_min = minval(probe%mu) - de
232
233 n_max = 0.0_dp
234 kts = 0.0_dp
235 !***HP loop
236 DO ispin = 1, nspin
237 DO ikp = 1, nkp
238 CALL hp_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), &
239 maxocc, fermi_max, occ(:, ikp, ispin), kts_kp)
240
241 kts = kts + kts_kp*wk(ikp) !entropic contribution
242 n_max = n_max + accurate_sum(occ(1:nmo, ikp, ispin))*wk(ikp)
243 END DO
244 END DO
245 !***HP loop
246
247 n_min = 0.0_dp
248 kts = 0.0_dp
249 !***HP loop
250 DO ispin = 1, nspin
251 DO ikp = 1, nkp
252 CALL hp_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), &
253 maxocc, fermi_min, occ(:, ikp, ispin), kts_kp)
254
255 kts = kts + kts_kp*wk(ikp) !entropic contribution
256 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 DO WHILE (abs(n_max - n_min) > n*epsocc)
263 iter = iter + 1
264 fermi_half = (fermi_max + fermi_min)/2.0_dp
265 n_half = 0.0_dp
266 kts = 0.0_dp
267 !***HP loop
268 DO ispin = 1, nspin
269 DO ikp = 1, nkp
270 CALL hp_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), &
271 maxocc, fermi_half, occ(:, ikp, ispin), kts_kp)
272
273 kts = kts + kts_kp*wk(ikp) !entropic contribution
274 n_half = n_half + accurate_sum(occ(1:nmo, ikp, ispin))*wk(ikp)
275 END DO
276 END DO
277 !***HP loop
278
279 h = fermi_half - fermi_min
280 IF (h .GT. n*epsocc*100) THEN
281 y0 = n_min - n
282 y1 = n_half - n
283 y2 = n_max - n
284
285 CALL three_point_zero(y0, y1, y2, h, delta_fermi)
286 fermi_fit = fermi_min + delta_fermi
287 n_fit = 0.0_dp
288 kts = 0.0_dp
289
290 !***HP loop
291 DO ispin = 1, nspin
292 DO ikp = 1, nkp
293 CALL hp_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), &
294 maxocc, fermi_fit, occ(:, ikp, ispin), kts_kp)
295
296 kts = kts + kts_kp*wk(ikp) !entropic contribution
297 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 IF (n_half < n) THEN
306 fermi_min = fermi_half
307 n_min = n_half
308 ELSE IF (n_half > n) THEN
309 fermi_max = fermi_half
310 n_max = n_half
311 ELSE
312 fermi_min = fermi_half
313 n_min = n_half
314 fermi_max = fermi_half
315 n_max = n_half
316 h = 0.0d0
317 END IF
318
319 !define 2nd bracker using fermi_fit
320 IF (h .GT. n*epsocc*100) THEN
321 IF (fermi_fit .GE. fermi_min .AND. fermi_fit .LE. fermi_max) THEN
322 IF (n_fit < n) THEN
323 fermi_min = fermi_fit
324 n_min = n_fit
325 ELSE IF (n_fit > n) THEN
326 fermi_max = fermi_fit
327 n_max = n_fit
328 END IF
329 END IF
330 END IF
331
332 IF (abs(n_max - n) < n*epsocc) THEN
333 fermi = fermi_max
334 EXIT
335 ELSE IF (abs(n_min - n) < n*epsocc) THEN
336 fermi = fermi_min
337 EXIT
338 END IF
339
340 IF (iter > bisect_max_iter) THEN
341 cpwarn("Maximum number of iterations reached while finding the Fermi energy")
342 EXIT
343 END IF
344 END DO
345
346!**************************************************************************************
347!3)calculate occupation numbers using the hairy probe formula
348!**************************************************************************************
349 n_now = 0.0_dp
350 kts = 0.0_dp
351
352 !***HP loop
353 DO ispin = 1, nspin
354 DO ikp = 1, nkp
355 CALL hp_occupancy(probe, coeff_squared(:, :, ikp, ispin), energies(:, ikp, ispin), &
356 maxocc, fermi, occ(:, ikp, ispin), kts_kp)
357
358 kts = kts + kts_kp*wk(ikp) !entropic contribution
359 n_now = n_now + accurate_sum(occ(1:nmo, ikp, ispin))*wk(ikp)
360 END DO
361 END DO
362 !***HP loop
363
364 DEALLOCATE (coeff_squared)
365
366 CALL timestop(handle)
367 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 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 nmo = SIZE(matrix, 2)
398 np = SIZE(probe)
399 kts = 0.0_dp
400 alpha = 1.0_dp
401 ! kTS is the entropic contribution to the electronic energy
402
403 mos2: DO imo = 1, nmo
404
405 sum_fermi_fun = 0.0_dp
406 sum_coeff = 0.0_dp
407
408 sum_fermi_fun_sol = 0.0_dp
409 sum_coeff_sol = 0.0_dp
410 fermi_fun_sol = 0.0_dp
411
412 probes2: DO ip = 1, np
413 IF (probe(ip)%alpha .LT. 1.0_dp) THEN
414 alpha = probe(ip)%alpha
415 !fermi distribution, solution probes
416 CALL fermi_distribution(energies(imo), fermi, probe(ip)%T, fermi_fun_sol)
417 !sum of coefficients, solution probes
418 sum_coeff_sol = sum_coeff_sol + sum(matrix(probe(ip)%first_ao:probe(ip)%last_ao, imo))
419 ELSE
420 c = sum(matrix(probe(ip)%first_ao:probe(ip)%last_ao, imo))
421 !bias probes
422 mu = fermi - probe(ip)%mu
423 !fermi distribution, main probes
424 CALL fermi_distribution(energies(imo), mu, probe(ip)%T, fermi_fun)
425 !sum fermi distribution * coefficients
426 sum_fermi_fun = sum_fermi_fun + (fermi_fun*c)
427 !sum cofficients
428 sum_coeff = sum_coeff + c
429 END IF
430 END DO probes2
431
432 sum_fermi_fun_sol = alpha*fermi_fun_sol*sum_coeff_sol
433 sum_coeff_sol = alpha*sum_coeff_sol
434 f = (sum_fermi_fun_sol + sum_fermi_fun)/(sum_coeff_sol + sum_coeff)
435 occupancy(imo) = f*maxocc
436
437 !entropy kTS= kT*[f ln f + (1-f) ln (1-f)]
438 IF (f .EQ. 0.0d0 .OR. f .EQ. 1.0d0) THEN
439 s = 0.0d0
440 ELSE
441 s = f*log(f) + (1.0d0 - f)*log(1.0d0 - f)
442 END IF
443 kts = kts + probe(np)%T*maxocc*s
444 END DO mos2
445
446 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 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 INTEGER, DIMENSION(:), POINTER :: nshell
472 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 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
486 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, nsgf=nsgf) !indexes for orbital symbols
487
488 !get iAO boundaries:
489 p_atom = SIZE(probe%atom_ids)
490 probe%first_ao = nsgf !nAO
491 probe%last_ao = 1
492 isgf = 0
493
494 atoms: DO iatom = 1, natom
495 NULLIFY (orb_basis_set)
496 !1) iatom is used to find the correct atom kind and its index (ikind) is the particle set
497 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 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 IF (ASSOCIATED(orb_basis_set)) THEN
504 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
505 nset=nset, nshell=nshell, l=l) !? ,cgf_symbol=bcgf_symbol )
506 END IF
507
508 !4) get iAO boundaries
509 sets: DO iset = 1, nset
510
511 shells: DO ishell = 1, nshell(iset)
512 lshell = l(ishell, iset)
513
514 isgf = isgf + nso(lshell)
515
516 boundaries: DO ii = 1, p_atom
517 IF (iatom .NE. probe%atom_ids(ii)) THEN
518 cycle boundaries
519 ELSE
520 !defines iAO boundaries***********
521 probe%first_ao = min(probe%first_ao, isgf)
522 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 IF (isgf .NE. nao) cpwarn("row count does not correspond to nAO, number of rows in mo_coeff")
531
532 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 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 IF (e > pot) THEN
554 arg = -(e - pot)/temp
555 ! exponential is smaller than 1
556 exponential = exp(arg)
557 exponential_plus_1 = exponential + 1.0_dp
558
559 one_minus_f = exponential/exponential_plus_1
560 f = 1.0_dp/exponential_plus_1
561 f_p = one_minus_f
562 ELSE
563 arg = (e - pot)/temp
564 ! exponential is smaller than 1
565 exponential = exp(arg)
566 exponential_plus_1 = exponential + 1.0_dp
567
568 f = 1.0_dp/exponential_plus_1
569 one_minus_f = exponential/exponential_plus_1
570 f_p = f
571 END IF
572
573 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 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 a = (y2 - 2.0_dp*y1 + y0)/(2.0_dp*h*h)
593 b = (4.0_dp*y1 - 3.0_dp*y0 - y2)/(2.0_dp*h)
594 c = y0
595
596 IF (abs(a) .LT. 1.0e-15_dp) THEN
597 delta = 0.0_dp
598 RETURN
599 END IF
600
601 d = sqrt(b*b - 4.0_dp*a*c)
602
603 u0 = (-b + d)/(2.0_dp*a)
604
605 IF (u0 .GE. 0.0_dp .AND. u0 .LE. 2.0_dp*h) THEN
606 delta = u0
607 !RETURN
608 ELSE
609 u0 = (-b - d)/(2.0_dp*a)
610 IF (u0 .GE. 0.0_dp .AND. u0 .LE. 2.0_dp*h) THEN
611 delta = u0
612 !RETURN
613 ELSE
614 IF (y1 .LT. 0.0_dp) delta = 2.0_dp*h
615 IF (y1 .GE. 0.0_dp) delta = 0.0_dp
616 END IF
617 END IF
618
619 END SUBROUTINE three_point_zero
620
621END MODULE hairy_probes
kahan_sum::accurate_sum
Definition kahan_sum.F:41
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Definition basis_set_types.F:624
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1009
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:894
hairy_probes
Definition hairy_probes.F:8
hairy_probes::ao_boundaries
subroutine, public ao_boundaries(probe, atomic_kind_set, qs_kind_set, particle_set, nao)
...
Definition hairy_probes.F:462
hairy_probes::probe_occupancy_kp
subroutine, public probe_occupancy_kp(occ, fermi, kts, energies, rcoeff, icoeff, maxocc, probe, n, wk)
subroutine to calculate occupation number and 'Fermi' level using the
Definition hairy_probes.F:194
hairy_probes::probe_occupancy
subroutine, public probe_occupancy(occ, fermi, kts, energies, coeff, maxocc, probe, n)
subroutine to calculate occupation number and 'Fermi' level using the
Definition hairy_probes.F:52
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:879
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_control_types::hairy_probes_type
Definition cp_control_types.F:41
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
particle_types::particle_type
Definition particle_types.F:35
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171