Added alternativ implementation of the broyden subroutine for continuous restarting

This is not yet fully functional and not yet activated in the code. To activate it
the broyden call in "mix" has to be replaced by a broyden2 call.

eigen/hsmt_extra.F90

@@ -135,6 +135,7 @@ CONTAINS
        ENDIF
        DO nn = 1,atoms%neq(n)
           na = na + 1
+          i_u_save = i_u
           IF (atoms%lnonsph(n).LT.0) CYCLE ntype_loop
           IF ((atoms%invsat(na).EQ.0) .OR. (atoms%invsat(na).EQ.1)) THEN
              IF (atoms%invsat(na).EQ.0) invsfct = 1
@@ -270,12 +271,10 @@ CONTAINS
                 ENDIF
              END IF
 
-
              IF (atoms%n_u.GT.0) THEN
                 nkvecprevatuTemp = nkvecprevatu
                 iilouTemp = iilou
                 locoluTemp = locolu
-                i_u_save = i_u
                 DO WHILE (i_u.LE.atoms%n_u)
                    IF (atoms%lda_u(i_u)%atomType.GT.n) EXIT
                    nkvecprevatuTemp = nkvecprevatu

io/broyd_io.F90

+!--------------------------------------------------------------------------------
+! 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.
+!--------------------------------------------------------------------------------
+
 MODULE m_broyd_io
 
 USE m_types
@@ -192,6 +198,247 @@ SUBROUTINE writeVVec(input,hybrid,vecLen,currentIter,dfivi,vVec)
 
 END SUBROUTINE writeVVec
 
