Branch data Line data Source code
1 : : !!****m*ABINIT/m_scdm_math
2 : : !! NAME
3 : : !! m_scdm_math
4 : : !!
5 : : !! FUNCTION
6 : : !! Math functions used by SCDM
7 : : !! (select columns of density matrix method) for getting wannier functions.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2005-2026 ABINIT group (hexu)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : :
22 : :
23 : : #include "abi_common.h"
24 : :
25 : : module m_scdm_math
26 : : use defs_basis
27 : : use m_abicore
28 : : use m_errors
29 : : implicit none
30 : : real(dp), parameter, public:: tpi = 2*PI
31 : : complex(dp), parameter, public:: tpi_im = cmplx(0.0_dp, tpi, kind = dp)
32 : :
33 : :
34 : : ! wrapper for LAPACK routine ZGEEV
35 : : type, public:: eigensolver
36 : : integer:: ndim = -1, lwork = 0
37 : : complex(dp), allocatable :: work(:)
38 : : real(dp), allocatable:: rwork(:)
39 : : contains
40 : : procedure:: run => eigensolver_run
41 : : procedure:: finalize => eigensolver_finalize
42 : : end type eigensolver
43 : :
44 : :
45 : : public:: complex_QRCP_piv_only
46 : : public:: real_svd
47 : : public:: complex_svd
48 : : public:: gaussian
49 : : public:: fermi
50 : : public:: insertion_sort_double
51 : : public:: build_Rgrid
52 : : private
53 : :
54 : : contains
55 : :
56 : : !
57 [ - + - + ]: 2 : subroutine complex_QRCP_Piv_only(A, Piv)
58 : : complex(dp), intent(in):: A(:, :)
59 : : integer, intent(inout):: Piv(:)
60 : 4 : real(dp):: rwork(size(A, 2)*2)
61 : 4 : complex(dp):: tau(min(size(A, 1), size(A, 2)))
62 : : integer:: m, n
63 : 2 : complex(dp), allocatable:: work(:)
64 : : complex(dp):: tmp_work(2)
65 : : integer:: lwork
66 : : integer:: info
67 : : EXTERNAL ZGEQP3
68 : :
69 : 2 : m = size(A, 1)
70 : 2 : n = size(A, 2)
71 [ + + ]: 62 : rwork(:)=0.0_dp
72 [ - + - - : 2 : call ZGEQP3(m, n, A, m, piv, tau, tmp_work, -1, rwork, info)
- - - - -
+ - - - -
- + - - -
- - + -
- ]
73 [ - + ]: 2 : if(info /= 0) then
74 : 0 : ABI_ERROR("Error in doing QRCP")
75 : : endif
76 [ + + ]: 32 : Piv(:) = 0
77 [ + + ]: 62 : rwork(:) = 0.0_DP
78 : 2 : lwork = INT(AINT(REAL(tmp_work(1))))
79 : :
80 [ - + + - : 6 : ABI_MALLOC(work, (lwork))
+ - + - -
+ - + ]
81 [ + + ]: 1026 : work(:) = (0.0_DP, 0.0_DP)
82 : :
83 [ - + - - : 2 : CALL ZGEQP3(m, n, A, m, piv, tau, work, lwork, rwork, info)
- - - - -
+ - - - -
- + - - -
- - + -
- ]
84 [ - + ]: 2 : if(info /= 0) then
85 : 0 : ABI_ERROR("Error in doing QRCP")
86 : : endif
87 [ + - ]: 2 : ABI_SFREE(work)
88 : 2 : end subroutine complex_QRCP_Piv_only
89 : :
90 : :
91 : :
92 : : ! wrapper for LAPACK routine DGESVD
93 : : ! A = U*S*VT
94 : : ! Note: A is overwritten by U
95 : : ! S is the singular values
96 : : ! VT is the transpose of V
97 [ # # # # : 0 : subroutine real_svd(A, U, S, VT)
# # # # ]
98 : : real(dp), intent(in):: A(:, :)
99 : : real(dp), intent(inout):: U(:, :), S(:), VT(:,:)
100 : : integer:: LWMAX
101 : : real(dp):: tmp(5)
102 : 0 : real(dp), allocatable:: WORK(:)
103 : : integer :: M, N
104 : : integer:: LDA, LDU, LDVT
105 : : integer:: INFO, LWORK
106 : : EXTERNAL DGESVD
107 : 0 : M = size(A, 1)
108 : 0 : N = size(A, 2)
109 : 0 : LDA = M
110 : 0 : LDU = M
111 : 0 : LDVT = N
112 : 0 : LWORK = -1
113 : 0 : LWMAX = max(size(A, 1), size(A, 2))*10
114 : : CALL DGESVD( 'All', 'All', M, N, A, LDA, S, U, LDU, VT, LDVT, &
115 [ # # # # : 0 : & tmp, LWORK, INFO )
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
116 : 0 : LWORK = MIN( LWMAX, INT( tmp( 1 ) ) )
117 : :
118 [ # # # # : 0 : ABI_MALLOC(work, (lwork))
# # # # #
# # # ]
119 : : CALL DGESVD( 'All', 'All', M, N, A, LDA, S, U, LDU, VT, LDVT, &
120 [ # # # # : 0 : WORK, LWORK, INFO )
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
121 [ # # ]: 0 : IF( INFO .GT. 0 ) THEN
122 : 0 : ABI_ERROR('The algorithm computing SVD failed to converge.')
123 : 0 : STOP
124 : : END IF
125 [ # # ]: 0 : ABI_SFREE(work)
126 : 0 : end subroutine real_svd
127 : :
128 : : ! wrapper for LAPACK routine ZGESVD
129 : : ! A = U*S*VT
130 : : ! Note: A is overwritten by U
131 : : ! S is the singular values
132 : : ! VT is the transpose of V
133 : : ! mode: A, or S, or N or O, see
134 : : ! https://netlib.org/lapack/explore-html/d3/da8/group__complex16_g_esing_gad6f0c85f3cca2968e1ef901d2b6014ee.html
135 [ - + - + : 192 : subroutine complex_svd(A, U, S, VT, mode)
- + - + ]
136 : : complex(dp), intent(in):: A(:, :)
137 : : complex(dp), intent(inout):: U(:, :), VT(:,:)
138 : : real(dp), intent(inout):: S(:)
139 : : character, intent(in):: mode ! A, or S, or N/O
140 : : integer:: LWMAX
141 : : complex(dp):: tmp(2)
142 : 192 : complex(dp), allocatable:: WORK(:)
143 : 192 : real(dp), allocatable:: rwork(:)
144 : : integer :: M, N
145 : : integer:: LDA, LDU, LDVT
146 : : integer:: INFO, LWORK
147 : : EXTERNAL ZGESVD
148 : :
149 : 192 : M = size(A, 1)
150 : 192 : N = size(A, 2)
151 : 192 : LWMAX = max(size(A, 1), size(A, 2))*10
152 : 192 : LDA = M
153 : 192 : LDU = M
154 : 192 : LDVT = N
155 : 192 : LWORK = -1
156 [ + - - + : 576 : ABI_MALLOC(rwork, (min(M, N)*6))
+ - + - +
- - + -
+ ]
157 : : CALL ZGESVD( mode, mode, M, N, A, LDA, S, U, LDU, VT, LDVT, &
158 [ - + - - : 192 : & tmp, LWORK, rwork, INFO )
- - - - -
+ - - - -
- + - - -
- - - - +
- - - - -
- - + - -
- - - + -
- - + - -
- - - + -
- - - ]
159 : 192 : LWORK = MIN( LWMAX, INT( tmp( 1 ) ) )
160 [ - + + - : 576 : ABI_MALLOC(work, (lwork))
+ - + - -
+ - + ]
161 : : CALL ZGESVD( mode, mode, M, N, A, LDA, S, U, LDU, VT, LDVT, &
162 [ - + - - : 192 : WORK, LWORK, rwork, INFO )
- - - - -
+ - - - -
- + - - -
- - - - +
- - - - -
- - + - -
- - - + -
- - + - -
- - - + -
- - - ]
163 [ - + ]: 192 : IF( INFO .GT. 0 ) THEN
164 : 0 : ABI_ERROR('The algorithm computing SVD failed to converge.')
165 : 0 : STOP
166 : : END IF
167 [ + - ]: 192 : ABI_SFREE(work)
168 [ + - ]: 192 : ABI_SFREE(rwork)
169 : 192 : end subroutine complex_svd
170 : :
171 : :
172 : 0 : function gaussian(x, mu, sigma) result(y)
173 : : real(dp), intent(in):: x, mu, sigma
174 : : real(dp):: y
175 : 0 : y = exp(-1.0 * (x-mu)**2/sigma**2)
176 : 0 : end function gaussian
177 : :
178 : : !===============================================================
179 : : ! Complementary error function.
180 : : !===============================================================
181 : : pure function erfcd(x) result(y)
182 : : real(dp), intent(in):: x
183 : : real(dp):: y
184 : : real(dp):: t, z
185 : : z = abs(x)
186 : : t = 1.0 / ( 1.0+0.5*z )
187 : :
188 : : y = t*exp( -z*z - 1.26551223+t * &
189 : : ( 1.00002368+t * ( 0.37409196+t * &
190 : : ( 0.09678418+t * (-0.18628806+t * &
191 : : ( 0.27886807+t * (-1.13520398+t * &
192 : : ( 1.48851587+t * (-0.82215223+t * 0.17087277 )))))))))
193 : : if ( x .lt. 0.0 ) y = 2.0-y
194 : :
195 : : end function erfcd
196 : :
197 : 1935 : function fermi(x, mu, sigma) result(y)
198 : : real(dp), intent(in):: x, mu, sigma
199 : : real(dp):: y
200 : 1935 : y = 0.5*erfc((x-mu) / sigma)
201 : 1935 : end function fermi
202 : :
203 : : ! wrapper for LAPACK routine ZHEEV
204 : : ! evecs are used as the input matrix which will be overwritten
205 : : ! evals are the eigenvalues
206 [ - + - + ]: 344 : subroutine eigensolver_run(self, evals, evecs)
207 : : class(eigensolver), intent(inout):: self
208 : : real(dp), intent(inout):: evals(:)
209 : : complex(dp), intent(inout):: evecs(:,:)
210 : : integer :: info
211 : : external ZHEEV
212 [ + - ]: 344 : if (self%ndim == -1) then
213 : 344 : self%ndim = size(evecs, 1)
214 : 344 : self%lwork = -1
215 [ - + - + ]: 344 : ABI_MALLOC(self%work, (1))
216 [ + - - + : 1032 : ABI_MALLOC(self%rwork, (3*size(evecs, 1)-2))
+ - + - +
- - + - +
- + ]
217 [ - + - - : 344 : call ZHEEV('V', 'U', self%ndim, evecs, self%ndim, evals, self%work, self%lwork, self%rwork, info)
- - - - -
+ - - - -
- + - - -
- - + -
- ]
218 : 344 : self%lwork = INT(self%work(1))
219 : 344 : ABI_SFREE(self%work)
220 [ - + + - : 1032 : ABI_MALLOC(self%work, (self%lwork))
+ - + - -
+ - + -
+ ]
221 [ # # ]: 0 : else if (self%ndim /= size(evecs, 1)) then
222 : 0 : ABI_ERROR("Eigensovler: The size of the evecs is not the same as previous one.")
223 : : end if
224 [ - + - - : 344 : call ZHEEV('V', 'U', self%ndim, evecs, self%ndim, evals, self%work, self%lwork, self%rwork, info)
- - - - -
+ - - - -
- + - - -
- - + -
- ]
225 : 344 : end subroutine eigensolver_run
226 : :
227 : : ! free the memory used by eigensolver
228 : 344 : subroutine eigensolver_finalize(self)
229 : : class(eigensolver), intent(inout):: self
230 : 344 : self%ndim = -1
231 : 344 : self%lwork = -1
232 [ + - ]: 344 : ABI_SFREE(self%work)
233 [ + - ]: 344 : ABI_SFREE(self%rwork)
234 : 344 : end subroutine eigensolver_finalize
235 : :
236 : : !----------------------------------------------------------------------
237 : : !> @brief insertion_sort_int: sort a array using insertion sort algorithm
238 : : !> it is a memory safe method but is generally slow.
239 : : !> @param[inout] a: the array to be sorted. and will output inplace
240 : : !> @param[inout] order (optional) the sorted index, it can be used to sort
241 : : !> other arrays so that the order in consistent.
242 : : !----------------------------------------------------------------------
243 [ # # ]: 0 : subroutine insertion_sort_double(a, order)
244 : : real(dp), intent(inout):: a(:)
245 : : integer, optional, intent(inout):: order(size(a))
246 : : integer:: n, i, j
247 : : real(dp):: v
248 : 0 : n = size(a)
249 [ # # ]: 0 : if (present(order)) then
250 [ # # ]: 0 : do i = 1, n
251 : 0 : order(i)=i
252 : : end do
253 : : end if
254 [ # # ]: 0 : do i = 2, n
255 : 0 : v = a(i)
256 : 0 : j = i-1
257 [ # # ]: 0 : do while(j >= 1 )
258 [ # # ]: 0 : if (a(j)<=v) exit
259 : 0 : a(j+1)=a(j)
260 [ # # ]: 0 : if(present(order)) order(j+1)=order(j)
261 : 0 : j = j-1
262 : : end do
263 : 0 : a(j+1)=v
264 [ # # ]: 0 : if(present(order)) order(j+1)=i
265 : : end do
266 : :
267 [ # # ]: 0 : end subroutine insertion_sort_double
268 : :
269 : :
270 : : ! -4/2-> -1, 4/2->2. Then (-1, 0, 1, 2)
271 : 126 : pure function div2(x) result(y)
272 : : integer, intent(in):: x
273 : : integer:: y
274 [ + - - + : 126 : if(mod(x, 2)==0 .and. x < 0) then
+ - - + +
- - + ]
275 : 63 : y = x/2+1
276 : : else
277 : 63 : y = x/2
278 : : endif
279 : : end function div2
280 : :
281 : : ! build the Rlist for the given kmesh
282 : : ! Rlist is a 3*n array, where n is the number of R vectors
283 : 3 : subroutine build_Rgrid(kmesh, Rlist)
284 : : integer, intent(in):: kmesh(3)
285 : : integer, allocatable, intent(inout):: Rlist(:,:)
286 : : integer:: n, i1, i2, i3, i
287 : :
288 : 3 : n = kmesh(1) * kmesh(2) *kmesh(3)
289 [ - + + - : 9 : ABI_MALLOC(Rlist, (3, n))
+ - + - -
+ - + -
+ ]
290 : 3 : i = 0
291 : : ! Note that C/Fortran integer division is "truncate towards 0" division,
292 : : ! whereas Python one is "floor" division.
293 : : ! For C/Fortran, the behavior for even and odd numbers is
294 : : ! not consistent and need special treatment in div2.
295 [ + - + - : 21 : do i3 = div2(-kmesh(3)), div2(kmesh(3))
+ + ]
296 [ + - + - : 87 : do i2 = div2(-kmesh(2)), div2(kmesh(2))
+ + ]
297 [ + - + - : 348 : do i1 = div2(-kmesh(1)), div2(kmesh(1))
+ + ]
298 : 192 : i = i+1
299 [ + + ]: 816 : Rlist(:, i) = [i1, i2, i3]
300 : : end do
301 : : end do
302 : : end do
303 : 3 : end subroutine build_Rgrid
304 : :
305 : :
306 : :
307 [ # # # # : 0 : end module m_scdm_math
# # # # #
# # # #
# ]
308 : :
309 : : !!***
|
