Merge branch 'develop' of iffgit.fz-juelich.de:fleur/fleur into develop

diagonalization/elpa_cpp.F90

@@ -41,6 +41,7 @@ IMPLICIT NONE
     CPP_REALCOMPLEX, ALLOCATABLE :: eigvec(:,:)
     INTEGER, EXTERNAL :: numroc, indxl2g  !SCALAPACK functions
 #ifdef CPP_ELPA_201705003    
+    INTEGER :: kernel
     CLASS(elpa_t),pointer :: elpa_obj
 
     err = elpa_init(20170403)
@@ -65,6 +66,10 @@ IMPLICIT NONE
     CALL elpa_obj%set("solver", ELPA_SOLVER_1STAGE)
 #endif
     err = elpa_obj%setup()
+    if (myid == 0) then
+       call elpa_obj%get("complex_kernel", kernel)
+       print *, "elpa uses " // elpa_int_value_to_string("complex_kernel", kernel) // " kernel"
+    endif
 #endif
 
 

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
@@ -7,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
 
@@ -27,18 +34,23 @@ SUBROUTINE readLastIterInAndDiffDen(hybrid,vecLen,nextIter,alpha,inDenVec,diffDe
       OPEN (57,file='broyd',form='unformatted',status='unknown')
    ENDIF
 
+   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
 
@@ -52,7 +64,11 @@ SUBROUTINE writeLastIterInAndDiffDen(hybrid,vecLen,nextIter,alpha,inDenVec,diffD
       OPEN (57,file='broyd',form='unformatted',status='unknown')
    ENDIF
 
+   IF(PRESENT(inDenVecMet)) THEN
+      WRITE(57) nextIter,alpha,diffDenVec,inDenVec,inDenVecMet,diffDenVecMet
+   ELSE
       WRITE(57) nextIter,alpha,diffDenVec,inDenVec
+   END IF
 
    CLOSE(57)
 
@@ -192,6 +208,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
@@ -55,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)
 
@@ -115,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:
@@ -149,12 +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,fmMet,smMet)
    END IF
 
    !initiatlize mixed density and extract it with brysh2 call
@@ -168,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
@@ -179,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
@@ -193,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'/),&
@@ -224,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/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)