TABLE OF CONTENTS


ABINIT/m_paw_tools [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_tools

FUNCTION

  This module contains miscelaneous routines used in the PAW context.

COPYRIGHT

 Copyright (C) 2018-2018 ABINIT group (FJ,MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_paw_tools
24 
25  use defs_basis
26  use defs_abitypes
27  use m_abicore
28  use m_errors
29  use m_xmpi
30 
31  use m_paral_atom,       only : get_my_atmtab, free_my_atmtab
32  use m_electronpositron, only : electronpositron_type,electronpositron_calctype,EP_POSITRON
33  use m_pawang,           only : pawang_type, mat_slm2ylm, mat_mlms2jmj
34  use m_pawtab,           only : pawtab_type
35  use m_paw_ij,           only : paw_ij_type, paw_ij_free, paw_ij_nullify, paw_ij_gather
36  use m_pawdij,           only : pawdij_print_dij
37  use m_pawrhoij,         only : pawrhoij_type, pawrhoij_free, pawrhoij_gather, pawrhoij_nullify
38  use m_paw_io,           only : pawio_print_ij
39  use m_paw_sphharm,      only : mat_mlms2jmj, mat_slm2ylm
40  use m_paw_correlations, only : setnoccmmp
41 
42  implicit none
43 
44  private
45 
46 !public procedures.
47  public :: chkpawovlp
48  public :: pawprt
49 
50 CONTAINS  !========================================================================================

m_paw_tools/chkpawovlp [ Functions ]

[ Top ] [ m_paw_tools ] [ Functions ]

NAME

 chkpawovlp

FUNCTION

 Verify that the PAW spheres are not overlapping

INPUTS

  natom=number of atoms in cell.
  ntypat=number of types of atoms in unit cell.
  pawovlp=percentage of voluminal overlap ratio allowed to continue execution
          (if negative value, execution always continues)
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data:
  rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).
  typat(natom)=type (integer) for each atom
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  (only checking)

NOTES

PARENTS

      bethe_salpeter,respfn,scfcv,screening,sigma,wfk_analyze

CHILDREN

      wrtout

SOURCE

 85 subroutine chkpawovlp(natom,ntypat,pawovlp,pawtab,rmet,typat,xred)
 86 
 87 
 88 !This section has been created automatically by the script Abilint (TD).
 89 !Do not modify the following lines by hand.
 90 #undef ABI_FUNC
 91 #define ABI_FUNC 'chkpawovlp'
 92 !End of the abilint section
 93 
 94  implicit none
 95 
 96 !Arguments ---------------------------------------------
 97 !scalars
 98  integer,intent(in) :: natom,ntypat
 99  real(dp) :: pawovlp
100 !arrays
101  integer,intent(in) :: typat(natom)
102  real(dp),intent(in) :: rmet(3,3),xred(3,natom)
103  type(pawtab_type),intent(in) :: pawtab(ntypat)
104 
105 !Local variables ---------------------------------------
106 !scalars
107  integer :: ia,ib,ii,t1,t2,t3
108  logical :: stop_on_error
109  real(dp) :: dd,dif1,dif2,dif3,ha,hb,norm2
110  real(dp) :: ratio_percent,va,vb,vv
111  character(len=750) :: message
112 !arrays
113  integer :: iamax(2),ibmax(2),iovl(2)
114  real(dp) :: norm2_min(2),r2cut(2),ratio_percent_max(2),rcuta(2),rcutb(2)
115 
116 
117 ! *************************************************************************
118 
119  DBG_ENTER("COLL")
120 
121  iamax(:)=-1;ibmax(:)=-1
122  norm2_min(:)=-1.d0;ratio_percent_max(:)=-1.d0
123  iovl(:)=0
124 
125 !Loop on "overlapping" atoms with the maximum overlap
126  do ia=1,natom
127 
128    rcuta(1)=pawtab(typat(ia))%rpaw
129    rcuta(2)=pawtab(typat(ia))%rshp
130 
131    do ib=ia,natom
132 
133      rcutb(1)=pawtab(typat(ib))%rpaw
134      rcutb(2)=pawtab(typat(ib))%rshp
135      r2cut(1)=(rcuta(1)+rcutb(1))**2
136      r2cut(2)=(rcuta(2)+rcutb(2))**2
137 
138 !    Visit the box and its first images:
139      do t3=-1,1
140        do t2=-1,1
141          do t1=-1,1
142 
143            dif1=xred(1,ia)-(xred(1,ib)+dble(t1))
144            dif2=xred(2,ia)-(xred(2,ib)+dble(t2))
145            dif3=xred(3,ia)-(xred(3,ib)+dble(t3))
146            norm2=sqnrm_pawovlp(dif1,dif2,dif3)
147 
148            do ii=1,2
149 
150              if(norm2>tol10.and.norm2<r2cut(ii)) then
151 
152                iovl(ii)=iovl(ii)+1
153 
154 !              Compute the overlap ratio:
155                dd=sqrt(norm2)
156                va=4._dp/3._dp*pi*rcuta(ii)**3
157                vb=4._dp/3._dp*pi*rcutb(ii)**3
158                ha=(rcutb(ii)**2-(dd-rcuta(ii))**2)/(two*dd)
159                hb=(rcuta(ii)**2-(dd-rcutb(ii))**2)/(two*dd)
160                vv=pi/3.d0*(ha**2*(three*rcuta(ii)-ha)+hb**2*(three*rcutb(ii)-hb))
161                ratio_percent=100._dp*min(vv/min(va,vb),one)
162                if (ratio_percent>ratio_percent_max(ii)) then
163                  ratio_percent_max(ii)=ratio_percent
164                  norm2_min(ii)=norm2
165                  iamax(ii)=ia;ibmax(ii)=ib
166                end if
167 
168              end if
169            end do
170          end do
171        end do
172      end do
173    end do
174  end do
175 
176  stop_on_error=(abs(pawovlp)<=tol6.or.(pawovlp>tol6.and.ratio_percent_max(1)>pawovlp))
177 
178 !Print adapted message with overlap value
179  if (iovl(1)+iovl(2)>0) then
180 
181    !ii=1: PAW augmentation regions overlap
182    !ii=2: compensation charges overlap
183    if (iovl(2)==0) ii=1
184    if (iovl(2)> 0) ii=2
185 
186    if (iovl(ii)>0) then
187 
188      if (ii==1) write(message,' (a)' ) 'PAW SPHERES ARE OVERLAPPING!'
189      if (ii==2) write(message, '(2a)' )'PAW COMPENSATION DENSITIES ARE OVERLAPPING !!!!'
190 
191      if (iovl(ii)==1) then
192        write(message, '(3a)' ) trim(message),ch10,&
193 &       '   There is one pair of overlapping atoms.'
194      else
195        write(message, '(3a,i5,a)' ) trim(message),ch10,&
196 &       '   There are ', iovl(1),' pairs of overlapping atoms.'
197      end if
198      write(message, '(3a,i3,a,i3,a)' ) trim(message),ch10,&
199      '   The maximum overlap percentage is obtained for the atoms ',iamax(ii),' and ',ibmax(ii),'.'
200      write(message, '(2a,2(a,i3),a,f9.5,a,2(a,i3,a,f9.5,a),a,f5.2,a)' ) trim(message),ch10,&
201 &     '    | Distance between atoms ',iamax(ii),' and ',ibmax(ii),' is  : ',sqrt(norm2_min(ii)),ch10,&
202 &     '    | PAW radius of the sphere around atom ',iamax(ii),' is: ',pawtab(typat(iamax(ii)))%rpaw,ch10,&
203 &     '    | PAW radius of the sphere around atom ',ibmax(ii),' is: ',pawtab(typat(ibmax(ii)))%rpaw,ch10,&
204 &     '    | This leads to a (voluminal) overlap ratio of ',ratio_percent_max(ii),' %'
205      if (ii==2) then
206        write(message, '(3a)' ) trim(message),ch10,&
207 &       'THIS IS DANGEROUS !, as PAW formalism assumes non-overlapping compensation densities.'
208      end if
209 
210      if (stop_on_error) then
211        MSG_ERROR_NOSTOP(message,ia) !ia is dummy
212      else
213        MSG_WARNING(message)
214      end if
215 
216    end if
217 
218 !  Print advice
219    if (stop_on_error) then
220      write(message, '(3a)' )&
221 &     '  Action: 1- decrease cutoff radius of PAW dataset',ch10,&
222 &     '    OR  2- ajust "pawovlp" input variable to allow overlap (risky)'
223      MSG_ERROR(message)
224    end if
225 
226 !  Print last message if execution continues:
227    if (pawovlp<=tol6) then
228      write(message, '(6a)' ) &
229 &     '       Results might be approximate,',ch10,&
230 &     '       and even inaccurate (if overlap is too big) !',ch10,&
231 &     '       Assume experienced user. Execution will continue.',ch10
232      call wrtout(std_out,message,'COLL')
233    else if (ratio_percent_max(1)<=pawovlp) then
234      write(message, '(8a)' ) &
235 &     '       Overlap ratio seems to be acceptable (less than value',ch10,&
236 &     '       of "pawovlp" input parameter): execution will continue.',ch10,&
237 &     '       But be aware that results might be approximate,',ch10,&
238 &     '       and even inaccurate (depending on your physical system) !',ch10
239      call wrtout(std_out,message,'COLL')
240    end if
241 
242  end if !iovl>0
243 
244  DBG_EXIT("COLL")
245 
246  contains
247 
248    function sqnrm_pawovlp(u1,u2,u3)
249 !squared norm of a vector
250 
251 !This section has been created automatically by the script Abilint (TD).
252 !Do not modify the following lines by hand.
253 #undef ABI_FUNC
254 #define ABI_FUNC 'sqnrm_pawovlp'
255 !End of the abilint section
256 
257    real(dp) :: sqnrm_pawovlp
258    real(dp),intent(in) :: u1,u2,u3
259 
260    sqnrm_pawovlp=rmet(1,1)*u1*u1+rmet(2,1)*u2*u1+rmet(3,1)*u3*u1&
261 &   +rmet(1,2)*u1*u2+rmet(2,2)*u2*u2+rmet(3,2)*u3*u2&
262 &   +rmet(1,3)*u1*u3+rmet(2,3)*u2*u3+rmet(3,3)*u3*u3
263 
264  end function sqnrm_pawovlp
265 
266 end subroutine chkpawovlp

m_paw_tools/pawprt [ Functions ]

[ Top ] [ m_paw_tools ] [ Functions ]

NAME

 pawprt

FUNCTION

 Print out data concerning PAW formalism
 (pseudopotential strength, augmentation occupancies...)
 To be called at the end of the SCF cycle

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (FJ,MT,BA)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

  dtset <type(dataset_type)>=all input variables in this dataset
   | enunit=parameter determining units of output energies
   | kptopt=option for the generation of k points
   | natom=number of atoms in cell
   | ntypat = number of atom types
   | pawprtvol= printing volume
   | pawspnorb=flag: 1 if spin-orbit coupling is activated
   | typat(natom)=type of each atom
  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data

OUTPUT

  (only printing)

PARENTS

      bethe_salpeter,outscfcv,screening,sigma

CHILDREN

      free_my_atmtab,get_my_atmtab,mat_mlms2jmj,mat_slm2ylm,paw_ij_free
      paw_ij_gather,paw_ij_nullify,pawio_print_ij,pawrhoij_free
      pawrhoij_gather,pawrhoij_nullify,setnoccmmp,wrtout,xmpi_comm_group
      xmpi_group_free,xmpi_group_translate_ranks

SOURCE

318 subroutine pawprt(dtset,my_natom,paw_ij,pawrhoij,pawtab,&
319 &                 electronpositron,& ! optional argument
320 &                 mpi_atmtab,comm_atom) ! optional arguments (parallelism)
321 
322 
323 !This section has been created automatically by the script Abilint (TD).
324 !Do not modify the following lines by hand.
325 #undef ABI_FUNC
326 #define ABI_FUNC 'pawprt'
327 !End of the abilint section
328 
329  implicit none
330 
331 !Arguments ------------------------------------
332 !scalars
333  integer,intent(in) :: my_natom
334  integer,optional,intent(in) :: comm_atom
335  type(dataset_type),intent(in) :: dtset
336  type(electronpositron_type),pointer,optional :: electronpositron
337 !arrays
338  integer,optional,target,intent(in) :: mpi_atmtab(:)
339  type(paw_ij_type),target,intent(inout) :: paw_ij(my_natom)
340  type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom)
341  type(pawtab_type),target,intent(in) :: pawtab(dtset%ntypat)
342 
343 !Local variables-------------------------------
344 !scalars
345  integer,parameter :: natmax=2
346  integer :: cplex_dij,group1,group2,iat,iatom,ierr,ii,im1,im2,ipositron,irhoij,ispden
347  integer :: i_unitfi,itypat,jrhoij,ll,llp,me_atom,my_comm_atom,natprt,ndij,nspden,nsppol
348  integer :: optsym,sz1,sz2,unitfi,unt
349  real(dp) :: mnorm,mx,my,mz,ntot,valmx,localm
350  logical,parameter :: debug=.false.
351  logical :: my_atmtab_allocated,paral_atom,useexexch,usepawu
352  type(pawang_type):: pawang_dum
353  character(len=7),parameter :: dspin1(6)=(/"up     ","down   ","up-up  ","dwn-dwn","up-dwn ","dwn-up "/)
354  character(len=8),parameter :: dspin2(6)=(/"up      ","down    ","dens (n)","magn (x)","magn (y)","magn (z)"/)
355  character(len=9),parameter :: dspin3(6)=(/"up       ","down     ","up-up    ","down-down","Re[up-dn]","Im[up-dn]"/)
356  character(len=500) :: msg0,msg
357 !arrays
358  integer :: idum(1)
359  integer :: idum1(0),idum3(0,0,0)
360  integer,allocatable :: jatom(:),opt_l_index(:)
361  integer,pointer :: my_atmtab(:)
362  real(dp) :: rdum2(0,0),rdum4(0,0,0,0)
363  real(dp),allocatable :: rhoijs(:,:)
364  complex(dpc),allocatable :: noccmmp_ylm(:,:,:),noccmmp_jmj(:,:),noccmmp_slm(:,:,:)
365  type(paw_ij_type), ABI_CONTIGUOUS pointer :: paw_ij_all(:)
366  type(pawrhoij_type),ABI_CONTIGUOUS pointer :: pawrhoij_all(:)
367 
368 ! *********************************************************************
369 
370  DBG_ENTER("COLL")
371 
372 !Set up parallelism over atoms
373  paral_atom=(present(comm_atom).and.(my_natom/=dtset%natom))
374  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
375  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
376  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,dtset%natom,my_natom_ref=my_natom)
377  me_atom=xmpi_comm_rank(my_comm_atom)
378 
379 !Continue only if comm_atom contains the master of the output comm
380  if (paral_atom) then
381    call xmpi_comm_group(abinit_comm_output,group1,ierr)
382    call xmpi_comm_group(my_comm_atom,group2,ierr)
383    call xmpi_group_translate_ranks(group1,1,(/0/),group2,idum,ierr)
384    call xmpi_group_free(group1)
385    call xmpi_group_free(group2)
386    if (idum(1)==xmpi_undefined) then
387      call free_my_atmtab(my_atmtab,my_atmtab_allocated)
388      return
389    end if
390  end if
391 
392 !Initializations
393  natprt=natmax;if (dtset%natom==1) natprt=1
394  if (dtset%pawprtvol<0) natprt=dtset%natom
395  ABI_ALLOCATE(jatom,(natprt))
396  if (natprt==1) then
397    jatom(1)=1
398  else if (natprt==2) then
399    jatom(1)=1;jatom(2)=dtset%natom
400  else if (natprt==dtset%natom) then
401    do iat=1,dtset%natom
402      jatom(iat)=iat
403    end do
404  else
405    MSG_BUG("invalid value of natprt!")
406  end if
407  usepawu=(count(pawtab(:)%usepawu>0)>0)
408  useexexch=(count(pawtab(:)%useexexch>0)>0)
409  ipositron=0
410  if (present(electronpositron)) then
411    if (associated(electronpositron)) ipositron=electronpositron%calctype
412  end if
413 
414 !Main title
415  write(msg, '(2a)' ) ch10,&
416 & ' ==== Results concerning PAW augmentation regions ===='
417  call wrtout(ab_out,msg,'COLL')
418  call wrtout(std_out,msg,'COLL')
419  msg=' '
420  call wrtout(ab_out,msg,'COLL')
421  call wrtout(std_out,msg,'COLL')
422 
423 !If atomic data are distributed, retrieve all Dij on master proc
424  if (paral_atom) then
425    if (me_atom==0) then
426      ABI_DATATYPE_ALLOCATE(paw_ij_all,(dtset%natom))
427      call paw_ij_nullify(paw_ij_all)
428    else
429      ABI_DATATYPE_ALLOCATE(paw_ij_all,(0))
430    end if
431    call paw_ij_gather(paw_ij,paw_ij_all,0,my_comm_atom)
432  else
433    paw_ij_all => paw_ij
434  end if
435 
436 !Print out pseudopotential strength
437 !----------------------------------
438  if (me_atom==0) then
439    do i_unitfi=1,2
440      unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
441      do unt=1,2
442        if (((unt==1).and.(dtset%enunit==0.or.dtset%enunit==2)).or.&
443 &          ((unt==2).and.(dtset%enunit==1.or.dtset%enunit==2))) then
444          if ((unt==1).and.(dtset%enunit==0.or.dtset%enunit==2)) then
445            write(msg,'(a)') ' Total pseudopotential strength Dij (hartree):'
446          else if ((unt==2).and.(dtset%enunit==1.or.dtset%enunit==2)) then
447            write(msg,'(a)') ' Total pseudopotential strength Dij (eV):'
448          end if
449          call wrtout(unitfi,msg,'COLL')
450          if (ipositron>0) then
451            if (electronpositron%has_pos_ham==0) then
452              write(msg,'(a)') ' -Note: these are the electronic Dij'
453            else
454              write(msg,'(a)') ' -Note: these are the positronic Dij'
455            end if
456            call wrtout(unitfi,msg,'COLL')
457          end if
458          valmx=100._dp;if (ipositron>0) valmx=-1._dp
459          do iat=1,natprt
460            iatom=jatom(iat)
461            call pawdij_print_dij(paw_ij_all(iatom)%dij,paw_ij_all(iatom)%cplex_dij,&
462 &                  paw_ij_all(iatom)%cplex_rf,iatom,dtset%natom,paw_ij_all(iatom)%nspden,&
463 &                  paw_ij_all(iatom)%nsppol,test_value=valmx,unit=unitfi,&
464 &                  Ha_or_eV=unt,opt_prtvol=dtset%pawprtvol)
465          end do
466        end if
467        msg=' '
468        call wrtout(unitfi,msg,'COLL')
469      end do
470    end do
471  end if
472  if (paral_atom.and.(.not.usepawu).and.(.not.useexexch)) then
473    call paw_ij_free(paw_ij_all)
474    ABI_DATATYPE_DEALLOCATE(paw_ij_all)
475  end if
476 
477 !If atomic data are distributed, retrieve all Rhoij on master proc
478  if (paral_atom) then
479    if (me_atom==0) then
480      ABI_DATATYPE_ALLOCATE(pawrhoij_all,(dtset%natom))
481    else
482      ABI_DATATYPE_ALLOCATE(pawrhoij_all,(0))
483    end if
484    call pawrhoij_nullify(pawrhoij_all)
485    call pawrhoij_gather(pawrhoij,pawrhoij_all,0,my_comm_atom,&
486 &   with_grhoij=.false.,with_lmnmix=.false.,&
487 &   with_rhoij_=.false.,with_rhoijres=.false.)
488  else
489    pawrhoij_all => pawrhoij
490  end if
491 
492 !Print out SYMMETRIZED occupancies of the partial waves
493 !------------------------------------------------------
494  if (me_atom==0) then
495    do i_unitfi=1,2
496      unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
497      write(msg,'(a)') ' Augmentation waves occupancies Rhoij:'
498      call wrtout(unitfi,msg,'COLL')
499      if (ipositron>0) then
500        if (electronpositron%particle==EP_POSITRON) then
501          write(msg,'(a)') ' -Note: these are the electronic Rhoij'
502        else
503          write(msg,'(a)') ' -Note: these are the positronic Rhoij'
504        end if
505        call wrtout(unitfi,msg,'COLL')
506      end if
507      if (dtset%pawspnorb>0.and.pawrhoij_all(1)%cplex==1.and.dtset%kptopt/=1.and.dtset%kptopt/=2) then
508        write(msg,'(6a)') ' pawprt: - WARNING:',ch10,&
509 &       '       Spin-orbit coupling is activated but only real part of Rhoij occupancies',ch10,&
510 &       '       has been computed; they could have an imaginary part (not printed here).'
511        call wrtout(unitfi,msg,'COLL')
512      end if
513      do iat=1,natprt
514        iatom=jatom(iat);nspden=pawrhoij_all(iatom)%nspden
515        optsym=2;if (pawrhoij_all(iatom)%cplex==2.and.dtset%nspinor==1) optsym=1
516        do ispden=1,nspden
517          valmx=-1._dp;if (ispden==1) valmx=25._dp
518          msg='' ; msg0=''
519          if (dtset%natom>1.or.nspden>1) write(msg0, '(a,i3)' ) ' Atom #',iatom
520          if (nspden==1) write(msg,'(a)')     trim(msg0)
521          if (nspden==2) write(msg,'(2a,i1)') trim(msg0),' - Spin component ',ispden
522          if (nspden==4) write(msg,'(3a)')    trim(msg0),' - Component ',dspin2(ispden+2*(nspden/4))
523          if (dtset%natom>1.or.nspden>1) then
524            call wrtout(unitfi,msg,'COLL')
525          end if
526          call pawio_print_ij(unitfi,pawrhoij_all(iatom)%rhoijp(:,ispden),pawrhoij_all(iatom)%nrhoijsel,&
527 &         pawrhoij_all(iatom)%cplex,pawrhoij_all(iatom)%lmn_size,-1,idum,1,dtset%pawprtvol,&
528 &         pawrhoij_all(iatom)%rhoijselect(:),valmx,1,opt_sym=optsym)
529        end do
530      end do
531      msg=' '
532      call wrtout(unitfi,msg,'COLL')
533    end do
534  end if
535 
536 !PAW+U or local exact-exchange: print out +U components of occupancies
537 !-------------------------------------------------------------------------------
538  if ((usepawu.or.useexexch).and.ipositron/=1.and.me_atom==0) then
539    do i_unitfi=1,2
540      unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
541      if(useexexch) write(msg,'(a)') &
542 &     ' "Local exact-exchange" part of augmentation waves occupancies Rhoij:'
543      if(usepawu) write(msg,'(a)') &
544 &     ' "PAW+U" part of augmentation waves occupancies Rhoij:'
545      call wrtout(unitfi,msg,'COLL')
546      valmx=-1._dp
547      do iatom=1,dtset%natom
548        nspden=pawrhoij_all(iatom)%nspden;itypat=dtset%typat(iatom)
549        cplex_dij=paw_ij_all(iatom)%cplex_dij
550        ll=-1;llp=-1
551        if (pawtab(itypat)%usepawu>0) ll=pawtab(itypat)%lpawu
552        if (pawtab(itypat)%useexexch>0) llp=pawtab(itypat)%lexexch
553        if (ll/=llp.and.ll/=-1.and.llp/=-1) then
554          MSG_BUG(" lpawu/=lexexch forbidden!")
555        end if
556        ABI_ALLOCATE(opt_l_index, (pawtab(itypat)%lmn_size))
557        ll=max(ll,llp)
558        if (ll>=0) then
559          optsym=2;if (pawrhoij_all(iatom)%cplex==2.and.dtset%nspinor==1) optsym=1
560          do ispden=1,nspden
561            msg='' ; msg0=''
562            write(msg0,'(a,i3,a,i1,a)') ' Atom #',iatom,' - L=',ll,' ONLY'
563            if (nspden==1) write(msg,'(a)')     trim(msg0)
564            if (nspden==2) write(msg,'(2a,i1)') trim(msg0),' - Spin component ',ispden
565            if (nspden==4) write(msg,'(3a)')    trim(msg0),' - Component ',dspin2(ispden+2*(nspden/4))
566            call wrtout(unitfi,msg,'COLL')
567 
568            opt_l_index = pawtab(itypat)%indlmn(1,1:pawtab(itypat)%lmn_size)
569            call pawio_print_ij(unitfi,pawrhoij_all(iatom)%rhoijp(:,ispden),pawrhoij_all(iatom)%nrhoijsel,&
570 &           pawrhoij_all(iatom)%cplex,pawrhoij_all(iatom)%lmn_size,ll,&
571 &           opt_l_index,1,dtset%pawprtvol,&
572 &           pawrhoij_all(iatom)%rhoijselect(:),valmx,1,opt_sym=optsym)
573 
574          end do
575          if (debug.and.nspden==4) then
576            sz1=paw_ij_all(iatom)%lmn2_size*cplex_dij
577            sz2=paw_ij_all(iatom)%ndij
578            ABI_ALLOCATE(rhoijs,(sz1,sz2))
579            do irhoij=1,pawrhoij_all(iatom)%nrhoijsel
580              jrhoij=cplex_dij*(irhoij-1)+1
581              rhoijs(jrhoij,1)=pawrhoij_all(iatom)%rhoijp(jrhoij,1)+pawrhoij_all(iatom)%rhoijp(jrhoij,4)
582              rhoijs(jrhoij+1,1)=pawrhoij_all(iatom)%rhoijp(jrhoij+1,1)+pawrhoij_all(iatom)%rhoijp(jrhoij+1,4)
583              rhoijs(jrhoij,2)=pawrhoij_all(iatom)%rhoijp(jrhoij,1)-pawrhoij_all(iatom)%rhoijp(jrhoij,4)
584              rhoijs(jrhoij+1,2)=pawrhoij_all(iatom)%rhoijp(jrhoij+1,1)-pawrhoij_all(iatom)%rhoijp(jrhoij+1,4)
585              rhoijs(jrhoij,3)=pawrhoij_all(iatom)%rhoijp(jrhoij,2)+pawrhoij_all(iatom)%rhoijp(jrhoij+1,3)
586              rhoijs(jrhoij+1,3)=pawrhoij_all(iatom)%rhoijp(jrhoij+1,2)-pawrhoij_all(iatom)%rhoijp(jrhoij,3)
587              rhoijs(jrhoij,4)=pawrhoij_all(iatom)%rhoijp(jrhoij ,2)-pawrhoij_all(iatom)%rhoijp(jrhoij+1,3)
588              rhoijs(jrhoij+1,4)=pawrhoij_all(iatom)%rhoijp(jrhoij+1,2)+pawrhoij_all(iatom)%rhoijp(jrhoij,3)
589            end do
590            do ispden=1,nspden
591              write(msg,'(3a)') trim(msg0),' - Component ',dspin1(ispden+2*(nspden/4))
592              opt_l_index = pawtab(itypat)%indlmn(1,1:pawtab(itypat)%lmn_size)
593 
594              call pawio_print_ij(unitfi,rhoijs(:,ispden),pawrhoij_all(iatom)%nrhoijsel,&
595 &             pawrhoij_all(iatom)%cplex,pawrhoij_all(iatom)%lmn_size,ll,&
596 &             opt_l_index,1,dtset%pawprtvol,&
597 &             pawrhoij_all(iatom)%rhoijselect(:),valmx,1,opt_sym=optsym)
598 
599              call wrtout(unitfi,"WARNING: half of the array is not correct",'COLL')
600            end do
601            ABI_DEALLOCATE(rhoijs)
602          end if
603        end if
604        ABI_DEALLOCATE(opt_l_index)
605      end do ! iatom
606      call wrtout(unitfi,' ','COLL')
607    end do
608  end if
609 
610 !PAW+U: print out occupations for correlated orbitals
611 !----------------------------------------------------
612  if (usepawu.and.ipositron/=1.and.me_atom==0) then
613    do i_unitfi=1,2
614      unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
615      write(msg,'(3a)') &
616 &     ' ---------- LDA+U DATA --------------------------------------------------- ',ch10
617      call wrtout(unitfi,msg,'COLL')
618      do iatom=1,dtset%natom
619        itypat=dtset%typat(iatom);ll=pawtab(itypat)%lpawu
620        nspden=paw_ij_all(iatom)%nspden;ndij=paw_ij_all(iatom)%ndij
621        cplex_dij=paw_ij_all(iatom)%cplex_dij
622        if ((ll>=0).and.(pawtab(itypat)%usepawu>0)) then
623          write(msg,fmt='(a,i5,a,i4,a)') " ====== For Atom ", iatom,&
624 &         ", occupations for correlated orbitals. lpawu =",ll,ch10
625          call wrtout(unitfi,msg,'COLL')
626          if(pawtab(itypat)%usepawu>=10) then
627            write(msg,fmt='(a)') "  (This is PAW atomic orbital occupations)"
628            call wrtout(unitfi,msg,'COLL')
629            write(msg,fmt='(a)') "  (For Wannier orbital occupations, refer to DFT+DMFT occupations above)"
630            call wrtout(unitfi,msg,'COLL')
631          end if
632          if(nspden==2) then
633            do ispden=1,nspden
634              write(msg,fmt='(a,i4,a,i3,a,f10.5)') " Atom", iatom,&
635 &             ". Occ. for lpawu and for spin",ispden," =",paw_ij_all(iatom)%nocctot(ispden)
636              call wrtout(unitfi,msg,'COLL')
637            end do
638            localm=paw_ij_all(iatom)%nocctot(2)-paw_ij_all(iatom)%nocctot(1)
639            write(msg,fmt='(a,i4,a,2x,f12.6)') " => On atom",iatom,&
640 &           ",  local Mag. for lpawu is  ",localm
641            call wrtout(unitfi,msg,'COLL')
642          end if
643          if(ndij==4) then
644            ntot=paw_ij_all(iatom)%nocctot(1)
645            mx=paw_ij_all(iatom)%nocctot(2)
646            my=paw_ij_all(iatom)%nocctot(3)
647            mz=paw_ij_all(iatom)%nocctot(4)
648            mnorm=sqrt(mx*mx+my*my+mz*mz)
649            write(msg,'(a,i4,a,2x,e15.8)') " => On atom",iatom,", for  lpawu, local Mag. x is  ",mx
650            call wrtout(unitfi,msg,'COLL')
651            write(msg,'(14x,a,2x,e15.8)') "               local Mag. y is  ",my
652            call wrtout(unitfi,msg,'COLL')
653            write(msg,'(14x,a,2x,e15.8)') "               local Mag. z is  ",mz
654            call wrtout(unitfi,msg,'COLL')
655            write(msg,'(14x,a,2x,e15.8)') "               norm of Mag. is  ",mnorm
656            call wrtout(unitfi,msg,'COLL')
657            write(msg,fmt='(8x,a,2x,f10.5)') " (along mag axis)    occ. for majority spin is = ",&
658 &           half*(ntot+mnorm)
659            call wrtout(unitfi,msg,'COLL')
660            write(msg,fmt='(8x,a,2x,f10.5)') " (along mag axis)    occ. for minority spin is = ",&
661 &           half*(ntot-mnorm)
662            call wrtout(unitfi,msg,'COLL')
663          end if
664          write(msg,'(3a)') ch10," == Occupation matrix for correlated orbitals:",ch10
665          call wrtout(unitfi,msg,'COLL')
666          do ispden=1,ndij
667            if (nspden==1.and.ndij/=4.and.(cplex_dij==1)) write(msg,fmt='(a)') " Up component only..."
668            if (nspden==2) write(msg,fmt='(a,i3)')" Occupation matrix for spin",ispden
669            if (ndij==4.or.(cplex_dij==2)) &
670 &           write(msg,fmt='(2a)')  " Occupation matrix for component ",trim(dspin1(ispden+2*(ndij/4)))
671            call wrtout(unitfi,msg,'COLL')
672            do im1=1,ll*2+1
673              if(cplex_dij==1)&
674 &             write(msg,'(12(1x,9(1x,f10.5)))') (paw_ij_all(iatom)%noccmmp(1,im1,im2,ispden),im2=1,ll*2+1)
675              if(cplex_dij==2)&
676 &             write(msg,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
677 &             (paw_ij_all(iatom)%noccmmp(:,im1,im2,ispden),im2=1,ll*2+1)
678              call wrtout(unitfi,msg,'COLL')
679            end do
680            write(msg,'(2a)') ch10,' '
681            call wrtout(unitfi,msg,'COLL')
682          end do
683 !        Transformation matrices: real->complex spherical harmonics
684          if(paw_ij_all(iatom)%ndij==4) then
685            ABI_ALLOCATE(noccmmp_ylm,(2*ll+1,2*ll+1,paw_ij_all(iatom)%ndij))
686            noccmmp_ylm=czero
687            ABI_ALLOCATE(noccmmp_slm,(2*ll+1,2*ll+1,paw_ij_all(iatom)%ndij))
688            noccmmp_slm=czero
689 !          Go from real notation for complex noccmmp to complex notation in noccmmp_slm
690            noccmmp_slm(:,:,:)=cmplx(paw_ij_all(iatom)%noccmmp(1,:,:,:),paw_ij_all(iatom)%noccmmp(2,:,:,:))
691            ii=std_out;if (unitfi==ab_out) ii=-1
692            call mat_slm2ylm(ll,noccmmp_slm,noccmmp_ylm,paw_ij_all(iatom)%ndij,&
693 &           1,1,dtset%pawprtvol,ii,'COLL') ! optspin=1 because up spin are first
694            do ispden=1,paw_ij_all(iatom)%ndij
695              write(msg,'(3a)') ch10,&
696 &             "== Occupation matrix in the complex harmonics basis for component ",&
697 &             trim(dspin1(ispden+2*(paw_ij_all(iatom)%ndij/4)))
698              call wrtout(unitfi,msg,'COLL')
699              do im1=1,ll*2+1
700                write(msg,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))') &
701 &               (noccmmp_ylm(im1,im2,ispden),im2=1,ll*2+1)
702                call wrtout(unitfi,msg,'COLL')
703              end do
704            end do
705            write(msg,'(a)') ch10
706            call wrtout(unitfi,msg,'COLL')
707            if (dtset%pawspnorb>0) then
708              ABI_ALLOCATE(noccmmp_jmj,(2*(2*ll+1),2*(2*ll+1)))
709              noccmmp_jmj=czero
710              ii=std_out;if (unitfi==ab_out) ii=-1
711              call mat_mlms2jmj(ll,noccmmp_ylm,noccmmp_jmj,paw_ij_all(iatom)%ndij,&
712 &             1,1,dtset%pawprtvol,-1,'COLL') !  optspin=1: up spin are first
713              write(msg,'(3a)') ch10,"== Occupation matrix in the J (= L-1/2, L+1/2) and M_J basis"
714              call wrtout(unitfi,msg,'COLL')
715              do im1=1,2*(ll*2+1)
716                write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') &
717 &               (noccmmp_jmj(im1,im2),im2=1,2*(ll*2+1))
718                call wrtout(unitfi,msg,'COLL')
719              end do
720              write(msg,'(a)') ch10
721              call wrtout(unitfi,msg,'COLL')
722              ABI_DEALLOCATE(noccmmp_jmj)
723            end if ! pawspnorb
724            ABI_DEALLOCATE(noccmmp_ylm)
725            ABI_DEALLOCATE(noccmmp_slm)
726          end if ! ndij==4
727        end if ! ((ll>=0).and.(pawtab(itypat)%usepawu>0))
728      end do
729    end do
730  end if
731 
732 !Exact exchange: print out occupations for correlated orbitals
733 !-------------------------------------------------------------
734  if (useexexch.and.ipositron/=1.and.me_atom==0) then
735    nspden=paw_ij_all(1)%nspden;nsppol=paw_ij_all(1)%nsppol;ndij=paw_ij_all(1)%ndij
736    do iatom=1,dtset%natom
737      itypat=dtset%typat(iatom);ll=pawtab(itypat)%lexexch
738      cplex_dij=paw_ij_all(iatom)%cplex_dij
739      if (ll>=0.and.pawtab(itypat)%useexexch>0) then
740        ABI_ALLOCATE(paw_ij_all(iatom)%noccmmp,(cplex_dij,2*ll+1,2*ll+1,ndij))
741        ABI_ALLOCATE(paw_ij_all(iatom)%nocctot,(nspden))
742      end if
743    end do
744    call setnoccmmp(1,0,rdum4,0,0,idum3,dtset%natom,dtset%natom,0,1,nsppol,0,dtset%ntypat,&
745 &   paw_ij_all,pawang_dum,dtset%pawprtvol,pawrhoij_all,pawtab,rdum2,idum1,dtset%typat,1,0)
746    do i_unitfi=1,2
747      unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
748      write(msg, '(3a)' ) &
749 &     ' ---------- Exact Exchange --------------------------------------------------- ',ch10
750      call wrtout(unitfi,msg,'COLL')
751      do iatom=1,dtset%natom
752        itypat=dtset%typat(iatom);ll=pawtab(itypat)%lexexch
753        cplex_dij=paw_ij_all(iatom)%cplex_dij
754        if ((ll>=0).and.(pawtab(itypat)%useexexch>0)) then
755          write(msg,fmt='(a,i5,a,i4,a)') " ====== For Atom",iatom,&
756 &         ", occupations for correlated orbitals. l =",ll,ch10
757          call wrtout(unitfi,msg,'COLL')
758          do ispden=1,ndij
759            if (nspden==1.and.ndij/=4) write(msg,fmt='(a)')   " Up component only..."
760            if (nspden==2) write(msg,fmt='(a,i3)')" Occupation matrix for spin",ispden
761            if (ndij==4) write(msg,fmt='(2a)')  " Occupation matrix for component ",&
762 &           trim(dspin2(ispden+2*(ndij/4)))
763            call wrtout(unitfi,msg,'COLL')
764            do im1=1,ll*2+1
765              if(cplex_dij==1)&
766 &             write(msg,'(12(1x,9(1x,f10.5)))')&
767 &             (paw_ij_all(iatom)%noccmmp(1,im1,im2,ispden),im2=1,ll*2+1)
768              if(cplex_dij==2)&
769 &             write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') &
770 &             (paw_ij_all(iatom)%noccmmp(:,im1,im2,ispden),im2=1,ll*2+1)
771              call wrtout(unitfi,msg,'COLL')
772            end do
773            call wrtout(unitfi,' ','COLL')
774          end do
775        end if
776      end do
777    end do
778    do iatom=1,dtset%natom
779      if (allocated(paw_ij_all(iatom)%noccmmp)) then
780        ABI_DEALLOCATE(paw_ij_all(iatom)%noccmmp)
781      end if
782      if (allocated(paw_ij_all(iatom)%nocctot)) then
783        ABI_DEALLOCATE(paw_ij_all(iatom)%nocctot)
784      end if
785    end do
786  end if
787 
788  msg=' '
789  call wrtout(ab_out,msg,'COLL')
790  call wrtout(std_out,msg,'COLL')
791 
792 !Destroy temporary stored atomic data
793  ABI_DEALLOCATE(jatom)
794  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
795  if (paral_atom) then
796    call pawrhoij_free(pawrhoij_all)
797    ABI_DATATYPE_DEALLOCATE(pawrhoij_all)
798    if (usepawu.or.useexexch) then
799      call paw_ij_free(paw_ij_all)
800      ABI_DATATYPE_DEALLOCATE(paw_ij_all)
801    end if
802  end if
803 
804  DBG_EXIT("COLL")
805 
806 end subroutine pawprt