LCOV - code coverage report
Current view: top level - shared/libpaw/src - m_pawcprj.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 53.5 % 1040 556
Test Date: 2026-09-21 19:39:32 Functions: 78.6 % 28 22

            Line data    Source code
       1              : !!****m* abinit/m_pawcprj
       2              : !! NAME
       3              : !!  m_pawcprj
       4              : !!
       5              : !! FUNCTION
       6              : !!  This module contains functions used to manipulate variables of
       7              : !!   structured datatype pawcprj_type.
       8              : !!   pawcprj_type variables are <p_lmn|Cnk> projected quantities,
       9              : !!   where |p_lmn> are non-local projectors
      10              : !!         |Cnk> are wave functions
      11              : !!
      12              : !! COPYRIGHT
      13              : !! Copyright (C) 2012-2026 ABINIT group (MT,JWZ)
      14              : !! This file is distributed under the terms of the
      15              : !! GNU General Public License, see ~abinit/COPYING
      16              : !! or http://www.gnu.org/copyleft/gpl.txt .
      17              : !!
      18              : !! NOTES
      19              : !!  FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
      20              : !!  please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
      21              : !!
      22              : !! SOURCE
      23              : 
      24              : #include "libpaw.h"
      25              : 
      26              : module m_pawcprj
      27              : 
      28              :  USE_DEFS
      29              :  USE_MSG_HANDLING
      30              :  USE_MPI_WRAPPERS
      31              :  USE_MEMORY_PROFILING
      32              : 
      33              :  use m_pawtab, only : pawtab_type
      34              : 
      35              :  implicit none
      36              : 
      37              :  private
      38              : !!***
      39              : 
      40              : !!****t* m_pawcprj/pawcprj_type
      41              : !! NAME
      42              : !! pawcprj_type
      43              : !!
      44              : !! FUNCTION
      45              : !! This structured datatype contains <p_lmn|Cnk> projected scalars and derivatives
      46              : !!             where |p_lmn> are non-local projectors for a given atom
      47              : !!                   |Cnk> is a wave function
      48              : !! Used only for PAW calculations.
      49              : !!
      50              : !! SOURCE
      51              : 
      52              :  type,public :: pawcprj_type
      53              : 
      54              : !Integer scalars
      55              : 
      56              :   integer :: ncpgr=0
      57              :    ! Number of gradients of cp=<p_lmn|Cnk>
      58              : 
      59              :   integer :: nlmn=0
      60              :    ! Number of (l,m,n) non-local projectors
      61              : 
      62              : !Real (real(dp)) arrays
      63              : 
      64              :   real(dp), allocatable :: cp (:,:)
      65              :    ! cp(2,nlmn)
      66              :    ! <p_lmn|Cnk> projected scalars for a given atom and wave function
      67              : 
      68              :   real(dp), allocatable :: dcp (:,:,:)
      69              :    ! dcp(2,ncpgr,nlmn)
      70              :    ! derivatives of <p_lmn|Cnk> projected scalars for a given atom and wave function
      71              : 
      72              :  end type pawcprj_type
      73              : 
      74              : !public procedures.
      75              :  public :: pawcprj_alloc          ! Allocation
      76              :  public :: pawcprj_free           ! Deallocation
      77              :  public :: pawcprj_set_zero       ! Set to zero all arrays in a cprj datastructure
      78              :  public :: pawcprj_copy           ! Copy a cprj datastructure into another
      79              :  public :: pawcprj_axpby          ! cprjy(:,:) <- alpha.cprjx(:,:)+beta.cprjy(:,:)
      80              :  public :: pawcprj_zaxpby         ! cprjy(:,:) <- alpha.cprjx(:,:)+beta.cprjy(:,:), alpha and beta are COMPLEX scalars
      81              :  public :: pawcprj_projbd         ! cprjy(:,:) <- alpha.cprjx(:,:)+beta.cprjy(:,:), alpha and beta are COMPLEX scalars
      82              :  public :: pawcprj_conjg          ! cprj(:,:) <- conjugate(cprj(:,:))
      83              :  public :: pawcprj_symkn          ! construct cprj from that at a symmetry related k point
      84              :  public :: pawcprj_lincom         ! Compute a LINear COMbination of cprj datastructure:
      85              :  public :: pawcprj_output         ! Output a cprj. Useful for debugging.
      86              :  public :: pawcprj_get            ! Read the cprj for a given k-point from memory or from a temporary file
      87              :  public :: pawcprj_put            ! Write the cprj for a given set of (n,k) into memory or into a temporary file
      88              :  public :: pawcprj_reorder        ! Change the order of a cprj datastructure
      89              :  public :: pawcprj_mpi_allgather  ! Perform MPI_ALLGATHER on a pawcprj_type inside a MPI communicator.
      90              :  public :: pawcprj_bcast          ! Broadcast a pawcprj_type from master to all nodes inside a MPI communicator.
      91              :  public :: pawcprj_transpose      ! Transpose a cprj datastructure FOR A GIVEN (K,SPIN)
      92              :  public :: pawcprj_gather_spin    ! Collect spin distributed cprjs.
      93              :  public :: pawcprj_mpi_exch       ! Exchange a pawcprj_type between two processors inside a MPI communicator.
      94              :  public :: pawcprj_mpi_send       ! Send a pawcprj_type inside a MPI communicator.
      95              :  public :: pawcprj_mpi_recv       ! Receive a pawcprj_type inside a MPI communicator.
      96              :  public :: pawcprj_mpi_sum        ! Perform MPI_SUM on a pawcprj_type inside a MPI communicator.
      97              :  public :: pawcprj_getdim         ! Returns the number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom.
      98              :  public :: paw_overlap            ! Compute the onsite contribution to the overlap between two states.
      99              :  public :: pawcprj_pack           ! Copy data from a cprj to a simple real buffer
     100              :  public :: pawcprj_unpack         ! Copy data from a simple real buffer to a cprj
     101              : !!***
     102              : 
     103              : CONTAINS
     104              : 
     105              : !===========================================================
     106              : !!***
     107              : 
     108              : !!****f* m_pawcprj/pawcprj_alloc
     109              : !! NAME
     110              : !! pawcprj_alloc
     111              : !!
     112              : !! FUNCTION
     113              : !! Allocation of a cprj datastructure
     114              : !!
     115              : !! INPUTS
     116              : !!  ncpgr=number of gradients to be allocated
     117              : !!  nlmn(:)=sizes of cprj%cp
     118              : !!
     119              : !! SIDE EFFECTS
     120              : !!  cprj(:,:) <type(pawcprj_type)>= cprj datastructure
     121              : !!
     122              : !! SOURCE
     123              : 
     124      4211900 :  subroutine pawcprj_alloc(cprj,ncpgr,nlmn)
     125              : 
     126              : !Arguments ------------------------------------
     127              : !scalars
     128              :  integer,intent(in) :: ncpgr
     129              : !arrays
     130              :  integer,intent(in) :: nlmn(:)
     131              :  type(pawcprj_type),intent(inout) :: cprj(:,:)
     132              : 
     133              : !Local variables-------------------------------
     134              : !scalars
     135              :  integer :: ii,jj,n1dim,n2dim,nn
     136              :  character(len=500) :: msg
     137              : 
     138              : ! *************************************************************************
     139              : 
     140      4211900 :  n1dim=size(cprj,dim=1);n2dim=size(cprj,dim=2);nn=size(nlmn,dim=1)
     141      4211900 :  if (nn/=n1dim) then
     142            0 :    write(msg,*) 'wrong sizes (pawcprj_alloc)! :',nn,n1dim
     143            0 :    LIBPAW_ERROR(msg)
     144              :  end if
     145              : 
     146     13374647 :  do jj=1,n2dim
     147     32548576 :    do ii=1,n1dim
     148     19173929 :      if (allocated(cprj(ii,jj)%cp)) then
     149            0 :        LIBPAW_DEALLOCATE(cprj(ii,jj)%cp)
     150              :      end if
     151     19173929 :      if (allocated(cprj(ii,jj)%dcp)) then
     152            0 :        LIBPAW_DEALLOCATE(cprj(ii,jj)%dcp)
     153              :      end if
     154     19173929 :      nn=nlmn(ii)
     155     19173929 :      cprj(ii,jj)%nlmn=nn
     156     57521787 :      LIBPAW_ALLOCATE(cprj(ii,jj)%cp,(2,nn))
     157    550538696 :      cprj(ii,jj)%cp=zero
     158     19173929 :      cprj(ii,jj)%ncpgr=ncpgr
     159     28336676 :      if (ncpgr>0) then
     160      8669316 :        LIBPAW_ALLOCATE(cprj(ii,jj)%dcp,(2,ncpgr,nn))
     161    110329299 :        cprj(ii,jj)%dcp=zero
     162              :      end if
     163              :    end do
     164              :  end do
     165              : 
     166      4211900 : end subroutine pawcprj_alloc
     167              : !!***
     168              : 
     169              : !----------------------------------------------------------------------
     170              : 
     171              : !!****f* m_pawcprj/pawcprj_free
     172              : !! NAME
     173              : !! pawcprj_free
     174              : !!
     175              : !! FUNCTION
     176              : !! Deallocation of a cprj datastructure
     177              : !!
     178              : !! SIDE EFFECTS
     179              : !!  cprj(:,:) <type(pawcprj_type)>= cprj datastructure
     180              : !!
     181              : !! SOURCE
     182              : 
     183      4220061 :  subroutine pawcprj_free(cprj)
     184              : 
     185              : !Arguments ------------------------------------
     186              : !scalars
     187              : !arrays
     188              :  type(pawcprj_type),intent(inout) :: cprj(:,:)
     189              : 
     190              : !Local variables-------------------------------
     191              : !scalars
     192              :  integer :: ii,jj,n1dim,n2dim
     193              : 
     194              : ! *************************************************************************
     195              : 
     196      4220061 :  n1dim=size(cprj,dim=1);n2dim=size(cprj,dim=2)
     197              : 
     198     13385528 :  do jj=1,n2dim
     199     32564609 :    do ii=1,n1dim
     200     19179081 :      if (allocated(cprj(ii,jj)%cp))  then
     201     19173033 :        LIBPAW_DEALLOCATE(cprj(ii,jj)%cp)
     202              :      end if
     203     28344548 :      if (allocated(cprj(ii,jj)%dcp))  then
     204      2167329 :        LIBPAW_DEALLOCATE(cprj(ii,jj)%dcp)
     205              :      end if
     206              :    end do
     207              :  end do
     208              : 
     209      4220061 : end subroutine pawcprj_free
     210              : !!***
     211              : 
     212              : !----------------------------------------------------------------------
     213              : 
     214              : !!****f* m_pawcprj/pawcprj_set_zero
     215              : !! NAME
     216              : !! pawcprj_set_zero
     217              : !!
     218              : !! FUNCTION
     219              : !! Set to zero all arrays in a cprj datastructure
     220              : !!
     221              : !! SIDE EFFECTS
     222              : !!  cprj(:,:) <type(pawcprj_type)>= cprj datastructure
     223              : !!
     224              : !! SOURCE
     225              : 
     226       621523 :  subroutine pawcprj_set_zero(cprj)
     227              : 
     228              : !Arguments ------------------------------------
     229              : !scalars
     230              : !arrays
     231              :  type(pawcprj_type),intent(inout) :: cprj(:,:)
     232              : 
     233              : !Local variables-------------------------------
     234              : !scalars
     235              :  integer :: ii,jj,n1dim,n2dim
     236              : 
     237              : ! *************************************************************************
     238              : 
     239       621523 :  n1dim=size(cprj,dim=1);n2dim=size(cprj,dim=2)
     240              : 
     241      1618199 :  do jj=1,n2dim
     242      3699733 :    do ii=1,n1dim
     243     56143592 :      if (cprj(ii,jj)%nlmn>0)  cprj(ii,jj)%cp(:,:)=zero
     244     24263918 :      if (cprj(ii,jj)%ncpgr>0) cprj(ii,jj)%dcp(:,:,:)=zero
     245              :    end do
     246              :  end do
     247              : 
     248       621523 : end subroutine pawcprj_set_zero
     249              : !!***
     250              : 
     251              : !----------------------------------------------------------------------
     252              : 
     253              : !!****f* m_pawcprj/pawcprj_copy
     254              : !! NAME
     255              : !! pawcprj_copy
     256              : !!
     257              : !! FUNCTION
     258              : !! Copy a cprj datastructure into another
     259              : !!
     260              : !! INPUTS
     261              : !!  icpgr= (optional argument) if present, only component icpgr of
     262              : !!         input cprj gradient is copied into output cprj
     263              : !!         Not used if cprj(:,:)%ncpgr<icpgr
     264              : !!         -1 only copy cp
     265              : !!  cprj_in(:,:) <type(pawcprj_type)>= input cprj datastructure
     266              : !!
     267              : !! OUTPUT
     268              : !!  cprj_out(:,:) <type(pawcprj_type)>= output cprj datastructure
     269              : !!
     270              : !! NOTES
     271              : !!  MG: What about an option to report a pointer to cprj_in?
     272              : !!
     273              : !! SOURCE
     274              : 
     275      5550184 :  subroutine pawcprj_copy(cprj_in,cprj_out,&
     276              : &                    icpgr) ! optional argument
     277              : 
     278              : !Arguments ------------------------------------
     279              : !scalars
     280              :  integer,intent(in),optional :: icpgr
     281              : !arrays
     282              :  type(pawcprj_type),intent(in) :: cprj_in(:,:)
     283              :  type(pawcprj_type),intent(inout) :: cprj_out(:,:)
     284              : 
     285              : !Local variables-------------------------------
     286              : !scalars
     287              :  integer :: ii,jj,kk,n1dim_in,n1dim_out,n2dim_in,n2dim_out,ncpgr_in,ncpgr_out,nlmn
     288              :  logical :: has_icpgr,copy_dcp
     289              :  character(len=500) :: msg
     290              : 
     291              : ! *************************************************************************
     292              : 
     293      5550184 :  n1dim_in=size(cprj_in,dim=1); n1dim_out=size(cprj_out,dim=1)
     294      5550184 :  n2dim_in=size(cprj_in,dim=2); n2dim_out=size(cprj_out,dim=2)
     295      5550184 :  ncpgr_in=cprj_in(1,1)%ncpgr;  ncpgr_out=cprj_out(1,1)%ncpgr
     296              : 
     297      5550184 :  if (n1dim_in/=n1dim_out) then
     298            0 :    write(msg,'(a,2(1x,i0))')" Error in pawcprj_copy: n1 wrong sizes ",n1dim_in,n1dim_out
     299            0 :    LIBPAW_ERROR(msg)
     300              :  end if
     301      5550184 :  if (n2dim_in/=n2dim_out) then
     302            0 :    write(msg,'(a,2(1x,i0))')" Error in pawcprj_copy: n2 wrong sizes ",n2dim_in,n2dim_out
     303            0 :    LIBPAW_ERROR(msg)
     304              :  end if
     305      5550184 :  if (ncpgr_in<ncpgr_out)  then
     306            0 :    write(msg,'(a,2(1x,i0))')" Error in pawcprj_copy: ncpgr wrong sizes ",ncpgr_in,ncpgr_out
     307            0 :    LIBPAW_ERROR(msg)
     308              :  end if
     309              : 
     310              : !Check if icgr is present and if dcp have to be copy
     311      5550184 :  has_icpgr=present(icpgr)
     312      5550184 :  copy_dcp = .TRUE.
     313      5550184 :  if(has_icpgr)then
     314       309168 :    copy_dcp = icpgr>=0
     315              :  end if
     316              : 
     317     20515720 :  do jj=1,n2dim_in
     318     50023004 :    do ii=1,n1dim_in
     319     29507284 :      nlmn=cprj_in(ii,jj)%nlmn
     320     29507284 :      cprj_out(ii,jj)%nlmn =nlmn
     321    282354026 :      do kk=1,nlmn
     322    743150902 :        cprj_out(ii,jj)%cp(1:2,kk)=cprj_in(ii,jj)%cp(1:2,kk)
     323              :      end do
     324              :    end do
     325              :  end do
     326              : 
     327      5550184 :  if (ncpgr_in>0.and.copy_dcp) then
     328      1491805 :    if (has_icpgr) has_icpgr=(ncpgr_out>0.and.icpgr>0.or.icpgr<=ncpgr_in)
     329              : 
     330              :    if (has_icpgr) then
     331       400896 :      do jj=1,n2dim_in
     332       838464 :        do ii=1,n1dim_in
     333       437568 :          nlmn=cprj_in(ii,jj)%nlmn
     334      4325520 :          do kk=1,nlmn
     335     11500080 :            cprj_out(ii,jj)%dcp(1:2,1,kk)=cprj_in(ii,jj)%dcp(1:2,icpgr,kk)
     336              :          end do
     337              :        end do
     338              :      end do
     339              :    else
     340      1291357 :      if (ncpgr_out>=ncpgr_in) then
     341      4383459 :        do jj=1,n2dim_in
     342      9879677 :          do ii=1,n1dim_in
     343      5496218 :            nlmn=cprj_in(ii,jj)%nlmn
     344     52835827 :            do kk=1,nlmn
     345   1113446776 :              cprj_out(ii,jj)%dcp(1:2,1:ncpgr_in,kk)=cprj_in(ii,jj)%dcp(1:2,1:ncpgr_in,kk)
     346              :            end do
     347              :          end do
     348              :        end do
     349              :      end if
     350              :    end if
     351              :  end if
     352              : 
     353      5550184 : end subroutine pawcprj_copy
     354              : !!***
     355              : 
     356              : !----------------------------------------------------------------------
     357              : 
     358              : !!****f* m_pawcprj/pawcprj_axpby
     359              : !! NAME
     360              : !! pawcprj_axpby
     361              : !!
     362              : !! FUNCTION
     363              : !! Apply AXPBY (blas-like) operation with 2 cprj datastructures:
     364              : !!  cprjy(:,:) <- alpha.cprjx(:,:)+beta.cprjy(:,:)
     365              : !!  alpha and beta are REAL scalars
     366              : !!
     367              : !! INPUTS
     368              : !!  alpha,beta= alpha,beta REAL factors
     369              : !!  cprjx(:,:) <type(pawcprj_type)>= input cprjx datastructure
     370              : !!
     371              : !! SIDE EFFECTS
     372              : !!  cprjy(:,:) <type(pawcprj_type)>= input/output cprjy datastructure
     373              : !!
     374              : !! SOURCE
     375              : 
     376      2148002 :  subroutine pawcprj_axpby(alpha,beta,cprjx,cprjy)
     377              : 
     378              : !Arguments ------------------------------------
     379              : !scalars
     380              :  real(dp),intent(in) :: alpha,beta
     381              : !arrays
     382              :  type(pawcprj_type),intent(in) :: cprjx(:,:)
     383              :  type(pawcprj_type),intent(inout) :: cprjy(:,:)
     384              : 
     385              : !Local variables-------------------------------
     386              : !scalars
     387              :  integer :: ii,jj,kk,n1dimx,n1dimy,n2dimx,n2dimy,ncpgrx,ncpgry,nlmn
     388              :  character(len=500) :: msg
     389              : 
     390              : ! *************************************************************************
     391              : 
     392      2148002 :  n1dimy=size(cprjy,dim=1);n2dimy=size(cprjy,dim=2);ncpgry=cprjy(1,1)%ncpgr
     393      2148002 :  if (abs(alpha)>tol16) then
     394      2147950 :    n1dimx=size(cprjx,dim=1);n2dimx=size(cprjx,dim=2);ncpgrx=cprjx(1,1)%ncpgr
     395      2147950 :    msg = ""
     396      2147950 :    if (n1dimx/=n1dimy) msg = TRIM(msg)//"Error in pawcprj_axpby: n1 wrong sizes !"//ch10
     397      2147950 :    if (n2dimx/=n2dimy) msg = TRIM(msg)//"Error in pawcprj_axpby: n2 wrong sizes !"//ch10
     398      2147950 :    if (ncpgrx/=ncpgry) msg = TRIM(msg)//"Error in pawcprj_axpby: ncpgr wrong sizes !"//ch10
     399      2147950 :    if (LEN_TRIM(msg) > 0) then
     400            0 :      LIBPAW_ERROR(msg)
     401              :    end if
     402              :  else
     403              :    n1dimx=0;n2dimx=0;ncpgrx=0
     404              :  end if
     405              : 
     406      2148002 :  if (abs(alpha)<=tol16) then
     407          104 :    do jj=1,n2dimy
     408          158 :      do ii=1,n1dimy
     409           54 :        nlmn=cprjy(ii,jj)%nlmn
     410          538 :        do kk=1,nlmn
     411         1350 :          cprjy(ii,jj)%cp(1:2,kk)=beta*cprjy(ii,jj)%cp(1:2,kk)
     412              :        end do
     413              :      end do
     414              :    end do
     415           52 :    if (ncpgry>0) then
     416            0 :      do jj=1,n2dimy
     417            0 :        do ii=1,n1dimy
     418            0 :          nlmn=cprjy(ii,jj)%nlmn
     419            0 :          do kk=1,nlmn
     420            0 :            cprjy(ii,jj)%dcp(1:2,1:ncpgry,kk)=beta*cprjy(ii,jj)%dcp(1:2,1:ncpgry,kk)
     421              :          end do
     422              :        end do
     423              :      end do
     424              :    end if
     425      2147950 :  else if (abs(beta)<=tol16) then
     426       562728 :    do jj=1,n2dimx
     427      1256172 :      do ii=1,n1dimx
     428       693444 :        nlmn=cprjx(ii,jj)%nlmn
     429       693444 :        cprjy(ii,jj)%nlmn=nlmn
     430      6746760 :        do kk=1,nlmn
     431     18009300 :          cprjy(ii,jj)%cp(1:2,kk)=alpha*cprjx(ii,jj)%cp(1:2,kk)
     432              :        end do
     433              :      end do
     434              :    end do
     435       281364 :    if (ncpgrx>0) then
     436            0 :      do jj=1,n2dimx
     437            0 :        do ii=1,n1dimx
     438            0 :          nlmn=cprjx(ii,jj)%nlmn
     439            0 :          do kk=1,nlmn
     440            0 :            cprjy(ii,jj)%dcp(1:2,1:ncpgrx,kk)=alpha*cprjx(ii,jj)%dcp(1:2,1:ncpgrx,kk)
     441              :          end do
     442              :        end do
     443              :      end do
     444              :    end if
     445              :  else  ! alpha/=0 and beta/=0
     446      4778950 :    do jj=1,n2dimx
     447     10982565 :      do ii=1,n1dimx
     448      6203615 :        nlmn=cprjx(ii,jj)%nlmn
     449      6203615 :        cprjy(ii,jj)%nlmn=nlmn
     450     62261539 :        do kk=1,nlmn
     451              :          cprjy(ii,jj)%cp(1:2,kk)=alpha*cprjx(ii,jj)%cp(1:2,kk) &
     452    165640295 : &         +beta *cprjy(ii,jj)%cp(1:2,kk)
     453              :        end do
     454              :      end do
     455              :    end do
     456      1866586 :    if (ncpgrx>0) then
     457            0 :      do jj=1,n2dimx
     458            0 :        do ii=1,n1dimx
     459            0 :          nlmn=cprjx(ii,jj)%nlmn
     460            0 :          do kk=1,nlmn
     461              :            cprjy(ii,jj)%dcp(1:2,1:ncpgrx,kk)=alpha*cprjx(ii,jj)%dcp(1:2,1:ncpgrx,kk) &
     462            0 : &           +beta *cprjy(ii,jj)%dcp(1:2,1:ncpgrx,kk)
     463              :          end do
     464              :        end do
     465              :      end do
     466              :    end if
     467              :  end if
     468              : 
     469      2148002 : end subroutine pawcprj_axpby
     470              : !!***
     471              : 
     472              : !----------------------------------------------------------------------
     473              : 
     474              : !!****f* m_pawcprj/pawcprj_zaxpby
     475              : !! NAME
     476              : !! pawcprj_zaxpby
     477              : !!
     478              : !! FUNCTION
     479              : !! Apply ZAXPBY (blas-like) operation with 2 cprj datastructures:
     480              : !!  cprjy(:,:) <- alpha.cprjx(:,:)+beta.cprjy(:,:)
     481              : !!  alpha and beta are COMPLEX scalars
     482              : !!
     483              : !! INPUTS
     484              : !!  alpha(2),beta(2)= alpha,beta COMPLEX factors
     485              : !!  cprjx(:,:) <type(pawcprj_type)>= input cprjx datastructure
     486              : !!
     487              : !! SIDE EFFECTS
     488              : !!  cprjy(:,:) <type(pawcprj_type)>= input/output cprjy datastructure
     489              : !!
     490              : !! SOURCE
     491              : 
     492      1107685 :  subroutine pawcprj_zaxpby(alpha,beta,cprjx,cprjy)
     493              : 
     494              : !Arguments ------------------------------------
     495              : !scalars
     496              :  real(dp),intent(in) :: alpha(2),beta(2)
     497              : !arrays
     498              :  type(pawcprj_type),intent(in) :: cprjx(:,:)
     499              :  type(pawcprj_type),intent(inout) :: cprjy(:,:)
     500              : 
     501              : !Local variables-------------------------------
     502              : !scalars
     503              :  integer :: ii,jj,kk,ll,n1dimx,n1dimy,n2dimx,n2dimy,ncpgrx,ncpgry,nlmn
     504              :  real(dp) :: cp1,cp2,norma,normb
     505              :  character(len=500) :: msg
     506              : 
     507              : ! *************************************************************************
     508              : 
     509      1107685 :  norma=alpha(1)**2+alpha(2)**2
     510      1107685 :  normb=beta(1) **2+beta(2) **2
     511      1107685 :  n1dimy=size(cprjy,dim=1);n2dimy=size(cprjy,dim=2);ncpgry=cprjy(1,1)%ncpgr
     512      1107685 :  if (norma>tol16*tol16) then
     513      1107606 :    n1dimx=size(cprjx,dim=1);n2dimx=size(cprjx,dim=2);ncpgrx=cprjx(1,1)%ncpgr
     514      1107606 :    msg = ""
     515      1107606 :    if (n1dimx/=n1dimy) msg = TRIM(msg)//"Error in pawcprj_zaxpby: n1 wrong sizes !"//ch10
     516      1107606 :    if (n2dimx/=n2dimy) msg = TRIM(msg)//"Error in pawcprj_zaxpby: n2 wrong sizes !"//ch10
     517      1107606 :    if (ncpgrx/=ncpgry) msg = TRIM(msg)//"Error in pawcprj_zaxpby: ncpgr wrong sizes !"//ch10
     518      1107606 :    if (LEN_TRIM(msg) > 0) then
     519            0 :      LIBPAW_ERROR(msg)
     520              :    end if
     521              :  end if
     522              : 
     523      1107685 :  if (norma<=tol16*tol16) then
     524          158 :    do jj=1,n2dimy
     525          438 :      do ii=1,n1dimy
     526          280 :        nlmn=cprjy(ii,jj)%nlmn
     527         2719 :        do kk=1,nlmn
     528         2360 :          cp1=beta(1)*cprjy(ii,jj)%cp(1,kk)-beta(2)*cprjy(ii,jj)%cp(2,kk)
     529         2360 :          cp2=beta(1)*cprjy(ii,jj)%cp(2,kk)+beta(2)*cprjy(ii,jj)%cp(1,kk)
     530         2360 :          cprjy(ii,jj)%cp(1,kk)=cp1
     531         2640 :          cprjy(ii,jj)%cp(2,kk)=cp2
     532              :        end do
     533              :      end do
     534              :    end do
     535           79 :    if (ncpgry>0) then
     536            0 :      do jj=1,n2dimy
     537            0 :        do ii=1,n1dimy
     538            0 :          nlmn=cprjy(ii,jj)%nlmn
     539            0 :          do kk=1,nlmn
     540            0 :            do ll=1,ncpgry
     541            0 :              cp1=beta(1)*cprjy(ii,jj)%dcp(1,ll,kk)-beta(2)*cprjy(ii,jj)%dcp(2,ll,kk)
     542            0 :              cp2=beta(1)*cprjy(ii,jj)%dcp(2,ll,kk)+beta(2)*cprjy(ii,jj)%dcp(1,ll,kk)
     543            0 :              cprjy(ii,jj)%dcp(1,ll,kk)=cp1
     544            0 :              cprjy(ii,jj)%dcp(2,ll,kk)=cp2
     545              :            end do
     546              :          end do
     547              :        end do
     548              :      end do
     549              :    end if
     550      1107606 :  else if (normb<=tol16*tol16) then
     551            0 :    do jj=1,n2dimx
     552            0 :      do ii=1,n1dimx
     553            0 :        nlmn=cprjx(ii,jj)%nlmn
     554            0 :        cprjy(ii,jj)%nlmn=nlmn
     555            0 :        do kk=1,nlmn
     556            0 :          cprjy(ii,jj)%cp(1,kk)=alpha(1)*cprjx(ii,jj)%cp(1,kk)-alpha(2)*cprjx(ii,jj)%cp(2,kk)
     557            0 :          cprjy(ii,jj)%cp(2,kk)=alpha(1)*cprjx(ii,jj)%cp(2,kk)+alpha(2)*cprjx(ii,jj)%cp(1,kk)
     558              :        end do
     559              :      end do
     560              :    end do
     561            0 :    if (ncpgrx>0) then
     562            0 :      do jj=1,n2dimx
     563            0 :        do ii=1,n1dimx
     564            0 :          nlmn=cprjx(ii,jj)%nlmn
     565            0 :          do kk=1,nlmn
     566              :            cprjy(ii,jj)%dcp(1,1:ncpgrx,kk)=alpha(1)*cprjx(ii,jj)%dcp(1,1:ncpgrx,kk) &
     567            0 : &           -alpha(2)*cprjx(ii,jj)%dcp(2,1:ncpgrx,kk)
     568              :            cprjy(ii,jj)%dcp(2,1:ncpgrx,kk)=alpha(1)*cprjx(ii,jj)%dcp(2,1:ncpgrx,kk) &
     569            0 : &           +alpha(2)*cprjx(ii,jj)%dcp(1,1:ncpgrx,kk)
     570              :          end do
     571              :        end do
     572              :      end do
     573              :    end if
     574              : ! else if (abs(beta(1)-one)<tol16.and.abs(beta(2))<tol16) then
     575              : !   do jj=1,n2dimx
     576              : !     do ii=1,n1dimx
     577              : !       nlmn=cprjx(ii,jj)%nlmn
     578              : !       cprjy(ii,jj)%nlmn =nlmn
     579              : !       do kk=1,nlmn
     580              : !         cp1=cprjy(ii,jj)%cp(1,kk)
     581              : !         cp2=cprjy(ii,jj)%cp(2,kk)
     582              : !         cp1=cp1+alpha(1)*cprjx(ii,jj)%cp(1,kk)-alpha(2)*cprjx(ii,jj)%cp(2,kk)
     583              : !         cp2=cp2+alpha(1)*cprjx(ii,jj)%cp(2,kk)+alpha(2)*cprjx(ii,jj)%cp(1,kk)
     584              : !         cprjy(ii,jj)%cp(1,kk)=cp1
     585              : !         cprjy(ii,jj)%cp(2,kk)=cp2
     586              : !       end do
     587              : !     end do
     588              : !   end do
     589              : !   if (ncpgrx>0) then
     590              : !     do jj=1,n2dimx
     591              : !       do ii=1,n1dimx
     592              : !         nlmn=cprjx(ii,jj)%nlmn
     593              : !         do kk=1,nlmn
     594              : !           do ll=1,ncpgrx
     595              : !             cp1=cprjy(ii,jj)%dcp(1,ll,kk)
     596              : !             cp2=cprjy(ii,jj)%dcp(2,ll,kk)
     597              : !             cp1=cp1+alpha(1)*cprjx(ii,jj)%dcp(1,ll,kk)-alpha(2)*cprjx(ii,jj)%dcp(2,ll,kk)
     598              : !             cp2=cp2+alpha(1)*cprjx(ii,jj)%dcp(2,ll,kk)+alpha(2)*cprjx(ii,jj)%dcp(1,ll,kk)
     599              : !             cprjy(ii,jj)%dcp(1,ll,kk)=cp1
     600              : !             cprjy(ii,jj)%dcp(2,ll,kk)=cp2
     601              : !           end do
     602              : !         end do
     603              : !       end do
     604              : !     end do
     605              : !   end if
     606              :  else
     607      2257512 :    do jj=1,n2dimx
     608      4605817 :      do ii=1,n1dimx
     609      2348305 :        nlmn=cprjx(ii,jj)%nlmn
     610      2348305 :        cprjy(ii,jj)%nlmn =nlmn
     611     24723531 :        do kk=1,nlmn
     612              :          cp1=alpha(1)*cprjx(ii,jj)%cp(1,kk)-alpha(2)*cprjx(ii,jj)%cp(2,kk) &
     613     21225320 : &         +beta(1) *cprjy(ii,jj)%cp(1,kk)-beta(2) *cprjy(ii,jj)%cp(2,kk)
     614              :          cp2=alpha(1)*cprjx(ii,jj)%cp(2,kk)+alpha(2)*cprjx(ii,jj)%cp(1,kk) &
     615     21225320 : &         +beta(1) *cprjy(ii,jj)%cp(2,kk)+beta(2) *cprjy(ii,jj)%cp(1,kk)
     616              : !         cp1=beta(1) *cprjy(ii,jj)%cp(1,kk)-beta(2) *cprjy(ii,jj)%cp(2,kk)
     617              : !         cp1=cp1+alpha(1)*cprjx(ii,jj)%cp(1,kk)-alpha(2)*cprjx(ii,jj)%cp(2,kk)
     618              : !         cp2=beta(1) *cprjy(ii,jj)%cp(2,kk)+beta(2) *cprjy(ii,jj)%cp(1,kk)
     619              : !         cp2=cp2+alpha(1)*cprjx(ii,jj)%cp(2,kk)+alpha(2)*cprjx(ii,jj)%cp(1,kk)
     620     21225320 :          cprjy(ii,jj)%cp(1,kk)=cp1
     621     23573625 :          cprjy(ii,jj)%cp(2,kk)=cp2
     622              :        end do
     623              :      end do
     624              :    end do
     625      1107606 :    if (ncpgrx>0) then
     626            0 :      do jj=1,n2dimx
     627            0 :        do ii=1,n1dimx
     628            0 :          nlmn=cprjx(ii,jj)%nlmn
     629            0 :          do kk=1,nlmn
     630            0 :            do ll=1,ncpgrx
     631              :              cp1=alpha(1)*cprjx(ii,jj)%dcp(1,ll,kk)-alpha(2)*cprjx(ii,jj)%dcp(2,ll,kk) &
     632            0 : &             +beta(1) *cprjy(ii,jj)%dcp(1,ll,kk)-beta(2) *cprjy(ii,jj)%dcp(2,ll,kk)
     633              :              cp2=alpha(1)*cprjx(ii,jj)%dcp(2,ll,kk)+alpha(2)*cprjx(ii,jj)%dcp(1,ll,kk) &
     634            0 : &             +beta(1) *cprjy(ii,jj)%dcp(2,ll,kk)+beta(2) *cprjy(ii,jj)%dcp(1,ll,kk)
     635              : !             cp1=beta(1) *cprjy(ii,jj)%dcp(1,ll,kk)-beta(2) *cprjy(ii,jj)%dcp(2,ll,kk)
     636              : !             cp1=cp1+alpha(1)*cprjx(ii,jj)%dcp(1,ll,kk)-alpha(2)*cprjx(ii,jj)%dcp(2,ll,kk)
     637              : !             cp2=beta(1) *cprjy(ii,jj)%dcp(2,ll,kk)+beta(2) *cprjy(ii,jj)%dcp(1,ll,kk)
     638              : !             cp2=cp2+alpha(1)*cprjx(ii,jj)%dcp(2,ll,kk)+alpha(2)*cprjx(ii,jj)%dcp(1,ll,kk)
     639            0 :              cprjy(ii,jj)%dcp(1,ll,kk)=cp1
     640            0 :              cprjy(ii,jj)%dcp(2,ll,kk)=cp2
     641              :            end do
     642              :          end do
     643              :        end do
     644              :      end do
     645              :    end if
     646              :  end if
     647              : 
     648      1107685 : end subroutine pawcprj_zaxpby
     649              : !!***
     650              : 
     651              : !----------------------------------------------------------------------
     652              : 
     653              : !!****f* m_pawcprj/pawcprj_projbd
     654              : !! NAME
     655              : !! pawcprj_projbd
     656              : !!
     657              : !! FUNCTION
     658              : !! Apply ZAXPBY (blas-like) operation with 2 cprj datastructures:
     659              : !!  cprjy(:,:) <- alpha.cprjx(:,:)+beta.cprjy(:,:)
     660              : !!  alpha and beta are COMPLEX scalars
     661              : !!
     662              : !! INPUTS
     663              : !!  alpha(2),beta(2)= alpha,beta COMPLEX factors
     664              : !!  cprjx(:,:) <type(pawcprj_type)>= input cprjx datastructure
     665              : !!
     666              : !! SIDE EFFECTS
     667              : !!  cprjy(:,:) <type(pawcprj_type)>= input/output cprjy datastructure
     668              : !!
     669              : !! SOURCE
     670              : 
     671            0 :  subroutine pawcprj_projbd(alpha,cprjx,cprjy)
     672              : 
     673              : !Arguments ------------------------------------
     674              : !scalars
     675              :  real(dp),intent(in) :: alpha(:,:)
     676              : !arrays
     677              :  type(pawcprj_type),intent(in) :: cprjx(:,:)
     678              :  type(pawcprj_type),intent(inout) :: cprjy(:,:)
     679              : 
     680              : !Local variables-------------------------------
     681              : !scalars
     682              :  integer :: ia,ii,jj,kk,ll,n1dima,n1dimx,n1dimy,n2dimx,n2dimy,n2dima,ncpgrx,ncpgry,nlmn
     683              :  real(dp) :: cp1,cp2,norma
     684              :  character(len=500) :: msg
     685              : 
     686              : ! *************************************************************************
     687              : 
     688            0 :  n1dimy=size(cprjy,dim=1);n2dimy=size(cprjy,dim=2);ncpgry=cprjy(1,1)%ncpgr
     689            0 :  n1dimx=size(cprjx,dim=1);n2dimx=size(cprjx,dim=2);ncpgrx=cprjx(1,1)%ncpgr
     690            0 :  n1dima=size(alpha,dim=1);n2dima=size(alpha,dim=2)
     691            0 :  msg = ""
     692            0 :  if (n1dima/=2) msg = TRIM(msg)//"Error in pawcprj_projbd: alpha n1 wrong sizes !"//ch10
     693            0 :  if (n1dimx/=n1dimy) msg = TRIM(msg)//"Error in pawcprj_projbd: n1 wrong sizes !"//ch10
     694            0 :  if (n2dimx/=n2dimy*n2dima) msg = TRIM(msg)//"Error in pawcprj_projbd: n2 wrong sizes !"//ch10
     695            0 :  if (ncpgrx/=ncpgry) msg = TRIM(msg)//"Error in pawcprj_projbd: ncpgr wrong sizes !"//ch10
     696            0 :  if (LEN_TRIM(msg) > 0) then
     697            0 :    LIBPAW_ERROR(msg)
     698              :  end if
     699              : 
     700            0 :  do ia=1,n2dima
     701            0 :    norma=alpha(1,ia)**2+alpha(2,ia)**2
     702            0 :    if (norma>tol16*tol16) then
     703            0 :      do jj=1,n2dimy
     704            0 :        do ii=1,n1dimx
     705            0 :          nlmn=cprjy(ii,jj)%nlmn
     706              :          cprjy(ii,jj)%nlmn =nlmn
     707            0 :          do kk=1,nlmn
     708            0 :            cp1=alpha(1,ia)*cprjx(ii,jj+(ia-1)*n2dimy)%cp(1,kk)-alpha(2,ia)*cprjx(ii,jj+(ia-1)*n2dimy)%cp(2,kk)
     709            0 :            cp2=alpha(1,ia)*cprjx(ii,jj+(ia-1)*n2dimy)%cp(2,kk)+alpha(2,ia)*cprjx(ii,jj+(ia-1)*n2dimy)%cp(1,kk)
     710            0 :            cprjy(ii,jj)%cp(1,kk)=cprjy(ii,jj)%cp(1,kk)+cp1
     711            0 :            cprjy(ii,jj)%cp(2,kk)=cprjy(ii,jj)%cp(2,kk)+cp2
     712              :          end do
     713              :        end do
     714              :      end do
     715            0 :      if (ncpgrx>0) then
     716            0 :        do jj=1,n2dimy
     717            0 :          do ii=1,n1dimx
     718            0 :            nlmn=cprjy(ii,jj)%nlmn
     719            0 :            do kk=1,nlmn
     720            0 :              do ll=1,ncpgrx
     721            0 :                cp1=alpha(1,ia)*cprjx(ii,jj+(ia-1)*n2dimy)%dcp(1,ll,kk)-alpha(2,ia)*cprjx(ii,jj+(ia-1)*n2dimy)%dcp(2,ll,kk)
     722            0 :                cp2=alpha(1,ia)*cprjx(ii,jj+(ia-1)*n2dimy)%dcp(2,ll,kk)+alpha(2,ia)*cprjx(ii,jj+(ia-1)*n2dimy)%dcp(1,ll,kk)
     723            0 :                cprjy(ii,jj)%dcp(1,ll,kk)=cprjy(ii,jj)%dcp(1,ll,kk)+cp1
     724            0 :                cprjy(ii,jj)%dcp(2,ll,kk)=cprjy(ii,jj)%dcp(2,ll,kk)+cp2
     725              :              end do
     726              :            end do
     727              :          end do
     728              :        end do
     729              :      end if
     730              :    end if
     731              :  end do
     732              : 
     733            0 : end subroutine pawcprj_projbd
     734              : !!***
     735              : 
     736              : !----------------------------------------------------------------------
     737              : 
     738              : !!****f* m_pawcprj/pawcprj_symkn
     739              : !! NAME
     740              : !! pawcprj_symkn
     741              : !!
     742              : !! FUNCTION
     743              : !! compute cprj for a given band and k point based on cprj at a symmetry-related
     744              : !! k point.
     745              : !!
     746              : !! INPUTS
     747              : !!  cprj_ikn (pawcprj_type) :: cprj for a single band and k point, typically a k point in the IBZ
     748              : !!  cprj_sym(4,nsym,natom) :: 1:3 shift, and 4 final atom, of symmetry isym operating on iatom
     749              : !!                            (S^{-1}(R - t) = r0 + L, see symatm.F90
     750              : !!  dimlmn(natom) :: ln dimension of each atom
     751              : !!  iband :: number of bands to treat, use -1 to treat all nband bands
     752              : !!  indlmn(6,lmnmax,ntypat) :: n,l,m dimensions for each atom type (see psps type)
     753              : !!  isym :: symmetry element used in current application
     754              : !!  itim :: 1 if time reversal also used, 0 else
     755              : !!  kpt(3) :: kpt vector used
     756              : !!  lmax :: max l value
     757              : !!  lmnmax :: max lmn value
     758              : !!  mband :: maximum number of bands
     759              : !!  natom :: number of atoms in cell
     760              : !!  nband :: number of bands in cprj_ikn
     761              : !!  nspinor :: number of spinors
     762              : !!  nsym :: total number of symmetry elements
     763              : !!  ntypat :: number of types of atoms
     764              : !!  typat(natom) :: type of each atom
     765              : !!  zarot(2*lmax+1,2*lmax+1,lmax+1,nsym) :: elements of rotation matrix for angular momentum states
     766              : !!                                          and symmetry operations. See m_paw_sphharm/setsym_ylm.
     767              : !!
     768              : !! OUTPUT
     769              : !!  cprj_fkn (pawcprj_type) :: cprj for a single band and k point where the k point is related to
     770              : !!    the input k point by a symmetry operation
     771              : !!
     772              : !! SIDE EFFECTS
     773              : !!
     774              : !! NOTES
     775              : !!  This routine is based on M. Giantomassi's doctoral dissertation, formula 7.77. It is not clear
     776              : !!  whether it is implemented correctly for nonsymmorphic symmetries.
     777              : !!
     778              : !! SOURCE
     779              : 
     780      1305428 :  subroutine pawcprj_symkn(cprj_fkn,cprj_ikn,cprj_sym,dimlmn,iband,indlmn,&
     781              : &                       isym,itim,kpt,lmax,lmnmax,mband,natom,nband,nspinor,nsym,ntypat,&
     782      1305428 : &                       typat,zarot)
     783              : 
     784              : !Arguments---------------------------
     785              : !scalars
     786              :  integer,intent(in) :: iband,isym,itim,lmax,lmnmax,mband
     787              :  integer,intent(in) :: natom,nband,nspinor,nsym,ntypat
     788              : 
     789              : !arrays
     790              :  integer,intent(in) :: cprj_sym(4,nsym,natom),dimlmn(natom)
     791              :  integer,intent(in) :: indlmn(6,lmnmax,ntypat),typat(natom)
     792              :  real(dp),intent(in) :: kpt(3)
     793              :  real(dp),intent(in) :: zarot(2*lmax+1,2*lmax+1,lmax+1,nsym)
     794              :  type(pawcprj_type),intent(in) :: cprj_ikn(natom,mband*nspinor)
     795              :  type(pawcprj_type),intent(inout) :: cprj_fkn(natom,mband*nspinor) !vz_i
     796              : 
     797              : !Local variables---------------------------
     798              : !scalars
     799              :  integer :: iatm,iatom, ibct, ibnd, ibsp, ibst, icpgr, iin, il, il0, im
     800              :  integer :: ilmn, iln, iln0, ilpm, indexi, ispinor, itypat, jatm,jatom, mm, nlmn
     801              :  real(dp) :: kdotL, phr, phi
     802              : !arrays
     803              :  real(dp) :: rl(3), t1(2), t2(2)
     804              : 
     805              : ! *************************************************************************
     806              : 
     807      1305428 :  if (iband == -1) then
     808      1302958 :    ibst = 1
     809      1302958 :    ibnd = nband
     810              :  else
     811              :    ibst = iband
     812              :    ibnd = iband
     813              :  end if
     814              : 
     815      3916284 :  do iatom = 1, natom
     816      2610856 :    iatm=iatom
     817      2610856 :    itypat = typat(iatom)
     818      2610856 :    nlmn = dimlmn(iatm)
     819      2610856 :    jatom = cprj_sym(4,isym,iatom)
     820      2610856 :    jatm=jatom
     821     10443424 :    rl(:) = cprj_sym(1:3,isym,iatom)
     822     10443424 :    kdotL = dot_product(rl,kpt)
     823      2610856 :    phr = cos(two_pi*kdotL)
     824      2610856 :    phi = sin(two_pi*kdotL)
     825              : 
     826      2610856 :    il0 = -1; iln0 = -1; indexi = 1
     827     24799952 :    do ilmn = 1, nlmn
     828              : 
     829     20883668 :      il = indlmn(1,ilmn,itypat)
     830     20883668 :      im = indlmn(2,ilmn,itypat)
     831     20883668 :      iin = indlmn(3,ilmn,itypat)
     832     20883668 :      iln = indlmn(5,ilmn,itypat)
     833     20883668 :      ilpm = 1 + il + im
     834     20883668 :      if (iln /= iln0) indexi = indexi + 2*il0 + 1
     835              : 
     836    104309320 :      do ibct = ibst, ibnd
     837              : 
     838    187734972 :        do ispinor = 1, nspinor
     839              : 
     840     83425652 :          ibsp = nspinor*(ibct-1) + ispinor
     841              : 
     842     83425652 :          t1(:) = zero
     843    291988192 :          do mm = 1, 2*il+1
     844    208562540 :            t1(1) = t1(1) + zarot(mm,ilpm,il+1,isym)*cprj_ikn(jatm,ibsp)%cp(1,indexi+mm)
     845    291988192 :            t1(2) = t1(2) + zarot(mm,ilpm,il+1,isym)*cprj_ikn(jatm,ibsp)%cp(2,indexi+mm)
     846              :          end do
     847     83425652 :          t2(1) = t1(1)*phr - t1(2)*phi
     848     83425652 :          t2(2) = t1(2)*phr + t1(1)*phi
     849              : 
     850     83425652 :          if (itim == 1) t2(2) = -t2(2)
     851              : 
     852     83425652 :          cprj_fkn(iatm,ibsp)%cp(1,ilmn) = t2(1)
     853     83425652 :          cprj_fkn(iatm,ibsp)%cp(2,ilmn) = t2(2)
     854              : 
     855              : ! do same transformations for gradients of cprj_ikn
     856              : ! note that ncpgr = 0 if no gradients present so this loop will not be executed
     857              : ! in this case
     858              : 
     859    339022084 :          do icpgr = 1, cprj_ikn(jatom,ibsp)%ncpgr
     860              :            t1(:) = zero
     861              : 
     862    602592960 :            do mm = 1, 2*il+1
     863    430422180 :              t1(1) = t1(1) + zarot(mm,ilpm,il+1,isym)*cprj_ikn(jatm,ibsp)%dcp(1,icpgr,indexi+mm)
     864    602592960 :              t1(2) = t1(2) + zarot(mm,ilpm,il+1,isym)*cprj_ikn(jatm,ibsp)%dcp(2,icpgr,indexi+mm)
     865              :            end do
     866              : 
     867    172170780 :            t2(1) = t1(1)*phr - t1(2)*phi
     868    172170780 :            t2(2) = t1(2)*phr + t1(1)*phi
     869              : 
     870    172170780 :            if (itim == 1) t2(2) = -t2(2)
     871              : 
     872    172170780 :            cprj_fkn(iatm,ibsp)%dcp(1,icpgr,ilmn) = t2(1)
     873    255596432 :            cprj_fkn(iatm,ibsp)%dcp(2,icpgr,ilmn) = t2(2)
     874              : 
     875              :          end do ! end loop over ncpgr
     876              : 
     877              :        end do ! end loop over nspinor
     878              : 
     879              :      end do ! end loop over bands
     880              : 
     881     23494524 :      il0 = il; iln0 = iln
     882              :    end do ! end loop over ilmn
     883              :  end do ! end loop over atoms
     884              : 
     885      1305428 :  end subroutine pawcprj_symkn
     886              : !!***
     887              : 
     888              : !----------------------------------------------------------------------
     889              : 
     890              : !!****f* m_pawcprj/pawcprj_conjg
     891              : !! NAME
     892              : !! pawcprj_conjg
     893              : !!
     894              : !! FUNCTION
     895              : !! conjugate a cprj datastructures:
     896              : !!  cprj(:,:) <- conjugate(cprj(:,:))
     897              : !!
     898              : !! INPUTS
     899              : !!
     900              : !! SIDE EFFECTS
     901              : !!  cprj(:,:) <type(pawcprj_type)>= input/output cprj datastructure
     902              : !!
     903              : !! SOURCE
     904              : 
     905            0 :  subroutine pawcprj_conjg(cprj)
     906              : 
     907              : !Arguments ------------------------------------
     908              : !scalars
     909              : !arrays
     910              :  type(pawcprj_type),intent(inout) :: cprj(:,:)
     911              : 
     912              : !Local variables-------------------------------
     913              : !scalars
     914              :  integer :: ii,jj,kk,n1dim,n2dim,ncpgr,nlmn
     915              : 
     916              : ! *************************************************************************
     917              : 
     918              : 
     919            0 :  n1dim=size(cprj,dim=1);n2dim=size(cprj,dim=2);ncpgr=cprj(1,1)%ncpgr
     920              : 
     921            0 :  do jj=1,n2dim
     922            0 :    do ii=1,n1dim
     923            0 :      nlmn=cprj(ii,jj)%nlmn
     924            0 :      do kk=1,nlmn
     925            0 :        cprj(ii,jj)%cp(2,kk)=-cprj(ii,jj)%cp(2,kk)
     926              :      end do
     927              :    end do
     928              :  end do
     929            0 :  if (ncpgr>0) then
     930            0 :    do jj=1,n2dim
     931            0 :      do ii=1,n1dim
     932            0 :        nlmn=cprj(ii,jj)%nlmn
     933            0 :        do kk=1,nlmn
     934            0 :          cprj(ii,jj)%dcp(2,1:ncpgr,kk)=-cprj(ii,jj)%dcp(2,1:ncpgr,kk)
     935              :        end do
     936              :      end do
     937              :    end do
     938              :  end if
     939              : 
     940            0 : end subroutine pawcprj_conjg
     941              : !!***
     942              : 
     943              : !----------------------------------------------------------------------
     944              : 
     945              : !!****f* m_pawcprj/pawcprj_lincom
     946              : !! NAME
     947              : !! pawcprj_lincom
     948              : !!
     949              : !! FUNCTION
     950              : !! Compute a LINear COMbination of cprj datastructure:
     951              : !!  cprj_out(:,:) <--- Sum_i [ alpha_i . cprj_i(:,:) ]
     952              : !!  alpha_i are COMPLEX scalars
     953              : !!
     954              : !! INPUTS
     955              : !!  alpha(2,nn)= alpha COMPLEX factors
     956              : !!  cprj_in(:,:) <type(pawcprj_type)>= input cprj_in datastructure
     957              : !!  nn= number of cprj involved in the linear combination
     958              : !!
     959              : !! OUTPUT
     960              : !!  cprj_out(:,:) <type(pawcprj_type)>= output cprj_out datastructure
     961              : !!
     962              : !! NOTES
     963              : !!  cprj_in and cprj_out must be dimensionned as cprj_in(n1,n2*nn) and cprj_in(n1,n2)
     964              : !!
     965              : !! SOURCE
     966              : 
     967       288348 :  subroutine pawcprj_lincom(alpha,cprj_in,cprj_out,nn)
     968              : 
     969              : !Arguments ------------------------------------
     970              : !scalars
     971              :  integer,intent(in) :: nn
     972              :  real(dp),intent(in) :: alpha(2,nn)
     973              : !arrays
     974              :  type(pawcprj_type),intent(in) :: cprj_in(:,:)
     975              :  type(pawcprj_type),intent(inout) :: cprj_out(:,:)
     976              : 
     977              : !Local variables-------------------------------
     978              : !scalars
     979              :  integer :: ii,in,jj,jn,kk,ll,n1in,n1out,n2in,n2out,ncpgrin,ncpgrout,nlmn
     980              :  real(dp) :: cp1,cp2
     981              :  character(len=500) :: msg
     982              : 
     983              : ! *************************************************************************
     984              : 
     985       288348 :  n1in=size(cprj_in,dim=1);n1out=size(cprj_out,dim=1)
     986       288348 :  n2in=size(cprj_in,dim=2);n2out=size(cprj_out,dim=2)
     987       288348 :  ncpgrin=cprj_in(1,1)%ncpgr;ncpgrout=cprj_out(1,1)%ncpgr
     988              : 
     989       288348 :  msg = ""
     990       288348 :  if (n1in/=n1out) msg = TRIM(msg)//"Bug in pawcprj_lincom: n1 wrong sizes!"//ch10
     991       288348 :  if (n2in/=n2out*nn) msg = TRIM(msg)//"Bug in pawcprj_lincom: n2 wrong sizes!"//ch10
     992       288348 :  if (ncpgrin/=ncpgrout) msg = TRIM(msg)//"Bug in pawcprj_lincom: ncpgr wrong sizes!"//ch10
     993       288348 :  if (LEN_TRIM(msg) > 0) then
     994            0 :    LIBPAW_ERROR(msg)
     995              :  end if
     996              : 
     997       576696 :  do jj=1,n2out
     998      1324932 :    do ii=1,n1out
     999       748236 :      nlmn=cprj_in(ii,jj)%nlmn
    1000       748236 :      cprj_out(ii,jj)%nlmn=nlmn
    1001     19379100 :      cprj_out(ii,jj)%cp(1:2,1:nlmn)=zero
    1002              :      jn=jj
    1003      6315696 :      do in=1,nn
    1004     50231448 :        do kk=1,nlmn
    1005              :          cp1=cprj_out(ii,jj)%cp(1,kk) &
    1006     44952336 : &         +alpha(1,in)*cprj_in(ii,jn)%cp(1,kk)-alpha(2,in)*cprj_in(ii,jn)%cp(2,kk)
    1007              :          cp2=cprj_out(ii,jj)%cp(2,kk) &
    1008     44952336 : &         +alpha(1,in)*cprj_in(ii,jn)%cp(2,kk)+alpha(2,in)*cprj_in(ii,jn)%cp(1,kk)
    1009     44952336 :          cprj_out(ii,jj)%cp(1,kk)=cp1
    1010     50231448 :          cprj_out(ii,jj)%cp(2,kk)=cp2
    1011              :        end do
    1012      6027348 :        jn=jn+n2out
    1013              :      end do
    1014              :    end do
    1015              :  end do
    1016              : 
    1017       288348 :  if (ncpgrin>0) then
    1018            0 :    do jj=1,n2out
    1019            0 :      do ii=1,n1out
    1020            0 :        nlmn=cprj_in(ii,jj)%nlmn
    1021            0 :        cprj_out(ii,jj)%dcp(1:2,1:ncpgrin,1:nlmn)=zero
    1022              :        jn=jj
    1023            0 :        do in=1,nn
    1024            0 :          do kk=1,nlmn
    1025            0 :            do ll=1,ncpgrin
    1026              :              cp1=cprj_out(ii,jj)%dcp(1,ll,kk) &
    1027              : &             +alpha(1,in)*cprj_in(ii,jn)%dcp(1,ll,kk) &
    1028            0 : &             -alpha(2,in)*cprj_in(ii,jn)%dcp(2,ll,kk)
    1029              :              cp2=cprj_out(ii,jj)%dcp(2,ll,kk) &
    1030              : &             +alpha(1,in)*cprj_in(ii,jn)%dcp(2,ll,kk) &
    1031            0 :              +alpha(2,in)*cprj_in(ii,jn)%dcp(1,ll,kk)
    1032            0 :              cprj_out(ii,jj)%dcp(1,ll,kk)=cp1
    1033            0 :              cprj_out(ii,jj)%dcp(2,ll,kk)=cp2
    1034              :            end do
    1035              :          end do
    1036            0 :          jn=jn+n2out
    1037              :        end do
    1038              :      end do
    1039              :    end do
    1040              :  end if
    1041              : 
    1042       288348 : end subroutine pawcprj_lincom
    1043              : !!***
    1044              : 
    1045              : !----------------------------------------------------------------------
    1046              : 
    1047              : !!****f* m_pawcprj/pawcprj_output
    1048              : !! NAME
    1049              : !! pawcprj_output
    1050              : !!
    1051              : !! FUNCTION
    1052              : !! Output a cprj. Useful for debugging.
    1053              : !!
    1054              : !! INPUTS
    1055              : !!  cprj(:,:) <type(pawcprj_type)>= cprj datastructure
    1056              : !!  prtgrads :: optional, 1 to print gradients also
    1057              : !!
    1058              : !! OUTPUT
    1059              : !!
    1060              : !! SOURCE
    1061              : 
    1062            0 :  subroutine pawcprj_output(cprj,prtgrads)
    1063              : 
    1064              : !Arguments ------------------------------------
    1065              : !scalar
    1066              : integer,optional :: prtgrads
    1067              : !arrays
    1068              :  type(pawcprj_type),intent(in) :: cprj(:,:)
    1069              : 
    1070              : !Local variables-------------------------------
    1071              : !scalar
    1072              :  integer :: ii,jj,kk,nlmn,n1dim,n2dim
    1073              :  logical :: gradoutput
    1074              : 
    1075              : ! *************************************************************************
    1076              : 
    1077            0 :  n1dim=size(cprj,dim=1)
    1078            0 :  n2dim=size(cprj,dim=2)
    1079            0 :  gradoutput = .FALSE.
    1080            0 :  if(present(prtgrads)) then
    1081            0 :    gradoutput = (prtgrads .EQ. 1)
    1082              :  end if
    1083              : 
    1084            0 :  write(std_out,'(a)')' pawcprj_output '
    1085              : 
    1086            0 :  do jj=1,n2dim
    1087            0 :    do ii=1,n1dim
    1088            0 :      write(std_out,'(a,i4,a,i4)')'atom ',ii,' band*k ',jj
    1089            0 :      nlmn=cprj(ii,jj)%nlmn
    1090            0 :      do kk=1,nlmn
    1091            0 :        write(std_out,'(2f12.8)')cprj(ii,jj)%cp(1,kk),cprj(ii,jj)%cp(2,kk)
    1092            0 :        if(gradoutput) then
    1093            0 :          write(std_out,'(6f12.8)')cprj(ii,jj)%dcp(1,1,kk),cprj(ii,jj)%dcp(2,1,kk),&
    1094            0 :            &cprj(ii,jj)%dcp(1,2,kk),cprj(ii,jj)%dcp(2,2,kk),&
    1095            0 :            &cprj(ii,jj)%dcp(1,3,kk),cprj(ii,jj)%dcp(2,3,kk)
    1096              :        end if
    1097              :      end do
    1098              :    end do
    1099              :  end do
    1100              : 
    1101            0 : end subroutine pawcprj_output
    1102              : !!***
    1103              : 
    1104              : !----------------------------------------------------------------------
    1105              : 
    1106              : !!****f* m_pawcprj/pawcprj_get
    1107              : !! NAME
    1108              : !! pawcprj_get
    1109              : !!
    1110              : !! FUNCTION
    1111              : !! Read the cprj_k for a given k-point from memory in cprj or from a temporary file
    1112              : !!
    1113              : !! INPUTS
    1114              : !!  atind(natom)=index table for atoms (see iorder below)
    1115              : !!  cprj(dimcp,nspinor*mband*mkmem*nsppol)=input cprj (used if mkmem/=0)
    1116              : !!  dimcp=first dimension of cprj_k,cprj arrays (1 or natom)
    1117              : !!  iband1=index of first band in cprj
    1118              : !!  ibg=shift in cprj array to locate current k-point
    1119              : !!  [icpgr]= (optional argument) if present, only component icpgr of
    1120              : !!           input cprj gradient is copied into output cprj
    1121              : !!           Not used if cprj(:,:)%ncpgr<icpgr (mkmem>0)
    1122              : !!                    or ncpgr(optional)<icpgr (mkmem=0)
    1123              : !!  ikpt=index of current k-point (only needed for the parallel distribution)
    1124              : !!  iorder=0 if cprj ordering does not change during reading
    1125              : !!         1 if cprj ordering changes during reading, depending on content of atind array:
    1126              : !!              - if atind=atindx  (type-sorted=>unsorted)
    1127              : !!              - if atind=atindx1 (unsorted=>type-sorted)
    1128              : !!  isppol=index of current spin component
    1129              : !!  mband=maximum number of bands
    1130              : !!  mkmem=number of k points which can fit in memory; set to 0 if use disk
    1131              : !!  [mpi_comm]=(optional argument) MPI communicator over (k-pts,bands,spins)
    1132              : !!             Must be used in association with proc_distrb argument
    1133              : !!  natom=number of atoms in cell
    1134              : !!  nband=number of bands to import (usually 1 or nband_k)
    1135              : !!  nband_k=total number of bands for this k-point
    1136              : !!  [ncpgr]=(optional argument) second dimension of cprj%dcp(2,ncpgr,nlmn
    1137              : !!          stored in memory (mkmem>0) or present on disk (mkmem=0))
    1138              : !!          needed only when optional argument icpgr is present
    1139              : !!  nspinor=number of spinorial components of the wavefunctions (on current proc)
    1140              : !!  nsppol=1 for unpolarized, 2 for spin-polarized
    1141              : !!  [proc_distrb(nkpt,nband,nsppol)]=(optional argument) processor distribution
    1142              : !!          Describe how cprj datastructures are distributed over processors
    1143              : !!          When present, mpicomm argument must be also present
    1144              : !!  uncp=unit number for cprj data (used if mkmem=0)
    1145              : !!
    1146              : !! OUTPUT
    1147              : !!  cprj_k(dimcp,nspinor*nband) <type(pawcprj_type)>= output cprj datastructure
    1148              : !!
    1149              : !! SOURCE
    1150              : 
    1151     10067736 :  subroutine pawcprj_get(atind,cprj_k,cprj,dimcp,iband1,ibg,ikpt,iorder,isppol,mband,&
    1152              : &                    mkmem,natom,nband,nband_k,nspinor,nsppol,uncp,&
    1153      5033868 : &                    icpgr,ncpgr,mpicomm,proc_distrb) ! optionals arguments
    1154              : 
    1155              : !Arguments ------------------------------------
    1156              : !scalars
    1157              :  integer,intent(in) :: dimcp,iband1,ibg,ikpt,iorder,isppol,mband,mkmem,natom
    1158              :  integer,intent(in) :: nband,nband_k,nspinor,nsppol,uncp
    1159              :  integer,intent(in),optional :: icpgr,mpicomm,ncpgr
    1160              : !arrays
    1161              :  integer,intent(in) :: atind(natom)
    1162              :  integer,intent(in),optional :: proc_distrb(:,:,:)
    1163              :  type(pawcprj_type),intent(in) :: cprj(dimcp,nspinor*mband*mkmem*nsppol)
    1164              :  type(pawcprj_type),intent(inout) :: cprj_k(dimcp,nspinor*nband)
    1165              : 
    1166              : !Local variables-------------------------------
    1167              : !scalars
    1168              :  integer :: iatm,iatom,ib,ibsp,icpgr_,isp,ispinor,jband,me,nband0,ncpgr_
    1169              :  logical :: has_distrb,has_icpgr
    1170              :  character(len=500) :: msg
    1171              : !arrays
    1172      5033868 :  real(dp),allocatable :: tmp(:,:,:)
    1173              : 
    1174              : ! *************************************************************************
    1175              : 
    1176      5033868 :  ncpgr_=cprj_k(1,1)%ncpgr;if (present(ncpgr)) ncpgr_=ncpgr
    1177      5033868 :  icpgr_=-1;if(present(icpgr)) icpgr_=icpgr
    1178      5033868 :  has_icpgr=(icpgr_>0.and.icpgr_<=ncpgr_)
    1179      5033868 :  if (present(icpgr).and.(.not.present(ncpgr))) then
    1180            0 :    msg='ncpgr must be present when icpgr is present (pawcprj_get)!'
    1181            0 :    LIBPAW_BUG(msg)
    1182              :  end if
    1183      5033868 :  if (has_icpgr.and.cprj_k(1,1)%ncpgr<1) then
    1184            0 :    msg='cprj_k%ncpgr not consistent with icpgr (pawcprj_get)!'
    1185            0 :    LIBPAW_BUG(msg)
    1186              :  end if
    1187              : 
    1188              : !MPI data
    1189      5033868 :  has_distrb=present(proc_distrb)
    1190      5033868 :  if (has_distrb) then
    1191      4266153 :    if (.not.present(mpicomm)) then
    1192            0 :      msg='mpicomm must be present when proc_distrb is present (pawcprj_get)!'
    1193            0 :      LIBPAW_BUG(msg)
    1194              :    end if
    1195      4266153 :    me=xmpi_comm_rank(mpicomm)
    1196              :  end if
    1197              : 
    1198      5033868 :  if (mkmem==0) then
    1199              : 
    1200            0 :    if (iband1==1) then
    1201            0 :      read(uncp) nband0
    1202            0 :      if (nband_k/=nband0) then
    1203            0 :        msg='_PAW file was not created with the right options (pawcprj_get)!'
    1204            0 :        LIBPAW_BUG(msg)
    1205              :      end if
    1206              :    end if
    1207              : 
    1208            0 :    isp=0;jband=iband1-1
    1209            0 :    do ib=1,nband
    1210            0 :      jband=jband+1
    1211            0 :      if (has_distrb) then
    1212            0 :        if (abs(proc_distrb(ikpt,jband,isppol)-me)/=0) then
    1213            0 :          isp=isp+nspinor
    1214            0 :          cycle
    1215              :        end if
    1216              :      end if
    1217            0 :      do ispinor=1,nspinor
    1218            0 :        isp=isp+1
    1219            0 :        if (iorder==0) then
    1220            0 :          if (ncpgr_==0) then
    1221            0 :            do iatom=1,dimcp
    1222            0 :              read(uncp) cprj_k(iatom,isp)%cp(:,:)
    1223              :            end do
    1224              :          else
    1225            0 :            if (has_icpgr) then
    1226            0 :              do iatom=1,dimcp
    1227            0 :                LIBPAW_ALLOCATE(tmp,(2,ncpgr_,cprj_k(iatom,1)%nlmn))
    1228            0 :                read(uncp) cprj_k(iatom,isp)%cp(:,:),tmp(:,:,:)
    1229            0 :                cprj_k(iatom,isp)%dcp(:,1,:)=tmp(:,icpgr_,:)
    1230            0 :                LIBPAW_DEALLOCATE(tmp)
    1231              :              end do
    1232              :            else
    1233            0 :              do iatom=1,dimcp
    1234            0 :                read(uncp) cprj_k(iatom,isp)%cp(:,:),cprj_k(iatom,isp)%dcp(:,:,:)
    1235              :              end do
    1236              :            end if
    1237              :          end if
    1238              :        else
    1239            0 :          if (ncpgr_==0) then
    1240            0 :            do iatom=1,dimcp
    1241            0 :              iatm=min(atind(iatom),dimcp)
    1242            0 :              read(uncp) cprj_k(iatm,isp)%cp(:,:)
    1243              :            end do
    1244              :          else
    1245            0 :            if (has_icpgr) then
    1246            0 :              do iatom=1,dimcp
    1247            0 :                iatm=min(atind(iatom),dimcp)
    1248            0 :                LIBPAW_ALLOCATE(tmp,(2,ncpgr_,cprj_k(iatm,1)%nlmn))
    1249            0 :                read(uncp) cprj_k(iatm,isp)%cp(:,:),tmp(:,:,:)
    1250            0 :                cprj_k(iatm,isp)%dcp(:,1,:)=tmp(:,icpgr_,:)
    1251            0 :                LIBPAW_DEALLOCATE(tmp)
    1252              :              end do
    1253              :            else
    1254            0 :              do iatom=1,dimcp
    1255            0 :                iatm=min(atind(iatom),dimcp)
    1256            0 :                read(uncp) cprj_k(iatm,isp)%cp(:,:),cprj_k(iatm,isp)%dcp(:,:,:)
    1257              :              end do
    1258              :            end if
    1259              :          end if
    1260              :        end if
    1261              :      end do
    1262              :    end do
    1263              : 
    1264              :  else
    1265              : 
    1266      5033868 :    isp=0;ibsp=ibg+nspinor*(iband1-1);jband=iband1-1
    1267     14620616 :    do ib=1,nband
    1268      9586748 :      jband=jband+1
    1269      9586748 :      if (has_distrb) then
    1270      8815472 :        if (abs(proc_distrb(ikpt,jband,isppol)-me)/=0) then
    1271            0 :          isp=isp+nspinor;ibsp=ibsp+nspinor
    1272            0 :          cycle
    1273              :        end if
    1274              :      end if
    1275     24444481 :      do ispinor=1,nspinor
    1276      9823865 :        isp=isp+1;ibsp=ibsp+1
    1277     19410613 :        if (iorder==0) then
    1278      7882361 :          if (ncpgr_==0) then
    1279     17889235 :            do iatom=1,dimcp
    1280    350674450 :              cprj_k(iatom,isp)%cp(:,:)=cprj(iatom,ibsp)%cp(:,:)
    1281              :            end do
    1282              :          else
    1283      2074560 :            if (has_icpgr) then
    1284      1791010 :              do iatom=1,dimcp
    1285     31196550 :                cprj_k(iatom,isp)%cp(:,:)   =cprj(iatom,ibsp)%cp(:,:)
    1286     31777522 :                cprj_k(iatom,isp)%dcp(:,1,:)=cprj(iatom,ibsp)%dcp(:,icpgr_,:)
    1287              :              end do
    1288              :            else
    1289      4509012 :              do iatom=1,dimcp
    1290     76849600 :                cprj_k(iatom,isp)%cp(:,:)   =cprj(iatom,ibsp)%cp(:,:)
    1291    602693188 :                cprj_k(iatom,isp)%dcp(:,:,:)=cprj(iatom,ibsp)%dcp(:,:,:)
    1292              :              end do
    1293              :            end if
    1294              :          end if
    1295              :        else
    1296      1941504 :          if (ncpgr_==0) then
    1297      5824512 :            do iatom=1,dimcp
    1298      3883008 :              iatm=min(atind(iatom),dimcp)
    1299     99016704 :              cprj_k(iatm,isp)%cp(:,:)=cprj(iatom,ibsp)%cp(:,:)
    1300              :            end do
    1301              :          else
    1302            0 :            if (has_icpgr) then
    1303            0 :              do iatom=1,dimcp
    1304            0 :                iatm=min(atind(iatom),dimcp)
    1305            0 :                cprj_k(iatm,isp)%cp(:,:)   =cprj(iatom,ibsp)%cp(:,:)
    1306            0 :                cprj_k(iatm,isp)%dcp(:,1,:)=cprj(iatom,ibsp)%dcp(:,icpgr_,:)
    1307              :              end do
    1308              :            else
    1309            0 :              do iatom=1,dimcp
    1310            0 :                iatm=min(atind(iatom),dimcp)
    1311            0 :                cprj_k(iatm,isp)%cp(:,:)   =cprj(iatom,ibsp)%cp(:,:)
    1312            0 :                cprj_k(iatm,isp)%dcp(:,:,:)=cprj(iatom,ibsp)%dcp(:,:,:)
    1313              :              end do
    1314              :            end if
    1315              :          end if
    1316              :        end if
    1317              :      end do
    1318              :    end do
    1319              : 
    1320              :  end if
    1321              : 
    1322      5033868 : end subroutine pawcprj_get
    1323              : !!***
    1324              : 
    1325              : !----------------------------------------------------------------------
    1326              : 
    1327              : !!****f* m_pawcprj/pawcprj_put
    1328              : !! NAME
    1329              : !! pawcprj_put
    1330              : !!
    1331              : !! FUNCTION
    1332              : !! Write cprj_k for a given set of (n,k) into memory in cprj, or into a temporary file
    1333              : !!
    1334              : !! INPUTS
    1335              : !!  atind(natom)=index table for atoms (see iorder below)
    1336              : !!  cprj_k(dimcp,nspinor*nband) <type(pawcprj_type)>= input cprj datastructure
    1337              : !!  dimcp=first dimension of cprj_k,cprjnk arrays (1 or natom)
    1338              : !!  iband1=index of first band in cprj
    1339              : !!  ibg=shift in cprj array to locate current k-point
    1340              : !!  ikpt=index of current k-point (only needed for the parallel distribution)
    1341              : !!  iorder=0 if cprj ordering does not change during reading
    1342              : !!         1 if cprj ordering changes during writing, depending on content of atind array:
    1343              : !!              - if atind=atindx  (type-sorted->unsorted)
    1344              : !!              - if atind=atindx1 (unsorted->type-sorted)
    1345              : !!  isppol=index of current spin component
    1346              : !!  mband=maximum number of bands
    1347              : !!  mkmem=number of k points which can fit in memory; set to 0 if use disk
    1348              : !!  [mpi_comm]=(optional argument) MPI communicator over (k-pts,bands,spins)
    1349              : !!             Must be used in association with proc_distrb argument
    1350              : !!  [mpi_comm_band]=(optional argument) MPI communicator over bands
    1351              : !!             Must be used in association with proc_distrb argument
    1352              : !!  natom=number of atoms in cell
    1353              : !!  nband=number of bands to export (usually 1, nband_k or nblockbd)
    1354              : !!  nband_k=total number of bands for this k-point
    1355              : !!  nlmn(dimcp)=array of dimensions of cprj_k,cprjnk datastructures
    1356              : !!  nspinor=number of spinorial components of the wavefunctions (on current proc)
    1357              : !!  nsppol=1 for unpolarized, 2 for spin-polarized
    1358              : !!  [proc_distrb(nkpt,nband,nsppol)]=(optional argument) processor distribution
    1359              : !!          Describe how cprj datastructures are distributed over processors
    1360              : !!          When present, mpicomm argument must be also present
    1361              : !!  [to_be_gathered]=(optional argument) TRUE if cprj_k arrays have to be
    1362              : !!                   gathered between procs (band-fft parallelism only)
    1363              : !!  uncp=unit number for cprj data (used if mkmem=0)
    1364              : !!
    1365              : !! SIDE EFFECTS
    1366              : !!  cprj(dimcp,nspinor*mband*mkmem*nsppol)=output cprj (used if mkmem/=0)
    1367              : !!
    1368              : !! SOURCE
    1369              : 
    1370      5025816 :  subroutine pawcprj_put(atind,cprj_k,cprj,dimcp,iband1,ibg,ikpt,iorder,isppol,mband,&
    1371      2512908 : &           mkmem,natom,nband,nband_k,nlmn,nspinor,nsppol,uncp,&
    1372      2512908 : &           mpicomm,mpi_comm_band,proc_distrb,to_be_gathered) ! Optional arguments
    1373              : 
    1374              : !Arguments ------------------------------------
    1375              : !scalars
    1376              :  integer,intent(in) :: iband1,ibg,ikpt,iorder,isppol,dimcp,mband,mkmem
    1377              :  integer,intent(in) :: natom,nband,nband_k,nspinor,nsppol,uncp
    1378              :  integer,intent(in),optional :: mpicomm,mpi_comm_band
    1379              :  logical,optional,intent(in) :: to_be_gathered
    1380              : !arrays
    1381              :  integer,intent(in) :: atind(natom),nlmn(dimcp)
    1382              :  integer,intent(in),optional :: proc_distrb(:,:,:)
    1383              :  type(pawcprj_type),intent(inout) :: cprj(dimcp,nspinor*mband*mkmem*nsppol)
    1384              :  type(pawcprj_type),intent(in) :: cprj_k(dimcp,nspinor*nband)
    1385              : 
    1386              : !Local variables-------------------------------
    1387              : !scalars
    1388              :  integer :: iatm,iatom,iband,ibsp,icpgr,ierr,ii,ilmn,isp,ispinor,jband,jj
    1389              :  integer :: lmndim,me,ncpgr,nproc_band
    1390              :  logical :: has_distrb,to_be_gathered_
    1391              :  character(len=500) :: msg
    1392              : !arrays
    1393      2512908 :  real(dp),allocatable :: buffer1(:),buffer2(:)
    1394              : ! *************************************************************************
    1395              : 
    1396      2512908 :  ncpgr=cprj_k(1,1)%ncpgr
    1397      2512908 :  to_be_gathered_=.false.;if (present(to_be_gathered)) to_be_gathered_=to_be_gathered
    1398              : 
    1399              : !MPI data
    1400      2512908 :  nproc_band=1;if (present(mpi_comm_band)) nproc_band=xmpi_comm_size(mpi_comm_band)
    1401      2512908 :  has_distrb=present(proc_distrb)
    1402      2512908 :  if (has_distrb) then
    1403      1577190 :    if (.not.present(mpicomm)) then
    1404            0 :      msg='mpicomm must be present when proc_distrb is present (pawcprj_put)!'
    1405            0 :      LIBPAW_BUG(msg)
    1406              :    end if
    1407      1577190 :    me=xmpi_comm_rank(mpicomm)
    1408              :  end if
    1409              : 
    1410      2512908 :  if (nproc_band==1.or.(.not.to_be_gathered_)) then
    1411              : 
    1412      2512908 :    if (mkmem==0) then
    1413              : 
    1414            0 :      if (iband1==1) write(uncp) nband_k
    1415              : 
    1416            0 :      isp=0;jband=iband1-1
    1417            0 :      do iband=1,nband
    1418            0 :        jband=jband+1
    1419            0 :        if (has_distrb) then
    1420            0 :          if (abs(proc_distrb(ikpt,jband,isppol)-me)/=0) then
    1421            0 :            isp=isp+nspinor
    1422            0 :            cycle
    1423              :          end if
    1424              :        end if
    1425            0 :        do ispinor=1,nspinor
    1426            0 :          isp=isp+1
    1427            0 :          if (iorder==0) then
    1428            0 :            do iatom=1,dimcp
    1429            0 :              if (ncpgr==0) then
    1430            0 :                write(uncp) cprj_k(iatom,isp)%cp(:,:)
    1431              :              else
    1432            0 :                write(uncp) cprj_k(iatom,isp)%cp(:,:),cprj_k(iatom,isp)%dcp(:,:,:)
    1433              :              end if
    1434              :            end do
    1435              :          else
    1436            0 :            do iatom=1,dimcp
    1437            0 :              iatm=min(atind(iatom),dimcp)
    1438            0 :              if (ncpgr==0) then
    1439            0 :                write(uncp) cprj_k(iatm,isp)%cp(:,:)
    1440              :              else
    1441            0 :                write(uncp) cprj_k(iatm,isp)%cp(:,:),cprj_k(iatm,isp)%dcp(:,:,:)
    1442              :              end if
    1443              :            end do
    1444              :          end if
    1445              :        end do
    1446              :      end do
    1447              : 
    1448              :    else
    1449              : 
    1450      2512908 :      isp=0;ibsp=ibg+nspinor*(iband1-1);jband=iband1-1
    1451      5195050 :      do iband=1,nband
    1452      2682142 :        jband=jband+1
    1453      2682142 :        if (has_distrb) then
    1454      1726456 :          if (abs(proc_distrb(ikpt,jband,isppol)-me)/=0) then
    1455            0 :            isp=isp+nspinor;ibsp=ibsp+nspinor
    1456            0 :            cycle
    1457              :          end if
    1458              :        end if
    1459      8094708 :        do ispinor=1,nspinor
    1460      2899658 :          isp=isp+1;ibsp=ibsp+1
    1461      5581800 :          if (iorder==0) then
    1462      7569322 :            do iatom=1,dimcp
    1463    176813092 :              cprj(iatom,ibsp)%cp(:,:)=cprj_k(iatom,isp)%cp(:,:)
    1464     40447326 :              if (ncpgr>0) cprj(iatom,ibsp)%dcp(:,:,:)=cprj_k(iatom,isp)%dcp(:,:,:)
    1465              :            end do
    1466              :          else
    1467      1869378 :            do iatom=1,dimcp
    1468      1243954 :              iatm=min(atind(iatom),dimcp)
    1469     31256260 :              cprj(iatom,ibsp)%cp(:,:)=cprj_k(iatm,isp)%cp(:,:)
    1470     10860610 :              if (ncpgr>0) cprj(iatom,ibsp)%dcp(:,:,:)=cprj_k(iatm,isp)%dcp(:,:,:)
    1471              :            end do
    1472              :          end if
    1473              :        end do
    1474              :      end do
    1475              : 
    1476              :    end if
    1477              : 
    1478              :  else ! np_band>1
    1479              : 
    1480            0 :    lmndim=2*sum(nlmn(1:dimcp))*(1+ncpgr)*nspinor
    1481            0 :    LIBPAW_ALLOCATE(buffer1,(lmndim))
    1482            0 :    LIBPAW_ALLOCATE(buffer2,(lmndim*nproc_band))
    1483            0 :    isp=0;ibsp=ibg+nspinor*(iband1-1)
    1484            0 :    do iband=1,nband  ! must be nblockbd for band-fft parallelism
    1485              :      jj=1
    1486            0 :      do ispinor=1,nspinor
    1487            0 :        isp=isp+1
    1488            0 :        do iatom=1,dimcp
    1489            0 :          if (iorder==0) then
    1490              :            iatm=iatom
    1491              :          else
    1492            0 :            iatm=min(atind(iatom),dimcp)
    1493              :          end if
    1494            0 :          do ilmn=1,nlmn(iatm)
    1495            0 :            buffer1(jj:jj+1)=cprj_k(iatm,isp)%cp(1:2,ilmn)
    1496            0 :            jj=jj+2
    1497              :          end do
    1498            0 :          if (ncpgr>0) then
    1499            0 :            do ilmn=1,nlmn(iatm)
    1500            0 :              do icpgr=1,ncpgr
    1501            0 :                buffer1(jj:jj+1)=cprj_k(iatm,isp)%dcp(1:2,icpgr,ilmn)
    1502            0 :                jj=jj+2
    1503              :              end do
    1504              :            end do
    1505              :          end if
    1506              :        end do !iatom
    1507              :      end do !ispinor
    1508            0 :      call xmpi_allgather(buffer1,lmndim,buffer2,mpi_comm_band,ierr)
    1509            0 :      jj=1
    1510            0 :      do ii=1,nproc_band
    1511            0 :        do ispinor=1,nspinor
    1512            0 :          ibsp=ibsp+1
    1513            0 :          do iatom=1,dimcp
    1514            0 :            if (iorder==0) then
    1515              :              iatm=iatom
    1516              :            else
    1517            0 :              iatm=min(atind(iatom),dimcp)
    1518              :            end if
    1519            0 :            do ilmn=1,nlmn(iatm)
    1520            0 :              cprj(iatom,ibsp)%cp(1:2,ilmn)=buffer2(jj:jj+1)
    1521            0 :              jj=jj+2
    1522              :            end do
    1523            0 :            if (ncpgr>0) then
    1524            0 :              do ilmn=1,nlmn(iatm)
    1525            0 :                do icpgr=1,ncpgr
    1526            0 :                  cprj(iatom,ibsp)%dcp(1:2,icpgr,ilmn)=buffer2(jj:jj+1)
    1527            0 :                  jj=jj+2
    1528              :                end do
    1529              :              end do
    1530              :            end if
    1531              :          end do !iatom
    1532              :        end do !ispinor
    1533              :      end do !ii=1,nproc_band
    1534              :    end do !iband
    1535            0 :    LIBPAW_DEALLOCATE(buffer1)
    1536            0 :    LIBPAW_DEALLOCATE(buffer2)
    1537              : 
    1538              :  end if ! mode_para=b, nband
    1539              : 
    1540      2512908 : end subroutine pawcprj_put
    1541              : !!***
    1542              : 
    1543              : !----------------------------------------------------------------------
    1544              : 
    1545              : !!****f* m_pawcprj/pawcprj_reorder
    1546              : !! NAME
    1547              : !! pawcprj_reorder
    1548              : !!
    1549              : !! FUNCTION
    1550              : !! Change the order of a cprj datastructure
    1551              : !!   From unsorted cprj to atom-sorted cprj (atm_indx=atindx)
    1552              : !!   From atom-sorted cprj to unsorted cprj (atm_indx=atindx1)
    1553              : !!
    1554              : !! INPUTS
    1555              : !!  atm_indx(natom)=index table for atoms
    1556              : !!   From unsorted cprj to atom-sorted cprj (atm_indx=atindx)
    1557              : !!   From atom-sorted cprj to unsorted cprj (atm_indx=atindx1)
    1558              : !!
    1559              : !! OUTPUT
    1560              : !!
    1561              : !! SIDE EFFECTS
    1562              : !!  cprj(:,:) <type(pawcprj_type)>= cprj datastructure
    1563              : !!
    1564              : !! SOURCE
    1565              : 
    1566          844 :  subroutine pawcprj_reorder(cprj,atm_indx)
    1567              : 
    1568              : !Arguments ------------------------------------
    1569              : !scalars
    1570              : !arrays
    1571              :  integer,intent(in) :: atm_indx(:)
    1572              :  type(pawcprj_type),intent(inout) :: cprj(:,:)
    1573              : 
    1574              : !Local variables-------------------------------
    1575              : !scalars
    1576              :  integer :: iexit,ii,jj,kk,n1atindx,n1cprj,n2cprj,ncpgr
    1577              :  character(len=100) :: msg
    1578              : !arrays
    1579          844 :  integer,allocatable :: nlmn(:)
    1580          844 :  type(pawcprj_type),allocatable :: cprj_tmp(:,:)
    1581              : 
    1582              : ! *************************************************************************
    1583              : 
    1584          844 :  n1cprj=size(cprj,dim=1);n2cprj=size(cprj,dim=2)
    1585          844 :  n1atindx=size(atm_indx,dim=1)
    1586          844 :  if (n1cprj==0.or.n2cprj==0.or.n1atindx<=1) return
    1587              : 
    1588          819 :  if (n1cprj/=n1atindx) then
    1589            0 :    msg='wrong sizes (pawcprj_reorder)!'
    1590            0 :    LIBPAW_BUG(msg)
    1591              :  end if
    1592              : 
    1593              : !Nothing to do when the atoms are already sorted
    1594          819 :  iexit=1;ii=0
    1595         2501 :  do while (iexit==1.and.ii<n1atindx)
    1596         1682 :    ii=ii+1
    1597         2501 :    if (atm_indx(ii)/=ii) iexit=0
    1598              :  end do
    1599          819 :  if (iexit==1) return
    1600              : 
    1601            0 :  LIBPAW_ALLOCATE(nlmn,(n1cprj))
    1602            0 :  do ii=1,n1cprj
    1603            0 :    nlmn(ii)=cprj(ii,1)%nlmn
    1604              :  end do
    1605            0 :  ncpgr=cprj(1,1)%ncpgr
    1606            0 :  LIBPAW_DATATYPE_ALLOCATE(cprj_tmp,(n1cprj,n2cprj))
    1607            0 :  call pawcprj_alloc(cprj_tmp,ncpgr,nlmn)
    1608            0 :  call pawcprj_copy(cprj,cprj_tmp)
    1609            0 :  call pawcprj_free(cprj)
    1610              : 
    1611            0 :  do jj=1,n2cprj
    1612            0 :    do ii=1,n1cprj
    1613            0 :      kk=atm_indx(ii)
    1614            0 :      cprj(kk,jj)%nlmn=nlmn(ii)
    1615            0 :      cprj(kk,jj)%ncpgr=ncpgr
    1616            0 :      LIBPAW_ALLOCATE(cprj(kk,jj)%cp,(2,nlmn(ii)))
    1617            0 :      cprj(kk,jj)%cp(:,:)=cprj_tmp(ii,jj)%cp(:,:)
    1618            0 :      if (ncpgr>0) then
    1619            0 :        LIBPAW_ALLOCATE(cprj(kk,jj)%dcp,(2,ncpgr,nlmn(ii)))
    1620            0 :        cprj(kk,jj)%dcp(:,:,:)=cprj_tmp(ii,jj)%dcp(:,:,:)
    1621              :      end if
    1622              :    end do
    1623              :  end do
    1624              : 
    1625            0 :  call pawcprj_free(cprj_tmp)
    1626            0 :  LIBPAW_DATATYPE_DEALLOCATE(cprj_tmp)
    1627            0 :  LIBPAW_DEALLOCATE(nlmn)
    1628              : 
    1629          844 : end subroutine pawcprj_reorder
    1630              : !!***
    1631              : 
    1632              : !----------------------------------------------------------------------
    1633              : 
    1634              : !!****f* m_pawcprj/pawcprj_mpi_exch
    1635              : !! NAME
    1636              : !! pawcprj_mpi_exch
    1637              : !!
    1638              : !! FUNCTION
    1639              : !! Exchange a pawcprj_type between two processors inside a MPI communicator.
    1640              : !!
    1641              : !! INPUTS
    1642              : !!  natom=Number of atoms (size of first dimension of Cprj_send and Cprj_recv).
    1643              : !!  n2dim=Size of the second dimension.
    1644              : !!  nlmn(natom)=Number of nlm partial waves for each atom.
    1645              : !!  Cprj_send= The datatype to be transmitted.
    1646              : !!  receiver=ID of the receiver in spaceComm.
    1647              : !!  sender=ID of the sender in spaceComm.
    1648              : !!  spaceComm=MPI Communicator.
    1649              : !!  mtag= message tag
    1650              : !!
    1651              : !! OUTPUT
    1652              : !!  ierr=Error status.
    1653              : !!  Cprj_recv=The datatype copied on proc. receiver.
    1654              : !!
    1655              : !! NOTES
    1656              : !!  If sender==receiver, Cprj_send is copied into Cprj_recv.
    1657              : !!  It should be easy to avoid this additional copy in the calling routine.
    1658              : !!
    1659              : !! SOURCE
    1660              : 
    1661           16 : subroutine pawcprj_mpi_exch(natom,n2dim,nlmn,ncpgr,Cprj_send,Cprj_recv,sender,receiver,spaceComm,mtag,ierr)
    1662              : 
    1663              : !Arguments ------------------------------------
    1664              : !scalars
    1665              :  integer,intent(in) :: mtag,natom,n2dim,ncpgr
    1666              :  integer,intent(in) :: sender,receiver,spaceComm
    1667              :  integer,intent(out) :: ierr
    1668              : !arrays
    1669              :  integer,intent(in) :: nlmn(natom)
    1670              :  type(pawcprj_type),intent(in) :: Cprj_send(:,:)
    1671              :  type(pawcprj_type),intent(inout) :: Cprj_recv(:,:)
    1672              : 
    1673              : !Local variables-------------------------------
    1674              : !scalars
    1675              :  integer :: iat,jj,t2dim,tcpgr,n1dim,nn
    1676              :  integer :: ntotcp,ipck,rank
    1677              :  character(len=500) :: msg
    1678              : !arrays
    1679           16 :  real(dp),allocatable :: buffer_cp(:,:),buffer_cpgr(:,:,:)
    1680              : 
    1681              : ! *************************************************************************
    1682              : 
    1683           16 :  n1dim=0
    1684           16 :  t2dim=0
    1685           16 :  tcpgr=0
    1686           16 :  ierr=0
    1687           16 :  if (sender==receiver) then
    1688            0 :    call pawcprj_copy(Cprj_send,Cprj_recv)
    1689              :    return
    1690              :  end if
    1691              : 
    1692           16 :  rank = xmpi_comm_rank(spaceComm)
    1693              : 
    1694           16 :  nn=size(nlmn,dim=1)
    1695           16 :  if (rank==sender) then
    1696            8 :    n1dim=size(Cprj_send,dim=1)
    1697            8 :    t2dim=size(Cprj_send,dim=2)
    1698            8 :    tcpgr=Cprj_send(1,1)%ncpgr
    1699              :  end if
    1700           16 :  if (rank==receiver) then
    1701            8 :    n1dim=size(Cprj_recv,dim=1)
    1702            8 :    t2dim=size(Cprj_recv,dim=2)
    1703            8 :    tcpgr=Cprj_recv(1,1)%ncpgr
    1704              :  end if
    1705           16 :  if (rank/=sender.and.rank/=receiver) then
    1706              :    write(msg,'(a,3i0)') &
    1707            0 : &   'rank is not equal to sender or receiver (pawcprj_mpi_exch): ',rank, sender, receiver
    1708            0 :    LIBPAW_BUG(msg)
    1709              :  end if
    1710              : 
    1711           32 :  ntotcp=n2dim*SUM(nlmn(:))
    1712              : 
    1713           48 :  LIBPAW_ALLOCATE(buffer_cp,(2,ntotcp))
    1714           16 :  if (ncpgr/=0)  then
    1715            0 :    LIBPAW_ALLOCATE(buffer_cpgr,(2,ncpgr,ntotcp))
    1716              :  end if
    1717              : 
    1718              : !=== Pack Cprj_send ===
    1719           16 :  if (rank==sender) then
    1720              :    ipck=0
    1721          168 :    do jj=1,n2dim
    1722          328 :      do iat=1,natom
    1723          160 :        nn=nlmn(iat)
    1724         8800 :        buffer_cp(:,ipck+1:ipck+nn)=Cprj_send(iat,jj)%cp(:,1:nn)
    1725          160 :        if (ncpgr/=0) buffer_cpgr(:,:,ipck+1:ipck+nn)=Cprj_send(iat,jj)%dcp(:,:,1:nn)
    1726          320 :        ipck=ipck+nn
    1727              :      end do
    1728              :    end do
    1729              :  end if
    1730              : 
    1731              : !=== Transmit data ===
    1732           16 :  call xmpi_exch(buffer_cp,2*ntotcp,sender,buffer_cp,receiver,spaceComm,2*mtag,ierr)
    1733           16 :  if (ncpgr/=0) then
    1734            0 :    call xmpi_exch(buffer_cpgr,2*ncpgr*ntotcp,sender,buffer_cpgr,receiver,spaceComm,2*mtag+1,ierr)
    1735              :  end if
    1736              : 
    1737              : !=== UnPack buffers into Cprj_recv ===
    1738           16 :  if (rank==receiver) then
    1739              :    ipck=0
    1740          168 :    do jj=1,n2dim
    1741          328 :      do iat=1,natom
    1742          160 :        nn=nlmn(iat)
    1743         8800 :        Cprj_recv(iat,jj)%cp(:,1:nn)=buffer_cp(:,ipck+1:ipck+nn)
    1744          160 :        if (ncpgr/=0) Cprj_recv(iat,jj)%dcp(:,:,1:nn)=buffer_cpgr(:,:,ipck+1:ipck+nn)
    1745          320 :        ipck=ipck+nn
    1746              :      end do
    1747              :    end do
    1748              :  end if
    1749              : 
    1750           16 :  LIBPAW_DEALLOCATE(buffer_cp)
    1751           16 :  if (ncpgr/=0)  then
    1752            0 :    LIBPAW_DEALLOCATE(buffer_cpgr)
    1753              :  end if
    1754              : 
    1755           16 : end subroutine pawcprj_mpi_exch
    1756              : !!***
    1757              : 
    1758              : !----------------------------------------------------------------------
    1759              : 
    1760              : !!****f* m_pawcprj/pawcprj_mpi_send
    1761              : !! NAME
    1762              : !! pawcprj_mpi_send
    1763              : !!
    1764              : !! FUNCTION
    1765              : !! Send a pawcprj_type inside a MPI communicator.
    1766              : !!
    1767              : !! INPUTS
    1768              : !!  natom=Number of atoms (size of first dimension of cprj_out).
    1769              : !!  n2dim=Size of the second dimension.
    1770              : !!  nlmn(natom)=Number of nlm partial waves for each atom.
    1771              : !!  ncpgr = number of gradients in cprj_out
    1772              : !!  cprj_out= The datatype to be transmitted.
    1773              : !!  receiver=ID of the receiver in spaceComm.
    1774              : !!  spaceComm=MPI Communicator.
    1775              : !!
    1776              : !! OUTPUT
    1777              : !!  ierr=Error status.
    1778              : !!
    1779              : !! NOTES
    1780              : !!   perhaps in general it is more efficient to use pawcprj_mpi_exch but it is
    1781              : !!   convenient for coding to have separate send and recieve routines.
    1782              : !!
    1783              : !! SOURCE
    1784              : 
    1785        40884 : subroutine pawcprj_mpi_send(natom,n2dim,nlmn,ncpgr,cprj_out,receiver,spaceComm,ierr)
    1786              : 
    1787              : !Arguments ------------------------------------
    1788              : !scalars
    1789              :  integer,intent(in) :: natom,n2dim,ncpgr
    1790              :  integer,intent(in) :: receiver,spaceComm
    1791              :  integer,intent(out) :: ierr
    1792              : !arrays
    1793              :  integer,intent(in) :: nlmn(natom)
    1794              :  type(pawcprj_type),intent(in) :: cprj_out(:,:)
    1795              : 
    1796              : !Local variables-------------------------------
    1797              : !scalars
    1798              :  integer :: iat,jj,t2dim,tcpgr,n1dim,nn
    1799              :  integer :: ntotcp,ipck,tag
    1800              :  character(len=100) :: msg
    1801              : !arrays
    1802        40884 :  real(dp),allocatable :: buffer_cp(:,:),buffer_cpgr(:,:,:)
    1803              : 
    1804              : ! *************************************************************************
    1805              : 
    1806        40884 :  n1dim=0
    1807        40884 :  t2dim=0
    1808        40884 :  tcpgr=0
    1809        40884 :  ierr=0
    1810              : 
    1811        40884 :  nn=size(nlmn,dim=1)
    1812        40884 :  n1dim=size(cprj_out,dim=1)
    1813        40884 :  t2dim=size(cprj_out,dim=2)
    1814        40884 :  tcpgr=cprj_out(1,1)%ncpgr
    1815              : 
    1816        40884 :  if (nn/=n1dim) then
    1817            0 :    msg='size mismatch in natom (pawcprj_mpi_send)!'
    1818            0 :    LIBPAW_BUG(msg)
    1819              :  end if
    1820        40884 :  if (t2dim/=n2dim) then
    1821            0 :    msg='size mismatch in dim=2 (pawcprj_mpi_send)!'
    1822            0 :    LIBPAW_BUG(msg)
    1823              :  end if
    1824        40884 :  if (tcpgr/=ncpgr) then
    1825            0 :    msg='size mismatch in ncpgr (pawcprj_mpi_send)!'
    1826            0 :    LIBPAW_BUG(msg)
    1827              :  end if
    1828              : 
    1829       122652 :  ntotcp=n2dim*SUM(nlmn(:))
    1830              : 
    1831       122652 :  LIBPAW_ALLOCATE(buffer_cp,(2,ntotcp))
    1832        40884 :  if (ncpgr/=0)  then
    1833       142080 :    LIBPAW_ALLOCATE(buffer_cpgr,(2,ncpgr,ntotcp))
    1834              :  end if
    1835              : 
    1836              : !=== Pack cprj_out ====
    1837              :  ipck=0
    1838       204420 :  do jj=1,n2dim
    1839       531492 :    do iat=1,natom
    1840       327072 :      nn=nlmn(iat)
    1841      8176800 :      buffer_cp(:,ipck+1:ipck+nn)=cprj_out(iat,jj)%cp(:,1:nn)
    1842     63978912 :      if (ncpgr/=0) buffer_cpgr(:,:,ipck+1:ipck+nn)=cprj_out(iat,jj)%dcp(:,:,1:nn)
    1843       490608 :      ipck=ipck+nn
    1844              :    end do
    1845              :  end do
    1846              : 
    1847              : !=== Transmit data ===
    1848        40884 :  tag = 2*ntotcp
    1849        40884 :  call xmpi_send(buffer_cp,receiver,tag,spaceComm,ierr)
    1850        40884 :  if (ncpgr/=0) then
    1851        35520 :    tag=tag*ncpgr
    1852        35520 :    call xmpi_send(buffer_cpgr,receiver,tag,spaceComm,ierr)
    1853              :  end if
    1854              : 
    1855              : !=== Clean up ===
    1856        40884 :  LIBPAW_DEALLOCATE(buffer_cp)
    1857        40884 :  if (ncpgr/=0)  then
    1858        35520 :    LIBPAW_DEALLOCATE(buffer_cpgr)
    1859              :  end if
    1860              : 
    1861        81768 : end subroutine pawcprj_mpi_send
    1862              : !!***
    1863              : 
    1864              : !----------------------------------------------------------------------
    1865              : 
    1866              : !!****f* m_pawcprj/pawcprj_mpi_recv
    1867              : !! NAME
    1868              : !! pawcprj_mpi_recv
    1869              : !!
    1870              : !! FUNCTION
    1871              : !! Receive a pawcprj_type inside a MPI communicator.
    1872              : !!
    1873              : !! INPUTS
    1874              : !!  natom=Number of atoms (size of first dimension of Cprj_in).
    1875              : !!  n2dim=Size of the second dimension.
    1876              : !!  nlmn(natom)=Number of nlm partial waves for each atom.
    1877              : !!  ncpgr = number of gradients in cprj_in
    1878              : !!  sender=ID of the sender in spaceComm.
    1879              : !!  spaceComm=MPI Communicator.
    1880              : !!
    1881              : !! OUTPUT
    1882              : !!  ierr=Error status.
    1883              : !!  cprj_in=The datatype copied on proc. receiver.
    1884              : !!
    1885              : !! NOTES
    1886              : !!   Perhaps in general it is more efficient to use pawcprj_mpi_exch but it is
    1887              : !!   convenient for coding to have separate send and receive routines.
    1888              : !!
    1889              : !! SOURCE
    1890              : 
    1891        40884 : subroutine pawcprj_mpi_recv(natom,n2dim,nlmn,ncpgr,cprj_in,sender,spaceComm,ierr)
    1892              : 
    1893              : !Arguments ------------------------------------
    1894              : !scalars
    1895              :  integer,intent(in) :: natom,n2dim,ncpgr
    1896              :  integer,intent(in) :: sender,spaceComm
    1897              :  integer,intent(out) :: ierr
    1898              : !arrays
    1899              :  integer,intent(in) :: nlmn(natom)
    1900              :  type(pawcprj_type),intent(inout) :: cprj_in(:,:)
    1901              : 
    1902              : !Local variables-------------------------------
    1903              : !scalars
    1904              :  integer :: iat,jj,t2dim,tcpgr,n1dim,nn
    1905              :  integer :: ntotcp,ipck,tag
    1906              :  character(len=100) :: msg
    1907              : !arrays
    1908        40884 :  real(dp),allocatable :: buffer_cp(:,:),buffer_cpgr(:,:,:)
    1909              : 
    1910              : ! *************************************************************************
    1911              : 
    1912        40884 :  n1dim=0
    1913        40884 :  t2dim=0
    1914        40884 :  tcpgr=0
    1915        40884 :  ierr=0
    1916              : 
    1917        40884 :  nn=size(nlmn,dim=1)
    1918        40884 :  n1dim=size(cprj_in,dim=1)
    1919        40884 :  t2dim=size(cprj_in,dim=2)
    1920        40884 :  tcpgr=cprj_in(1,1)%ncpgr
    1921              : 
    1922        40884 :  if (nn/=n1dim) then
    1923            0 :    msg='size mismatch in natom (pawcprj_mpi_recv)!'
    1924            0 :    LIBPAW_BUG(msg)
    1925              :  end if
    1926        40884 :  if (t2dim/=n2dim) then
    1927            0 :    msg='size mismatch in dim=2 (pawcprj_mpi_recv)!'
    1928            0 :    LIBPAW_BUG(msg)
    1929              :  end if
    1930        40884 :  if (tcpgr/=ncpgr) then
    1931            0 :    msg='size mismatch in ncpgr (pawcprj_mpi_recv)!'
    1932            0 :    LIBPAW_BUG(msg)
    1933              :  end if
    1934              : 
    1935       122652 :  ntotcp=n2dim*SUM(nlmn(:))
    1936              : 
    1937       122652 :  LIBPAW_ALLOCATE(buffer_cp,(2,ntotcp))
    1938        40884 :  if (ncpgr/=0)  then
    1939       142080 :    LIBPAW_ALLOCATE(buffer_cpgr,(2,ncpgr,ntotcp))
    1940              :  end if
    1941              : 
    1942              : !=== Receive data ===
    1943        40884 :  tag = 2*ntotcp
    1944        40884 :  call xmpi_recv(buffer_cp,sender,tag,spaceComm,ierr)
    1945        40884 :  if (ncpgr/=0) then
    1946        35520 :    tag=tag*ncpgr
    1947        35520 :    call xmpi_recv(buffer_cpgr,sender,tag,spaceComm,ierr)
    1948              :  end if
    1949              : 
    1950              : !=== UnPack buffers into cprj_in ===
    1951              :  ipck=0
    1952       204420 :  do jj=1,n2dim
    1953       531492 :    do iat=1,natom
    1954       327072 :      nn=nlmn(iat)
    1955      8176800 :      cprj_in(iat,jj)%cp(:,1:nn)=buffer_cp(:,ipck+1:ipck+nn)
    1956     63978912 :      if (ncpgr/=0) cprj_in(iat,jj)%dcp(:,:,1:nn)=buffer_cpgr(:,:,ipck+1:ipck+nn)
    1957       490608 :      ipck=ipck+nn
    1958              :    end do
    1959              :  end do
    1960              : 
    1961              : !=== Clean up ===
    1962        40884 :  LIBPAW_DEALLOCATE(buffer_cp)
    1963        40884 :  if (ncpgr/=0)  then
    1964        35520 :    LIBPAW_DEALLOCATE(buffer_cpgr)
    1965              :  end if
    1966              : 
    1967        81768 : end subroutine pawcprj_mpi_recv
    1968              : !!***
    1969              : 
    1970              : !----------------------------------------------------------------------
    1971              : 
    1972              : !!****f* m_pawcprj/pawcprj_mpi_sum
    1973              : !! NAME
    1974              : !! pawcprj_mpi_sum
    1975              : !!
    1976              : !! FUNCTION
    1977              : !! Perform MPI_SUM on a pawcprj_type inside a MPI communicator.
    1978              : !!
    1979              : !! INPUTS
    1980              : !!  spaceComm=MPI Communicator.
    1981              : !!
    1982              : !! SIDE EFFECTS
    1983              : !!  cprj=the cprj datastructure
    1984              : !!  ierr=Error status.
    1985              : !!
    1986              : !! SOURCE
    1987              : 
    1988       880950 : subroutine pawcprj_mpi_sum(cprj,spaceComm,ierr)
    1989              : 
    1990              : !Arguments ------------------------------------
    1991              : !scalars
    1992              :  integer,intent(in) :: spaceComm
    1993              :  integer,intent(out) :: ierr
    1994              : !arrays
    1995              :  type(pawcprj_type),intent(inout) :: cprj(:,:)
    1996              : 
    1997              : !Local variables-------------------------------
    1998              : !scalars
    1999              :  integer,parameter :: maxBytes=100*1024*1024 ! 100 MBytes
    2000              :  integer :: ii,ipck,jj,ncpgr,nlmn,nn,n1dim,n2dim,n2dim1,n2dim2,sizeBytes,step
    2001              :  logical,parameter :: save_memory=.true.
    2002              : !arrays
    2003       880950 :  real(dp),allocatable :: buffer_cprj(:,:,:)
    2004              : 
    2005              : ! *************************************************************************
    2006              : 
    2007       880950 :  if (xmpi_comm_size(spaceComm)<2) return
    2008              : 
    2009        12000 :  n1dim=size(cprj,1);n2dim=size(cprj,2)
    2010       132000 :  nlmn=sum(cprj(:,:)%nlmn)
    2011       132000 :  ncpgr=maxval(cprj(:,:)%ncpgr)
    2012              : 
    2013        12000 :  step=n2dim
    2014              :  if (save_memory) then
    2015        12000 :    sizeBytes=2*(1+ncpgr)*nlmn *8
    2016        12000 :    step=n2dim/max(1,sizeBytes/maxBytes)
    2017        12000 :    if (step==0) step=1
    2018              :  end if
    2019              : 
    2020        12000 :  do n2dim1=1,n2dim,step
    2021              : 
    2022        12000 :    n2dim2=min(n2dim1+step-1,n2dim)
    2023       132000 :    nlmn=sum(cprj(:,n2dim1:n2dim2)%nlmn)
    2024        60000 :    LIBPAW_ALLOCATE(buffer_cprj,(2,1+ncpgr,nlmn))
    2025              : 
    2026      3468000 :    ipck=0 ; buffer_cprj=zero
    2027        24000 :    do jj=n2dim1,n2dim2
    2028       132000 :      do ii=1,n1dim
    2029       108000 :        nn=cprj(ii,jj)%nlmn
    2030      2700000 :        buffer_cprj(:,1,ipck+1:ipck+nn)=cprj(ii,jj)%cp(:,1:nn)
    2031       108000 :        if (cprj(ii,jj)%ncpgr/=0) buffer_cprj(:,2:1+ncpgr,ipck+1:ipck+nn)=cprj(ii,jj)%dcp(:,1:ncpgr,1:nn)
    2032       120000 :        ipck=ipck+nn
    2033              :      end do
    2034              :    end do
    2035              : 
    2036        12000 :    call xmpi_sum(buffer_cprj,spaceComm,ierr)
    2037              : 
    2038        12000 :    ipck=0
    2039        24000 :    do jj=n2dim1,n2dim2
    2040       132000 :      do ii=1,n1dim
    2041       108000 :        nn=cprj(ii,jj)%nlmn
    2042      2700000 :        cprj(ii,jj)%cp(:,1:nn)=buffer_cprj(:,1,ipck+1:ipck+nn)
    2043       108000 :        if (cprj(ii,jj)%ncpgr/=0) cprj(ii,jj)%dcp(:,1:ncpgr,1:nn)=buffer_cprj(:,2:1+ncpgr,ipck+1:ipck+nn)
    2044       120000 :        ipck=ipck+nn
    2045              :      end do
    2046              :    end do
    2047              : 
    2048        24000 :    LIBPAW_DEALLOCATE(buffer_cprj)
    2049              : 
    2050              :  end do
    2051              : 
    2052        12000 : end subroutine pawcprj_mpi_sum
    2053              : !!***
    2054              : 
    2055              : !----------------------------------------------------------------------
    2056              : 
    2057              : !!****f* m_pawcprj/pawcprj_mpi_allgather
    2058              : !! NAME
    2059              : !! pawcprj_mpi_allgather
    2060              : !!
    2061              : !! FUNCTION
    2062              : !! Perform MPI_ALLGATHER on a pawcprj_type inside a MPI communicator.
    2063              : !!
    2064              : !! INPUTS
    2065              : !!  cprj_loc= The cprj on the local proc being all-gathered
    2066              : !!  natom=Number of atoms (size of first dimension of cprj_loc).
    2067              : !!  n2dim=Size of the second dimension of cprj_loc.
    2068              : !!  n2std=Stride of n2 dimension
    2069              : !!   if n2std=1, cprj_loc(:,1) is on proc 0, cprj_loc(:,2) is on proc 1, cprj_loc(:,3) is on proc 2, etc.
    2070              : !!   if n2std>1, cprj_loc(:,1:n2std) are on proc 0, cprj_loc(:,n2std+1,2*n2std) are on proc 1, etc.
    2071              : !!  nlmn(natom)=Number of nlm partial waves for each atom.
    2072              : !!  ncpgr = number of gradients in cprj_loc
    2073              : !!  nproc=number of processors being gathered
    2074              : !!  spaceComm=MPI Communicator.
    2075              : !!  [rank_ordered]= optional, default=FALSE
    2076              : !!                  TRUE: second dimension of gathered datastructure is rank-ordered
    2077              : !!                  FALSE: second dimension of gathered datastructure is not rank-ordered
    2078              : !!
    2079              : !! OUTPUT
    2080              : !!  cprj_gat=the gathered cprjs
    2081              : !!  ierr=Error status.
    2082              : !!
    2083              : !! SOURCE
    2084              : 
    2085        30042 : subroutine pawcprj_mpi_allgather(cprj_loc,cprj_gat,natom,n2dim,n2std,nlmn,ncpgr,nproc,spaceComm,ierr,&
    2086              : &                                rank_ordered)
    2087              : 
    2088              : !Arguments ------------------------------------
    2089              : !scalars
    2090              :  integer,intent(in) :: natom,n2dim,n2std,ncpgr,nproc,spaceComm
    2091              :  integer,intent(out) :: ierr
    2092              :  logical,optional,intent(in) :: rank_ordered
    2093              : !arrays
    2094              :  integer,intent(in) :: nlmn(natom)
    2095              :  type(pawcprj_type),intent(in) :: cprj_loc(:,:)
    2096              :  type(pawcprj_type),intent(inout) :: cprj_gat(:,:)
    2097              : 
    2098              : !Local variables-------------------------------
    2099              : !scalars
    2100              :  integer :: iat,ii,jj,t2dim,tcpgr,tg2dim,n1dim,nn
    2101              :  integer :: ntotcp,ibuf,ipck,iproc
    2102              :  logical :: rank_ordered_
    2103              :  character(len=100) :: msg
    2104              : !arrays
    2105        30042 :  real(dp),allocatable :: buffer_cpgr(:,:,:),buffer_cpgr_all(:,:,:)
    2106              : 
    2107              : ! *************************************************************************
    2108              : 
    2109        30042 :  n1dim=0
    2110        30042 :  t2dim=0
    2111        30042 :  tg2dim=0
    2112        30042 :  tcpgr=0
    2113        30042 :  ierr=0
    2114              : 
    2115        30042 :  nn=size(nlmn,dim=1)
    2116        30042 :  n1dim=size(cprj_loc,dim=1)
    2117        30042 :  t2dim=size(cprj_loc,dim=2)
    2118        30042 :  tg2dim=size(cprj_gat,dim=2)
    2119        30042 :  tcpgr=cprj_loc(1,1)%ncpgr
    2120              : 
    2121        30042 :  if (nn/=n1dim) then
    2122            0 :    msg='size mismatch in natom (pawcprj_mpi_allgather)!'
    2123            0 :    LIBPAW_BUG(msg)
    2124              :  end if
    2125        30042 :  if (t2dim/=n2dim) then
    2126            0 :    msg='size mismatch in dim=2 (pawcprj_mpi_allgather)!'
    2127            0 :    LIBPAW_BUG(msg)
    2128              :  end if
    2129        30042 :  if (tg2dim/=n2dim*nproc) then
    2130            0 :    msg='size mismatch in dim=2 (pawcprj_mpi_allgather)!'
    2131            0 :    LIBPAW_BUG(msg)
    2132              :  end if
    2133        30042 :  if (tcpgr/=ncpgr) then
    2134            0 :    msg='size mismatch in ncpgr (pawcprj_mpi_allgather)!'
    2135            0 :    LIBPAW_BUG(msg)
    2136              :  end if
    2137        30042 :  if (mod(n2dim,n2std)/=0) then
    2138            0 :    msg='n2std should divide n2dim (pawcprj_mpi_allgather)!'
    2139            0 :    LIBPAW_BUG(msg)
    2140              :  end if
    2141              : 
    2142        30042 :  rank_ordered_=.false.;if(present(rank_ordered)) rank_ordered_=rank_ordered
    2143              : 
    2144        90126 :  ntotcp=n2dim*SUM(nlmn(:))
    2145       120168 :  LIBPAW_ALLOCATE(buffer_cpgr,(2,1+ncpgr,ntotcp))
    2146       120168 :  LIBPAW_ALLOCATE(buffer_cpgr_all,(2,1+ncpgr,nproc*ntotcp))
    2147              : 
    2148              : !=== Pack cprj_loc ====
    2149        30042 :  ipck=0
    2150       150250 :  do jj=1,n2dim
    2151       390666 :    do iat=1,natom
    2152       240416 :      nn=nlmn(iat)
    2153      6010400 :      buffer_cpgr(:,1,ipck+1:ipck+nn)=cprj_loc(iat,jj)%cp(:,1:nn)
    2154      9211168 :      if (ncpgr/=0) buffer_cpgr(:,2:1+ncpgr,ipck+1:ipck+nn)=cprj_loc(iat,jj)%dcp(:,:,1:nn)
    2155       360624 :      ipck=ipck+nn
    2156              :    end do
    2157              :  end do
    2158              : 
    2159              : !=== allgather data ===
    2160        30042 :  call xmpi_allgather(buffer_cpgr,2*(ncpgr+1)*ntotcp,buffer_cpgr_all,spaceComm,ierr)
    2161              : 
    2162              : !=== unpack gathered data into cprj(natom,n2dim*nproc)
    2163              : !=== second dimension is rank-ordered if rank_ordered_=true
    2164        30042 :  ipck=0
    2165        93444 :  do iproc=1,nproc
    2166       231084 :    do jj=1,n2dim/n2std
    2167       454650 :      do ii=1,n2std
    2168       253608 :        if (rank_ordered_) then
    2169        98024 :          ibuf=(iproc-1)*n2dim+(jj-1)*n2std+ii
    2170              :        else
    2171       155584 :          ibuf=(iproc+(jj-1)*nproc-1)*n2std+ii
    2172              :        end if
    2173       898464 :        do iat=1,natom
    2174       507216 :          nn=nlmn(iat)
    2175     12680400 :          cprj_gat(iat,ibuf)%cp(:,1:nn)=buffer_cpgr_all(:,1,ipck+1:ipck+nn)
    2176       507216 :          if (ncpgr/=0) cprj_gat(iat,ibuf)%dcp(:,1:ncpgr,1:nn)=&
    2177     20761200 : &            buffer_cpgr_all(:,2:1+ncpgr,ipck+1:ipck+nn)
    2178       760824 :          ipck=ipck+nn
    2179              :        end do
    2180              :      end do
    2181              :    end do
    2182              :  end do
    2183              : 
    2184              : !=== Clean up ===
    2185        30042 :  LIBPAW_DEALLOCATE(buffer_cpgr)
    2186        30042 :  LIBPAW_DEALLOCATE(buffer_cpgr_all)
    2187              : 
    2188        30042 : end subroutine pawcprj_mpi_allgather
    2189              : !!***
    2190              : 
    2191              : !----------------------------------------------------------------------
    2192              : 
    2193              : !!****f* m_pawcprj/pawcprj_bcast
    2194              : !! NAME
    2195              : !! pawcprj_bcast
    2196              : !!
    2197              : !! FUNCTION
    2198              : !! Broadcast a pawcprj_type from master to all nodes inside a MPI communicator.
    2199              : !!
    2200              : !! INPUTS
    2201              : !!  natom=Number of atoms (size of the first dimension of Cprj).
    2202              : !!  n2dim=Size of the second dimension of Cprj.
    2203              : !!  ncpgr=Number of gradients that have to be cast. It is a bit redundant but, it can be used to
    2204              : !!   broad cast only the %cp"s without caring about the gradients. Just set it to 0 but be careful!
    2205              : !!  nlmn(natom)=Number of nlm partial waves for each atom.
    2206              : !!  master=ID of the sending node in spaceComm.
    2207              : !!  spaceComm=MPI Communicator.
    2208              : !!
    2209              : !! OUTPUT
    2210              : !!  ierr=Error status.
    2211              : !!  Cprj(natom,n2dim)<pawcprj_type>=The datatype to be transmitted by master and received by the others nodes.
    2212              : !!
    2213              : !! SOURCE
    2214              : 
    2215         2470 : subroutine pawcprj_bcast(Cprj,natom,n2dim,nlmn,ncpgr,master,spaceComm,ierr)
    2216              : 
    2217              : !Arguments ------------------------------------
    2218              : !scalars
    2219              :  integer,intent(in) :: natom,n2dim,ncpgr,master,spaceComm
    2220              :  integer,intent(out) :: ierr
    2221              : !arrays
    2222              :  integer,intent(in) :: nlmn(natom)
    2223              :  type(pawcprj_type),intent(inout) :: Cprj(natom,n2dim)
    2224              : 
    2225              : !Local variables-------------------------------
    2226              : !scalars
    2227              :  integer :: iat,jj,n1dim,nn
    2228              :  integer :: ntotcp,ipck,rank,nprocs
    2229              :  character(len=100) :: msg
    2230              : !arrays
    2231         2470 :  real(dp),allocatable :: buffer_cp(:,:),buffer_cpgr(:,:,:)
    2232              : 
    2233              : ! *************************************************************************
    2234              : 
    2235         2470 :  ierr=0
    2236         2470 :  nprocs = xmpi_comm_size(spaceComm)
    2237         2470 :  if (nprocs==1) return
    2238              : 
    2239            0 :  rank = xmpi_comm_rank(spaceComm)
    2240              : 
    2241            0 :  nn=size(nlmn,dim=1)
    2242            0 :  n1dim=size(Cprj,dim=1)
    2243              :  if (nn/=n1dim) then
    2244              :    msg='size mismatch in natom (pawcprj_bcast)!'
    2245              :    LIBPAW_BUG(msg)
    2246              :  end if
    2247              : 
    2248            0 :  ntotcp=n2dim*SUM(nlmn(:))
    2249              : 
    2250            0 :  LIBPAW_ALLOCATE(buffer_cp,(2,ntotcp))
    2251            0 :  if (ncpgr/=0)  then
    2252            0 :    LIBPAW_ALLOCATE(buffer_cpgr,(2,ncpgr,ntotcp))
    2253              :  end if
    2254              : 
    2255              : !=== Master packs Cprj ===
    2256              : !Write a routine to pack/unpack?
    2257            0 :  if (rank==master) then
    2258              :    ipck=0
    2259            0 :    do jj=1,n2dim
    2260            0 :      do iat=1,natom
    2261            0 :        nn=nlmn(iat)
    2262            0 :        buffer_cp(:,ipck+1:ipck+nn)=Cprj(iat,jj)%cp(:,1:nn)
    2263            0 :        if (ncpgr/=0) buffer_cpgr(:,:,ipck+1:ipck+nn)=Cprj(iat,jj)%dcp(:,:,1:nn)
    2264            0 :        ipck=ipck+nn
    2265              :      end do
    2266              :    end do
    2267              :  end if
    2268              : 
    2269              : !=== Transmit data ===
    2270            0 :  call xmpi_bcast(buffer_cp,master,spaceComm,ierr)
    2271            0 :  if (ncpgr/=0) then
    2272            0 :    call xmpi_bcast(buffer_cpgr,master,spaceComm,ierr)
    2273              :  end if
    2274              : 
    2275              : !=== UnPack the received buffer ===
    2276            0 :  if (rank/=master) then
    2277              :    ipck=0
    2278            0 :    do jj=1,n2dim
    2279            0 :      do iat=1,natom
    2280            0 :        nn=nlmn(iat)
    2281            0 :        Cprj(iat,jj)%cp(:,1:nn)=buffer_cp(:,ipck+1:ipck+nn)
    2282            0 :        if (ncpgr/=0) Cprj(iat,jj)%dcp(:,:,1:nn)=buffer_cpgr(:,:,ipck+1:ipck+nn)
    2283            0 :        ipck=ipck+nn
    2284              :      end do
    2285              :    end do
    2286              :  end if
    2287              : 
    2288            0 :  LIBPAW_DEALLOCATE(buffer_cp)
    2289            0 :  if (ncpgr/=0)  then
    2290            0 :    LIBPAW_DEALLOCATE(buffer_cpgr)
    2291              :  end if
    2292              : 
    2293            0 : end subroutine pawcprj_bcast
    2294              : !!***
    2295              : 
    2296              : !----------------------------------------------------------------------
    2297              : 
    2298              : !!****f* m_pawcprj/pawcprj_transpose
    2299              : !! NAME
    2300              : !! pawcprj_transpose
    2301              : !!
    2302              : !! FUNCTION
    2303              : !! Transpose a cprj datastructure FOR A GIVEN (K,SPIN)
    2304              : !! in order to change the parallel distribution from atom to band (or the contrary).
    2305              : !! At input, cprj is distributed over bands (or atoms); at output, it is distributed over atoms (or bands)
    2306              : !!
    2307              : !! INPUTS
    2308              : !!  cprjin(n1indim,n2indim)<pawcprj_type>=the input cprj datastructure
    2309              : !!  cprj_bandpp=number of bands to be treated simultaneoulsy by a processor
    2310              : !!  natom=number of atoms in cell
    2311              : !!  nband=number of bands
    2312              : !!  nspinor=number of spinorial components
    2313              : !!  spaceComm=MPI Communicator.
    2314              : !!
    2315              : !! OUTPUT
    2316              : !!  cprjout(n1outdim,n2outdim)<pawcprj_type>=the output cprj datastructure with another distribution
    2317              : !!
    2318              : !! NOTES
    2319              : !!  On the dimensions:
    2320              : !!   To transfer cprj from band distribution to atom distribution, dimensions should be:
    2321              : !!    n1indim =natom       n2indim =nband/nproc*nspinor
    2322              : !!    n1outdim=natom/nproc n2outdim=nband*nspinor
    2323              : !!   To transfer cprj from atom distribution to band distribution, dimensions should be:
    2324              : !!    n1indim =natom       n2indim =nband/nproc*nspinor
    2325              : !!    n1outdim=natom/nproc n2outdim=nband*nspinor
    2326              : !!
    2327              : !! SOURCE
    2328              : 
    2329            0 :  subroutine pawcprj_transpose(cprjin,cprjout,cprj_bandpp,natom,nband,nspinor,spaceComm)
    2330              : 
    2331              : !Arguments-------------------------------------
    2332              : !scalars
    2333              :  integer :: cprj_bandpp,natom,nband,nspinor,spaceComm
    2334              : !arrays
    2335              :  type(pawcprj_type),intent(in) ::  cprjin(:,:)
    2336              :  type(pawcprj_type),intent(out) :: cprjout(:,:)
    2337              : 
    2338              : !Local variables-------------------------------
    2339              : !scalars
    2340              :  integer :: bpp,buf_indx
    2341              :  integer :: iashft,iatom,iatom_max_sd,iatom_max_rc,iatom_1,iatom_2,iatm1_sd,iatm1_rc,iatm2_sd,iatm2_rc
    2342              :  integer :: ib,iband,iband_1,iband_2,iband_shift,iblock_atom,iblock_band,ibshft
    2343              :  integer :: ierr,ip,ispinor,me,nba,nbb,nbnp_sd,nbnp_rc,ncpgr,nlmn,np
    2344              :  integer :: rbufsize,sbufsize,size11,size12,size21,size22,transpose_mode
    2345              :  character(len=100) :: msg
    2346              : !arrays
    2347            0 :  integer,allocatable :: cprjsz_atom(:),cprjsz_block(:,:)
    2348            0 :  integer,allocatable,target :: count_atom(:),count_band(:),displ_atom(:),displ_band(:)
    2349            0 :  integer,pointer :: scount(:),sdispl(:),rcount(:),rdispl(:)
    2350            0 :  real(dp),allocatable :: rbuf(:),sbuf(:)
    2351              : 
    2352              : ! *************************************************************************
    2353              : 
    2354              : !MPI data
    2355            0 :  me = xmpi_comm_rank(spaceComm)
    2356            0 :  np = xmpi_comm_size(spaceComm)
    2357              : 
    2358              : !Nothing to do if nprocs=1
    2359            0 :  if (np==1) then
    2360            0 :    call pawcprj_copy(cprjin,cprjout)
    2361            0 :    return
    2362              :  end if
    2363              : 
    2364              : !Compute bloc sizes
    2365            0 :  bpp=cprj_bandpp
    2366            0 :  nba=natom/np;if (mod(natom,np)/=0) nba=nba+1
    2367            0 :  nbb=nband/(np*bpp)
    2368              : 
    2369              : !Check sizes, select direction of transposition
    2370            0 :  transpose_mode=0
    2371            0 :  size11=size(cprjin,1);size12=size(cprjin,2)
    2372            0 :  size21=size(cprjout,1);size22=size(cprjout,2)
    2373              :  if (size11==natom.and.size12==nbb*bpp*nspinor.and.&
    2374            0 : & size21==nba.and.size22==nband*nspinor) then
    2375              :    transpose_mode=1
    2376              :  else if (size11==nba.and.size12==nband*nspinor.and.&
    2377            0 : &   size21==natom.and.size22==nbb*bpp*nspinor) then
    2378              :  else
    2379            0 :    msg='wrong cprjin/cprjout sizes (pawcprj_transpose)!'
    2380            0 :    LIBPAW_BUG(msg)
    2381              :  end if
    2382              : 
    2383              : !Compute size of atom bloc (wr to cprj)
    2384            0 :  LIBPAW_ALLOCATE(cprjsz_atom,(natom))
    2385            0 :  LIBPAW_ALLOCATE(cprjsz_block,(np,nba))
    2386            0 :  cprjsz_atom=0;cprjsz_block=0
    2387            0 :  if (transpose_mode==1) then
    2388            0 :    do iatom=1,natom
    2389            0 :      cprjsz_atom(iatom)=2*cprjin(iatom,1)%nlmn*(1+cprjin(iatom,1)%ncpgr)
    2390              :    end do
    2391              :  else
    2392            0 :    do iblock_atom=1,nba
    2393            0 :      iatom=(iblock_atom-1)*np+1+me
    2394            0 :      if (iatom<=natom) cprjsz_atom(iatom)=2*cprjin(iblock_atom,1)%nlmn*(1+cprjin(iblock_atom,1)%ncpgr)
    2395              :    end do
    2396            0 :    call xmpi_sum(cprjsz_atom,spaceComm,ierr)
    2397              :  end if
    2398            0 :  do iblock_atom=1,nba
    2399            0 :    iashft=(iblock_atom-1)*np
    2400            0 :    iatom_1=iashft+1;iatom_2=iashft+np
    2401            0 :    if (iatom_1>natom) cycle
    2402            0 :    if (iatom_2>natom) iatom_2=natom
    2403            0 :    do iatom=iatom_1,iatom_2
    2404            0 :      cprjsz_block(iatom-iashft,iblock_atom)=cprjsz_atom(iatom)+2  ! +2 for nlmn et ncpgr
    2405              :    end do
    2406              :  end do
    2407            0 :  LIBPAW_DEALLOCATE(cprjsz_atom)
    2408              : 
    2409              : !Allocations for MPI_ALLTOALL
    2410            0 :  LIBPAW_ALLOCATE(count_atom,(np))
    2411            0 :  LIBPAW_ALLOCATE(displ_atom,(np))
    2412            0 :  LIBPAW_ALLOCATE(count_band,(np))
    2413            0 :  LIBPAW_ALLOCATE(displ_band,(np))
    2414              : 
    2415              : !Loop on blocks of bands
    2416            0 :  do iblock_band=1,nbb !(note: np divides nband)
    2417            0 :    ibshft=(iblock_band-1)*np*bpp
    2418            0 :    iband_1=ibshft+1;iband_2=ibshft+np*bpp
    2419            0 :    if (iband_1>nband.or.iband_2>nband) cycle ! for security
    2420              : 
    2421              : !  Loop on blocks of atoms
    2422            0 :    do iblock_atom=1,nba
    2423            0 :      iashft=(iblock_atom-1)*np
    2424            0 :      iatom_1=iashft+1;iatom_2=iashft+np
    2425            0 :      if (iatom_1>natom) cycle
    2426            0 :      if (iatom_2>natom) iatom_2=natom
    2427              : 
    2428              : !    Computation of displacements and sizes of blocks when data are band-distributed
    2429            0 :      count_band(1)=cprjsz_block(1,iblock_atom)*nspinor*bpp;displ_band(1)=0
    2430            0 :      do ip=2,np
    2431            0 :        count_band(ip)=cprjsz_block(ip,iblock_atom)*nspinor*bpp
    2432            0 :        displ_band(ip)=displ_band(ip-1)+count_band(ip-1)
    2433              :      end do
    2434              : 
    2435              : !    Computation of displacements and sizes of blocks when data are atom-distributed
    2436            0 :      count_atom(1)=cprjsz_block(1+me,iblock_atom)*bpp*nspinor;displ_atom(1)=0
    2437            0 :      do ip=2,np
    2438            0 :        count_atom(ip)=count_atom(1)
    2439            0 :        displ_atom(ip)=displ_atom(ip-1)+count_atom(ip-1)
    2440              :      end do
    2441              : 
    2442              : !    According to transposition mode, select
    2443              : !    - displacements and sizes of blocks
    2444              : !    - shifts in arrays
    2445            0 :      if (transpose_mode==1) then
    2446            0 :        scount => count_band ; sdispl => displ_band
    2447            0 :        rcount => count_atom ; rdispl => displ_atom
    2448            0 :        nbnp_sd=bpp;nbnp_rc=np*bpp
    2449            0 :        iatm1_sd=iatom_1;iatm2_sd=iatom_2
    2450            0 :        iatm1_rc=iblock_atom;iatm2_rc=iatm1_rc
    2451            0 :        iatom_max_sd=iatom_2;iatom_max_rc=iashft+1+me
    2452              :      else
    2453            0 :        scount => count_atom ; sdispl => displ_atom
    2454            0 :        rcount => count_band ; rdispl => displ_band
    2455            0 :        nbnp_sd=np*bpp;nbnp_rc=bpp
    2456            0 :        iatm1_sd=iblock_atom;iatm2_sd=iatm1_sd
    2457            0 :        iatm1_rc=iatom_1;iatm2_rc=iatom_2
    2458            0 :        iatom_max_sd=iashft+1+me;iatom_max_rc=iatom_2
    2459              :      end if
    2460              : 
    2461              : !    Allocation of buffers
    2462            0 :      sbufsize=sdispl(np)+scount(np)
    2463            0 :      rbufsize=rdispl(np)+rcount(np)
    2464            0 :      LIBPAW_ALLOCATE(sbuf,(sbufsize))
    2465            0 :      LIBPAW_ALLOCATE(rbuf,(rbufsize))
    2466              : 
    2467              : !    Coying of input cprj to buffer for sending
    2468            0 :      buf_indx=0
    2469            0 :      iband_shift=(iblock_band-1)*nbnp_sd-1
    2470            0 :      if (iatom_max_sd<=natom) then
    2471            0 :        do iatom=iatm1_sd,iatm2_sd
    2472            0 :          do ib=1,nbnp_sd
    2473            0 :            iband=(iband_shift+ib)*nspinor
    2474            0 :            do ispinor=1,nspinor
    2475            0 :              iband=iband+1
    2476            0 :              nlmn=cprjin(iatom,iband)%nlmn;ncpgr=cprjin(iatom,iband)%ncpgr
    2477            0 :              sbuf(buf_indx+1)=dble(nlmn) ;buf_indx=buf_indx+1
    2478            0 :              sbuf(buf_indx+1)=dble(ncpgr);buf_indx=buf_indx+1
    2479            0 :              sbuf(buf_indx+1:buf_indx+2*nlmn)=reshape(cprjin(iatom,iband)%cp(1:2,1:nlmn),(/2*nlmn/))
    2480            0 :              buf_indx=buf_indx+2*nlmn
    2481            0 :              if (ncpgr>0) then
    2482            0 :                sbuf(buf_indx+1:buf_indx+2*ncpgr*nlmn)=reshape(cprjin(iatom,iband)%dcp(1:2,1:ncpgr,1:nlmn),(/2*ncpgr*nlmn/))
    2483            0 :                buf_indx=buf_indx+2*ncpgr*nlmn
    2484              :              end if
    2485              :            end do
    2486              :          end do
    2487              :        end do
    2488              :      end if
    2489            0 :      if (buf_indx/=sbufsize) then
    2490            0 :        msg='wrong buffer size for sending (pawcprj_transpose)!'
    2491            0 :        LIBPAW_BUG(msg)
    2492              :      end if
    2493              : 
    2494              : !    Main call to MPI_ALLTOALL
    2495            0 :      call xmpi_alltoallv(sbuf,scount,sdispl,rbuf,rcount,rdispl,spaceComm,ierr)
    2496              : 
    2497              : !    Retrieving of output cprj for received buffer
    2498            0 :      buf_indx=0
    2499            0 :      iband_shift=(iblock_band-1)*nbnp_rc-1
    2500            0 :      if (iatom_max_rc<=natom) then
    2501            0 :        do iatom=iatm1_rc,iatm2_rc
    2502            0 :          do ib=1,nbnp_rc
    2503            0 :            iband=(iband_shift+ib)*nspinor
    2504            0 :            do ispinor=1,nspinor
    2505            0 :              iband=iband+1
    2506            0 :              nlmn =int(rbuf(buf_indx+1));buf_indx=buf_indx+1
    2507            0 :              ncpgr=int(rbuf(buf_indx+1));buf_indx=buf_indx+1
    2508            0 :              cprjout(iatom,iband)%nlmn=nlmn;cprjout(iatom,iband)%ncpgr=ncpgr
    2509            0 :              cprjout(iatom,iband)%cp(1:2,1:nlmn)=reshape(rbuf(buf_indx+1:buf_indx+2*nlmn),(/2,nlmn/))
    2510            0 :              buf_indx=buf_indx+2*nlmn
    2511            0 :              if (ncpgr>0) then
    2512            0 :                cprjout(iatom,iband)%dcp(1:2,1:ncpgr,1:nlmn)=reshape(rbuf(buf_indx+1:buf_indx+2*nlmn*ncpgr),(/2,ncpgr,nlmn/))
    2513            0 :                buf_indx=buf_indx+2*nlmn*ncpgr
    2514              :              end if
    2515              :            end do
    2516              :          end do
    2517              :        end do
    2518              :      else
    2519            0 :        cprjout(iatom,iband)%nlmn=0;cprjout(iatom,iband)%ncpgr=0
    2520              :      end if
    2521            0 :      if (buf_indx/=rbufsize) then
    2522            0 :        msg='wrong buffer size for receiving (pawcprj_transpose)!'
    2523            0 :        LIBPAW_BUG(msg)
    2524              :      end if
    2525              : 
    2526              : !    Deallocation of buffers
    2527            0 :      LIBPAW_DEALLOCATE(sbuf)
    2528            0 :      LIBPAW_DEALLOCATE(rbuf)
    2529              : 
    2530              : !    End of loops
    2531              :    end do ! do iblock_atom
    2532              :  end do ! do iblock_atom
    2533              : 
    2534              : !Free memory
    2535            0 :  LIBPAW_DEALLOCATE(count_atom)
    2536            0 :  LIBPAW_DEALLOCATE(displ_atom)
    2537            0 :  LIBPAW_DEALLOCATE(count_band)
    2538            0 :  LIBPAW_DEALLOCATE(displ_band)
    2539            0 :  LIBPAW_DEALLOCATE(cprjsz_block)
    2540            0 :  nullify(scount,rcount,sdispl,rdispl)
    2541              : 
    2542            0 :  end subroutine pawcprj_transpose
    2543              : !!***
    2544              : 
    2545              : !----------------------------------------------------------------------
    2546              : 
    2547              : !!****f* m_pawcprj/pawcprj_gather_spin
    2548              : !! NAME
    2549              : !! pawcprj_gather_spin
    2550              : !!
    2551              : !! FUNCTION
    2552              : !!
    2553              : !! INPUTS
    2554              : !!  cprj(:,:)=the input cprj datastructure
    2555              : !!  n2size=number of cprj datastructures to be gathered (second dim)
    2556              : !!  nspinor : number of spinorial component (on current proc)
    2557              : !!  nspinortot : total number of spinorial component
    2558              : !!
    2559              : !! OUTPUT
    2560              : !!  cprj_gat(:,:) = the cprj containing all nspinor componants
    2561              : !!
    2562              : !! NOTES
    2563              : !! The cprj has been built like the following:
    2564              : !!   loop on nsppol
    2565              : !!   loop on k point
    2566              : !!   loop over band or block of band
    2567              : !! These quantities were build only if treated by the current proc
    2568              : !! the inner quantities being nspinor
    2569              : !!
    2570              : !! SOURCE
    2571          464 :  subroutine pawcprj_gather_spin(cprj,cprj_gat,natom,n2size,nspinor,nspinortot,&
    2572              : &                            spaceComm_spin,ierr)
    2573              : 
    2574              : !Arguments ------------------------------------
    2575              : !scalars
    2576              :  integer,intent(in) :: natom,nspinor,nspinortot,n2size
    2577              :  integer,intent(in) :: spaceComm_spin
    2578              :  integer,intent(out) :: ierr
    2579              : !arrays
    2580              :  type(pawcprj_type),intent(in) :: cprj(:,:)
    2581              :  type(pawcprj_type),intent(inout) :: cprj_gat(:,:)
    2582              : 
    2583              : !Local variables-------------------------------
    2584              : !scalars
    2585              :  integer :: i1,iatom,ibsp,icpgr,ilmn,isp,ispinor,jj,lmndim,n2dim,n2dim_gat,ncpgr
    2586              :  character(len=100) :: msg
    2587              : !arrays
    2588          928 :  integer :: nlmn(natom)
    2589          464 :  real(dp),allocatable :: buffer1(:),buffer2(:)
    2590              : 
    2591              : ! *************************************************************************
    2592              : 
    2593          464 :  n2dim    =size(cprj,dim=2)
    2594          464 :  n2dim_gat=size(cprj_gat,dim=2)
    2595          464 :  if (n2dim_gat/=(nspinortot/nspinor)*n2dim) then
    2596            0 :    msg='wrong dims (pawcprj_gather_spin)!'
    2597            0 :    LIBPAW_BUG(msg)
    2598              :  end if
    2599              : 
    2600         1392 :  do iatom=1,natom
    2601         1392 :    nlmn(iatom)=size(cprj(iatom,1)%cp(1,:))
    2602              :  end do
    2603          464 :  ncpgr=cprj(1,1)%ncpgr
    2604         1392 :  lmndim=2*n2size*sum(nlmn(1:natom))*(1+ncpgr)
    2605         1392 :  LIBPAW_ALLOCATE(buffer1,(lmndim))
    2606         1392 :  LIBPAW_ALLOCATE(buffer2,(lmndim*nspinortot))
    2607              : 
    2608          464 :  isp=0;ibsp=0
    2609          464 :  jj=1
    2610         8432 :  do i1=1,n2size
    2611         7968 :    isp=isp+1
    2612        24368 :    do iatom=1,natom
    2613       302784 :      do ilmn=1,nlmn(iatom)
    2614       860544 :        buffer1(jj:jj+1)=cprj(iatom,isp)%cp(1:2,ilmn)
    2615       302784 :        jj=jj+2
    2616              :      end do
    2617        23904 :      if (ncpgr>0) then
    2618            0 :        do ilmn=1,nlmn(iatom)
    2619            0 :          do icpgr=1,ncpgr
    2620            0 :            buffer1(jj:jj+1)=cprj(iatom,isp)%dcp(1:2,icpgr,ilmn)
    2621            0 :            jj=jj+2
    2622              :          end do
    2623              :        end do
    2624              :      end if
    2625              :    end do
    2626              :  end do
    2627              : 
    2628          464 :  call xmpi_allgather(buffer1,lmndim,buffer2,spaceComm_spin,ierr)
    2629              : 
    2630          464 :  jj=1
    2631         1392 :  do ispinor=1,nspinortot
    2632        17328 :    do i1 =1,n2size
    2633        15936 :      ibsp=(i1-1)*nspinortot + ispinor
    2634        48736 :      do iatom=1,natom
    2635       605568 :        do ilmn=1,nlmn(iatom)
    2636      1721088 :          cprj_gat(iatom,ibsp)%cp(1:2,ilmn)=buffer2(jj:jj+1)
    2637       605568 :          jj=jj+2
    2638              :        end do
    2639        47808 :        if (ncpgr>0) then
    2640            0 :          do ilmn=1,nlmn(iatom)
    2641            0 :            do icpgr=1,ncpgr
    2642            0 :              cprj_gat(iatom,ibsp)%dcp(1:2,icpgr,ilmn)=buffer2(jj:jj+1)
    2643            0 :              jj=jj+2
    2644              :            end do
    2645              :          end do
    2646              :        end if
    2647              :      end do
    2648              :    end do
    2649              :  end do
    2650              : 
    2651          464 :  LIBPAW_DEALLOCATE(buffer1)
    2652          464 :  LIBPAW_DEALLOCATE(buffer2)
    2653              : 
    2654          464 :  end subroutine pawcprj_gather_spin
    2655              : !!***
    2656              : 
    2657              : !----------------------------------------------------------------------
    2658              : 
    2659              : !!****f* m_pawcprj/pawcprj_getdim
    2660              : !! NAME
    2661              : !! pawcprj_getdim
    2662              : !!
    2663              : !! FUNCTION
    2664              : !!  Helper function returning the number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom.
    2665              : !!  Used to initialize the dimensioning array that is passed to the pawcprj_alloc routines when the
    2666              : !!  pawcprj_type structure is allocated and initialized.
    2667              : !!
    2668              : !! INPUTS
    2669              : !! natom=number of atoms in the unit cell
    2670              : !! nattyp(ntypat)=number of atoms of each type
    2671              : !! ntypat=number of atom types
    2672              : !! typat(natom-= type of each atom
    2673              : !! Pawtab(ntypat)<pawtab_type>=PAW tabulated starting data.
    2674              : !! sort_mode(len=*)=String defining the sorting of the atoms in the Cprj arrays.
    2675              : !!   Two modes are possible:
    2676              : !!   -- "O[rdered]", if atoms are sorted by atom type.
    2677              : !!   -- "R[andom]", if atoms are sorted randomly i.e. according the values of typat specified in the input file.
    2678              : !!
    2679              : !! OUTPUT
    2680              : !!  dimcprj(natom)=Number of nlm elements in the <p_{lmn}^i|\psi> matrix elements for i=1,...,natom.
    2681              : !!
    2682              : !! SOURCE
    2683              : 
    2684        42424 : subroutine pawcprj_getdim(dimcprj,natom,nattyp,ntypat,typat,Pawtab,sort_mode)
    2685              : 
    2686              : !Arguments ------------------------------------
    2687              :  integer,intent(in) :: natom,ntypat
    2688              :  character(len=*),intent(in) :: sort_mode
    2689              : !arrays
    2690              :  integer,intent(in) :: nattyp(:),typat(natom)
    2691              :  integer,intent(inout) :: dimcprj(natom)
    2692              :  type(Pawtab_type),intent(in) :: Pawtab(ntypat)
    2693              : 
    2694              : !Local variables-------------------------------
    2695              :  integer :: iatom,itypat
    2696              :  character(len=500) :: msg
    2697              : 
    2698              : ! *************************************************************************
    2699              : 
    2700              :  SELECT CASE (sort_mode(1:1))
    2701              : 
    2702              :  CASE ("o","O") ! Ordered by atom-type
    2703              : 
    2704              :   iatom=0
    2705        92188 :   do itypat=1,ntypat
    2706       149091 :    dimcprj(iatom+1:iatom+nattyp(itypat))=Pawtab(itypat)%lmn_size
    2707        92188 :    iatom=iatom+nattyp(itypat)
    2708              :   end do
    2709              : 
    2710              :  CASE ("r","R") ! Randomly ordered (typat from input file)
    2711              : 
    2712        12525 :   do iatom=1,natom
    2713         7939 :    itypat=typat(iatom)
    2714        12525 :    dimcprj(iatom)=Pawtab(itypat)%lmn_size
    2715              :   end do
    2716              : 
    2717              :  CASE DEFAULT
    2718            0 :   msg='Wrong value for sort_mode: '//TRIM(sort_mode)
    2719        42424 :   LIBPAW_ERROR(msg)
    2720              :  END SELECT
    2721              : 
    2722        42424 : end subroutine pawcprj_getdim
    2723              : !!***
    2724              : 
    2725              : !----------------------------------------------------------------------
    2726              : 
    2727              : !!****f* m_pawcprj/paw_overlap
    2728              : !! NAME
    2729              : !! paw_overlap
    2730              : !!
    2731              : !! FUNCTION
    2732              : !!  Helper function returning the onsite contribution to the overlap between two states.
    2733              : !!
    2734              : !! INPUTS
    2735              : !!   spinor_comm= (optional) communicator over spinorial components
    2736              : !!   typat(:)=The type of each atom.
    2737              : !!   Pawtab(ntypat)<type(pawtab_type)>=paw tabulated starting data.
    2738              : !!   cprj1,cprj2<pawcprj_type>
    2739              : !!     Projected wave functions <Proj_i|Cnk> with all NL projectors for the left and the right wavefunction,respectively.
    2740              : !!
    2741              : !! OUTPUT
    2742              : !!
    2743              : !! SOURCE
    2744              : 
    2745         2310 : function paw_overlap(cprj1,cprj2,typat,pawtab,spinor_comm) result(onsite)
    2746              : 
    2747              : !Arguments ------------------------------------
    2748              : !scalars
    2749              :  integer,intent(in),optional :: spinor_comm
    2750              : !arrays
    2751              :  integer,intent(in) :: typat(:)
    2752              :  real(dp) :: onsite(2)
    2753              :  type(pawcprj_type),intent(in) :: cprj1(:,:),cprj2(:,:)
    2754              :  type(pawtab_type),intent(in) :: pawtab(:)
    2755              : 
    2756              : !Local variables-------------------------------
    2757              : !scalars
    2758              :  integer :: iatom,ilmn,itypat,j0lmn,jlmn,klmn,natom,nspinor,isp
    2759              :  real(dp) :: sij
    2760              :  character(len=500) :: msg
    2761              : !arrays
    2762              : 
    2763              : ! *************************************************************************
    2764              : 
    2765         2310 :  natom=SIZE(typat)
    2766              : 
    2767         2310 :  if (SIZE(cprj1,DIM=1)/=SIZE(cprj2,DIM=1) .or. SIZE(cprj1,DIM=1)/=natom) then
    2768            0 :    write(msg,'(a,3i4)')' Wrong size in typat, cprj1, cprj2 : ',natom,SIZE(cprj1),SIZE(cprj2)
    2769            0 :    LIBPAW_ERROR(msg)
    2770              :  end if
    2771              : 
    2772         2310 :  nspinor = SIZE(cprj1,DIM=2)
    2773              : 
    2774         6930 :  onsite=zero
    2775         4620 :  do iatom=1,natom
    2776         2310 :    itypat=typat(iatom)
    2777        46200 :    do jlmn=1,pawtab(itypat)%lmn_size
    2778        41580 :      j0lmn=jlmn*(jlmn-1)/2
    2779       438900 :      do ilmn=1,jlmn
    2780       395010 :        klmn=j0lmn+ilmn
    2781       395010 :        sij=pawtab(itypat)%sij(klmn); if (jlmn==ilmn) sij=sij*half
    2782       436590 :        if (ABS(sij)>tol16) then
    2783       124740 :          do isp=1,nspinor
    2784              : 
    2785              :            onsite(1)=onsite(1) + sij*(                                 &
    2786              : &            cprj1(iatom,isp)%cp(1,ilmn) * cprj2(iatom,isp)%cp(1,jlmn) &
    2787              : &           +cprj1(iatom,isp)%cp(2,ilmn) * cprj2(iatom,isp)%cp(2,jlmn) &
    2788              : &           +cprj1(iatom,isp)%cp(1,jlmn) * cprj2(iatom,isp)%cp(1,ilmn) &
    2789              : &           +cprj1(iatom,isp)%cp(2,jlmn) * cprj2(iatom,isp)%cp(2,ilmn) &
    2790        62370 : &           )
    2791              : 
    2792              :            onsite(2)=onsite(2) + sij*(                                 &
    2793              : &           cprj1(iatom,isp)%cp(1,ilmn) * cprj2(iatom,isp)%cp(2,jlmn)  &
    2794              : &           -cprj1(iatom,isp)%cp(2,ilmn) * cprj2(iatom,isp)%cp(1,jlmn) &
    2795              : &           +cprj1(iatom,isp)%cp(1,jlmn) * cprj2(iatom,isp)%cp(2,ilmn) &
    2796              : &           -cprj1(iatom,isp)%cp(2,jlmn) * cprj2(iatom,isp)%cp(1,ilmn) &
    2797       124740 : &           )
    2798              :          end do
    2799              :        end if
    2800              :      end do
    2801              :    end do
    2802              :  end do
    2803              : 
    2804         2310 :  if (present(spinor_comm)) then
    2805         2310 :    call xmpi_sum(onsite,spinor_comm,isp)
    2806              :  end if
    2807              : 
    2808         2310 : end function paw_overlap
    2809              : !!***
    2810              : 
    2811              : !----------------------------------------------------------------------
    2812              : 
    2813              : !!****f* m_pawcprj/pawcprj_pack
    2814              : !! NAME
    2815              : !! pawcprj_pack
    2816              : !!
    2817              : !! FUNCTION
    2818              : !! Pack structured data into a simple buffer
    2819              : !!
    2820              : !! INPUTS
    2821              : !!  nlmn(natom)=Number of nlm partial waves for each atom.
    2822              : !!  ncpgr = number of gradients in cprj_out
    2823              : !!  cprj= The datatype to be packed.
    2824              : !!
    2825              : !! OUTPUT
    2826              : !!  buffer = the data packed, dim : (2, n2dim*sum(nlmn))
    2827              : !!  [buffer_gr] = if present the gradient data packed, dim : (2, ncpgr, n2dim*sum(nlmn))
    2828              : !!
    2829              : !! SOURCE
    2830              : 
    2831          600 : subroutine pawcprj_pack(nlmn,cprj,buffer,buffer_gr)
    2832              : 
    2833              : !Arguments ------------------------------------
    2834              : !scalars
    2835              : !arrays
    2836              :  integer,intent(in) :: nlmn(:)
    2837              :  type(pawcprj_type),intent(in) :: cprj(:,:)
    2838              :  real(dp),intent(out) :: buffer(:,:)
    2839              :  real(dp),intent(out),optional :: buffer_gr(:,:,:)
    2840              : 
    2841              : !Local variables-------------------------------
    2842              : !scalars
    2843              :  integer :: natom,n2buffer,ncpgr,n2dim
    2844              :  integer :: iat,jj,n1dim,nn
    2845              :  integer :: ipck
    2846              :  character(len=100) :: msg
    2847              : !arrays
    2848              : 
    2849              : ! *************************************************************************
    2850              : 
    2851          600 :  natom=size(nlmn,dim=1)
    2852          600 :  n2buffer=size(buffer,dim=2)
    2853          600 :  n1dim=size(cprj,dim=1)
    2854          600 :  n2dim=size(cprj,dim=2)
    2855              : 
    2856          600 :  if (natom/=n1dim) then
    2857            0 :    msg='size mismatch in natom (pawcprj_pack)!'
    2858            0 :    LIBPAW_BUG(msg)
    2859              :  end if
    2860         3600 :  if (n2dim*SUM(nlmn)/=n2buffer) then
    2861            0 :    msg='size mismatch in dim=2 (pawcprj_pack)!'
    2862            0 :    LIBPAW_BUG(msg)
    2863              :  end if
    2864          600 :  ncpgr=0
    2865          600 :  if (present(buffer_gr)) then
    2866            0 :    ncpgr=size(buffer_gr,dim=2)
    2867              :  end if
    2868              : 
    2869              : !=== Pack cprj ====
    2870          600 :  ipck=0
    2871         1200 :  do jj=1,n2dim
    2872         4200 :    do iat=1,natom
    2873         3000 :      nn=nlmn(iat)
    2874       118200 :      buffer(:,ipck+1:ipck+nn)=cprj(iat,jj)%cp(:,1:nn)
    2875         3000 :      if (ncpgr/=0) then
    2876            0 :        buffer_gr(:,:,ipck+1:ipck+nn)=cprj(iat,jj)%dcp(:,:,1:nn)
    2877              :      end if
    2878         3600 :      ipck=ipck+nn
    2879              :    end do
    2880              :  end do
    2881              : 
    2882          600 : end subroutine pawcprj_pack
    2883              : !!***
    2884              : 
    2885              : !----------------------------------------------------------------------
    2886              : 
    2887              : !!****f* m_pawcprj/pawcprj_unpack
    2888              : !! NAME
    2889              : !! pawcprj_unpack
    2890              : !!
    2891              : !! FUNCTION
    2892              : !! Unpack structured data from a simple buffer
    2893              : !!
    2894              : !! INPUTS
    2895              : !!  nlmn(natom)=Number of nlm partial waves for each atom.
    2896              : !!  ncpgr = number of gradients in cprj_in
    2897              : !!  buffer = the data to be unpacked, dim : (2, n2dim*sum(nlmn))
    2898              : !!  [buffer_gr] = if present the gradient data to be unpacked, dim : (2, ncpgr, n2dim*sum(nlmn))
    2899              : !!
    2900              : !! OUTPUT
    2901              : !!  cprj=The datatype unpacked
    2902              : !!
    2903              : !! SOURCE
    2904              : 
    2905         1508 : subroutine pawcprj_unpack(nlmn,cprj,buffer,buffer_gr)
    2906              : 
    2907              : !Arguments ------------------------------------
    2908              : !scalars
    2909              : !arrays
    2910              :  integer,intent(in) :: nlmn(:)
    2911              :  real(dp),intent(in) :: buffer(:,:)
    2912              :  real(dp),intent(in),optional :: buffer_gr(:,:,:)
    2913              :  type(pawcprj_type),intent(inout) :: cprj(:,:)
    2914              : 
    2915              : !Local variables-------------------------------
    2916              : !scalars
    2917              :  integer :: natom,n2buffer,ncpgr,n2dim
    2918              :  integer :: iat,jj,n1dim,nn
    2919              :  integer :: ipck
    2920              :  character(len=100) :: msg
    2921              : !arrays
    2922              : 
    2923              : ! *************************************************************************
    2924              : 
    2925         1508 :  natom=size(nlmn,dim=1)
    2926         1508 :  n2buffer=size(buffer,dim=2)
    2927         1508 :  n1dim=size(cprj,dim=1)
    2928         1508 :  n2dim=size(cprj,dim=2)
    2929              : 
    2930         1508 :  if (natom/=n1dim) then
    2931            0 :    msg='size mismatch in natom (pawcprj_unpack)!'
    2932            0 :    LIBPAW_BUG(msg)
    2933              :  end if
    2934         9048 :  if (n2dim*SUM(nlmn)/=n2buffer) then
    2935            0 :    msg='size mismatch in dim=2 (pawcprj_unpack)!'
    2936            0 :    LIBPAW_BUG(msg)
    2937              :  end if
    2938         1508 :  ncpgr=0
    2939         1508 :  if (present(buffer_gr)) then
    2940            0 :    ncpgr=size(buffer_gr,dim=2)
    2941              :  end if
    2942              : 
    2943              : !=== Unpack buffers into cprj ===
    2944         1508 :  ipck=0
    2945         3016 :  do jj=1,n2dim
    2946        10556 :    do iat=1,natom
    2947         7540 :      nn=nlmn(iat)
    2948       297076 :      cprj(iat,jj)%cp(:,1:nn)=buffer(:,ipck+1:ipck+nn)
    2949         7540 :      if (ncpgr/=0) then
    2950            0 :        cprj(iat,jj)%dcp(:,:,1:nn)=buffer_gr(:,:,ipck+1:ipck+nn)
    2951              :      end if
    2952         9048 :      ipck=ipck+nn
    2953              :    end do
    2954              :  end do
    2955              : 
    2956         1508 : end subroutine pawcprj_unpack
    2957              : 
    2958            0 : end module m_pawcprj
    2959              : !!***
        

Generated by: LCOV version 2.3-1