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 222256 : subroutine pawio_print_ij(unit,a_ij,adim,cplex,ndim,opt_l,opt_l_index, &
91 111128 : & opt_pack,opt_prtvol,pack2ij,test_value,Ha_or_eV, &
92 77235 : & 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 222256 : real(dp) :: tabmax(cplex),tabmin(cplex)
121 111128 : 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 111128 : mode_paral_='COLL';if (present(mode_paral)) mode_paral_=mode_paral
135 111128 : use_asym=present(asym_ij)
136 111128 : if (present(opt_sym)) then
137 111127 : optsym=opt_sym
138 : else
139 : optsym=2
140 : end if
141 111128 : l_force_print=.false.; if (present(force_print)) l_force_print=.true.
142 :
143 : !Define size of square matrix
144 111128 : if (opt_prtvol>=0) then
145 62861 : maxprt=maxprt_default
146 : else
147 48267 : maxprt=ndim
148 : end if
149 111128 : nmin=min(ndim,maxprt)
150 :
151 128168 : if (opt_l>=0) nmin=count(opt_l_index(:)==opt_l)
152 555640 : LIBPAW_ALLOCATE(prtab,(cplex,nmin,nmin))
153 111128 : dplex=cplex-1
154 :
155 : !Eventually unpack input matrix(es)
156 333384 : LIBPAW_ALLOCATE(b_ij,(cplex*ndim*(ndim+1)/2))
157 111128 : if (opt_pack==0) then
158 15614722 : b_ij=a_ij
159 17813 : else if (opt_pack==1) then
160 2493939 : b_ij=zero
161 1074814 : do klmn=1,adim
162 1057001 : klmn1=cplex*klmn-dplex
163 1057001 : klmn2=cplex*pack2ij(klmn)-dplex
164 2756424 : b_ij(klmn2:klmn2+dplex)=a_ij(klmn1:klmn1+dplex)
165 : end do
166 : end if
167 111128 : if (opt_prtvol<0.and.opt_l<0) then
168 47975 : if (cplex==1) then
169 3203740 : tabmax(1)=maxval(abs(b_ij))
170 3203740 : 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 111128 : if (use_asym) then
183 154470 : LIBPAW_ALLOCATE(bsym_ij,(cplex*ndim*(ndim+1)/2))
184 77235 : if (opt_pack==0) then
185 13156834 : 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 77235 : 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 111128 : jlmn1=0
212 1579656 : do jlmn=1,ndim
213 1468528 : if (opt_l<0) then
214 1452428 : jlmn1=jlmn;if (jlmn1>nmin) cycle
215 16100 : else if (opt_l_index(jlmn)==opt_l) then
216 7888 : jlmn1=jlmn1+1
217 : else
218 : cycle
219 : end if
220 1275000 : ilmn1=0;j0lmn=jlmn*(jlmn-1)/2
221 10776657 : do ilmn=1,jlmn
222 9390529 : if (opt_l<0) then
223 : ilmn1=ilmn
224 115128 : else if (opt_l_index(ilmn)==opt_l) then
225 40768 : ilmn1=ilmn1+1
226 : else
227 : cycle
228 : end if
229 9316169 : klmn=j0lmn+ilmn
230 10784697 : if (cplex==1) then
231 5326713 : prtab(1,ilmn1,jlmn1)=b_ij(klmn)
232 5326713 : if (use_asym) then
233 3668377 : prtab(1,jlmn1,ilmn1)=fact_re(optsym)*bsym_ij(klmn)
234 : else
235 1658336 : prtab(1,jlmn1,ilmn1)=fact_re(optsym)*b_ij(klmn)
236 : end if
237 : else
238 3989456 : klmn=2*klmn
239 11968368 : prtab(1:2,ilmn1,jlmn1)=b_ij(klmn-1:klmn)
240 3989456 : if (use_asym) then
241 2973024 : prtab(1,jlmn1,ilmn1)=fact_re(optsym)*bsym_ij(klmn-1)
242 2973024 : prtab(2,jlmn1,ilmn1)=fact_im(optsym)*bsym_ij(klmn )
243 : else
244 1016432 : prtab(1,jlmn1,ilmn1)=fact_re(optsym)*b_ij(klmn-1)
245 1016432 : prtab(2,jlmn1,ilmn1)=fact_im(optsym)*b_ij(klmn )
246 : end if
247 : end if
248 : end do
249 : end do
250 111128 : LIBPAW_DEALLOCATE(b_ij)
251 :
252 111128 : if (use_asym) then
253 77235 : LIBPAW_DEALLOCATE(bsym_ij)
254 : end if
255 :
256 : ! Test if the matrix contains high values if required
257 111128 : nhigh=0
258 111128 : if (test_value>zero) then
259 36776 : testval=test_value
260 : !;if (Ha_or_eV==2) testval=testval*Ha_eV
261 6507721 : nhigh=count(abs(prtab(:,:,:))>=testval)
262 : end if
263 :
264 111128 : if (opt_prtvol/=0 .or. nhigh>0 .or. l_force_print) then
265 :
266 76286 : if (present(title)) call wrtout(unit,title,mode_paral_)
267 :
268 76286 : 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 76286 : if (cplex==2) then
277 25451 : write(msg,'(3x,a)') '=== REAL PART:'
278 25451 : call wrtout(unit,msg,mode_paral_)
279 : end if
280 :
281 228858 : LIBPAW_ALLOCATE(out_arr,(nmin))
282 76286 : if (ndim<=maxprt.or.opt_l>=0) then
283 843254 : do ilmn=1,nmin
284 13810848 : out_arr = prtab(1,1:nmin,ilmn)
285 784532 : write(msg,fmt=10) out_arr
286 843254 : call wrtout(unit,msg,mode_paral_)
287 : end do
288 : else
289 228332 : do ilmn=1,nmin
290 2950752 : out_arr = prtab(1,1:nmin,ilmn)
291 210768 : write(msg,fmt=11) out_arr,' ...'
292 228332 : call wrtout(unit,msg,mode_paral_)
293 : end do
294 17564 : write(msg,'(3x,a,i2,a)') '... only ',maxprt,' components have been written...'
295 17564 : call wrtout(unit,msg,mode_paral_)
296 : end if
297 76286 : if (opt_prtvol<0.and.opt_l<0) then
298 47975 : write(msg,'(3x,2(a,es9.2))') 'max. value= ',tabmax(1),', min. value= ',tabmin(1)
299 47975 : call wrtout(unit,msg,mode_paral_)
300 : end if
301 :
302 76286 : if (cplex==2) then
303 25451 : write(msg,'(3x,a)') '=== IMAGINARY PART:'
304 25451 : call wrtout(unit,msg,mode_paral_)
305 25451 : if (ndim<=maxprt.or.opt_l>=0) then
306 355952 : do ilmn=1,nmin
307 6536464 : out_arr = prtab(2,1:nmin,ilmn)
308 336000 : write(msg,fmt=10) out_arr
309 355952 : 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 25451 : 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 76286 : LIBPAW_DEALLOCATE(out_arr)
326 :
327 : ! Print warning if the matrix has high values
328 76286 : if (test_value>zero) then
329 11479 : 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 111128 : LIBPAW_DEALLOCATE(prtab)
345 :
346 : !DEBUG
347 : !write(std_out,*)' pawio_print_ij : exit '
348 : !ENDDEBUG
349 :
350 222256 : end subroutine pawio_print_ij
351 : !!***
352 :
353 : end module m_paw_io
354 : !!***
|
