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!------------------------------------------------------------- 9! 10 module meshsubs 11 12 implicit none 13 14 public :: initMesh ! initialises the mesh 15 public :: initAtomMesh ! " atom info on the mesh 16 public :: sameMeshAndAtoms ! have the mesh or atoms changed? 17 public :: PartialCoreOnMesh ! partial core density on the mesh 18 public :: NeutralAtomOnMesh ! neutral atom potential on mesh 19 public :: LocalChargeOnMesh ! local (positive ion) charge 20 public :: ModelCoreChargeOnMesh ! model core charge for Bader 21 public :: PhiOnMesh ! Orbital values on the mesh 22 public :: setupExtMesh ! Sets up extended mesh 23 public :: distriPhiOnMesh ! Sets up data distribution 24 25 private 26 27 CONTAINS 28 29 subroutine InitMesh( na, cell, norb, iaorb, iphorb, isa, 30 & rmax, G2max, G2mesh, nsc, nmpl, 31 & nm, nml, ntm, ntml, ntpl, dvol) 32 33! Modules 34 use precision, only : dp, i8b 35 use parallel, only : Node, Nodes 36 use parallelsubs, only : HowManyMeshPerNode, GlobalToLocalMesh 37 use moreMeshSubs, only : initMeshDistr, setMeshDistr 38 use moreMeshSubs, only : UNIFORM 39 use alloc, only : re_alloc, de_alloc 40 use mesh, only : ne, mop, nmeshg, nmsc, nsm, nsp 41 use mesh, only : cmesh, rcmesh, xdsp, nmuc, nusc 42 use siesta_cml 43 use cellsubs, only : reclat ! Finds reciprocal unit cell 44 use cellsubs, only : volcel ! Finds unit cell volume 45#ifdef MPI 46 use mpi_siesta 47#endif 48 49! Passed arguments 50 integer, intent(in):: na ! Number of atoms in supercell 51 real(dp),intent(in):: cell(3,3) ! Unit cell vectors 52 integer, intent(in):: norb ! Total number of basis 53 ! orbitals in supercell 54 integer, intent(in):: iaorb(norb) ! Atom to which orbitals belong 55 integer, intent(in):: iphorb(norb) ! Orbital index (within atom) 56 ! of each orbital 57 integer, intent(in):: isa(na) ! Species index of each atom 58 real(dp),intent(in):: rmax ! Maximum orbital radius 59 integer, intent(in):: nsc(3) ! Number of unit-cells in each 60 ! supercell direction 61 real(dp),intent(inout):: G2max ! Effective planewave cutoff 62 ! On input : Value required 63 ! On output: Value used, which 64 ! may be larger 65 real(dp),intent(out):: G2mesh ! Planewave cutoff of mesh used 66 ! (same as G2max on output) 67 integer, intent(out):: nmpl ! Number of mesh points stored 68 ! by my processor 69 integer, intent(out):: nm(3) ! Number of mesh divisions of 70 ! each unit cell vector 71 integer, intent(out):: nml(3) ! Sizes of my processor's box 72 ! of mesh points 73 integer, intent(out):: ntm(3) ! Mesh divisions of each unit 74 ! cell vector, incl. subpoints 75 integer, intent(out):: ntml(3) ! Sizes of my processor's mesh 76 ! box, including subpoints 77 integer, intent(out):: ntpl ! Mesh points stored by my 78 ! processor, incl. subpoints 79!!OLD integer, intent(out):: ntopl ! Total number of nonzero 80 ! orbital points stored by my 81 ! processor 82 real(dp),intent(out):: dvol ! Volume per mesh (super)point 83 84C ---------------------------------------------------------------------- 85C Internal variables and arrays: 86C ---------------------------------------------------------------------- 87C real*8 dx(3) : Vector from atom to mesh sub-point 88C real*8 dxp(3) : Vector from atom to mesh point 89C integer i : General-purpose index 90C integer ia : Looping variable for number of atoms 91C integer i1,i2,i3 : Mesh indexes in each mesh direction 92C integer is : Species index 93C integer isp : Sub-Point index 94C integer j : General-porpose index 95C integer j1,j2,j3 : Mesh indexes in each mesh direction 96C integer nep : Number of extended-mesh points 97C integer nmp : Number of mesh points in unit cell 98C real*8 pldist : Distance between mesh planes 99C real*8 r : Distance between atom and mesh point 100C real*8 vecmod : Vector modulus 101C real*8 volume : Unit cell volume 102C logical within : Is a mesh point within orbital range? 103C ---------------------------------------------------------------------- 104C Units : 105C ---------------------------------------------------------------------- 106C 107C Energies in Rydbergs 108C Distances in Bohr 109C 110 111C Local variables 112 integer :: i, ia, i1, i2, i3, io, iphi, is, isp, 113 & ity, j, nep, j1, j2, j3, indi 114 integer(i8b) :: ntp 115 real(dp) :: dx(3), dxp(3), pldist, r, 116 & volume 117 logical :: within 118C Functions 119 real(dp), external :: dismin 120!--------------------------------------------------------------------- BEGIN 121#ifdef DEBUG 122 call write_debug( ' PRE InitMesh' ) 123#endif 124 125! ---------------------------------------------------------------------- 126! Mesh initialization 127! ---------------------------------------------------------------------- 128 ! determines ntm 129 call get_mesh_sizes(cell,G2max,nsm,ntm,G2mesh) 130 131! Save some numbers 132! nmeshg in module array 133 nm(:) = ntm(:) / nsm 134 nmsc(:) = nm(:) * nsc(:) 135 nmeshg(:) = ntm(:) ! Total number of mesh points, incl. subpoints 136 nmuc(:) = nm(:) ! Mesh points in each unit cell vector 137 nusc(:) = nsc(:) ! Unit cells in each supercell direction 138 139! Find number of mesh points in unit cell. 140! Note needed integer*8 ! 141!!AG** Check proper kind name 142 143 ntp = product(INT(ntm,kind=i8b)) 144 145C Create the first mesh distribution 146 call initMeshDistr( oDistr=UNIFORM, nm=nm ) 147 148C Find and sets local number of Mesh points of each kind 149! Sets nml, nmpl, ntml and ntpl 150 call setMeshDistr( UNIFORM, nsm, nsp, nml, nmpl, ntml, ntpl ) 151 152C Find volume of unit cell and of mesh cell 153 volume = volcel( cell ) 154 dvol = volume / ntp 155 156C Output current mesh dimensions and cut-off 157 if (Node.eq.0) then 158 write(6,'(/,a,3(i6,a),i12)') 'InitMesh: MESH =', 159 & ntm(1),' x',ntm(2),' x',ntm(3),' =', ntp 160 write(6,'(a,3(i6,a),i12)') 'InitMesh: (bp) =', 161 & nm(1),' x',nm(2),' x',nm(3),' =', product(INT(nm,8)) 162 write(6,'(a,2f10.3,a)') 163 & 'InitMesh: Mesh cutoff (required, used) =', 164 & G2max, G2mesh, ' Ry' 165 endif 166 if (cml_p) then 167 call cmlStartPropertyList(mainXML) 168 call cmlAddProperty(xf=mainXML, value=ntm, 169 . dictref='siesta:ntm', title='Mesh', 170 . units='cmlUnits:countable') 171 call cmlAddProperty(xf=mainXML, value=G2max, 172 . units='siestaUnits:Ry', 173 . dictref='siesta:g2max', title='Requested Cut-Off') 174 call cmlAddProperty(xf=mainXML, value=G2mesh, 175 . units='siestaUnits:Ry', 176 . dictref='siesta:g2mesh', title='Actual Cut-Off') 177 call cmlEndPropertyList(mainXML) 178 endif 179 G2max = G2mesh 180 181C Find mesh-cell vectors 182 do i = 1,3 183 do j = 1,3 184 cmesh(j,i) = cell(j,i) / nm(i) 185 enddo 186 enddo 187 188C Find reciprocal mesh-cell vectors (not multiplied by 2*pi) 189 call reclat( cmesh, rcmesh, 0 ) 190 191! Find number ne(:) of extended-mesh intervals for each cell vector. 192 do i = 1,3 193 ! pldist is the distance between mesh planes 194 pldist = 1.0_dp / sqrt(dot_product(rcmesh(:,i),rcmesh(:,i))) 195 ! Find number of planes spanned by rmax 196 ne(i) = rmax / pldist 197 enddo ! i 198 199! For an atom at x=0, ne=rmax/pldist is the last mesh point within rmax. 200! But an atom at mesh cell ix=0 can be almost at x=dx. Therefore we have 201! to go up to ne+1. And subpoints in mesh cell ix=-ne-1 can be almost at 202! x=-ne*dx. Therefore we have to go up to -ne-1 to the left. Thus, we 203! just increase ne and forget about these two effects from now on. 204 ne(:) = ne(:) + 1 205 206C Find sub-points 207!$OMP parallel default(shared), private(i1,i2,i3,i,isp) 208 209!$OMP do 210 do i3 = 0, nsm-1 211 isp = i3 * nsm ** 2 212 do i2 = 0, nsm-1 213 do i1 = 0, nsm-1 214 isp = isp + 1 215 do i = 1,3 216 xdsp(i,isp) = ( cmesh(i,1) * i1 + 217 & cmesh(i,2) * i2 + 218 & cmesh(i,3) * i3 ) / nsm 219 enddo 220 enddo 221 enddo 222 enddo 223!$OMP end do 224 225! Find number of points within rmax (orbital points) 226!$OMP single ! Only used for OpenMP 227 mop = 0 228!$OMP end single 229!$OMP do private(dxp,within,dx,r), reduction(+:mop) 230 do i3 = -ne(3), ne(3) ! Loop over neighbor (super) points 231 do i2 = -ne(2), ne(2) 232 do i1 = -ne(1), ne(1) 233 dxp(:) = cmesh(:,1) * i1 + ! (Super) point coordinates 234 . cmesh(:,2) * i2 + 235 . cmesh(:,3) * i3 236 ! Find if any subpoint can be within rmax of an atom that 237 ! is in the origin's mesh cell 238 within = .false. 239 do isp = 1,nsp ! Loop over sub-points 240 dx(1:3) = dxp(1:3) + xdsp(1:3,isp) ! Subpoint coords. 241 ! Find distance from point to the mesh cell at the origin 242 ! (since the atom might be at any point of the mesh cell) 243 r = dismin( cmesh, dx ) 244 if ( r .lt. rmax ) within = .true. 245 enddo ! isp 246 if ( within ) then 247 mop = mop + 1 ! Total number of points within rmax 248 endif ! (within) 249 enddo ! i1 250 enddo ! i2 251 enddo ! i3 252!$OMP end do nowait 253 254!$OMP end parallel 255 256C Create extended mesh arrays for the first data distribution 257 call setupExtMesh( UNIFORM, rmax ) 258 259#ifdef DEBUG 260 call write_debug( ' POS InitMesh' ) 261#endif 262!----------------------------------------------------------------------- END 263 end subroutine InitMesh 264!------------------------------------ 265 266 subroutine get_mesh_sizes(cell,G2max,nsm,ntm,RealCutoff) 267 ! Computes the mesh dimensions ntm for a given cell 268 ! and mesh cutoff, making them multiples of the big-point 269 ! sampling size nsm. 270 271 use precision, only : dp 272 use cellsubs, only : reclat ! Finds reciprocal unit cell 273 use m_chkgmx, only : chkgmx ! Checks planewave cutoff of a mesh 274 use fft1d, only : nfft ! Finds allowed value for 1-D FFT 275 use fdf 276 use sys, only : die, message 277 278 real(dp), intent(in) :: cell(3,3) 279 real(dp), intent(in) :: G2max 280 real(dp), intent(out) :: RealCutoff 281 integer, intent(in) :: nsm 282 integer, intent(out) :: ntm(3) 283 284 real(dp) :: rcell(3,3), vecmod 285 integer :: i, itmp 286 real(dp), parameter :: k0(3) = (/ 0.0, 0.0, 0.0 /) 287 288 type(block_fdf) :: bfdf 289 type(parsed_line), pointer :: pline 290 291 ! Find reciprocal cell vectors (multiplied by 2*pi) 292 call reclat( cell, rcell, 1 ) 293 294! Find number of mesh intervals for each cell vector. 295! The reciprocal vectors of the mesh unit cell (cell/ntm) 296! are rcell*ntm, and must be larger than 2*G2max 297 do i = 1,3 298 vecmod = sqrt(dot_product(rcell(:,i),rcell(:,i))) 299 ntm(i) = 2 * sqrt(G2max) / vecmod + 1 300 enddo 301 302! Impose that mesh cut-off is large enough 303 impose_cutoff: do 304 305 ! Require that ntm is suitable for FFTs and a multiple of nsm 306 do i = 1,3 307 impose_fft: do 308 ! nfft increases ntm to the next integer suitable for FFTs 309 call nfft( ntm(i) ) 310 if ( mod( ntm(i), nsm )==0 ) exit impose_fft 311 ntm(i) = ntm(i) + 1 312 end do impose_fft 313 enddo ! i 314 315 ! Check that effective cut-off is large enough as for non-right 316 ! angled unit cells this is not guaranteed to be the case. 317 ! If cut-off needs to be larger, increase ntm and try again. 318 ! chkgmx decreases G2mesh to the right cutoff value of the mesh 319 RealCutoff = huge(1.0_dp) 320 call chkgmx( k0, rcell, ntm, RealCutoff ) 321 if (RealCutoff >= G2max) then 322 exit impose_cutoff 323 else 324 ntm(1:3) = ntm(1:3) + 1 325 end if 326 327 end do impose_cutoff 328 329 330! Let the user manually decide the mesh-box size 331! This option may be useful when testing Mesh.SubDivisions 332! for performance. 333! NOTE/TODO: manual specification is not checked with respect 334! to nfft (allowed integers). I don't know if this should 335! forced or not? 336 if ( fdf_block('Mesh.Sizes',bfdf) ) then 337 338 if ( fdf_bline(bfdf, pline) ) then 339 ntm(1) = fdf_bintegers(pline,1) 340 ntm(2) = fdf_bintegers(pline,2) 341 ntm(3) = fdf_bintegers(pline,3) 342 else 343 call die('initmesh: ERROR in Mesh.Sizes block') 344 endif 345 call fdf_bclose(bfdf) 346 347 ! Calculate the actual meshcutoff from the grid 348 RealCutoff = huge(1.0_dp) 349 call chkgmx( k0, rcell, ntm, RealCutoff ) 350 351 else if ( fdf_islist('Mesh.Sizes') ) then 352 353 i = -1 354 call fdf_list('Mesh.Sizes', i, ntm) 355 if ( i /= 3 ) 356 & call die('initmesh: ERROR in Mesh.Sizes list') 357 ! Do the actual read of the mesh sizes 358 call fdf_list('Mesh.Sizes', i, ntm) 359 360 ! Calculate the actual meshcutoff from the grid 361 RealCutoff = huge(1.0_dp) 362 call chkgmx( k0, rcell, ntm, RealCutoff ) 363 364 end if 365 366! Issue a warning if the mesh sizes are smaller than 367! the requested mesh cutoff 368 if ( RealCutoff < G2max ) then 369 call message("WARNING","Read mesh sizes give a smaller" 370 $ // " cutoff than requested") 371 end if 372 373! Check whether the resulting mesh sizes are multiples 374! of nsm. 375 if ( any(mod(ntm,nsm) /= 0) ) then 376 call die("initmesh: Read mesh sizes are not multiples of nsm") 377 end if 378 379! Check that the numbers are divisible with the FFT algorithms 380 do i = 1, 3 381 itmp = ntm(i) 382 call nfft( itmp ) 383 if ( itmp /= ntm(i) ) then 384 call die("initmesh: Read mesh sizes are not compatible " 385 & // "with the FFT algorithm: % [2, 3, 5] == 0") 386 end if 387 end do 388 389 end subroutine get_mesh_sizes 390 391!-------------------------- 392 393 subroutine SameMeshAndAtoms(na, xa, ucell, rmax, G2max, G2mesh, 394 & samesh, samexa) 395C 396C Checks whether anything has changed that requires the mesh to be 397C reinitialised or quantities relating to the atoms relative to 398C the mesh. 399C 400C ---------------------------------------------------------------------- 401C Input : 402C ---------------------------------------------------------------------- 403C integer na : Number of atoms in the supercell 404C real*8 xa(3,na) : Coordinates of the atoms in the supercell 405C real*8 ucell(3,3) : Current unit cell vectors 406C real*8 rmax : Maximum orbital radius 407C real*8 G2max : Requested mesh cut-off 408C real*8 G2mesh : Mesh cut-off from last call 409C ---------------------------------------------------------------------- 410C Output : 411C ---------------------------------------------------------------------- 412C logical samesh : If .true. then the mesh must be reinitialised 413C logical samexa : If .true. then atom related quantities need 414C : to be recalculated 415C ---------------------------------------------------------------------- 416C Internal : 417C ---------------------------------------------------------------------- 418C integer lastna : Number of atoms from last call 419C real*8 lastxa(3,na) : Position of atoms from last call 420C real*8 lastra : Rmax from last call 421C real*8 lastc(3,3) : Unit cell from last call 422C ---------------------------------------------------------------------- 423 use precision, only : dp 424 use alloc, only : re_alloc, de_alloc 425 implicit none 426C Passed arguments 427 integer, intent(in) :: na 428 real(dp), intent(in) :: G2max, G2mesh 429 real(dp), intent(in) :: rmax 430 real(dp), intent(in) :: ucell(3,3) 431 real(dp), intent(in) :: xa(3,na) 432 logical, intent(out) :: samesh, samexa 433C Local variables 434C Saved 435 integer, save :: lastna = 0 436 real(dp), save :: tiny = 1.0e-12_dp, 437 & lastc(3,3) = 0.743978657912656e50_dp, 438 & lastra = 0.0_dp 439C non-saved 440 integer :: i, ia, j 441 real(dp), pointer, save :: lastxa(:,:) 442 443!------------------------------------------------------------------------- BEGIN 444C If number of atoms has changed then deallocate lastxa 445 if (na.ne.lastna) then 446 call re_alloc(lastxa, 1, 3, 1, na, 'lastxa', 'SameMeshAndAtoms') 447 lastxa(1:3,1:na) = 0.0_dp 448 lastxa(1,1) = 0.987654321273567e50_dp 449 endif 450 451C Find if mesh has to be changed due to unit cell 452 samesh = .true. 453 do i = 1,3 454 do j = 1,3 455 if ( ucell(j,i) .ne. lastc(j,i) ) samesh = .false. 456 lastc(j,i) = ucell(j,i) 457 enddo 458 enddo 459 460C Find if mesh has to be changed due to unit cell 461 if ( G2max .gt. G2mesh * (1.0_dp + tiny) ) samesh = .false. 462 463C Find if mesh has to be changed due to rmax 464 if (rmax .ne. lastra) samesh = .false. 465 lastra = rmax 466 467C Find if atoms have moved having checked the number of atoms first 468 samexa = (na.eq.lastna) 469 if (samexa) then 470 do ia = 1,na 471 do i = 1,3 472 if ( xa(i,ia) .ne. lastxa(i,ia) ) samexa = .false. 473 enddo 474 enddo 475 endif 476 477C Copy the number of atoms and coordinates to save for next call 478 lastna = na 479 lastxa(1:3,1:na) = xa(1:3,1:na) 480 481C If cell has changed then it is necessary to reinitialise coordinates 482 if (.not.samesh) samexa = .false. 483 484!--------------------------------------------------------------------------- END 485 end subroutine SameMeshAndAtoms 486 487! subroutine InitAtomMesh 488! ---------------------------------------------------------------------- 489! Initialises the atom information relating to the mesh for a given data 490! distribution. It creates and initialise ipa. This arrays is allocated 491! inside a meshDisType structure and can be accesed by pointers of the 492! mesh module. 493! 494! Only necesary for those data distributions who needs to use the 495! extended mesh. 496! 497! ---------------------------------------------------------------------- 498! Input : 499! ---------------------------------------------------------------------- 500! integer distr : Used data distribution 501! integer na : number of atoms 502! integer xa(3,na) : Atomic coordinates 503! ---------------------------------------------------------------------- 504! Output : 505! ---------------------------------------------------------------------- 506! All output quantities are in the mesh and moreMeshSubs modules. 507! The ipa array is stored in moreMeshSubs but accessed using the mesh 508! module pointer. 509! 510!*********************************************************************** 511 subroutine InitAtomMesh( distr, na, xa ) 512 use precision, only: dp 513 use mesh, only: ipa ! atoms' mesh points within my Extended mesh 514 use mesh, only: dxa ! atom's position relative to mesh point 515 use mesh, only: cmesh ! Mesh cell vectors 516 use mesh, only: meshLim ! My processor's box of mesh points 517 use mesh, only: ne ! Points in rmax along each lat. vector 518 use mesh, only: nmsc ! Mesh points in each supercell vector 519 use mesh, only: rcmesh ! Reciprocal mesh cell vectors 520 use mesh, only: iatfold ! Supercell vector that keeps track of the 521 ! of the folding of the atomic coordinates 522 ! in the mesh 523 use alloc,only: re_alloc ! (Re)allocation routine 524 use moreMeshSubs, only : allocIpaDistr 525 implicit none 526! Input parameters 527 integer, intent(in) :: na ! Number of atoms in supercell 528 real(dp), intent(in) :: xa(3,na) ! Atomic coordinates 529 integer, intent(in) :: distr ! Used data distribution 530! Local variables 531 character(len=*),parameter:: myName = 'InitAtomMesh ' 532 integer, save :: lastna = 0 ! Number of atoms on last call 533 integer :: ia, ixa(3), jsc(3), jxa(3), nem(3), 534 & j1, j2, j3, myBox(2,3), myExtBox(2,3) 535 integer :: ixabeffold(3) 536 real(dp) :: cxa(3) 537 538#ifdef DEBUG 539 call write_debug( ' PRE InitAtomMesh' ) 540#endif 541! (Re)allocate atom-mesh arrays 542 call allocIpaDistr( distr, na ) 543 if (na.ne.lastna) then 544 call re_alloc( dxa, 1,3, 1,na, myName//'dxa' ) 545 call re_alloc( iatfold, 1,3, 1,na, myName//'iatfold' ) 546 endif 547 lastna = na 548 549 myBox(:,:) = meshLim(:,:) - 1 550 551! Add 'wings extensions' to myBox of mesh points, containing all points 552! within the orbital spheres that may intersect myBox. Thus, wings must 553! be one diameter wide. 554 myExtBox(1,:) = myBox(1,:) - 2*ne(:) 555 myExtBox(2,:) = myBox(2,:) + 2*ne(:) 556 557! Find number of extended-box points. 558 nem(:) = myExtBox(2,:) - myExtBox(1,:) + 1 559 560! Find atomic positions relative to mesh 561 do ia = 1,na 562 563! Find atomic coordinates in mesh basis vectors 564 cxa(:) = matmul( xa(:,ia), rcmesh(:,:) ) 565 ixabeffold(:) = floor(cxa(:)) 566 ! fix: Make sure we consider equivalent positions 567 ! inside the cell, even if some atoms are outside. 568 cxa(:) = modulo( cxa(:), real(nmsc(:),dp) ) 569 570! Find indexes of mesh cell in which atom is 571 ixa(:) = floor(cxa(:)) 572 iatfold(:,ia) = ( ixa(:) - ixabeffold(:) ) / nmsc(:) 573 574! Find atom position within mesh cell 575 cxa(:) = cxa(:) - ixa(:) ! Mesh coordinates 576 dxa(:,ia) = matmul( cmesh(:,:), cxa(:) ) ! Cartesian coords 577 578! Find atom's mesh index within myExtBox, but only if the 579! (periodic) atomic sphere (of width ne) can intersect myBox. 580! Otherwise, make ipa=0. nem is the width of myExtBox. 581! nmsc is the width of the supercell. 582! Notice that, even if several periodic images of the same 583! atom may fit within myExtBox, only one must be considered, 584! since the points within its sphere, once folded into myBox, 585! will be equivalent to those of any other periodic image. 586 ipa(ia) = 0 587 sc: do j3 = -1,1 ! Loop on neighbor supercells 588 do j2 = -1,1 589 do j1 = -1,1 590 jsc(:) = (/j1,j2,j3/) 591 jxa(:) = ixa(:) + jsc(:)*nmsc(:) 592 if ( all(jxa(:)>=myBox(1,:)-ne(:)) .and. 593 & all(jxa(:)<=myBox(2,:)+ne(:)) ) then 594 jxa(:) = jxa(:) - myExtBox(1,:) 595 ipa(ia) = 1 + jxa(1) + nem(1)*jxa(2) 596 & + nem(1)*nem(2)*jxa(3) 597 598 iatfold(1,ia) = iatfold(1,ia) + j1 599 iatfold(2,ia) = iatfold(2,ia) + j2 600 iatfold(3,ia) = iatfold(3,ia) + j3 601 602 exit sc 603 end if 604 end do ! j1 605 end do ! j2 606 end do sc ! j3 607 enddo ! ia 608 609#ifdef DEBUG 610 call write_debug( ' POS InitAtomMesh' ) 611#endif 612 end subroutine InitAtomMesh 613 614 subroutine PartialCoreOnMesh( na, isa, ntpl, rhopcc, indxua, 615 & nsd, dvol, volume, Vscf, Vaux, 616 & fal, stressl, Forces, Stress ) 617C 618C Finds the partial-core-correction energy density on the mesh 619C 620C ---------------------------------------------------------------------- 621C Input : 622C ---------------------------------------------------------------------- 623C integer na : Number of atoms in supercell 624C integer isa(na) : Species index of all atoms in supercell 625C integer ntpl : Number of mesh Total Points in unit cell 626C : (including subpoints) locally 627C integer indxua : Index of equivalent atom in unit cell 628C integer nsd : Number of diagonal spin values (1 or 2) 629C real*8 dvol : Mesh-cell volume 630C real*8 volume : Unit cell volume 631C real*4 Vscf(ntpl) : Hartree potential of SCF density 632C real*4 Vaux(ntpl) : Auxiliary potential array 633C logical Forces : Are the forces to be calculated? 634C logical Stress : Are the stresses to be calculated? 635C ---------------------------------------------------------------------- 636C Output : (non-gradient case) 637C ---------------------------------------------------------------------- 638C real*4 rhopcc(ntpl) : Partial-core-correction density for xc 639C ---------------------------------------------------------------------- 640C Output : (gradient case) 641C ---------------------------------------------------------------------- 642C real*8 fal(3,:) : Local copy of atomic forces 643C real*8 stressl(3,3) : Local copy of stress tensor 644C ---------------------------------------------------------------------- 645 use precision, only: dp, grid_p 646 use atmfuncs, only: rcore, chcore_sub 647 use mesh, only: dxa, idop, indexp, ipa, mop, nsp, xdop, xdsp 648 implicit none 649C Passed arguments 650 integer, intent(in) :: na, isa(na), ntpl 651 integer, intent(in) :: indxua(*) 652 integer, intent(in) :: nsd 653 real(dp), intent(in) :: dvol, volume 654 real(grid_p), intent(in) :: Vscf(ntpl,*) 655 real(grid_p), intent(in) :: Vaux(ntpl) 656 logical, intent(in) :: Forces 657 logical, intent(in) :: Stress 658 real(grid_p), intent(out) :: rhopcc(ntpl) 659 real(dp), intent(inout) :: fal(3,*) 660 real(dp), intent(inout) :: stressl(3,3) 661C Internal variables 662 integer :: i, ia, iop, ip, ip0, is, isp, 663 & ispin, iua 664 real(dp) :: dfa(3), dx(3), grrho(3), r, ra, 665 & rhop, vxc 666 real(dp), save :: tiny 667 logical :: Gradients 668 data tiny / 1.0e-12_dp / 669 670!------------------------------------------------------------------------- BEGIN 671#ifdef DEBUG 672 call write_debug( ' PRE PartialCoreOnMesh' ) 673#endif 674C Find out whether this is a gradient run based on the arguments passed 675 Gradients = ( Forces .or. Stress ) 676 677C Initialise array to zero 678 if (.not.Gradients) rhopcc(1:ntpl) = 0.0_dp 679 680 do ia = 1,na 681 if (ipa(ia)==0) cycle ! Atomic sphere does not intersec my Box 682 iua = indxua(ia) 683 is = isa(ia) 684 ra = rcore( is ) 685 if (ra .gt. tiny) then 686C Loop over mesh points inside rmax 687 do iop = 1,mop 688 ip0 = indexp( ipa(ia) + idop(iop) ) 689 if (ip0 .gt. 0) then 690C Loop over sub-points 691 do isp = 1,nsp 692 dx(1:3) = xdop(1:3,iop) + xdsp(1:3,isp) - dxa(1:3,ia) 693 r = sqrt(dot_product(dx,dx)) 694C if (r .lt. ra .and. r .gt. tiny) then 695 if (r .lt. ra) then 696 ip = isp + nsp * (ip0 - 1) 697 call chcore_sub( is, dx, rhop, grrho ) 698 if (Gradients) then 699C Calculate gradients of PCC 700 do ispin = 1,nsd 701 vxc = Vscf(ip,ispin) - Vaux(ip) 702 do i = 1,3 703 dfa(i) = dvol * vxc * grrho(i) / nsd 704 if (Forces) fal(i,iua) = fal(i,iua) + dfa(i) 705 if (Stress) stressl(1:3,i) = stressl(1:3,i) + 706 & dx(1:3)*dfa(i)/volume 707 enddo 708 enddo 709 else 710C Calculate density due to PCC 711 rhopcc(ip) = rhopcc(ip) + rhop 712 endif 713 endif 714 enddo 715 endif 716 enddo 717 endif 718 enddo 719#ifdef DEBUG 720 call write_debug( ' POS PartialCoreOnMesh' ) 721#endif 722!--------------------------------------------------------------------------- END 723 end subroutine PartialCoreOnMesh 724 725 subroutine NeutralAtomOnMesh( na, isa, ntpl, vna, 726 & indxua, dvol, volume, 727 & drho, fal, stressl, 728 & Forces, Stress ) 729C 730C Finds the potential due to the neutral atoms on the mesh 731C 732C ---------------------------------------------------------------------- 733C Input : 734C ---------------------------------------------------------------------- 735C integer na : Number of atoms in supercell 736C integer isa(na) : Species index of all atoms in supercell 737C integer ntpl : Number of mesh Total Points in unit cell 738C : (including subpoints) locally 739C integer indxua : Index of equivalent atom in unit cell 740C real*8 dvol : Mesh-cell volume 741C real*8 volume : Unit cell volume 742C real drho(ntpl) : SCF density difference 743C logical Forces : Should the forces be calculated? 744C logical Stress : Should the stress be calculated? 745C ---------------------------------------------------------------------- 746C Input/Output : (Output for non-gradient call, must be kept otherwise) 747C ---------------------------------------------------------------------- 748C real vna(ntpl) : Sum of neutral-atom potentials 749C ---------------------------------------------------------------------- 750C Output : (gradient call) 751C ---------------------------------------------------------------------- 752C real*8 fal(3,:) : Local copy of atomic forces 753C real*8 stressl(3,3) : Local copy of stress tensor 754C ---------------------------------------------------------------------- 755 use precision, only: dp, grid_p 756 use atmfuncs, only: rcut, phiatm 757 use mesh, only: dxa, idop, indexp, ipa, mop, nsp, xdop, xdsp 758 implicit none 759C Passed arguments 760 integer, intent(in) :: na, isa(na), ntpl 761 integer, intent(in) :: indxua(*) 762 real(dp), intent(in) :: dvol, volume 763 real(grid_p), intent(in) :: drho(*) 764 logical, intent(in) :: Forces, Stress 765 real(grid_p), intent(inout) :: vna(ntpl) 766 real(dp), intent(inout) :: fal(3,*) 767 real(dp), intent(inout) :: stressl(3,3) 768C Internal variables 769 integer :: i, ia, iop, ip, ip0, is, isp, iua 770 real(dp) :: dx(3), grva(3), r, ra, va 771 logical :: Gradients 772 773#ifdef DEBUG 774 call write_debug( ' PRE NeutralAtomOnMesh' ) 775#endif 776C Check whether forces and/or stress has been 777C requested based on arguments 778 Gradients = (Forces.or.Stress) 779 780C Initialise array to zero if not computing gradients 781C If we are computing gradients, the current value of vna 782C must be kept. This is particularly important for grid-cell sampling... 783 784 if (.not.Gradients) vna(1:ntpl) = 0.0_dp 785 786 do ia = 1,na 787 if (ipa(ia)==0) cycle ! Atomic sphere does not intersec my Box 788 iua = indxua(ia) 789 is = isa(ia) 790 ra = rcut( is, 0 ) 791 792C Loop over mesh points inside rmax 793 do iop = 1,mop 794 ip0 = indexp( ipa(ia) + idop(iop) ) 795 if (ip0 .gt. 0) then 796C Loop over sub-points 797 do isp = 1,nsp 798 dx(1:3) = xdop(1:3,iop) + xdsp(1:3,isp) - dxa(1:3,ia) 799 r = sqrt(dot_product(dx,dx)) 800 if (r .lt. ra) then 801 ip = isp + nsp * (ip0 - 1) 802 call phiatm( is, 0, dx, va, grva ) 803 if (Gradients) then 804 do i = 1,3 805 grva(i) = dvol * grva(i) * drho(ip) 806 if (Forces) fal(i,iua) = fal(i,iua) + grva(i) 807 if (Stress) stressl(1:3,i) = stressl(1:3,i) + 808 & dx(1:3) * grva(i) / volume 809 enddo 810 else 811 vna(ip) = vna(ip) + va 812 endif 813 endif 814 enddo 815 endif 816 enddo 817 enddo 818#ifdef DEBUG 819 call write_debug( ' POS NeutralAtomOnMesh' ) 820#endif 821!--------------------------------------------------------------------------- END 822 end subroutine NeutralAtomOnMesh 823 824 subroutine LocalChargeOnMesh( na, isa, ntpl, Chlocal, indxua ) 825C 826C Finds the diffuse ionic charge, whose electrostatic potential is equal 827C to Vlocal on the mesh 828C 829C ---------------------------------------------------------------------- 830C Input : 831C ---------------------------------------------------------------------- 832C integer na : Number of atoms in supercell 833C integer isa(na) : Species index of all atoms in supercell 834C integer ntpl : Number of mesh Total Points in unit cell 835C : (including subpoints) locally 836C integer indxua : Index of equivalent atom in unit cell 837C ---------------------------------------------------------------------- 838C Output : 839C ---------------------------------------------------------------------- 840C real Chlocal(ntpl) : Sum of diffuse ionic charge densities 841C ---------------------------------------------------------------------- 842 use precision, only: dp, grid_p 843 use atmfuncs, only: rcut, psch 844 use atm_types, only: species 845 use radial, only: rad_func 846 use mesh, only: dxa, idop, indexp, ipa, mop, nsp, xdop, xdsp 847 implicit none 848C Passed arguments 849 integer, intent(in) :: na, isa(na), ntpl 850 integer, intent(in) :: indxua(*) 851 real(grid_p), intent(out) :: Chlocal(ntpl) 852C Internal variables 853 integer :: ia, iop, ip, ip0, is, isp, iua 854 real(dp) :: dx(3), grpsch(3), r, ra, pschloc 855 type(rad_func), pointer :: func 856 857!------------------------------------------------------------------------- BEGIN 858#ifdef DEBUG 859 call write_debug( ' PRE LocalChargeOnMesh' ) 860#endif 861C Initialise array to zero 862 Chlocal(1:ntpl) = 0.0_grid_p 863 864 do ia = 1,na 865 if (ipa(ia)==0) cycle ! Atomic sphere does not intersec my Box 866 iua = indxua(ia) 867 is = isa(ia) 868 func => species(is)%chlocal 869 ra = func%cutoff 870C Loop over mesh points inside rmax 871 do iop = 1,mop 872 ip0 = indexp( ipa(ia) + idop(iop) ) 873 if (ip0 .gt. 0) then 874C Loop over sub-points 875 do isp = 1,nsp 876 dx(1:3) = xdop(1:3,iop) + xdsp(1:3,isp) - dxa(1:3,ia) 877 r = sqrt(dot_product(dx,dx)) 878 if (r .lt. ra) then 879 ip = isp + nsp * (ip0 - 1) 880 call psch( is, dx, pschloc, grpsch ) 881 Chlocal(ip) = Chlocal(ip) + pschloc 882 endif 883 enddo 884 endif 885 enddo 886 enddo 887#ifdef DEBUG 888 call write_debug( ' POS LocalChargeOnMesh' ) 889#endif 890!--------------------------------------------------------------------------- END 891 end subroutine LocalChargeOnMesh 892 893 subroutine ModelCoreChargeOnMesh( na, isa, ntpl, 894 $ ModelCore, indxua ) 895C 896C Constructs a model total core charge, useful for Bader analysis 897C 898C ---------------------------------------------------------------------- 899C Input : 900C ---------------------------------------------------------------------- 901C integer na : Number of atoms in supercell 902C integer isa(na) : Species index of all atoms in supercell 903C integer ntpl : Number of mesh Total Points in unit cell 904C : (including subpoints) locally 905C integer indxua : Index of equivalent atom in unit cell 906C ---------------------------------------------------------------------- 907C Output : 908C ---------------------------------------------------------------------- 909C real ModelCore(ntpl) : Sum of model core charge densities 910C ---------------------------------------------------------------------- 911 912 use precision, only : dp, grid_p 913 use atmfuncs, only: izofis, zvalfis 914 use atm_types, only : species 915 use mesh, only: dxa, idop, indexp, ipa, mop, nsp, xdop, xdsp 916 use fdf 917 use parallel, only: ionode 918 919 implicit none 920 921C 922C Passed arguments 923C 924 integer, intent(in) :: na, isa(na), ntpl 925 integer, intent(in) :: indxua(*) 926 real(grid_p), intent(out) :: ModelCore(ntpl) 927 928C 929C Internal variables 930C 931 integer 932 . ia, iop, ip, ip0, is, isp, iua 933 934 real(dp) dx(3), ra, charge_ratio, r 935 real(dp) core_charge 936 real(dp), parameter :: pi = 3.141592653589793_dp 937 938 ! We want a total core charge of (Z-Zval), localized in a 939 ! region of around 1.0 bohr. 940 ! Note that we put 1 electron (extra!) for Hydrogen atoms, 941 ! to provide the necessary "cuspy" behavior for the reference 942 ! file for Bader analysis. For H, the radius is 0.6 bohr 943 real(dp) :: radius_standard, radius_h 944 945 ! The model function is N*f(r/ra), where 946 ! f(x) = exp( - (sinh(2.0_dp*x)/sinh(1.0_dp))**2 ) 947 ! f(1) is zero to within 10^-5 948 ! For normalization, we need \int_{0,1} { x*x*f(x) }, which is: 949 real(dp), parameter :: Integral = 0.0411548_dp 950 951 ! so that N above is Q/(4*pi*Integral*ra**3) 952 953 ! Note that this is quite a small confinement region, so 954 ! you might need a high cutoff (of the order of 300-500 Ry) 955 ! to represent it on the grid. It might be better to use 956 ! Denchar for this, rather than performing a Siesta calculation 957 ! (even if it reads the converged DM) 958 959 960 ! Initialise array to zero 961 ModelCore(1:ntpl) = 0.0_grid_p 962 963 radius_standard = fdf_get("bader-core-radius-standard", 964 $ 1.0_dp,"Bohr") 965 radius_h = fdf_get("bader-core-radius-hydrogen", 966 $ 0.6_dp,"Bohr") 967 968 if (ionode) then 969 write(6,"(a,2f6.3)") "Bader Analysis core-charge setup. " // 970 $ "Radii (standard, H): ", 971 $ radius_standard, radius_h 972 endif 973 974 do ia = 1,na 975 if (ipa(ia)==0) cycle ! Atomic sphere does not intersect my Box 976 is = isa(ia) 977 if (izofis(is) == 1) then 978 ra = radius_h 979 charge_ratio = 1.0_dp / 980 $ (4.0_dp*pi*Integral*ra**3.0_dp) 981 else 982 ra = radius_standard 983 charge_ratio = (izofis(is)-zvalfis(is)) / 984 $ (4.0_dp*pi*Integral*ra**3.0_dp) 985 endif 986 987 988C Loop over mesh points inside rmax 989 do iop = 1,mop 990 ip0 = indexp( ipa(ia) + idop(iop) ) 991 if (ip0 .gt. 0) then 992 993C Loop over sub-points 994 do isp = 1,nsp 995 dx(1:3) = xdop(1:3,iop) + xdsp(1:3,isp) - dxa(1:3,ia) 996 r = sqrt(dot_product(dx,dx)) 997 if (r .lt. ra) then 998 ip = isp + nsp * (ip0 - 1) 999 core_charge = charge_ratio*model(r/ra) 1000 ModelCore(ip) = ModelCore(ip) + core_charge 1001 endif 1002 enddo 1003 1004 endif 1005 enddo 1006 1007 enddo 1008 1009 CONTAINS 1010 1011 function model(x) 1012 real(dp), intent(in) :: x 1013 real(dp) :: model 1014 1015 model = exp( - (sinh(2.0_dp*x)/sinh(1.0_dp))**2 ) 1016 end function model 1017 1018 end subroutine ModelCoreChargeOnMesh 1019! 1020 subroutine PhiOnMesh( nmpl, norb, iaorb, iphorb, isa, numphi ) 1021C 1022C Calculates the values of the orbitals at the mesh points 1023C 1024C ---------------------------------------------------------------------- 1025C Input : 1026C ---------------------------------------------------------------------- 1027C integer nmpl : Number of mesh points in unit cell locally 1028C integer norb : Total number of basis orbitals in supercell 1029C integer iaorb(norb) : Atom to which each orbital belongs 1030C integer iphorb(norb) : Orbital index (within atom) 1031C integer isa(na) : Species index of all atoms in supercell 1032C ---------------------------------------------------------------------- 1033C Output : 1034C ---------------------------------------------------------------------- 1035C All output quantities are in the module meshphi 1036C ---------------------------------------------------------------------- 1037 1038C Modules 1039 use precision, only : dp, grid_p 1040 use atmfuncs, only : rcut, phiatm 1041 use mesh, only : dxa, idop, indexp, ipa, mop, nsp, xdop, 1042 & xdsp, meshLim 1043 use meshphi, only : directphi, nphi, phi, lstpht, listp2, endpht 1044 use meshphi, only : resetMeshPhi 1045 use fdf, only : fdf_get 1046 use parallel, only : IONode 1047 use alloc, only : re_alloc 1048 1049 implicit none 1050C Passed arguments 1051 integer, intent(in) :: nmpl, norb, iaorb(norb), 1052 & iphorb(norb), isa(*) 1053 1054 ! This array was set in distriPhiOnMesh, and holds 1055 ! the pre-computed count of orbitals touching a given 1056 ! point. It is re-used and overwritten here. 1057 integer, intent(inout) :: numphi(nmpl) 1058 1059C Local variables 1060 integer :: ia, io, iop, ip, ip0, iphi, is, isp, n, 1061 & nlist, nliste, nmax 1062 1063 logical :: within 1064 logical, save :: firsttime = .true. 1065 real(dp) :: dxsp(3,nsp), grphi(3), phip, r2o, 1066 & r2sp(nsp) 1067 1068!--------------------------------------------------------------------- BEGIN 1069#ifdef DEBUG 1070 call write_debug( ' PRE PhiOnMesh' ) 1071#endif 1072 1073 call timer( 'PHION', 1 ) 1074C On first call set the logical DirectPhi 1075 if (firsttime) then 1076 DirectPhi = fdf_get('DirectPhi', .false. ) 1077 firsttime = .false. 1078 else 1079 ! reset the data structures in module meshphi 1080 call resetMeshPhi( ) 1081 endif 1082 1083C Allocate memory related to nmpl 1084 nullify( endpht, phi, lstpht, listp2 ) 1085 call re_alloc( endpht, 0, nmpl, 'endpht', 'PhiOnMesh' ) 1086 1087C Initialise pointer array 1088 endpht(0) = 0 1089 do ip = 1,nmpl 1090 endpht(ip) = endpht(ip-1) + numphi(ip) 1091 enddo 1092 1093C Allocate phi if this is not a direct calculation 1094 nlist = endpht(nmpl) 1095 if (DirectPhi) then 1096 nphi = 1 1097 else 1098 nphi = nlist 1099 if (IONode) then 1100 write(6,"(a,i20)") 1101 $ "PhiOnMesh: Number of (b)points on node 0 = ", nmpl 1102 write(6,"(a,i20)") 1103 $ "PhiOnMesh: nlist on node 0 = ", nlist 1104 endif 1105 endif 1106 1107C Add an extra margin of error to nlist to minimise reallocations 1108 nliste = 1.01*nlist 1109 1110C Adjust dimensions of phi if necessary 1111 if (associated(phi)) then 1112 nmax = size(phi,2) 1113 else 1114 nmax = 0 1115 end if 1116 if (nphi.gt.nmax) then 1117 if (DirectPhi) then 1118 call re_alloc( phi, 1, nsp, 1, nphi, 'phi', 'PhiOnMesh' ) 1119 else 1120 call re_alloc( phi, 1, nsp, 1, nliste, 'phi', 'PhiOnMesh' ) 1121 endif 1122 endif 1123 1124C Adjust dimensions of list arrays if necessary 1125 if (associated(lstpht) .and. associated(listp2)) then 1126 nmax = min( size(lstpht), size(listp2) ) 1127 else 1128 nmax = 0 1129 end if 1130 if (nlist.gt.nmax) then 1131 call re_alloc( lstpht, 1, nliste, 'lstpht', 'PhiOnMesh' ) 1132 call re_alloc( listp2, 1, nliste, 'listp2', 'PhiOnMesh' ) 1133 endif 1134 1135C Find indexes and values of atomic orbitals at mesh points 1136 numphi = 0 1137 do io = 1,norb 1138 ia = iaorb(io) 1139 if (ipa(ia)==0) cycle ! Atomic sphere does not intersect my Box 1140 iphi = iphorb(io) 1141 is = isa(ia) 1142 r2o = rcut(is,iphi)**2 1143 do iop = 1, mop 1144 ip0 = indexp( ipa(ia) + idop(iop) ) 1145 if (ip0 .gt. 0) then 1146 within = .false. 1147 do isp = 1,nsp 1148 dxsp(1:3,isp) = xdop(1:3,iop)+xdsp(1:3,isp)-dxa(1:3,ia) 1149 r2sp(isp) = sum( dxsp(1:3,isp)**2 ) 1150 if (r2sp(isp) .lt. r2o) within = .true. 1151 enddo 1152 if (within) then 1153 numphi(ip0) = numphi(ip0) + 1 1154 n = endpht(ip0-1) + numphi(ip0) 1155 lstpht(n) = io 1156 listp2(n) = iop 1157 if (.not.DirectPhi) then 1158 do isp = 1,nsp 1159 if (r2sp(isp) .lt. r2o) then 1160 call phiatm( is, iphi, dxsp(1,isp), 1161 & phip, grphi ) 1162 phi(isp,n) = phip 1163 else 1164 phi(isp,n) = 0.0_grid_p 1165 endif 1166 enddo 1167 endif 1168 endif 1169 endif 1170 enddo 1171 enddo 1172 1173 call timer( 'PHION', 2 ) 1174#ifdef DEBUG 1175 call write_debug( ' POS PhiOnMesh' ) 1176#endif 1177!------------------------------------------------------------------------ END 1178 end subroutine PhiOnMesh 1179 1180 subroutine setupExtMesh( distr, rmax ) 1181! ---------------------------------------------------------------------- 1182! Setup the extended mesh variables for a given data distribution. 1183! It creates and initialise indexp, xdop and idop. These arrays are 1184! allocated inside a meshDisType structure and can be accesed by pointers 1185! of the mesh module. 1186! 1187! Only necesary for those data distributions who needs to use the 1188! extended mesh. 1189! ---------------------------------------------------------------------- 1190! Input : 1191! ---------------------------------------------------------------------- 1192! integer distr : distribution number 1193! real(dp) rmax : Maximum orbital radius 1194! 1195! ---------------------------------------------------------------------- 1196! Output : 1197! ---------------------------------------------------------------------- 1198! All output quantities are in the mesh and moreMeshSubs modules. 1199! Data is stored in moreMeshSubs but accessed using the mesh module 1200! pointers. 1201! 1202!*********************************************************************** 1203 use precision, only : dp 1204 use parallel, only : IONode 1205 use mesh, only : meshLim ! My processor's box of mesh points 1206 use mesh, only : ne ! Points in rmax along each lat. vector 1207 use mesh, only : nem ! Extended-mesh divisions in each direction 1208 use mesh, only : nmsc ! Mesh points in each supercell vector 1209 use mesh, only : mop ! Accumulated num. of orbital points 1210 use mesh, only : nsp ! Number of sub-points of each mesh point 1211 use mesh, only : cmesh ! Mesh cell vectors 1212 use mesh, only : xdsp ! Vector to mesh sub-points 1213 use mesh, only : indexp ! ranslation from extended to 1214 ! normal mesh index 1215 use mesh, only : idop ! Mesh-index span of points 1216 ! within an atomic sphere 1217 use mesh, only : xdop ! Coordinates of (super)points 1218 ! within an atomic sphere 1219 use moreMeshSubs, only : allocExtMeshDistr 1220 1221 implicit none 1222C Input parameters 1223 integer :: distr 1224 real(dp) :: rmax 1225C Local parameters 1226 integer :: i1, i2, i3, j1, j2, j3, k1, k2, k3, 1227 & i, j, k, nep, boxWidth(3), extWidth(3), 1228 & isp, myBox(2,3), myExtBox(2,3) 1229 real(dp) :: r, dxp(3), dx(3) 1230 logical :: within 1231C Functions 1232 real(dp) :: dismin 1233 1234 ! Correspondence between JMS and RG box descriptors 1235 myBox(:,:) = meshLim(:,:) - 1 1236 1237! Add 'wings extensions' to myBox of mesh points, containing all points 1238! within the orbital spheres that may intersect myBox. Thus, wings must 1239! be one diameter wide. 1240 myExtBox(1,:) = myBox(1,:) - 2*ne(:) 1241 myExtBox(2,:) = myBox(2,:) + 2*ne(:) 1242 1243! Find number of extended-box points. 1244 nem(:) = myExtBox(2,:) - myExtBox(1,:) + 1 1245 nep = nem(1) * nem(2) * nem(3) 1246 1247 if (IONode) then 1248 write(6,'(a,3(i6,a),i12)') 'ExtMesh (bp) on 0 =', 1249 & nem(1),' x',nem(2),' x',nem(3),' =', nep 1250 endif 1251 1252! Allocate memory for indexp(nep), idop(mop) and xdop(3,mop) of 1253! the current data distribution 1254 call allocExtMeshDistr( distr, nep, mop ) 1255 1256! Find relationship between extended and normal box points 1257 boxWidth(:) = myBox(2,:) - myBox(1,:) + 1 1258 extWidth(:) = myExtBox(2,:) - myExtBox(1,:) + 1 1259 do i3 = myExtBox(1,3),myExtBox(2,3) 1260 do i2 = myExtBox(1,2),myExtBox(2,2) 1261 do i1 = myExtBox(1,1),myExtBox(2,1) 1262 1263 ! Find periodic indexes within first supercell 1264 j1 = modulo( i1, nmsc(1) ) 1265 j2 = modulo( i2, nmsc(2) ) 1266 j3 = modulo( i3, nmsc(3) ) 1267 1268 ! Find indexes relative to box origins 1269 j1 = j1 - myBox(1,1) 1270 j2 = j2 - myBox(1,2) 1271 j3 = j3 - myBox(1,3) 1272 k1 = i1 - myExtBox(1,1) 1273 k2 = i2 - myExtBox(1,2) 1274 k3 = i3 - myExtBox(1,3) 1275 1276 ! Find combined mesh indexes. 1277 j = 1 + j1 + boxWidth(1)*j2 + boxWidth(1)*boxWidth(2)*j3 1278 k = 1 + k1 + extWidth(1)*k2 + extWidth(1)*extWidth(2)*k3 1279 1280 ! Find myExtBox -> myBox index translation 1281 if (j1>=0 .and. j1<boxWidth(1) .and. 1282 . j2>=0 .and. j2<boxWidth(2) .and. 1283 . j3>=0 .and. j3<boxWidth(3)) then 1284 indexp(k) = j ! Point is within myBox 1285 else 1286 indexp(k) = 0 ! Point is outside myBox 1287 endif ! (j>=myBox(1) and j<=myBox(2)) 1288 1289 enddo ! i1 1290 enddo ! i2 1291 enddo ! i3 1292 1293! This loop is done for each distribution, and replicates 1294! the one in InitMesh, except that now we actually fill 1295! the data in the arrays. Note that mop does not depend 1296! on the distribution. 1297 1298! Find points within rmax (orbital points) 1299 mop = 0 1300 do i3 = -ne(3), ne(3) 1301 do i2 = -ne(2), ne(2) 1302 do i1 = -ne(1), ne(1) 1303 dxp(:) = cmesh(:,1) * i1 + 1304 & cmesh(:,2) * i2 + 1305 & cmesh(:,3) * i3 1306 within = .false. 1307 do isp = 1,nsp 1308 dx(1:3) = dxp(1:3) + xdsp(1:3,isp) 1309 r = dismin( cmesh, dx ) 1310 if ( r .lt. rmax ) within = .true. 1311 enddo 1312 if ( within ) then 1313 mop = mop + 1 1314 ! Store index-distance and vector-distance to point. 1315 idop(mop) = i1 + nem(1)*i2 + nem(1)*nem(2)*i3 1316 xdop(1:3,mop) = dxp(1:3) 1317 endif ! (within) 1318 enddo ! i1 1319 enddo ! i2 1320 enddo ! i3 1321 1322 end subroutine setupExtMesh 1323 1324!------------------------------------------------------------------ 1325 subroutine distriPhiOnMesh( nm, nmpl, norb, iaorb, 1326 & iphorb, isa, numphi, 1327 & need_gradients_in_xc ) 1328 1329C Computes the number of orbitals that intersects every mesh point 1330C and calls initMeshDistr in order to compute a new data distribution 1331C based in the number of orbital intersections. 1332C ================================================================== 1333C 1334C INPUT: 1335C integer nm(3) : Number of Mesh divisions of each cell vector 1336C integer nmpl : Local number of mesh points 1337C integer norb : Total number of basis orbitals in supercell 1338C integer iaorb(norb) : Atom to which each orbital belongs 1339C integer iphorb(norb) : Orbital index (within atom) of each orbital 1340C integer isa(na) : Species index of all atoms in supercell 1341C logical need_gradients_in_xc : 1342C OUTPUT: 1343C The output values are in the module moreMeshSubs (New limits of the mesh 1344C and the information needed to transfer data between data 1345C distributions) 1346C 1347C BEHAVIOR: 1348C The amount of computation in every point of the mesh depends on the 1349C number of orbitals that intersecs in that mesh point. This can lead 1350C to unbalanced distributions in parallel executions. We want to 1351C compute a new data distribution, where the load is balanced among 1352C the several processes. 1353C 1354C In order to compute the number of orbitals that intersects in every 1355C mesh point, we loop in all orbitals of every atom and we increase the 1356C value of numphi in every point of the mesh that is inside of the 1357C radius of the current orbital. 1358C 1359C Finally, we call initMeshDistr that computes the new data distributions 1360C that will be used inside setup_hamiltonian. We are going to have 3 1361C different data distributions: 1362C - Uniform: Same amount of mesh in every process. Used in Poison 1363C - Quadratic: The load in every point of the mesh is proporcional 1364C to the square of the number of orbitals that intersects 1365C in that point. Used in rhoofd and vmat. 1366C - Lineal: We try to distribute those points of the mesh that have 1367C at least one orbital intersection. Used in cellxc. 1368C 1369C ================================================================== 1370 use mesh, only: mop, ipa, idop, nsp, xdop, xdsp, dxa, 1371 & indexp, meshLim 1372 use precision, only: dp, grid_p 1373 use atmfuncs, only: rcut 1374 use parallel, only: Nodes, node 1375#ifdef MPI 1376 use moreMeshSubs, only: initMeshExtencil 1377#endif 1378 use moreMeshSubs, only: initMeshDistr 1379 use moreMeshSubs, only: allocASynBuffer 1380 use moreMeshSubs, only: UNIFORM, QUADRATIC, LINEAR 1381 use alloc, only: re_alloc, de_alloc 1382 1383 implicit none 1384C Passed arguments 1385 integer, intent(in) :: nm(3), nmpl, norb, iaorb(norb), 1386 & iphorb(norb), isa(*) 1387 integer, intent(out) :: numphi(nmpl) 1388 logical, intent(in) :: need_gradients_in_xc 1389C Local variables 1390 integer :: ii, io, ia, iphi, is, iop, ip0, isp 1391 real(dp) :: r2o, dxsp(3,nsp), r2sp(nsp) 1392 logical :: within 1393 integer, pointer :: wload(:) 1394!------------------------------------------------------------------------- BEGIN 1395#ifdef DEBUG 1396 call write_debug( ' PRE distriPhiOnMesh' ) 1397#endif 1398 call timer( 'REMESH', 1 ) 1399 1400 do ii= 1, nmpl 1401 numphi(ii) = 0 1402 enddo 1403 1404C Find number of atomic orbitals at mesh points 1405 do io = 1,norb 1406 ia = iaorb(io) 1407 if (ipa(ia)==0) cycle ! Atomic sphere does not intersec my Box 1408 iphi = iphorb(io) 1409 is = isa(ia) 1410 r2o = rcut(is,iphi)**2 1411C Loop over mesh points inside rmax 1412 do iop = 1, mop 1413 ip0 = indexp( ipa(ia) + idop(iop) ) 1414 if (ip0 .gt. 0) then 1415C Loop over sub-points to find if point is within range 1416 within = .false. 1417 do isp = 1,nsp 1418 dxsp(1:3,isp) = xdop(1:3,iop)+xdsp(1:3,isp)-dxa(1:3,ia) 1419 r2sp(isp) = sum( dxsp(1:3,isp)**2 ) 1420 if (r2sp(isp) .lt. r2o) within = .true. 1421 enddo 1422C If within range, add one to number of point orbitals 1423 if (within) numphi(ip0) = numphi(ip0) + 1 1424 endif 1425 enddo 1426 enddo 1427 1428#ifdef MPI 1429 if (nodes.gt.1) then 1430 nullify( wload ) 1431 call re_alloc( wload, 1, nmpl, 'wload', 'distriPhiOnMesh' ) 1432 do ip0= 1, nmpl 1433 wload(ip0) = numphi(ip0)*(numphi(ip0)+1) + 1 1434 enddo 1435 call initMeshDistr( UNIFORM, QUADRATIC, nm, wload ) 1436 1437 do ip0= 1, nmpl 1438! The computational work of a mesh point in cellxc is almost 0 when 1439! when the charge density is zero. We use a proportion of 400/1. 1440 if (numphi(ip0).ne.0) then 1441 wload(ip0) = 400 1442 else 1443 wload(ip0) = 1 1444 endif 1445 enddo 1446 call initMeshDistr( UNIFORM, LINEAR, nm, wload ) 1447 if (need_gradients_in_xc) call initMeshExtencil( LINEAR, nm ) 1448 1449C Free local memory 1450 call de_alloc( wload, 'wload', 'distriPhiOnMesh' ) 1451 1452C Allocate memory for asynchronous communications. 1453C It does nothing if we use synchronous communications. 1454 call allocASynBuffer( LINEAR ) 1455 endif 1456 1457#endif 1458 call timer( 'REMESH', 2 ) 1459#ifdef DEBUG 1460 call write_debug( ' POS distriPhiOnMesh' ) 1461#endif 1462!--------------------------------------------------------------------------- END 1463 end subroutine distriPhiOnMesh 1464 1465 end module meshsubs 1466