Move metric calls from the broyden2 to the mix subroutine

This implies the reduction of required metric calls in the broyden code path
from 4 to 2. This is now equivalent to the number of metric calls in the
old broyden code path.

io/broyd_io.F90

@@ -13,13 +13,14 @@ IMPLICIT NONE
 
 CONTAINS
 
-SUBROUTINE readLastIterInAndDiffDen(hybrid,vecLen,nextIter,alpha,inDenVec,diffDenVec)
+SUBROUTINE readLastIterInAndDiffDen(hybrid,vecLen,nextIter,alpha,inDenVec,diffDenVec,inDenVecMet,diffDenVecMet)
 
    TYPE(t_hybrid), INTENT(IN)  :: hybrid
    INTEGER,        INTENT(IN)  :: vecLen
    INTEGER,        INTENT(OUT) :: nextIter
    REAL,           INTENT(OUT) :: alpha
    REAL,           INTENT(OUT) :: inDenVec(vecLen), diffDenVec(vecLen)
+   REAL,OPTIONAL,  INTENT(OUT) :: inDenVecMet(vecLen), diffDenVecMet(vecLen)
 
    INTEGER             :: mode
 
@@ -33,18 +34,23 @@ SUBROUTINE readLastIterInAndDiffDen(hybrid,vecLen,nextIter,alpha,inDenVec,diffDe
       OPEN (57,file='broyd',form='unformatted',status='unknown')
    ENDIF
 
-   READ(57) nextIter,alpha,diffDenVec,inDenVec
+   IF(PRESENT(inDenVecMet)) THEN
+      READ(57) nextIter,alpha,diffDenVec,inDenVec,inDenVecMet,diffDenVecMet
+   ELSE
+      READ(57) nextIter,alpha,diffDenVec,inDenVec
+   END IF
 
    CLOSE(57)
 
 END SUBROUTINE readLastIterInAndDiffDen
 
-SUBROUTINE writeLastIterInAndDiffDen(hybrid,vecLen,nextIter,alpha,inDenVec,diffDenVec)
+SUBROUTINE writeLastIterInAndDiffDen(hybrid,vecLen,nextIter,alpha,inDenVec,diffDenVec,inDenVecMet,diffDenVecMet)
 
    TYPE(t_hybrid), INTENT(IN) :: hybrid
    INTEGER,        INTENT(IN) :: nextIter,vecLen
    REAL,           INTENT(IN) :: alpha
    REAL,           INTENT(IN) :: inDenVec(vecLen), diffDenVec(vecLen)
+   REAL,OPTIONAL,  INTENT(IN) :: inDenVecMet(vecLen), diffDenVecMet(vecLen)
 
    INTEGER             :: mode
 
@@ -58,7 +64,11 @@ SUBROUTINE writeLastIterInAndDiffDen(hybrid,vecLen,nextIter,alpha,inDenVec,diffD
       OPEN (57,file='broyd',form='unformatted',status='unknown')
    ENDIF
 
-   WRITE(57) nextIter,alpha,diffDenVec,inDenVec
+   IF(PRESENT(inDenVecMet)) THEN
+      WRITE(57) nextIter,alpha,diffDenVec,inDenVec,inDenVecMet,diffDenVecMet
+   ELSE
+      WRITE(57) nextIter,alpha,diffDenVec,inDenVec
+   END IF
 
    CLOSE(57)
 

main/mix.F90

@@ -56,11 +56,11 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
    INTEGER i,imap,js
    INTEGER mmap,mmaph,nmaph,nmap,mapmt,mapvac,mapvac2
    INTEGER iofl,n_u_keep
-   LOGICAL l_exist,l_ldaU, l_densityMatrixPresent
+   LOGICAL l_exist,l_ldaU, l_densityMatrixPresent, l_pot
 
    !Local Arrays
    REAL dist(6)
-   REAL, ALLOCATABLE :: sm(:), fsm(:)
+   REAL, ALLOCATABLE :: sm(:), fsm(:), fmMet(:), smMet(:)
    CHARACTER(LEN=20) :: attributes(2)
    COMPLEX           :: n_mmpTemp(-3:3,-3:3,MAX(1,atoms%n_u),input%jspins)
 
@@ -116,6 +116,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
    ! LDA+U (end)
 
    ALLOCATE (sm(mmap),fsm(mmap))
+   ALLOCATE (smMet(mmap),fmMet(mmap))
    dist(:) = 0.0
 
    !determine type of mixing:
@@ -150,14 +151,26 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
    !store the difference fsm - sm in fsm
    fsm(:nmap) = fsm(:nmap) - sm(:nmap)
 
+   l_pot = .FALSE.
+   ! Apply metric w to fsm and store in fmMet:  w |fsm>
+   CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
+               mmap,nmaph,mapmt,mapvac2,fsm,fmMet,l_pot)
+
    !mixing of the densities
    IF (input%imix.EQ.0) THEN
       CALL stmix(atoms,input,noco, nmap,nmaph,fsm, sm)
    ELSE
       CALL broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
                    hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
+
+!     Replace the broyden call above by the commented metric and broyden2 calls
+!     below to switch on the continuous restart of the Broyden method.
+!      ! Apply metric w to sm and store in smMet:  w |sm>
+!      CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
+!                  mmap,nmaph,mapmt,mapvac2,sm,smMet,l_pot)
+!
 !      CALL broyden2(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
-!                    hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
+!                    hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm,fmMet,smMet)
    END IF
 
    !initiatlize mixed density and extract it with brysh2 call
