broyden2.F90 14.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 106 107 108 109
    dNVec = 0.0
    dFVec = 0.0
    vVec  = 0.0
    dNMet = 0.0
    dFMet = 0.0

    mit = 0
    l_exist = initBroydenHistory(input,hybrid,nmap) ! returns true if there already exists a Broyden history
    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 131
    END IF

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

    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>  
144 145
!       CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!                   mmap,nmaph,mapmt,mapvac2,dNVec,dNMet,l_pot)
146 147

       ! Apply metric w to dFVec and store in dFMet:  w |dFVec>  
148 149
!       CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!                   mmap,nmaph,mapmt,mapvac2,dFVec,dFMet,l_pot)
150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 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

       ! 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-

       iread = MIN(mit-1,input%maxiter+1)
       historyLength = iread

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

       DO it = 2, iread
          CALL readDeltaNVec(input,hybrid,nmap,it-mit,mit,deltaN_i)
          CALL readDeltaFVec(input,hybrid,nmap,it-mit,mit,deltaF_i)

          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)

       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)

       ! 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
          CALL readDeltaNVec(input,hybrid,nmap,it-mit,mit,deltaN_i)
          CALL readDeltaFVec(input,hybrid,nmap,it-mit,mit,deltaF_i)

          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>
349 350
!       CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
!                   mmap,nmaph,mapmt,mapvac2,fm,FMet,l_pot)
351

352
       vFMetProd = CPP_BLAS_sdot(nmap,vVec,1,fmMet,1)
353 354 355 356 357 358

       ! 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

359
    DEALLOCATE (fmLast,smLast,fmLastMet,smLastMet,dNVec,dFVec,vVec)
360 361 362

  END SUBROUTINE broyden2
END MODULE m_broyden2