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