TABLE OF CONTENTS


ABINIT/optics_paw [ Functions ]

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

COPYRIGHT

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

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,int_ang,nderiv_gen,pawcprj_alloc,pawcprj_free,pawcprj_get
      pawcprj_mpi_allgather,pawrad_deducer0,simp_gen,timab,wffclose,wffopen
      xmpi_exch,xmpi_sum,xmpi_sum_master

SOURCE

 57 #if defined HAVE_CONFIG_H
 58 #include "config.h"
 59 #endif
 60 
 61 #include "abi_common.h"
 62 
 63  subroutine optics_paw(atindx1,cg,cprj,dimcprj,dtfil,dtset,eigen0,gprimd,hdr,kg,&
 64 &               mband,mcg,mcprj,mkmem,mpi_enreg,mpsang,mpw,natom,nkpt,npwarr,nsppol,&
 65 &               pawrad,pawtab)
 66 
 67  use defs_basis
 68  use defs_abitypes
 69  use m_xmpi
 70  use m_errors
 71  use m_wffile
 72  use m_profiling_abi
 73  use m_hdr
 74 
 75  use m_io_tools,  only : get_unit
 76  use m_pawrad,    only : pawrad_type, pawrad_deducer0, simp_gen, nderiv_gen
 77  use m_pawtab,    only : pawtab_type
 78  use m_pawcprj,   only : pawcprj_type, pawcprj_alloc, pawcprj_get, &
 79 &                        pawcprj_free,pawcprj_mpi_allgather
 80 
 81 !This section has been created automatically by the script Abilint (TD).
 82 !Do not modify the following lines by hand.
 83 #undef ABI_FUNC
 84 #define ABI_FUNC 'optics_paw'
 85  use interfaces_18_timing
 86  use interfaces_32_util
 87  use interfaces_65_paw, except_this_one => optics_paw
 88 !End of the abilint section
 89 
 90  implicit none
 91 
 92 !Arguments ------------------------------------
 93 !scalars
 94  integer,intent(in) :: mband,mcg,mcprj,mkmem,mpsang,mpw,natom,nkpt,nsppol
 95  type(MPI_type),intent(in) :: mpi_enreg
 96  type(datafiles_type),intent(in) :: dtfil
 97  type(dataset_type),intent(in) :: dtset
 98  type(hdr_type),intent(inout) :: hdr
 99 !arrays
