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 : !!***
|
