TABLE OF CONTENTS


ABINIT/ctocprj [ Functions ]

[ Top ] [ Functions ]

NAME

 ctocprj

FUNCTION

  Compute all <Proj_i|Cnk> for every wave function |Cnk> expressed in reciprocal space.
  |Proj_i> are non-local projectors (for each atom and each l,m,n)
  Can also compute derivatives of <Proj_i|Cnk> wrt to several parameters

INPUTS

  atindx(natom)=index table for atoms
  cg(2,mcg)=planewave coefficients of wavefunctions
  choice: chooses derivatives to compute:
          =1 => no derivatives
          =2 => 1st derivatives with respect to atomic position(s)
          =3 => 1st derivatives with respect to strain(s)
          =23=> 1st derivatives with respect to strain(s) and atm pos
          =4 => 2nd derivatives with respect to atomic pos.
          =24=> 1st and 2nd derivatives with respect to atomic pos.
          =5 => derivatives with respect to k wavevector
          =6 => 2nd derivatives with respect to strain and atm. pos.
  gmet(3,3)=reciprocal space metric tensor in bohr**-2
  gprimd(3,3)=dimensional reciprocal space primitive translations
  iatom= if <=0, cprj=<p_i|Cnk> are computed for all atoms 1...natom
         if  >0  cprj=<p_i|Cnk> are computed only for atom with index iatom
  idir=direction of the derivative, i.e. dir. of - atom to be moved  in the case choice=2
                                                 - strain component  in the case choice=3
                                                 - k point direction in the case choice=5
       Compatible only with choice=2,3,5; if idir=0, all derivatives are computed
  iorder_cprj=0 if output cprj=<p_i|Cnk> are sorted by atom type
                (first all elements of atom type 1, followed by those of atom type 2 and so on).
              1 if output cprj=<p_i|Cnk> are sorted according to
                the variable typat in the main input file
  istwfk(nkpt)=option parameter that describes the storage of wfs
  kg(3,mpw*mkmem)=reduced planewave coordinates
  kpt(3,nkpt)=reduced coordinates of k points.
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mgfft=maximum size of 1D FFTs
  mkmem=number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  mpsang=1+maximum angular momentum for nonlocal pseudopotentials
  mpw=maximum dimensioned size of npw
  natom=number of atoms in cell
  nattyp(ntypat)= # atoms of each type
  nband(nkpt*nsppol)=number of bands at this k point for that spin polarization
  ncprj=1st dim. of cprj array (natom if iatom<=0, 1 if iatom>0)
  ngfft(18)=contain all needed information about 3D FFT, see ~ABINIT/Infos/vargs.htm#ngfft
  nkpt=number of k points
  nloalg(3)=governs the choice of the algorithm for nonlocal operator
  npwarr(nkpt)=number of planewaves in basis at this k point
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  ntypat=number of types of atoms in unit cell
  paral_kgb= 1 if kpt-band-FFT is activated
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rmet(3,3)=real space metric (bohr**2)
  tim_ctocprj=timing code of the calling routine
  typat(natom)= types of atoms
  uncp=unit number for <P_lmn|Cnk> data (if used)
  xred(3,natom)=reduced dimensionless atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang)=real spherical harmonics for each G and k point
  ylmgr(npw*mkmem,nylmgr,mpsang*mpsang*useylmgr)=gradients of real spherical harmonics wrt (k+G)

