TABLE OF CONTENTS


ABINIT/m_paw_optics [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_optics

FUNCTION

  This module contains several routines related to conductivity:
    optical conductivity, X spectroscopy, linear susceptibility, ...

COPYRIGHT

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

SOURCE

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 MODULE m_paw_optics
25 
26  use defs_basis
27  use defs_abitypes
28  use m_xmpi
29  use m_errors
30  use m_wffile
31  use m_abicore
32  use m_hdr
33 
34  use m_time,         only : timab
35  use m_io_tools,     only : open_file,get_unit
36  use m_pawpsp,       only : pawpsp_read_corewf
37  use m_pawrad,       only : pawrad_type, pawrad_deducer0, simp_gen
38  use m_pawtab,       only : pawtab_type
39  use m_pawcprj,      only : pawcprj_type, pawcprj_alloc, pawcprj_get, &
40 &                           pawcprj_free,pawcprj_mpi_allgather
41  use m_paw_onsite,   only : pawnabla_init,pawnabla_core_init
42  use m_mpinfo,       only : destroy_mpi_enreg,nullify_mpi_enreg,initmpi_seq,proc_distrb_cycle
43  use m_numeric_tools,only : kramerskronig
44  use m_geometry,     only : metric
45  use m_hide_lapack,  only : matrginv
46 
47  implicit none
48 
49  private
50 
51 !public procedures.
52  public :: optics_paw
53  public :: optics_paw_core
54  public :: linear_optics_paw
55 
56 CONTAINS  !========================================================================================

m_paw_optics/linear_optics_paw [ Functions ]

[ Top ] [ m_paw_optics ] [ Functions ]

NAME

 linear_optics_paw

FUNCTION

 This program computes the elements of the optical frequency dependent
 linear susceptiblity using matrix elements <-i Nabla> obtained from a
 PAW ground state calculation. It uses formula 17 from Gadoc et al,
 Phys. Rev. B 73, 045112 (2006) [[cite:Gajdo2006]] together with a scissors correction. It uses
 a Kramers-Kronig transform to compute the real part from the imaginary part, and
 it will work on all types of unit cells. It outputs all tensor elements of
 both the real and imaginary parts.

INPUTS

  filnam: base of file names to read data from
  mpi_enreg: mpi set up variable, not used in this code

OUTPUT

  _real and _imag output files

NOTES

  This routine is not tested

PARENTS

      conducti

CHILDREN

      destroy_mpi_enreg,hdr_free,hdr_io,hdr_read_from_fname,initmpi_seq
      kramerskronig,matrginv,metric,wffopen

SOURCE

 856  subroutine linear_optics_paw(filnam,filnam_out,mpi_enreg_seq)
 857 
 858 
 859 !This section has been created automatically by the script Abilint (TD).
 860 !Do not modify the following lines by hand.
 861 #undef ABI_FUNC
 862 #define ABI_FUNC 'linear_optics_paw'
 863 !End of the abilint section
 864 
 865  implicit none
 866 
 867 !Arguments -----------------------------------
 868 !scalars
 869  character(len=fnlen),intent(in) :: filnam,filnam_out
 870  type(MPI_type),intent(inout) :: mpi_enreg_seq
 871 
 872 !Local variables-------------------------------
 873  integer,parameter :: master=0
 874  integer :: iomode,bantot,bdtot_index,fform1,headform
 875  integer :: iband,ierr,ii,ikpt,iom,iout,isppol,isym,jband,jj,me,mband
 876  integer :: method,mom,nband_k,nkpt,nspinor,nsppol,nsym,occopt,only_check
 877  integer :: rdwr,spaceComm,inpunt,reunt,imunt,wfunt
 878  integer,allocatable :: nband(:),symrel(:,:,:)
 879  real(dp) :: del,dom,fij,gdelta,omin,omax,paijpbij(2),mbpt_sciss,wij,ucvol
 880  real(dp) :: diffwp, diffwm
 881  real(dp) :: e2rot(3,3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),rprimdinv(3,3),symd(3,3),symdinv(3,3)
 882  real(dp),allocatable :: e1(:,:,:),e2(:,:,:,:),epsilon_tot(:,:,:,:),eigen0(:),eig0_k(:)
 883  real(dp),allocatable :: kpts(:,:),occ(:),occ_k(:),oml1(:),wtk(:)
 884  complex,allocatable :: eps_work(:)
 885  character(len=fnlen) :: filnam1,filnam_gen
 886  character(len=500) :: msg
 887  type(hdr_type) :: hdr
 888  type(wffile_type) :: wff1
 889 !arrays
 890  real(dp),allocatable :: psinablapsi(:,:,:,:)
 891 
 892 ! *********************************************************************************
 893 
 894  DBG_ENTER("COLL")
 895 
 896 !* Fake MPI_type for the sequential part.
 897 !This routine should not be parallelized as communicating gbig and other
 898 !tables takes more time than recalculating them in sequential.
 899  call initmpi_seq(MPI_enreg_seq)
 900 
 901 !write(std_out,'(a)')' Give the name of the output file ...'
 902 !read(std_in, '(a)') filnam_out
 903 !write(std_out,'(a)')' The name of the output file is :',filnam_out
 904 
 905 !Read data file
 906  if (open_file(filnam,msg,newunit=inpunt,form='formatted') /= 0 ) then
 907    MSG_ERROR(msg)
 908  end if
 909 
 910  rewind(inpunt)
 911  read(inpunt,*)
 912  read(inpunt,'(a)')filnam_gen       ! generic name for the files
 913  filnam1=trim(filnam_gen)//'_OPT' ! nabla matrix elements file
 914 
 915 !Open the Wavefunction and optic files
 916 !These default values are typical of sequential use
 917  iomode=IO_MODE_FORTRAN ; spaceComm=xmpi_comm_self; me=0
 918 
 919 ! Read the header of the optic files
 920  call hdr_read_from_fname(hdr, filnam1, fform1, spaceComm)
 921  call hdr_free(hdr)
 922  if (fform1 /= 610) then
 923    MSG_ERROR("Abinit8 requires an OPT file with fform = 610")
 924  end if
 925 
 926 !Open the conducti optic files
 927  wfunt = get_unit()
 928  call WffOpen(iomode,spaceComm,filnam1,ierr,wff1,master,me,wfunt)
 929 
 930 !Read the header from Ground state file
 931  rdwr=1
 932  call hdr_io(fform1,hdr,rdwr,wff1)
 933 
 934 !Extract info from the header
 935  headform=hdr%headform
 936  bantot=hdr%bantot
 937  nkpt=hdr%nkpt
 938  ABI_ALLOCATE(kpts,(3,nkpt))
 939  ABI_ALLOCATE(wtk,(nkpt))
 940  kpts(:,:)=hdr%kptns(:,:)
 941  wtk(:)=hdr%wtk(:)
 942  nspinor=hdr%nspinor
 943  nsppol=hdr%nsppol
 944  occopt=hdr%occopt
 945  rprimd(:,:)=hdr%rprimd(:,:)
 946  rprimdinv(:,:) = rprimd(:,:)
 947  call matrginv(rprimdinv,3,3) ! need the inverse of rprimd to symmetrize the tensors
 948  ABI_ALLOCATE(nband,(nkpt*nsppol))
 949  ABI_ALLOCATE(occ,(bantot))
 950  occ(1:bantot)=hdr%occ(1:bantot)
 951  nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
 952  nsym=hdr%nsym
 953  ABI_ALLOCATE(symrel,(3,3,nsym))
 954  symrel(:,:,:)=hdr%symrel(:,:,:)
 955 
 956 !Get mband, as the maximum value of nband(nkpt)
 957  mband=maxval(nband(:))
 958 
 959 !get ucvol etc.
 960  iout = -1
 961  call metric(gmet,gprimd,iout,rmet,rprimd,ucvol)
 962 
 963  write(std_out,*)
 964  write(std_out,'(a,3f10.5,a)' )' rprimd(bohr)      =',rprimd(1:3,1)
 965  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(1:3,2)
 966  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(1:3,3)
 967  write(std_out,*)
 968  write(std_out,'(a,3f10.5,a)' )' rprimdinv         =',rprimdinv(1:3,1)
 969  write(std_out,'(a,3f10.5,a)' )'                    ',rprimdinv(1:3,2)
 970  write(std_out,'(a,3f10.5,a)' )'                    ',rprimdinv(1:3,3)
 971  write(std_out,'(a,2i8)')      ' nkpt,mband        =',nkpt,mband
 972 
 973 !get eigen0
 974  ABI_ALLOCATE(eigen0,(mband*nkpt*nsppol))
 975  read(wfunt)(eigen0(iband),iband=1,mband*nkpt*nsppol)
 976 
 977  read(inpunt,*)mbpt_sciss
 978  read(inpunt,*)dom,omin,omax,mom
 979  close(inpunt)
 980 
 981  ABI_ALLOCATE(oml1,(mom))
 982  ABI_ALLOCATE(e1,(3,3,mom))
 983  ABI_ALLOCATE(e2,(2,3,3,mom))
 984  ABI_ALLOCATE(epsilon_tot,(2,3,3,mom))
 985  ABI_ALLOCATE(eps_work,(mom))
 986  del=(omax-omin)/(mom-1)
 987  do iom=1,mom
 988    oml1(iom)=omin+dble(iom-1)*del
 989  end do
 990  write(std_out,'(a,i8,4f10.5,a)')' npts,omin,omax,width,mbpt_sciss      =',mom,omin,omax,dom,mbpt_sciss,' Ha'
 991 
 992  ABI_ALLOCATE(psinablapsi,(2,3,mband,mband))
 993 
 994 !loop over spin components
 995  do isppol=1,nsppol
 996    bdtot_index = 0
 997 !  loop over k points
 998    do ikpt=1,nkpt
 999 !
1000 !    number of bands for this k point
1001      nband_k=nband(ikpt+(isppol-1)*nkpt)
1002      ABI_ALLOCATE(eig0_k,(nband_k))
1003      ABI_ALLOCATE(occ_k,(nband_k))
1004 !    eigenvalues for this k-point
1005      eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
1006 !    occupation numbers for this k-point
1007      occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
1008 !    values of -i*nabla matrix elements for this k point
1009      psinablapsi=zero
1010      read(wfunt)((psinablapsi(1:2,1,iband,jband),iband=1,nband_k),jband=1,nband_k)
1011      read(wfunt)((psinablapsi(1:2,2,iband,jband),iband=1,nband_k),jband=1,nband_k)
1012      read(wfunt)((psinablapsi(1:2,3,iband,jband),iband=1,nband_k),jband=1,nband_k)
1013 
1014 !    occupation numbers for k-point
1015      occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
1016 !    accumulate e2 for this k point, Eq. 17 from PRB 73, 045112 (2006) [[cite:Gajdo2006]]
1017      do iband = 1, nband_k
1018        do jband = 1, nband_k
1019          fij = occ_k(iband) - occ_k(jband) !occ number difference
1020          wij = eig0_k(iband) - eig0_k(jband) !energy difference
1021          if (abs(fij) > zero) then ! only consider states of differing occupation numbers
1022            do ii = 1, 3
1023              do jj = 1, 3
1024                paijpbij(1) = psinablapsi(1,ii,iband,jband)*psinablapsi(1,jj,iband,jband) + &
1025 &               psinablapsi(2,ii,iband,jband)*psinablapsi(2,jj,iband,jband)
1026                paijpbij(2) = psinablapsi(2,ii,iband,jband)*psinablapsi(1,jj,iband,jband) - &
1027 &               psinablapsi(1,ii,iband,jband)*psinablapsi(2,jj,iband,jband)
1028                do iom = 1, mom
1029 !                original version
1030 !                diffw = wij + mbpt_sciss - oml1(iom) ! apply scissors term here
1031 !                gdelta = exp(-diffw*diffw/(4.0*dom*dom))/(2.0*dom*sqrt(pi)) ! delta fnc resolved as Gaussian
1032 !                e2(1,ii,jj,iom) = e2(1,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(1)*gdelta/(oml1(iom)*oml1(iom))
1033 !                e2(2,ii,jj,iom) = e2(2,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(2)*gdelta/(oml1(iom)*oml1(iom))
1034                  diffwm = wij - mbpt_sciss + oml1(iom) ! apply scissors term here
1035                  diffwp = wij + mbpt_sciss - oml1(iom) ! apply scissors term here
1036                  gdelta = exp(-diffwp*diffwp/(4.0*dom*dom))/(2.0*dom*sqrt(pi))
1037                  e2(1,ii,jj,iom) = e2(1,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(1)*gdelta/(wij*wij)
1038                  e2(2,ii,jj,iom) = e2(2,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(2)*gdelta/(wij*wij)
1039                end do ! end loop over spectral points
1040              end do ! end loop over jj = 1, 3
1041            end do ! end loop over ii = 1, 3
1042          end if ! end selection on fij /= 0
1043        end do ! end loop over jband
1044      end do ! end loop over iband
1045 
1046      ABI_DEALLOCATE(eig0_k)
1047      ABI_DEALLOCATE(occ_k)
1048      bdtot_index=bdtot_index+nband_k
1049    end do ! end loop over k points
1050  end do ! end loop over spin polarizations
1051 
1052 !here apply nsym symrel transformations to reconstruct full tensor from IBZ part
1053  epsilon_tot(:,:,:,:) = zero
1054  do isym = 1, nsym
1055    symd(:,:)=matmul(rprimd(:,:),matmul(symrel(:,:,isym),rprimdinv(:,:)))
1056    symdinv(:,:)=symd(:,:)
1057    call matrginv(symdinv,3,3)
1058    do iom = 1, mom
1059      e2rot(:,:)=matmul(symdinv(:,:),matmul(e2(1,:,:,iom),symd(:,:)))
1060      epsilon_tot(2,:,:,iom) = epsilon_tot(2,:,:,iom)+e2rot(:,:)/nsym
1061    end do
1062  end do
1063 
1064 !generate e1 from e2 via KK transforma
1065  method=0 ! use naive integration ( = 1 for simpson)
1066  only_check=0 ! compute real part of eps in kk routine
1067  do ii = 1, 3
1068    do jj = 1, 3
1069      eps_work(:) = cmplx(0.0,epsilon_tot(2,ii,jj,:))
1070      call kramerskronig(mom,oml1,eps_work,method,only_check)
1071      epsilon_tot(1,ii,jj,:) = real(eps_work(:))
1072      if (ii /= jj) epsilon_tot(1,ii,jj,:) = epsilon_tot(1,ii,jj,:)- 1.0
1073    end do ! end loop over jj
1074  end do ! end loop over ii
1075 
1076  if (open_file(trim(filnam_out)//'_imag',msg,newunit=reunt,form='formatted') /= 0) then
1077    MSG_ERROR(msg)
1078  end if
1079 
1080  if (open_file(trim(filnam_out)//'_real',msg,unit=imunt,form='formatted') /= 0) then
1081    MSG_ERROR(msg)
1082  end if
1083 
1084  write(reunt,'(a12,6a13)')' # Energy/Ha ','eps_2_xx','eps_2_yy','eps_2_zz',&
1085 & 'eps_2_yz','eps_2_xz','eps_2_xy'
1086  write(imunt,'(a12,6a13)')' # Energy/Ha ','eps_1_xx','eps_1_yy','eps_1_zz',&
1087 & 'eps_1_yz','eps_1_xz','eps_1_xy'
1088 
1089  do iom = 1, mom
1090    write(reunt,'(ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4)') oml1(iom),' ',&
1091 &   epsilon_tot(2,1,1,iom),' ',epsilon_tot(2,2,2,iom),' ',epsilon_tot(2,3,3,iom),' ',&
1092 &   epsilon_tot(2,2,3,iom),' ',epsilon_tot(2,1,3,iom),' ',epsilon_tot(2,1,2,iom)
1093    write(imunt,'(ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4)') oml1(iom),' ',&
1094 &   epsilon_tot(1,1,1,iom),' ',epsilon_tot(1,2,2,iom),' ',epsilon_tot(1,3,3,iom),' ',&
1095 &   epsilon_tot(1,2,3,iom),' ',epsilon_tot(1,1,3,iom),' ',epsilon_tot(1,1,2,iom)
1096  end do
1097 
1098  close(reunt)
1099  close(imunt)
1100 
1101  ABI_DEALLOCATE(nband)
1102  ABI_DEALLOCATE(oml1)
1103  ABI_DEALLOCATE(e2)
1104  ABI_DEALLOCATE(e1)
1105  ABI_DEALLOCATE(occ)
1106  ABI_DEALLOCATE(psinablapsi)
1107  ABI_DEALLOCATE(eigen0)
1108  ABI_DEALLOCATE(wtk)
1109  ABI_DEALLOCATE(kpts)
1110 
1111  call hdr_free(hdr)
1112  call destroy_mpi_enreg(MPI_enreg_seq)
1113 
1114  DBG_EXIT("COLL")
1115 
1116  end subroutine linear_optics_paw

m_paw_optics/optics_paw [ Functions ]

[ Top ] [ m_paw_optics ] [ Functions ]

NAME

 optics_paw

FUNCTION

 Compute matrix elements need for optical conductivity (in the PAW context) and store them in a file
  Matrix elements = <Phi_i|Nabla|Phi_j>

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
  cg(2,mcg)=planewave coefficients of wavefunctions.
  cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk>
                                          and each |p_lmn> non-local projector
  dimcprj(natom)=array of dimensions of array cprj (not ordered)
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  gprimd(3,3)=dimensional reciprocal space primitive translations
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  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.
  nkpt=number of k points.
  npwarr(nkpt)=number of planewaves in basis at this k point
  nsppol=1 for unpolarized, 2 for spin-polarized
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data:
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data

OUTPUT

  (only writing in a file)

SIDE EFFECTS

NOTES

PARENTS

      outscfcv

CHILDREN

      hdr_io,setnabla_ylm,pawcprj_alloc,pawcprj_free,pawcprj_get
      pawcprj_mpi_allgather,pawrad_deducer0,simp_gen,timab,wffclose,wffopen
      xmpi_exch,xmpi_sum,xmpi_sum_master

SOURCE

110  subroutine optics_paw(atindx1,cg,cprj,dimcprj,dtfil,dtset,eigen0,gprimd,hdr,kg,&
111 &               mband,mcg,mcprj,mkmem,mpi_enreg,mpsang,mpw,natom,nkpt,npwarr,nsppol,&
112 &               pawrad,pawtab)
113 
114 
115 !This section has been created automatically by the script Abilint (TD).
116 !Do not modify the following lines by hand.
117 #undef ABI_FUNC
118 #define ABI_FUNC 'optics_paw'
119 !End of the abilint section
120 
121  implicit none
122 
123 !Arguments ------------------------------------
124 !scalars
125  integer,intent(in) :: mband,mcg,mcprj,mkmem,mpsang,mpw,natom,nkpt,nsppol
126  type(MPI_type),intent(in) :: mpi_enreg
127  type(datafiles_type),intent(in) :: dtfil
128  type(dataset_type),intent(in) :: dtset
129  type(hdr_type),intent(inout) :: hdr
130 !arrays
131  integer,intent(in) :: atindx1(natom),dimcprj(natom)
132  integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt)
133  real(dp),intent(in) :: eigen0(mband*nkpt*nsppol)
134  real(dp),intent(in) :: gprimd(3,3)
135  real(dp),intent(inout) :: cg(2,mcg)
136  type(pawcprj_type),target,intent(inout) :: cprj(natom,mcprj)
137  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
138  type(pawtab_type),target,intent(inout) :: pawtab(dtset%ntypat)
139 
140 !Local variables-------------------------------
141 !scalars
142  integer :: iomode,bdtot_index,cplex,etiq,fformopt,iatom,ib,ibg,ibsp
143  integer :: icg,ierr,ikg,ikpt,ilmn,ount
144  integer :: iorder_cprj,ipw,ispinor,isppol,istwf_k,itypat,iwavef
145  integer :: jb,jbsp,jlmn,jwavef,lmn_size,mband_cprj
146  integer :: me,me_kpt,my_nspinor,nband_k,nband_cprj_k,npw_k,sender
147  integer :: spaceComm_band,spaceComm_bandfftspin,spaceComm_fft,spaceComm_k,spaceComm_spin,spaceComm_w
148  logical :: already_has_nabla,cprj_paral_band,mykpt
149  real(dp) :: cgnm1,cgnm2,cpnm1,cpnm2,intg
150  character(len=500) :: message
151 !arrays
152  integer :: tmp_shape(3)
153  integer,allocatable :: kg_k(:,:)
154  real(dp) :: kpoint(3),tsec(2)
155  real(dp),allocatable :: kpg_k(:,:),psinablapsi(:,:,:,:),tnm(:,:,:,:)
156  type(pawcprj_type),pointer :: cprj_k(:,:),cprj_k_loc(:,:)
157  type(wffile_type) :: wff1
158 
159 ! ************************************************************************
160 
161  DBG_ENTER("COLL")
162 
163 !Compatibility tests
164  ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
165 
166 !----------------------------------------------------------------------------------
167 !1- Computation of on-site contribution: <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j>
168 !----------------------------------------------------------------------------------
169 
170  already_has_nabla=all(pawtab(:)%has_nabla==2)
171  call pawnabla_init(mpsang,dtset%ntypat,pawrad,pawtab)
172 
173 !----------------------------------------------------------------------------------
174 !2- Computation of <psi_n|-i.nabla|psi_m> for each k
175 !----------------------------------------------------------------------------------
176 
177 !Init parallelism
178  spaceComm_w=mpi_enreg%comm_cell
179  me=xmpi_comm_rank(spaceComm_w)
180  if (mpi_enreg%paral_kgb==1) then
181    spaceComm_k=mpi_enreg%comm_kpt
182    spaceComm_fft=mpi_enreg%comm_fft
183    spaceComm_band=mpi_enreg%comm_band
184    spaceComm_spin=mpi_enreg%comm_spinor
185    spaceComm_bandfftspin=mpi_enreg%comm_bandspinorfft
186    me_kpt=mpi_enreg%me_kpt
187  else
188    spaceComm_k=spaceComm_w
189    spaceComm_band=0;spaceComm_fft=0;spaceComm_spin=xmpi_comm_self;spaceComm_bandfftspin=0
190    me_kpt=me
191  end if
192  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
193 
194 !Determine if cprj datastructure is distributed over bands
195  mband_cprj=mcprj/(my_nspinor*mkmem*nsppol)
196  cprj_paral_band=(mband_cprj<mband)
197 
198 !PAW file
199  iorder_cprj=0
200 
201 !Initialize main variables
202  ABI_ALLOCATE(psinablapsi,(2,3,mband,mband))
203  psinablapsi=zero
204 
205 !open _OPT file for proc 0
206  iomode=IO_MODE_FORTRAN_MASTER
207  fformopt=610
208  ount = get_unit()
209  call WffOpen(iomode,spaceComm_k,dtfil%fnameabo_app_opt,ierr,wff1,0,me,ount)
210  call hdr_io(fformopt,hdr,2,wff1)
211  if (me==0) then
212    write(ount)(eigen0(ib),ib=1,mband*nkpt*nsppol)
213  end if
214 
215 !LOOP OVER SPINS
216  ibg=0;icg=0
217  bdtot_index=0
218  do isppol=1,nsppol
219 
220 !  LOOP OVER k POINTS
221    ikg=0
222    do ikpt=1,nkpt
223 
224      nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
225      etiq=ikpt+(isppol-1)*nkpt
226 
227 !    Select k-points for current proc
228      mykpt=.true.
229      mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt))
230      if (mykpt) then
231 
232 !      Allocations depending on k-point
233        kpoint(:)=dtset%kptns(:,ikpt)
234        istwf_k=dtset%istwfk(ikpt)
235        npw_k=npwarr(ikpt)
236        cplex=2;if (istwf_k>1) cplex=1
237        ABI_ALLOCATE(kg_k,(3,npw_k))
238        ABI_ALLOCATE(kpg_k,(npw_k*dtset%nspinor,3))
239 
240 !      Extract cprj for this k-point according to mkmem
241        nband_cprj_k=nband_k;if (cprj_paral_band) nband_cprj_k=nband_k/mpi_enreg%nproc_band
242        if (mkmem*nsppol/=1) then
243          ABI_DATATYPE_ALLOCATE(cprj_k_loc,(natom,my_nspinor*nband_cprj_k))
244          call pawcprj_alloc(cprj_k_loc,0,dimcprj)
245          call pawcprj_get(atindx1,cprj_k_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isppol,&
246 &         mband_cprj,mkmem,natom,nband_cprj_k,nband_cprj_k,my_nspinor,nsppol,dtfil%unpaw,&
247 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
248        else
249          cprj_k_loc => cprj
250        end if
251 
252 !      if cprj are distributed over bands, gather them (because we need to mix bands)
253        if (cprj_paral_band) then
254          ABI_DATATYPE_ALLOCATE(cprj_k,(natom,my_nspinor*nband_k))
255          call pawcprj_alloc(cprj_k,0,dimcprj)
256          call pawcprj_mpi_allgather(cprj_k_loc,cprj_k,natom,my_nspinor*nband_cprj_k,dimcprj,0,&
257 &         mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.false.)
258        else
259          cprj_k => cprj_k_loc
260        end if
261 
262 !      Get G-vectors for this k-point.
263        kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
264        ikg=ikg+npw_k
265 
266 !      Calculation of k+G in cartesian coordinates
267        do ipw=1,npw_k
268          kpg_k(ipw,1)=(kpoint(1)+kg_k(1,ipw))*gprimd(1,1)&
269 &         +(kpoint(2)+kg_k(2,ipw))*gprimd(1,2)&
270 &         +(kpoint(3)+kg_k(3,ipw))*gprimd(1,3)
271          kpg_k(ipw,2)=(kpoint(1)+kg_k(1,ipw))*gprimd(2,1)&
272 &         +(kpoint(2)+kg_k(2,ipw))*gprimd(2,2)&
273 &         +(kpoint(3)+kg_k(3,ipw))*gprimd(2,3)
274          kpg_k(ipw,3)=(kpoint(1)+kg_k(1,ipw))*gprimd(3,1)&
275 &         +(kpoint(2)+kg_k(2,ipw))*gprimd(3,2)&
276 &         +(kpoint(3)+kg_k(3,ipw))*gprimd(3,3)
277        end do !ipw
278        kpg_k=two_pi*kpg_k
279        if (dtset%nspinor==2) kpg_k(npw_k+1:2*npw_k,1:3)=kpg_k(1:npw_k,1:3)
280        ABI_DEALLOCATE(kg_k)
281 
282 !      2-A Computation of <psi_tild_n|-i.nabla|psi_tild_m>
283 !      ----------------------------------------------------------------------------------
284 !      Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates
285 
286        ABI_ALLOCATE(tnm,(2,3,nband_k,nband_k))
287        tnm=zero
288 
289 !      Loops on bands
290        do jb=1,nband_k
291          jwavef=(jb-1)*npw_k*my_nspinor+icg
292          if (xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then
293            tmp_shape = shape(mpi_enreg%proc_distrb)
294            if (ikpt > tmp_shape(1)) then
295              message = ' optics_paw : ikpt out of bounds '
296              MSG_ERROR(message)
297            end if
298 !          MJV 6/12/2008: looks like mpi_enreg may not be completely initialized here
299            if (abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-me_kpt)/=0) cycle
300          end if
301          do ib=1,jb
302            iwavef=(ib-1)*npw_k*my_nspinor+icg
303 
304 !          Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates
305            if (cplex==1) then
306              do ipw=1,npw_k*my_nspinor
307                cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)
308                tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3)
309              end do
310            else
311              do ipw=1,npw_k*my_nspinor
312                cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)+cg(2,ipw+iwavef)*cg(2,ipw+jwavef)
313                cgnm2=cg(1,ipw+iwavef)*cg(2,ipw+jwavef)-cg(2,ipw+iwavef)*cg(1,ipw+jwavef)
314                tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3)
315                tnm(2,1:3,ib,jb)=tnm(2,1:3,ib,jb)+cgnm2*kpg_k(ipw,1:3)
316              end do
317            end if
318 
319 !          Second half of the (n,m) matrix
320            if (ib/=jb) then
321              tnm(1,1:3,jb,ib)= tnm(1,1:3,ib,jb)
322              tnm(2,1:3,jb,ib)=-tnm(2,1:3,ib,jb)
323            end if
324 
325          end do ! ib
326        end do ! jb
327 
328 !      ATTEMPT: USING BLAS (not successful)
329 !      allocate(dg(2,npw_k*my_nspinor,nband_k),tnm_tmp(2,nband_k,nband_k))
330 !      do ii=1,3
331 !      do ib=1,nband_k
332 !      iwavef=icg+(ib-1)*npw_k*my_nspinor
333 !      do ipw=1,my_nspinor*npw_k
334 !      dg(:,ipw,nband_k)=cg(:,ipw+iwavef)*kpg_k(ipw,ii)
335 !      end do
336 !      end do
337 !      call zgemm('C','N',nband_k,nband_k,npw_k*my_nspinor,one,dg,npw_k*my_nspinor,&
338 !      &                cg(:,1+icg:npw_k*my_nspinor*nband_k+icg),npw_k*my_nspinor,zero,tnm_tmp,nband_k)
339 !      tnm(:,ii,:,:)=tnm_tmp(:,:,:)
340 !      deallocate(dg,tnm_tmp)
341 !      end do
342 
343 !      Reduction in case of parallelism
344        if (mpi_enreg%paral_kgb == 1) then
345          call timab(48,1,tsec)
346          call xmpi_sum_master(tnm,0,spaceComm_bandfftspin,ierr)
347          call timab(48,2,tsec)
348        end if
349 
350        psinablapsi(:,:,:,:)=tnm(:,:,:,:)
351 
352 !      2-B Computation of <psi_n|p_i><p_j|psi_m>(<phi_i|-i.nabla|phi_j>-<tphi_i|-i.nabla|tphi_j>)
353 !      ----------------------------------------------------------------------------------
354 
355        tnm=zero
356 
357 !      Loops on bands
358        do jb=1,nband_k
359 
360          if (mpi_enreg%paral_kgb==1) then
361            if (mod(jb-1,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle
362            if (abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-me_kpt)/=0) cycle
363          end if
364          do ib=1,jb
365            ibsp=(ib-1)*my_nspinor;jbsp=(jb-1)*my_nspinor
366 
367            if (cplex==1) then
368              do ispinor=1,my_nspinor
369                ibsp=ibsp+1;jbsp=jbsp+1
370                do iatom=1,natom
371                  itypat=dtset%typat(iatom)
372                  lmn_size=pawtab(itypat)%lmn_size
373                  do jlmn=1,lmn_size
374                    do ilmn=1,lmn_size
375                      cpnm1=cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn)
376                      tnm(2,:,ib,jb)=tnm(2,:,ib,jb)+cpnm1*pawtab(itypat)%nabla_ij(:,ilmn,jlmn)
377                    end do !ilmn
378                  end do !jlmn
379                end do !iatom
380              end do !ispinor
381            else
382              do ispinor=1,my_nspinor
383                ibsp=ibsp+1;jbsp=jbsp+1
384                do iatom=1,natom
385                  itypat=dtset%typat(iatom)
386                  lmn_size=pawtab(itypat)%lmn_size
387                  do jlmn=1,lmn_size
388                    do ilmn=1,lmn_size
389                      cpnm1=(cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn) &
390 &                     +cprj_k(iatom,ibsp)%cp(2,ilmn)*cprj_k(iatom,jbsp)%cp(2,jlmn))
391                      cpnm2=(cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(2,jlmn) &
392 &                     -cprj_k(iatom,ibsp)%cp(2,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn))
393                      tnm(1,:,ib,jb)=tnm(1,:,ib,jb)+cpnm2*pawtab(itypat)%nabla_ij(:,ilmn,jlmn)
394                      tnm(2,:,ib,jb)=tnm(2,:,ib,jb)-cpnm1*pawtab(itypat)%nabla_ij(:,ilmn,jlmn)
395                    end do !ilmn
396                  end do !jlmn
397                end do !iatom
398              end do !ispinor
399            end if
400 
401 !          Second half of the (n,m) matrix
402            if (ib/=jb) then
403              tnm(1,1:3,jb,ib)=-tnm(1,1:3,ib,jb)
404              tnm(2,1:3,jb,ib)= tnm(2,1:3,ib,jb)
405            end if
406 
407 !          End loops on bands
408          end do ! jb
409        end do !ib
410 
411 !      Reduction in case of parallelism
412        if (mpi_enreg%paral_kgb==1) then
413          call timab(48,1,tsec)
414          call xmpi_sum(tnm,mpi_enreg%comm_bandspinor,ierr)
415          call timab(48,2,tsec)
416        else
417 !        In that case parallelisation on nspinor not implemented yet;put this line just in case
418          call xmpi_sum(tnm,spacecomm_spin,ierr)
419        end if
420 
421        psinablapsi(:,:,:,:)=psinablapsi(:,:,:,:)+tnm(:,:,:,:)
422        ABI_DEALLOCATE(tnm)
423 
424        if (mkmem/=0) then
425          ibg = ibg +       my_nspinor*nband_cprj_k
426          icg = icg + npw_k*my_nspinor*nband_k
427        end if
428 
429        if (cprj_paral_band) then
430          call pawcprj_free(cprj_k)
431          ABI_DATATYPE_DEALLOCATE(cprj_k)
432        end if
433        if (mkmem*nsppol/=1) then
434          call pawcprj_free(cprj_k_loc)
435          ABI_DATATYPE_DEALLOCATE(cprj_k_loc)
436        end if
437        ABI_DEALLOCATE(kpg_k)
438 
439        if (me==0) then
440          write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k)
441          write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k)
442          write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k)
443        elseif (mpi_enreg%me_band==0.and.mpi_enreg%me_fft==0) then
444          call xmpi_exch(psinablapsi,etiq,me_kpt,psinablapsi,0,spaceComm_k,ierr)
445        end if
446 
447      elseif (me==0) then
448        sender=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
449        call xmpi_exch(psinablapsi,etiq,sender,psinablapsi,0,spaceComm_k,ierr)
450        write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k)
451        write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k)
452        write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k)
453      end if ! mykpt
454 
455      bdtot_index=bdtot_index+nband_k
456 !    End loop on spin,kpt
457    end do ! ikpt
458  end do !isppol
459 
460 !Close file
461  call WffClose(wff1,ierr)
462 
463 !Datastructures deallocations
464  ABI_DEALLOCATE(psinablapsi)
465  if (.not.already_has_nabla) then
466    do itypat=1,dtset%ntypat
467      if (allocated(pawtab(itypat)%nabla_ij)) then
468        ABI_DEALLOCATE(pawtab(itypat)%nabla_ij)
469        pawtab(itypat)%has_nabla=0
470      end if
471    end do
472  end if
473 
474  DBG_EXIT("COLL")
475 
476  end subroutine optics_paw