@@ -171,8 +184,8 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
    !calculate the distance of charge densities...
 
    !induce metric in fsm use sm as an output array: |sm> = w |fsm>
-   CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
-               mmap,nmaph,mapmt,mapvac2,fsm, sm)
+!   CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
+!               mmap,nmaph,mapmt,mapvac2,fsm, sm)
 
    !calculate the charge density distance for each spin
    IF(hybrid%l_calhf) THEN
@@ -182,7 +195,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
    END IF
 
    DO js = 1,input%jspins
-      dist(js) = CPP_BLAS_sdot(nmaph,fsm(nmaph*(js-1)+1),1, sm(nmaph*(js-1)+1),1)
+      dist(js) = CPP_BLAS_sdot(nmaph,fsm(nmaph*(js-1)+1),1, fmMet(nmaph*(js-1)+1),1)
 
       attributes = ''
       WRITE(attributes(1),'(i0)') js
@@ -196,14 +209,14 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
          WRITE ( 6,FMT=7900) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
       END IF
    END DO
-   IF (noco%l_noco) dist(6) = CPP_BLAS_sdot((nmap-2*nmaph), fsm(nmaph*2+1),1,sm(nmaph*2+1),1)
+   IF (noco%l_noco) dist(6) = CPP_BLAS_sdot((nmap-2*nmaph), fsm(nmaph*2+1),1,fmMet(nmaph*2+1),1)
    IF (noco%l_noco) WRITE (6,FMT=7900) 3,inDen%iter,1000*SQRT(ABS(dist(6)/cell%vol))
 
    !calculate the distance of total charge and spin density
    !|rho/m(o) - rho/m(i)| = |rh1(o) -rh1(i)|+ |rh2(o) -rh2(i)| +/_
    !                        +/_2<rh2(o) -rh2(i)|rh1(o) -rh1(i)>
    IF (input%jspins.EQ.2) THEN
-      dist(3) = CPP_BLAS_sdot(nmaph,fsm,1,sm(nmaph+1),1)
+      dist(3) = CPP_BLAS_sdot(nmaph,fsm,1,fmMet(nmaph+1),1)
       dist(4) = dist(1) + dist(2) + 2.0e0*dist(3)
       dist(5) = dist(1) + dist(2) - 2.0e0*dist(3)
       CALL writeXMLElementFormPoly('overallChargeDensity',(/'distance'/),&
@@ -227,7 +240,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
       !(e.g. when calculating non-magnetic systems with jspins=2).
    END IF
    results%last_distance=maxval(1000*SQRT(ABS(dist/cell%vol)))
-   DEALLOCATE (sm,fsm)
+   DEALLOCATE (sm,fsm,smMet,fmMet)
    CALL closeXMLElement('densityConvergence')
 
    !fix charge of the new density

mix/broyden2.F90

@@ -29,7 +29,7 @@ MODULE m_broyden2
   !################################################################
 CONTAINS
   SUBROUTINE broyden2(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
-                      hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fm,sm,lpot)
+                      hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fm,sm,fmMet,smMet,lpot)
 
 #include"cpp_double.h"
 