OUTPUT

  cprj(ncprj,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors
                                       Usually ncprj=natom

PARENTS

      dfpt_looppert,extrapwf,forstr,scfcv,update_e_field_vars,vtorho

CHILDREN

      getcprj,mkffnl,mkkpg,pawcprj_alloc,pawcprj_free,pawcprj_mpi_sum
      pawcprj_put,pawcprj_set_zero,ph1d3d,strconv,xmpi_allgather
      xmpi_allgatherv,xmpi_alltoallv

SOURCE

451  subroutine ctocprj(atindx,cg,choice,cprj,gmet,gprimd,iatom,idir,&
452 & iorder_cprj,istwfk,kg,kpt,mcg,mcprj,mgfft,mkmem,mpi_enreg,mpsang,&
453 & mpw,natom,nattyp,nband,ncprj,ngfft,nkpt,nloalg,npwarr,nspinor,&
454 & nsppol,ntypat,paral_kgb,ph1d,psps,rmet,typat,ucvol,uncp,xred,ylm,ylmgr)
455 
456 
457 !This section has been created automatically by the script Abilint (TD).
458 !Do not modify the following lines by hand.
459 #undef ABI_FUNC
460 #define ABI_FUNC 'ctocprj'
461 !End of the abilint section
462 
463  implicit none
464 
465 !Arguments -------------------------------
466 !scalars
467  integer,intent(in) :: choice,iatom,idir,iorder_cprj,mcg,mcprj,mgfft,mkmem,mpsang,mpw
468  integer,intent(in) :: natom,ncprj,nkpt,nspinor,nsppol,ntypat,paral_kgb,uncp
469  real(dp),intent(in) :: ucvol
470  type(MPI_type),intent(in) :: mpi_enreg
471  type(pseudopotential_type),target,intent(in) :: psps
472 !arrays
473  integer,intent(in) :: istwfk(nkpt),nband(nkpt*nsppol)
474  integer,intent(in) :: ngfft(18),nloalg(3),npwarr(nkpt),kg(3,mpw*mkmem),typat(natom)
475  integer,intent(in),target :: atindx(natom),nattyp(ntypat)
476  real(dp),intent(in) :: cg(2,mcg)
477  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt(3,nkpt),rmet(3,3)
478  real(dp),intent(in) :: xred(3,natom),ylm(:,:),ylmgr(:,:,:)
479  real(dp),intent(in),target :: ph1d(2,3*(2*mgfft+1)*natom)
480  type(pawcprj_type),intent(inout) :: cprj(ncprj,mcprj)
481 
482 !Local variables-------------------------------
483 !scalars
484  integer :: blocksz,cg_bandpp,counter,cpopt,cprj_bandpp,dimffnl,ia,iatm,iatom1,iatom2
485  integer :: iband_max,iband_min,iband_start,ibg,ibgb,iblockbd,ibp,icg,icgb,icp1,icp2
486  integer :: ider,idir0,iend,ierr,ig,ii,ikg,ikpt,ilm,ipw,isize,isppol,istart,istwf_k,itypat,iwf1,iwf2,jdir
487  integer :: matblk,mband_cprj,me_distrb,my_nspinor,n1,n1_2p1,n2,n2_2p1,n3,n3_2p1,kk,nlmn
488  integer :: nband_k,nband_cprj_k,nblockbd,ncpgr,nkpg,npband_bandfft,npws,npw_k,npw_nk,ntypat0
489  integer :: shift1,shift1b,shift2,shift2b,shift3,shift3b
490  integer :: spaceComm,spaceComm_band,spaceComm_fft,useylmgr
491  logical :: cg_band_distributed,cprj_band_distributed,one_atom
492  real(dp) :: arg
493  character(len=500) :: msg
494 !arrays
495  integer,allocatable :: bufsize(:),bufsize_wf(:),bufdisp(:),bufdisp_wf(:)
496  integer,allocatable :: dimlmn(:),kg_k(:,:),kg_k_loc(:,:)
497  integer,allocatable :: npw_block(:),npw_disp(:)
498  integer,pointer :: atindx_atm(:),indlmn_atm(:,:,:),nattyp_atm(:),pspso_atm(:)
499  real(dp) :: kpoint(3),work(6)
500  real(dp),allocatable :: cwavef(:,:),cwavef_tmp(:,:)
501  real(dp),allocatable :: ffnl(:,:,:,:),ffnl_npw(:,:,:,:),ffnl_tmp(:,:,:,:),ffnl_tmp_npw(:,:,:,:)
502  real(dp),allocatable :: kpg_k(:,:)
503  real(dp),allocatable :: ph3d(:,:,:),ph3d_npw(:,:,:),ph3d_tmp(:,:,:),ph3d_tmp_npw(:,:,:)
504  real(dp),allocatable :: phkxred(:,:),ylm_k(:,:),ylmgr_k(:,:,:)
505  real(dp),ABI_CONTIGUOUS pointer :: ekb_atm(:,:),ffspl_atm(:,:,:,:),ph1d_atm(:,:)
506  type(pawcprj_type),allocatable :: cwaveprj(:,:)
507 
508 ! *********************************************************************
509 
510  DBG_ENTER('COLL')
511 
512 !Nothing to do if current MPI process does treat kpoints or plane-waves
513  if (mcg==0.or.mcprj==0) return
514 
515 !Preliminary tests
516  if (psps%useylm==0) then
517    msg='Not available for useylm=0!'
518    MSG_ERROR(msg)
519  end if
520  if ((choice<1.or.choice>6).and.choice/=23.and.choice/=24) then
521    msg='Bad choice!'
522    MSG_BUG(msg)
523  end if
524  if (idir>0.and.choice/=2.and.choice/=3.and.choice/=5) then
525    msg='Does not support idir>0 for that choice!'
526    MSG_BUG(msg)
527  end if
528  if (size(ylm)/=mpw*mkmem*mpsang*mpsang) then
529    msg='Wrong size for Ylm!'
530    MSG_BUG(msg)
531  end if
532  useylmgr=merge(0,1,(size(ylmgr)==0))
533  if (useylmgr==0.and.(choice==3.or.choice==5.or.choice==23)) then
534    msg=' Ylm gradients have to be in memory for choice=3, 5, or 23!'
535    MSG_BUG(msg)
536  end if
537 
538 !Init parallelism
539  if (paral_kgb==1) then
540    me_distrb=mpi_enreg%me_kpt
541    spaceComm=mpi_enreg%comm_kpt
542    spaceComm_fft=mpi_enreg%comm_fft
543    npband_bandfft=mpi_enreg%nproc_band
544    cg_bandpp=mpi_enreg%bandpp
545    cprj_bandpp=mpi_enreg%bandpp
546    spaceComm_band=mpi_enreg%comm_band
547    cg_band_distributed=.true.
548    cprj_band_distributed=(mcprj/=mcg/mpw)
549  else
550    me_distrb=mpi_enreg%me_kpt
551    spaceComm=mpi_enreg%comm_cell
552    spaceComm_fft=xmpi_comm_self
553    npband_bandfft=1
554    cg_bandpp=1;cprj_bandpp=1
555    cg_band_distributed=.false.
556    cprj_band_distributed=.false.
557    if (mpi_enreg%paralbd==1) then
558      spaceComm_band=mpi_enreg%comm_band
559    else
560      spaceComm_band=xmpi_comm_self
561    end if
562  end if
563  if (cg_bandpp/=cprj_bandpp) then
564    MSG_BUG('cg_bandpp must be equal to cprj_bandpp !')
565  end if
566 
567 !Manage parallelization over spinors
568  my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
569 !Check sizes for cprj (distribution is tricky)
570  one_atom=(iatom>0)
571  if (one_atom.and.ncprj/=1) then
572    MSG_BUG('Bad value for ncprj dimension (should be 1) !')
573  end if
574  if (.not.one_atom.and.ncprj/=natom) then
575    MSG_BUG('Bad value for ncprj dimension (should be natom) !')
576  end if
577 
578 !Initialize some variables
579  mband_cprj=mcprj/(my_nspinor*mkmem*nsppol)
580  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
581  n1_2p1=2*n1+1;n2_2p1=2*n2+1;n3_2p1=2*n3+1
582  ibg=0;icg=0;cpopt=0
583  ider=0;idir0=0;istart=idir;iend=idir
584  if (choice==3.or.choice==5.or.choice==23) ider=1
585  if (idir>0) then
586    if (choice==3) idir0=-idir
587    if (choice==5) idir0=idir
588  else
589 !   if (choice==23) idir0=-7
590    if (choice==3) idir0=-7
591    if (choice==5) idir0=4
592  end if
593  if (idir0==0.or.idir0==4) then
594    dimffnl=1+3*ider
595  else if (idir0/=-7) then
596    dimffnl=1+ider
597  else
598    dimffnl=1+6*ider
599    if(choice==3)then
600      istart=ider;iend=6*ider
601    end if
602  end if
603  nkpg=0
604  if (choice==3.or.choice==2.or.choice==23) nkpg=3*nloalg(3)
605  if (choice==4.or.choice==24) nkpg=9*nloalg(3)
606 
607 !Set number of gradients for <p_i|Cnk>
608  ncpgr=0
609  if (idir==0) then
610    if (choice==2) ncpgr=3
611    if (choice==3) ncpgr=6
612    if (choice==23)ncpgr=9
613    if (choice==4) ncpgr=6
614    if (choice==24)ncpgr=9
615    if (choice==5) ncpgr=3
616    if (choice==6) ncpgr=63
617  else
618    ncpgr=1
619  end if
620 !Test cprj gradients dimension (just to be sure)
621  if (cprj(1,1)%ncpgr/=ncpgr) then
622    MSG_BUG('cprj are badly allocated !')
623  end if
624 
625 
626 !Extract data for treated atom(s)
627  if (one_atom) then
628    iatom1=iatom;iatom2=iatom
629    ntypat0=1;itypat=typat(iatom)
630    ABI_ALLOCATE(nattyp_atm,(ntypat0))
631    nattyp_atm(1)=1
632    ABI_ALLOCATE(atindx_atm,(ntypat0))
633    atindx_atm(1)=atindx(iatom)
634    ABI_ALLOCATE(ph1d_atm,(2,(n1_2p1+n2_2p1+n3_2p1)*ntypat0))
635    shift1=(atindx(iatom)-1)*n1_2p1
636    shift2=(atindx(iatom)-1)*n2_2p1+natom*n1_2p1
637    shift3=(atindx(iatom)-1)*n3_2p1+natom*(n1_2p1+n2_2p1)
638    shift1b=0;shift2b=n1_2p1;shift3b=n1_2p1+n2_2p1
639    ph1d_atm(:,shift1b+1:shift1b+n1_2p1)=ph1d(:,shift1+1:shift1+n1_2p1)
640    ph1d_atm(:,shift2b+1:shift2b+n2_2p1)=ph1d(:,shift2+1:shift2+n2_2p1)
641    ph1d_atm(:,shift3b+1:shift3b+n3_2p1)=ph1d(:,shift3+1:shift3+n3_2p1)
642    ABI_ALLOCATE(ekb_atm,(psps%dimekb,ntypat0))
643    ABI_ALLOCATE(indlmn_atm,(6,psps%lmnmax,ntypat0))
644    ABI_ALLOCATE(ffspl_atm,(psps%mqgrid_ff,2,psps%lnmax,ntypat0))
645    ABI_ALLOCATE(pspso_atm,(ntypat0))
646    ekb_atm(:,1)=psps%ekb(:,itypat)
647    indlmn_atm(:,:,1)=psps%indlmn(:,:,itypat)
648    ffspl_atm(:,:,:,1)=psps%ffspl(:,:,:,itypat)
649    pspso_atm(1)=psps%pspso(itypat)
650  else
651    iatom1=1;iatom2=natom
652    ntypat0=ntypat
653    atindx_atm => atindx
654    nattyp_atm => nattyp
655    ph1d_atm => ph1d
656    ekb_atm => psps%ekb
657    indlmn_atm => psps%indlmn
658    ffspl_atm => psps%ffspl
659    pspso_atm => psps%pspso
660  end if
661 
662 !Dimensioning and allocation of <p_i|Cnk>
663  ABI_ALLOCATE(dimlmn,(ncprj))
664  dimlmn=0  ! Type-sorted cprj
665  if (one_atom) then
666    itypat=typat(iatom)
667    dimlmn(ia+1:ia+nattyp(itypat))=count(indlmn_atm(3,:,itypat)>0)
668  else
669    ia=0
670    do itypat=1,ntypat0
671      dimlmn(ia+1:ia+nattyp(itypat))=count(indlmn_atm(3,:,itypat)>0)
672      ia=ia+nattyp(itypat)
673    end do
674  end if
675  ABI_DATATYPE_ALLOCATE(cwaveprj,(ncprj,my_nspinor*cprj_bandpp))
676  call pawcprj_alloc(cwaveprj,ncpgr,dimlmn)
677 
678 !Additional statements if band-fft parallelism
679  if (npband_bandfft>1) then
680    ABI_ALLOCATE(npw_block,(npband_bandfft))
681    ABI_ALLOCATE(npw_disp,(npband_bandfft))
682    ABI_ALLOCATE(bufsize,(npband_bandfft*cg_bandpp))
683    ABI_ALLOCATE(bufdisp,(npband_bandfft*cg_bandpp))
684    ABI_ALLOCATE(bufsize_wf,(npband_bandfft*cg_bandpp))
685    ABI_ALLOCATE(bufdisp_wf,(npband_bandfft*cg_bandpp))
686  end if
687 
688 !Set output datastructure to zero
689  call pawcprj_set_zero(cprj)
690 
691 !LOOP OVER SPINS
692  do isppol=1,nsppol
693    ikg=0
694 
695 !  BIG FAT k POINT LOOP
696    do ikpt=1,nkpt
697      counter=100*ikpt+isppol
698 
699 !    Select k point to be treated by this proc
700      nband_k=nband(ikpt+(isppol-1)*nkpt)
701      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
702 
703 !    Retrieve k-point
704      kpoint(:)=kpt(:,ikpt)
705      istwf_k=istwfk(ikpt)
706 
707 !    Retrieve number of plane waves
708      npw_k=npwarr(ikpt)
709      if (npband_bandfft>1) then
710 !      Special treatment for band-fft //
711        call xmpi_allgather(npw_k,npw_block,spaceComm_band,ierr)
712        npw_nk=sum(npw_block);npw_disp(1)=0
713        do ii=2,npband_bandfft
714          npw_disp(ii)=npw_disp(ii-1)+npw_block(ii-1)
715        end do
716      else
717        npw_nk=npw_k
718      end if
719 
720 !    Retrieve (k+G) points and spherical harmonics
721      ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang))
722      ABI_ALLOCATE(ylmgr_k,(npw_k,3,mpsang*mpsang*useylmgr))
723      ABI_ALLOCATE(kg_k,(3,npw_nk))
724      if (npband_bandfft>1) then
725 !      Special treatment for band-fft //
726        ABI_ALLOCATE(kg_k_loc,(3,npw_k))
727        kg_k_loc(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
728        bufsize(:)=3*npw_block(:);bufdisp(:)=3*npw_disp(:)
729        call xmpi_allgatherv(kg_k_loc,3*npw_k,kg_k,bufsize,bufdisp,spaceComm_band,ierr)
730      else
731        kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
732      end if
733      do ilm=1,mpsang*mpsang
734        ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
735        if (useylmgr>0) ylmgr_k(1:npw_k,1:3,ilm)=ylmgr(1+ikg:npw_k+ikg,1:3,ilm)
736      end do
737 
738 !    Compute (k+G) vectors
739      ABI_ALLOCATE(kpg_k,(npw_nk,nkpg))
740      if (nkpg>0) then
741        call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_nk)
742      end if
743 !    Allocate and compute the arrays phkxred and ph3d
744      ABI_ALLOCATE(phkxred,(2,ncprj))
745      do ia=iatom1,iatom2
746        iatm=min(atindx_atm(ia),ncprj)
747        arg=two_pi*(kpoint(1)*xred(1,ia)+kpoint(2)*xred(2,ia)+kpoint(3)*xred(3,ia))
748        phkxred(1,iatm)=cos(arg);phkxred(2,iatm)=sin(arg)
749      end do
750      matblk=ncprj;if (nloalg(2)<=0) matblk=0
751      ABI_ALLOCATE(ph3d,(2,npw_nk,matblk))
752      if (matblk>0)then
753 !      Here, precomputation of ph3d
754        if (npband_bandfft>1) then
755 !        Special treatment for band-fft //
756          ABI_ALLOCATE(ph3d_tmp,(2,npw_k,matblk))
757          call ph1d3d(1,ncprj,kg_k_loc,matblk,ncprj,npw_k,n1,n2,n3,phkxred,ph1d_atm,ph3d_tmp)
758          ABI_ALLOCATE(ph3d_tmp_npw,(2,matblk,npw_k))
759          ABI_ALLOCATE(ph3d_npw,(2,matblk,npw_nk))
760          isize=2*matblk;bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:)
761          do ipw=1,npw_k
762            ph3d_tmp_npw(:,:,ipw)=ph3d_tmp(:,ipw,:)
763          end do
764          call xmpi_allgatherv(ph3d_tmp_npw,isize*npw_k,ph3d_npw,bufsize,bufdisp,spaceComm_band,ierr)
765          do ipw=1,npw_nk
766            ph3d(:,ipw,:)=ph3d_npw(:,:,ipw)
767          end do
768          ABI_DEALLOCATE(ph3d_npw)
769          ABI_DEALLOCATE(ph3d_tmp_npw)
770          ABI_DEALLOCATE(ph3d_tmp)
771        else
772          call ph1d3d(1,ncprj,kg_k,matblk,ncprj,npw_k,n1,n2,n3,phkxred,ph1d_atm,ph3d)
773        end if
774      else if (npband_bandfft>1) then
775        MSG_ERROR('Band-fft parallelism +nloag(1)<0 forbidden !')
776      end if
777 
778 !    Compute nonlocal form factors ffnl at all (k+G)
779      ABI_ALLOCATE(ffnl,(npw_nk,dimffnl,psps%lmnmax,ntypat0))
780      if (npband_bandfft>1) then
781 !      Special treatment for band-fft //
782        ABI_ALLOCATE(ffnl_tmp,(npw_k,dimffnl,psps%lmnmax,ntypat0))
783        call mkffnl(psps%dimekb,dimffnl,ekb_atm,ffnl_tmp,ffspl_atm,&
784 &       gmet,gprimd,ider,idir0,indlmn_atm,kg_k_loc,kpg_k,kpoint,psps%lmnmax,&
785 &       psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat0,&
786 &       pspso_atm,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
787        ABI_ALLOCATE(ffnl_tmp_npw,(dimffnl,psps%lmnmax,ntypat0,npw_k))
788        ABI_ALLOCATE(ffnl_npw,(dimffnl,psps%lmnmax,ntypat0,npw_nk))
789        isize=dimffnl*psps%lmnmax*ntypat0
790        bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:)
791        do ipw=1,npw_k
792          ffnl_tmp_npw(:,:,:,ipw)=ffnl_tmp(ipw,:,:,:)
793        end do
794        call xmpi_allgatherv(ffnl_tmp_npw,isize*npw_k,ffnl_npw,bufsize,bufdisp,spaceComm_band,ierr)
795        do ipw=1,npw_nk
796          ffnl(ipw,:,:,:)=ffnl_npw(:,:,:,ipw)
797        end do
798        ABI_DEALLOCATE(ffnl_npw)
799        ABI_DEALLOCATE(ffnl_tmp_npw)
800        ABI_DEALLOCATE(ffnl_tmp)
801      else
802        call mkffnl(psps%dimekb,dimffnl,ekb_atm,ffnl,ffspl_atm,&
803 &       gmet,gprimd,ider,idir0,indlmn_atm,kg_k,kpg_k,kpoint,psps%lmnmax,&
804 &       psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat0,&
805 &       pspso_atm,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
806      end if
807 
808 !    No more need of kg_g_tmp
809      if (npband_bandfft>1)  then
810        ABI_DEALLOCATE(kg_k_loc)
811      end if
812 
813 !    Allocate arrays for a wave-function (or a block of WFs)
814      ABI_ALLOCATE(cwavef,(2,npw_nk*my_nspinor*cg_bandpp))
815      if (npband_bandfft>1) then
816        isize=2*my_nspinor*cg_bandpp;bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:)
817        isize=2*my_nspinor*npw_k*cg_bandpp;bufsize_wf(:)=isize
818        do ii=1,npband_bandfft*cg_bandpp
819          bufdisp_wf(ii)=(ii-1)*isize
820        end do
821      end if
822 
823 !    Loop over bands or blocks of bands
824      icgb=icg ; ibgb=ibg ; iband_start=1
825      blocksz=npband_bandfft*cg_bandpp
826      nblockbd=nband_k/blocksz
827      nband_cprj_k=merge(nband_k/npband_bandfft,nband_k,cprj_band_distributed)
828      do iblockbd=1,nblockbd
829        iband_min=1+(iblockbd-1)*blocksz
830        iband_max=iblockbd*blocksz
831 
832        if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min,iband_max,isppol,me_distrb)) then
833          if (.not.cg_band_distributed) icgb=icgb+npw_k*my_nspinor*blocksz
834          if (.not.cprj_band_distributed) ibgb=ibgb+my_nspinor*blocksz
835          cycle
836        end if
837 
838 !      Extract wavefunction information
839 !      Special treatment for band-fft parallelism
840        if (npband_bandfft>1) then
841          !Transpose WF to get them in "FFT" representation
842          ABI_ALLOCATE(cwavef_tmp,(2,npw_k*my_nspinor*blocksz))
843          do ig=1,npw_k*my_nspinor*blocksz
844            cwavef_tmp(1,ig)=cg(1,ig+icgb)
845            cwavef_tmp(2,ig)=cg(2,ig+icgb)
846          end do
847          call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef,bufsize,bufdisp,spaceComm_band,ierr)
848          ABI_DEALLOCATE(cwavef_tmp)
849          !Reorder WF according to cg_bandpp and/or spinor
850          if (cg_bandpp>1.or.my_nspinor>1) then
851            ABI_ALLOCATE(cwavef_tmp,(2,npw_nk*my_nspinor*blocksz))
852            do ig=1,npw_nk*my_nspinor*blocksz
853              cwavef_tmp(:,ig)=cwavef(:,ig)
854            end do
855            shift1=0
856            do iwf2=1,cg_bandpp
857              do ig=1,my_nspinor
858                shift2=0
859                do iwf1=1,npband_bandfft
860                  npws=npw_block(iwf1)
861                  ipw=shift2+(iwf2-1)*my_nspinor*npws+(ig-1)*npws
862                  cwavef(:,shift1+1:shift1+npws)=cwavef_tmp(:,ipw+1:ipw+npws)
863                  shift1=shift1+npws ; shift2=shift2+cg_bandpp*my_nspinor*npws
864                end do
865              end do
866            end do
867            ABI_DEALLOCATE(cwavef_tmp)
868          end if
869        else
870          do ig=1,npw_k*my_nspinor*cg_bandpp
871            cwavef(1,ig)=cg(1,ig+icgb)
872            cwavef(2,ig)=cg(2,ig+icgb)
873          end do
874        end if
875 
876 !      Compute scalar product of wavefunction with all NL projectors
877        do ibp=1,cg_bandpp   ! Note: we suppose cp_bandpp=cprj_bandpp
878          iwf1=1+(ibp-1)*npw_nk*my_nspinor;iwf2=ibp*npw_nk*my_nspinor
879          icp1=1+(ibp-1)*my_nspinor;icp2=ibp*my_nspinor
880          do jdir=istart,iend
881            call getcprj(choice,cpopt,cwavef(:,iwf1:iwf2),cwaveprj(:,icp1:icp2),&
882 &           ffnl,jdir,indlmn_atm,istwf_k,kg_k,kpg_k,kpoint,psps%lmnmax,&
883 &           mgfft,mpi_enreg,ncprj,nattyp_atm,ngfft,nloalg,&
884 &           npw_nk,my_nspinor,ntypat0,phkxred,ph1d_atm,ph3d,ucvol,psps%useylm)
885          end do
886        end do
887 !      Export cwaveprj to big array cprj
888        call pawcprj_put(atindx_atm,cwaveprj,cprj,ncprj,iband_start,ibgb,ikpt,iorder_cprj,isppol,&
889 &       mband_cprj,mkmem,natom,cprj_bandpp,nband_cprj_k,dimlmn,my_nspinor,nsppol,uncp,&
890 &       mpi_comm_band=spaceComm_band,to_be_gathered=(cg_band_distributed.and.(.not.cprj_band_distributed)))
891 
892        iband_start=iband_start+merge(cg_bandpp,blocksz,cprj_band_distributed)
893 
894 !      End loop over bands
895        icgb=icgb+npw_k*my_nspinor*blocksz
896      end do
897 
898 !    Shift array memory (if mkmem/=0)
899      if (mkmem/=0) then
900        ibg=ibg+my_nspinor*nband_cprj_k
901        icg=icg+my_nspinor*nband_k*npw_k
902        ikg=ikg+npw_k
903      end if
904 
905 !    End big k point loop
906      ABI_DEALLOCATE(ffnl)
907      ABI_DEALLOCATE(ph3d)
908      ABI_DEALLOCATE(phkxred)
909      ABI_DEALLOCATE(kg_k)
910      ABI_DEALLOCATE(kpg_k)
911      ABI_DEALLOCATE(ylm_k)
912      ABI_DEALLOCATE(ylmgr_k)
913      ABI_DEALLOCATE(cwavef)
914    end do
915 !  End loop over spins
916  end do
917 
918  if ((iatom<=0).and.(choice==23)) then
919    do iatom1=1,ncprj
920      do ii=1,mcprj
921        nlmn=cprj(iatom1,ii)%nlmn
922        do kk=1,nlmn
923          work(1:6)=cprj(iatom1,ii)%dcp(1,1:6,kk)
924          call strconv(work,gprimd,work)
925          cprj(iatom1,ii)%dcp(1,1:6,kk)=work(1:6)
926          work(1:6)=cprj(iatom1,ii)%dcp(2,1:6,kk)
927          call strconv(work,gprimd,work)
928          cprj(iatom1,ii)%dcp(2,1:6,kk)=work(1:6)
929        end do
930      end do
931    end do
932  end if
933 
934 !If needed, gather computed scalars
935  if (.not.cg_band_distributed) then
936    call pawcprj_mpi_sum(cprj,spaceComm_band,ierr)
937  end if
938 
939 !Deallocate temporary storage
940  if (one_atom)  then
941    ABI_DEALLOCATE(atindx_atm)
942    ABI_DEALLOCATE(nattyp_atm)
943    ABI_DEALLOCATE(ph1d_atm)
944    ABI_DEALLOCATE(ekb_atm)
945    ABI_DEALLOCATE(indlmn_atm)
946    ABI_DEALLOCATE(ffspl_atm)
947    ABI_DEALLOCATE(pspso_atm)
948  end if
949  nullify(atindx_atm,nattyp_atm,ph1d_atm,ekb_atm,indlmn_atm,ffspl_atm,pspso_atm)
950  call pawcprj_free(cwaveprj)
951  ABI_DATATYPE_DEALLOCATE(cwaveprj)
952  ABI_DEALLOCATE(dimlmn)
953  if (npband_bandfft>1) then
954    ABI_DEALLOCATE(npw_block)
955    ABI_DEALLOCATE(npw_disp)
956    ABI_DEALLOCATE(bufsize)
957    ABI_DEALLOCATE(bufdisp)
958    ABI_DEALLOCATE(bufsize_wf)
959    ABI_DEALLOCATE(bufdisp_wf)
960  end if
961 
962  DBG_EXIT('COLL')
963 
964  end subroutine ctocprj

