1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 subroutine transition_rate( ng, psi, ek, efermi, 9 . ekloc, efermiloc, temp, 10 . smooth, wmin, wmax, 11 . Sr, Aux, Aux2, numh, listhptr, listh, 12 . indxuo, no, nuo, nuotot, xij, 13 . maxnh, nbands, kpoint, matrix, 14 . intraband ) 15C ********************************************************************* 16C Finds the matrix element for the dipolar transition between to states 17C Written by DSP. August 1999 18C Restyled for f90 version by JDG. June 2004 19C Modified by DSP. April 2010 20C **************************** INPUT ********************************** 21C integer ng : first dimension of psi, Aux and Aux2 22C real*8 psi(ng,nuotot,nuo) : Wavefunctions in current k point 23C real*8 ek(nuotot) : Eigenvalues 24C real*8 temp : Electronic temperature 25C real*8 efermi : Fermi level 26C real*8 ekloc(nuotot) : Eigenvalues (possibly modified by escissor) 27C real*8 efermiloc : Fermi level (possibly modified by escissor) 28C real*8 wmin : minimum transition energy required 29C real*8 wmax : maximum transition energy required 30C real*8 smooth : artificial width given to transitions 31C real*8 Sr(maxnh) : Position operator matrix elements (sparse) 32C real*8 Aux(ng,nuotot,nuo) : Auxiliary space 33C real*8 Aux2(ng,nuotot,nuo) : Auxiliary space 34C integer numh(nuo) : Number of nonzero elements of each row 35C of hamiltonian matrix 36C integer listhptr(nuo) : Pointer to start of row in listh 37C integer listh(maxnh) : Nonzero hamiltonian-matrix element 38C column indexes for each matrix row 39C integer indxuo(no) : Index of equivalent orbital in unit cell 40C Unit cell orbitals must be the first in 41C orbital lists, i.e. indxuo.le.nuo, with 42C nuo the number of orbitals in unit cell 43C integer no : Number of basis orbitals in supercell 44C integer nuo : Number of orbitals in the cell (locally) 45C integer nuotot : Number of orbitals in the cell (globally) 46C integer maxnh : Maximum dimension of listh 47C real*8 xij(3,maxnh) : Vectors between orbital centers (sparse) 48C (not used if only gamma point) 49C integer nocc : number of occupied states 50C real*8 kpoint(3) : Current kpoint 51C real*8 dk(3) : Vector joining the previous and current 52C kpoint 53C character matrix*1 : 'R' or 'P' for position or momentum operator 54C *************************** OUTPUT **************************** 55C real*8 Aux(2,nuotot,nuo) : matrix elements of the dipolar transition 56C real*8 intraband(2,nuo) : weight of intraband transitions (if a band 57C crosses the Fermi energy 58C *************************** UNITS *********************************** 59C Lengths in atomic units (Bohr). 60C k vectors in reciprocal atomic units. 61C Energies in Rydbergs. 62C ********************************************************************* 63 64 use precision 65 use parallel, only : BlockSize, Node, Nodes 66 use parallelsubs, only : LocalToGlobalOrb 67 use m_fermid, only : stepf 68 use sys, only : die 69 use alloc, only : re_alloc, de_alloc 70#ifdef MPI 71 use mpi_siesta 72#endif 73 74 implicit none 75 76C Passed variables 77 integer nuo, nuotot, maxnh, no, 78 . listh(maxnh), numh(nuo), listhptr(nuo), 79 . indxuo(no), nbands, ng 80 81 real(dp) 82 . xij(3,maxnh), Sr(maxnh), 83 . psi(ng,nuotot,nuo), kpoint(3), Aux(ng,nuotot,nuo), 84 . Aux2(ng,nuotot,nuo), ek(nuotot), temp, efermi, 85 . intraband(ng,nuo), smooth, wmin, wmax, ekloc(nuotot), 86 . efermiloc, ediff 87 88 character 89 . matrix*1 90 91C Internal variables 92 integer 93 . ind, iuo, juo, j, ie, iie, iio, jje, je, jo, 94 . BNodei, Bnodej, BTest, ig 95 real(dp) 96 . kxij, skxij, ckxij, pipj1, pipj2, 97 . f1, f2, intrb(2), intrb2(2) 98 99#ifdef MPI 100 integer :: MPIerror 101 real(dp), pointer :: AuxLocal2(:,:) 102#endif 103 real(dp), pointer :: AuxLocal(:,:) 104 real(dp), pointer :: psibandi(:,:) 105 real(dp), pointer :: psibandj(:,:) 106 real(dp), save :: ediffmin = 1.0d-3 107 real(dp), save :: tiny = 1.0d-9 108 109 110C Start timer 111 call timer('transrate',1) 112 113C Check input matrix 114 if(matrix.ne.'P'.and.matrix.ne.'R') 115 $ call die('transrate: matrix only can take values R or P') 116 117C Initialise matrix elements to zero 118 do iuo = 1,nuo 119 do juo = 1,nuotot 120 do ig=1,ng 121 Aux(ig,juo,iuo) = 0.0_dp 122 Aux2(ig,juo,iuo) = 0.0_dp 123 enddo 124 enddo 125 do ig=1,ng 126 intraband(ig,iuo)=0.0_dp 127 enddo 128 enddo 129 130C Compute matrix elements 131 do iuo = 1,nuo 132 do j = 1,numh(iuo) 133 ind = listhptr(iuo) + j 134 jo = listh(ind) 135 juo = indxuo(jo) 136 if(ng.gt.1) then 137 kxij = kpoint(1) * xij(1,ind) + 138 . kpoint(2) * xij(2,ind) + 139 . kpoint(3) * xij(3,ind) 140 ckxij = dcos(kxij) 141 skxij = dsin(kxij) 142 Aux2(1,juo,iuo) = Aux2(1,juo,iuo) + Sr(ind)*ckxij 143 Aux2(2,juo,iuo) = Aux2(2,juo,iuo) + Sr(ind)*skxij 144 else 145 Aux2(1,juo,iuo) = Aux2(1,juo,iuo) + Sr(ind) 146 endif 147 enddo 148 enddo 149 150C Allocate workspace array 151#ifdef MPI 152 nullify( AuxLocal2 ) 153 call re_alloc( AuxLocal2, 1, ng, 1, nuotot, 154 & name = 'AuxLocal2', 155 & routine = 'transition_rate' ) 156#endif 157 nullify( AuxLocal ) 158 call re_alloc( AuxLocal, 1, ng, 1, nuotot, 159 & name = 'AuxLocal', 160 & routine = 'transition_rate' ) 161 nullify( psibandi ) 162 call re_alloc( psibandi, 1, ng, 1, nuotot, 163 & name = 'psibandi', 164 & routine = 'transition_rate' ) 165 nullify( psibandj ) 166 call re_alloc( psibandj, 1, ng, 1, nuotot, 167 & name = 'psibandj', 168 & routine = 'transition_rate' ) 169 170 BNodei = 0 171 iie = 0 172 do ie = 1,nbands 173 f1 = 2.0d0*stepf((ekloc(ie)-efermiloc)/temp) 174 if (Node.eq.BNodei) then 175 iie = iie + 1 176 do j = 1,nuotot 177 do ig=1,ng 178 psibandi(ig,j) = psi(ig,j,iie) 179 enddo 180 enddo 181 endif 182#ifdef MPI 183 call MPI_Bcast(psibandi(1,1),ng*nuotot,MPI_double_precision, 184 . BNodei,MPI_Comm_World,MPIerror) 185#endif 186 187C Compute the intra-band contribution to include a Drude-like term 188C for metals 189 if(matrix.eq.'P') then 190 intrb(:)=0.0_dp 191 ediff=ek(ie)-efermi 192 if(dabs(ediff).lt.2.0d0*smooth) then 193 do iuo = 1,nuo 194 call LocalToGlobalOrb(iuo,Node,Nodes,iio) 195 do juo = 1,nuotot 196 if(ng.eq.2) then 197 pipj1 = psibandi(1,iio)*psibandi(1,juo) + 198 . psibandi(2,iio)*psibandi(2,juo) 199 pipj2 = psibandi(1,iio)*psibandi(2,juo) - 200 . psibandi(2,iio)*psibandi(1,juo) 201 202C The factor 1/2 is necessary because the in subroutine phirphi_opt 203C we have added a factor of 2, to be canceled by the energy denominator 204C in Ry (rather than in Ha). For the calculation of the 205C Drude term we do not have the energy denominator 206 intrb(1) = intrb(1) 207 . + 0.5_dp * (pipj1*Aux2(1,juo,iuo) 208 . - pipj2*Aux2(2,juo,iuo)) 209 210 intrb(2) = intrb(2) 211 . + 0.5_dp * (pipj1*Aux2(2,juo,iuo) 212 . + pipj2*Aux2(1,juo,iuo)) 213 else 214 pipj1 = psibandi(1,iio)*psibandi(1,juo) 215 216 intrb(1) = intrb(1) 217 . + 0.5_dp* (pipj1*Aux2(1,juo,iuo)) 218 intrb(2) = 0.0_dp 219 endif 220 enddo 221 enddo 222 endif 223#ifdef MPI 224 call MPI_Reduce(intrb,intrb2,2, 225 . MPI_double_precision,MPI_sum,BNodei,MPI_Comm_World,MPIerror) 226 if (Node.eq.BNodei) then 227 intraband(1:ng,iie) = intrb2(1:ng) 228 endif 229#else 230 if (Node.eq.BNodei) then 231 intraband(1:ng,iie) = intrb(1:ng) 232 endif 233#endif 234 endif 235 236 AuxLocal(1:ng,1:nuotot) = 0.0d0 237 238 BNodej = 0 239 jje = 0 240 do je = 1,nbands 241 if (Node.eq.BNodej) then 242 jje = jje + 1 243 endif 244 if (dabs(ek(ie)-ek(je)).gt.ediffmin) then 245 f2 = 2.0d0*stepf((ekloc(je)-efermiloc)/temp) 246 247 if (f1*(2.0d0-f2).gt.tiny) then 248 249 ediff = ekloc(je) - ekloc(ie) 250 if ((ediff.ge.max(wmin-2.0d0*smooth,0.0d0)).and. 251 . (ediff.le.wmax+2.0d0*smooth)) then 252 253 if (Node.eq.BNodej) then 254 do j = 1,nuotot 255 do ig=1,ng 256 psibandj(ig,j) = psi(ig,j,jje) 257 enddo 258 enddo 259 endif 260#ifdef MPI 261 call MPI_Bcast(psibandj(1,1),ng*nuotot, 262 $ MPI_double_precision,BNodej,MPI_Comm_World,MPIerror) 263#endif 264 265 do iuo = 1,nuo 266 call LocalToGlobalOrb(iuo,Node,Nodes,iio) 267 do juo = 1,nuotot 268 if(ng.eq.2) then 269 pipj1 = psibandi(1,iio)*psibandj(1,juo) + 270 . psibandi(2,iio)*psibandj(2,juo) 271 pipj2 = psibandi(1,iio)*psibandj(2,juo) - 272 . psibandi(2,iio)*psibandj(1,juo) 273 274 AuxLocal(1,je) = AuxLocal(1,je) 275 . + pipj1*Aux2(1,juo,iuo) 276 . - pipj2*Aux2(2,juo,iuo) 277 278 AuxLocal(2,je) = AuxLocal(2,je) 279 . + pipj1*Aux2(2,juo,iuo) 280 . + pipj2*Aux2(1,juo,iuo) 281 else 282 pipj1 = psibandi(1,iio)*psibandj(1,juo) 283 284 AuxLocal(1,je) = AuxLocal(1,je) 285 . + pipj1*Aux2(1,juo,iuo) 286 endif 287 288 enddo 289 enddo 290 if (matrix.eq.'P') then 291 do ig=1,ng 292 AuxLocal(ig,je) = AuxLocal(ig,je)/(ek(je)-ek(ie)) 293 enddo 294 endif 295 endif 296 endif 297 endif 298 BTest = je/BlockSize 299 if (BTest*BlockSize.eq.je) then 300 BNodej = BNodej + 1 301 if (BNodej .gt. Nodes-1) BNodej = 0 302 endif 303 enddo 304#ifdef MPI 305 call MPI_Reduce(AuxLocal(1,1),AuxLocal2(1,1),ng*nuotot, 306 . MPI_double_precision,MPI_sum,BNodei,MPI_Comm_World,MPIerror) 307 if (Node.eq.BNodei) then 308 Aux(1:ng,1:nuotot,iie) = AuxLocal2(1:ng,1:nuotot) 309 endif 310#else 311 if (Node.eq.BNodei) then 312 Aux(1:ng,1:nuotot,iie) = AuxLocal(1:ng,1:nuotot) 313 endif 314#endif 315 BTest = ie/BlockSize 316 if (BTest*BlockSize.eq.ie) then 317 BNodei = BNodei + 1 318 if (BNodei .gt. Nodes-1) BNodei = 0 319 endif 320 enddo 321 322C Free workspace array 323 call de_alloc( psibandj, 'psibandj', 'transition_rate') 324 call de_alloc( psibandi, 'psibandi', 'transition_rate') 325 call de_alloc( AuxLocal, 'AuxLocal', 'transition_rate') 326#ifdef MPI 327 call de_alloc( AuxLocal2, 'AuxLocal2', 'transition_rate') 328#endif 329 330C Stop timer 331 call timer('transrate',2) 332 333 end 334