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