TABLE OF CONTENTS
ABINIT/m_optics_vloc [ Modules ]
NAME
m_optics_vloc
FUNCTION
COPYRIGHT
Copyright (C) 2010-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_optics_vloc 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_wffile 28 use m_xmpi 29 use m_hdr 30 use m_dtset 31 use m_dtfil 32 33 use defs_abitypes, only : MPI_type 34 use m_time, only : timab 35 use m_io_tools, only : get_unit 36 use m_mpinfo, only : proc_distrb_cycle 37 38 implicit none 39 40 private
ABINIT/optics_vloc [ Functions ]
NAME
optics_vloc
FUNCTION
Compute matrix elements need for optical conductivity in a LOCAL potential and store them in a file. Matrix elements = <Phi_i|Nabla|Phi_j>
INPUTS
cg(2,mcg)=planewave coefficients of wavefunctions. dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset fildata= name of the output file 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 mkmem =number of k points treated by this node. mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw. nkpt=number of k points. npwarr(nkpt)=number of planewaves in basis at this k point nsppol=1 for unpolarized, 2 for spin-polarized
OUTPUT
(only writing in a file)
SIDE EFFECTS
NOTES
SOURCE
83 subroutine optics_vloc(cg,dtfil,dtset,eigen0,gprimd,hdr,kg,mband,mcg,mkmem,mpi_enreg,mpw,& 84 & nkpt,npwarr,nsppol) 85 !Arguments ------------------------------------ 86 !scalars 87 integer,intent(in) :: mband,mcg,mkmem,mpw,nkpt,nsppol 88 type(datafiles_type),intent(in) :: dtfil 89 type(dataset_type),intent(in) :: dtset 90 type(MPI_type),intent(in) :: mpi_enreg 91 type(hdr_type),intent(inout) :: hdr 92 !arrays 93 integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt) 94 real(dp),intent(in) :: gprimd(3,3) 95 real(dp),intent(in) :: eigen0(mband*nkpt*nsppol) 96 real(dp),intent(inout) :: cg(2,mcg) 97 98 !Local variables------------------------------- 99 !scalars 100 integer :: iomode,bdtot_index,cplex,etiq,fformopt,ib,icg,ierr,ikg,ikpt 101 integer :: ipw,isppol,istwf_k,iwavef,jb,jwavef 102 integer :: me,me_kpt,my_nspinor,nband_k,npw_k,sender,ount,pnp_size 103 integer :: spaceComm_band,spaceComm_bandfftspin,spaceComm_fft,spaceComm_k 104 logical :: mykpt 105 real(dp) :: cgnm1,cgnm2 106 character(len=500) :: msg 107 !arrays 108 integer :: tmp_shape(3) 109 integer,allocatable :: kg_k(:,:) 110 real(dp) :: kpoint(3),tsec(2) 111 real(dp),allocatable :: eig0_k(:),kpg_k(:,:) 112 real(dp),allocatable :: psinablapsi(:,:,:,:),tnm(:,:,:,:) 113 type(wffile_type) :: wff1 114 ! ************************************************************************ 115 116 DBG_ENTER("COLL") 117 118 !Compatibility tests 119 120 !---------------------------------------------------------------------------------- 121 !2- Computation of <psi_n|-i.nabla|psi_m> for each k 122 !---------------------------------------------------------------------------------- 123 124 !Init parallelism 125 if (mpi_enreg%paral_kgb==1) then 126 spaceComm_k=mpi_enreg%comm_kpt 127 spaceComm_fft=mpi_enreg%comm_fft 128 spaceComm_band=mpi_enreg%comm_band 129 spaceComm_bandfftspin=mpi_enreg%comm_bandspinorfft 130 me_kpt=mpi_enreg%me_kpt 131 else 132 spaceComm_band=0;spaceComm_fft=0;spaceComm_bandfftspin=0 133 spaceComm_k=mpi_enreg%comm_cell 134 me=xmpi_comm_rank(spaceComm_k) 135 me_kpt=me 136 end if 137 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 138 139 !Initialize main variables 140 ABI_MALLOC(psinablapsi,(2,3,mband,mband)) 141 pnp_size=size(psinablapsi) 142 psinablapsi=zero 143 144 iomode= IO_MODE_FORTRAN_MASTER 145 fformopt=612 146 ount = get_unit() 147 call WffOpen(iomode,spaceComm_k,dtfil%fnameabo_app_opt,ierr,wff1,0,me,ount) 148 call hdr_io(fformopt,hdr,2,wff1) 149 !if (me == 0) then 150 !call hdr_fort_write(hdr, wff1%unwff, fformopt,ierr) 151 !ABI_CHECK(ierr /= 0, "hdr_fort_write returned ierr = 0") 152 !end if 153 154 !LOOP OVER SPINS 155 icg=0 156 do isppol=1,nsppol 157 158 ! LOOP OVER k POINTS 159 ikg=0 160 do ikpt=1,nkpt 161 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 162 etiq=ikpt+(isppol-1)*nkpt 163 if (me==0) then 164 ABI_MALLOC(eig0_k,(nband_k)) 165 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 166 end if 167 168 ! Select my k-points 169 mykpt=.true. 170 mykpt=(.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt))) 171 if (mykpt) then 172 173 ! Allocations depending on k-point 174 kpoint(:)=dtset%kptns(:,ikpt) 175 istwf_k=dtset%istwfk(ikpt) 176 npw_k=npwarr(ikpt) 177 cplex=2;if (istwf_k>1) cplex=1 178 ABI_MALLOC(kg_k,(3,npw_k)) 179 ABI_MALLOC(kpg_k,(npw_k*dtset%nspinor,3)) 180 181 ! Get G-vectors for this k-point 182 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 183 ikg=ikg+npw_k 184 185 ! Calculation of k+G in cartesian coordinates 186 do ipw=1,npw_k 187 kpg_k(ipw,1)=(kpoint(1)+kg_k(1,ipw))*gprimd(1,1)& 188 & +(kpoint(2)+kg_k(2,ipw))*gprimd(1,2)& 189 & +(kpoint(3)+kg_k(3,ipw))*gprimd(1,3) 190 kpg_k(ipw,2)=(kpoint(1)+kg_k(1,ipw))*gprimd(2,1)& 191 & +(kpoint(2)+kg_k(2,ipw))*gprimd(2,2)& 192 & +(kpoint(3)+kg_k(3,ipw))*gprimd(2,3) 193 kpg_k(ipw,3)=(kpoint(1)+kg_k(1,ipw))*gprimd(3,1)& 194 & +(kpoint(2)+kg_k(2,ipw))*gprimd(3,2)& 195 & +(kpoint(3)+kg_k(3,ipw))*gprimd(3,3) 196 end do !ipw 197 kpg_k=two_pi*kpg_k 198 if (dtset%nspinor==2) kpg_k(npw_k+1:2*npw_k,1:3)=kpg_k(1:npw_k,1:3) 199 ABI_FREE(kg_k) 200 201 ! 2-A Computation of <psi_tild_n|-i.nabla|psi_tild_m> 202 ! ---------------------------------------------------------------------------------- 203 ! Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates 204 205 ABI_MALLOC(tnm,(2,3,nband_k,nband_k)) 206 tnm=zero 207 208 ! Loops on bands 209 do jb=1,nband_k 210 jwavef=(jb-1)*npw_k*my_nspinor+icg 211 if (mpi_enreg%paral_kgb/=1) then 212 tmp_shape = shape(mpi_enreg%proc_distrb) 213 if (ikpt > tmp_shape(1)) then 214 msg=' ikpt out of bounds ' 215 ABI_BUG(msg) 216 end if 217 if (abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-me_kpt)/=0) cycle 218 end if 219 do ib=1,jb 220 iwavef=(ib-1)*npw_k*my_nspinor+icg 221 222 ! Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates 223 if (cplex==1) then 224 do ipw=1,npw_k*my_nspinor 225 cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef) 226 tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3) 227 end do 228 else 229 do ipw=1,npw_k*my_nspinor 230 cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)+cg(2,ipw+iwavef)*cg(2,ipw+jwavef) 231 cgnm2=cg(1,ipw+iwavef)*cg(2,ipw+jwavef)-cg(2,ipw+iwavef)*cg(1,ipw+jwavef) 232 tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3) 233 tnm(2,1:3,ib,jb)=tnm(2,1:3,ib,jb)+cgnm2*kpg_k(ipw,1:3) 234 end do 235 end if 236 237 ! Second half of the (n,m) matrix 238 if (ib/=jb) then 239 tnm(1,1:3,jb,ib)= tnm(1,1:3,ib,jb) 240 tnm(2,1:3,jb,ib)=-tnm(2,1:3,ib,jb) 241 end if 242 243 end do ! ib 244 end do ! jb 245 246 ! Reduction in case of parallelism 247 if (mpi_enreg%paral_kgb == 1) then 248 call timab(48,1,tsec) 249 call xmpi_sum_master(tnm,0,spaceComm_bandfftspin,ierr) 250 call timab(48,2,tsec) 251 end if 252 253 psinablapsi(:,:,:,:)=tnm(:,:,:,:) 254 255 ABI_FREE(tnm) 256 257 if (mkmem/=0) then 258 icg = icg + npw_k*my_nspinor*nband_k 259 end if 260 261 ABI_FREE(kpg_k) 262 263 if (me==0) then 264 write(ount)(eig0_k(ib),ib=1,nband_k) 265 write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k) 266 write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k) 267 write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k) 268 elseif (mpi_enreg%me_band==0.and.mpi_enreg%me_fft==0) then 269 call xmpi_exch(psinablapsi,pnp_size,me_kpt,psinablapsi,0,spaceComm_k,etiq,ierr) 270 end if 271 272 elseif (me==0) then 273 sender=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol)) 274 call xmpi_exch(psinablapsi,pnp_size,sender,psinablapsi,0,spaceComm_k,etiq,ierr) 275 write(ount)(eig0_k(ib),ib=1,nband_k) 276 write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k) 277 write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k) 278 write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k) 279 end if ! mykpt 280 281 bdtot_index=bdtot_index+nband_k 282 if (me==0) then 283 ABI_FREE(eig0_k) 284 end if 285 ! End loop on spin,kpt 286 end do ! ikpt 287 end do !isppol 288 289 !Close file 290 call WffClose(wff1,ierr) 291 292 !Datastructures deallocations 293 ABI_FREE(psinablapsi) 294 295 DBG_EXIT("COLL") 296 297 end subroutine optics_vloc