ABINIT/getcprj [ Functions ]

[ Top ] [ Functions ]

NAME

 getcprj

FUNCTION

  Compute <Proj_i|Cnk> for one wave function |Cnk> expressed in reciprocal space.
  Compute also derivatives of <Proj_i|Cnk>.
  |Proj_i> are non-local projectors (for each atom and each l,m,n)

INPUTS

  choice=chooses possible output:
    In addition to projected wave function:
    choice=1 => nothing else
          =2 => 1st gradients with respect to atomic position(s)
          =3 => 1st gradients with respect to strain(s)
          =23=> 1st gradients with respect to strain(s) and atm pos
          =4 => 2nd derivatives with respect to atomic pos.
          =24=> 1st and 2nd derivatives with respect to atomic pos.
          =5 => 1st gradients with respect to k wavevector
          =6 => 2nd derivatives with respect to strain and atm. pos.
  cpopt=1 if <Proj_i|Cnk> are already in memory; see below (side effects).
  cwavef(2,nspinor*npw_k)=input cmplx wavefunction coefficients <G|Cnk>
  ffnl(npw_k,dimffnl,lmnmax,ntypat)=nonlocal form factors to be used for the application of the nl operator
  idir=direction of the derivative, i.e. dir. of - atom to be moved  in the case choice=2
                                                 - strain component  in the case choice=3
                                                 - k point direction in the case choice=5
       Compatible only with choice=2,3,5; if idir=0, all derivatives are computed
  indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=lmn
  istwf_k=option parameter that describes the storage of wfs
  kg_k(3,npw_k)=reduced planewave coordinates
  kpg(npw_k,npk)=(k+G) components and related data
  kpoint(3)=k point in terms of recip. translations
  lmnmax=max. number of (l,m,n) components over all types of atoms
  mgfft=maximum size of 1D FFTs
  mpi_enreg=informations about MPI parallelization
  natom=number of atoms in cell
  nattyp(ntypat)=number of atoms of each type
  ngfft(18)=contain all needed information about 3D FFT, see ~ABINIT/Infos/vargs.htm#ngfft
  nloalg(3)=governs the choice of the algorithm for nonlocal operator
  npw_k=number of planewaves for given k point
  nspinor=number of spinorial components of the wavefunctions (on current proc)
  ntypat=number of types of atoms in unit cell
  phkxred(2,natom)=phase factors exp(2 pi kpoint.xred)
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
  ph3d(2,npw_k,natom)=3D structure factors, for each atom and plane wave
                      only used if nloalg(2)>0
  ucvol= unit cell volume
  useylm=governs the way the nonlocal operator is to be applied

