TABLE OF CONTENTS


ABINIT/m_paw_io [ Modules ]

[ Top ] [ 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