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! This code segment has been improved or fully created by:
9! Nick Papior Andersen, 2014, nickpapior@gmail.com
10!
11module m_mesh_node
12
13! Module for retaining information about the mesh.
14! It is used by the Hartree module and the bias module.
15!
16! Created and copyrighted by: Nick Papior Andersen, 2014
17! The use of this program is allowed for not-for-profit research only.
18! Copy or disemination of all or part of this package is not
19! permitted without prior and explicit authorization by the author.
20
21  use precision, only : dp
22
23  implicit none
24
25  public
26  save
27
28  ! The offset for the current node
29  real(dp) :: offset_r(3) = 0._dp
30  ! the voxel-vectors for each sub-mesh element
31  real(dp) :: dL(3,3) = 0._dp
32  ! the voxel length along each direction
33  real(dp) :: dMesh(3) = 0._dp
34
35  ! offsets for the local node
36  integer :: meshl(3) = 0
37  integer :: offset_i(3) = 0
38
39contains
40
41  subroutine init_mesh_node( ucell , meshG , meshLim , nsm)
42
43    use intrinsic_missing, only : VNORM
44    use parallel, only : Node, Nodes, ProcessorY, IONode
45    use units, only: Ang
46#ifdef MPI
47    use mpi_siesta
48#endif
49
50    ! The unit cell
51    real(dp), intent(in) :: ucell(3,3)
52    ! Number of mesh divisions of each lattice vector
53    integer, intent(in) :: meshG(3), meshLim(2,3)
54    ! Number of fine points per big point (see iogrid_netcdf)
55    integer, intent(in) :: nsm
56
57    ! Number of big division points
58    integer :: nm(3)
59    ! Processor specifics
60    integer :: ProcessorZ, blocY, blocZ, nremY, nremZ, idx(3)
61    ! dimension tracking of the divisions
62    integer :: iniX, iniY, iniZ, dimX, dimY, dimZ
63    ! Local node dimensionality of the grid
64    integer :: ldimX, ldimY, ldimZ
65    ! Loop stuff
66    integer :: node_Y, node_Z, cur_Node
67    integer :: i
68#ifdef MPI
69    real(dp) :: ll(3)
70    integer :: ix, iy, iz
71    integer :: MPIerror, status(MPI_STATUS_SIZE)
72#endif
73
74    ! We calculate the spacing in each direction
75    ! Notice that we now have:
76    ! dL(1,1) = dX for stepping in the x-direction
77    ! dL(1,2) = dX for stepping in the y-direction
78    ! dL(1,3) = dX for stepping in the z-direction
79    do i = 1 , 3
80       ldimX = max(meshG(i),1)
81       ! The dimension stepping in each direction.
82       dL(:,i) = ucell(:,i) / ldimX
83    end do
84    ! The voxel box-size in Cartesian coordinates
85    ! is calculated by adding all three vectors
86    dMesh = matmul(dL,(/1._dp,1._dp,1._dp/))
87
88    ! For nodes == 1 we have no offset
89    ! (also some of the arrays are not initialized, which
90    !  could lead to errors)
91    if ( Nodes == 1 ) then
92       meshl = meshG ! same as meshLim(2,:) * nsm
93       return
94    end if
95
96    ! Now we need to calculate the offset of the local node
97
98    ! Calculate the number of big-points
99    nm(1:3) = meshG(1:3) / nsm
100
101    ! In order to check that we do it correctly...
102    dimX  = nm(1) * nsm
103    ldimX = dimX
104
105    ! Calculate number of z-processor divisions
106    ProcessorZ = Nodes / ProcessorY
107
108    ! Calculate the block-division in the Y-direction
109    blocY = nm(2) / ProcessorY
110    ! Calculate the block-division in the Z-direction
111    blocZ = nm(3) / ProcessorZ
112    ! If there are any remaining mesh-points, they will be added one
113    ! to each of the processors, in sequence
114    nremY = nm(2) - blocY*ProcessorY
115    nremZ = nm(3) - blocZ*ProcessorZ
116
117    ! Initialize the loop construct variables
118    cur_Node = 0
119    ! Notice that we start from zero as we would like to calculate the
120    ! offset
121    iniX = 0
122    iniY = 0
123    do node_Y = 1, ProcessorY
124
125       ! Initialize the dimension size in the Y-direction
126       dimY = blocY
127       ! The if there were too many points to be perfectly divisable
128       ! we will add them.
129       if ( node_Y <= nremY ) dimY = dimY + 1
130       ! Extend into the mesh size
131       dimY = dimY * nsm
132
133       ! Initialize the Z-division (notice we index from zero
134       ! as we would like to calculate the offset)
135       iniZ = 0
136       do node_Z = 1, ProcessorZ
137          dimZ = blocZ
138          if ( node_Z <= nremZ ) dimZ = dimZ + 1
139          dimZ = dimZ * nsm ! For fine points
140
141          if ( Node == cur_Node ) then
142             ! Calculate the offset in the [YZ]-direction for the processor
143             ! We know that iniX == 0, so no need, but we have it for
144             ! consistency
145             offset_r(:) = iniX * dL(:,1) + iniY * dL(:,2) + iniZ * dL(:,3)
146             ldimY = dimY
147             ldimZ = dimZ
148          end if
149
150          iniZ = iniZ + dimZ
151          cur_Node = cur_Node + 1
152       end do
153       iniY = iniY + dimY
154    end do
155
156    ! Find quantities in mesh coordinates
157    meshl(1) = (meshLim(2,1) - meshLim(1,1)+1)*nsm
158    if ( ldimX /= meshl(1) ) &
159         call die('Incorrect number of A divisions found.')
160    meshl(2) = (meshLim(2,2) - meshLim(1,2)+1)*nsm
161    if ( ldimY /= meshl(2) ) &
162         call die('Incorrect number of B divisions found.')
163    meshl(3) = (meshLim(2,3) - meshLim(1,3)+1)*nsm
164    if ( ldimZ /= meshl(3) ) &
165         call die('Incorrect number of C divisions found.')
166
167    ! Calculate starting point for grid
168    offset_i(:) = (meshLim(1,:) - 1)*nsm
169
170#ifdef MESH_DEBUG
171    ! We will print out information about the boxes that each processor
172    ! has..
173    if ( IONode ) then
174       write(*,'(t2,a)') 'Printing offsets and corners in Ang'
175       write(*,'(t3,a,3(tr1,f10.5))') &
176               'dL-x: ',dL(:,1)/Ang
177       write(*,'(t3,a,3(tr1,f10.5))') &
178               'dL-y: ',dL(:,2)/Ang
179       write(*,'(t3,a,3(tr1,f10.5))') &
180               'dL-z: ',dL(:,3)/Ang
181       write(*,'(t3,a,3(tr1,f10.5))') &
182               'dMesh:',dMesh(:)/Ang
183    end if
184#ifdef MPI
185    do i = 0 , Nodes - 1
186       if ( Node == i ) then
187          if ( IONode ) then
188             ll = offset_r
189             ix = ldimX
190             iy = ldimY
191             iz = ldimZ
192             idx = offset_i
193          else
194             call MPI_Send(offset_r,3,MPI_Double_precision, &
195                  0,0,MPI_Comm_World,MPIerror)
196             call MPI_Send(ldimX,1,MPI_Integer, &
197                  0,1,MPI_Comm_World,MPIerror)
198             call MPI_Send(ldimY,1,MPI_Integer, &
199                  0,2,MPI_Comm_World,MPIerror)
200             call MPI_Send(ldimZ,1,MPI_Integer, &
201                  0,3,MPI_Comm_World,MPIerror)
202             call MPI_Send(offset_i,3,MPI_Integer, &
203                  0,4,MPI_Comm_World,MPIerror)
204          end if
205       else if ( IONode ) then
206          call MPI_Recv(ll,3,MPI_Double_precision, &
207               i,0,MPI_Comm_World,status,MPIerror)
208          call MPI_Recv(ix,1,MPI_Integer, &
209               i,1,MPI_Comm_World,status,MPIerror)
210          call MPI_Recv(iy,1,MPI_Integer, &
211               i,2,MPI_Comm_World,status,MPIerror)
212          call MPI_Recv(iz,1,MPI_Integer, &
213               i,3,MPI_Comm_World,status,MPIerror)
214          call MPI_Recv(idx,3,MPI_Integer, &
215               i,4,MPI_Comm_World,status,MPIerror)
216       end if
217       if ( IONode ) then
218          print *, i, ix,iy,iz
219          write(*,'(t3,a,3(tr1,i8))') &
220               'Indices:',idx(:)
221          write(*,'(t3,a,3(tr1,f12.5))') &
222               'Lower-left:',ll(:)/Ang
223          write(*,'(t3,a,3(tr1,f12.5))') &
224               'Upper-x   :',(ll(:)+ix * dL(:,1))/Ang
225          write(*,'(t3,a,3(tr1,f12.5))') &
226               'Upper-y   :',(ll(:)+iy * dL(:,2))/Ang
227          write(*,'(t3,a,3(tr1,f12.5))') &
228               'Upper-z   :',(ll(:)+iz * dL(:,3))/Ang
229       end if
230    end do
231#else
232    if ( IONode ) then
233       write(*,'(t3,a,3(tr1,i8))') &
234            'Indices:',idx(:)
235       write(*,'(t3,a,3(tr1,f12.5))') &
236            'Lower-left:',offset_r(:)/Ang
237       write(*,'(t3,a,3(tr1,f12.5))') &
238            'Upper-x   :',(offset_r(:)+ldimX * dL(:,1))/Ang
239       write(*,'(t3,a,3(tr1,f12.5))') &
240            'Upper-y   :',(offset_r(:)+ldimY * dL(:,2))/Ang
241       write(*,'(t3,a,3(tr1,f12.5))') &
242            'Upper-z   :',(offset_r(:)+ldimZ * dL(:,3))/Ang
243    end if
244#endif
245#endif
246
247  end subroutine init_mesh_node
248
249  elemental subroutine mesh_correct_idx(mesh,idx)
250    integer, intent(in) :: mesh
251    integer, intent(inout) :: idx
252    ! negative "supercell"
253    do while ( idx <= 0 )
254       idx = idx + mesh
255    end do
256    ! positive "supercell"
257    do while ( mesh < idx )
258       idx = idx - mesh
259    end do
260
261  end subroutine mesh_correct_idx
262
263end module m_mesh_node
264
265