broyden2.F90 15.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
!--------------------------------------------------------------------------------
! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! This is a continuously restartable implementation of the original
!!! broyden subroutine. It constructs the required u and v vectors for
!!! the relevant iterations on the fly in terms of linear combinations
!!! of the deltaN_i and deltaF_i vectors. It is slower than the original
!!! broyden subroutine so if the continuous restartability is not needed
!!! the old routine should be used.
!!!
!!!                                   GM'2018
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

MODULE m_broyden2
  USE m_juDFT
  !################################################################
  !     IMIX = 3 : BROYDEN'S FIRST METHOD
  !     IMIX = 5 : BROYDEN'S SECOND METHOD
  !     IMIX = 7 : GENERALIZED ANDERSEN METHOD
  !     sm   : input charge density of iteration m
  !            afterwards update rho(m+1)
  !     fm   : output minus input charge density of iteration m
  !################################################################
CONTAINS
  SUBROUTINE broyden2(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
32
                      hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fm,sm,fmMet,smMet,lpot)
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

#include"cpp_double.h"

    USE m_metric
    USE m_types
    USE m_broyd_io

    IMPLICIT NONE

    TYPE(t_oneD),INTENT(IN)    :: oneD
    TYPE(t_input),INTENT(IN)   :: input
    TYPE(t_vacuum),INTENT(IN)  :: vacuum
    TYPE(t_noco),INTENT(IN)    :: noco
    TYPE(t_sym),INTENT(IN)     :: sym
    TYPE(t_stars),INTENT(IN)   :: stars
    TYPE(t_cell),INTENT(IN)    :: cell
    TYPE(t_sphhar),INTENT(IN)  :: sphhar
    TYPE(t_atoms),INTENT(IN)   :: atoms
    TYPE(t_hybrid),INTENT(IN)  :: hybrid

    ! Scalar Arguments
    INTEGER, INTENT (IN)        :: mmap,nmap
    INTEGER, INTENT (IN)        :: mapmt,mapvac2
    LOGICAL,OPTIONAL,INTENT(IN) :: lpot

    ! Array Arguments
    REAL,    INTENT (IN)    :: fm(nmap) 
    REAL,    INTENT (INOUT) :: sm(nmap)
61 62
    REAL,    INTENT (IN)    :: fmMet(nmap) 
    REAL,    INTENT (INOUT) :: smMet(nmap)
63 64 65 66 67 68 69

    ! Local Scalars
    INTEGER         :: i,j,it,k,nit,iread,nmaph, mit, historyLength
    REAL            :: vFMetProd,alphan,coeff,vNorm
    LOGICAL         :: l_pot, l_exist

    ! Local Arrays
70 71
    REAL, ALLOCATABLE :: fmLast(:),smLast(:),fmLastMet(:),smLastMet(:)
    REAL, ALLOCATABLE :: uVec(:),vVec(:)
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
    REAL, ALLOCATABLE :: dNVec(:),dFVec(:),dNMet(:), dFMet(:), FMet(:)
    REAL, ALLOCATABLE :: deltaN_i(:), deltaF_i(:)
    REAL, ALLOCATABLE :: dNdNMat(:,:), dNdFMat(:,:), dFdNMat(:,:), dFdFMat(:,:)
    REAL, ALLOCATABLE :: dNdNLast(:), dNdFLast(:), dFdNLast(:), dFdFLast(:)
    REAL, ALLOCATABLE :: uDNTableau(:,:), vDNTableau(:,:), wDNTableau(:,:)
    REAL, ALLOCATABLE :: uDFTableau(:,:), vDFTableau(:,:), wDFTableau(:,:)

    ! External Functions
    REAL CPP_BLAS_sdot
    EXTERNAL CPP_BLAS_sdot

    ! External Subroutines
    EXTERNAL CPP_BLAS_saxpy,CPP_BLAS_sscal

    l_pot = .FALSE.
    IF (PRESENT(lpot)) l_pot = lpot

