Branch data Line data Source code
1 : : !!****m* ABINIT/m_paw_io
2 : : !! NAME
3 : : !! m_paw_io
4 : : !!
5 : : !! FUNCTION
6 : : !! PAW I/O related operations
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2012-2026 ABINIT group (MT, TR)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! NOTES
15 : : !! FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
16 : : !! please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
17 : : !!
18 : : !! SOURCE
19 : :
20 : : #include "libpaw.h"
21 : :
22 : : module m_paw_io
23 : :
24 : : USE_DEFS
25 : : USE_MSG_HANDLING
26 : : USE_MEMORY_PROFILING
27 : :
28 : : implicit none
29 : :
30 : : private
31 : :
32 : : public :: pawio_print_ij
33 : : !!***
34 : :
35 : : CONTAINS !===========================================================
36 : :
37 : : !!****f* m_paw_io/pawio_print_ij
38 : : !! NAME
39 : : !! pawio_print_ij
40 : : !!
41 : : !! FUNCTION
42 : : !! Print ij_ square matrixes in a "suitable" format.
43 : : !! Data are "energy-like" in Hartree units.
44 : : !! Devoted to the printing of rhoij, Dij -like PAW matrixes.
45 : : !!
46 : : !! INPUTS
47 : : !! a_ij(cplex*adim)= input square matrix
48 : : !! asym_ij(cplex*adim)= -OPTIONAL ARGUMENT-
49 : : !! When present, A(j,i) is deduced from asym_ij
50 : : !! instead of a_ij
51 : : !! adim= dimension of array a_ij:
52 : : !! adim=ndim*(ndim+1)/2 if opt_pack= 0
53 : : !! adim=number of non-zero values of a_ij if opt_pack=+1
54 : : !! cplex=1 if a_ij is real, 2 if it is complex
55 : : !! [mode_paral]= --optional argument, default='COLL'--
56 : : !! 'COLL' if all procs are calling the routine with the same message to be written once only.
57 : : !! 'PERS' if the procs are calling the routine with different messages each to be written,
58 : : !! or if one proc is calling the routine
59 : : !! ndim= dimension of input square matrix
60 : : !! opt_l= if <0 all parts of a_ij are printed
61 : : !! if >=0 only parts of a_ij corresponding to li=lj=opt_l are printed
62 : : !! opt_l_index(ndim)= array giving l quantum number for each 1<=ilmn<=ndim
63 : : !! not used if opt_l<0
64 : : !! opt_pack= 0 if a_ij is given as A(j(j-1)/2+i), i<=j
65 : : !! +1 if a_ij is given as A(j(j-1)/2+i) and is in "packed storage"
66 : : !! (i.e. only non-zero values are stored)
67 : : !! opt_prtvol= >0 if up to 12 components of _ij matrix have to be printed
68 : : !! <0 if all components of ij_ matrix have to be printed
69 : : !! =0 if _ij matrix is not to be printed
70 : : !! (except if the warning on high values is raised - see test_value)
71 : : !! opt_sym= -OPTIONAL ARGUMENT- (default if not present: opt_sym=2)
72 : : !! Define the symmetry of a_ij matrix:
73 : : !! opt_sym=1 : A(j,i)= A(i,j)
74 : : !! opt_sym=2 : A(j,i)= Conjg[A(i,j)]
75 : : !! opt_sym=3 : A(j,i)=-A(i,j)
76 : : !! opt_sym=4 : A(j,i)=-Conjg[A(i,j)]
77 : : !! When asym_ij argument is present, A[i,j] is taken from it.
78 : : !! pack2ij(adim)= gives the (i,j) index of of packed value of rhoij
79 : : !! used only if opt_packed=+1
80 : : !! test_value= (real number) if positive, print a warning when the
81 : : !! magnitude of a_ij is greater than opt_test
82 : : !! No test when test_value<0
83 : : !! title = if given then print this string before the matrix
84 : : !! force_print = if puts to true prints in any case even if opt_prtvol is 0
85 : : !! unit=the unit number for output
86 : : !! Ha_or_eV= 1: output in hartrees, 2: output in eV
87 : : !!
88 : : !! SOURCE
89 : :
90 : 225628 : subroutine pawio_print_ij(unit,a_ij,adim,cplex,ndim,opt_l,opt_l_index, &
91 : 112814 : & opt_pack,opt_prtvol,pack2ij,test_value,Ha_or_eV, &
92 : 78407 : & mode_paral,opt_sym,asym_ij,title,force_print) ! Optional arguments
93 : :
94 : : !Arguments ---------------------------------------------
95 : : !scalars
96 : : integer,intent(in) :: adim,cplex,ndim,opt_l,opt_pack,opt_prtvol,unit,Ha_or_eV
97 : : integer,intent(in),optional :: opt_sym
98 : : character(len=*),optional,intent(in) :: mode_paral
99 : : real(dp),intent(in) :: test_value
100 : : character(len=2000),optional,intent(in) :: title
101 : : logical,optional,intent(in) :: force_print
102 : : !arrays
103 : : integer,intent(in) :: opt_l_index(ndim*min(1+opt_l,1)),pack2ij(adim*opt_pack)
104 : : real(dp),intent(in) :: a_ij(cplex*adim)
105 : : real(dp),intent(in),optional :: asym_ij(cplex*adim)
106 : :
107 : : !Local variables ---------------------------------------
108 : : ! Adjust format bellow according to maxprt
109 : : !scalars
110 : : integer,parameter :: maxprt_default=12
111 : : integer :: dplex,ilmn,ilmn1,j0lmn,jlmn,jlmn1,klmn,klmn1,klmn2,maxprt,nhigh
112 : : integer :: nmin,optsym
113 : : real(dp) :: testval
114 : : logical :: use_asym
115 : : character(len=4) :: mode_paral_
116 : : character(len=500) :: msg=''
117 : : logical :: l_force_print
118 : : !arrays
119 : : real(dp),parameter :: fact_re(4)=(/one,one,-one,-one/),fact_im(4)=(/one,-one,-one,one/)
120 : 225628 : real(dp) :: tabmax(cplex),tabmin(cplex)
121 : 112814 : real(dp),allocatable :: b_ij(:),bsym_ij(:),prtab(:,:,:),out_arr(:)
122 : :
123 : : ! *************************************************************************
124 : :
125 : : 10 format(100(1x,f9.5))
126 : : 11 format(12(1x,f9.5),a) !Change this format according to variable "maxprt"
127 : :
128 : :
129 : : !DEBUG
130 : : !write(std_out,*)' pawio_print_ij : enter '
131 : : !ENDDEBUG
132 : :
133 : : !Optional arguments
134 [ + + - + ]: 112814 : mode_paral_='COLL';if (present(mode_paral)) mode_paral_=mode_paral
135 : 112814 : use_asym=present(asym_ij)
136 [ + + ]: 112814 : if (present(opt_sym)) then
137 : 112811 : optsym=opt_sym
138 : : else
139 : : optsym=2
140 : : end if
141 [ + + ]: 112814 : l_force_print=.false.; if (present(force_print)) l_force_print=.true.
142 : :
143 : : !Define size of square matrix
144 [ + + ]: 112814 : if (opt_prtvol>=0) then
145 : 64161 : maxprt=maxprt_default
146 : : else
147 : 48653 : maxprt=ndim
148 : : end if
149 : 112814 : nmin=min(ndim,maxprt)
150 : :
151 [ + + + + : 129966 : if (opt_l>=0) nmin=count(opt_l_index(:)==opt_l)
+ + ]
152 [ - + + - : 564070 : LIBPAW_ALLOCATE(prtab,(cplex,nmin,nmin))
- + + - -
+ + - + -
+ - - + -
+ ]
153 : 112814 : dplex=cplex-1
154 : :
155 : : !Eventually unpack input matrix(es)
156 [ + - - + : 338442 : LIBPAW_ALLOCATE(b_ij,(cplex*ndim*(ndim+1)/2))
+ - + - +
- - + -
+ ]
157 [ + + ]: 112814 : if (opt_pack==0) then
158 [ + - - + : 15693883 : b_ij=a_ij
- - + + ]
159 [ + - ]: 18133 : else if (opt_pack==1) then
160 [ + + ]: 2513663 : b_ij=zero
161 [ + + ]: 1089146 : do klmn=1,adim
162 : 1071013 : klmn1=cplex*klmn-dplex
163 : 1071013 : klmn2=cplex*pack2ij(klmn)-dplex
164 [ + + ]: 2785472 : b_ij(klmn2:klmn2+dplex)=a_ij(klmn1:klmn1+dplex)
165 : : end do
166 : : end if
167 [ + + + + ]: 112814 : if (opt_prtvol<0.and.opt_l<0) then
168 [ + + ]: 48361 : if (cplex==1) then
169 [ + - - + : 3229103 : tabmax(1)=maxval(abs(b_ij))
- - + + ]
170 [ + - - + : 3229103 : tabmin(1)=minval(abs(b_ij))
- - + + ]
171 : : else
172 [ + + + + ]: 110376 : tabmax(1:2)=zero;tabmin(1:2)=1.d20
173 [ + + ]: 3060972 : do klmn=1,size(b_ij)/cplex
174 : 3042576 : klmn2=2*klmn
175 : 3042576 : tabmax(1)=max(tabmax(1),b_ij(klmn2-1))
176 : 3042576 : tabmin(1)=min(tabmin(1),b_ij(klmn2-1))
177 : 3042576 : tabmax(2)=max(tabmax(2),b_ij(klmn2 ))
178 : 3060972 : tabmin(2)=min(tabmin(2),b_ij(klmn2 ))
179 : : end do
180 : : end if
181 : : end if
182 [ + + ]: 112814 : if (use_asym) then
183 [ - + + - : 156814 : LIBPAW_ALLOCATE(bsym_ij,(cplex*ndim*(ndim+1)/2))
- + - + ]
184 [ + - ]: 78407 : if (opt_pack==0) then
185 [ + - - + : 13223208 : bsym_ij=asym_ij
- - + + ]
186 [ # # ]: 0 : else if (opt_pack==1) then
187 [ # # ]: 0 : bsym_ij=zero
188 [ # # ]: 0 : do klmn=1,adim
189 : 0 : klmn1=cplex*klmn-dplex
190 : 0 : klmn2=cplex*pack2ij(klmn)-dplex
191 [ # # ]: 0 : bsym_ij(klmn2:klmn2+dplex)=asym_ij(klmn1:klmn1+dplex)
192 : : end do
193 : : end if
194 [ + + + - ]: 78407 : if (opt_prtvol<0.and.opt_l<0) then
195 [ + + ]: 24800 : if (cplex==1) then
196 [ + - - + : 1941652 : tabmax(1)=max(tabmax(1),maxval(abs(bsym_ij)))
- - + + ]
197 [ + - - + : 1941652 : tabmin(1)=min(tabmin(1),minval(abs(bsym_ij)))
- - + + ]
198 : : else
199 [ + + ]: 244484 : do klmn=1,ndim
200 : 231568 : klmn2=2*klmn
201 : 231568 : tabmax(1)=max(tabmax(1),bsym_ij(klmn2-1))
202 : 231568 : tabmin(1)=min(tabmin(1),bsym_ij(klmn2-1))
203 : 231568 : tabmax(2)=max(tabmax(2),bsym_ij(klmn2 ))
204 : 244484 : tabmin(2)=min(tabmin(2),bsym_ij(klmn2 ))
205 : : end do
206 : : end if
207 : : end if
208 : : end if
209 : :
210 : : !Transfer triangular matrix to rectangular one
211 : 112814 : jlmn1=0
212 [ + + ]: 1597050 : do jlmn=1,ndim
213 [ + + ]: 1484236 : if (opt_l<0) then
214 [ + + ]: 1468032 : jlmn1=jlmn;if (jlmn1>nmin) cycle
215 [ + + ]: 16204 : else if (opt_l_index(jlmn)==opt_l) then
216 : 7928 : jlmn1=jlmn1+1
217 : : else
218 : : cycle
219 : : end if
220 : 1290258 : ilmn1=0;j0lmn=jlmn*(jlmn-1)/2
221 [ + + ]: 10873796 : do ilmn=1,jlmn
222 [ + + ]: 9470724 : if (opt_l<0) then
223 : : ilmn1=ilmn
224 [ + + ]: 115568 : else if (opt_l_index(ilmn)==opt_l) then
225 : 40888 : ilmn1=ilmn1+1
226 : : else
227 : : cycle
228 : : end if
229 : 9396044 : klmn=j0lmn+ilmn
230 [ + + ]: 10880280 : if (cplex==1) then
231 : 5396796 : prtab(1,ilmn1,jlmn1)=b_ij(klmn)
232 [ + + ]: 5396796 : if (use_asym) then
233 : 3709789 : prtab(1,jlmn1,ilmn1)=fact_re(optsym)*bsym_ij(klmn)
234 : : else
235 : 1687007 : prtab(1,jlmn1,ilmn1)=fact_re(optsym)*b_ij(klmn)
236 : : end if
237 : : else
238 : 3999248 : klmn=2*klmn
239 [ + + ]: 11997744 : prtab(1:2,ilmn1,jlmn1)=b_ij(klmn-1:klmn)
240 [ + + ]: 3999248 : if (use_asym) then
241 : 2981952 : prtab(1,jlmn1,ilmn1)=fact_re(optsym)*bsym_ij(klmn-1)
242 : 2981952 : prtab(2,jlmn1,ilmn1)=fact_im(optsym)*bsym_ij(klmn )
243 : : else
244 : 1017296 : prtab(1,jlmn1,ilmn1)=fact_re(optsym)*b_ij(klmn-1)
245 : 1017296 : prtab(2,jlmn1,ilmn1)=fact_im(optsym)*b_ij(klmn )
246 : : end if
247 : : end if
248 : : end do
249 : : end do
250 [ - + ]: 112814 : LIBPAW_DEALLOCATE(b_ij)
251 : :
252 [ + + ]: 112814 : if (use_asym) then
253 [ - + ]: 78407 : LIBPAW_DEALLOCATE(bsym_ij)
254 : : end if
255 : :
256 : : ! Test if the matrix contains high values if required
257 : 112814 : nhigh=0
258 [ + + ]: 112814 : if (test_value>zero) then
259 : 37784 : testval=test_value
260 : : !;if (Ha_or_eV==2) testval=testval*Ha_eV
261 [ + + + + : 6675773 : nhigh=count(abs(prtab(:,:,:))>=testval)
+ + + + ]
262 : : end if
263 : :
264 [ + + + + ]: 112814 : if (opt_prtvol/=0 .or. nhigh>0 .or. l_force_print) then
265 : :
266 [ + + ]: 76920 : if (present(title)) call wrtout(unit,title,mode_paral_)
267 : :
268 [ + + ]: 76920 : if (Ha_or_eV==2) then
269 [ + + + + : 70056 : prtab=prtab*Ha_eV
+ + ]
270 [ - + - - ]: 464 : if (opt_prtvol<0.and.opt_l<0) then
271 [ # # ]: 0 : tabmax=tabmax*Ha_eV
272 [ # # ]: 0 : tabmin=tabmin*Ha_eV
273 : : end if
274 : : end if
275 : :
276 [ + + ]: 76920 : if (cplex==2) then
277 : 25499 : write(msg,'(3x,a)') '=== REAL PART:'
278 : 25499 : call wrtout(unit,msg,mode_paral_)
279 : : end if
280 : :
281 [ - + + - : 230760 : LIBPAW_ALLOCATE(out_arr,(nmin))
+ - + - -
+ - + ]
282 [ + + + + ]: 76920 : if (ndim<=maxprt.or.opt_l>=0) then
283 [ + + ]: 849690 : do ilmn=1,nmin
284 [ + - - + : 13881648 : out_arr = prtab(1,1:nmin,ilmn)
- - + + ]
285 : 790358 : write(msg,fmt=10) out_arr
286 : 849690 : call wrtout(unit,msg,mode_paral_)
287 : : end do
288 : : else
289 [ + + ]: 228644 : do ilmn=1,nmin
290 [ + - - + : 2954784 : out_arr = prtab(1,1:nmin,ilmn)
- - + + ]
291 : 211056 : write(msg,fmt=11) out_arr,' ...'
292 : 228644 : call wrtout(unit,msg,mode_paral_)
293 : : end do
294 : 17588 : write(msg,'(3x,a,i2,a)') '... only ',maxprt,' components have been written...'
295 : 17588 : call wrtout(unit,msg,mode_paral_)
296 : : end if
297 [ + + + + ]: 76920 : if (opt_prtvol<0.and.opt_l<0) then
298 : 48361 : write(msg,'(3x,2(a,es9.2))') 'max. value= ',tabmax(1),', min. value= ',tabmin(1)
299 : 48361 : call wrtout(unit,msg,mode_paral_)
300 : : end if
301 : :
302 [ + + ]: 76920 : if (cplex==2) then
303 : 25499 : write(msg,'(3x,a)') '=== IMAGINARY PART:'
304 : 25499 : call wrtout(unit,msg,mode_paral_)
305 [ + + + + ]: 25499 : if (ndim<=maxprt.or.opt_l>=0) then
306 [ + + ]: 356384 : do ilmn=1,nmin
307 [ + - - + : 6540304 : out_arr = prtab(2,1:nmin,ilmn)
- - + + ]
308 : 336384 : write(msg,fmt=10) out_arr
309 : 356384 : call wrtout(unit,msg,mode_paral_)
310 : : end do
311 : : else
312 [ + + ]: 71487 : do ilmn=1,nmin
313 [ + - - + : 923832 : out_arr = prtab(2,1:nmin,ilmn)
- - + + ]
314 : 65988 : write(msg,fmt=11) out_arr,' ...'
315 : 71487 : call wrtout(unit,msg,mode_paral_)
316 : : end do
317 : 5499 : write(msg,'(3x,a,i2,a)') '... only ',maxprt,' components have been written...'
318 : 5499 : call wrtout(unit,msg,mode_paral_)
319 : : end if
320 [ + + + + ]: 25499 : if (opt_prtvol<0.and.opt_l<0) then
321 : 18396 : write(msg,'(3x,2(a,es9.2))') 'max. value= ',tabmax(2),', min. value= ',tabmin(2)
322 : 18396 : call wrtout(unit,msg,mode_paral_)
323 : : end if
324 : : end if
325 [ - + ]: 76920 : LIBPAW_DEALLOCATE(out_arr)
326 : :
327 : : ! Print warning if the matrix has high values
328 [ + + ]: 76920 : if (test_value>zero) then
329 [ + + ]: 11687 : if (nhigh>0) then
330 [ + + ]: 1306 : if (Ha_or_eV==2) testval=testval*Ha_eV
331 : : write(msg,'(5a,i3,a,f6.1,7a)')&
332 : 1306 : & ' pawio_print_ij: WARNING -',ch10,&
333 : 1306 : & ' The matrix seems to have high value(s) !',ch10,&
334 : 1306 : & ' (',nhigh,' components have a value greater than ',testval,').',ch10,&
335 : 1306 : & ' It can cause instabilities during SCF convergence.',ch10,&
336 : 1306 : & ' Action: you should check your atomic dataset (psp file)',ch10,&
337 : 2612 : & ' and look for "high" projector functions...'
338 : 1306 : call wrtout(unit,msg,mode_paral_)
339 : : end if
340 : : end if
341 : :
342 : : end if
343 : :
344 [ - + ]: 112814 : LIBPAW_DEALLOCATE(prtab)
345 : :
346 : : !DEBUG
347 : : !write(std_out,*)' pawio_print_ij : exit '
348 : : !ENDDEBUG
349 : :
350 [ + + - + : 225628 : end subroutine pawio_print_ij
- + ]
351 : : !!***
352 : :
353 : : end module m_paw_io
354 : : !!***
|
