TABLE OF CONTENTS
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>
COPYRIGHT
Copyright (C) 2010-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
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
PARENTS
outscfcv
CHILDREN
hdr_io,timab,wffclose,wffopen,xmpi_exch,xmpi_sum_master
SOURCE
48 #if defined HAVE_CONFIG_H 49 #include "config.h" 50 #endif 51 52 #include "abi_common.h" 53 54 subroutine optics_vloc(cg,dtfil,dtset,eigen0,gprimd,hdr,kg,mband,mcg,mkmem,mpi_enreg,mpw,& 55 & nkpt,npwarr,nsppol) 56 57 use defs_basis 58 use defs_abitypes 59 use m_profiling_abi 60 use m_errors 61 use m_wffile 62 use m_xmpi 63 use m_hdr 64 65 use m_io_tools, only : get_unit 66 67 !This section has been created automatically by the script Abilint (TD). 68 !Do not modify the following lines by hand. 69 #undef ABI_FUNC 70 #define ABI_FUNC 'optics_vloc' 71 use interfaces_18_timing 72 use interfaces_32_util 73 !End of the abilint section 74 75 implicit none 76 77 !Arguments ------------------------------------ 78 !scalars 79 integer,intent(in) :: mband,mcg,mkmem,mpw,nkpt,nsppol 80 type(datafiles_type),intent(in) :: dtfil 81 type(dataset_type),intent(in) :: dtset 82 type(MPI_type),intent(in) :: mpi_enreg 83 type(hdr_type),intent(inout) :: hdr 84 !arrays 85 integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt) 86 real(dp),intent(in) :: gprimd(3,3) 87 real(dp),intent(in) :: eigen0(mband*nkpt*nsppol) 88 real(dp),intent(inout) :: cg(2,mcg) 89 90 !Local variables------------------------------- 91 !scalars 92 integer :: iomode,bdtot_index,cplex,etiq,fformopt,ib,icg,ierr,ikg,ikpt 93 integer :: ipw,isppol,istwf_k,iwavef,jb,jwavef 94 integer :: me,me_kpt,my_nspinor,nband_k,npw_k,sender,ount 95 integer :: spaceComm_band,spaceComm_bandfftspin,spaceComm_fft,spaceComm_k 96 logical :: mykpt 97 real(dp) :: cgnm1,cgnm2 98 character(len=500) :: msg 99 !arrays 100 integer :: tmp_shape(3) 101 integer,allocatable :: kg_k(:,:) 102 real(dp) :: kpoint(3),tsec(2) 103 real(dp),allocatable :: eig0_k(:),kpg_k(:,:) 104 real(dp),allocatable :: psinablapsi(:,:,:,:),tnm(:,:,:,:) 105 type(wffile_type) :: wff1 106 ! ************************************************************************ 107 108 DBG_ENTER("COLL") 109 110 !Compatibility tests 111 112 !---------------------------------------------------------------------------------- 113 !2- Computation of <psi_n|-i.nabla|psi_m> for each k 114 !---------------------------------------------------------------------------------- 115 116 !Init parallelism 117 if (mpi_enreg%paral_kgb==1) then 118 spaceComm_k=mpi_enreg%comm_kpt 119 spaceComm_fft=mpi_enreg%comm_fft 120 spaceComm_band=mpi_enreg%comm_band 121 spaceComm_bandfftspin=mpi_enreg%comm_bandspinorfft 122 me_kpt=mpi_enreg%me_kpt 123 else 124 spaceComm_band=0;spaceComm_fft=0;spaceComm_bandfftspin=0 125 spaceComm_k=mpi_enreg%comm_cell 126 me=xmpi_comm_rank(spaceComm_k) 127 me_kpt=me 128 end if 129 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 130 131 !Initialize main variables 132 ABI_ALLOCATE(psinablapsi,(2,3,mband,mband)) 133 psinablapsi=zero 134 135 iomode= IO_MODE_FORTRAN_MASTER 136 fformopt=612 137 ount = get_unit() 138 call WffOpen(iomode,spaceComm_k,dtfil%fnameabo_app_opt,ierr,wff1,0,me,ount) 139 call hdr_io(fformopt,hdr,2,wff1) 140 !if (me == 0) then 141 !call hdr_fort_write(hdr, wff1%unwff, fformopt,ierr) 142 !ABI_CHECK(ierr /= 0, "hdr_fort_write returned ierr = 0") 143 !end if 144 145 !LOOP OVER SPINS 146 icg=0 147 do isppol=1,nsppol 148 149 ! LOOP OVER k POINTS 150 ikg=0 151 do ikpt=1,nkpt 152 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 153 etiq=ikpt+(isppol-1)*nkpt 154 if (me==0) then 155 ABI_ALLOCATE(eig0_k,(nband_k)) 156 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 157 end if 158 159 ! Select my k-points 160 mykpt=.true. 161 mykpt=(.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt))) 162 if (mykpt) then 163 164 ! Allocations depending on k-point 165 kpoint(:)=dtset%kptns(:,ikpt) 166 istwf_k=dtset%istwfk(ikpt) 167 npw_k=npwarr(ikpt) 168 cplex=2;if (istwf_k>1) cplex=1 169 ABI_ALLOCATE(kg_k,(3,npw_k)) 170 ABI_ALLOCATE(kpg_k,(npw_k*dtset%nspinor,3)) 171 172 ! Get G-vectors for this k-point 173 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 174 ikg=ikg+npw_k 175 176 ! Calculation of k+G in cartesian coordinates 177 do ipw=1,npw_k 178 kpg_k(ipw,1)=(kpoint(1)+kg_k(1,ipw))*gprimd(1,1)& 179 & +(kpoint(2)+kg_k(2,ipw))*gprimd(1,2)& 180 & +(kpoint(3)+kg_k(3,ipw))*gprimd(1,3) 181 kpg_k(ipw,2)=(kpoint(1)+kg_k(1,ipw))*gprimd(2,1)& 182 & +(kpoint(2)+kg_k(2,ipw))*gprimd(2,2)& 183 & +(kpoint(3)+kg_k(3,ipw))*gprimd(2,3) 184 kpg_k(ipw,3)=(kpoint(1)+kg_k(1,ipw))*gprimd(3,1)& 185 & +(kpoint(2)+kg_k(2,ipw))*gprimd(3,2)& 186 & +(kpoint(3)+kg_k(3,ipw))*gprimd(3,3) 187 end do !ipw 188 kpg_k=two_pi*kpg_k 189 if (dtset%nspinor==2) kpg_k(npw_k+1:2*npw_k,1:3)=kpg_k(1:npw_k,1:3) 190 ABI_DEALLOCATE(kg_k) 191 192 ! 2-A Computation of <psi_tild_n|-i.nabla|psi_tild_m> 193 ! ---------------------------------------------------------------------------------- 194 ! Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates 195 196 ABI_ALLOCATE(tnm,(2,3,nband_k,nband_k)) 197 tnm=zero 198 199 ! Loops on bands 200 do jb=1,nband_k 201 jwavef=(jb-1)*npw_k*my_nspinor+icg 202 if (mpi_enreg%paral_kgb/=1) then 203 tmp_shape = shape(mpi_enreg%proc_distrb) 204 if (ikpt > tmp_shape(1)) then 205 msg=' ikpt out of bounds ' 206 MSG_BUG(msg) 207 end if 208 if (abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-me_kpt)/=0) cycle 209 end if 210 do ib=1,jb 211 iwavef=(ib-1)*npw_k*my_nspinor+icg 212 213 ! Computation of (C_nk^*)*C_mk*(k+g) in cartesian coordinates 214 if (cplex==1) then 215 do ipw=1,npw_k*my_nspinor 216 cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef) 217 tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3) 218 end do 219 else 220 do ipw=1,npw_k*my_nspinor 221 cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)+cg(2,ipw+iwavef)*cg(2,ipw+jwavef) 222 cgnm2=cg(1,ipw+iwavef)*cg(2,ipw+jwavef)-cg(2,ipw+iwavef)*cg(1,ipw+jwavef) 223 tnm(1,1:3,ib,jb)=tnm(1,1:3,ib,jb)+cgnm1*kpg_k(ipw,1:3) 224 tnm(2,1:3,ib,jb)=tnm(2,1:3,ib,jb)+cgnm2*kpg_k(ipw,1:3) 225 end do 226 end if 227 228 ! Second half of the (n,m) matrix 229 if (ib/=jb) then 230 tnm(1,1:3,jb,ib)= tnm(1,1:3,ib,jb) 231 tnm(2,1:3,jb,ib)=-tnm(2,1:3,ib,jb) 232 end if 233 234 end do ! ib 235 end do ! jb 236 237 ! Reduction in case of parallelism 238 if (mpi_enreg%paral_kgb == 1) then 239 call timab(48,1,tsec) 240 call xmpi_sum_master(tnm,0,spaceComm_bandfftspin,ierr) 241 call timab(48,2,tsec) 242 end if 243 244 psinablapsi(:,:,:,:)=tnm(:,:,:,:) 245 246 ABI_DEALLOCATE(tnm) 247 248 if (mkmem/=0) then 249 icg = icg + npw_k*my_nspinor*nband_k 250 end if 251 252 ABI_DEALLOCATE(kpg_k) 253 254 if (me==0) then 255 write(ount)(eig0_k(ib),ib=1,nband_k) 256 write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k) 257 write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k) 258 write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k) 259 elseif (mpi_enreg%me_band==0.and.mpi_enreg%me_fft==0) then 260 call xmpi_exch(psinablapsi,etiq,me_kpt,psinablapsi,0,spaceComm_k,ierr) 261 end if 262 263 elseif (me==0) then 264 sender=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol)) 265 call xmpi_exch(psinablapsi,etiq,sender,psinablapsi,0,spaceComm_k,ierr) 266 write(ount)(eig0_k(ib),ib=1,nband_k) 267 write(ount)((psinablapsi(1:2,1,ib,jb),ib=1,nband_k),jb=1,nband_k) 268 write(ount)((psinablapsi(1:2,2,ib,jb),ib=1,nband_k),jb=1,nband_k) 269 write(ount)((psinablapsi(1:2,3,ib,jb),ib=1,nband_k),jb=1,nband_k) 270 end if ! mykpt 271 272 bdtot_index=bdtot_index+nband_k 273 if (me==0) then 274 ABI_DEALLOCATE(eig0_k) 275 end if 276 ! End loop on spin,kpt 277 end do ! ikpt 278 end do !isppol 279 280 !Close file 281 call WffClose(wff1,ierr) 282 283 !Datastructures deallocations 284 ABI_DEALLOCATE(psinablapsi) 285 286 DBG_EXIT("COLL") 287 288 end subroutine optics_vloc