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!> Utilities to dea with real-space grid functions
10module m_gridfunc
11
12  implicit          none
13
14  integer, parameter, private :: dp = selected_real_kind(10,100)
15  integer, parameter, private :: sp = selected_real_kind(5,10)
16  integer, parameter          :: grid_p = sp  ! NOTE this
17
18  type, public :: gridfunc_t
19     !> Lattice vectors (by columns, in bohr)
20     real(dp)              ::  cell(3,3)
21     !> Possible shift for the origin of the box, useful in some
22     !> situations (sub-boxes, for example) (in bohr)
23     real(dp)              ::  origin(3) = [0.0_dp, 0.0_dp, 0.0_dp]
24     !> The 'grid' format can in principle hold functions that
25     !> are not periodic (for example, values in a sub-box)
26     logical               ::  is_periodic(3) = [ .true., .true., .true.]
27     !> Number of points along each cell direction
28     integer               ::  n(3) = [ 0,0,0 ]
29     !> Extra dimension, typically the spin index
30     integer               ::  nspin
31     !> Array holding the values of the grid function
32     real(grid_p), allocatable ::  val(:,:,:,:)
33  end type gridfunc_t
34
35      public :: write_gridfunc, read_gridfunc, clean_gridfunc
36#ifdef CDF
37      public :: read_gridfunc_netcdf, write_gridfunc_netcdf
38#endif
39      public :: get_planar_average, grid_p, monoclinic_z
40      public :: get_values_at_point
41      public :: monoclinic
42
43      private
44
45 CONTAINS
46
47   subroutine clean_gridfunc(gf)
48     type(gridfunc_t), intent(inout) :: gf
49
50     if (allocated(gf%val)) then
51        deallocate(gf%val)
52     endif
53     gf%n(:) = 0
54     gf%cell(:,:) = 0.0_dp
55     gf%origin(:) = 0.0_dp
56     gf%is_periodic(:) = .true.
57
58   end subroutine clean_gridfunc
59
60   subroutine write_gridfunc(gf,fname)
61     type(gridfunc_t), intent(in) :: gf
62     character(len=*), intent(in) :: fname
63
64     integer :: n(3)
65     integer :: isp, ix, iy, iz
66     integer :: iu
67
68     n = gf%n
69
70     call get_lun(iu)
71     open(unit=iu,file=fname,form="unformatted",status="unknown", &
72          position="rewind",action="write")
73     write(iu) gf%cell, gf%origin, gf%is_periodic
74     write(iu) n, gf%nspin
75     do isp=1,gf%nspin
76        do iz=1,n(3)
77           do iy=1,n(2)
78              write(iu) (gf%val(ix,iy,iz,isp),ix=1,n(1))
79           enddo
80        enddo
81     enddo
82     close( iu )
83
84   end subroutine write_gridfunc
85
86   !> Evaluates a 'grid function' at an arbitrary point in the box.
87   !> The function is computed using linear interpolation.
88   subroutine get_values_at_point(gf, xfrac, vals)
89     !> The object representing the grid function
90     type(gridfunc_t), intent(in) :: gf
91     !> The fractional coordinates of the point
92     real(dp), intent(in) :: xfrac(3)
93     !> Array holding the values (with the 'spin' dimension)
94     real(grid_p), allocatable :: vals(:)
95
96     integer  ::  n(3), lo(3), hi(3)
97     real(dp) ::  r(3), x(3), y(3)
98     integer  ::  i, j, k
99     real(dp) ::  nk
100
101     n(:) = gf%n(:)
102
103!           Find the right 3D "grid cube" and the reduced coordinates
104!           of the point in it. The double mod assures that negative
105!           numbers are well treated (the idea is to bring the coordinates
106!           to the [0,n(k)) interval)
107!
108     do k = 1, 3
109        nk = real(n(k),kind=dp)
110        r(k) =  modulo(n(k)*xfrac(k),nk)
111        lo(k) = int(r(k))
112        hi(k) = mod ( lo(k)+1, n(k) )
113        x(k) = r(k) - lo(k)
114        y(k) = 1 - x(k)
115     enddo
116
117     ! Switch to 1-based array convention
118     lo(:) = lo(:) + 1
119     hi(:) = hi(:) + 1
120
121!      compute charge density by linear interpolation
122
123     vals(:) = gf%val(lo(1),lo(2),lo(3),:) * y(1) * y(2) * y(3) + &
124          gf%val(lo(1),lo(2),hi(3),:) * y(1) * y(2) * x(3) + &
125          gf%val(lo(1),hi(2),lo(3),:) * y(1) * x(2) * y(3) + &
126          gf%val(lo(1),hi(2),hi(3),:) * y(1) * x(2) * x(3) + &
127          gf%val(hi(1),lo(2),lo(3),:) * x(1) * y(2) * y(3) + &
128          gf%val(hi(1),lo(2),hi(3),:) * x(1) * y(2) * x(3) + &
129          gf%val(hi(1),hi(2),lo(3),:) * x(1) * x(2) * y(3) + &
130          gf%val(hi(1),hi(2),hi(3),:) * x(1) * x(2) * x(3)
131
132   end subroutine get_values_at_point
133
134
135   !> Computes a simple-minded 'planar average' of a grid
136   !> function. This is a mathematical average. For a more
137   !> physical average, see the 'macroave' utility
138   subroutine get_planar_average(gf, dim, val)
139     !> The object representing the grid function
140     type(gridfunc_t), intent(in) :: gf
141     !> The dimension that indexes the average values.
142     !> These are obtained by summing over the two other dimensions.
143     integer, intent(in) :: dim
144     !> Array holding the averages (second dimension is the 'spin'))
145     real(grid_p), allocatable :: val(:,:)
146
147     integer :: isp, ix, iy, iz, n(3), spin_dim
148     real(grid_p) :: sum
149
150     n = gf%n
151     spin_dim = size(gf%val,dim=4)
152
153     allocate(val(n(dim),spin_dim))
154
155     do isp=1,gf%nspin
156        select case (dim)
157        case (3)
158           do iz=1,n(3)
159              sum = 0.0_grid_p
160              do iy=1,n(2)
161                 do ix=1,n(1)
162                    sum = sum + gf%val(ix,iy,iz,isp)
163                 enddo
164              enddo
165              val(iz,isp) = sum/(n(1)*n(2))
166           enddo
167        case (2)
168           do iy=1,n(2)
169              sum = 0.0_grid_p
170              do iz=1,n(3)
171                 do ix=1,n(1)
172                    sum = sum + gf%val(ix,iy,iz,isp)
173                 enddo
174              enddo
175              val(iy,isp) = sum/(n(1)*n(3))
176           enddo
177        case (1)
178           do ix=1,n(1)
179              sum = 0.0_grid_p
180              do iz=1,n(3)
181                 do iy=1,n(2)
182                    sum = sum + gf%val(ix,iy,iz,isp)
183                 enddo
184              enddo
185              val(ix,isp) = sum/(n(2)*n(3))
186           enddo
187        end select
188
189     enddo
190
191   end subroutine get_planar_average
192
193   subroutine read_gridfunc(fname,gf)
194     type(gridfunc_t), intent(inout) :: gf
195     character(len=*), intent(in) :: fname
196
197     integer :: n(3)
198     integer :: isp, ix, iy, iz
199     integer :: iu, iostat
200
201     call clean_gridfunc(gf)
202
203     call get_lun(iu)
204     open(unit=iu,file=fname,form="unformatted",status="old", &
205          position="rewind",action="read")
206
207     read(iu,iostat=iostat) gf%cell, gf%origin, gf%is_periodic
208     if (iostat /= 0) then
209        backspace(iu)
210        read(iu) gf%cell
211        gf%origin(:) = 0.0_dp
212        gf%is_periodic(:) = .true.
213     endif
214
215     read(iu) gf%n, gf%nspin
216     n = gf%n
217     allocate(gf%val(n(1),n(2),n(3),gf%nspin))
218
219     do isp=1,gf%nspin
220        do iz=1,n(3)
221           do iy=1,n(2)
222              read(iu) (gf%val(ix,iy,iz,isp),ix=1,n(1))
223           enddo
224        enddo
225     enddo
226     close( iu )
227
228   end subroutine read_gridfunc
229
230#ifdef CDF
231
232subroutine write_gridfunc_netcdf(gf,filename,description)
233use netcdf
234
235implicit none
236
237character(len=*), intent(in) :: filename
238type(gridfunc_t), intent(in)    :: gf
239character(len=*), intent(in), optional :: description
240
241integer  :: ncid
242integer  :: xyz_id, step_id, abc_id, spin_id
243integer  :: cell_id, n1_id, n2_id, n3_id, gridfunc_id
244integer  :: origin_id, is_periodic_id
245
246integer   ::    nspin  ! Number of spins
247integer   ::    n(3), int_from_bool(3), i
248integer   ::   ispin, iostat, ix, iy, iz, iret
249
250!-----------------------------------------------------
251
252call check( nf90_create(filename,NF90_CLOBBER,ncid))
253       iret = nf90_def_dim(ncid,'xyz',3,xyz_id)
254       call check(iret)
255       iret = nf90_def_dim(ncid,'abc',3,abc_id)
256       call check(iret)
257       iret = nf90_def_dim(ncid,'spin',gf%nspin,spin_id)
258       call check(iret)
259
260       n(:) = gf%n(:)
261
262       iret = nf90_def_dim(ncid,'n1',n(1),n1_id)
263       call check(iret)
264       iret = nf90_def_dim(ncid,'n2',n(2),n2_id)
265       call check(iret)
266       iret = nf90_def_dim(ncid,'n3',n(3),n3_id)
267       call check(iret)
268
269       iret = nf90_def_var(ncid,'cell',nf90_float,(/xyz_id,abc_id/),cell_id)
270       call check(iret)
271       iret = nf90_put_att(ncid,cell_id,'Description', &
272               "Cell vectors in Bohr: xyz, abc")
273       call check(iret)
274
275       iret = nf90_def_var(ncid,'origin',nf90_float,(/xyz_id/),origin_id)
276       call check(iret)
277       iret = nf90_put_att(ncid,origin_id,'Description', &
278               "Origin of coordinates: xyz")
279       call check(iret)
280
281       iret = nf90_def_var(ncid,'is_periodic',nf90_int,(/abc_id/),is_periodic_id)
282       call check(iret)
283       iret = nf90_put_att(ncid,cell_id,'Description', &
284               "Periodic along cell vectors (1:yes, 0:no): abc")
285       call check(iret)
286
287       iret = nf90_def_var(ncid,'gridfunc',nf90_float,(/n1_id,n2_id,n3_id,spin_id/),gridfunc_id)
288       call check(iret)
289       if (present(description)) then
290          iret = nf90_put_att(ncid,gridfunc_id,'Description', &
291               trim(description))
292       else
293          iret = nf90_put_att(ncid,gridfunc_id,'Description', &
294               "Grid function -- ")
295       endif
296       call check(iret)
297
298       iret = nf90_enddef(ncid)
299       call check(iret)
300!
301       iret = nf90_put_var(ncid, cell_id, gf%cell, start = (/1, 1 /), count = (/3, 3/) )
302       call check(iret)
303       iret = nf90_put_var(ncid, origin_id, gf%origin, start = [1], count = [3] )
304       call check(iret)
305       int_from_bool(:) = 0
306       do i = 1, 3
307          if (gf%is_periodic(i)) int_from_bool(i) = 1
308       enddo
309       iret = nf90_put_var(ncid, is_periodic_id, int_from_bool, start = [1], count = [3] )
310       call check(iret)
311
312      iret = nf90_put_var(ncid, gridfunc_id, gf%val, start = (/1, 1, 1, 1 /), &
313           count = (/n(1), n(2), n(3), nspin/) )
314      call check(iret)
315
316     call check( nf90_close(ncid) )
317
318end subroutine write_gridfunc_netcdf
319!--------
320
321subroutine read_gridfunc_netcdf(filename,gf)
322use netcdf
323
324implicit none
325
326character(len=*), intent(in) :: filename
327type(gridfunc_t), intent(inout)    :: gf
328
329integer  :: ncid
330integer  :: xyz_id, step_id, abc_id, spin_id
331integer  :: cell_id, n1_id, n2_id, n3_id, gridfunc_id
332integer  :: origin_id, is_periodic_id
333
334integer   ::    nspin  ! Number of spins
335integer   ::    n(3), int_from_bool(3), i
336integer   ::   ispin, iostat, ix, iy, iz, iret
337
338logical   ::   has_origin, has_is_periodic
339real(dp)  ::    cell(3,3)
340!-----------------------------------------------------
341
342call clean_gridfunc(gf)
343
344call check( nf90_open(filename,NF90_NOWRITE,ncid))
345       call check( nf90_inq_dimid(ncid,'spin',spin_id) )
346       call check( nf90_inquire_dimension(ncid, dimid=spin_id, len=nspin) )
347
348       call check( nf90_inq_dimid(ncid,'n1',n1_id) )
349       call check( nf90_inquire_dimension(ncid, dimid=n1_id, len=n(1)) )
350       call check( nf90_inq_dimid(ncid,'n2',n2_id) )
351       call check( nf90_inquire_dimension(ncid, dimid=n2_id, len=n(2)) )
352       call check( nf90_inq_dimid(ncid,'n3',n3_id) )
353       call check( nf90_inquire_dimension(ncid, dimid=n3_id, len=n(3)) )
354
355       call check( nf90_inq_varid(ncid, "cell", cell_id) )
356
357       iret =  nf90_inq_varid(ncid, "origin", origin_id)
358       has_origin = (iret == nf90_noerr)
359       iret =  nf90_inq_varid(ncid, "is_periodic", is_periodic_id)
360       has_is_periodic = (iret == nf90_noerr)
361
362       call check( nf90_inq_varid(ncid, "gridfunc", gridfunc_id) )
363
364       iret = nf90_get_var(ncid, cell_id, cell, start = (/1, 1 /), &
365            count = (/3, 3/) )
366
367       if (has_origin) then
368          iret = nf90_get_var(ncid, origin_id, gf%origin, start = [1], count = [3] )
369          call check(iret)
370       else
371          gf%origin(:) = 0.0_dp
372       endif
373
374       if (has_is_periodic) then
375          iret = nf90_get_var(ncid, is_periodic_id, int_from_bool, start = [1], count = [3] )
376          call check(iret)
377          do i = 1, 3
378             gf%is_periodic(i) = (int_from_bool(i) == 1)
379          enddo
380       else
381          gf%is_periodic(:) = .true.
382       endif
383
384   gf%n(:) = n(:)
385   gf%cell = cell
386
387   allocate(gf%val(n(1),n(2),n(3),nspin))
388
389   iret = nf90_get_var(ncid, gridfunc_id, gf%val, start = (/1, 1, 1, 1 /), &
390        count = (/n(1), n(2), n(3), nspin/) )
391   call check(iret)
392
393   call check( nf90_close(ncid) )
394
395end subroutine read_gridfunc_netcdf
396
397subroutine check(code)
398use netcdf, only: nf90_noerr, nf90_strerror
399integer, intent(in) :: code
400if (code /= nf90_noerr) then
401  print *, "netCDF error: " // NF90_STRERROR(code)
402  STOP
403endif
404end subroutine check
405
406#endif /* CDF */
407
408 subroutine get_lun(lun)
409   integer, intent(out) :: lun
410
411   logical :: busy
412   do lun = 1, 99
413      inquire(unit=lun,opened=busy)
414      if (.not. busy) RETURN
415   enddo
416   call die("Cannot get free unit")
417 end subroutine get_lun
418
419 subroutine die(str)
420   character(len=*), intent(in) :: str
421   write(0,"(a)") str
422   stop
423 end subroutine die
424
425 !> Physically, it should not matter whether the 'c' axis points along
426 !> z, since rotations are irrelevant. However, some programs are
427 !> hard-wired for the z-direction. Hence the name. See the
428 !> 'monoclinic' function for the rotationally invariant
429 !> version.
430
431 function monoclinic_z(cell)
432   real(dp), intent(in) :: cell(3,3)
433   logical monoclinic_z
434
435   real(dp), parameter :: tol = 1.0e-8_dp
436
437   ! Check that the 'c' vector is along z, and 'a' and 'b' in the XY plane
438
439   monoclinic_z =  (abs(CELL(3,1)) < tol   &
440               .and. abs(CELL(3,2)) < tol  &
441               .and. abs(CELL(1,3)) < tol  &
442               .and. abs(CELL(2,3)) < tol )
443
444   end function monoclinic_z
445
446   !> Checks that the lattice vector of index 'dim' is orthogonal to
447   !> the other two.
448
449   function monoclinic(cell,dim)
450     real(dp), intent(in) :: cell(3,3)
451     integer, intent(in)  :: dim
452     logical monoclinic
453
454     real(dp), parameter :: tol = 1.0e-8_dp
455
456   ! Check that the 'dim' vector is orthogonal to the other two
457     monoclinic = .true.
458     if (abs(dot_product(cell(:,dim),cell(:,next(dim,1)))) > tol) then
459        monoclinic = .false.
460     endif
461     if (abs(dot_product(cell(:,dim),cell(:,next(dim,2)))) > tol) then
462        monoclinic = .false.
463     endif
464
465   CONTAINS
466     !> generates cyclic neighbors of dimension 'i' with shift 'm'
467     !> For example:
468     !>   x --> y:   next(1,1) = 2
469     !>   y --> x:   next(2,2) = 1
470     function next(i,m) result(ind)
471       integer, intent(in) :: i, m
472       integer :: ind
473
474       ind = i + m
475       if (ind > 3) then
476          ind = ind - 3
477       endif
478     end function next
479   end function monoclinic
480
481 end module m_gridfunc
482
483