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 fully created by: 9! Nick Papior Andersen, 2013, nickpapior@gmail.com 10! 11 12! When dealing with high bias one often realizes that 13! the integral becomes rather murky. 14! One can then do a monitor run and check which energy ranges 15! account for the most noise in the signal. 16! This allows the user to customize the energy contour 17! segments to optimize the convergence 18! It will probably only be really useful for large bias, but should be rather 19! easy to use. 20 21module m_monitor 22 23 use precision, only : dp 24 25 implicit none 26 27contains 28 29 subroutine read_monitor(bName,block_dist,sparse,N,orb_list) 30 31 use class_OrbitalDistribution 32 use class_Sparsity 33 use parallel, only : IONode, Node, Nodes 34 use fdf 35 use alloc 36 use geom_helper, only : UCORB 37#ifdef MPI 38 use mpi_siesta 39#endif 40 41 ! the block for reading in the monitor list 42 character(len=*), intent(in) :: bName 43 ! The distribution for the sparsity-pattern 44 type(OrbitalDistribution), intent(inout) :: block_dist 45 ! the sparse pattern we wish to monitor 46 type(Sparsity), intent(inout) :: sparse 47 integer, intent(inout) :: N 48 integer, pointer :: orb_list(:,:) 49 50 integer, pointer :: tmp(:,:) => null() 51#ifdef MPI 52 integer, pointer :: tmpMPI(:,:) => null() 53#endif 54 55 type(block_fdf) :: bfdf 56 type(parsed_line), pointer :: pline => null() 57 58 integer, pointer :: l_ncol(:) => null() 59 integer, pointer :: l_ptr(:) => null() 60 integer, pointer :: l_col(:) => null() 61 62 integer :: iM 63 integer :: no_u, no_local 64 integer :: io, j, ptr, ic, jc 65 66 character(len=50) :: c 67#ifdef MPI 68 integer :: MPIerror 69#endif 70 71 ! This will read in the monitored orbitals 72 if ( N /= 0 ) return 73 74 ! if the block does not exist simply return 75 if ( .not. fdf_block(trim(bName),bfdf) ) return 76 77 ! The block we read is formatted like this: 78 79 ! %block <name> 80 ! # Onsite orbital 81 ! 200 82 ! # Connection orbital (if it exists) 83 ! # if it doesn't, nothing will happen 84 ! 200 210 85 ! %endblock <name> 86 87 ! First we read the number of entries in the block 88 89 N = 0 90 do while ( fdf_bline(bfdf,pline) ) 91 if ( fdf_bnintegers(pline) > 0 ) then 92 N = N + 1 93 end if 94 end do 95 96 ! if it is empty, return 97 if ( N == 0 ) return 98 99 ! allocate space 100 call re_alloc(tmp,1,3,1,N, routine='monitor') 101 tmp(3,:) = -1 102 103 ! read in the actual numbers 104 ! first rewind... 105 call fdf_brewind(bfdf) 106 107 N = 0 108 do while ( fdf_bline(bfdf,pline) ) 109 if ( fdf_bnintegers(pline) > 0 ) then 110 N = N + 1 111 112 tmp(1,N) = fdf_bintegers(pline,1) 113 if ( fdf_bnintegers(pline) > 1 ) then 114 tmp(2,N) = fdf_bintegers(pline,2) 115 else 116 ! onsite... 117 tmp(2,N) = tmp(1,N) 118 end if 119 end if 120 end do 121 122 ! done reading 123 124 ! find the indices in the sparse pattern 125 call attach(sparse,nrows=no_local,nrows_g=no_u, & 126 n_col=l_ncol,list_ptr=l_ptr,list_col=l_col) 127 128 do io = 1 , no_local 129 130 ! Get the global position 131 ic = index_local_to_global(block_dist,io,Node) 132 133 ! quick loop if not requested 134 if ( all(tmp(1,:) /= ic) ) cycle 135 136 ! Loop number of entries in the row... 137 do j = 1 , l_ncol(io) 138 139 ! The column 140 jc = ucorb(l_col(l_ptr(io) + j),no_u) 141 142 do iM = 1 , N 143 if ( tmp(1,iM) == ic .and. & 144 tmp(2,iM) == jc ) then 145 tmp(3,iM) = l_ptr(io) + j 146 end if 147 end do 148 149 end do 150 end do 151 152 ! remove any dead references 153 iM = count(tmp(3,:) > 0 ) 154 if ( iM == 0 ) then 155 ! everything the user has requested is non-existing 156 call de_alloc(tmp,routine='monitor') 157 N = 0 158 return 159 end if 160 if ( associated(orb_list) ) & 161 call de_alloc(orb_list,routine='monitor') 162 nullify(orb_list) 163 call re_alloc(orb_list,1,iM,1,4, & 164 routine='monitor') 165 iM = 0 166 do io = 1 , N 167 if ( tmp(3,io) == -1 ) cycle 168 iM = iM + 1 169 ! copy over orb_list 170 orb_list(iM,1) = tmp(1,io) 171 orb_list(iM,2) = tmp(2,io) 172 orb_list(iM,3) = 1 173 orb_list(iM,4) = tmp(3,io) 174 end do 175 call de_alloc(tmp,routine='monitor') 176 nullify(tmp) 177 178 return 179 180 ! TODO finalize this to be able to deal with a distributed 181 ! orbital 182#ifdef MPI 183 ! retrieve the number of orbitals on each node 184 !call re_alloc(tmpMPI,0,Nodes-1) 185 !tmpMPI = 0 186 !tmpMPI(Node) = N 187 ! The IONode should have the globalized array 188 !call MPI_Reduce(N,iM,1,MPI_Integer,MPI_Sum,0, MPI_Comm_World, MPIerror) 189 !call MPI_Reduce(MPI_InPlace,tmpMPI,Nodes,MPI_Integer, & 190 ! MPI_Sum,0, MPI_Comm_World, MPIerror) 191 192 !if ( IONode ) then 193 ! call re_alloc(tmp,1,2,1,iM,routine='monitor') 194 !end if 195 196#endif 197 198 end subroutine read_monitor 199 200 201 function fname_monitor(iB,iC,slabel,suffix,basename) result(fname) 202 integer, intent(in) :: iB, iC 203 character(len=*), intent(in), optional :: slabel, suffix, basename 204 character(len=200) :: fname 205 ! Initialize the files 206 if ( present(basename) ) then 207 if ( iB == iC ) then 208 write(fname,'(2a,i0)') trim(basename), & 209 '_',iB 210 else 211 write(fname,'(a,2(a,i0))') trim(basename), & 212 '_',iB,'_',iC 213 end if 214 else 215 if ( .not. & 216 ( present(slabel) .and. present(suffix) ) ) & 217 call die('Error in filename input') 218 if ( iB == iC ) then 219 write(fname,'(4a,i0)') trim(slabel),'.', & 220 trim(suffix),'_',iB 221 else 222 write(fname,'(3a,2(a,i0))') trim(slabel),'.', & 223 trim(suffix),'_',iB,'_',iC 224 end if 225 end if 226 end function fname_monitor 227 228end module m_monitor 229 230 231 232 233 234 235 236 237 238 239