TABLE OF CONTENTS
ABINIT/optics_paw [ 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