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      module m_broyden_mixing
9
10      use precision, only: dp
11
12      use m_broyddj, only: broyden_t, broyden_init, broyden_is_setup
13      use m_broyddj, only: broyden_reset, broyden_step
14      use fdf
15      use alloc,     only: re_alloc, de_alloc
16
17      use parallel, only: ionode
18      use m_mpi_utils, only: globalize_sum, globalize_max
19
20      public :: broyden_mixing
21      private
22
23      CONTAINS
24
25      subroutine broyden_mixing(iscf,mix_scf1,nbasis,maxnd,numd,
26     .                   listdptr,nspin,alpha,nkick,alpha_kick,
27     $                   dmnew,dmold,dmax)
28
29C ************************** INPUT **************************************
30
31C integer iscf               : Current SCF iteration
32C integer nbasis             : Number of atomic orbitals stored locally
33C integer maxnd              : First dimension of D.M., and
34C                              maximum number of nonzero elements of D.M.
35C integer numd(:)            : Control vector of D.M.
36C                              (number of nonzero elements of each row)
37C integer nspin              : Spin polarization (1=unpolarized, 2=polarized,
38C                               4=Non-collinear)
39C real*8 alpha               : Mixing parameter (for linear mixing)
40C integer nkick              : A kick is given every nkick iterations
41C real*8 alpha_kick          : Mixing parameter for kicks
42C logical mix_scf1           : Mix on first iteration?
43C ********************* INPUT AND OUTPUT*********************************
44C real*8 dmnew(maxnd)        : Density Matrix
45C                              Input: d.m. output in current SCF step
46C                              Output: d.m. input for next SCF iteration
47C real*8 dmold(maxnd)        : Density matrix
48C                              Input: d.m. input in current SCF step
49C                              Output: d.m. input for next SCF iteration
50C ************************** OUTPUT *************************************
51C real*8 dmax                : Maximum change of a DM element between
52C                              input and output
53
54      implicit none
55
56      integer, intent(in) :: iscf,maxnd,nbasis,nspin,nkick
57      logical, intent(in) :: mix_scf1
58      integer, intent(in) ::  numd(:), listdptr(:)
59      real(dp), intent(in) :: alpha, alpha_kick
60      real(dp), intent(inout) :: dmnew(maxnd,nspin),
61     $                           dmold(maxnd,nspin)
62      real(dp), intent(out) ::  dmax
63
64
65      logical, save           :: initialization_done = .false.
66
67
68      integer ::   i0,i,is,j, numel,k, global_numel
69
70      real(dp), dimension(:), pointer  :: rold, rnew, rdiff
71
72      type(broyden_t), save ::  br
73
74      real(dp), save :: jinv0
75      integer, save  :: maxit
76      logical, save  :: cycle_on_maxit, variable_weight
77      logical, save  :: do_not_mix, broyden_debug
78
79      real(dp) :: global_dnorm, global_dmax,  dnorm, diff, weight
80
81      numel = nspin * sum(numd(1:nbasis))
82      call Globalize_sum(numel,global_numel)
83
84      if (.not. initialization_done) then
85
86        if (ionode) then
87          print *, "Broyden: No of relevant DM elements: ", global_numel
88        endif
89        maxit = fdf_integer("DM.Number.Broyden",5)
90        cycle_on_maxit =
91     $          fdf_boolean("DM.Broyden.Cycle.On.Maxit",.true.)
92        variable_weight =
93     $          fdf_boolean("DM.Broyden.Variable.Weight",.true.)
94        broyden_debug =
95     $          fdf_boolean("DM.Broyden.Debug",.false.)
96
97        jinv0 = fdf_double("DM.Broyden.Initial.Mixing",alpha)
98        if (ionode) then
99          print *, "maxit for broyden: ", maxit
100          print *, "cycle on maxit: ", cycle_on_maxit
101          print *, "variable weight: ", variable_weight
102          print *, "initial alpha: ", jinv0
103        endif
104
105        call broyden_init(br,broyden_debug)
106
107        initialization_done = .true.
108
109      endif
110
111      do_not_mix = (iscf == 1 .and. .not. mix_scf1)
112      if (kick_is_due(iscf,nkick) .or. do_not_mix) then
113
114          ! Kick without saving any history for later
115         if (broyden_debug .and. ionode) then
116            if (do_not_mix) then
117               print *, "No mix in first iteration"
118            else
119               print *, "Kick"
120            endif
121         endif
122
123        ! Linear mixing with alpha_kick (or no mixing)
124
125         dmax = 0.0_dp
126         dnorm = 0.0_dp
127         do is = 1,nspin
128            do i = 1,nbasis
129               do j = 1,numd(i)
130                  k = listdptr(i) + j
131                  diff = dmnew(k,is) - dmold(k,is)
132                  dmax = max(dmax, abs(diff))
133                  dnorm = dnorm + diff**2
134                  if (.not. do_not_mix) then
135                     dmnew(k,is) = dmold(k,is) + alpha_kick*diff
136                  endif
137                  dmold(k,is) = dmnew(k,is)
138               enddo
139            enddo
140         enddo
141         call Globalize_sum(dnorm,global_dnorm)
142         call Globalize_max(dmax,global_dmax)
143!
144         global_dnorm = sqrt(global_dnorm)
145         dmax = global_dmax
146         if (broyden_debug .and. ionode)
147     $                   print *, "Dnorm = ", global_dnorm
148
149         if (broyden_is_setup(br)) then
150            if (broyden_debug .and. ionode)
151     $           print *, "Resetting history after kick or 1st."
152            call broyden_reset(br,numel,maxit,cycle_on_maxit,
153     $           jinv0,0.01_dp)
154         endif
155
156         RETURN
157
158      endif       ! Linear mixing
159!
160!     Broyden section
161!
162      nullify( rold )
163      call re_alloc( rold, 1, numel, name='rold',
164     &               routine='broyden_mixing' )
165      nullify( rnew )
166      call re_alloc( rnew, 1, numel, name='rnew',
167     &               routine='broyden_mixing' )
168      nullify( rdiff )
169      call re_alloc( rdiff, 1, numel, name='rdiff',
170     &               routine='broyden_mixing' )
171!
172!          Copy input to auxiliary arrays
173!          (memory will be saved by inlining the whole thing
174!           later on)
175!
176           i0 = 0
177           dmax = 0.0_dp
178           dnorm = 0.0_dp
179           do is = 1,nspin
180             do i = 1,nbasis
181               do j = 1,numd(i)
182                 i0 = i0 + 1
183                 k = listdptr(i) + j
184                 rold(i0) = dmold(k,is)
185                 rdiff(i0) = dmnew(k,is) -  dmold(k,is)
186                 dmax = max(dmax, abs(rdiff(i0)))
187                 dnorm = dnorm + rdiff(i0)**2
188               enddo
189             enddo
190           enddo
191           call Globalize_sum(dnorm,global_dnorm)
192           call Globalize_max(dmax,global_dmax)
193!
194           global_dnorm = sqrt(global_dnorm)
195           dmax = global_dmax
196           if (broyden_debug .and. ionode)
197     $                  print *, "Dnorm = ", global_dnorm
198!
199
200           if (iscf == 1) then
201              if (broyden_debug .and. ionode) then
202              ! Forget about history upon new scf cycle
203                 print *, "Resetting history for new SCF cycle."
204                 print *, "Broyden: No of relevant DM elements: ",
205     $                global_numel
206              endif
207              call broyden_reset(br,numel,maxit,cycle_on_maxit,
208     $             jinv0,0.01_dp)
209
210           endif
211
212           if (.not. broyden_is_setup(br)) then
213              if (broyden_debug .and. ionode) then
214                 print *, "Broyden: No of relevant DM elements: ",
215     $                global_numel
216              endif
217              call broyden_reset(br,numel,maxit,cycle_on_maxit,
218     $             jinv0,0.01_dp)
219
220           endif
221
222!          Broyden step
223!
224           if (variable_weight) then
225!
226!            Heuristic weight 1 < w < 300
227!
228              weight = exp(1.0_dp/(global_dmax+0.20))
229              if (broyden_debug .and. ionode)
230     $                     print *, "weight: ", weight
231           else
232              weight = 1.0_dp
233           endif
234
235           call broyden_step(br,rold,rdiff,w=weight,newx=rnew)
236
237!
238!         Copy back the results
239!
240           i0 = 0
241           do is = 1,nspin
242             do i = 1,nbasis
243               do j = 1,numd(i)
244                 i0 = i0 + 1
245                 k = listdptr(i)+j
246                 dmnew(k,is) = rnew(i0)
247                 dmold(k,is) = dmnew(k,is)
248               enddo
249             enddo
250           enddo
251
252           call de_alloc( rold, name='rold' )
253           call de_alloc( rnew, name='rnew' )
254           call de_alloc( rdiff, name='rdiff' )
255
256      end subroutine broyden_mixing
257
258!---------------------------------------------------------
259      function kick_is_due(n,step) result(due)
260      logical :: due
261      integer, intent(in) :: n, step
262
263      due = .false.
264      if (step .le. 0) RETURN
265      due = (modulo(n,step) == 0)
266
267      end function kick_is_due
268!---------------------------------------------------------
269
270      end module m_broyden_mixing
271
272