89 90
    ALLOCATE (fmLast(mmap),smLast(mmap),fmLastMet(mmap),smLastMet(mmap))
    ALLOCATE (dNVec(mmap),dFVec(mmap),uVec(mmap),vVec(mmap))
91 92 93 94
    ALLOCATE (deltaN_i(mmap),deltaF_i(mmap))
    ALLOCATE (dNdNLast(input%maxiter),dFdFLast(input%maxiter))
    ALLOCATE (dNdFLast(input%maxiter),dFdNLast(input%maxiter))

95
    ALLOCATE (dNMet(mmap), dFMet(mmap), FMet(mmap))
96

97 98
    fmLast = 0.0
    smLast = 0.0
99 100 101 102 103 104 105
    dNVec = 0.0
    dFVec = 0.0
    vVec  = 0.0
    dNMet = 0.0
    dFMet = 0.0

    mit = 0
106
    l_exist = initBroydenHistory2(input,hybrid,nmap) ! returns true if there already exists a Broyden history
107 108 109
    IF(.NOT.l_exist) mit = 1

    IF (mit.NE.1) THEN
110 111
       ! load input charge density (smLast) and difference of 
       ! in and out charge densities (fmLast) from previous iteration (m-1)
112

113
       CALL readLastIterInAndDiffDen(hybrid,nmap,mit,alphan,smLast(:nmap),fmLast(:nmap),smLastMet(:nmap),fmLastMet(:nmap))
114 115 116 117 118 119
       IF (ABS(input%alpha-alphan) > 0.0001) THEN
          WRITE (6,*) 'mixing parameter has been changed; reset'
          WRITE (6,*) 'broyden algorithm or set alpha to',alphan
          CALL juDFT_error("mixing parameter (input) changed", calledby ="broyden")
       END IF

120 121 122 123 124 125 126
       ! generate F_m   - F_(m-1)  ... smLast
       !      and rho_m - rho_(m-1) .. fmLast
       dNVec(1:nmap) = sm(1:nmap) - smLast(1:nmap)
       dFVec(1:nmap) = fm(1:nmap) - fmLast(1:nmap)

       dNMet(1:nmap) = smMet(1:nmap) - smLastMet(1:nmap)
       dFMet(1:nmap) = fmMet(1:nmap) - fmLastMet(1:nmap)
127 128 129 130
    END IF

    ! save F_m and rho_m for next iteration
    nit = mit +1
131
    CALL writeLastIterInAndDiffDen(hybrid,nmap,nit,input%alpha,sm,fm,smMet,fmMet)
132 133 134 135 136 137 138 139 140 141 142

    IF (mit.EQ.1) THEN 
       !     update for rho for mit=1 is straight mixing
       !     sm = sm + alpha*fm
       CALL CPP_BLAS_saxpy(nmap,input%alpha,fm,1,sm,1)
    ELSE

       CALL writeDeltaNVec(input,hybrid,nmap,mit,dNVec)
       CALL writeDeltaFVec(input,hybrid,nmap,mit,dFVec)

       ! Apply metric w to dNVec and store in dNMet:  w |dNVec>  
143 144
!       CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!                   mmap,nmaph,mapmt,mapvac2,dNVec,dNMet,l_pot)
145 146

       ! Apply metric w to dFVec and store in dFMet:  w |dFVec>  
147 148
!       CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!                   mmap,nmaph,mapmt,mapvac2,dFVec,dFMet,l_pot)
149 150 151 152

       ! Extend overlap matrices <delta n(i) | delta n(j)>, <delta F(i) | delta F(j)>,
       !                         <delta n(i) | delta F(j)>, <delta F(i) | delta n(j)> -start-

153
       iread = MIN(mit-1,input%maxiter)
154 155 156 157 158 159 160
       historyLength = iread

       dNdNLast = 0.0
       dFdFLast = 0.0
       dNdFLast = 0.0
       dFdNLast = 0.0

