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! 8module m_rhoofdsp 9 10 implicit none 11 12 private 13 public :: rhoofdsp 14 15contains 16 subroutine rhoofdsp( no, nml, nmpl, maxnd, numd, listdptr, & 17 listd, nspin, Dscf, rhoscf, nuo, nuotot, & 18 iaorb, iphorb, isa, q ) 19! ******************************************************************** 20! Finds the SCF density at the mesh points from the density matrix. 21! Written by P.Ordejon and J.M.Soler. May'95. 22! Re-ordered so that mesh is the outer loop and the orbitals are 23! handled as lower-half triangular. J.D.Gale and J.M.Soler, Feb'99 24! Version of rhoofd that optionally uses a direct algorithm to save 25! memory. Modified by J.D.Gale, November'99 26! Spiral version written by V. M. Garcia-Suarez. June 2002. 27! Fixed for parallel runs by Nick Papior, June 2020. 28! *********************** InpUT ************************************** 29! integer no : Number of basis orbitals 30! integer nml : Local numer of mesh divisions in each cell vector 31! integer nmpl : Local numer of mesh divisions 32! integer maxnd : First dimension of listD and Dscf, and 33! maximum number of nonzero elements in 34! any row of Dscf 35! integer numd(nuo) : Number of nonzero elemts in each row of Dscf 36! integer listdptr(nuo) : Pointer to start of rows in listd 37! integer listd(maxnd) : List of nonzero elements in each row of Dscf 38! integer nspin : Number of spin components 39! real*8 Dscf(maxnd) : Rows of Dscf that are non-zero 40! integer nuo : Number of orbitals in unit cell locally 41! integer nuotot : Number of orbitals in unit cell in total 42! integer iaorb(*) : Pointer to atom to which orbital belongs 43! integer iphorb(*) : Orbital index within each atom 44! integer isa(*) : Species index of all atoms 45! real*8 q(3) : Wave vector for spiral configuration. 46! *********************** OUTPUT ************************************** 47! real rhoscf(nsp,nmpl,nspin): SCF density at mesh points 48! ********************************************************************* 49 50! Modules 51 use precision, only: dp, grid_p 52 use parallel, only: Nodes 53 use atmfuncs, only: rcut, all_phi 54 use atm_types, only: nsmax=>nspecies 55 use atomlist, only: indxuo 56 use listsc_module, only: listsc 57 use mesh, only: nsp, dxa, xdop, xdsp, cmesh, nsm, meshLim 58 use meshdscf 59 use meshphi 60 use sys, only : die 61 use alloc, only: re_alloc, de_alloc 62 63! Argument types and dimensions 64 integer :: & 65 no, nml(3), nmpl, nspin, maxnd, nuo, nuotot, iaorb(*), & 66 iphorb(*), isa(*), numd(nuo), listdptr(nuo), listd(maxnd) 67 68 real(grid_p), intent(out) :: rhoscf(nsp,nmpl,nspin) 69 70 real(dp) :: Dscf(maxnd,nspin), q(3) 71 72 external :: memory, timer, ipack 73 74 ! Internal variables and arrays 75 integer, parameter :: minloc = 1000 ! Min buffer size for local copy of Dscf 76 integer, parameter :: maxoa = 100 ! Max # of orbitals per atom 77 78 integer :: triang, & 79 i, ia, ic, ii, ijl, il, imp, ind, ispin, io, ix, & 80 iop, ip, iphi, is, isp, iu, iul, j, jc, jl, & 81 last, lasta, lastop, maxloc, maxloc2, maxndl, nc, nphiloc, & 82 iii(3) 83 84 integer, pointer :: ilc(:) => null() 85 integer, pointer :: ilocal(:) => null() 86 integer, pointer :: iorb(:) => null() 87 real(dp), pointer :: phia(:,:) => null() 88 89 real(dp), pointer :: Clocal(:,:) => null() 90 real(dp), pointer :: Dlocal(:,:) => null() 91 92 logical :: ParallelLocal 93 94 real(dp) :: Cij, Dij(4), & 95 dxsp(3), r2cut(nsmax), r2sp, & 96 xr(3), Rdi, qRdi, cqRdi(nsp), sqRdi(nsp) 97 98 ! Start time counter 99 call timer('rhoofdsp',1) 100 101 ! Find atomic cutoff radiae 102 r2cut(:) = 0.0_dp 103 do i = 1,nuotot 104 ia = iaorb(i) 105 is = isa(ia) 106 io = iphorb(i) 107 r2cut(is) = max( r2cut(is), rcut(is,io)**2 ) 108 end do 109 110 ! Set algorithm logical 111 ParallelLocal = (Nodes > 1) 112 113 ! Find value of maxloc 114 maxloc2 = maxval(endpht(1:nmpl)-endpht(0:nmpl-1)) 115 maxloc = maxloc2 + minloc 116 maxloc = min( maxloc, no ) 117 triang = (maxloc+1)*(maxloc+2)/2 118 if ( ParallelLocal ) then 119 if ( nrowsDscfL > 0 ) then 120 maxndl = listdlptr(nrowsDscfL) + numdl(nrowsDscfL) 121 else 122 maxndl = 1 123 end if 124 end if 125 126 ! Allocate local memory 127 nullify(ilocal, ilc, iorb, phia, Clocal, Dlocal, DscfL) 128 129 call re_alloc( ilocal, 1, no, 'ilocal', 'rhoofdsp' ) 130 call re_alloc( ilc, 1, maxloc2, 'ilc', 'rhoofdsp' ) 131 call re_alloc( iorb, 0, maxloc, 'iorb', 'rhoofdsp' ) 132 133 call re_alloc( Dlocal, 1, nspin, 1, triang, 'Dlocal', 'rhoofdsp' ) 134 call re_alloc( Clocal, 1, nsp, 1, maxloc2, 'Clocal', 'rhoofdsp' ) 135 if ( DirectPhi ) then 136 call re_alloc(phia, 1, maxoa, 1, nsp, "phia", "rhoofdsp") 137 end if 138 139 if ( ParallelLocal ) then 140 call re_alloc( DscfL, 1, maxndl, 1, nspin, 'DscfL', 'rhoofdsp' ) 141 142 ! Redistribute Dscf to DscfL form 143 call matrixOtoM( maxnd, numd, listdptr, maxndl, nuo, & 144 nspin, Dscf, DscfL ) 145 endif 146 147 ! Initializations 148 rhoscf(:,:,:) = 0.0_grid_p 149 Dlocal(:,:) = 0.0_dp 150 ilocal(:) = 0 151 iorb(:) = 0 152 last = 0 153 154 ! Loop over grid points 155 do ip = 1,nmpl 156 157 ! Find point coordinates 158 call ipack(-1,3,nml,iii,ip) 159 ! Correct for processor offset 160 iii(:) = iii(:) + meshLim(1,:) - 1 161 do ix = 1,3 162 xr(ix) = iii(1)*cmesh(ix,1) + iii(2)*cmesh(ix,2) + iii(3)*cmesh(ix,3) 163 end do 164 165 ! Find number of nonzero orbitals at this point 166 nc = endpht(ip) - endpht(ip-1) 167 168 ! iorb(il)>0 means that row il of Dlocal must not be overwritten 169 ! iorb(il)=0 means that row il of Dlocal is empty 170 ! iorb(il)<0 means that row il of Dlocal contains a valid row of 171 ! Dscf, but which is not required at this point 172 do ic = 1,nc 173 imp = endpht(ip-1) + ic 174 i = lstpht(imp) 175 il = ilocal(i) 176 if ( il > 0 ) iorb(il) = i 177 end do 178 179 ! Look for required rows of Dscf not yet stored in Dlocal 180 do ic = 1,nc 181 imp = endpht(ip-1) + ic 182 i = lstpht(imp) 183 if ( ilocal(i) == 0 ) then 184 185 ! Look for an available row in Dlocal 186 do il = 1,maxloc 187 ! last runs circularly over rows of Dlocal 188 last = last + 1 189 if (last > maxloc) last = 1 190 if (iorb(last) <= 0) goto 10 191 enddo 192 call die('rhoofdsp: no slot available in Dlocal') 19310 continue 194 195 ! Copy row i of Dscf into row last of Dlocal 196 j = abs(iorb(last)) 197 if ( j /= 0 ) ilocal(j) = 0 198 ilocal(i) = last 199 iorb(last) = i 200 il = last 201 iu = indxuo(i) 202 if (ParallelLocal) then 203 iul = NeedDscfL(iu) 204 if (i .eq. iu) then 205 do ii = 1, numdl(iul) 206 ind = listdlptr(iul)+ii 207 j = listdl(ind) 208 jl = ilocal(j) 209 ijl = idx_ijl(il,jl) 210 do ispin = 1,nspin 211 Dlocal(ispin,ijl) = DscfL(ind,ispin) 212 enddo 213 enddo 214 else 215 do ii = 1, numdl(iul) 216 ind = listdlptr(iul)+ii 217 j = listsc( i, iu, listdl(ind) ) 218 jl = ilocal(j) 219 ijl = idx_ijl(il,jl) 220 do ispin = 1,nspin 221 Dlocal(ispin,ijl) = DscfL(ind,ispin) 222 enddo 223 enddo 224 endif 225 else 226 if (i .eq. iu) then 227 do ii = 1, numd(iu) 228 ind = listdptr(iu)+ii 229 j = listd(ind) 230 jl = ilocal(j) 231 ijl = idx_ijl(il,jl) 232 do ispin = 1,nspin 233 Dlocal(ispin,ijl) = Dscf(ind,ispin) 234 enddo 235 enddo 236 else 237 do ii = 1, numd(iu) 238 ind = listdptr(iu)+ii 239 j = listsc( i, iu, listd(ind) ) 240 jl = ilocal(j) 241 ijl = idx_ijl(il,jl) 242 do ispin = 1,nspin 243 Dlocal(ispin,ijl) = Dscf(ind,ispin) 244 enddo 245 enddo 246 endif 247 endif 248 endif 249 enddo 250 251 lasta = 0 252 lastop = 0 253 254 ! Generate or retrieve phi values for all orbitals up to nc 255 do ic = 1,nc 256 imp = endpht(ip-1) + ic 257 i = lstpht(imp) 258 il = ilocal(i) 259 iu = indxuo(i) 260 ia = iaorb(i) 261 is = isa(ia) 262 iop = listp2(imp) 263 ilc(ic) = il 264 265 if (DirectPhi) then 266 if (ia.ne.lasta .or. iop.ne.lastop) then 267 lasta = ia 268 lastop = iop 269 do isp = 1,nsp 270 dxsp(:) = xdsp(:,isp) + xdop(:,iop) - dxa(:,ia) 271 r2sp = sum(dxsp**2) 272 if (r2sp.lt.r2cut(is)) then 273 call all_phi( is, +1, dxsp, maxoa, nphiloc, phia(:,isp) ) 274 else 275 phia(:,isp) = 0.0_dp 276 endif 277 enddo 278 endif 279 iphi = iphorb(i) 280 Clocal(:,ic) = phia(iphi,:) 281 else 282 Clocal(:,ic) = phi(:,imp) 283 endif 284 285 enddo 286 287 ! Loop on first orbital of mesh point 288 do ic = 1,nc 289 imp = endpht(ip-1) + ic 290 i = lstpht(imp) 291 il = ilocal(i) 292 iu = indxuo(i) 293 ia = iaorb(i) 294 is = isa(ia) 295 iop = listp2(imp) 296 297 ! Calculate spiral phase 298 do isp = 1, nsp 299 qRdi = 0._dp 300 do ix = 1, 3 301 Rdi = xr(ix) + xdsp(ix,isp) - xdop(ix,iop) + dxa(ix,ia) 302 qRdi = qRdi + q(ix) * Rdi 303 end do 304 cqRdi(isp) = cos(qRdi) 305 sqRdi(isp) = sin(qRdi) 306 end do 307 308 ! Loop on second orbital of mesh point, here only jc < ic 309 do jc = 1, ic - 1 310 jl = ilc(jc) 311 312 ijl = idx_ijl(il,jl) 313 314 Dij(1) = 2*Dlocal(1,ijl) 315 Dij(2) = 2*Dlocal(2,ijl) 316 Dij(3) = Dlocal(3,ijl) 317 Dij(4) = Dlocal(4,ijl) 318 319 ! Loop over sub-points 320 do isp = 1, nsp 321 Cij = Clocal(isp,ic) * Clocal(isp,jc) 322 rhoscf(isp,ip,1) = rhoscf(isp,ip,1) + Dij(1)*Cij 323 rhoscf(isp,ip,2) = rhoscf(isp,ip,2) + Dij(2)*Cij 324 rhoscf(isp,ip,3) = rhoscf(isp,ip,3) + & 325 (Dij(3)*cqRdi(isp) + Dij(4)*sqRdi(isp))*Cij 326 rhoscf(isp,ip,4) = rhoscf(isp,ip,4) + & 327 (Dij(4)*cqRdi(isp) - Dij(3)*sqRdi(isp))*Cij 328 end do 329 330 enddo 331 332 ! onsite orbital 333 ijl = idx_ijl(il,il) 334 Dij(1) = Dlocal(1,ijl) 335 Dij(2) = Dlocal(2,ijl) 336 Dij(3) = Dlocal(3,ijl) 337 Dij(4) = Dlocal(4,ijl) 338 do isp = 1, nsp 339 Cij = Clocal(isp,ic) * Clocal(isp,ic) 340 rhoscf(isp,ip,1) = rhoscf(isp,ip,1) + Dij(1)*Cij 341 rhoscf(isp,ip,2) = rhoscf(isp,ip,2) + Dij(2)*Cij 342 rhoscf(isp,ip,3) = rhoscf(isp,ip,3) + & 343 (Dij(3)*cqRdi(isp) + Dij(4)*sqRdi(isp))*Cij 344 rhoscf(isp,ip,4) = rhoscf(isp,ip,4) + & 345 (Dij(4)*cqRdi(isp) - Dij(3)*sqRdi(isp))*Cij 346 end do 347 348 ! Loop on second orbital of mesh point (here only jc > ic) 349 do jc = ic + 1, nc 350 jl = ilc(jc) 351 ijl = idx_ijl(il,jl) 352 Dij(3) = Dlocal(3,ijl) 353 Dij(4) = Dlocal(4,ijl) 354 do isp = 1, nsp 355 Cij = Clocal(isp,ic) * Clocal(isp,jc) 356 rhoscf(isp,ip,3) = rhoscf(isp,ip,3) + & 357 (Dij(3)*cqRdi(isp) + Dij(4)*sqRdi(isp))*Cij 358 rhoscf(isp,ip,4) = rhoscf(isp,ip,4) + & 359 (Dij(4)*cqRdi(isp) - Dij(3)*sqRdi(isp))*Cij 360 end do 361 end do 362 363 enddo 364 365 ! Restore iorb for next point 366 do imp = 1+endpht(ip-1), endpht(ip) 367 i = lstpht(imp) 368 il = ilocal(i) 369 iorb(il) = -i 370 enddo 371 372 enddo 373 374 ! Free local memory 375 call de_alloc( ilocal, 'ilocal', 'rhoofdsp' ) 376 call de_alloc( ilc,'ilc', 'rhoofdsp' ) 377 call de_alloc( iorb,'iorb', 'rhoofdsp' ) 378 call de_alloc( Clocal, 'Clocal', 'rhoofdsp' ) 379 call de_alloc( Dlocal, 'Dlocal', 'rhoofdsp' ) 380 if ( DirectPhi ) then 381 call de_alloc( phia, 'phia', 'vmatsp' ) 382 end if 383 384 if (ParallelLocal) then 385 call de_alloc( DscfL, 'DscfL', 'rhoofdsp' ) 386 endif 387 388 call timer('rhoofdsp',2) 389 390 contains 391 392 ! In any case will the compiler most likely inline this 393 ! small routine. So it should not pose any problem. 394 pure function idx_ijl(i,j) result(ij) 395 integer, intent(in) :: i,j 396 integer :: ij 397 if ( i > j ) then 398 ij = i * (i + 1)/2 + j + 1 399 else 400 ij = j * (j + 1)/2 + i + 1 401 end if 402 end function idx_ijl 403 404 end subroutine rhoofdsp 405 406end module m_rhoofdsp 407