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 setspatial(na,xa,vecs,rcut,lspatialok) 9C 10C Sets up the spatial decomposition of the atoms into boxes 11C 12C Based on the equivalent routine from GULP, except that it 13C now uses fractional coordinates to handle low symmetry cells. 14C 15C On entry : 16C 17C na = number of atoms in cell 18C xa(3,na) = Cartesian coordinates of atoms 19C vecs(3,3) = cell vectors 20C rcut = maximum potential cut-off radius 21C Node = local Node number 22C Nodes = total number of processors 23C 24C On exit : 25C 26C lspatialok = if .true. then the cell can be decomposed 27C into boxes 28C 29C Set in modules: 30C 31C nspcellat(ncell) = number of atoms per spatial cell 32C nspcellatptr(maxspcellat, ncell) = pointer to global atom number from atom number in cell 33C nspcellatptrcell = pointer to image (1-27) of atom in cell 34C spmin(3) = min spatial cell length in x, y, z 35C spcell = 36C nspcell(3) = no. spatial cells in x, y, z 37C maxspcell = 38C maxspcellat = maxval(nspcellat) + a bit extra 39C 40C nspcellat(ncell) = number of atoms per cell 41C nspcellatptr(natom,ncell) = pointer to atom from atom number in cell 42C nspcellatptrcell(natom,ncell) = pointer to image (1-27) of atom in cell 43 44C 45C 46C Julian Gale, NRI, Curtin University, March 2004 47C 48 use parallel, only : Node 49 use precision, only : dp 50 use spatial, only : nspcellat, nspcellatptr, nspcellatptrcell 51 . , nbufferx, nbuffery, nbufferz 52 . , spmin, spcell, nspcell 53 . , maxspcell, maxspcellat 54 use alloc, only : re_alloc, de_alloc 55 implicit none 56C 57C Passed variables 58C 59 integer, intent(in) :: na 60 logical, intent(out) :: lspatialok 61 real(dp), intent(in) :: xa(3,na) 62 real(dp), intent(in) :: vecs(3,3) 63 real(dp), intent(in) :: rcut 64C 65C Local variables 66C 67 integer :: i 68 integer :: ii 69 integer :: ind 70 integer :: ixvec1cell(27) 71 integer :: iyvec1cell(27) 72 integer :: izvec1cell(27) 73 integer :: maxxy 74 integer :: maxx 75 integer :: n 76 integer :: n2bufferplus1x 77 integer :: n2bufferplus1y 78 integer :: n2bufferplus1z 79 integer :: np 80 integer :: npc 81 integer :: nx 82 integer :: ny 83 integer :: nz 84 integer :: nspcelltot 85 logical :: ldebug 86 logical, save :: firstcall = .true. 87 real(dp) :: a 88 real(dp) :: b 89 real(dp) :: c 90 real(dp) :: alpha 91 real(dp) :: beta 92 real(dp) :: gamma 93 real(dp) :: degtorad 94 real(dp) :: small 95 real(dp) :: xdivision 96 real(dp) :: ydivision 97 real(dp) :: zdivision 98 real(dp), pointer :: xyzfrac(:,:) 99 real(dp) :: xi 100 real(dp) :: yi 101 real(dp) :: zi 102 real(dp) :: xvec1cell(27) 103 real(dp) :: yvec1cell(27) 104 real(dp) :: zvec1cell(27) 105 real(dp), dimension(:), pointer :: xinbox 106 real(dp), dimension(:), pointer :: yinbox 107 real(dp), dimension(:), pointer :: zinbox 108C 109C Set the following according to whether debugging is required 110C 111c ldebug = .false. 112 ldebug = (Node.eq.0) 113C 114C On first call initialise spatial decomposition arrays 115C 116 if (firstcall) then 117 nullify(nspcellat) 118 nullify(nspcellatptr) 119 nullify(nspcellatptrcell) 120 nullify(xinbox) 121 nullify(yinbox) 122 nullify(zinbox) 123 firstcall = .false. 124 endif 125C 126C Allocate global memory 127C 128 call re_alloc(xinbox,1,na,name='xinbox') 129 call re_alloc(yinbox,1,na,name='yinbox') 130 call re_alloc(zinbox,1,na,name='zinbox') 131C 132C Allocate local workspace 133C 134 nullify(xyzfrac) 135 call re_alloc(xyzfrac,1,3,1,na,name="xyzfrac") 136C 137 small = 0.000001d0 138 degtorad = 4.0d0*atan(1.0d0)/180.0d0 139C 140C Find cell parameters 141C 142 call uncell(vecs,a,b,c,alpha,beta,gamma,degtorad) 143C 144C Check that cut-off is greater than zero 145C 146 lspatialok = .true. 147 if (rcut.lt.1.0d-10) lspatialok = .false. 148C 149C If cell is not compatible return 150C 151 if (.not.lspatialok) return 152C 153C Copy current atomic coordinates into box coordinate arrays 154C 155 xinbox(1:na) = xa(1,1:na) 156 yinbox(1:na) = xa(2,1:na) 157 zinbox(1:na) = xa(3,1:na) 158C 159C Find cell lengths in each direction 160C 161 spmin(1:3) = - small 162 spcell(1) = a 163 spcell(2) = b 164 spcell(3) = c 165C 166C Find numbers of cells in each direction - subtract a small amount for safety 167C Add 2*nbuffer to number of cells in each direction to allow for buffer on each side 168C 169 nspcell(1) = 2*nbufferx + (spcell(1) - 1.0d-6)/rcut 170 nspcell(2) = 2*nbuffery + (spcell(2) - 1.0d-6)/rcut 171 nspcell(3) = 2*nbufferz + (spcell(3) - 1.0d-6)/rcut 172C 173C Check that minimum number of cells in any direction is at least 2*nbuffer + 1 174C 175 n2bufferplus1x = 2*nbufferx + 1 176 n2bufferplus1y = 2*nbuffery + 1 177 n2bufferplus1z = 2*nbufferz + 1 178 nspcell(1) = max(nspcell(1),n2bufferplus1x) 179 nspcell(2) = max(nspcell(2),n2bufferplus1y) 180 nspcell(3) = max(nspcell(3),n2bufferplus1z) 181C 182C Set up cell images 183C 184 call cellimagelist(vecs,xvec1cell,yvec1cell,zvec1cell, 185 . ixvec1cell,iyvec1cell,izvec1cell) 186C 187C Create an array of fractional coordinates to decide which box each atom belongs to 188C and place Cartesian coordinates within central cell 189C 190 call cart2frac(na,xinbox,yinbox,zinbox,vecs,xyzfrac) 191C 192C Calculate the total number of cells and check that arrays are allocated to this size 193C 194 nspcelltot = nspcell(1)*nspcell(2)*nspcell(3) 195 if (nspcelltot.gt.maxspcell) then 196 maxspcell = nspcelltot 197 call re_alloc(nspcellat,1,maxspcell,name='nspcellat') 198 endif 199C 200C Initialise atom counters 201C 202 nspcellat(1:nspcelltot) = 0 203C 204C Find which box each atom belongs in. This is described by the arrays: 205C 206C nspcellat(ncell) = number of atoms per cell 207C nspcellatptr(natom,ncell) = pointer to atom from atom number in cell 208C nspcellatptrcell(natom,ncell) = pointer to image (1-27) of atom in cell 209C 210C Have to loop over cell images in order to fill buffer regions too 211C 212 xdivision = nspcell(1) - 2*nbufferx 213 ydivision = nspcell(2) - 2*nbuffery 214 zdivision = nspcell(3) - 2*nbufferz 215 maxxy = nspcell(1)*nspcell(2) 216 maxx = nspcell(1) 217 do i = 1,na 218 do ii = 1,27 219 nx = (xyzfrac(1,i) + ixvec1cell(ii) - spmin(1))*xdivision + 220 . 2*nbufferx 221 if (nx.ge.1.and.nx.le.nspcell(1)) then 222 ny = (xyzfrac(2,i) + iyvec1cell(ii) - spmin(2))*ydivision + 223 . 2*nbuffery 224 if (ny.ge.1.and.ny.le.nspcell(2)) then 225 nz = (xyzfrac(3,i) + izvec1cell(ii) - spmin(3))*zdivision+ 226 . 2*nbufferz 227 if (nz.ge.1.and.nz.le.nspcell(3)) then 228 ind = (nz-1)*maxxy + (ny-1)*maxx + nx 229 nspcellat(ind) = nspcellat(ind) + 1 230 if (nspcellat(ind).gt.maxspcellat) then 231 maxspcellat = nspcellat(ind) + 50 232 call re_alloc(nspcellatptr,1,maxspcellat, 233 . 1,maxspcell,name='nspcellatptr') 234 call re_alloc(nspcellatptrcell,1,maxspcellat, 235 . 1,maxspcell,name='nspcellatptrcell') 236 endif 237 nspcellatptr(nspcellat(ind),ind) = i 238 nspcellatptrcell(nspcellat(ind),ind) = ii 239 endif 240 endif 241 endif 242 enddo 243 enddo 244C 245C Debugging information 246C 247 if (ldebug) then 248 write(6,'(/,'' Spatial decomposition data : '',/)') 249 write(6,'('' Cutoff = '',f8.4,'' Angstroms'',/)') 250 . rcut*0.529177d0 251 write(6,'('' Spatial decomposition cells (containing atoms) :'' 252 . ,/)') 253 write(6,'(6x,''Ncells'',4x,''Cell min'',4x,''Cell length'')') 254 write(6,'('' x : '',i6,2(3x,f9.4))') nspcell(1),spmin(1), 255 . spcell(1) 256 write(6,'('' y : '',i6,2(3x,f9.4))') nspcell(2),spmin(2), 257 . spcell(2) 258 write(6,'('' z : '',i6,2(3x,f9.4))') nspcell(3),spmin(3), 259 . spcell(3) 260 write(6,'(/)') 261 write(6,'(''----------------------------------------------'', 262 . ''----------------------------------'')') 263 write(6,'('' Cell No. of atoms Atom No. x '', 264 . '' y z'')') 265 write(6,'(''----------------------------------------------'', 266 . ''----------------------------------'')') 267 do i = 1,nspcelltot 268 if (nspcellat(i).gt.0) then 269 np = nspcellatptr(1,i) 270 npc = nspcellatptrcell(1,i) 271 xi = xinbox(np) + xvec1cell(npc) 272 yi = yinbox(np) + yvec1cell(npc) 273 zi = zinbox(np) + zvec1cell(npc) 274 write(6,'(i6,4x,i8,4x,i8,3(1x,f12.4))') i,nspcellat(i), 275 . np,xi,yi,zi 276 do n = 2,nspcellat(i) 277 np = nspcellatptr(n,i) 278 npc = nspcellatptrcell(n,i) 279 xi = xinbox(np) + xvec1cell(npc) 280 yi = yinbox(np) + yvec1cell(npc) 281 zi = zinbox(np) + zvec1cell(npc) 282 write(6,'(22x,i8,3(1x,f12.4))') np,xi,yi,zi 283 enddo 284 write(6,'(''--------------------------------------------'', 285 . ''------------------------------------'')') 286 endif 287 enddo 288 write(6,'(/)') 289 endif 290C 291C Deallocate local workspace 292C 293 call de_alloc(xyzfrac,name='xyzfrac') 294 call de_alloc(xinbox,name='xinbox') 295 call de_alloc(yinbox,name='yinbox') 296 call de_alloc(zinbox,name='zinbox') 297C 298 return 299 end 300