161 162 163 164 165
!       WRITE(1400,*) '========================================'
!       WRITE(1400,*) '========================================'
!       WRITE(1400,*) '========================================'
!       WRITE(1400,*) 'mit: ', mit
!       WRITE(1400,*) 'iread, historyLength: ', iread, historyLength
166
       DO it = 2, iread
167 168 169 170 171 172
          CALL readDeltaNVec(input,hybrid,nmap,it-iread-1,mit,deltaN_i)
          CALL readDeltaFVec(input,hybrid,nmap,it-iread-1,mit,deltaF_i)

!          WRITE(1400,'(4i7)') it,mit,it-iread-1,mit+(it-iread-1)
!          WRITE(1400,'(a,5f15.8)') 'deltaN_i: ', deltaN_i(1), deltaN_i(2), deltaN_i(3), deltaN_i(4), deltaN_i(5)
!          WRITE(1400,'(a,5f15.8)') 'deltaF_i: ', deltaF_i(1), deltaF_i(2), deltaF_i(3), deltaF_i(4), deltaF_i(5)
173 174 175 176 177 178 179 180 181 182 183 184

          dNdNLast(it-1) = CPP_BLAS_sdot(nmap,deltaN_i,1,dNMet,1)
          dFdFLast(it-1) = CPP_BLAS_sdot(nmap,deltaF_i,1,dFMet,1)
          dNdFLast(it-1) = CPP_BLAS_sdot(nmap,deltaN_i,1,dFMet,1)
          dFdNLast(it-1) = CPP_BLAS_sdot(nmap,deltaF_i,1,dNMet,1)
       END DO

       dNdNLast(historyLength) = CPP_BLAS_sdot(nmap,dNVec,1,dNMet,1)
       dFdFLast(historyLength) = CPP_BLAS_sdot(nmap,dFVec,1,dFMet,1)
       dNdFLast(historyLength) = CPP_BLAS_sdot(nmap,dNVec,1,dFMet,1)
       dFdNLast(historyLength) = CPP_BLAS_sdot(nmap,dFVec,1,dNMet,1)

185 186 187 188 189
!       WRITE(1400,*) 'last overlaps:'
!       DO i = 1, historyLength
!          WRITE(1400,'(i7,4f20.13)') i,dNdNLast(i),dFdFLast(i),dNdFLast(i),dFdNLast(i)
!       END DO

190 191 192 193 194 195 196 197 198 199 200
       CALL writeBroydenOverlapExt(input,hybrid,mit,historyLength,&
                                   dNdNLast,dFdFLast,dNdFLast,dFdNLast)

       ALLOCATE (dNdNMat(historyLength,historyLength))
       ALLOCATE (dFdFMat(historyLength,historyLength))
       ALLOCATE (dNdFMat(historyLength,historyLength))
       ALLOCATE (dFdNMat(historyLength,historyLength))

       CALL readBroydenOverlaps(input,hybrid,mit,historyLength,&
                                dNdNMat,dFdFMat,dNdFMat,dFdNMat)

201 202 203 204 205 206 207 208
!       WRITE(1400,*) 'all overlaps'
!       DO i = 1, historyLength
!          DO j = 1, historyLength
!             WRITE(1400,'(2i7,4f20.13)') i,j,dNdNMat(j,i),dFdFMat(j,i),dNdFMat(j,i),dFdNMat(j,i)
!          END DO
!       END DO
!       WRITE(1400,*) '-----------------------------'