m_paw_optics/optics_paw_core [ Functions ]

[ Top ] [ m_paw_optics ] [ Functions ]

NAME

 optics_paw_core

FUNCTION

 Compute matrix elements need for X spectr. (in the PAW context) and store them in a file
  Matrix elements = <Phi_core|Nabla|Phi_j>

COPYRIGHT

 Copyright (C) 2005-2018 ABINIT group (SM,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 .

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
  cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
  dimcprj(natom)=array of dimensions of array cprj (not ordered)
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  filpsp(ntypat)=name(s) of the pseudopotential file(s)
  mband=maximum number of bands
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mkmem =number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  mpsang =1+maximum angular momentum for nonlocal pseudopotentials
  natom=number of atoms in cell.
  nkpt=number of k points.
  nsppol=1 for unpolarized, 2 for spin-polarized
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data:
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data

OUTPUT

  (only writing in a file)

PARENTS

      outscfcv

CHILDREN

      hdr_io,setnabla_ylm,pawcprj_alloc,pawcprj_free,pawcprj_get
      pawcprj_mpi_allgather,pawpsp_read_corewf,pawrad_deducer0,simp_gen,timab
      wffclose,wffopen,xmpi_exch,xmpi_sum_master

SOURCE

526  subroutine optics_paw_core(atindx1,cprj,dimcprj,dtfil,dtset,eigen0,filpsp,hdr,&
527 &               mband,mcprj,mkmem,mpi_enreg,mpsang,natom,nkpt,nsppol,pawrad,pawtab)
528 
529 
530 !This section has been created automatically by the script Abilint (TD).
531 !Do not modify the following lines by hand.
532 #undef ABI_FUNC
533 #define ABI_FUNC 'optics_paw_core'
534 !End of the abilint section
535 
536  implicit none
537 
538 !Arguments ------------------------------------
539 !scalars
540  integer,intent(in) :: mband,mcprj,mkmem,mpsang,natom,nkpt,nsppol
541  type(MPI_type),intent(in) :: mpi_enreg
542  type(datafiles_type),intent(in) :: dtfil
543  type(dataset_type),intent(in) :: dtset
544  type(hdr_type),intent(inout) :: hdr
545 !arrays
546  integer,intent(in) :: atindx1(natom),dimcprj(natom)
547  character(len=fnlen),intent(in) :: filpsp(dtset%ntypat)
548  real(dp),intent(in) :: eigen0(mband*nkpt*nsppol)
549  type(pawcprj_type),target,intent(inout) :: cprj(natom,mcprj)
550  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
551  type(pawtab_type),target,intent(inout) :: pawtab(dtset%ntypat)
552 
553 !Local variables-------------------------------
554 !scalars
555  integer :: basis_size,bdtot_index,cplex,etiq,iatom,ib,ibg
556  integer :: ierr,ikpt,il,ilmn,iln,ount
557  integer :: iorder_cprj,ispinor,isppol,istwf_k,itypat
558  integer :: jb,jbsp,jl,jlmn,lmn_size,lmncmax,mband_cprj
559  integer :: me,me_kpt,mesh_size,my_nspinor,nband_cprj_k,nband_k,nphicor
560  integer :: sender,spaceComm_bandspin,spaceComm_k,spaceComm_w
561  integer :: iomode,fformopt
562  logical :: already_has_nabla,cprj_paral_band,ex,mykpt
563  character(len=fnlen) :: filecore
564  real(dp) :: cpnm1,cpnm2,intg
565 !arrays
566  integer,allocatable :: indlmn_core(:,:),lcor(:),ncor(:)
567  real(dp) :: tsec(2)
568  real(dp),allocatable :: energy_cor(:),phi_cor(:,:),psinablapsi(:,:,:,:,:),tnm(:,:,:,:,:)
569  type(pawcprj_type),pointer :: cprj_k(:,:),cprj_k_loc(:,:)
570  type(wffile_type) :: wff1
571 
572 ! ************************************************************************
573 
574  DBG_ENTER("COLL")
575 
576 !Compatibility tests
577  ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
578 
579 !------------------------------------------------------------------------------------------------
580 !0-Reading core wavefunctions
581 !------------------------------------------------------------------------------------------------
582 
583 !Note: core WF is read for itypat=1
584  filecore=trim(filpsp(1))//'.corewf'
585  inquire(file=filecore,exist=ex)
586  if (ex) then
587    !Use <filepsp>.corewf
588    call pawpsp_read_corewf(energy_cor,indlmn_core,lcor,lmncmax,ncor,nphicor,pawrad(1),phi_cor,&
589 &   filename=filecore)
590  else
591    !Use default name
592    call pawpsp_read_corewf(energy_cor,indlmn_core,lcor,lmncmax,ncor,nphicor,pawrad(1),phi_cor)
593  end if
594 
595 !----------------------------------------------------------------------------------
596 !1-Computation of phipphj=<phi_i|nabla|phi_core>
597 !----------------------------------------------------------------------------------
598 
599  already_has_nabla=all(pawtab(:)%has_nabla==3)
600  call pawnabla_core_init(mpsang,dtset%ntypat,pawrad,pawtab,phi_cor,indlmn_core)
601 
602 !----------------------------------------------------------------------------------
603 !2- Computation of <psi_n|p_i>(<phi_i|-i.nabla|phi_core>)
604 !----------------------------------------------------------------------------------
605 
606 !Init parallelism
607  spaceComm_w=mpi_enreg%comm_cell
608  me=xmpi_comm_rank(spaceComm_w)
609  if (mpi_enreg%paral_kgb==1) then
610    spaceComm_k=mpi_enreg%comm_kpt
611    spaceComm_bandspin=mpi_enreg%comm_bandspinor
612    me_kpt=mpi_enreg%me_kpt
613  else
614    spaceComm_k=spaceComm_w
615    spaceComm_bandspin=0
616    me_kpt=me
617  end if
618  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
619 
620 !Determine if cprj datastructure is distributed over bands
621  mband_cprj=mcprj/(my_nspinor*mkmem*nsppol)
622  cprj_paral_band=(mband_cprj<mband)
623 
624 !Prepare temporary PAW file if mkmem==0
625  iorder_cprj=0
626 !Open _OPT2 file for proc 0
627  iomode=IO_MODE_FORTRAN_MASTER
628  fformopt=611
629  ount = get_unit()
630  call WffOpen(iomode,spaceComm_k,dtfil%fnameabo_app_opt2,ierr,wff1,0,me,ount)
631  call hdr_io(fformopt,hdr,2,wff1)
632  if (me==0) then
633    write(ount)(eigen0(ib),ib=1,mband*nkpt*nsppol)
634    write(ount) nphicor
635    do iln=1,nphicor
636      write(ount) ncor(iln),lcor(iln),half*energy_cor(iln)
637    end do
638  end if
639 
640  ABI_DEALLOCATE(ncor)
641  ABI_DEALLOCATE(lcor)
642  ABI_DEALLOCATE(phi_cor)
643  ABI_DEALLOCATE(energy_cor)
644 
645  ABI_ALLOCATE(psinablapsi,(2,3,mband,nphicor,natom))
646 
647 !LOOP OVER SPINS
648  ibg=0
649  bdtot_index=0
650  do isppol=1,nsppol
651 
652 !  LOOP OVER k POINTS
653    do ikpt=1,nkpt
654      nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
655      etiq=ikpt+(isppol-1)*nkpt
656      psinablapsi=zero
657 
658      mykpt=.true.
659      mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt))
660      if (mykpt) then
661 
662 !      Constants depending on k-point
663        istwf_k=dtset%istwfk(ikpt)
664        cplex=2;if (istwf_k>1) cplex=1
665 
666 !      Extract cprj for this k-point.
667        nband_cprj_k=nband_k;if (cprj_paral_band) nband_cprj_k=nband_k/mpi_enreg%nproc_band
668        if (mkmem*nsppol/=1) then
669          ABI_DATATYPE_ALLOCATE(cprj_k_loc,(natom,my_nspinor*nband_cprj_k))
670          call pawcprj_alloc(cprj_k_loc,0,dimcprj)
671          call pawcprj_get(atindx1,cprj_k_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isppol,&
672 &         mband_cprj,mkmem,natom,nband_cprj_k,nband_cprj_k,my_nspinor,nsppol,dtfil%unpaw,&
673 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
674        else
675          cprj_k_loc => cprj
676        end if
677 
678 !      if cprj are distributed over bands, gather them (because we need to mix bands)
679        if (cprj_paral_band) then
680          ABI_DATATYPE_ALLOCATE(cprj_k,(natom,my_nspinor*nband_k))
681          call pawcprj_alloc(cprj_k,0,dimcprj)
682          call pawcprj_mpi_allgather(cprj_k_loc,cprj_k,natom,my_nspinor*nband_cprj_k,dimcprj,0,&
683 &         mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.false.)
684        else
685          cprj_k => cprj_k_loc
686        end if
687 
688        ABI_ALLOCATE(tnm,(2,3,nband_k,nphicor,natom))
689        tnm=zero
690 
691 !      Loops on bands
692        do jb=1,nband_k
693 
694          if (mpi_enreg%paral_kgb==1) then
695            if (mod(jb-1,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle
696          elseif (xmpi_paral==1) then
697            if (abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-me_kpt)/=0) cycle
698          end if
699          jbsp=(jb-1)*my_nspinor
700 
701          if (cplex==1) then
702            do ispinor=1,my_nspinor
703              jbsp=jbsp+1
704              do iatom=1,natom
705                itypat=dtset%typat(iatom)
706                lmn_size=pawtab(itypat)%lmn_size
707                do jlmn=1,lmn_size
708                  do ilmn=1,lmncmax
709                    ib=indlmn_core(5,ilmn)
710                    cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn)
711                    tnm(2,:,jb,ib,iatom)=tnm(2,:,jb,ib,iatom)+cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
712                  end do !ilmn
713                end do !jlmn
714              end do !iatom
715            end do !ispinor
716          else
717            do ispinor=1,my_nspinor
718              jbsp=jbsp+1
719              do iatom=1,natom
720                itypat=dtset%typat(iatom)
721                lmn_size=pawtab(itypat)%lmn_size
722                do jlmn=1,lmn_size
723                  do ilmn=1,lmncmax
724                    ib=indlmn_core(5,ilmn)
725                    cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn)
726                    cpnm2=cprj_k(iatom,jbsp)%cp(2,jlmn)
727                    tnm(1,:,jb,ib,iatom)=tnm(1,:,jb,ib,iatom)+cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
728                    tnm(2,:,jb,ib,iatom)=tnm(2,:,jb,ib,iatom)+cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
729                  end do !ilmn
730                end do !jlmn
731              end do !iatom
732            end do !ispinor
733          end if
734 
735 !        End loops on bands
736        end do ! jb
737 
738 !      Reduction in case of parallelism
739        if (mpi_enreg%paral_kgb==1) then
740          call timab(48,1,tsec)
741          call xmpi_sum_master(tnm,0,spaceComm_bandspin,ierr)
742          call timab(48,2,tsec)
743        end if
744 
745        psinablapsi(:,:,:,:,:)=psinablapsi(:,:,:,:,:)+tnm(:,:,:,:,:)
746        ABI_DEALLOCATE(tnm)
747 
748        if (mkmem/=0) ibg = ibg + my_nspinor*nband_cprj_k
749 
750        if (cprj_paral_band) then
751          call pawcprj_free(cprj_k)
752          ABI_DATATYPE_DEALLOCATE(cprj_k)
753        end if
754        if (mkmem*nsppol/=1) then
755          call pawcprj_free(cprj_k_loc)
756          ABI_DATATYPE_DEALLOCATE(cprj_k_loc)
757        end if
758 
759        if (me==0) then
760          do iatom=1,natom
761            write(ount) ((psinablapsi(1:2,1,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
762            write(ount) ((psinablapsi(1:2,2,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
763            write(ount) ((psinablapsi(1:2,3,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
764          end do
765 !DEBUG
766 !         do iatom=1,natom
767 !           write(138,*) ((psinablapsi(1:2,1,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
768 !           write(138,*) ((psinablapsi(1:2,2,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
769 !           write(138,*) ((psinablapsi(1:2,3,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
770 !         end do
771 !DEBUG
772 
773        elseif (mpi_enreg%me_band==0.and.mpi_enreg%me_fft==0) then
774          call xmpi_exch(psinablapsi,etiq,me_kpt,psinablapsi,0,spaceComm_k,ierr)
775        end if
776 
777      elseif (me==0) then
778        sender=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
779        call xmpi_exch(psinablapsi,etiq,sender,psinablapsi,0,spaceComm_k,ierr)
780        do iatom=1,natom
781          write(ount) ((psinablapsi(1:2,1,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
782          write(ount) ((psinablapsi(1:2,2,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
783          write(ount) ((psinablapsi(1:2,3,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
784        end do
785 
786 !DEBUG
787 !      do iatom=1,natom
788 !         write(138,*) ((psinablapsi(1:2,1,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
789 !         write(138,*) ((psinablapsi(1:2,2,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
790 !         write(138,*) ((psinablapsi(1:2,3,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor)
791 !       end do
792 !DEBUG
793 
794      end if ! mykpt
795 
796      bdtot_index=bdtot_index+nband_k
797 
798 !    End loop on spin,kpt
799    end do ! ikpt
800  end do !isppol
801 
802 !Close file
803  call WffClose(wff1,ierr)
804 
805 !Datastructures deallocations
806  ABI_DEALLOCATE(indlmn_core)
807  ABI_DEALLOCATE(psinablapsi)
808  if (.not.already_has_nabla) then
809    do itypat=1,dtset%ntypat
810      if (allocated(pawtab(itypat)%nabla_ij)) then
811        ABI_DEALLOCATE(pawtab(itypat)%nabla_ij)
812        pawtab(itypat)%has_nabla=0
813      end if
814    end do
815  end if
816 
817  DBG_EXIT("COLL")
818 
819  end subroutine optics_paw_core