TABLE OF CONTENTS
ABINIT/m_paw_io [ Modules ]
NAME
m_paw_io
FUNCTION
PAW I/O related operations
COPYRIGHT
Copyright (C) 2012-2022 ABINIT group (MT, TR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
FOR DEVELOPPERS: in order to preserve the portability of libPAW library, please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
SOURCE
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
m_paw_io/pawio_print_ij [ Functions ]
[ Top ] [ m_paw_io ] [ Functions ]
NAME
pawio_print_ij
FUNCTION
Print ij_ square matrixes in a "suitable" format. Data are "energy-like" in Hartree units. Devoted to the printing of rhoij, Dij -like PAW matrixes.
INPUTS
a_ij(cplex*adim)= input square matrix asym_ij(cplex*adim)= -OPTIONAL ARGUMENT- When present, A(j,i) is deduced from asym_ij instead of a_ij adim= dimension of array a_ij: adim=ndim*(ndim+1)/2 if opt_pack= 0 adim=number of non-zero values of a_ij if opt_pack=+1 cplex=1 if a_ij is real, 2 if it is complex [mode_paral]= --optional argument, default='COLL'-- 'COLL' if all procs are calling the routine with the same message to be written once only. 'PERS' if the procs are calling the routine with different messages each to be written, or if one proc is calling the routine ndim= dimension of input square matrix opt_l= if <0 all parts of a_ij are printed if >=0 only parts of a_ij corresponding to li=lj=opt_l are printed opt_l_index(ndim)= array giving l quantum number for each 1<=ilmn<=ndim not used if opt_l<0 opt_pack= 0 if a_ij is given as A(j(j-1)/2+i), i<=j +1 if a_ij is given as A(j(j-1)/2+i) and is in "packed storage" (i.e. only non-zero values are stored) opt_prtvol= >=0 if up to 12 components of _ij matrix have to be printed <0 if all components of ij_ matrix have to be printed opt_sym= -OPTIONAL ARGUMENT- (default if not present: opt_sym=2) Define the symmetry of a_ij matrix: opt_sym=1 : A(j,i)= A(i,j) opt_sym=2 : A(j,i)= Conjg[A(i,j)] opt_sym=3 : A(j,i)=-A(i,j) opt_sym=4 : A(j,i)=-Conjg[A(i,j)] When asym_ij argument is present, A[i,j] is taken from it. pack2ij(adim)= gives the (i,j) index of of packed value of rhoij used only if opt_packed=+1 test_value= (real number) if positive, print a warning when the magnitude of a_ij is greater than opt_test No test when test_value<0 unit=the unit number for output Ha_or_eV= 1: output in hartrees, 2: output in eV
SOURCE
86 subroutine pawio_print_ij(unit,a_ij,adim,cplex,ndim,opt_l,opt_l_index, & 87 & opt_pack,opt_prtvol,pack2ij,test_value,Ha_or_eV, & 88 & mode_paral,opt_sym,asym_ij) ! Optional arguments 89 90 !Arguments --------------------------------------------- 91 !scalars 92 integer,intent(in) :: adim,cplex,ndim,opt_l,opt_pack,opt_prtvol,unit,Ha_or_eV 93 integer,intent(in),optional :: opt_sym 94 character(len=*),optional,intent(in) :: mode_paral 95 real(dp),intent(in) :: test_value 96 !arrays 97 integer,intent(in) :: opt_l_index(ndim*min(1+opt_l,1)),pack2ij(adim*opt_pack) 98 real(dp),intent(in) :: a_ij(cplex*adim) 99 real(dp),intent(in),optional :: asym_ij(cplex*adim) 100 101 !Local variables --------------------------------------- 102 ! Adjust format bellow according to maxprt 103 !scalars 104 integer,parameter :: maxprt_default=12 105 integer :: dplex,ilmn,ilmn1,j0lmn,jlmn,jlmn1,klmn,klmn1,klmn2,maxprt,nhigh 106 integer :: nmin,optsym 107 real(dp) :: testval 108 logical :: use_asym 109 character(len=4) :: mode_paral_ 110 character(len=500) :: msg='' 111 !arrays 112 real(dp),parameter :: fact_re(4)=(/one,one,-one,-one/),fact_im(4)=(/one,-one,-one,one/) 113 real(dp) :: tabmax(cplex),tabmin(cplex) 114 real(dp),allocatable :: b_ij(:),bsym_ij(:),prtab(:,:,:) 115 116 ! ************************************************************************* 117 118 10 format(100(1x,f9.5)) 119 11 format(12(1x,f9.5),a) !Change this format according to variable "maxprt" 120 121 122 !DEBUG 123 !write(std_out,*)' pawio_print_ij : enter ' 124 !ENDDEBUG 125 126 !Optional arguments 127 mode_paral_='COLL';if (present(mode_paral)) mode_paral_=mode_paral 128 use_asym=present(asym_ij) 129 if (present(opt_sym)) then 130 optsym=opt_sym 131 else 132 optsym=2 133 end if 134 135 !Define size of square matrix 136 if (opt_prtvol>=0) then 137 maxprt=maxprt_default 138 else 139 maxprt=ndim 140 end if 141 nmin=min(ndim,maxprt) 142 143 if (opt_l>=0) nmin=count(opt_l_index(:)==opt_l) 144 LIBPAW_ALLOCATE(prtab,(cplex,nmin,nmin)) 145 dplex=cplex-1 146 147 !Eventually unpack input matrix(es) 148 LIBPAW_ALLOCATE(b_ij,(cplex*ndim*(ndim+1)/2)) 149 if (opt_pack==0) then 150 b_ij=a_ij 151 else if (opt_pack==1) then 152 b_ij=zero 153 do klmn=1,adim 154 klmn1=cplex*klmn-dplex 155 klmn2=cplex*pack2ij(klmn)-dplex 156 b_ij(klmn2:klmn2+dplex)=a_ij(klmn1:klmn1+dplex) 157 end do 158 end if 159 if (opt_prtvol<0.and.opt_l<0) then 160 if (cplex==1) then 161 tabmax(1)=maxval(abs(b_ij)) 162 tabmin(1)=minval(abs(b_ij)) 163 else 164 tabmax(1:2)=zero;tabmin(1:2)=1.d20 165 do klmn=1,size(b_ij)/cplex 166 klmn2=2*klmn 167 tabmax(1)=max(tabmax(1),b_ij(klmn2-1)) 168 tabmin(1)=min(tabmin(1),b_ij(klmn2-1)) 169 tabmax(2)=max(tabmax(2),b_ij(klmn2 )) 170 tabmin(2)=min(tabmin(2),b_ij(klmn2 )) 171 end do 172 end if 173 end if 174 if (use_asym) then 175 LIBPAW_ALLOCATE(bsym_ij,(cplex*ndim*(ndim+1)/2)) 176 if (opt_pack==0) then 177 bsym_ij=asym_ij 178 else if (opt_pack==1) then 179 bsym_ij=zero 180 do klmn=1,adim 181 klmn1=cplex*klmn-dplex 182 klmn2=cplex*pack2ij(klmn)-dplex 183 bsym_ij(klmn2:klmn2+dplex)=asym_ij(klmn1:klmn1+dplex) 184 end do 185 end if 186 if (opt_prtvol<0.and.opt_l<0) then 187 if (cplex==1) then 188 tabmax(1)=max(tabmax(1),maxval(abs(bsym_ij))) 189 tabmin(1)=min(tabmin(1),minval(abs(bsym_ij))) 190 else 191 do klmn=1,ndim 192 klmn2=2*klmn 193 tabmax(1)=max(tabmax(1),bsym_ij(klmn2-1)) 194 tabmin(1)=min(tabmin(1),bsym_ij(klmn2-1)) 195 tabmax(2)=max(tabmax(2),bsym_ij(klmn2 )) 196 tabmin(2)=min(tabmin(2),bsym_ij(klmn2 )) 197 end do 198 end if 199 end if 200 end if 201 202 !Transfer triangular matrix to rectangular one 203 jlmn1=0 204 do jlmn=1,ndim 205 if (opt_l<0) then 206 jlmn1=jlmn;if (jlmn1>nmin) cycle 207 else if (opt_l_index(jlmn)==opt_l) then 208 jlmn1=jlmn1+1 209 else 210 cycle 211 end if 212 ilmn1=0;j0lmn=jlmn*(jlmn-1)/2 213 do ilmn=1,jlmn 214 if (opt_l<0) then 215 ilmn1=ilmn 216 else if (opt_l_index(ilmn)==opt_l) then 217 ilmn1=ilmn1+1 218 else 219 cycle 220 end if 221 klmn=j0lmn+ilmn 222 if (cplex==1) then 223 prtab(1,ilmn1,jlmn1)=b_ij(klmn) 224 if (use_asym) then 225 prtab(1,jlmn1,ilmn1)=fact_re(optsym)*bsym_ij(klmn) 226 else 227 prtab(1,jlmn1,ilmn1)=fact_re(optsym)*b_ij(klmn) 228 end if 229 else 230 klmn=2*klmn 231 prtab(1:2,ilmn1,jlmn1)=b_ij(klmn-1:klmn) 232 if (use_asym) then 233 prtab(1,jlmn1,ilmn1)=fact_re(optsym)*bsym_ij(klmn-1) 234 prtab(2,jlmn1,ilmn1)=fact_im(optsym)*bsym_ij(klmn ) 235 else 236 prtab(1,jlmn1,ilmn1)=fact_re(optsym)*b_ij(klmn-1) 237 prtab(2,jlmn1,ilmn1)=fact_im(optsym)*b_ij(klmn ) 238 end if 239 end if 240 end do 241 end do 242 LIBPAW_DEALLOCATE(b_ij) 243 244 if (use_asym) then 245 LIBPAW_DEALLOCATE(bsym_ij) 246 end if 247 248 if (Ha_or_eV==2) then 249 prtab=prtab*Ha_eV 250 if (opt_prtvol<0.and.opt_l<0) then 251 tabmax=tabmax*Ha_eV 252 tabmin=tabmin*Ha_eV 253 end if 254 end if 255 256 if (cplex==2) then 257 write(msg,'(3x,a)') '=== REAL PART:' 258 call wrtout(unit,msg,mode_paral_) 259 end if 260 261 if (ndim<=maxprt.or.opt_l>=0) then 262 do ilmn=1,nmin 263 write(msg,fmt=10) prtab(1,1:nmin,ilmn) 264 call wrtout(unit,msg,mode_paral_) 265 end do 266 else 267 do ilmn=1,nmin 268 write(msg,fmt=11) prtab(1,1:nmin,ilmn),' ...' 269 call wrtout(unit,msg,mode_paral_) 270 end do 271 write(msg,'(3x,a,i2,a)') '... only ',maxprt,' components have been written...' 272 call wrtout(unit,msg,mode_paral_) 273 end if 274 if (opt_prtvol<0.and.opt_l<0) then 275 write(msg,'(3x,2(a,es9.2))') 'max. value= ',tabmax(1),', min. value= ',tabmin(1) 276 call wrtout(unit,msg,mode_paral_) 277 end if 278 279 if (cplex==2) then 280 write(msg,'(3x,a)') '=== IMAGINARY PART:' 281 call wrtout(unit,msg,mode_paral_) 282 if (ndim<=maxprt.or.opt_l>=0) then 283 do ilmn=1,nmin 284 write(msg,fmt=10) prtab(2,1:nmin,ilmn) 285 call wrtout(unit,msg,mode_paral_) 286 end do 287 else 288 do ilmn=1,nmin 289 write(msg,fmt=11) prtab(2,1:nmin,ilmn),' ...' 290 call wrtout(unit,msg,mode_paral_) 291 end do 292 write(msg,'(3x,a,i2,a)') '... only ',maxprt,' components have been written...' 293 call wrtout(unit,msg,mode_paral_) 294 end if 295 if (opt_prtvol<0.and.opt_l<0) then 296 write(msg,'(3x,2(a,es9.2))') 'max. value= ',tabmax(2),', min. value= ',tabmin(2) 297 call wrtout(unit,msg,mode_paral_) 298 end if 299 end if 300 301 if (test_value>zero) then 302 testval=test_value;if (Ha_or_eV==2) testval=testval*Ha_eV 303 nhigh=0;nhigh=count(abs(prtab(:,:,:))>=testval) 304 if (nhigh>0) then 305 write(msg,'(5a,i3,a,f6.1,7a)')& 306 & ' pawio_print_ij: WARNING -',ch10,& 307 & ' The matrix seems to have high value(s) !',ch10,& 308 & ' (',nhigh,' components have a value greater than ',testval,').',ch10,& 309 & ' It can cause instabilities during SCF convergence.',ch10,& 310 & ' Action: you should check your atomic dataset (psp file)',ch10,& 311 & ' and look for "high" projector functions...' 312 call wrtout(unit,msg,mode_paral_) 313 end if 314 end if 315 316 LIBPAW_DEALLOCATE(prtab) 317 318 !DEBUG 319 !write(std_out,*)' pawio_print_ij : exit ' 320 !ENDDEBUG 321 322 end subroutine pawio_print_ij