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