209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349
       ! Extend overlap matrices <delta n(i) | delta n(j)>, <delta F(i) | delta F(j)>,
       !                         <delta n(i) | delta F(j)>, <delta F(i) | delta n(j)> -end-

       ! Construct u_i, v_i tableaus (prefactors for delta n(i) and delta F(i) vectors) -start-

       ALLOCATE (uDNTableau(historyLength,historyLength)) !first index: iteration of DN, second index: iteration of u
       ALLOCATE (vDNTableau(historyLength,historyLength)) !first index: iteration of DN, second index: iteration of v
       ALLOCATE (uDFTableau(historyLength,historyLength)) !first index: iteration of DF, second index: iteration of u
       ALLOCATE (vDFTableau(historyLength,historyLength)) !first index: iteration of DF, second index: iteration of v

       uDNTableau = 0.0
       vDNTableau = 0.0
       uDFTableau = 0.0
       vDFTableau = 0.0

       IF (input%imix.EQ.3) THEN ! Broyden's first method

          DO i = 1, historyLength
             uDNTableau(i,i) = 1.0
             uDFTableau(i,i) = input%alpha
             vDNTableau(i,i) = input%alpha

             DO j = 1, i-1
                ! Calculations for u_i
                coeff = 0.0
                DO k = 1, historyLength
                   coeff = coeff + vDNTableau(k,j)*dNdFMat(k,i)
                   coeff = coeff + vDFTableau(k,j)*dFdFMat(k,i)
                END DO
                DO k = 1, historyLength
                   uDNTableau(k,i) = uDNTableau(k,i) - coeff*uDNTableau(k,j)
                   uDFTableau(k,i) = uDFTableau(k,i) - coeff*uDFTableau(k,j)
                END DO

                ! Calculations for v_i (numerator)
                coeff = 0.0
                DO k = 1, historyLength
                   coeff = coeff + uDNTableau(k,j)*dNdNMat(k,i) ! This is in agreement with the old implementation.
                   coeff = coeff + uDFTableau(k,j)*dFdNMat(k,i)
!                   coeff = coeff + vDNTableau(k,j)*dNdNMat(k,i) ! This is in agreement with RP's thesis.
!                   coeff = coeff + vDFTableau(k,j)*dFdNMat(k,i)
                END DO
                DO k = 1, historyLength
                   vDNTableau(k,i) = vDNTableau(k,i) - coeff*vDNTableau(k,j) ! This is in agreement with the old implementation.
                   vDFTableau(k,i) = vDFTableau(k,i) - coeff*vDFTableau(k,j)