+
+
+
+
+
+
+
+
+SUBROUTINE readDeltaNVec(input,hybrid,vecLen,relIter,currentIter,deltaNVec)
+
+   TYPE(t_input),  INTENT(IN)  :: input
+   TYPE(t_hybrid), INTENT(IN)  :: hybrid
+   INTEGER,        INTENT(IN)  :: vecLen, relIter, currentIter
+   REAL,           INTENT(OUT) :: deltaNVec(vecLen)
+
+   INTEGER             :: mode, npos
+   INTEGER*8           :: recLen
+
+   CALL getIOMode(mode)
+
+   ! At the moment broyden IO is mode independent
+
+   recLen=(vecLen+1)*8
+
+   IF (hybrid%l_calhf) THEN
+      OPEN (59,file='hf_broyd_DN',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ELSE
+      OPEN (59,file='broyd_DN',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ENDIF
+
+   npos=currentIter+relIter-1
+   IF (currentIter.GT.input%maxiter+1) npos = MOD(currentIter-2,input%maxiter)+1
+
+   READ(59,rec=npos) deltaNVec(:vecLen)
+
+   CLOSE(59)
+
+END SUBROUTINE readDeltaNVec
+
+SUBROUTINE writeDeltaNVec(input,hybrid,vecLen,currentIter,deltaNVec)
+
+   TYPE(t_input),  INTENT(IN) :: input
+   TYPE(t_hybrid), INTENT(IN) :: hybrid
+   INTEGER,        INTENT(IN) :: vecLen, currentIter
+   REAL,           INTENT(IN) :: deltaNVec(vecLen)
+
+   INTEGER             :: mode, npos
+   INTEGER*8           :: recLen
+
+   CALL getIOMode(mode)
+
+   ! At the moment broyden IO is mode independent
+
+   recLen=(vecLen+1)*8
+
+   IF (hybrid%l_calhf) THEN
+      OPEN (59,file='hf_broyd_DN',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ELSE
+      OPEN (59,file='broyd_DN',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ENDIF
+
+   npos=currentIter-1
+   IF (currentIter.GT.input%maxiter+1) npos = MOD(currentIter-2,input%maxiter)+1
+
+   WRITE(59,rec=npos) deltaNVec(:vecLen)
+
+   CLOSE(59)
+
+END SUBROUTINE writeDeltaNVec
+
+SUBROUTINE readDeltaFVec(input,hybrid,vecLen,relIter,currentIter,deltaFVec)
+
+   TYPE(t_input),  INTENT(IN)  :: input
+   TYPE(t_hybrid), INTENT(IN)  :: hybrid
+   INTEGER,        INTENT(IN)  :: vecLen, relIter, currentIter
+   REAL,           INTENT(OUT) :: deltaFVec(vecLen)
+
+   INTEGER             :: mode, npos
+   INTEGER*8           :: recLen
+
+   CALL getIOMode(mode)
+
+   ! At the moment broyden IO is mode independent
+
+   recLen=(vecLen+1)*8
+
+   IF (hybrid%l_calhf) THEN
+      OPEN (59,file='hf_broyd_DF',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ELSE
+      OPEN (59,file='broyd_DF',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ENDIF
+
+   npos=currentIter+relIter-1
+   IF (currentIter.GT.input%maxiter+1) npos = MOD(currentIter-2,input%maxiter)+1
+
+   READ(59,rec=npos) deltaFVec(:vecLen)
+
+   CLOSE(59)
+
+END SUBROUTINE readDeltaFVec
+
+SUBROUTINE writeDeltaFVec(input,hybrid,vecLen,currentIter,deltaFVec)
+
+   TYPE(t_input),  INTENT(IN) :: input
+   TYPE(t_hybrid), INTENT(IN) :: hybrid
+   INTEGER,        INTENT(IN) :: vecLen, currentIter
+   REAL,           INTENT(IN) :: deltaFVec(vecLen)
+
+   INTEGER             :: mode, npos
+   INTEGER*8           :: recLen
+
+   CALL getIOMode(mode)
+
+   ! At the moment broyden IO is mode independent
+
+   recLen=(vecLen+1)*8
+
+   IF (hybrid%l_calhf) THEN
+      OPEN (59,file='hf_broyd_DF',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ELSE
+      OPEN (59,file='broyd_DF',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ENDIF
+
+   npos=currentIter-1
+   IF (currentIter.GT.input%maxiter+1) npos = MOD(currentIter-2,input%maxiter)+1
+
+   WRITE(59,rec=npos) deltaFVec(:vecLen)
+
+   CLOSE(59)
+
+END SUBROUTINE writeDeltaFVec
+
+
+SUBROUTINE writeBroydenOverlapExt(input,hybrid,currentIter,historyLength,&
+                                  dNdNLast,dFdFLast,dNdFLast,dFdNLast)
+
+   TYPE(t_input),  INTENT(IN) :: input
+   TYPE(t_hybrid), INTENT(IN) :: hybrid
+   INTEGER,        INTENT(IN) :: currentIter, historyLength
+   REAL,           INTENT(IN) :: dNdNLast(input%maxIter)
+   REAL,           INTENT(IN) :: dFdFLast(input%maxIter)
+   REAL,           INTENT(IN) :: dNdFLast(input%maxIter)
+   REAL,           INTENT(IN) :: dFdNLast(input%maxIter)
+
+   INTEGER                    :: recLen, npos
+
+   recLen = 8*4*input%maxIter    ! sizeOfReal*numberOfVectors*vectorLength
+   recLen = recLen + 2*8         ! storage for currentIter, historyLength
+
+   IF (hybrid%l_calhf) THEN
+      OPEN (59,file='hf_broydOvlp',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ELSE
+      OPEN (59,file='broydOvlp',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ENDIF
+
+   npos=currentIter-1
+   IF (currentIter.GT.input%maxiter+1) npos = MOD(currentIter-2,input%maxiter)+1
+
+   WRITE(59,rec=npos) currentIter, historyLength,&
+                      dNdNLast(:input%maxIter), dFdFLast(:input%maxIter), &
+                      dNdFLast(:input%maxIter), dFdNLast(:input%maxIter)
+
+   CLOSE(59)
+
+END SUBROUTINE writeBroydenOverlapExt
+
+SUBROUTINE readBroydenOverlaps(input,hybrid,currentIter,historyLength,&
+                               dNdNMat,dFdFMat,dNdFMat,dFdNMat)
+
+   TYPE(t_input),  INTENT(IN)    :: input
+   TYPE(t_hybrid), INTENT(IN)    :: hybrid
+   INTEGER,        INTENT(IN)    :: currentIter, historyLength
+   REAL,           INTENT(INOUT) :: dNdNMat(historyLength,historyLength)
+   REAL,           INTENT(INOUT) :: dFdFMat(historyLength,historyLength)
+   REAL,           INTENT(INOUT) :: dNdFMat(historyLength,historyLength)
+   REAL,           INTENT(INOUT) :: dFdNMat(historyLength,historyLength)
+
+   INTEGER                    :: i, j
+   INTEGER                    :: recLen, npos, iter
+   INTEGER                    :: recIter, recHistLen
+
+   REAL                       :: dNdNLast(input%maxIter)
+   REAL                       :: dFdFLast(input%maxIter)
+   REAL                       :: dNdFLast(input%maxIter)
+   REAL                       :: dFdNLast(input%maxIter)
+
+   recLen = 8*4*input%maxIter    ! sizeOfReal*numberOfVectors*vectorLength
+   recLen = recLen + 2*8         ! storage for currentIter, historyLength
+
+   IF (hybrid%l_calhf) THEN
+      OPEN (59,file='hf_broydOvlp',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ELSE
+      OPEN (59,file='broydOvlp',access='direct',&
+            recl=recLen,form='unformatted',status='unknown')
+   ENDIF
+
+   dNdNMat = 0.0
+   dFdFMat = 0.0
+   dNdFMat = 0.0
+   dFdNMat = 0.0
+
+   DO i = 1, historyLength
+
+      iter = currentIter - historyLength + i
+      npos=iter-1
+      IF (iter.GT.input%maxiter+1) npos = MOD(iter-2,input%maxiter)+1
+
+      READ(59,rec=npos) recIter, recHistLen,&
+                        dNdNLast(:input%maxIter), dFdFLast(:input%maxIter), &
+                        dNdFLast(:input%maxIter), dFdNLast(:input%maxIter)
+
+      DO j = 1, recHistLen
+         IF ((j-recHistLen+i).LE.0) CYCLE
+         dNdNMat(j-recHistLen+i,i) = dNdNLast(j)
+         dFdFMat(j-recHistLen+i,i) = dFdFLast(j)
+         dNdFMat(j-recHistLen+i,i) = dNdFLast(j)
+         dFdNMat(j-recHistLen+i,i) = dFdNLast(j)
+      END DO
+   END DO
+
+   CLOSE(59)
+
+END SUBROUTINE readBroydenOverlaps
+
+
+
+
+
+
+
 SUBROUTINE resetBroydenHistory()
 
    INTEGER           :: i

main/mix.F90

@@ -27,6 +27,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
    USE m_brysh1
    USE m_stmix
    USE m_broyden
+   USE m_broyden2
    USE m_brysh2
    USE m_metric
    USE m_qfix
@@ -155,6 +156,8 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
    ELSE
       CALL broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
                    hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
+!      CALL broyden2(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
+!                    hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
    END IF
 
    !initiatlize mixed density and extract it with brysh2 call

mix/CMakeLists.txt

@@ -3,6 +3,7 @@ mix/metr_z0.f
 )
 set(fleur_F90 ${fleur_F90}
 mix/broyden.F90
+mix/broyden2.F90
 mix/brysh1.f90
 mix/brysh2.f90
 mix/metric.f90

mix/broyden.F90

@@ -43,7 +43,7 @@ CONTAINS
     REAL,    INTENT (INOUT) :: sm(nmap)
 
     ! Local Scalars
-    INTEGER         :: i,it,k,nit,npos,iread,nmaph, mit
+    INTEGER         :: i,it,k,nit,iread,nmaph, mit
     REAL            :: bm,dfivi,fmvm,smnorm,vmnorm,alphan
     LOGICAL         :: l_pot, l_exist
     REAL, PARAMETER :: one=1.0, zero=0.0
@@ -109,13 +109,14 @@ CONTAINS
        !     |vi> = w|vi> 
        !     loop to generate um : um = sm1 + alpha*fm1 - \sum <fm1|w|vi> ui
        um(:nmap) = input%alpha * fm1(:nmap) + sm1(:nmap)
+
        iread = MIN(mit-1,input%maxiter+1)
        DO it = 2, iread
           CALL readUVec(input,hybrid,nmap,it-mit,mit,ui)
           CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)
 
           am(it) = CPP_BLAS_sdot(nmap,vi,1,fm1,1)
-          ! calculate um(:) = -am(it)*ui(:) + um
+          ! calculate um(:) = -am(it)*ui(:) + um(:)
           CALL CPP_BLAS_saxpy(nmap,-am(it),ui,1,um,1)
           WRITE(6,FMT='(5x,"<vi|w|Fm> for it",i2,5x,f10.6)')it,am(it) 
        END DO
@@ -134,6 +135,7 @@ CONTAINS
 
           ! generate vm = alpha*sm1  - \sum <ui|w|sm1> vi
           vm(:) = input%alpha * fm1(:)
+
           DO it = 2,iread
              CALL readUVec(input,hybrid,nmap,it-mit,mit,ui)
              CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)
@@ -155,7 +157,7 @@ CONTAINS
           !      broyden's second method
           !****************************************
 
-          ! multiply fm1 with metric matrix and store in vm:  w |fm1>  
+          ! multiply fm1 with metric matrix and store in vm:  w |fm1>
           CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
                       mmap,nmaph,mapmt,mapvac2,fm1,vm,l_pot)
 
@@ -192,16 +194,13 @@ CONTAINS
 
        ! write um,vm and dfivi on file broyd.?
 
-       npos=mit-1
-       IF (mit.GT.input%maxiter+1) npos = MOD(mit-2,input%maxiter)+1
-
        CALL writeUVec(input,hybrid,nmap,mit,um)
        CALL writeVVec(input,hybrid,nmap,mit,dfivi,vm)
 
        ! update rho(m+1)
        ! calculate <fm|w|vm>
        fmvm = CPP_BLAS_sdot(nmap,vm,1,fm,1)
-       ! calculate sm(:) = (1.0-fmvm)*ui(:) + sm
+       ! calculate sm(:) = (1.0-fmvm)*um(:) + sm
        CALL CPP_BLAS_saxpy(nmap,one-fmvm,um,1,sm,1)
     END IF
 

mix/broyden2.F90

+!--------------------------------------------------------------------------------
+! 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,&
+                      hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fm,sm,lpot)
+
+#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)
+
+    ! Local Scalars
+    INTEGER         :: i,j,it,k,nit,iread,nmaph, mit, historyLength
+    REAL            :: vFMetProd,alphan,coeff,vNorm
+    LOGICAL         :: l_pot, l_exist
+
+    ! Local Arrays
+    REAL, ALLOCATABLE :: fm1(:),sm1(:),uVec(:),vVec(:)
+    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
+
+    ALLOCATE (fm1(mmap),sm1(mmap),dNVec(mmap),dFVec(mmap),uVec(mmap),vVec(mmap))
+    ALLOCATE (deltaN_i(mmap),deltaF_i(mmap))
+    ALLOCATE (dNdNLast(input%maxiter),dFdFLast(input%maxiter))
+    ALLOCATE (dNdFLast(input%maxiter),dFdNLast(input%maxiter))
+
+    ALLOCATE (dnMet(mmap), dFMet(mmap), FMet(mmap))
+
+    fm1 = 0.0
+    sm1 = 0.0
+    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
+       ! load input charge density (sm1) and difference of 
+       ! in and out charge densities (fm1) from previous iteration (m-1)
+
+       CALL readLastIterInAndDiffDen(hybrid,nmap,mit,alphan,sm1(:nmap),fm1(:nmap))
+       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
+
+       ! generate F_m   - F_(m-1)  ... sm1
+       !      and rho_m - rho_(m-1) .. fm1
+       dNVec(1:nmap) = sm(1:nmap) - sm1(1:nmap)
+       dFVec(1:nmap) = fm(1:nmap) - fm1(1:nmap)
+    END IF
+
+    ! save F_m and rho_m for next iteration
+    nit = mit +1
+    IF (nit > input%maxiter+1) nit = 1
+    CALL writeLastIterInAndDiffDen(hybrid,nmap,nit,input%alpha,sm,fm)
+
+    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>  
+       CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
+                   mmap,nmaph,mapmt,mapvac2,dNVec,dNMet,l_pot)
+
+       ! Apply metric w to dFVec and store in dFMet:  w |dFVec>  
+       CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
+                   mmap,nmaph,mapmt,mapvac2,dFVec,dFMet,l_pot)
+
+       ! 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)
+