root / src / Step_GEDIIS_All.f90 @ 1
Historique | Voir | Annoter | Télécharger (21,02 ko)
| 1 | 1 | equemene | ! Geom = input parameter vector (Geometry), Grad = input gradient vector, HEAT is Energy(Geom) |
|---|---|---|---|
| 2 | 1 | equemene | SUBROUTINE Step_GEDIIS_All(NGeomF,IGeom,Step,Geom,Grad,HEAT,Hess,NCoord,allocation_flag,Tangent) |
| 3 | 1 | equemene | !SUBROUTINE Step_GEDIIS(Geom_new,Geom,Grad,HEAT,Hess,NCoord,FRST) |
| 4 | 1 | equemene | use Io_module |
| 5 | 1 | equemene | use Path_module, only : Nom, Atome, OrderInv, indzmat, Pi, Nat, Vfree |
| 6 | 1 | equemene | IMPLICIT NONE |
| 7 | 1 | equemene | |
| 8 | 1 | equemene | INTEGER(KINT) :: NGeomF,IGeom |
| 9 | 1 | equemene | INTEGER(KINT), INTENT(IN) :: NCoord |
| 10 | 1 | equemene | REAL(KREAL) :: Geom(NCoord), Grad(NCoord), Hess(NCoord*NCoord), Step(NCoord), Geom_new(NCoord) |
| 11 | 1 | equemene | REAL(KREAL) :: HEAT ! HEAT= Energy |
| 12 | 1 | equemene | LOGICAL :: allocation_flag |
| 13 | 1 | equemene | REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord) |
| 14 | 1 | equemene | |
| 15 | 1 | equemene | ! MRESET = maximum number of iterations. |
| 16 | 1 | equemene | INTEGER(KINT), PARAMETER :: MRESET=15, M2=(MRESET+1)*(MRESET+1) !M2 = 256 |
| 17 | 1 | equemene | REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:,:), GradSet(:,:) ! NGeomF,MRESET*NCoord |
| 18 | 1 | equemene | REAL(KREAL), ALLOCATABLE, SAVE :: GSAVE(:,:) !NGeomF,NCoord |
| 19 | 1 | equemene | REAL(KREAL), ALLOCATABLE, SAVE :: ESET(:,:) |
| 20 | 1 | equemene | REAL(KREAL) :: ESET_tmp(MRESET), B(M2), BS(M2), BST(M2), B_tmp(M2) ! M2=256 |
| 21 | 1 | equemene | LOGICAL :: DEBUG, PRINT, ci_lt_zero |
| 22 | 1 | equemene | INTEGER(KINT), ALLOCATABLE, SAVE :: MSET(:) ! mth Iteration |
| 23 | 1 | equemene | LOGICAL, ALLOCATABLE, SAVE :: FRST(:) |
| 24 | 1 | equemene | REAL(KREAL) :: ci(MRESET), ci_tmp(MRESET) ! MRESET = maximum number of iterations. |
| 25 | 1 | equemene | INTEGER(KINT) :: NGEDIIS, MPLUS, INV, ITERA, MM, cis_zero, IXX, JXX, MSET_minus_cis_zero |
| 26 | 1 | equemene | INTEGER(KINT) :: I,J,K, JJ, KJ, JNV, II, IONE, IJ, INK,ITmp, IX, JX, KX, MSET_C_cis_zero |
| 27 | 1 | equemene | INTEGER(KINT) :: current_size_B_mat, MyPointer, Iat, Isch, NFree, Idx |
| 28 | 1 | equemene | REAL(KREAL) :: XMax, XNorm, S, DET, THRES, tmp, ER_star, ER_star_tmp, Norm |
| 29 | 1 | equemene | REAL(KREAL), PARAMETER :: eps=1e-12 |
| 30 | 1 | equemene | REAL(KREAL), PARAMETER :: crit=1e-8 |
| 31 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord |
| 32 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: HFree(:) ! NFree*NFree |
| 33 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree |
| 34 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: Grad_free(:), Step_free(:) ! NFree |
| 35 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: Geom_free(:), Geom_new_free(:) ! NFree |
| 36 | 1 | equemene | REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet_free(:,:), GradSet_free(:,:) |
| 37 | 1 | equemene | |
| 38 | 1 | equemene | DEBUG=.TRUE. |
| 39 | 1 | equemene | PRINT=.FALSE. |
| 40 | 1 | equemene | |
| 41 | 1 | equemene | IF (PRINT) WRITE(*,'(/,'' BEGIN Step_GEDIIS_ALL '')') |
| 42 | 1 | equemene | |
| 43 | 1 | equemene | ! Initialization |
| 44 | 1 | equemene | IF (allocation_flag) THEN |
| 45 | 1 | equemene | ! allocation_flag will be set to False in SPACE_GEDIIS, so no need to modify it here |
| 46 | 1 | equemene | IF (ALLOCATED(GeomSet)) THEN |
| 47 | 1 | equemene | IF (PRINT) WRITE(*,'(/,'' In allocation_flag, GEDIIS_ALL Dealloc '')') |
| 48 | 1 | equemene | DEALLOCATE(GeomSet,GradSet,GSave,GeomSet_free,GradSet_free) |
| 49 | 1 | equemene | RETURN |
| 50 | 1 | equemene | ELSE |
| 51 | 1 | equemene | IF (PRINT) WRITE(*,'(/,'' In allocation_flag, GEDIIS_ALL Alloc '')') |
| 52 | 1 | equemene | ALLOCATE(GeomSet(NGeomF,MRESET*NCoord),GradSet(NGeomF,MRESET*NCoord),GSAVE(NGeomF,NCoord)) |
| 53 | 1 | equemene | ALLOCATE(GeomSet_free(NGeomF,MRESET*NCoord),GradSet_free(NGeomF,MRESET*NCoord)) |
| 54 | 1 | equemene | ALLOCATE(MSET(NGeomF),FRST(NGeomF),ESET(NGeomF,MRESET)) |
| 55 | 1 | equemene | DO I=1,NGeomF |
| 56 | 1 | equemene | FRST(I) = .TRUE. |
| 57 | 1 | equemene | GeomSet(I,:) = 0.d0 |
| 58 | 1 | equemene | GradSet(I,:) = 0.d0 |
| 59 | 1 | equemene | GSAVE(I,:)=0.d0 |
| 60 | 1 | equemene | GeomSet_free(I,:) = 0.d0 |
| 61 | 1 | equemene | GradSet_free(I,:) = 0.d0 |
| 62 | 1 | equemene | END DO |
| 63 | 1 | equemene | MSET(:)=0 |
| 64 | 1 | equemene | END IF |
| 65 | 1 | equemene | allocation_flag = .FALSE. |
| 66 | 1 | equemene | END IF ! IF (allocation_flag) THEN |
| 67 | 1 | equemene | |
| 68 | 1 | equemene | ! ADDED FROM HERE: |
| 69 | 1 | equemene | Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord), as of now, an Identity matrix. |
| 70 | 1 | equemene | ! we orthogonalize Vfree to the tangent vector of this geom only if Tangent/=0.d0 |
| 71 | 1 | equemene | Norm=sqrt(dot_product(Tangent,Tangent)) |
| 72 | 1 | equemene | IF (Norm.GT.eps) THEN |
| 73 | 1 | equemene | ALLOCATE(Tanf(NCoord)) |
| 74 | 1 | equemene | |
| 75 | 1 | equemene | ! We normalize Tangent |
| 76 | 1 | equemene | Tangent=Tangent/Norm |
| 77 | 1 | equemene | |
| 78 | 1 | equemene | ! We convert Tangent into Vfree only displacements. This is useless for now (2007.Apr.23) |
| 79 | 1 | equemene | ! as Vfree=Id matrix but it will be usefull as soon as we introduce constraints. |
| 80 | 1 | equemene | DO I=1,NCoord |
| 81 | 1 | equemene | Tanf(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Tangent) |
| 82 | 1 | equemene | END DO |
| 83 | 1 | equemene | Tangent=0.d0 |
| 84 | 1 | equemene | DO I=1,NCoord |
| 85 | 1 | equemene | Tangent=Tangent+Tanf(I)*Vfree(:,I) |
| 86 | 1 | equemene | END DO |
| 87 | 1 | equemene | ! first we subtract Tangent from vfree |
| 88 | 1 | equemene | DO I=1,NCoord |
| 89 | 1 | equemene | Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent) |
| 90 | 1 | equemene | Vfree(:,I)=Vfree(:,I)-Norm*Tangent |
| 91 | 1 | equemene | END DO |
| 92 | 1 | equemene | |
| 93 | 1 | equemene | Idx=0. |
| 94 | 1 | equemene | ! Schmidt orthogonalization of the Vfree vectors |
| 95 | 1 | equemene | DO I=1,NCoord |
| 96 | 1 | equemene | ! We subtract the first vectors, we do it twice as the Schmidt procedure is not numerically stable. |
| 97 | 1 | equemene | DO Isch=1,2 |
| 98 | 1 | equemene | DO J=1,Idx |
| 99 | 1 | equemene | Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,J),(/NCoord/))) |
| 100 | 1 | equemene | Vfree(:,I)=Vfree(:,I)-Norm*Vfree(:,J) |
| 101 | 1 | equemene | END DO |
| 102 | 1 | equemene | END DO |
| 103 | 1 | equemene | Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,I),(/NCoord/))) |
| 104 | 1 | equemene | IF (Norm.GE.crit) THEN |
| 105 | 1 | equemene | Idx=Idx+1 |
| 106 | 1 | equemene | Vfree(:,Idx)=Vfree(:,I)/sqrt(Norm) |
| 107 | 1 | equemene | END IF |
| 108 | 1 | equemene | END DO |
| 109 | 1 | equemene | |
| 110 | 1 | equemene | IF (Idx/= NCoord-1) THEN |
| 111 | 1 | equemene | WRITE(*,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom |
| 112 | 1 | equemene | WRITE(IOOut,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom |
| 113 | 1 | equemene | STOP |
| 114 | 1 | equemene | END IF |
| 115 | 1 | equemene | |
| 116 | 1 | equemene | DEALLOCATE(Tanf) |
| 117 | 1 | equemene | NFree=Idx |
| 118 | 1 | equemene | ELSE ! Tangent =0, matches IF (Norm.GT.eps) THEN |
| 119 | 1 | equemene | if (debug) WRITE(*,*) "Tangent=0, using full displacement" |
| 120 | 1 | equemene | NFree=NCoord |
| 121 | 1 | equemene | END IF !IF (Norm.GT.eps) THEN |
| 122 | 1 | equemene | |
| 123 | 1 | equemene | if (debug) WRITE(*,*) 'DBG Step_GEDIIS_All, IGeom, NFree=', IGeom, NFree |
| 124 | 1 | equemene | |
| 125 | 1 | equemene | ! We now calculate the new step |
| 126 | 1 | equemene | ! we project the hessian onto the free vectors |
| 127 | 1 | equemene | ALLOCATE(HFree(NFree*NFree),Htmp(NCoord,NFree),Grad_free(NFree)) |
| 128 | 1 | equemene | ALLOCATE(Geom_free(NFree),Step_free(NFree),Geom_new_free(NFree)) |
| 129 | 1 | equemene | DO J=1,NFree |
| 130 | 1 | equemene | DO I=1,NCoord |
| 131 | 1 | equemene | Htmp(I,J)=0.d0 |
| 132 | 1 | equemene | DO K=1,NCoord |
| 133 | 1 | equemene | Htmp(I,J)=Htmp(I,J)+Hess(((I-1)*NCoord)+K)*Vfree(K,J) |
| 134 | 1 | equemene | END DO |
| 135 | 1 | equemene | END DO |
| 136 | 1 | equemene | END DO |
| 137 | 1 | equemene | DO J=1,NFree |
| 138 | 1 | equemene | DO I=1,NFree |
| 139 | 1 | equemene | HFree(I+((J-1)*NFree))=0.d0 |
| 140 | 1 | equemene | DO K=1,NCoord |
| 141 | 1 | equemene | HFree(I+((J-1)*NFree))=HFree(I+((J-1)*NFree))+Vfree(K,I)*Htmp(K,J) |
| 142 | 1 | equemene | END DO |
| 143 | 1 | equemene | END DO |
| 144 | 1 | equemene | END DO |
| 145 | 1 | equemene | |
| 146 | 1 | equemene | DO I=1,NFree |
| 147 | 1 | equemene | Grad_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Grad) |
| 148 | 1 | equemene | Geom_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Geom) |
| 149 | 1 | equemene | END DO |
| 150 | 1 | equemene | !ADDED ENDS HERE.*********************************************** |
| 151 | 1 | equemene | |
| 152 | 1 | equemene | ! SPACE_GEDIIS SIMPLY LOADS THE CURRENT VALUES OF Geom AND Grad INTO THE ARRAYS GeomSet |
| 153 | 1 | equemene | ! AND GradSet, MSET is set to zero and then 1 in SPACE_GEDIIS_All at first iteration. |
| 154 | 1 | equemene | CALL SPACE_GEDIIS_All(NGeomF,IGeom,MRESET,MSET,Geom,Grad,HEAT,NCoord,GeomSet,GradSet,ESET,FRST) |
| 155 | 1 | equemene | |
| 156 | 1 | equemene | IF (PRINT) WRITE(*,'(/,'' GEDIIS after SPACE_GEDIIS_ALL '')') |
| 157 | 1 | equemene | |
| 158 | 1 | equemene | DO J=1,MSet(IGeom) |
| 159 | 1 | equemene | DO K=1,NFree |
| 160 | 1 | equemene | GradSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),& |
| 161 | 1 | equemene | GradSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord)) |
| 162 | 1 | equemene | GeomSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),& |
| 163 | 1 | equemene | GeomSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord)) |
| 164 | 1 | equemene | END DO |
| 165 | 1 | equemene | END DO |
| 166 | 1 | equemene | |
| 167 | 1 | equemene | ! INITIALIZE SOME VARIABLES AND CONSTANTS: |
| 168 | 1 | equemene | NGEDIIS = MSET(IGeom) !MSET=mth iteration |
| 169 | 1 | equemene | MPLUS = MSET(IGeom) + 1 |
| 170 | 1 | equemene | MM = MPLUS * MPLUS |
| 171 | 1 | equemene | |
| 172 | 1 | equemene | ! CONSTRUCT THE GEDIIS MATRIX: |
| 173 | 1 | equemene | ! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)> |
| 174 | 1 | equemene | JJ=0 |
| 175 | 1 | equemene | INV=-NFree |
| 176 | 1 | equemene | DO I=1,MSET(IGeom) |
| 177 | 1 | equemene | INV=INV+NFree |
| 178 | 1 | equemene | JNV=-NFree |
| 179 | 1 | equemene | DO J=1,MSET(IGeom) |
| 180 | 1 | equemene | JNV=JNV+NFree |
| 181 | 1 | equemene | JJ = JJ + 1 |
| 182 | 1 | equemene | B(JJ)=0.D0 |
| 183 | 1 | equemene | DO K=1, NFree |
| 184 | 1 | equemene | B(JJ) = B(JJ) + (((GradSet_free(IGeom,INV+K)-GradSet_free(IGeom,JNV+K))* & |
| 185 | 1 | equemene | (GeomSet_free(IGeom,INV+K)-GeomSet_free(IGeom,JNV+K)))/2.D0) |
| 186 | 1 | equemene | END DO |
| 187 | 1 | equemene | END DO |
| 188 | 1 | equemene | END DO |
| 189 | 1 | equemene | |
| 190 | 1 | equemene | ! The following shifting is required to correct indices of B_ij elements in the GEDIIS matrix. |
| 191 | 1 | equemene | ! The correction is needed because the last coloumn of the matrix contains all 1 and one zero. |
| 192 | 1 | equemene | DO I=MSET(IGeom)-1,1,-1 |
| 193 | 1 | equemene | DO J=MSET(IGeom),1,-1 |
| 194 | 1 | equemene | B(I*MSET(IGeom)+J+I) = B(I*MSET(IGeom)+J) |
| 195 | 1 | equemene | END DO |
| 196 | 1 | equemene | END DO |
| 197 | 1 | equemene | |
| 198 | 1 | equemene | ! For the last row and last column of GEDIIS matrix: |
| 199 | 1 | equemene | DO I=1,MPLUS |
| 200 | 1 | equemene | B(MPLUS*I) = 1.D0 |
| 201 | 1 | equemene | B(MPLUS*MSET(IGeom)+I) = 1.D0 |
| 202 | 1 | equemene | END DO |
| 203 | 1 | equemene | B(MM) = 0.D0 |
| 204 | 1 | equemene | |
| 205 | 1 | equemene | DO I=1, MPLUS |
| 206 | 1 | equemene | !WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS)) |
| 207 | 1 | equemene | END DO |
| 208 | 1 | equemene | |
| 209 | 1 | equemene | ! ELIMINATE ERROR VECTORS WITH THE LARGEST NORM: |
| 210 | 1 | equemene | 80 CONTINUE |
| 211 | 1 | equemene | DO I=1,MM !MM = (MSET(IGeom)+1) * (MSET(IGeom)+1) |
| 212 | 1 | equemene | BS(I) = B(I) !just a copy of the original B (GEDIIS) matrix |
| 213 | 1 | equemene | END DO |
| 214 | 1 | equemene | |
| 215 | 1 | equemene | IF (NGEDIIS .NE. MSET(IGeom)) THEN |
| 216 | 1 | equemene | DO II=1,MSET(IGeom)-NGEDIIS |
| 217 | 1 | equemene | XMAX = -1.D10 |
| 218 | 1 | equemene | ITERA = 0 |
| 219 | 1 | equemene | DO I=1,MSET(IGeom) |
| 220 | 1 | equemene | XNORM = 0.D0 |
| 221 | 1 | equemene | INV = (I-1) * MPLUS |
| 222 | 1 | equemene | DO J=1,MSET(IGeom) |
| 223 | 1 | equemene | XNORM = XNORM + ABS(B(INV + J)) |
| 224 | 1 | equemene | END DO |
| 225 | 1 | equemene | IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN |
| 226 | 1 | equemene | XMAX = XNORM |
| 227 | 1 | equemene | ITERA = I |
| 228 | 1 | equemene | IONE = INV + I |
| 229 | 1 | equemene | ENDIF |
| 230 | 1 | equemene | END DO |
| 231 | 1 | equemene | |
| 232 | 1 | equemene | DO I=1,MPLUS |
| 233 | 1 | equemene | INV = (I-1) * MPLUS |
| 234 | 1 | equemene | DO J=1,MPLUS |
| 235 | 1 | equemene | JNV = (J-1) * MPLUS |
| 236 | 1 | equemene | IF (J.EQ.ITERA) B(INV + J) = 0.D0 |
| 237 | 1 | equemene | B(JNV + I) = B(INV + J) |
| 238 | 1 | equemene | END DO |
| 239 | 1 | equemene | END DO |
| 240 | 1 | equemene | B(IONE) = 1.0D0 |
| 241 | 1 | equemene | END DO |
| 242 | 1 | equemene | END IF ! matches IF (NGEDIIS .NE. MSET(IGeom)) THEN |
| 243 | 1 | equemene | |
| 244 | 1 | equemene | ! SCALE GEDIIS MATRIX BEFORE INVERSION: |
| 245 | 1 | equemene | DO I=1,MPLUS |
| 246 | 1 | equemene | II = MPLUS * (I-1) + I ! B(II)=diagonal elements of B matrix |
| 247 | 1 | equemene | GSAVE(IGeom,I) = 1.D0 / DSQRT(1.D-20+DABS(B(II))) |
| 248 | 1 | equemene | !Print *, 'GSAVE(',IGeom,',',I,')=', GSAVE(IGeom,I)
|
| 249 | 1 | equemene | END DO |
| 250 | 1 | equemene | GSAVE(IGeom,MPLUS) = 1.D0 |
| 251 | 1 | equemene | DO I=1,MPLUS |
| 252 | 1 | equemene | DO J=1,MPLUS |
| 253 | 1 | equemene | IJ = MPLUS * (I-1) + J |
| 254 | 1 | equemene | B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J) |
| 255 | 1 | equemene | END DO |
| 256 | 1 | equemene | END DO |
| 257 | 1 | equemene | |
| 258 | 1 | equemene | ! INVERT THE GEDIIS MATRIX B: |
| 259 | 1 | equemene | DO I=1, MPLUS |
| 260 | 1 | equemene | !WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS)) |
| 261 | 1 | equemene | END DO |
| 262 | 1 | equemene | |
| 263 | 1 | equemene | CALL MINV(B,MPLUS,DET) ! matrix inversion. |
| 264 | 1 | equemene | |
| 265 | 1 | equemene | DO I=1, MPLUS |
| 266 | 1 | equemene | !WRITE(*,'(10(1X,F20.16))') B((I-1)*MPLUS+1:I*(MPLUS)) |
| 267 | 1 | equemene | END DO |
| 268 | 1 | equemene | |
| 269 | 1 | equemene | DO I=1,MPLUS |
| 270 | 1 | equemene | DO J=1,MPLUS |
| 271 | 1 | equemene | IJ = MPLUS * (I-1) + J |
| 272 | 1 | equemene | B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J) |
| 273 | 1 | equemene | END DO |
| 274 | 1 | equemene | END DO |
| 275 | 1 | equemene | |
| 276 | 1 | equemene | ! COMPUTE THE NEW INTERPOLATED PARAMETER VECTOR (Geometry): |
| 277 | 1 | equemene | ci=0.d0 |
| 278 | 1 | equemene | ci_tmp=0.d0 |
| 279 | 1 | equemene | |
| 280 | 1 | equemene | ci_lt_zero= .FALSE. |
| 281 | 1 | equemene | DO I=1, MSET(IGeom) |
| 282 | 1 | equemene | DO J=1, MSET(IGeom) ! B matrix is read column-wise |
| 283 | 1 | equemene | ci(I)=ci(I)+B((J-1)*(MPLUS)+I)*ESET(IGeom,J) !ESET is energy set. |
| 284 | 1 | equemene | END DO |
| 285 | 1 | equemene | ci(I)=ci(I)+B((MPLUS-1)*(MPLUS)+I) |
| 286 | 1 | equemene | !Print *, 'NO ci < 0 yet, c(',I,')=', ci(I)
|
| 287 | 1 | equemene | IF((ci(I) .LT. 0.0D0) .OR. (ci(I) .GT. 1.0D0)) THEN |
| 288 | 1 | equemene | ci_lt_zero=.TRUE. |
| 289 | 1 | equemene | EXIT |
| 290 | 1 | equemene | END IF |
| 291 | 1 | equemene | END DO !matches DO I=1, MSET(IGeom) |
| 292 | 1 | equemene | |
| 293 | 1 | equemene | IF (ci_lt_zero) Then |
| 294 | 1 | equemene | cis_zero = 0 |
| 295 | 1 | equemene | ER_star = 0.D0 |
| 296 | 1 | equemene | ER_star_tmp = 1e32 |
| 297 | 1 | equemene | |
| 298 | 1 | equemene | ! 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. |
| 299 | 1 | equemene | JJ=0 |
| 300 | 1 | equemene | INV=-NFree |
| 301 | 1 | equemene | DO IX=1,MSET(IGeom) |
| 302 | 1 | equemene | INV=INV+NFree |
| 303 | 1 | equemene | JNV=-NFree |
| 304 | 1 | equemene | DO JX=1,MSET(IGeom) |
| 305 | 1 | equemene | JNV=JNV+NFree |
| 306 | 1 | equemene | JJ = JJ + 1 |
| 307 | 1 | equemene | BST(JJ)=0.D0 |
| 308 | 1 | equemene | DO KX=1, NFree |
| 309 | 1 | equemene | BST(JJ) = BST(JJ) + (((GradSet_free(IGeom,INV+KX)-GradSet_free(IGeom,JNV+KX))* & |
| 310 | 1 | equemene | (GeomSet_free(IGeom,INV+KX)-GeomSet_free(IGeom,JNV+KX)))/2.D0) |
| 311 | 1 | equemene | END DO |
| 312 | 1 | equemene | END DO |
| 313 | 1 | equemene | END DO |
| 314 | 1 | equemene | |
| 315 | 1 | equemene | 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 |
| 316 | 1 | equemene | !Print *, 'Entering into DO I=1, (2**MSET(IGeom))-2 loop, MSET(IGeom)=', MSET(IGeom), ', I=', I |
| 317 | 1 | equemene | ci(:)=1.D0 |
| 318 | 1 | equemene | ! find out which cis are zero in I: |
| 319 | 1 | equemene | DO IX=1, MSET(IGeom) |
| 320 | 1 | equemene | JJ=IAND(I, 2**(IX-1)) |
| 321 | 1 | equemene | IF(JJ .EQ. 0) Then |
| 322 | 1 | equemene | ci(IX)=0.D0 |
| 323 | 1 | equemene | END IF |
| 324 | 1 | equemene | END DO |
| 325 | 1 | equemene | |
| 326 | 1 | equemene | ci_lt_zero = .FALSE. |
| 327 | 1 | equemene | ! 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. |
| 328 | 1 | equemene | DO IX=1, MSET(IGeom)*MSET(IGeom) |
| 329 | 1 | equemene | B(IX) = BST(IX) !just a copy of the original B (GEDIIS) matrix |
| 330 | 1 | equemene | END DO |
| 331 | 1 | equemene | |
| 332 | 1 | equemene | ! Removal of KXth row and KXth column in order to accomodate cI to be zero: |
| 333 | 1 | equemene | current_size_B_mat=MSET(IGeom) |
| 334 | 1 | equemene | cis_zero = 0 |
| 335 | 1 | equemene | ! The bits of I (index of the upper loop 'DO I=1, (2**MSET(IGeom))-2'), gives which cis are zero. |
| 336 | 1 | equemene | DO KX=1, MSET(IGeom) ! searching for each bit of I (index of the upper loop 'DO I=1, (2**MSET(IGeom))-2') |
| 337 | 1 | equemene | IF (ci(KX) .EQ. 0.D0) Then !remove KXth row and KXth column |
| 338 | 1 | equemene | cis_zero = cis_zero + 1 |
| 339 | 1 | equemene | |
| 340 | 1 | equemene | ! First row removal: (B matrix is read column-wise) |
| 341 | 1 | equemene | JJ=0 |
| 342 | 1 | equemene | DO IX=1,current_size_B_mat ! columns reading |
| 343 | 1 | equemene | DO JX=1,current_size_B_mat ! rows reading |
| 344 | 1 | equemene | IF (JX .NE. KX) Then |
| 345 | 1 | equemene | JJ = JJ + 1 |
| 346 | 1 | equemene | B_tmp(JJ) = B((IX-1)*current_size_B_mat+JX) |
| 347 | 1 | equemene | END IF |
| 348 | 1 | equemene | END DO |
| 349 | 1 | equemene | END DO |
| 350 | 1 | equemene | |
| 351 | 1 | equemene | DO IX=1,current_size_B_mat*(current_size_B_mat-1) |
| 352 | 1 | equemene | B(IX) = B_tmp(IX) |
| 353 | 1 | equemene | END DO |
| 354 | 1 | equemene | |
| 355 | 1 | equemene | ! Now column removal: |
| 356 | 1 | equemene | JJ=0 |
| 357 | 1 | equemene | DO IX=1,KX-1 ! columns reading |
| 358 | 1 | equemene | DO JX=1,current_size_B_mat-1 ! rows reading |
| 359 | 1 | equemene | JJ = JJ + 1 |
| 360 | 1 | equemene | B_tmp(JJ) = B(JJ) |
| 361 | 1 | equemene | END DO |
| 362 | 1 | equemene | END DO |
| 363 | 1 | equemene | |
| 364 | 1 | equemene | DO IX=KX+1,current_size_B_mat |
| 365 | 1 | equemene | DO JX=1,current_size_B_mat-1 |
| 366 | 1 | equemene | JJ = JJ + 1 |
| 367 | 1 | equemene | B_tmp(JJ) = B(JJ+current_size_B_mat-1) |
| 368 | 1 | equemene | END DO |
| 369 | 1 | equemene | END DO |
| 370 | 1 | equemene | |
| 371 | 1 | equemene | DO IX=1,(current_size_B_mat-1)*(current_size_B_mat-1) |
| 372 | 1 | equemene | B(IX) = B_tmp(IX) |
| 373 | 1 | equemene | END DO |
| 374 | 1 | equemene | current_size_B_mat = current_size_B_mat - 1 |
| 375 | 1 | equemene | END IF ! matches IF (ci(KX) .EQ. 0.D0) Then !remove |
| 376 | 1 | equemene | END DO ! matches DO KX=1, MSET(IGeom) |
| 377 | 1 | equemene | |
| 378 | 1 | equemene | ! The following shifting is required to correct indices of B_ij elements in the GEDIIS matrix. |
| 379 | 1 | equemene | ! The correction is needed because the last coloumn and row of the matrix contains all 1 and one zero. |
| 380 | 1 | equemene | DO IX=MSET(IGeom)-cis_zero-1,1,-1 |
| 381 | 1 | equemene | DO JX=MSET(IGeom)-cis_zero,1,-1 |
| 382 | 1 | equemene | B(IX*(MSET(IGeom)-cis_zero)+JX+IX) = B(IX*(MSET(IGeom)-cis_zero)+JX) |
| 383 | 1 | equemene | END DO |
| 384 | 1 | equemene | END DO |
| 385 | 1 | equemene | |
| 386 | 1 | equemene | ! for last row and last column of GEDIIS matrix |
| 387 | 1 | equemene | DO IX=1,MPLUS-cis_zero |
| 388 | 1 | equemene | B((MPLUS-cis_zero)*IX) = 1.D0 |
| 389 | 1 | equemene | B((MPLUS-cis_zero)*(MSET(IGeom)-cis_zero)+IX) = 1.D0 |
| 390 | 1 | equemene | END DO |
| 391 | 1 | equemene | B((MPLUS-cis_zero) * (MPLUS-cis_zero)) = 0.D0 |
| 392 | 1 | equemene | |
| 393 | 1 | equemene | DO IX=1, MPLUS |
| 394 | 1 | equemene | !WRITE(*,'(10(1X,F20.4))') B((IX-1)*MPLUS+1:IX*(MPLUS)) |
| 395 | 1 | equemene | END DO |
| 396 | 1 | equemene | |
| 397 | 1 | equemene | ! ELIMINATE ERROR VECTORS WITH THE LARGEST NORM: |
| 398 | 1 | equemene | IF (NGEDIIS .NE. MSET(IGeom)) THEN |
| 399 | 1 | equemene | JX = min(MSET(IGeom)-NGEDIIS,MSET(IGeom)-cis_zero-1) |
| 400 | 1 | equemene | DO II=1,JX |
| 401 | 1 | equemene | XMAX = -1.D10 |
| 402 | 1 | equemene | ITERA = 0 |
| 403 | 1 | equemene | DO IX=1,MSET(IGeom)-cis_zero |
| 404 | 1 | equemene | XNORM = 0.D0 |
| 405 | 1 | equemene | INV = (IX-1) * (MPLUS-cis_zero) |
| 406 | 1 | equemene | DO J=1,MSET(IGeom)-cis_zero |
| 407 | 1 | equemene | XNORM = XNORM + ABS(B(INV + J)) |
| 408 | 1 | equemene | END DO |
| 409 | 1 | equemene | IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN |
| 410 | 1 | equemene | XMAX = XNORM |
| 411 | 1 | equemene | ITERA = IX |
| 412 | 1 | equemene | IONE = INV + IX |
| 413 | 1 | equemene | ENDIF |
| 414 | 1 | equemene | END DO |
| 415 | 1 | equemene | |
| 416 | 1 | equemene | DO IX=1,MPLUS-cis_zero |
| 417 | 1 | equemene | INV = (IX-1) * (MPLUS-cis_zero) |
| 418 | 1 | equemene | DO J=1,MPLUS-cis_zero |
| 419 | 1 | equemene | JNV = (J-1) * (MPLUS-cis_zero) |
| 420 | 1 | equemene | IF (J.EQ.ITERA) B(INV + J) = 0.D0 |
| 421 | 1 | equemene | B(JNV + IX) = B(INV + J) |
| 422 | 1 | equemene | END DO |
| 423 | 1 | equemene | END DO |
| 424 | 1 | equemene | B(IONE) = 1.0D0 |
| 425 | 1 | equemene | END DO |
| 426 | 1 | equemene | END IF ! matches IF (NGEDIIS .NE. MSET(IGeom)) THEN |
| 427 | 1 | equemene | |
| 428 | 1 | equemene | ! SCALE GEDIIS MATRIX BEFORE INVERSION: |
| 429 | 1 | equemene | DO IX=1,MPLUS-cis_zero |
| 430 | 1 | equemene | II = (MPLUS-cis_zero) * (IX-1) + IX ! B(II)=diagonal elements of B matrix |
| 431 | 1 | equemene | GSAVE(IGeom,IX) = 1.D0 / DSQRT(1.D-20+DABS(B(II))) |
| 432 | 1 | equemene | END DO |
| 433 | 1 | equemene | GSAVE(IGeom,MPLUS-cis_zero) = 1.D0 |
| 434 | 1 | equemene | DO IX=1,MPLUS-cis_zero |
| 435 | 1 | equemene | DO JX=1,MPLUS-cis_zero |
| 436 | 1 | equemene | IJ = (MPLUS-cis_zero) * (IX-1) + JX |
| 437 | 1 | equemene | B(IJ) = B(IJ) * GSAVE(IGeom,IX) * GSAVE(IGeom,JX) |
| 438 | 1 | equemene | END DO |
| 439 | 1 | equemene | END DO |
| 440 | 1 | equemene | |
| 441 | 1 | equemene | ! INVERT THE GEDIIS MATRIX B: |
| 442 | 1 | equemene | CALL MINV(B,MPLUS-cis_zero,DET) ! matrix inversion. |
| 443 | 1 | equemene | |
| 444 | 1 | equemene | DO IX=1,MPLUS-cis_zero |
| 445 | 1 | equemene | DO JX=1,MPLUS-cis_zero |
| 446 | 1 | equemene | IJ = (MPLUS-cis_zero) * (IX-1) + JX |
| 447 | 1 | equemene | B(IJ) = B(IJ) * GSAVE(IGeom,IX) * GSAVE(IGeom,JX) |
| 448 | 1 | equemene | END DO |
| 449 | 1 | equemene | END DO |
| 450 | 1 | equemene | |
| 451 | 1 | equemene | DO IX=1, MPLUS |
| 452 | 1 | equemene | !WRITE(*,'(10(1X,F20.4))') B((IX-1)*MPLUS+1:IX*(MPLUS)) |
| 453 | 1 | equemene | END DO |
| 454 | 1 | equemene | |
| 455 | 1 | equemene | ! ESET is rearranged to handle zero cis and stored in ESET_tmp: |
| 456 | 1 | equemene | JJ=0 |
| 457 | 1 | equemene | DO IX=1, MSET(IGeom) |
| 458 | 1 | equemene | IF (ci(IX) .NE. 0) Then |
| 459 | 1 | equemene | JJ=JJ+1 |
| 460 | 1 | equemene | ESET_tmp(JJ) = ESET(IGeom,IX) |
| 461 | 1 | equemene | END IF |
| 462 | 1 | equemene | END DO |
| 463 | 1 | equemene | |
| 464 | 1 | equemene | ! DETERMINATION OF nonzero cis: |
| 465 | 1 | equemene | MyPointer=1 |
| 466 | 1 | equemene | DO IX=1, MSET(IGeom)-cis_zero |
| 467 | 1 | equemene | tmp = 0.D0 |
| 468 | 1 | equemene | DO J=1, MSET(IGeom)-cis_zero ! B matrix is read column-wise |
| 469 | 1 | equemene | tmp=tmp+B((J-1)*(MPLUS-cis_zero)+IX)*ESET_tmp(J) |
| 470 | 1 | equemene | END DO |
| 471 | 1 | equemene | tmp=tmp+B((MPLUS-cis_zero-1)*(MPLUS-cis_zero)+IX) |
| 472 | 1 | equemene | IF((tmp .LT. 0.0D0) .OR. (tmp .GT. 1.0D0)) THEN |
| 473 | 1 | equemene | ci_lt_zero=.TRUE. |
| 474 | 1 | equemene | EXIT |
| 475 | 1 | equemene | ELSE |
| 476 | 1 | equemene | DO JX=MyPointer,MSET(IGeom) |
| 477 | 1 | equemene | IF (ci(JX) .NE. 0) Then |
| 478 | 1 | equemene | ci(JX) = tmp |
| 479 | 1 | equemene | MyPointer=JX+1 |
| 480 | 1 | equemene | EXIT |
| 481 | 1 | equemene | END IF |
| 482 | 1 | equemene | END DO |
| 483 | 1 | equemene | END IF |
| 484 | 1 | equemene | END DO !matches DO I=1, MSET(IGeom)-cis_zero |
| 485 | 1 | equemene | !Print *, 'Local set of cis, first 10:, MSET(IGeom)=', MSET(IGeom), ', I of (2**MSET(IGeom))-2=', I |
| 486 | 1 | equemene | !WRITE(*,'(10(1X,F20.4))') ci(1:MSET(IGeom)) |
| 487 | 1 | equemene | !Print *, 'Local set of cis ends:****************************************' |
| 488 | 1 | equemene | |
| 489 | 1 | equemene | ! new set of cis determined based on the lower energy (ER_star): |
| 490 | 1 | equemene | IF (.NOT. ci_lt_zero) Then |
| 491 | 1 | equemene | Call Energy_GEDIIS(MRESET,MSET(IGeom),ci,GeomSet_free(IGeom,:),GradSet_free(IGeom,:),ESET(IGeom,:),NFree,ER_star) |
| 492 | 1 | equemene | IF (ER_star .LT. ER_star_tmp) Then |
| 493 | 1 | equemene | ci_tmp=ci |
| 494 | 1 | equemene | ER_star_tmp = ER_star |
| 495 | 1 | equemene | END IF |
| 496 | 1 | equemene | END IF ! matches IF (.NOT. ci_lt_zero) Then |
| 497 | 1 | equemene | 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 |
| 498 | 1 | equemene | ci = ci_tmp |
| 499 | 1 | equemene | END IF! matches IF (ci_lt_zero) Then |
| 500 | 1 | equemene | |
| 501 | 1 | equemene | Print *, 'Final set of cis, first 10:***********************************' |
| 502 | 1 | equemene | WRITE(*,'(10(1X,F20.4))') ci(1:MSET(IGeom)) |
| 503 | 1 | equemene | Print *, 'Final set of cis ends:****************************************' |
| 504 | 1 | equemene | Geom_new_free(:) = 0.D0 |
| 505 | 1 | equemene | DO I=1, MSET(IGeom) |
| 506 | 1 | equemene | Geom_new_free(:) = Geom_new_free(:) + (ci(I)*GeomSet_free(IGeom,(I-1)*NFree+1:I*NFree)) !MPLUS=MSET(IGeom)+1 |
| 507 | 1 | equemene | ! R_(N+1)=R*+DeltaR: |
| 508 | 1 | equemene | DO J=1, NFree |
| 509 | 1 | equemene | tmp=0.D0 |
| 510 | 1 | equemene | DO K=1,NFree |
| 511 | 1 | equemene | ! this can be commented: |
| 512 | 1 | equemene | !tmp=tmp+HFree((J-1)*NFree+K)*GradSet_free(IGeom,(I-1)*NFree+K) ! If Hinv=.False., then we need to invert Hess |
| 513 | 1 | equemene | END DO |
| 514 | 1 | equemene | Geom_new_free(J) = Geom_new_free(J) - (ci(I)*tmp) |
| 515 | 1 | equemene | END DO |
| 516 | 1 | equemene | END DO |
| 517 | 1 | equemene | |
| 518 | 1 | equemene | Step_free(:) = Geom_new_free(:) - Geom_free(:) |
| 519 | 1 | equemene | |
| 520 | 1 | equemene | XNORM = SQRT(DOT_PRODUCT(Step_free,Step_free)) |
| 521 | 1 | equemene | IF (PRINT) THEN |
| 522 | 1 | equemene | WRITE (6,'(/10X,''DEVIATION IN X '',F10.4,8X,''DETERMINANT '',G9.3)') XNORM, DET |
| 523 | 1 | equemene | !WRITE(*,'(10X,''GEDIIS COEFFICIENTS'')') |
| 524 | 1 | equemene | !WRITE(*,'(10X,5F12.5)') (B(MPLUS*MSET(IGeom)+I),I=1,MSET(IGeom)) |
| 525 | 1 | equemene | ENDIF |
| 526 | 1 | equemene | |
| 527 | 1 | equemene | ! THE FOLLOWING TOLERENCES FOR XNORM AND DET ARE SOMEWHAT ARBITRARY! |
| 528 | 1 | equemene | THRES = MAX(10.D0**(-NFree), 1.D-25) |
| 529 | 1 | equemene | IF (XNORM.GT.2.D0 .OR. DABS(DET) .LT. THRES) THEN |
| 530 | 1 | equemene | IF (PRINT)THEN |
| 531 | 1 | equemene | WRITE(*,*) "THE GEDIIS MATRIX IS ILL CONDITIONED" |
| 532 | 1 | equemene | WRITE(*,*) " - PROBABLY, VECTORS ARE LINEARLY DEPENDENT - " |
| 533 | 1 | equemene | WRITE(*,*) "THE GEDIIS STEP WILL BE REPEATED WITH A SMALLER SPACE" |
| 534 | 1 | equemene | END IF |
| 535 | 1 | equemene | DO K=1,MM |
| 536 | 1 | equemene | B(K) = BS(K) ! why this is reverted? Because "IF (NGEDIIS .GT. 0) GO TO 80", see below |
| 537 | 1 | equemene | END DO |
| 538 | 1 | equemene | NGEDIIS = NGEDIIS - 1 |
| 539 | 1 | equemene | IF (NGEDIIS .GT. 0) GO TO 80 |
| 540 | 1 | equemene | IF (PRINT) WRITE(*,'(10X,''NEWTON-RAPHSON STEP TAKEN'')') |
| 541 | 1 | equemene | Geom_new_free(:) = Geom_free(:) ! Geom_new is set to original Geom, thus Step = Geom(:) - Geom_new(:)=zero, the whole |
| 542 | 1 | equemene | ! new update to Geom_new is discarded, since XNORM.GT.2.D0 .OR. DABS(DET) .LT. THRES |
| 543 | 1 | equemene | END IF ! matches IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN |
| 544 | 1 | equemene | |
| 545 | 1 | equemene | !****************************************************************************************************************** |
| 546 | 1 | equemene | Geom_new_free(:) = 0.D0 |
| 547 | 1 | equemene | DO I=1, MSET(IGeom) |
| 548 | 1 | equemene | Geom_new_free(:) = Geom_new_free(:) + (ci(I)*GeomSet_free(IGeom,(I-1)*NFree+1:I*NFree)) !MPLUS=MSET(IGeom)+1 |
| 549 | 1 | equemene | ! R_(N+1)=R*+DeltaR: |
| 550 | 1 | equemene | DO J=1, NFree |
| 551 | 1 | equemene | tmp=0.D0 |
| 552 | 1 | equemene | DO K=1,NFree |
| 553 | 1 | equemene | tmp=tmp+HFree((J-1)*NFree+K)*GradSet_free(IGeom,(I-1)*NFree+K) ! If Hinv=.False., then we need to invert Hess |
| 554 | 1 | equemene | END DO |
| 555 | 1 | equemene | Geom_new_free(J) = Geom_new_free(J) - (ci(I)*tmp) |
| 556 | 1 | equemene | END DO |
| 557 | 1 | equemene | END DO |
| 558 | 1 | equemene | |
| 559 | 1 | equemene | Step_free(:) = Geom_new_free(:) - Geom_free(:) |
| 560 | 1 | equemene | !****************************************************************************************************************** |
| 561 | 1 | equemene | Step = 0.d0 |
| 562 | 1 | equemene | DO I=1,NFree |
| 563 | 1 | equemene | Step = Step + Step_free(I)*Vfree(:,I) |
| 564 | 1 | equemene | END DO |
| 565 | 1 | equemene | |
| 566 | 1 | equemene | DEALLOCATE(Hfree,Htmp,Grad_free,Step_free,Geom_free,Geom_new_free) |
| 567 | 1 | equemene | |
| 568 | 1 | equemene | IF (PRINT) WRITE(*,'(/,'' END Step_GEDIIS_ALL '',/)') |
| 569 | 1 | equemene | |
| 570 | 1 | equemene | END SUBROUTINE Step_GEDIIS_All |