SIDE EFFECTS

  cwaveprj(natom,nspinor) <type(pawcprj_type)>=projected input wave function <Proj_i|Cnk> with all NL projectors
                                (and derivatives)
                                if cpopt=1 the projected scalars have been already been computed and
                                           only derivatives are computed here
                                if cpopt=0 the projected scalars and derivatives are computed here

TODO

  Spin-orbit

PARENTS

      cgwf,ctocprj,debug_tools,dfpt_accrho,dfpt_nstpaw,ks_ddiago,m_wfd
      rf2_init

CHILDREN

      mkkpg,opernla_ylm,ph1d3d

SOURCE

124  subroutine getcprj(choice,cpopt,cwavef,cwaveprj,ffnl,&
125 &                   idir,indlmn,istwf_k,kg_k,kpg,kpoint,lmnmax,mgfft,mpi_enreg,&
126 &                   natom,nattyp,ngfft,nloalg,npw_k,nspinor,ntypat,&
127 &                   phkxred,ph1d,ph3d,ucvol,useylm)
128 
129 
130 !This section has been created automatically by the script Abilint (TD).
131 !Do not modify the following lines by hand.
132 #undef ABI_FUNC
133 #define ABI_FUNC 'getcprj'
134 !End of the abilint section
135 
136  implicit none
137 
138 !Arguments -------------------------------
139 !scalars
140  integer,intent(in) :: choice,cpopt,idir,istwf_k,lmnmax
141  integer,intent(in) :: mgfft,natom,npw_k,nspinor,ntypat,useylm
142  real(dp),intent(in) :: ucvol
143  type(MPI_type),intent(in) :: mpi_enreg
144 !arrays
145  integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg_k(3,npw_k),nattyp(ntypat)
146  integer,intent(in) :: ngfft(18),nloalg(3)
147  real(dp),intent(in) :: cwavef(2,npw_k*nspinor)
148  real(dp),intent(in),target :: ffnl(:,:,:,:),kpg(:,:),ph3d(:,:,:)
149  real(dp),intent(in) :: kpoint(3),ph1d(2,3*(2*mgfft+1)*natom),phkxred(2,natom)
150  type(pawcprj_type),intent(inout) :: cwaveprj(natom,nspinor)
151 
152 !Local variables-------------------------------
153 !scalars
154  integer :: choice_,cplex,dimffnl,ia,ia1,ia2,ia3,ia4,iatm,ic,ii,ilmn,ishift,ispinor,itypat
155  integer :: jc,matblk,mincat,nd2gxdt,ndgxdt,nincat,nkpg,nkpg_,nlmn,signs
156 !arrays
157  integer,allocatable :: cplex_dgxdt(:),cplex_d2gxdt(:),indlmn_typ(:,:)
158  real(dp),allocatable :: d2gxdt(:,:,:,:,:),dgxdt(:,:,:,:,:),ffnl_typ(:,:,:)
159  real(dp),allocatable :: gx(:,:,:,:)
160  real(dp), pointer :: kpg_(:,:),ph3d_(:,:,:)
161 
162 ! *********************************************************************
163 
164  DBG_ENTER('COLL')
165 
166 !Nothing to do in that case
167  if (cpopt==1.and.choice==1) return
168 
169 !Not available for useylm=0
170  if (useylm==0) then
171    MSG_ERROR('Not available for useylm=0 !')
172  end if
173 
174 !Error on bad choice
175  if ((choice<1.or.choice>6).and.choice/=23.and.choice/=24) then
176    MSG_BUG('Does not presently support this choice !')
177  end if
178 
179 !Error on bad idir
180  if (idir>0.and.choice/=2.and.choice/=3.and.choice/=5) then
181    MSG_BUG('Does not support idir>0 for this choice')
182  end if
183 
184 !Error on sizes
185  nkpg=size(kpg,2)
186  if (nkpg>0) then
187    if( (choice==2.and.nkpg<3) .or. &
188 &   ((choice==4.or.choice==24).and.nkpg<9) .or. &
189 &   ((choice==6.or.choice==3.or.choice==23).and.nkpg<3) ) then
190      MSG_BUG('Incorrect size for kpg array !')
191    end if
192  end if
193  if (size(ffnl,1)/=npw_k.or.size(ffnl,3)/=lmnmax) then
194    MSG_BUG('Incorrect size for ffnl!')
195  end if
196  if (size(ph3d)>0) then
197    if (size(ph3d,2)/=npw_k) then
198      MSG_BUG('Incorrect size for ph3d!')
199    end if
200  end if
201 
202 !Define dimensions of projected scalars
203  dimffnl=size(ffnl,2)
204  ndgxdt=0;nd2gxdt=0
205  if (idir==0) then
206    if (choice==2) ndgxdt=3
207    if (choice==3) ndgxdt=6
208    if (choice==23) ndgxdt=9
209    if (choice==4) nd2gxdt=6
210    if (choice==24) then
211      ndgxdt=3;nd2gxdt=6
212    end if
213    if (choice==5) ndgxdt=3
214    if (choice==6) then
215      ndgxdt=9;nd2gxdt=54
216    end if
217  else
218    ndgxdt=1
219  end if
220  if(cwaveprj(1,1)%ncpgr<ndgxdt+nd2gxdt) then
221    MSG_BUG('Incorrect size for ncpgr')
222  end if
223 !Eventually re-compute (k+G) vectors (and related data)
224  if (nkpg==0) then
225    nkpg_=0
226    if (choice==4.or.choice==24) nkpg_=9
227    if (choice==2.or.choice==3.or.choice==23) nkpg_=3
228    ABI_ALLOCATE(kpg_,(npw_k,nkpg_))
229    if (nkpg_>0) then
230      call mkkpg(kg_k,kpg_,kpoint,nkpg_,npw_k)
231    end if
232  else
233    nkpg_=nkpg
234    kpg_ => kpg
235  end if
236 
237 !Some other dims
238  mincat=min(NLO_MINCAT,maxval(nattyp))
239  cplex=2;if (istwf_k>1) cplex=1
240  choice_=choice;if (cpopt==1) choice_=-choice
241  signs=1;if (idir>0) signs=2
242 !Eventually allocate temporary array for ph3d
243  if (nloalg(2)<=0) then
244    matblk=mincat
245    ABI_ALLOCATE(ph3d_,(2,npw_k,matblk))
246  else
247    matblk=size(ph3d,3)
248    ph3d_ => ph3d
249  end if
250 
251 !Loop over atom types
252  ia1=1;iatm=0
253  do itypat=1,ntypat
254    ia2=ia1+nattyp(itypat)-1;if (ia2<ia1) cycle
255    nlmn=count(indlmn(3,:,itypat)>0)
256 
257 !  Retrieve some data for this type of atom
258    ABI_ALLOCATE(indlmn_typ,(6,nlmn))
259    ABI_ALLOCATE(ffnl_typ,(npw_k,dimffnl,nlmn))
260    indlmn_typ(:,1:nlmn)=indlmn(:,1:nlmn,itypat)
261    ffnl_typ(:,:,1:nlmn)=ffnl(:,:,1:nlmn,itypat)
262 
263 !  Loop on blocks of atoms inside type
264    do ia3=ia1,ia2,mincat
265      ia4=min(ia2,ia3+mincat-1);nincat=ia4-ia3+1
266 !     Prepare the phase factors if they were not already computed
267      if (nloalg(2)<=0) then
268        call ph1d3d(ia3,ia4,kg_k,matblk,natom,npw_k,ngfft(1),ngfft(2),ngfft(3),&
269 &       phkxred,ph1d,ph3d_)
270      end if
271 
272 !    Allocate memory for projected scalars
273      ABI_ALLOCATE(gx,(cplex,nlmn,nincat,nspinor))
274      ABI_ALLOCATE(dgxdt,(cplex,ndgxdt,nlmn,nincat,nspinor))
275      ABI_ALLOCATE(d2gxdt,(cplex,nd2gxdt,nlmn,nincat,nspinor))
276      ABI_ALLOCATE(cplex_dgxdt,(ndgxdt))
277      ABI_ALLOCATE(cplex_d2gxdt,(nd2gxdt))
278 
279 !    Retrieve eventually <p_i|c> coeffs
280      if (cpopt==1) then
281        do ispinor=1,nspinor
282          do ia=1,nincat
283            gx(1:cplex,1:nlmn,ia,ispinor)=cwaveprj(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)
284          end do
285        end do
286      end if
287 
288 !    Compute <p_i|c> scalars (and derivatives) for this block of atoms
289      call opernla_ylm(choice_,cplex,cplex_dgxdt,cplex_d2gxdt,dimffnl,d2gxdt,dgxdt,ffnl_typ,gx,&
290 &     ia3,idir,indlmn_typ,istwf_k,kpg_,matblk,mpi_enreg,nd2gxdt,ndgxdt,nincat,nkpg_,nlmn,&
291 &     nloalg,npw_k,nspinor,ph3d_,signs,ucvol,cwavef)
292 
293 !    Transfer result to output variable cwaveprj
294      if (cpopt==0) then
295        do ispinor=1,nspinor
296          do ia=1,nincat
297            cwaveprj(iatm+ia,ispinor)%nlmn=nlmn
298            cwaveprj(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)=gx(1:cplex,1:nlmn,ia,ispinor)
299            if(cplex==1) cwaveprj(iatm+ia,ispinor)%cp(2,1:nlmn)=zero
300          end do
301        end do
302      end if
303      if (cpopt>=0.and.choice>1) then
304        ishift=0
305        if ((idir>0).and.(cwaveprj(1,1)%ncpgr>ndgxdt)) ishift=idir-1
306        if(cplex==2)then
307          do ispinor=1,nspinor
308            do ia=1,nincat
309 !             cwaveprj(iatm+ia,ispinor)%ncpgr=ndgxdt+nd2gxdt
310              if (ndgxdt>0) cwaveprj(iatm+ia,ispinor)%dcp(1:2,1+ishift:ndgxdt+ishift,1:nlmn)=&
311 &             dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor)
312              if (nd2gxdt>0)cwaveprj(iatm+ia,ispinor)%dcp(1:2,ndgxdt+1+ishift:ndgxdt+nd2gxdt+ishift,1:nlmn)=&
313 &             d2gxdt(1:2,1:nd2gxdt,1:nlmn,ia,ispinor)
314            end do
315          end do
316        else
317 !        cplex_dgxdt(i)  = 1 if dgxdt(1,i,:,:)  is real, 2 if it is pure imaginary
318 !        cplex_d2gxdt(i) = 1 if d2gxdt(1,i,:,:) is real, 2 if it is pure imaginary
319          do ispinor=1,nspinor
320            do ia=1,nincat
321 !             cwaveprj(iatm+ia,ispinor)%ncpgr=ndgxdt+nd2gxdt
322              if (ndgxdt>0) then
323                do ilmn =1,nlmn
324                  do ii = 1,ndgxdt
325                    ic = cplex_dgxdt(ii) ; jc = 3 - ic
326                    cwaveprj(iatm+ia,ispinor)%dcp(ic,ii+ishift,ilmn)=dgxdt(1,ii,ilmn,ia,ispinor)
327                    cwaveprj(iatm+ia,ispinor)%dcp(jc,ii+ishift,ilmn)=zero
328                  end do
329                end do
330              end if
331              if (nd2gxdt>0) then
332                do ilmn =1,nlmn
333                  do ii = 1,nd2gxdt
334                    ic = cplex_d2gxdt(ii) ; jc = 3 - ic
335                    cwaveprj(iatm+ia,ispinor)%dcp(ic,ndgxdt+ii+ishift,ilmn)=d2gxdt(1,ii,ilmn,ia,ispinor)
336                    cwaveprj(iatm+ia,ispinor)%dcp(jc,ndgxdt+ii+ishift,ilmn)=zero
337                  end do
338                end do
339              end if
340            end do
341          end do
342        end if
343      end if
344 
345 !    End loop inside block of atoms
346      iatm=iatm+nincat
347      ABI_DEALLOCATE(gx)
348      ABI_DEALLOCATE(dgxdt)
349      ABI_DEALLOCATE(d2gxdt)
350      ABI_DEALLOCATE(cplex_dgxdt)
351      ABI_DEALLOCATE(cplex_d2gxdt)
352    end do
353 
354 !  End loop over atom types
355    ia1=ia2+1
356    ABI_DEALLOCATE(indlmn_typ)
357    ABI_DEALLOCATE(ffnl_typ)
358  end do
359 
360  if (nkpg==0) then
361    ABI_DEALLOCATE(kpg_)
362  end if
363  if (nloalg(2)<=0) then
364    ABI_DEALLOCATE(ph3d_)
365  end if
366 
367  DBG_EXIT('COLL')
368 
369  end subroutine getcprj

ABINIT/m_cgprj [ Modules ]

[ Top ] [ Modules ]

NAME

  m_cgprj

FUNCTION

   Routines to compute <Proj_i|Cnk> with |Cnk> expressed in reciprocal space.

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_cgprj
28 
29  use defs_basis
30  use m_abicore
31  use m_errors
32  use m_xmpi
33 
34  use defs_abitypes, only : MPI_type
35  use defs_datatypes, only : pseudopotential_type
36  use m_kg,       only : ph1d3d, mkkpg
37  use m_geometry, only : strconv
38  use m_mkffnl,   only : mkffnl
39  use m_mpinfo,   only : proc_distrb_cycle
40  use m_pawcprj,  only : pawcprj_type, pawcprj_alloc, pawcprj_put, pawcprj_free, &
41                         pawcprj_set_zero, pawcprj_mpi_sum
42  use m_opernla_ylm, only : opernla_ylm
43 
44  implicit none
45 
46  private