@@ -58,6 +58,8 @@ CONTAINS
     ! Array Arguments
     REAL,    INTENT (IN)    :: fm(nmap) 
     REAL,    INTENT (INOUT) :: sm(nmap)
+    REAL,    INTENT (IN)    :: fmMet(nmap) 
+    REAL,    INTENT (INOUT) :: smMet(nmap)
 
     ! Local Scalars
     INTEGER         :: i,j,it,k,nit,iread,nmaph, mit, historyLength
@@ -65,7 +67,8 @@ CONTAINS
     LOGICAL         :: l_pot, l_exist
 
     ! Local Arrays
-    REAL, ALLOCATABLE :: fm1(:),sm1(:),uVec(:),vVec(:)
+    REAL, ALLOCATABLE :: fmLast(:),smLast(:),fmLastMet(:),smLastMet(:)
+    REAL, ALLOCATABLE :: uVec(:),vVec(:)
     REAL, ALLOCATABLE :: dNVec(:),dFVec(:),dNMet(:), dFMet(:), FMet(:)
     REAL, ALLOCATABLE :: deltaN_i(:), deltaF_i(:)
     REAL, ALLOCATABLE :: dNdNMat(:,:), dNdFMat(:,:), dFdNMat(:,:), dFdFMat(:,:)
@@ -83,15 +86,16 @@ CONTAINS
     l_pot = .FALSE.
     IF (PRESENT(lpot)) l_pot = lpot
 
-    ALLOCATE (fm1(mmap),sm1(mmap),dNVec(mmap),dFVec(mmap),uVec(mmap),vVec(mmap))
+    ALLOCATE (fmLast(mmap),smLast(mmap),fmLastMet(mmap),smLastMet(mmap))
+    ALLOCATE (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))
+    ALLOCATE (dNMet(mmap), dFMet(mmap), FMet(mmap))
 
-    fm1 = 0.0
-    sm1 = 0.0
+    fmLast = 0.0
+    smLast = 0.0
     dNVec = 0.0
     dFVec = 0.0
     vVec  = 0.0
@@ -103,26 +107,29 @@ CONTAINS
     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)
+       ! load input charge density (smLast) and difference of 
+       ! in and out charge densities (fmLast) from previous iteration (m-1)
 
-       CALL readLastIterInAndDiffDen(hybrid,nmap,mit,alphan,sm1(:nmap),fm1(:nmap))
+       CALL readLastIterInAndDiffDen(hybrid,nmap,mit,alphan,smLast(:nmap),fmLast(:nmap),smLastMet(:nmap),fmLastMet(: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)
+       ! 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)
     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)
+    CALL writeLastIterInAndDiffDen(hybrid,nmap,nit,input%alpha,sm,fm,smMet,fmMet)
 
     IF (mit.EQ.1) THEN 
        !     update for rho for mit=1 is straight mixing
@@ -134,12 +141,12 @@ CONTAINS
        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)
+!       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)
+!       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-
@@ -339,17 +346,17 @@ CONTAINS
        ! Construct the final uVec and vVec -end-
 
        ! Apply metric w to dFVec and store in FMet:  w |dFVec>
-       CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
-                   mmap,nmaph,mapmt,mapvac2,fm,FMet,l_pot)
+!       CALL metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
+!                   mmap,nmaph,mapmt,mapvac2,fm,FMet,l_pot)
 
-       vFMetProd = CPP_BLAS_sdot(nmap,vVec,1,FMet,1)
+       vFMetProd = CPP_BLAS_sdot(nmap,vVec,1,fmMet,1)
 
        ! 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
 
-    DEALLOCATE (fm1,sm1,dNVec,dFVec,vVec)
+    DEALLOCATE (fmLast,smLast,fmLastMet,smLastMet,dNVec,dFVec,vVec)
 
   END SUBROUTINE broyden2
 END MODULE m_broyden2