TABLE OF CONTENTS
ABINIT/optics_paw_core [ 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,int_ang,nderiv_gen,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
47 #if defined HAVE_CONFIG_H 48 #include "config.h" 49 #endif 50 51 #include "abi_common.h" 52 53 subroutine optics_paw_core(atindx1,cprj,dimcprj,dtfil,dtset,eigen0,filpsp,hdr,& 54 & mband,mcprj,mkmem,mpi_enreg,mpsang,natom,nkpt,nsppol,pawrad,pawtab) 55 56 57 use defs_basis 58 use defs_abitypes 59 use m_profiling_abi 60 use m_xmpi 61 use m_errors 62 use m_wffile 63 use m_hdr 64 65 use m_io_tools, only : get_unit 66 use m_pawpsp, only : pawpsp_read_corewf 67 use m_pawtab, only : pawtab_type 68 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, & 69 & pawcprj_free, pawcprj_mpi_allgather 70 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free, pawrad_ifromr, & 71 & pawrad_deducer0, simp_gen, nderiv_gen, bound_deriv 72 use m_splines, only : spline,splint 73 74 !This section has been created automatically by the script Abilint (TD). 75 !Do not modify the following lines by hand. 76 #undef ABI_FUNC 77 #define ABI_FUNC 'optics_paw_core' 78 use interfaces_18_timing 79 use interfaces_32_util 80 use interfaces_65_paw, except_this_one => optics_paw_core 81 !End of the abilint section 82 83 implicit none 84 85 !Arguments ------------------------------------ 86 !scalars 87 integer,intent(in) :: mband,mcprj,mkmem,mpsang,natom,nkpt,nsppol 88 type(MPI_type),intent(in) :: mpi_enreg 89 type(datafiles_type),intent(in) :: dtfil 90 type(dataset_type),intent(in) :: dtset 91 type(hdr_type),intent(inout) :: hdr 92 !arrays 93 integer,intent(in) :: atindx1(natom),dimcprj(natom) 94 character(len=fnlen),intent(in) :: filpsp(dtset%ntypat) 95 real(dp),intent(in) :: eigen0(mband*nkpt*nsppol) 96 type(pawcprj_type),target,intent(inout) :: cprj(natom,mcprj) 97 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat) 98 type(pawtab_type),target,intent(in) :: pawtab(dtset%ntypat) 99 100 !Local variables------------------------------- 101 !scalars 102 integer :: basis_size,bdtot_index,cplex,etiq,iatom,ib,ibg 103 integer :: ierr,ij_size,ikpt,il,ilm,ilmn,iln,ount 104 integer :: iorder_cprj,ispinor,isppol,istwf_k,itypat 105 integer :: jb,jbsp,jl,jlm,jlmn,jln,lmn_size,lmncmax,mband_cprj 106 integer :: me,me_kpt,mesh_size,my_nspinor,nband_cprj_k,nband_k,nphicor 107 integer :: sender,spaceComm_bandspin,spaceComm_k,spaceComm_w 108 integer :: iomode,fformopt 109 logical :: cprj_paral_band,ex,mykpt 110 character(len=fnlen) :: filecore 111 real(dp) :: cpnm1,cpnm2,intg 112 !arrays 113 integer,allocatable :: indlmn_core(:,:),lcor(:),ncor(:) 114 integer,ABI_CONTIGUOUS pointer :: indlmn(:,:) 115 real(dp) :: ang_phipphj(mpsang**2,mpsang**2,8) 116 real(dp) :: tsec(2) 117 real(dp),allocatable :: dphi(:),energy_cor(:),ff(:),int1(:,:),int2(:,:) 118 real(dp),allocatable :: rad(:),phi_cor(:,:),psinablapsi(:,:,:,:,:) 119 real(dp),allocatable :: tnm(:,:,:,:,:) 120 type(coeff3_type), allocatable :: phipphj(:) 121 type(pawcprj_type),pointer :: cprj_k(:,:),cprj_k_loc(:,:) 122 type(wffile_type) :: wff1 123 124 ! ************************************************************************ 125 126 DBG_ENTER("COLL") 127 128 !Compatibility tests 129 ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!") 130 131 !------------------------------------------------------------------------------------------------ 132 !0-Reading core wavefunctions 133 !------------------------------------------------------------------------------------------------ 134 135 !Note: core WF is read for itypat=1 136 filecore=trim(filpsp(1))//'.corewf' 137 inquire(file=filecore,exist=ex) 138 if (ex) then 139 !Use <filepsp>.corewf 140 call pawpsp_read_corewf(energy_cor,indlmn_core,lcor,lmncmax,ncor,nphicor,pawrad(1),phi_cor,& 141 & filename=filecore) 142 else 143 !Use default name 144 call pawpsp_read_corewf(energy_cor,indlmn_core,lcor,lmncmax,ncor,nphicor,pawrad(1),phi_cor) 145 end if 146 147 !---------------------------------------------------------------------------------- 148 !1-Computation of phipphj=<phi_i|nabla|phi_core> 149 !---------------------------------------------------------------------------------- 150 151 !1-A Integration of the angular part : all angular integrals have been 152 !computed outside Abinit and tabulated for each (l,m) value 153 !---------------------------------------------------------------------------------- 154 155 call int_ang(ang_phipphj,mpsang) 156 157 ABI_DATATYPE_ALLOCATE(phipphj,(dtset%ntypat)) 158 159 !We consider the impurity to be the first atom 160 !loop on atoms type 161 do itypat=1,dtset%ntypat 162 163 mesh_size=pawtab(itypat)%mesh_size 164 lmn_size=pawtab(itypat)%lmn_size 165 basis_size=pawtab(itypat)%basis_size 166 ij_size=lmn_size*lmn_size 167 indlmn => pawtab(itypat)%indlmn 168 169 ABI_ALLOCATE(ff,(mesh_size)) 170 ABI_ALLOCATE(rad,(mesh_size)) 171 ABI_ALLOCATE(int2,(lmn_size,lmncmax)) 172 ABI_ALLOCATE(int1,(lmn_size,lmncmax)) 173 ABI_ALLOCATE(dphi,(mesh_size)) 174 ABI_ALLOCATE(phipphj(itypat)%value,(3,lmn_size,lmncmax)) 175 176 rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size) 177 178 ! 1-B Computation of int1=\int phi phi_core /r dr 179 ! ---------------------------------------------------------------------------------- 180 do jln=1,nphicor 181 do iln=1,basis_size 182 ff(2:mesh_size)=(pawtab(itypat)%phi(2:mesh_size,iln)*phi_cor(2:mesh_size,jln))/rad(2:mesh_size) 183 call pawrad_deducer0(ff,mesh_size,pawrad(itypat)) 184 call simp_gen(intg,ff,pawrad(itypat)) 185 int1(iln,jln)=intg 186 end do 187 end do 188 189 ! 1-C Computation of int2=\int phi/r d/dr(phi_core/r) - phi phi_core/r dr 190 ! ---------------------------------------------------------------------------------- 191 do jln=1,nphicor 192 ff(1:mesh_size)=phi_cor(1:mesh_size,jln) 193 call nderiv_gen(dphi,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)*phi_cor(2:mesh_size,jln)/ & 198 & rad(2:mesh_size) 199 call pawrad_deducer0(ff,mesh_size,pawrad(itypat)) 200 call simp_gen(intg,ff,pawrad(itypat)) 201 int2(iln,jln)=intg 202 end do 203 end do 204 205 ! 1-D Integration of the radial part 206 ! ---------------------------------------------------------------------------------- 207 do jlmn=1,lmncmax 208 jlm=indlmn_core(4,jlmn) 209 jl=indlmn_core(5,jlmn) 210 do ilmn=1,lmn_size 211 ilm=indlmn(4,ilmn) 212 il=indlmn(5,ilmn) 213 phipphj(itypat)%value(1,ilmn,jlmn)= int2(il,jl)*ang_phipphj(ilm,jlm,1)& 214 & + int1(il,jl)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3)) 215 phipphj(itypat)%value(2,ilmn,jlmn)= int2(il,jl)*ang_phipphj(ilm,jlm,4)& 216 & + int1(il,jl)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6)) 217 phipphj(itypat)%value(3,ilmn,jlmn)= int2(il,jl)*ang_phipphj(ilm,jlm,7)& 218 & + int1(il,jl)*ang_phipphj(ilm,jlm,8) 219 end do 220 end do 221 222 ABI_DEALLOCATE(ff) 223 ABI_DEALLOCATE(rad) 224 ABI_DEALLOCATE(int2) 225 ABI_DEALLOCATE(int1) 226 ABI_DEALLOCATE(dphi) 227 228 ! end loop on atoms type 229 end do 230 231 !---------------------------------------------------------------------------------- 232 !2- Computation of <psi_n|p_i>(<phi_i|-i.nabla|phi_core>) 233 !---------------------------------------------------------------------------------- 234 235 !Init parallelism 236 spaceComm_w=mpi_enreg%comm_cell 237 me=xmpi_comm_rank(spaceComm_w) 238 if (mpi_enreg%paral_kgb==1) then 239 spaceComm_k=mpi_enreg%comm_kpt 240 spaceComm_bandspin=mpi_enreg%comm_bandspinor 241 me_kpt=mpi_enreg%me_kpt 242 else 243 spaceComm_k=spaceComm_w 244 spaceComm_bandspin=0 245 me_kpt=me 246 end if 247 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 248 249 !Determine if cprj datastructure is distributed over bands 250 mband_cprj=mcprj/(my_nspinor*mkmem*nsppol) 251 cprj_paral_band=(mband_cprj<mband) 252 253 !Prepare temporary PAW file if mkmem==0 254 iorder_cprj=0 255 !Open _OPT2 file for proc 0 256 iomode=IO_MODE_FORTRAN_MASTER 257 fformopt=611 258 ount = get_unit() 259 call WffOpen(iomode,spaceComm_k,dtfil%fnameabo_app_opt2,ierr,wff1,0,me,ount) 260 call hdr_io(fformopt,hdr,2,wff1) 261 if (me==0) then 262 write(ount)(eigen0(ib),ib=1,mband*nkpt*nsppol) 263 write(ount) nphicor 264 do iln=1,nphicor 265 write(ount) ncor(iln),lcor(iln),half*energy_cor(iln) 266 end do 267 end if 268 269 ABI_DEALLOCATE(ncor) 270 ABI_DEALLOCATE(lcor) 271 ABI_DEALLOCATE(phi_cor) 272 ABI_DEALLOCATE(energy_cor) 273 274 ABI_ALLOCATE(psinablapsi,(2,3,mband,nphicor,natom)) 275 276 !LOOP OVER SPINS 277 ibg=0 278 bdtot_index=0 279 do isppol=1,nsppol 280 281 ! LOOP OVER k POINTS 282 do ikpt=1,nkpt 283 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 284 etiq=ikpt+(isppol-1)*nkpt 285 psinablapsi=zero 286 287 mykpt=.true. 288 mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) 289 if (mykpt) then 290 291 ! Constants depending on k-point 292 istwf_k=dtset%istwfk(ikpt) 293 cplex=2;if (istwf_k>1) cplex=1 294 295 ! Extract cprj for this k-point. 296 nband_cprj_k=nband_k;if (cprj_paral_band) nband_cprj_k=nband_k/mpi_enreg%nproc_band 297 if (mkmem*nsppol/=1) then 298 ABI_DATATYPE_ALLOCATE(cprj_k_loc,(natom,my_nspinor*nband_cprj_k)) 299 call pawcprj_alloc(cprj_k_loc,0,dimcprj) 300 call pawcprj_get(atindx1,cprj_k_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isppol,& 301 & mband_cprj,mkmem,natom,nband_cprj_k,nband_cprj_k,my_nspinor,nsppol,dtfil%unpaw,& 302 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 303 else 304 cprj_k_loc => cprj 305 end if 306 307 ! if cprj are distributed over bands, gather them (because we need to mix bands) 308 if (cprj_paral_band) then 309 ABI_DATATYPE_ALLOCATE(cprj_k,(natom,my_nspinor*nband_k)) 310 call pawcprj_alloc(cprj_k,0,dimcprj) 311 call pawcprj_mpi_allgather(cprj_k_loc,cprj_k,natom,my_nspinor*nband_cprj_k,dimcprj,0,& 312 & mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.false.) 313 else 314 cprj_k => cprj_k_loc 315 end if 316 317 ABI_ALLOCATE(tnm,(2,3,nband_k,nphicor,natom)) 318 tnm=zero 319 320 ! Loops on bands 321 do jb=1,nband_k 322 323 if (mpi_enreg%paral_kgb==1) then 324 if (mod(jb-1,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle 325 elseif (xmpi_paral==1) then 326 if (abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-me_kpt)/=0) cycle 327 end if 328 jbsp=(jb-1)*my_nspinor 329 330 if (cplex==1) then 331 do ispinor=1,my_nspinor 332 jbsp=jbsp+1 333 do iatom=1,natom 334 itypat=dtset%typat(iatom) 335 lmn_size=pawtab(itypat)%lmn_size 336 do jlmn=1,lmn_size 337 do ilmn=1,lmncmax 338 ib=indlmn_core(5,ilmn) 339 cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn) 340 tnm(2,:,jb,ib,iatom)=tnm(2,:,jb,ib,iatom)+cpnm1*phipphj(itypat)%value(:,jlmn,ilmn) 341 end do !ilmn 342 end do !jlmn 343 end do !iatom 344 end do !ispinor 345 else 346 do ispinor=1,my_nspinor 347 jbsp=jbsp+1 348 do iatom=1,natom 349 itypat=dtset%typat(iatom) 350 lmn_size=pawtab(itypat)%lmn_size 351 do jlmn=1,lmn_size 352 do ilmn=1,lmncmax 353 ib=indlmn_core(5,ilmn) 354 cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn) 355 cpnm2=cprj_k(iatom,jbsp)%cp(2,jlmn) 356 tnm(1,:,jb,ib,iatom)=tnm(1,:,jb,ib,iatom)+cpnm1*phipphj(itypat)%value(:,jlmn,ilmn) 357 tnm(2,:,jb,ib,iatom)=tnm(2,:,jb,ib,iatom)+cpnm2*phipphj(itypat)%value(:,jlmn,ilmn) 358 end do !ilmn 359 end do !jlmn 360 end do !iatom 361 end do !ispinor 362 end if 363 364 ! End loops on bands 365 end do ! jb 366 367 ! Reduction in case of parallelism 368 if (mpi_enreg%paral_kgb==1) then 369 call timab(48,1,tsec) 370 call xmpi_sum_master(tnm,0,spaceComm_bandspin,ierr) 371 call timab(48,2,tsec) 372 end if 373 374 psinablapsi(:,:,:,:,:)=psinablapsi(:,:,:,:,:)+tnm(:,:,:,:,:) 375 ABI_DEALLOCATE(tnm) 376 377 if (mkmem/=0) ibg = ibg + my_nspinor*nband_cprj_k 378 379 if (cprj_paral_band) then 380 call pawcprj_free(cprj_k) 381 ABI_DATATYPE_DEALLOCATE(cprj_k) 382 end if 383 if (mkmem*nsppol/=1) then 384 call pawcprj_free(cprj_k_loc) 385 ABI_DATATYPE_DEALLOCATE(cprj_k_loc) 386 end if 387 388 if (me==0) then 389 do iatom=1,natom 390 write(ount) ((psinablapsi(1:2,1,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 391 write(ount) ((psinablapsi(1:2,2,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 392 write(ount) ((psinablapsi(1:2,3,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 393 end do 394 !DEBUG 395 ! do iatom=1,natom 396 ! write(138,*) ((psinablapsi(1:2,1,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 397 ! write(138,*) ((psinablapsi(1:2,2,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 398 ! write(138,*) ((psinablapsi(1:2,3,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 399 ! end do 400 !DEBUG 401 402 elseif (mpi_enreg%me_band==0.and.mpi_enreg%me_fft==0) then 403 call xmpi_exch(psinablapsi,etiq,me_kpt,psinablapsi,0,spaceComm_k,ierr) 404 end if 405 406 elseif (me==0) then 407 sender=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol)) 408 call xmpi_exch(psinablapsi,etiq,sender,psinablapsi,0,spaceComm_k,ierr) 409 do iatom=1,natom 410 write(ount) ((psinablapsi(1:2,1,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 411 write(ount) ((psinablapsi(1:2,2,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 412 write(ount) ((psinablapsi(1:2,3,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 413 end do 414 415 !DEBUG 416 ! do iatom=1,natom 417 ! write(138,*) ((psinablapsi(1:2,1,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 418 ! write(138,*) ((psinablapsi(1:2,2,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 419 ! write(138,*) ((psinablapsi(1:2,3,ib,jb,iatom),ib=1,nband_k),jb=1,nphicor) 420 ! end do 421 !DEBUG 422 423 end if ! mykpt 424 425 bdtot_index=bdtot_index+nband_k 426 427 ! End loop on spin,kpt 428 end do ! ikpt 429 end do !isppol 430 431 !Close file 432 call WffClose(wff1,ierr) 433 434 !Datastructures deallocations 435 do itypat=1,dtset%ntypat 436 ABI_DEALLOCATE(phipphj(itypat)%value) 437 end do 438 ABI_DATATYPE_DEALLOCATE(phipphj) 439 ABI_DEALLOCATE(indlmn_core) 440 ABI_DEALLOCATE(psinablapsi) 441 442 DBG_EXIT("COLL") 443 444 end subroutine optics_paw_core