!                   vDNTableau(k,i) = vDNTableau(k,i) - coeff*uDNTableau(k,j) ! This is in agreement with RP's thesis.
!                   vDFTableau(k,i) = vDFTableau(k,i) - coeff*uDFTableau(k,j)
                END DO
             END DO

             ! Calculations for v_i (denominator)
             vNorm = dNdNMat(i,i)
             DO k = 1, historyLength
                vNorm = vNorm - uDNTableau(k,i)*dNdNMat(k,i)
                vNorm = vNorm - uDFTableau(k,i)*dFdNMat(k,i)
             END DO

             vNorm = -vNorm ! Note: This is in agreement with the old implementation.
                            !       The equations in RP's thesis don't have this sign.

             DO k = 1, historyLength
                vDNTableau(k,i) = vDNTableau(k,i) / vNorm
                vDFTableau(k,i) = vDFTableau(k,i) / vNorm
             END DO

          END DO
          
       ELSE IF (input%imix.EQ.5) THEN ! Broyden's second method

          DO i = 1, historyLength
             uDNTableau(i,i) = 1.0
             uDFTableau(i,i) = input%alpha
             vDFTableau(i,i) = 1.0 / dFdFMat(i,i)

             DO j = 1, i-1
                coeff = dFdFMat(j,i)*vDFTableau(j,j)
                DO k = 1, historyLength
                   uDNTableau(k,i) = uDNTableau(k,i) - coeff*uDNTableau(k,j)
                   uDFTableau(k,i) = uDFTableau(k,i) - coeff*uDFTableau(k,j)
                END DO
             END DO
          END DO

       ELSE IF (input%imix.EQ.7) THEN ! generalized anderson method

          ALLOCATE (wDNTableau(historyLength,historyLength)) !first index: iteration of DN, second index: iteration of w
          ALLOCATE (wDFTableau(historyLength,historyLength)) !first index: iteration of DF, second index: iteration of w

          wDNTableau = 0.0
          wDFTableau = 0.0

          DO i = 1, historyLength

             uDNTableau(i,i) = 1.0
             uDFTableau(i,i) = input%alpha
             wDFTableau(i,i) = 1.0

             DO j = 1, i-1
                ! Calculations for u_i and w_i
                coeff = 0.0
                DO k = 1, historyLength
                   coeff = coeff + vDNTableau(k,j)*dNdFMat(k,i)
                   coeff = coeff + vDFTableau(k,j)*dFdFMat(k,i)
                END DO
                DO k = 1, historyLength
                   ! Calculations for u_i
                   uDNTableau(k,i) = uDNTableau(k,i) - coeff*uDNTableau(k,j)
                   uDFTableau(k,i) = uDFTableau(k,i) - coeff*uDFTableau(k,j)

                   ! Calculations for w_i
                   wDNTableau(k,i) = wDNTableau(k,i) - coeff*wDNTableau(k,j)
                   wDFTableau(k,i) = wDFTableau(k,i) - coeff*wDFTableau(k,j)
                END DO

             END DO

             ! Calculations for v_i
             vNorm = 0.0
             DO k = 1, historyLength
                vNorm = vNorm + wDNTableau(k,i)*dNdFMat(k,i)
                vNorm = vNorm + wDFTableau(k,i)*dFdFMat(k,i)
             END DO

             DO k = 1, historyLength
                vDNTableau(k,i) = wDNTableau(k,i) / vNorm
                vDFTableau(k,i) = wDFTableau(k,i) / vNorm
             END DO

          END DO

          DEALLOCATE (wDNTableau, wDFTableau)

       END IF

       ! Construct u_i, v_i tableau (prefactors for delta n(i) and delta F(i) vectors) -end-

       ! Construct the final uVec and vVec -start-
       uVec = 0.0
       vVec = 0.0

       DO it = 2, iread
350 351
          CALL readDeltaNVec(input,hybrid,nmap,it-iread-1,mit,deltaN_i)
          CALL readDeltaFVec(input,hybrid,nmap,it-iread-1,mit,deltaF_i)
352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369

          DO k = 1, nmap
             uVec(k) = uVec(k) + uDNTableau(it-1,historyLength)*deltaN_i(k)
             uVec(k) = uVec(k) + uDFTableau(it-1,historyLength)*deltaF_i(k)
             vVec(k) = vVec(k) + vDNTableau(it-1,historyLength)*deltaN_i(k)
             vVec(k) = vVec(k) + vDFTableau(it-1,historyLength)*deltaF_i(k)
          END DO
       END DO

       DO k = 1, nmap
          uVec(k) = uVec(k) + uDNTableau(historyLength,historyLength)*dNVec(k)
          uVec(k) = uVec(k) + uDFTableau(historyLength,historyLength)*dFVec(k)
          vVec(k) = vVec(k) + vDNTableau(historyLength,historyLength)*dNVec(k)
          vVec(k) = vVec(k) + vDFTableau(historyLength,historyLength)*dFVec(k)
       END DO
       ! Construct the final uVec and vVec -end-

       ! Apply metric w to dFVec and store in FMet:  w |dFVec>
370 371
!       CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!                   mmap,nmaph,mapmt,mapvac2,fm,FMet,l_pot)
372

373
       vFMetProd = CPP_BLAS_sdot(nmap,vVec,1,fmMet,1)
374 375 376 377 378 379

       ! update rho(m+1)
       ! calculate sm(:) = (1.0-vFMetProd)*uVec(:) + sm
       CALL CPP_BLAS_saxpy(nmap,1.0-vFMetProd,uVec,1,sm,1)
    END IF

380
    DEALLOCATE (fmLast,smLast,fmLastMet,smLastMet,dNVec,dFVec,vVec)
381 382 383

  END SUBROUTINE broyden2
END MODULE m_broyden2