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