Branch data Line data Source code
1 : : !!****m* ABINIT/m_paw_lmn
2 : : !! NAME
3 : : !! m_paw_lmn
4 : : !!
5 : : !! FUNCTION
6 : : !! This module provides tools to calculate tables commonly used to iterate
7 : : !! over the the (l,m,n) channels of the PAW partial waves.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2008-2026 ABINIT group (MG)
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 : : #include "libpaw.h"
18 : :
19 : : MODULE m_paw_lmn
20 : :
21 : : USE_DEFS
22 : : USE_MSG_HANDLING
23 : : USE_MPI_WRAPPERS
24 : : USE_MEMORY_PROFILING
25 : :
26 : : implicit none
27 : :
28 : : private
29 : :
30 : : ! Public procedures.
31 : : public :: ilm2lm ! Returns the value of l and m in from the index ilm
32 : : public :: make_indlmn ! Calculates the indlmn(6,lmn_size) table giving l,m,n,lm,ln,spin for i=lmn.
33 : : public :: make_indklmn ! Calculates the indklmn(8,lmn2_size) table giving
34 : : ! klm, kln, abs(il-jl), (il+jl), ilm and jlm, ilmn and jlmn for each symmetric klmn=(ilmn,jlmn)
35 : : public :: make_kln2ln ! Calculates the kln2ln(6,ln2_size) table giving
36 : : ! il, jl ,in, jn, iln, jln for each symmetric kln=(iln,jln)
37 : : public :: make_klm2lm ! Calculates the klm2lm(6,lm2_size) table giving
38 : : ! il, jl ,im, jm, ilm, jlm for each symmetric klm=(ilm,jlm)
39 : : public :: klmn2ijlmn ! Calculates ilmn and jlmn from klmn.
40 : : public :: make_indln ! Calculates indln(2,ln_size) giving l and n for i=ln
41 : : public :: uppert_index ! The sequential index of an element in the upper triangle of a matrix
42 : :
43 : : CONTAINS !========================================================================================
44 : : !!***
45 : :
46 : : !!****f* m_paw_lmn/ilm2lm
47 : : !! NAME
48 : : !! ilm2lm
49 : : !!
50 : : !! FUNCTION
51 : : !! Returns the value of l and m in from the index ilm
52 : : !!
53 : : !! INPUTS
54 : : !! ilm=The contracted index for (l,m).
55 : : !!
56 : : !! OUTPUT
57 : : !! ll=The angular momentun defined in [0,1,2,3,4].
58 : : !! mm=The magnetic number in the interval [-l,....+l].
59 : : !!
60 : : !! SOURCE
61 : :
62 : 0 : subroutine ilm2lm(ilm,ll,mm)
63 : :
64 : : !Arguments ------------------------------------
65 : : !scalars
66 : : integer,intent(in) :: ilm
67 : : integer,intent(out) :: ll,mm
68 : :
69 : : !Local variables-------------------------------
70 : : integer :: ii
71 : : ! *********************************************************************
72 : :
73 [ # # ]: 0 : if (ilm<1) then
74 : 0 : LIBPAW_ERROR("Wrong ilm")
75 : : end if
76 : :
77 : 0 : ll = -1
78 [ # # ]: 0 : do ii=0,100
79 [ # # ]: 0 : if ( (ii+1)**2 >= ilm) then
80 : 0 : ll = ii
81 : 0 : EXIT
82 : : end if
83 : : end do
84 : :
85 : 0 : mm = ilm - ll**2 -ll-1
86 : :
87 [ # # ]: 0 : if (ll==-1) then
88 : 0 : LIBPAW_ERROR("l>100 not programmed!")
89 : : end if
90 : :
91 : 0 : end subroutine ilm2lm
92 : : !!***
93 : :
94 : : !----------------------------------------------------------------------
95 : :
96 : : !!****f* m_paw_lmn/make_indlmn
97 : : !! NAME
98 : : !! make_indlmn
99 : : !!
100 : : !! FUNCTION
101 : : !! Performs the setup of the indlmn table for PAW calculations (indices for (l,m,n) basis)
102 : : !!
103 : : !! INPUTS
104 : : !! ln_size= Total number of nl components
105 : : !! lmn_size= Second dimension in indlmn. Total number of (l,m,n) components for this atom.
106 : : !! orbitals(ln_size)=Give the value of l for each element of the augmented basis set.
107 : : !!
108 : : !! OUTPUT
109 : : !! indlmn(6,lmn_size)=array giving l,m,n,lm,ln,s for i=lmn
110 : : !! indlmn(8,lmn_size)=array giving l,m,sign of kappa,ln,spinor,2*j,2*m_j in case of dirac relativism
111 : : !!
112 : : !! SOURCE
113 : :
114 : 12 : subroutine make_indlmn(ln_size,lmn_size,orbitals,indlmn,kappa)
115 : :
116 : : !Arguments ------------------------------------
117 : : !scalars
118 : : integer,intent(in) :: ln_size,lmn_size
119 : : integer,intent(in) :: orbitals(ln_size)
120 : : !scalars
121 : : integer,allocatable,intent(inout) :: indlmn(:,:)
122 : : integer, intent(in), optional :: kappa(ln_size)
123 : :
124 : : !Local variables ------------------------------
125 : : !scalars
126 : : integer :: ilmn,ib,il,iln,ilm,kappa_sign,spinor,i2j,i2mj,im,ilmn_ws
127 : : !arrays
128 : 12 : integer,allocatable :: nprj(:)
129 : :
130 : : !************************************************************************
131 : :
132 [ + + ]: 12 : if(.not.present(kappa)) then
133 [ + + - + : 58 : LIBPAW_ALLOCATE(nprj,(0:MAXVAL(orbitals)))
+ - + - +
- - + -
+ ]
134 [ - + + - : 30 : LIBPAW_ALLOCATE(indlmn,(6,lmn_size))
+ - + - -
+ - + -
+ ]
135 [ + + + + ]: 332 : indlmn=0
136 [ + + ]: 28 : ilmn=0; iln=0; nprj=0
137 [ + + ]: 38 : do ib=1,ln_size
138 : 28 : il=orbitals(ib)
139 : 28 : nprj(il)=nprj(il)+1
140 : 28 : iln=iln+1
141 [ + + ]: 74 : do ilm=1,2*il+1
142 : 46 : indlmn(1,ilmn+ilm)=il ! l
143 : 46 : indlmn(2,ilmn+ilm)=ilm-(il+1) ! m
144 : 46 : indlmn(3,ilmn+ilm)=nprj(il) ! n
145 : 46 : indlmn(4,ilmn+ilm)=il*il+ilm ! lm index
146 : 46 : indlmn(5,ilmn+ilm)=iln ! ln index
147 : 74 : indlmn(6,ilmn+ilm)=1 ! spin (not yet used!)
148 : : end do
149 : 38 : ilmn=ilmn+2*il+1
150 : : end do
151 [ - + ]: 10 : LIBPAW_DEALLOCATE(nprj)
152 : : else
153 [ - + + - : 6 : LIBPAW_ALLOCATE(indlmn,(9,lmn_size))
+ - + - -
+ - + -
+ ]
154 [ + + + + ]: 402 : indlmn=0;ilmn=0;iln=0;ilmn_ws=0
155 [ + + ]: 18 : do ib=1,2*ln_size
156 : 16 : iln=iln+modulo(ib,2)
157 : 16 : il=orbitals(iln)
158 : 16 : kappa_sign=sign(1,kappa(iln)) ! sgn(kappa)=+1 or -1
159 : 16 : spinor=2-modulo(ib,2) ! spinor= 1 or 2
160 : 16 : i2j=2*il-kappa_sign ! j=l-sgn(kappa)/2 = l-1/2 or l+1/2
161 [ + + ]: 56 : do ilm=1,i2j+1
162 : 40 : ilmn_ws=ilmn_ws+1
163 : : !mj= -j,...,j
164 : 40 : i2mj=-i2j+2*(ilm-1) ! 2m_j= -jc ... +jc
165 : 40 : im=(i2mj-3+2*spinor)/2 ! m=m_j-1/2 (spinor=1) or m_j+1/2 (spinor=2)
166 [ + + ]: 56 : if(abs(im)<=il) then
167 : : !Valid value for sph. harm., i.e. abs(m)<=l
168 : 28 : indlmn(1,ilmn+ilm)=il !l
169 : 28 : indlmn(2,ilmn+ilm)=im !m
170 : 28 : indlmn(3,ilmn+ilm)=kappa_sign !sign of kappa
171 : 28 : indlmn(4,ilmn+ilm)=il*il+im+il+1 !lm
172 : 28 : indlmn(5,ilmn+ilm)=iln !ln also includes the two spinor values here
173 : 28 : indlmn(6,ilmn+ilm)=spinor !spinor index (1 up, 2 down)
174 : 28 : indlmn(7,ilmn+ilm)=i2j !2*j (times 2 to make it an integer)
175 : 28 : indlmn(8,ilmn+ilm)=i2mj !2*m_j (times 2 to make it an integer)
176 : 28 : indlmn(9,ilmn+ilm)=ilmn_ws ! ilmn without spinor
177 : : else
178 : : !Invalid value for sph. harm. ; will be multiplied by zero later
179 : 12 : indlmn(1,ilmn+ilm)=-1 !Invalid value that should be checked later
180 [ + + ]: 108 : indlmn(2:9,ilmn+ilm)=-1 ; indlmn(3,ilmn+ilm)=0
181 : : endif
182 : : end do
183 : 16 : ilmn=ilmn+i2j+1
184 [ + + ]: 18 : if(modulo(ib,2)==1) ilmn_ws=ilmn_ws-i2j-1
185 : : end do
186 : : endif
187 : :
188 : 12 : end subroutine make_indlmn
189 : : !!***
190 : :
191 : : !----------------------------------------------------------------------
192 : :
193 : : !!****f* m_paw_lmn/make_indklmn
194 : : !! NAME
195 : : !! make_indklmn
196 : : !!
197 : : !! FUNCTION
198 : : !! Performs the setup of the indklmn table for PAW calculations.
199 : : !! Compute the indklmn indexes giving klm, kln, abs(il-jl) and (il+jl), ilm and jlm, ilmn and jlmn
200 : : !! for each klmn=(ilmn,jlmn) with jlmn >= ilmn
201 : : !!
202 : : !! INPUTS
203 : : !! lcutdens=Maximum l for densities/potentials moments computations
204 : : !! lmn_size=Number of (l,m,n) elements for the PAW basis set.
205 : : !! lmn2_size=Number of elements in the symmetric basis set: lmn2_size=lmn_size*(lmn_size+1)/2
206 : : !! indlmn(6,lmn_size)=Array giving l,m,n,lm,ln,spin for i=lmn.
207 : : !!
208 : : !! OUTPUT
209 : : !! indklmn(6,lmn2_size)=Array giving klm, kln, abs(il-jl), (il+jl), ilm and jlm, ilmn and jlmn
210 : : !! for each klmn=(ilmn,jlmn). Note: ilmn=(il,im,in) and ilmn<=jlmn
211 : : !! klm_diag(lmn2_size)=1 il==jl and im==jm, 0 otherwise.
212 : : !!
213 : : !! SOURCE
214 : :
215 : 12 : subroutine make_indklmn(lcutdens,lmn_size,lmn2_size,indlmn,indklmn,klm_diag)
216 : :
217 : : !Arguments ------------------------------------
218 : : !scalars
219 : : integer,intent(in) :: lmn2_size,lmn_size
220 : : integer,intent(in) :: lcutdens
221 : : !scalars
222 : : integer,intent(in) :: indlmn(6,lmn_size)
223 : : integer,intent(out) :: indklmn(8,lmn2_size)
224 : : integer,intent(out) :: klm_diag(lmn2_size)
225 : :
226 : : !Local variables ------------------------------
227 : : !scalars
228 : : integer :: i0lm,i0ln,il,ilm,ilmn,iln
229 : : integer :: j0lm,j0lmn,j0ln,jl,jlm,jlmn,jln,klmn
230 : :
231 : : !************************************************************************
232 : :
233 [ + + ]: 584 : klm_diag=0
234 : :
235 [ + + ]: 98 : do jlmn=1,lmn_size
236 : 86 : jl= indlmn(1,jlmn); jlm=indlmn(4,jlmn); jln=indlmn(5,jlmn)
237 : 86 : j0lmn=jlmn*(jlmn-1)/2
238 : 86 : j0lm =jlm *(jlm -1)/2
239 : 86 : j0ln =jln *(jln -1)/2
240 [ + + ]: 670 : do ilmn=1,jlmn
241 : 572 : il=indlmn(1,ilmn); ilm=indlmn(4,ilmn); iln=indlmn(5,ilmn)
242 : 572 : klmn=j0lmn+ilmn
243 [ + + ]: 572 : if (ilm<=jlm) then
244 : 483 : indklmn(1,klmn)=j0lm+ilm !klm
245 : : else
246 : 89 : i0lm=ilm*(ilm-1)/2
247 : 89 : indklmn(1,klmn)=i0lm+jlm
248 : : end if
249 [ + + ]: 572 : if (iln<=jln) then
250 : 452 : indklmn(2,klmn)=j0ln+iln !kln
251 : : else
252 : 120 : i0ln=iln*(iln-1)/2
253 : 120 : indklmn(2,klmn)=i0ln+jln
254 : : end if
255 : : !MG This is not safe, what happens if lcutdens < |il-jl|?
256 : 572 : indklmn(3,klmn)=MIN(ABS(il-jl),lcutdens) ! abs(li-lj) NB this is a l-value, not an index >=1.
257 : 572 : indklmn(4,klmn)=MIN(il+jl,lcutdens) ! abs(li+lj) NB this is a l-value, not an index >=1.
258 : :
259 : 572 : indklmn(5,klmn)=ilm ! ilm
260 : 572 : indklmn(6,klmn)=jlm ! jlm
261 : :
262 : 572 : indklmn(7,klmn)=ilmn ! ilmn
263 : 572 : indklmn(8,klmn)=jlmn ! jlmn
264 : :
265 : :
266 [ + + ]: 658 : if (ilm==jlm) klm_diag(klmn)=1
267 : : end do
268 : : end do
269 : :
270 : 12 : end subroutine make_indklmn
271 : : !!***
272 : :
273 : : !----------------------------------------------------------------------
274 : :
275 : : !!****f* m_paw_lmn/make_kln2ln
276 : : !! NAME
277 : : !! make_kln2ln
278 : : !!
279 : : !! FUNCTION
280 : : !! Performs the setup of the kln2ln table for PAW calculations.
281 : : !!
282 : : !! INPUTS
283 : : !! lmn_size=Number of (l,m,n) elements for the PAW basis set.
284 : : !! lmn2_size=Number of elements in the symmetric basis set: lmn2_size=lmn_size*(lmn_size+1)/2
285 : : !! ln2_size=Number of symmetric (l,n) channels i.e. ln_size*(ln_size+1)/2
286 : : !! indlmn(6,lmn_size)=Array giving l,m,n,lm,ln,spin for i=lmn.
287 : : !! indklmn(8,lmn2_size)=Array giving klm, kln, abs(il-jl), (il+jl), ilm and jlm, ilmn and jlmn
288 : : !! for each klmn=(ilmn,jlmn). Note: ilmn=(il,im,in) and ilmn<=jlmn
289 : : !!
290 : : !! OUTPUT
291 : : !! kln2ln(6,ln2_size)=Table giving il, jl ,in, jn, iln, jln for each kln=(iln,jln)
292 : : !! where iln=(il,in) and iln<=jln. NB: kln2ln is an application and not a bijection
293 : : !!
294 : : !! SOURCE
295 : :
296 : 6 : subroutine make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
297 : :
298 : : !Arguments ------------------------------------
299 : : !scalars
300 : : integer,intent(in) :: lmn2_size,lmn_size,ln2_size
301 : : !arrays
302 : : integer,intent(in) :: indlmn(6,lmn_size)
303 : : integer,intent(in) :: indklmn(8,lmn2_size)
304 : : integer,intent(out) :: kln2ln(6,ln2_size)
305 : :
306 : : !Local variables ------------------------------
307 : : !scalars
308 : : integer :: il,in,ilmn,iln
309 : : integer :: jl,jn,jlmn,jln
310 : : integer :: klmn,kln,kln_old
311 : :
312 : : !************************************************************************
313 : :
314 [ + + + + ]: 580 : kln2ln = 0
315 : :
316 : 6 : kln_old = -1
317 [ + + ]: 492 : do klmn=1,lmn2_size
318 : 486 : kln = indklmn(2,klmn)
319 [ + + ]: 492 : if (kln /= kln_old) then
320 : 254 : kln_old = kln
321 : :
322 : 254 : call klmn2ijlmn(klmn,lmn_size,ilmn,jlmn)
323 : :
324 : 254 : il= indlmn(1,ilmn); iln=indlmn(5,ilmn)
325 : 254 : jl= indlmn(1,jlmn); jln=indlmn(5,jlmn)
326 : : !$in = indklmn(9 ,klmn) !=indlmn(3,ilmn)
327 : : !$jn = indklmn(10,klmn) !=indlmn(3,jlmn)
328 : :
329 : 254 : in = indlmn(3,ilmn)
330 : 254 : jn = indlmn(3,jlmn)
331 : :
332 : : ! Shift il to have the index instead of the value of l.
333 : 254 : kln2ln(1,kln) = il +1
334 : 254 : kln2ln(2,kln) = jl +1
335 : 254 : kln2ln(3,kln) = in
336 : 254 : kln2ln(4,kln) = jn
337 : 254 : kln2ln(5,kln) = iln
338 : 254 : kln2ln(6,kln) = jln
339 : : end if
340 : : end do !klmn
341 : :
342 : 6 : end subroutine make_kln2ln
343 : : !!***
344 : :
345 : : !----------------------------------------------------------------------
346 : :
347 : : !!****f* m_paw_lmn/make_klm2lm
348 : : !! NAME
349 : : !! make_klm2lm
350 : : !!
351 : : !! FUNCTION
352 : : !! Performs the setup of the klm2lm table for PAW calculations.
353 : : !!
354 : : !! INPUTS
355 : : !! lmn_size=Number of (l,m,n) elements for the PAW basis set.
356 : : !! lmn2_size=Number of elements in the symmetric basis set: lmn2_size=lmn_size*(lmn_size+1)/2
357 : : !! lm2_size)=Number of (l.m) elements in the symmetric basis set.
358 : : !! indlmn(6,lmn_size)=Array giving l,m,n,lm,ln,spin for i=lmn.
359 : : !! indklmn(8,lmn2_size)=Array giving klm, kln, abs(il-jl), (il+jl), ilm and jlm
360 : : !! for each klmn=(ilmn,jlmn). Note: ilmn=(il,im,in) and ilmn<=jlmn
361 : : !!
362 : : !! OUTPUT
363 : : !! klm2lm(6,lm2_size)=Table giving il, jl ,im, jm, ilm, jlm for each klm=(ilm,jlm)
364 : : !! where ilm=(il,im) and ilm<=jlm. NB: klm2lm is an application and not a bijection.
365 : : !!
366 : : !! NOTES
367 : : !! klm2lm can be calculated easily if we assume that all (l,m) channels
368 : : !! are ordered by increasing l and m. This is the standard convention
369 : : !! used in most of the PAW datasets. This routines, howevever, works
370 : : !! works also in the unlikely case in with (l,m) are not ordered.
371 : : !!
372 : : !! SOURCE
373 : :
374 : 3 : subroutine make_klm2lm(lmn_size,lmn2_size,lm2_size,indlmn,indklmn,klm2lm)
375 : :
376 : : !Arguments ------------------------------------
377 : : !scalars
378 : : integer,intent(in) :: lmn2_size,lmn_size,lm2_size
379 : : !arrays
380 : : integer,intent(in) :: indlmn(6,lmn_size)
381 : : integer,intent(in) :: indklmn(8,lmn2_size)
382 : : integer,intent(out) :: klm2lm(6,lm2_size)
383 : :
384 : : !Local variables ------------------------------
385 : : !scalars
386 : : integer :: il,ilm,ilmn,jl,jlm,jlmn,im,jm
387 : : integer :: klmn,klm,klm_old !,iklm
388 : :
389 : : !************************************************************************
390 : :
391 [ + + + + ]: 458 : klm2lm = 0
392 : 3 : klm_old = -1
393 [ + + ]: 246 : do klmn=1,lmn2_size
394 : 243 : klm = indklmn(1,klmn)
395 [ + + ]: 246 : if (klm /= klm_old) then
396 : 209 : klm_old = klm
397 : :
398 : 209 : ilm = indklmn(5,klmn)
399 : 209 : jlm = indklmn(6,klmn)
400 : :
401 : 209 : call klmn2ijlmn(klmn,lmn_size,ilmn,jlmn)
402 : :
403 : 209 : il=indlmn(1,ilmn); im=indlmn(2,ilmn)
404 : 209 : jl=indlmn(1,jlmn); jm=indlmn(2,jlmn)
405 : :
406 : : !shift to have the index instead of l or m.
407 : 209 : klm2lm(1,klm) = il +1 ! il
408 : 209 : klm2lm(2,klm) = jl +1 ! jl
409 : 209 : klm2lm(3,klm) = im + il +1 ! im
410 : 209 : klm2lm(4,klm) = jm + jl +1 ! jm
411 : 209 : klm2lm(5,klm) = ilm ! ilm
412 : 209 : klm2lm(6,klm) = jlm ! jlm
413 : : end if
414 : : end do !klmn
415 : :
416 [ + + + + : 458 : if (ANY(klm2lm==0)) then
+ - - + ]
417 : 0 : LIBPAW_BUG("check klm2lm")
418 : : end if
419 : :
420 : : !DEBUG
421 : : #if 0
422 : : write(std_out,*)"Debugging make_klm2lm:"
423 : : do iklm=1,lm2_size
424 : : il = klm2lm(1,iklm) -1 ! li
425 : : jl = klm2lm(2,iklm) -1 ! lj
426 : : im = klm2lm(3,iklm) -il -1 ! mi
427 : : jm = klm2lm(4,iklm) -jl -1 ! mj
428 : : ilm = klm2lm(5,iklm) ! ilm
429 : : jlm = klm2lm(6,iklm) ! jlm
430 : : write(std_out,'(i3,2(a,2i3,a))')"iklm ",iklm," l m (",il,im,")"," l m (",jl,jm,")"
431 : : end do
432 : : #endif
433 : :
434 : 3 : end subroutine make_klm2lm
435 : : !!***
436 : :
437 : : !----------------------------------------------------------------------
438 : :
439 : : !!****f* m_paw_lmn/klmn2ijlmn
440 : : !! NAME
441 : : !! klmn2ijlmn
442 : : !!
443 : : !! FUNCTION
444 : : !! Find ilmn and jlmn from klmn and lmn_size.
445 : : !!
446 : : !! INPUTS
447 : : !! lmn_size=Number of (l,m,n) elements for the PAW basis set.
448 : : !! klmn=The index corresponding to (ilmn,jlmn) in packed form.
449 : : !!
450 : : !! OUTPUT
451 : : !! jlmn, ilmn=The two symmetrix indices corresponding to klmn. NB: jlmn >= ilmn
452 : : !!
453 : : !! SOURCE
454 : :
455 : 463 : subroutine klmn2ijlmn(klmn,lmn_size,ilmn,jlmn)
456 : :
457 : : !Arguments ------------------------------------
458 : : !scalars
459 : : integer,intent(in) :: klmn,lmn_size
460 : : integer,intent(out) :: ilmn,jlmn
461 : :
462 : : !Local variables ------------------------------
463 : : !scalars
464 : : integer :: ii,jj,k0
465 : : !************************************************************************
466 : :
467 : 463 : ilmn=-1; jlmn=-1
468 : :
469 [ + - ]: 4531 : do jj=1,lmn_size
470 : 4531 : k0=jj*(jj-1)/2
471 [ + + ]: 31759 : do ii=1,jj
472 [ + + ]: 31759 : if (klmn==ii+k0) then
473 : 463 : ilmn = ii
474 : 463 : jlmn = jj; RETURN
475 : : end if
476 : : end do
477 : : end do
478 : :
479 : 0 : LIBPAW_BUG("Not able to found ilmn and jlmn")
480 : :
481 : : end subroutine klmn2ijlmn
482 : : !!***
483 : :
484 : : !----------------------------------------------------------------------
485 : :
486 : : !!****f* m_paw_lmn/make_indln
487 : : !! NAME
488 : : !! make_indln
489 : : !!
490 : : !! FUNCTION
491 : : !! Performs the setup of the indln table for PAW calculations.
492 : : !! Compute the indln indexes giving ilmn=(l,m,n)
493 : : !!
494 : : !! INPUTS
495 : : !! lmn_size=Number of (l,m,n) elements for the PAW basis set.
496 : : !! ln_size=Number of (l,n) elements
497 : : !! indlmn(6,lmn_size)=Array giving l,m,n,lm,ln,spin for i=lmn.
498 : : !!
499 : : !! OUTPUT
500 : : !! indln(2,ln_size)=Array giving l and n for i=ln
501 : : !!
502 : : !! SOURCE
503 : :
504 : 3 : subroutine make_indln(lmn_size,ln_size,indlmn,indln)
505 : :
506 : : !Arguments ------------------------------------
507 : : !scalars
508 : : integer,intent(in) :: lmn_size,ln_size
509 : : !arrays
510 : : integer,intent(in) :: indlmn(6,lmn_size)
511 : : integer,intent(out) :: indln(2,ln_size)
512 : :
513 : : !Local variables ------------------------------
514 : : !scalars
515 : : integer :: ilmn,ll,nn,ll_old,nn_old,ii
516 : : !************************************************************************
517 : :
518 : 3 : ii=0; ll_old=-1; nn_old=-1
519 [ + + ]: 37 : do ilmn=1,lmn_size
520 : 34 : ll=indlmn(1,ilmn)
521 : 34 : nn=indlmn(3,ilmn)
522 [ + + ]: 37 : if (ll/=ll_old.or.nn/=nn_old) then
523 : 14 : ll_old=ll
524 : 14 : nn_old=nn
525 : 14 : ii=ii+1
526 : 14 : indln(1,ii)=ll
527 : 14 : indln(2,ii)=nn
528 : : end if
529 : : end do
530 : :
531 [ - + ]: 3 : if(ii/=ln_size) LIBPAW_ERROR("ii/=ln_size")
532 : :
533 : 3 : end subroutine make_indln
534 : : !!***
535 : :
536 : : !----------------------------------------------------------------------
537 : :
538 : : !!****f* m_paw_lmn/uppert_index
539 : : !! NAME
540 : : !! uppert_index
541 : : !!
542 : : !! FUNCTION
543 : : !! Helper function returning the sequential index of an element in the upper triangle of a matrix
544 : : !! given the row-index ii and the column-index jj. If ii>jj the index of the element a_{jj,ii} is returned.
545 : : !!
546 : : !! INPUTS
547 : : !! ii=Row index
548 : : !! jj=column index
549 : : !!
550 : : !! SOURCE
551 : :
552 : 0 : function uppert_index(ii,jj)
553 : :
554 : : !Arguments ------------------------------------
555 : : !scalars
556 : : integer,intent(in) :: ii,jj
557 : : integer :: uppert_index
558 : : !scalars
559 : :
560 : : !************************************************************************
561 : :
562 : 0 : if (jj>=jj) then
563 : 0 : uppert_index = ii + jj*(jj-1)/2
564 : : else
565 : : uppert_index = jj + ii*(ii-1)/2
566 : : end if
567 : :
568 : 0 : end function uppert_index
569 : : !!***
570 : :
571 : : END MODULE m_paw_lmn
572 : : !!***
|