100  integer,intent(in) :: atindx1(natom),dimcprj(natom)
101  integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt)
102  real(dp),intent(in) :: eigen0(mband*nkpt*nsppol)
103  real(dp),intent(in) :: gprimd(3,3)
104  real(dp),intent(inout) :: cg(2,mcg)
105  type(pawcprj_type),target,intent(inout) :: cprj(natom,mcprj)
106  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
107  type(pawtab_type),target,intent(in) :: pawtab(dtset%ntypat)
108 
109 !Local variables-------------------------------
110 !scalars
111  integer :: iomode,basis_size,bdtot_index,cplex,etiq,fformopt,iatom,ib,ibg,ibsp
112  integer :: icg,ierr,ij_size,ikg,ikpt,il,ilm,ilmn,iln,ount
113  integer :: iorder_cprj,ipw,ispinor,isppol,istwf_k,itypat,iwavef
114  integer :: jb,jbsp,jl,jlm,jlmn,jln,jwavef,lmn_size,mband_cprj
115  integer :: me,me_kpt,mesh_size,my_nspinor,nband_k,nband_cprj_k,npw_k,sender
116  integer :: spaceComm_band,spaceComm_bandfftspin,spaceComm_fft,spaceComm_k,spaceComm_spin,spaceComm_w
117  logical :: cprj_paral_band,mykpt
118  real(dp) :: cgnm1,cgnm2,cpnm1,cpnm2,intg
119  character(len=500) :: message
120 !arrays
121  integer :: tmp_shape(3)
122  integer,allocatable :: kg_k(:,:)
123  integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
124  real(dp) :: ang_phipphj(mpsang**2,mpsang**2,8),kpoint(3)
125  real(dp) :: tsec(2)
126  real(dp),allocatable :: dphi(:),dtphi(:)
127  real(dp),allocatable ::ff(:),int1(:,:),int2(:,:)
128  real(dp),allocatable :: kpg_k(:,:),phidphj(:,:),psinablapsi(:,:,:,:)
129  real(dp),allocatable :: rad(:),tnm(:,:,:,:),tphidtphj(:,:)
130  type(coeff3_type), allocatable :: phipphj(:)
131  type(pawcprj_type),pointer :: cprj_k(:,:),cprj_k_loc(:,:)
132  type(wffile_type) :: wff1
133 
134 ! ************************************************************************
135 
136  DBG_ENTER("COLL")
137 
138 !Compatibility tests
139  ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
140 
141 !----------------------------------------------------------------------------------
142 !1- Computation of phipphj=<phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j>
143 !----------------------------------------------------------------------------------
144 
145 !1-A Integration of the angular part : all angular integrals have been
146 !computed outside Abinit and tabulated for each (l,m) value
147 !----------------------------------------------------------------------------------
148 
149  call int_ang(ang_phipphj,mpsang)
150 
151  ABI_DATATYPE_ALLOCATE(phipphj,(dtset%ntypat))
152 
153 !loop on atoms type
154  do itypat=1,dtset%ntypat
155 
156    mesh_size=pawtab(itypat)%mesh_size
157    lmn_size=pawtab(itypat)%lmn_size
158    basis_size=pawtab(itypat)%basis_size
159    ij_size=lmn_size*lmn_size
160 
161    ABI_ALLOCATE(ff,(mesh_size))
162    ABI_ALLOCATE(rad,(mesh_size))
163    ABI_ALLOCATE(int2,(lmn_size,lmn_size))
164    ABI_ALLOCATE(int1,(lmn_size,lmn_size))
165    ABI_ALLOCATE(dphi,(mesh_size))
166    ABI_ALLOCATE(dtphi,(mesh_size))
167    ABI_ALLOCATE(phidphj,(mesh_size,ij_size))
168    ABI_ALLOCATE(tphidtphj,(mesh_size,ij_size))
169    ABI_ALLOCATE(phipphj(itypat)%value,(3,lmn_size,lmn_size))
170 
171    indlmn => pawtab(itypat)%indlmn
172    rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size)
173 
174 !  1-B  Computation of int1=\int phi phj /r dr - \int tphi tphj /r dr
175 !  ----------------------------------------------------------------------------------
176    do jln=1,basis_size
177      do iln=1,basis_size
178        ff(2:mesh_size)=(pawtab(itypat)%phi(2:mesh_size,iln)*pawtab(itypat)%phi(2:mesh_size,jln)&
179 &       -pawtab(itypat)%tphi(2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln))/rad(2:mesh_size)
180        call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
181        call simp_gen(intg,ff,pawrad(itypat))
182        int1(iln,jln)=intg
183      end do
184    end do
185 
186 !  1-C Computation of int2=\int phi/r d/dr(phj/r) r^2dr - \int tphi/r d/dr(tphj/r)r^2 dr
187 !  ----------------------------------------------------------------------------------
188    do jln=1,basis_size
189 
190      ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,jln)
191      call nderiv_gen(dphi,ff,pawrad(itypat))
192      ff(1:mesh_size)=pawtab(itypat)%tphi(1:mesh_size,jln)
193      call nderiv_gen(dtphi,ff,pawrad(itypat))
194 
195      do iln=1,basis_size
196        ff(2:mesh_size)=pawtab(itypat)%phi(2:mesh_size,iln)*dphi(2:mesh_size) &
197 &       -pawtab(itypat)%phi (2:mesh_size,iln)*pawtab(itypat)%phi(2:mesh_size,jln)/ &
198 &       rad(2:mesh_size)-(pawtab(itypat)%tphi(2:mesh_size,iln)*dtphi(2:mesh_size) &
199 &       -pawtab(itypat)%tphi (2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln)/rad(2:mesh_size))
200        call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
201        call simp_gen(intg,ff,pawrad(itypat))
202        int2(iln,jln)=intg
203      end do
204    end do
205 
206 !  1-D Integration of the radial part
207 !  ----------------------------------------------------------------------------------
208    do jlmn=1,lmn_size
209      jlm=indlmn(4,jlmn)
210      jl=indlmn(5,jlmn)
211      do ilmn=1,lmn_size
212        ilm=indlmn(4,ilmn)
213        il=indlmn(5,ilmn)
214        phipphj(itypat)%value(1,ilmn,jlmn)= int2(il,jl)*ang_phipphj(ilm,jlm,1)&
215 &       + int1(il,jl)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3))
216        phipphj(itypat)%value(2,ilmn,jlmn)= int2(il,jl)*ang_phipphj(ilm,jlm,4)&
217 &       + int1(il,jl)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6))
218        phipphj(itypat)%value(3,ilmn,jlmn)= int2(il,jl)*ang_phipphj(ilm,jlm,7)&
219 &       + int1(il,jl)*ang_phipphj(ilm,jlm,8)
220      end do
221    end do
222 
223    ABI_DEALLOCATE(ff)
224    ABI_DEALLOCATE(rad)
225    ABI_DEALLOCATE(int2)
226    ABI_DEALLOCATE(int1)
227    ABI_DEALLOCATE(dphi)
228    ABI_DEALLOCATE(dtphi)
229    ABI_DEALLOCATE(phidphj)
230    ABI_DEALLOCATE(tphidtphj)
231 
232 !  end loop on atoms type
233  end do
234 
235 !----------------------------------------------------------------------------------
236 !2- Computation of <psi_n|-i.nabla|psi_m> for each k
237 !----------------------------------------------------------------------------------
238 
239 !Init parallelism
240  spaceComm_w=mpi_enreg%comm_cell
241  me=xmpi_comm_rank(spaceComm_w)
242  if (mpi_enreg%paral_kgb==1) then
243    spaceComm_k=mpi_enreg%comm_kpt
244    spaceComm_fft=mpi_enreg%comm_fft
245    spaceComm_band=mpi_enreg%comm_band
246    spaceComm_spin=mpi_enreg%comm_spinor
247    spaceComm_bandfftspin=mpi_enreg%comm_bandspinorfft
248    me_kpt=mpi_enreg%me_kpt
249  else
250    spaceComm_k=spaceComm_w
251    spaceComm_band=0;spaceComm_fft=0;spaceComm_spin=xmpi_comm_self;spaceComm_bandfftspin=0
252    me_kpt=me
253  end if
254  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
255 
256 !Determine if cprj datastructure is distributed over bands
257  mband_cprj=mcprj/(my_nspinor*mkmem*nsppol)
258  cprj_paral_band=(mband_cprj<mband)
259 
260 !PAW file
261  iorder_cprj=0
262 
263 !Initialize main variables
264  ABI_ALLOCATE(psinablapsi,(2,3,mband,mband))
265  psinablapsi=zero
266 
267 !open _OPT file for proc 0
268  iomode=IO_MODE_FORTRAN_MASTER
269  fformopt=610
270  ount = get_unit()
271  call WffOpen(iomode,spaceComm_k,dtfil%fnameabo_app_opt,ierr,wff1,0,me,ount)
272  call hdr_io(fformopt,hdr,2,wff1)
273  if (me==0) then
274    write(ount)(eigen0(ib),ib=1,mband*nkpt*nsppol)
275  end if
276 
277 !LOOP OVER SPINS
278  ibg=0;icg=0
279  bdtot_index=0
280  do isppol=1,nsppol
281 
282 !  LOOP OVER k POINTS
283    ikg=0
284    do ikpt=1,nkpt
285 
286      nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
287      etiq=ikpt+(isppol-1)*nkpt
288 
289 !    Select k-points for current proc
290      mykpt=.true.
291      mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt))
292      if (mykpt) then
293 
294 !      Allocations depending on k-point
295        kpoint(:)=dtset%kptns(:,ikpt)
296        istwf_k=dtset%istwfk(ikpt)
297        npw_k=npwarr(ikpt)
298        cplex=2;if (istwf_k>1) cplex=1
299        ABI_ALLOCATE(kg_k,(3,npw_k))
300        ABI_ALLOCATE(kpg_k,(npw_k*dtset%nspinor,3))
301 
302 !      Extract cprj for this k-point according to mkmem
303        nband_cprj_k=nband_k;if (cprj_paral_band) nband_cprj_k=nband_k/mpi_enreg%nproc_band
304        if (mkmem*nsppol/=1) then
305          ABI_DATATYPE_ALLOCATE(cprj_k_loc,(natom,my_nspinor*nband_cprj_k))
306          call pawcprj_alloc(cprj_k_loc,0,dimcprj)
307          call pawcprj_get(atindx1,cprj_k_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isppol,&
308 &         mband_cprj,mkmem,natom,nband_cprj_k,nband_cprj_k,my_nspinor,nsppol,dtfil%unpaw,&
309 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
310        else
311          cprj_k_loc => cprj
312        end if
313 
314 !      if cprj are distributed over bands, gather them (because we need to mix bands)
315        if (cprj_paral_band) then
316          ABI_DATATYPE_ALLOCATE(cprj_k,(natom,my_nspinor*nband_k))
317          call pawcprj_alloc(cprj_k,0,dimcprj)
318          call pawcprj_mpi_allgather(cprj_k_loc,cprj_k,natom,my_nspinor*nband_cprj_k,dimcprj,0,&
319 &         mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.false.)
320        else
321          cprj_k => cprj_k_loc
322        end if
323 
324 !      Get G-vectors for this k-point.
325        kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
326        ikg=ikg+npw_k
327 
328 !      Calculation of k+G in cartesian coordinates
329        do ipw=1,npw_k
330          kpg_k(ipw,1)=(kpoint(1)+kg_k(1,ipw))*gprimd(1,1)&
331 &         +(kpoint(2)+kg_k(2,ipw))*gprimd(1,2)&
332 &         +(kpoint(3)+kg_k(3,ipw))*gprimd(1,3)
333          kpg_k(ipw,2)=(kpoint(1)+kg_k(1,ipw))*gprimd(2,1)&
334 &         +(kpoint(2)+kg_k(2,ipw))*gprimd(2,2)&
335 &         +(kpoint(3)+kg_k(3,ipw))*gprimd(2,3)
336          kpg_k(ipw,3)=(kpoint(1)+kg_k(1,ipw))*gprimd(3,1)&
337 &         +(kpoint(2)+kg_k(2,ipw))*gprimd(3,2)&
338 &         +(kpoint(3)+kg_k(3,ipw))*gprimd(3,3)
339        end do !ipw
340        kpg_k=two_pi*kpg_k
341        if (dtset%nspinor==2) kpg_k(npw_k+1:2*npw_k,1:3)=kpg_k(1:npw_k,1:3)
342        ABI_DEALLOCATE(kg_k)
343 
344 !      2-A Computation of <psi_tild_n|-i.nabla|psi_tild_m>
345 !      ----------------------------------------------------------------------------------
346 !      Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates
347 
348        ABI_ALLOCATE(tnm,(2,3,nband_k,nband_k))
349        tnm=zero
350 
351 !      Loops on bands
352        do jb=1,nband_k
353          jwavef=(jb-1)*npw_k*my_nspinor+icg
354          if (xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then
355            tmp_shape = shape(mpi_enreg%proc_distrb)
356            if (ikpt > tmp_shape(1)) then
357              message = ' optics_paw : ikpt out of bounds '
358              MSG_ERROR(message)
359            end if
360 !          MJV 6/12/2008: looks like mpi_enreg may not be completely initialized here
361            if (abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-me_kpt)/=0) cycle
362          end if
363          do ib=1,jb
364            iwavef=(ib-1)*npw_k*my_nspinor+icg
365 
366 !          Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates
367            if (cplex==1) then
368              do ipw=1,npw_k*my_nspinor
369                cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)
370                tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3)
371              end do
372            else
373              do ipw=1,npw_k*my_nspinor
374                cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)+cg(2,ipw+iwavef)*cg(2,ipw+jwavef)
375                cgnm2=cg(1,ipw+iwavef)*cg(2,ipw+jwavef)-cg(2,ipw+iwavef)*cg(1,ipw+jwavef)
376                tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3)
377                tnm(2,1:3,ib,jb)=tnm(2,1:3,ib,jb)+cgnm2*kpg_k(ipw,1:3)
378              end do
379            end if
380 
381 !          Second half of the (n,m) matrix
382            if (ib/=jb) then
383              tnm(1,1:3,jb,ib)= tnm(1,1:3,ib,jb)
384              tnm(2,1:3,jb,ib)=-tnm(2,1:3,ib,jb)
385            end if
386 
387          end do ! ib
388        end do ! jb
389 
390 !      ATTEMPT: USING BLAS (not successful)
391 !      allocate(dg(2,npw_k*my_nspinor,nband_k),tnm_tmp(2,nband_k,nband_k))
392 !      do ii=1,3
393 !      do ib=1,nband_k
394 !      iwavef=icg+(ib-1)*npw_k*my_nspinor
395 !      do ipw=1,my_nspinor*npw_k
396 !      dg(:,ipw,nband_k)=cg(:,ipw+iwavef)*kpg_k(ipw,ii)
397 !      end do
398 !      end do
399 !      call zgemm('C','N',nband_k,nband_k,npw_k*my_nspinor,one,dg,npw_k*my_nspinor,&
400 !      &                cg(:,1+icg:npw_k*my_nspinor*nband_k+icg),npw_k*my_nspinor,zero,tnm_tmp,nband_k)
401 !      tnm(:,ii,:,:)=tnm_tmp(:,:,:)
402 !      deallocate(dg,tnm_tmp)
403 !      end do
404 
405 !      Reduction in case of parallelism
406        if (mpi_enreg%paral_kgb == 1) then
407          call timab(48,1,tsec)
408          call xmpi_sum_master(tnm,0,spaceComm_bandfftspin,ierr)
409          call timab(48,2,tsec)
410        end if
411 
412        psinablapsi(:,:,:,:)=tnm(:,:,:,:)
413 
414 !      2-B Computation of <psi_n|p_i><p_j|psi_m>(<phi_i|-i.nabla|phi_j>-<tphi_i|-i.nabla|tphi_j>)
415 !      ----------------------------------------------------------------------------------
416 
417        tnm=zero
418 
419 !      Loops on bands
420        do jb=1,nband_k
421 
422          if (mpi_enreg%paral_kgb==1) then
423            if (mod(jb-1,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle
424            if (abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-me_kpt)/=0) cycle
425          end if
426          do ib=1,jb
427            ibsp=(ib-1)*my_nspinor;jbsp=(jb-1)*my_nspinor
428 
429            if (cplex==1) then
430              do ispinor=1,my_nspinor
431                ibsp=ibsp+1;jbsp=jbsp+1
432                do iatom=1,natom
433                  itypat=dtset%typat(iatom)
434                  lmn_size=pawtab(itypat)%lmn_size
435                  do jlmn=1,lmn_size
436                    do ilmn=1,lmn_size
437                      cpnm1=cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn)
438                      tnm(2,:,ib,jb)=tnm(2,:,ib,jb)+cpnm1*phipphj(itypat)%value(:,ilmn,jlmn)
439                    end do !ilmn
440                  end do !jlmn
441                end do !iatom
442              end do !ispinor
443            else
444              do ispinor=1,my_nspinor
445                ibsp=ibsp+1;jbsp=jbsp+1
446                do iatom=1,natom
447                  itypat=dtset%typat(iatom)
448                  lmn_size=pawtab(itypat)%lmn_size
449                  do jlmn=1,lmn_size
450                    do ilmn=1,lmn_size
451                      cpnm1=(cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn) &
452 &                     +cprj_k(iatom,ibsp)%cp(2,ilmn)*cprj_k(iatom,jbsp)%cp(2,jlmn))
453                      cpnm2=(cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(2,jlmn) &
454 &                     -cprj_k(iatom,ibsp)%cp(2,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn))
455                      tnm(1,:,ib,jb)=tnm(1,:,ib,jb)+cpnm2*phipphj(itypat)%value(:,ilmn,jlmn)
456                      tnm(2,:,ib,jb)=tnm(2,:,ib,jb)-cpnm1*phipphj(itypat)%value(:,ilmn,jlmn)
457                    end do !ilmn
458                  end do !jlmn
459                end do !iatom
460              end do !ispinor
461            end if
462 
463 !          Second half of the (n,m) matrix
464            if (ib/=jb) then
465              tnm(1,1:3,jb,ib)=-tnm(1,1:3,ib,jb)
466              tnm(2,1:3,jb,ib)= tnm(2,1:3,ib,jb)
467            end if
468 
469 !          End loops on bands
470          end do ! jb
471        end do !ib
472 
473 !      Reduction in case of parallelism
474        if (mpi_enreg%paral_kgb==1) then
475          call timab(48,1,tsec)
476          call xmpi_sum(tnm,mpi_enreg%comm_bandspinor,ierr)
477          call timab(48,2,tsec)
478        else
479 !        In that case parallelisation on nspinor not implemented yet;put this line just in case
480          call xmpi_sum(tnm,spacecomm_spin,ierr)
481        end if
482 
483        psinablapsi(:,:,:,:)=psinablapsi(:,:,:,:)+tnm(:,:,:,:)
484        ABI_DEALLOCATE(tnm)
485 
486        if (mkmem/=0) then
487          ibg = ibg +       my_nspinor*nband_cprj_k
488          icg = icg + npw_k*my_nspinor*nband_k
489        end if
490 
491        if (cprj_paral_band) then
492          call pawcprj_free(cprj_k)
493          ABI_DATATYPE_DEALLOCATE(cprj_k)
494        end if
495        if (mkmem*nsppol/=1) then
496          call pawcprj_free(cprj_k_loc)
497          ABI_DATATYPE_DEALLOCATE(cprj_k_loc)
498        end if
499        ABI_DEALLOCATE(kpg_k)
500 
501        if (me==0) then
502          write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k)
503          write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k)
504          write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k)
505        elseif (mpi_enreg%me_band==0.and.mpi_enreg%me_fft==0) then
506          call xmpi_exch(psinablapsi,etiq,me_kpt,psinablapsi,0,spaceComm_k,ierr)
507        end if
508 
509      elseif (me==0) then
510        sender=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
511        call xmpi_exch(psinablapsi,etiq,sender,psinablapsi,0,spaceComm_k,ierr)
512        write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k)
513        write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k)
514        write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k)
515      end if ! mykpt
516 
517      bdtot_index=bdtot_index+nband_k
518 !    End loop on spin,kpt
519    end do ! ikpt
520  end do !isppol
521 
522 !Close file
523  call WffClose(wff1,ierr)
524 
525 !Datastructures deallocations
526  do itypat=1,dtset%ntypat
527    ABI_DEALLOCATE(phipphj(itypat)%value)
528  end do
529  ABI_DATATYPE_DEALLOCATE(phipphj)
530  ABI_DEALLOCATE(psinablapsi)
531 
532  DBG_EXIT("COLL")
533 
534  end subroutine optics_paw