Statistiques
| Révision :

## root / src / Step_GEDIIS_All.f90 @ 8

Historique | Voir | Annoter | Télécharger (21,64 ko)

 1  ! 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)   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   INTEGER(KINT) :: I, J, K, JJ, JNV, II, IONE, IJ, IX, JX, KX   INTEGER(KINT) :: current_size_B_mat, MyPointer, Isch, NFree, Idx   REAL(KREAL) :: XMax, XNorm, 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   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 , 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 , 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