     ! Geom = input parameter vector (Geometry), Grad = input gradient vector, HEAT is Energy(Geom)
      SUBROUTINE Step_GEDIIS_All(NGeomF,IGeom,Step,Geom,Grad,HEAT,Hess,NCoord,allocation_flag,Tangent)
      !SUBROUTINE Step_GEDIIS(Geom_new,Geom,Grad,HEAT,Hess,NCoord,FRST)
    use Io_module
    use Path_module, only : Nom, Atome, OrderInv, indzmat, Pi, Nat, Vfree
      IMPLICIT NONE

      INTEGER(KINT) :: NGeomF,IGeom
      INTEGER(KINT), INTENT(IN) :: NCoord
      REAL(KREAL) :: Geom(NCoord), Grad(NCoord), Hess(NCoord*NCoord), Step(NCoord), Geom_new(NCoord)
    REAL(KREAL) :: HEAT ! HEAT= Energy
    LOGICAL :: allocation_flag
    REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord)
    
      ! MRESET = maximum number of iterations.
      INTEGER(KINT), PARAMETER :: MRESET=15, M2=(MRESET+1)*(MRESET+1) !M2 = 256
      REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:,:), GradSet(:,:) ! NGeomF,MRESET*NCoord
      REAL(KREAL), ALLOCATABLE, SAVE :: GSAVE(:,:) !NGeomF,NCoord
    REAL(KREAL), ALLOCATABLE, SAVE :: ESET(:,:)
    REAL(KREAL) :: ESET_tmp(MRESET), B(M2), BS(M2), BST(M2), B_tmp(M2) ! M2=256
      LOGICAL :: DEBUG, PRINT, ci_lt_zero
      INTEGER(KINT), ALLOCATABLE, SAVE :: MSET(:) ! mth Iteration
    LOGICAL, ALLOCATABLE, SAVE :: FRST(:)
    REAL(KREAL) :: ci(MRESET), ci_tmp(MRESET) ! MRESET = maximum number of iterations.
      INTEGER(KINT) :: NGEDIIS, MPLUS, INV, ITERA, MM, cis_zero, IXX, JXX, MSET_minus_cis_zero
      INTEGER(KINT) :: I,J,K, JJ, KJ, JNV, II, IONE, IJ, INK,ITmp, IX, JX, KX, MSET_C_cis_zero
    INTEGER(KINT) :: current_size_B_mat, MyPointer, Iat, Isch, NFree, Idx
      REAL(KREAL) :: XMax, XNorm, S, DET, THRES, tmp, ER_star, ER_star_tmp, Norm
    REAL(KREAL), PARAMETER :: eps=1e-12
      REAL(KREAL), PARAMETER :: crit=1e-8
    REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord
    REAL(KREAL), ALLOCATABLE :: HFree(:) ! NFree*NFree
    REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree
    REAL(KREAL), ALLOCATABLE :: Grad_free(:), Step_free(:) ! NFree
    REAL(KREAL), ALLOCATABLE :: Geom_free(:), Geom_new_free(:) ! NFree
    REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet_free(:,:), GradSet_free(:,:)

      DEBUG=.TRUE.
      PRINT=.FALSE.
    
      IF (PRINT)  WRITE(*,'(/,''      BEGIN Step_GEDIIS_ALL   '')')
    
      ! Initialization
      IF (allocation_flag) THEN
      ! allocation_flag will be set to False in SPACE_GEDIIS, so no need to modify it here
         IF (ALLOCATED(GeomSet)) THEN
            IF (PRINT)  WRITE(*,'(/,''    In allocation_flag, GEDIIS_ALL  Dealloc  '')')
            DEALLOCATE(GeomSet,GradSet,GSave,GeomSet_free,GradSet_free)
            RETURN
         ELSE
            IF (PRINT)  WRITE(*,'(/,''     In allocation_flag,  GEDIIS_ALL  Alloc  '')')
            ALLOCATE(GeomSet(NGeomF,MRESET*NCoord),GradSet(NGeomF,MRESET*NCoord),GSAVE(NGeomF,NCoord))
      ALLOCATE(GeomSet_free(NGeomF,MRESET*NCoord),GradSet_free(NGeomF,MRESET*NCoord))
      ALLOCATE(MSET(NGeomF),FRST(NGeomF),ESET(NGeomF,MRESET))
      DO I=1,NGeomF
         FRST(I) = .TRUE.
         GeomSet(I,:) = 0.d0
         GradSet(I,:) = 0.d0
         GSAVE(I,:)=0.d0
         GeomSet_free(I,:) = 0.d0
         GradSet_free(I,:) = 0.d0
      END DO
      MSET(:)=0      
         END IF
     allocation_flag = .FALSE.
      END IF ! IF (allocation_flag) THEN
    
    ! ADDED FROM HERE:
    Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord), as of now, an Identity matrix.
      ! we orthogonalize Vfree to the tangent vector of this geom only if Tangent/=0.d0
      Norm=sqrt(dot_product(Tangent,Tangent))
      IF (Norm.GT.eps) THEN
         ALLOCATE(Tanf(NCoord))

         ! We normalize Tangent
         Tangent=Tangent/Norm

         ! We convert Tangent into Vfree only displacements. This is useless for now (2007.Apr.23)
     ! as Vfree=Id matrix but it will be usefull as soon as we introduce constraints.
         DO I=1,NCoord
            Tanf(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Tangent)
         END DO
         Tangent=0.d0
         DO I=1,NCoord
            Tangent=Tangent+Tanf(I)*Vfree(:,I)
         END DO
         ! first we subtract Tangent from vfree
         DO I=1,NCoord
            Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent)
            Vfree(:,I)=Vfree(:,I)-Norm*Tangent
     END DO

         Idx=0.
         ! Schmidt orthogonalization of the Vfree vectors
         DO I=1,NCoord
            ! We subtract the first vectors, we do it twice as the Schmidt procedure is not numerically stable.
            DO Isch=1,2
               DO J=1,Idx
                  Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,J),(/NCoord/)))
                  Vfree(:,I)=Vfree(:,I)-Norm*Vfree(:,J)
               END DO
            END DO
            Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,I),(/NCoord/)))
            IF (Norm.GE.crit) THEN
               Idx=Idx+1
               Vfree(:,Idx)=Vfree(:,I)/sqrt(Norm)
            END IF
         END DO
       
         IF (Idx/= NCoord-1) THEN
            WRITE(*,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom
            WRITE(IOOut,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom
            STOP
         END IF
     
         DEALLOCATE(Tanf)
         NFree=Idx
      ELSE ! Tangent =0, matches IF (Norm.GT.eps) THEN
         if (debug) WRITE(*,*) "Tangent=0, using full displacement"
         NFree=NCoord
      END IF !IF (Norm.GT.eps) THEN
  
      if (debug) WRITE(*,*) 'DBG Step_GEDIIS_All, IGeom, NFree=', IGeom, NFree

      ! We now calculate the new step
      ! we project the hessian onto the free vectors
      ALLOCATE(HFree(NFree*NFree),Htmp(NCoord,NFree),Grad_free(NFree))
    ALLOCATE(Geom_free(NFree),Step_free(NFree),Geom_new_free(NFree))
      DO J=1,NFree
        DO I=1,NCoord
           Htmp(I,J)=0.d0
           DO K=1,NCoord
              Htmp(I,J)=Htmp(I,J)+Hess(((I-1)*NCoord)+K)*Vfree(K,J)
           END DO
        END DO
      END DO
      DO J=1,NFree
        DO I=1,NFree
       HFree(I+((J-1)*NFree))=0.d0
           DO K=1,NCoord
        HFree(I+((J-1)*NFree))=HFree(I+((J-1)*NFree))+Vfree(K,I)*Htmp(K,J)
           END DO
        END DO
      END DO

      DO I=1,NFree
         Grad_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Grad)
     Geom_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Geom)
      END DO
      !ADDED ENDS HERE.***********************************************
    
      ! SPACE_GEDIIS SIMPLY LOADS THE CURRENT VALUES OF Geom AND Grad INTO THE ARRAYS GeomSet
      ! AND GradSet, MSET is set to zero and then 1 in SPACE_GEDIIS_All at first iteration.
      CALL SPACE_GEDIIS_All(NGeomF,IGeom,MRESET,MSET,Geom,Grad,HEAT,NCoord,GeomSet,GradSet,ESET,FRST)

      IF (PRINT)  WRITE(*,'(/,''       GEDIIS after SPACE_GEDIIS_ALL   '')')
    
    DO J=1,MSet(IGeom)
         DO K=1,NFree
            GradSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&
               GradSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
        GeomSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&
                GeomSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
         END DO
    END DO

      ! INITIALIZE SOME VARIABLES AND CONSTANTS:
      NGEDIIS = MSET(IGeom) !MSET=mth iteration
      MPLUS = MSET(IGeom) + 1
      MM = MPLUS * MPLUS

      ! CONSTRUCT THE GEDIIS MATRIX:
      ! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)>
      JJ=0
      INV=-NFree
      DO I=1,MSET(IGeom)
         INV=INV+NFree
         JNV=-NFree
         DO J=1,MSET(IGeom)
            JNV=JNV+NFree
            JJ = JJ + 1
            B(JJ)=0.D0
      DO K=1, NFree
         B(JJ) = B(JJ) + (((GradSet_free(IGeom,INV+K)-GradSet_free(IGeom,JNV+K))* &
                 (GeomSet_free(IGeom,INV+K)-GeomSet_free(IGeom,JNV+K)))/2.D0)
      END DO
         END DO
      END DO

     ! The following shifting is required to correct indices of B_ij elements in the GEDIIS matrix.
   ! The correction is needed because the last coloumn of the matrix contains all 1 and one zero.
      DO I=MSET(IGeom)-1,1,-1
         DO J=MSET(IGeom),1,-1
            B(I*MSET(IGeom)+J+I) = B(I*MSET(IGeom)+J) 
         END DO
    END DO
   
      ! For the last row and last column of GEDIIS matrix:
      DO I=1,MPLUS
         B(MPLUS*I) = 1.D0
         B(MPLUS*MSET(IGeom)+I) = 1.D0
      END DO
      B(MM) = 0.D0
    
    DO I=1, MPLUS
       !WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS))
    END DO
        
      ! ELIMINATE ERROR VECTORS WITH THE LARGEST NORM:
   80 CONTINUE
      DO I=1,MM !MM = (MSET(IGeom)+1) * (MSET(IGeom)+1)
         BS(I) = B(I) !just a copy of the original B (GEDIIS) matrix
      END DO
    
      IF (NGEDIIS .NE. MSET(IGeom)) THEN
        DO II=1,MSET(IGeom)-NGEDIIS
           XMAX = -1.D10
           ITERA = 0
           DO I=1,MSET(IGeom)
              XNORM = 0.D0
              INV = (I-1) * MPLUS
              DO J=1,MSET(IGeom)
                 XNORM = XNORM + ABS(B(INV + J))
              END DO
              IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN
                 XMAX = XNORM
                 ITERA = I
                 IONE = INV + I
              ENDIF
           END DO
     
           DO I=1,MPLUS
              INV = (I-1) * MPLUS
              DO J=1,MPLUS
                 JNV = (J-1) * MPLUS
                 IF (J.EQ.ITERA) B(INV + J) = 0.D0
                 B(JNV + I) = B(INV + J)
              END DO
       END DO
           B(IONE) = 1.0D0
        END DO
    END IF ! matches IF (NGEDIIS .NE. MSET(IGeom)) THEN
    
      ! SCALE GEDIIS MATRIX BEFORE INVERSION:
      DO I=1,MPLUS
         II = MPLUS * (I-1) + I ! B(II)=diagonal elements of B matrix
         GSAVE(IGeom,I) = 1.D0 / DSQRT(1.D-20+DABS(B(II)))
     !Print *, 'GSAVE(',IGeom,',',I,')=', GSAVE(IGeom,I) 
      END DO
      GSAVE(IGeom,MPLUS) = 1.D0
      DO I=1,MPLUS
         DO J=1,MPLUS
            IJ = MPLUS * (I-1) + J
            B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J)
         END DO
      END DO
  
     ! INVERT THE GEDIIS MATRIX B:
    DO I=1, MPLUS
       !WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS))
    END DO
        
      CALL MINV(B,MPLUS,DET) ! matrix inversion.
    
    DO I=1, MPLUS
       !WRITE(*,'(10(1X,F20.16))') B((I-1)*MPLUS+1:I*(MPLUS))
    END DO
    
      DO I=1,MPLUS
         DO J=1,MPLUS
            IJ = MPLUS * (I-1) + J
            B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J)
         END DO
      END DO
    
      ! COMPUTE THE NEW INTERPOLATED PARAMETER VECTOR (Geometry):    
    ci=0.d0
    ci_tmp=0.d0
  
    ci_lt_zero= .FALSE.
    DO I=1, MSET(IGeom)
     DO J=1, MSET(IGeom) ! B matrix is read column-wise
        ci(I)=ci(I)+B((J-1)*(MPLUS)+I)*ESET(IGeom,J) !ESET is energy set.                            
     END DO
     ci(I)=ci(I)+B((MPLUS-1)*(MPLUS)+I)
     !Print *, 'NO ci < 0 yet, c(',I,')=', ci(I)
     IF((ci(I) .LT. 0.0D0) .OR. (ci(I) .GT. 1.0D0)) THEN
       ci_lt_zero=.TRUE.
       EXIT
     END IF
      END DO !matches DO I=1, MSET(IGeom)
    
    IF (ci_lt_zero) Then
       cis_zero = 0
         ER_star = 0.D0
     ER_star_tmp = 1e32
     
     ! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)>, Full B matrix created first and then rows and columns are removed. 
         JJ=0
         INV=-NFree
         DO IX=1,MSET(IGeom)
            INV=INV+NFree
            JNV=-NFree
            DO JX=1,MSET(IGeom)
               JNV=JNV+NFree
               JJ = JJ + 1
               BST(JJ)=0.D0
           DO KX=1, NFree
              BST(JJ) = BST(JJ) + (((GradSet_free(IGeom,INV+KX)-GradSet_free(IGeom,JNV+KX))* &
                    (GeomSet_free(IGeom,INV+KX)-GeomSet_free(IGeom,JNV+KX)))/2.D0)
           END DO
            END DO
       END DO  
    
     DO I=1, (2**MSET(IGeom))-2 ! all (2**MSET(IGeom))-2 combinations of cis, except the one where all cis are .GT. 0 and .LT. 1
         !Print *, 'Entering into DO I=1, (2**MSET(IGeom))-2  loop, MSET(IGeom)=', MSET(IGeom), ', I=', I
         ci(:)=1.D0
         ! find out which cis are zero in I:
       DO IX=1, MSET(IGeom)
          JJ=IAND(I, 2**(IX-1))
        IF(JJ .EQ. 0) Then
          ci(IX)=0.D0
          END IF
       END DO
       
       ci_lt_zero = .FALSE.
       ! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)>, Full B matrix created first and then rows and columns are removed. 
       DO IX=1, MSET(IGeom)*MSET(IGeom)
                B(IX) = BST(IX) !just a copy of the original B (GEDIIS) matrix
             END DO
       
             ! Removal of KXth row and KXth column in order to accomodate cI to be zero:
       current_size_B_mat=MSET(IGeom)
       cis_zero = 0
       ! The bits of I (index of the upper loop 'DO I=1, (2**MSET(IGeom))-2'), gives which cis are zero.
       DO KX=1, MSET(IGeom) ! searching for each bit of I (index of the upper loop 'DO I=1, (2**MSET(IGeom))-2')
          IF (ci(KX) .EQ. 0.D0) Then !remove KXth row and KXth column
           cis_zero = cis_zero + 1
           
             ! First row removal: (B matrix is read column-wise)
             JJ=0
                   DO IX=1,current_size_B_mat ! columns reading
                      DO JX=1,current_size_B_mat ! rows reading
                 IF (JX .NE. KX) Then
                     JJ = JJ + 1
                     B_tmp(JJ) = B((IX-1)*current_size_B_mat+JX)
                 END IF
              END DO   
             END DO
       
             DO IX=1,current_size_B_mat*(current_size_B_mat-1)
                B(IX) = B_tmp(IX)
             END DO
       
             ! Now column removal:
             JJ=0   
                   DO IX=1,KX-1 ! columns reading
                      DO JX=1,current_size_B_mat-1 ! rows reading
                 JJ = JJ + 1
                 B_tmp(JJ) = B(JJ)
              END DO
             END DO        
        
                   DO IX=KX+1,current_size_B_mat
                      DO JX=1,current_size_B_mat-1
                 JJ = JJ + 1
                    B_tmp(JJ) = B(JJ+current_size_B_mat-1)
              END DO
             END DO
              
             DO IX=1,(current_size_B_mat-1)*(current_size_B_mat-1)
                B(IX) = B_tmp(IX)
             END DO
           current_size_B_mat = current_size_B_mat - 1
        END IF ! matches IF (ci(KX) .EQ. 0.D0) Then !remove
           END DO ! matches DO KX=1, MSET(IGeom)

       ! The following shifting is required to correct indices of B_ij elements in the GEDIIS matrix.
       ! The correction is needed because the last coloumn and row of the matrix contains all 1 and one zero.
       DO IX=MSET(IGeom)-cis_zero-1,1,-1
        DO JX=MSET(IGeom)-cis_zero,1,-1
           B(IX*(MSET(IGeom)-cis_zero)+JX+IX) = B(IX*(MSET(IGeom)-cis_zero)+JX) 
        END DO
       END DO
       
       ! for last row and last column of GEDIIS matrix 
       DO IX=1,MPLUS-cis_zero
        B((MPLUS-cis_zero)*IX) = 1.D0
        B((MPLUS-cis_zero)*(MSET(IGeom)-cis_zero)+IX) = 1.D0
       END DO
       B((MPLUS-cis_zero) * (MPLUS-cis_zero)) = 0.D0
         
           DO IX=1, MPLUS
              !WRITE(*,'(10(1X,F20.4))') B((IX-1)*MPLUS+1:IX*(MPLUS))
           END DO
              
       ! ELIMINATE ERROR VECTORS WITH THE LARGEST NORM:
             IF (NGEDIIS .NE. MSET(IGeom)) THEN
          JX = min(MSET(IGeom)-NGEDIIS,MSET(IGeom)-cis_zero-1)
                DO II=1,JX
                   XMAX = -1.D10
                   ITERA = 0
                   DO IX=1,MSET(IGeom)-cis_zero
                      XNORM = 0.D0
                      INV = (IX-1) * (MPLUS-cis_zero)
                      DO J=1,MSET(IGeom)-cis_zero
                         XNORM = XNORM + ABS(B(INV + J))
                      END DO
                      IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN
                         XMAX = XNORM
                         ITERA = IX
                         IONE = INV + IX
                      ENDIF
           END DO
     
                   DO IX=1,MPLUS-cis_zero
                      INV = (IX-1) * (MPLUS-cis_zero)
                      DO J=1,MPLUS-cis_zero
                         JNV = (J-1) * (MPLUS-cis_zero)
                         IF (J.EQ.ITERA) B(INV + J) = 0.D0
                         B(JNV + IX) = B(INV + J)
            END DO
               END DO
                   B(IONE) = 1.0D0
          END DO
           END IF ! matches IF (NGEDIIS .NE. MSET(IGeom)) THEN

       ! SCALE GEDIIS MATRIX BEFORE INVERSION:
       DO IX=1,MPLUS-cis_zero
        II = (MPLUS-cis_zero) * (IX-1) + IX ! B(II)=diagonal elements of B matrix
        GSAVE(IGeom,IX) = 1.D0 / DSQRT(1.D-20+DABS(B(II)))
       END DO
       GSAVE(IGeom,MPLUS-cis_zero) = 1.D0
       DO IX=1,MPLUS-cis_zero
        DO JX=1,MPLUS-cis_zero
           IJ = (MPLUS-cis_zero) * (IX-1) + JX
           B(IJ) = B(IJ) * GSAVE(IGeom,IX) * GSAVE(IGeom,JX)
        END DO
       END DO

       ! INVERT THE GEDIIS MATRIX B:
       CALL MINV(B,MPLUS-cis_zero,DET) ! matrix inversion.
    
       DO IX=1,MPLUS-cis_zero
        DO JX=1,MPLUS-cis_zero
           IJ = (MPLUS-cis_zero) * (IX-1) + JX
           B(IJ) = B(IJ) * GSAVE(IGeom,IX) * GSAVE(IGeom,JX)
        END DO
       END DO
       
           DO IX=1, MPLUS
              !WRITE(*,'(10(1X,F20.4))') B((IX-1)*MPLUS+1:IX*(MPLUS))
           END DO
              
             ! ESET is rearranged to handle zero cis and stored in ESET_tmp:
       JJ=0
       DO IX=1, MSET(IGeom)
        IF (ci(IX) .NE. 0) Then
           JJ=JJ+1
           ESET_tmp(JJ) = ESET(IGeom,IX)
        END IF
       END DO
       
       ! DETERMINATION OF nonzero cis:
       MyPointer=1
            DO IX=1, MSET(IGeom)-cis_zero
          tmp = 0.D0
            DO J=1, MSET(IGeom)-cis_zero ! B matrix is read column-wise
           tmp=tmp+B((J-1)*(MPLUS-cis_zero)+IX)*ESET_tmp(J)
        END DO
            tmp=tmp+B((MPLUS-cis_zero-1)*(MPLUS-cis_zero)+IX)
            IF((tmp .LT. 0.0D0) .OR. (tmp .GT. 1.0D0)) THEN
               ci_lt_zero=.TRUE.
               EXIT
        ELSE
           DO JX=MyPointer,MSET(IGeom)
              IF (ci(JX) .NE. 0) Then
                 ci(JX) = tmp
               MyPointer=JX+1
               EXIT
            END IF
           END DO
            END IF
             END DO !matches DO I=1, MSET(IGeom)-cis_zero
           !Print *, 'Local set of cis, first 10:, MSET(IGeom)=', MSET(IGeom), ', I of (2**MSET(IGeom))-2=', I
       !WRITE(*,'(10(1X,F20.4))') ci(1:MSET(IGeom))
           !Print *, 'Local set of cis ends:****************************************'

             ! new set of cis determined based on the lower energy (ER_star):
       IF (.NOT. ci_lt_zero) Then
                Call Energy_GEDIIS(MRESET,MSET(IGeom),ci,GeomSet_free(IGeom,:),GradSet_free(IGeom,:),ESET(IGeom,:),NFree,ER_star)
        IF (ER_star .LT. ER_star_tmp) Then
           ci_tmp=ci
           ER_star_tmp = ER_star
        END IF
             END IF ! matches IF (.NOT. ci_lt_zero) Then
          END DO !matches DO I=1, (2**K)-2 ! all (2**K)-2 combinations of cis, except the one where all cis are .GT. 0 and .LT. 1 
      ci = ci_tmp
    END IF! matches IF (ci_lt_zero) Then
    
    Print *, 'Final set of cis, first 10:***********************************'
    WRITE(*,'(10(1X,F20.4))') ci(1:MSET(IGeom))
    Print *, 'Final set of cis ends:****************************************'
    Geom_new_free(:) = 0.D0
    DO I=1, MSET(IGeom)    
         Geom_new_free(:) = Geom_new_free(:) + (ci(I)*GeomSet_free(IGeom,(I-1)*NFree+1:I*NFree)) !MPLUS=MSET(IGeom)+1
     ! R_(N+1)=R*+DeltaR:
     DO J=1, NFree
      tmp=0.D0
      DO K=1,NFree
         ! this can be commented:
         !tmp=tmp+HFree((J-1)*NFree+K)*GradSet_free(IGeom,(I-1)*NFree+K) ! If Hinv=.False., then we need to invert Hess
      END DO
      Geom_new_free(J) = Geom_new_free(J) - (ci(I)*tmp)
     END DO
    END DO
    
    Step_free(:) =  Geom_new_free(:) - Geom_free(:)
    
      XNORM = SQRT(DOT_PRODUCT(Step_free,Step_free))
      IF (PRINT) THEN
         WRITE (6,'(/10X,''DEVIATION IN X '',F10.4,8X,''DETERMINANT '',G9.3)') XNORM, DET
         !WRITE(*,'(10X,''GEDIIS COEFFICIENTS'')')
         !WRITE(*,'(10X,5F12.5)') (B(MPLUS*MSET(IGeom)+I),I=1,MSET(IGeom))
      ENDIF

      ! THE FOLLOWING TOLERENCES FOR XNORM AND DET ARE SOMEWHAT ARBITRARY!
      THRES = MAX(10.D0**(-NFree), 1.D-25)
      IF (XNORM.GT.2.D0 .OR. DABS(DET) .LT. THRES) THEN
         IF (PRINT)THEN
            WRITE(*,*) "THE GEDIIS MATRIX IS ILL CONDITIONED"     
            WRITE(*,*) " - PROBABLY, VECTORS ARE LINEARLY DEPENDENT - "
            WRITE(*,*) "THE GEDIIS STEP WILL BE REPEATED WITH A SMALLER SPACE"
         END IF
         DO K=1,MM
      B(K) = BS(K) ! why this is reverted? Because "IF (NGEDIIS .GT. 0) GO TO 80", see below
         END DO
         NGEDIIS = NGEDIIS - 1
         IF (NGEDIIS .GT. 0) GO TO 80
         IF (PRINT) WRITE(*,'(10X,''NEWTON-RAPHSON STEP TAKEN'')')
         Geom_new_free(:) = Geom_free(:) ! Geom_new is set to original Geom, thus Step = Geom(:) - Geom_new(:)=zero, the whole 
                           ! new update to Geom_new is discarded, since XNORM.GT.2.D0 .OR. DABS(DET) .LT. THRES
      END IF ! matches IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN
    
    !******************************************************************************************************************
    Geom_new_free(:) = 0.D0
    DO I=1, MSET(IGeom)    
         Geom_new_free(:) = Geom_new_free(:) + (ci(I)*GeomSet_free(IGeom,(I-1)*NFree+1:I*NFree)) !MPLUS=MSET(IGeom)+1
     ! R_(N+1)=R*+DeltaR:
     DO J=1, NFree
      tmp=0.D0
      DO K=1,NFree
         tmp=tmp+HFree((J-1)*NFree+K)*GradSet_free(IGeom,(I-1)*NFree+K) ! If Hinv=.False., then we need to invert Hess
      END DO
      Geom_new_free(J) = Geom_new_free(J) - (ci(I)*tmp)
     END DO
    END DO
    
    Step_free(:) =  Geom_new_free(:) - Geom_free(:)
   !******************************************************************************************************************
    Step = 0.d0
      DO I=1,NFree
     Step = Step + Step_free(I)*Vfree(:,I)
      END DO    
    
      DEALLOCATE(Hfree,Htmp,Grad_free,Step_free,Geom_free,Geom_new_free)
    
      IF (PRINT)  WRITE(*,'(/,''       END Step_GEDIIS_ALL   '',/)')

      END SUBROUTINE Step_GEDIIS_All

