root / src / Step_GEDIIS.f90 @ 1
Historique | Voir | Annoter | Télécharger (15,13 ko)
| 1 | 1 | equemene | ! Geom = input parameter vector (Geometry), Grad = input gradient vector. |
|---|---|---|---|
| 2 | 1 | equemene | ! HEAT is Energy(Geom) |
| 3 | 1 | equemene | SUBROUTINE Step_GEDIIS(Geom_new,Geom,Grad,HEAT,Hess,NCoord,FRST) |
| 4 | 1 | equemene | |
| 5 | 1 | equemene | use Io_module |
| 6 | 1 | equemene | use Path_module, only : Nom, Atome, OrderInv, indzmat, Pi, Nat |
| 7 | 1 | equemene | |
| 8 | 1 | equemene | IMPLICIT NONE |
| 9 | 1 | equemene | |
| 10 | 1 | equemene | INTEGER(KINT) :: NCoord |
| 11 | 1 | equemene | REAL(KREAL) :: Geom_new(NCoord), Grad(NCoord), Hess(NCoord*NCoord) |
| 12 | 1 | equemene | REAL(KREAL), INTENT(IN) :: Geom(NCoord) |
| 13 | 1 | equemene | REAL(KREAL) :: HEAT ! HEAT= Energy |
| 14 | 1 | equemene | LOGICAL :: FRST |
| 15 | 1 | equemene | |
| 16 | 1 | equemene | ! MRESET = maximum number of iterations. |
| 17 | 1 | equemene | INTEGER(KINT), PARAMETER :: MRESET=15, M2=(MRESET+1)*(MRESET+1) !M2 = 256 |
| 18 | 1 | equemene | REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:), GradSet(:) ! MRESET*NCoord |
| 19 | 1 | equemene | REAL(KREAL), ALLOCATABLE, SAVE :: DX(:), GSAVE(:) !NCoord |
| 20 | 1 | equemene | REAL(KREAL), SAVE :: ESET(MRESET) |
| 21 | 1 | equemene | REAL(KREAL) :: ESET_tmp(MRESET), B(M2),BS(M2),BST(M2), B_tmp(M2) ! M2=256 |
| 22 | 1 | equemene | LOGICAL DEBUG, PRINT, ci_lt_zero |
| 23 | 1 | equemene | INTEGER(KINT), SAVE :: MSET ! mth Iteration |
| 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 |
| 28 | 1 | equemene | REAL(KREAL) :: XMax, XNorm, S, DET, THRES, tmp, ER_star, ER_star_tmp |
| 29 | 1 | equemene | |
| 30 | 1 | equemene | DEBUG=.TRUE. |
| 31 | 1 | equemene | PRINT=.FALSE. |
| 32 | 1 | equemene | |
| 33 | 1 | equemene | IF (PRINT) WRITE(*,'(/,'' BEGIN GEDIIS '')') |
| 34 | 1 | equemene | |
| 35 | 1 | equemene | ! Initialization |
| 36 | 1 | equemene | IF (FRST) THEN |
| 37 | 1 | equemene | ! FRST will be set to False in SPACE_GEDIIS, so no need to modify it here |
| 38 | 1 | equemene | IF (ALLOCATED(GeomSet)) THEN |
| 39 | 1 | equemene | IF (PRINT) WRITE(*,'(/,'' In FRST, GEDIIS Dealloc '')') |
| 40 | 1 | equemene | DEALLOCATE(GeomSet,GradSet,DX,GSave) |
| 41 | 1 | equemene | RETURN |
| 42 | 1 | equemene | ELSE |
| 43 | 1 | equemene | IF (PRINT) WRITE(*,'(/,'' In FRST, GEDIIS Alloc '')') |
| 44 | 1 | equemene | ALLOCATE(GeomSet(MRESET*NCoord),GradSet(MRESET*NCoord),DX(NCoord),GSAVE(NCoord)) |
| 45 | 1 | equemene | END IF |
| 46 | 1 | equemene | END IF ! IF (FRST) THEN |
| 47 | 1 | equemene | |
| 48 | 1 | equemene | ! SPACE_GEDIIS SIMPLY LOADS THE CURRENT VALUES OF Geom AND Grad INTO THE ARRAYS GeomSet |
| 49 | 1 | equemene | ! AND GradSet, MSET is set to zero and then 1 in SPACE_GEDIIS at first iteration. |
| 50 | 1 | equemene | CALL SPACE_GEDIIS(MRESET,MSET,Geom,Grad,HEAT,NCoord,GeomSet,GradSet,ESET,FRST) |
| 51 | 1 | equemene | |
| 52 | 1 | equemene | IF (PRINT) WRITE(*,'(/,'' GEDIIS after SPACE_GEDIIS '')') |
| 53 | 1 | equemene | |
| 54 | 1 | equemene | ! INITIALIZE SOME VARIABLES AND CONSTANTS: |
| 55 | 1 | equemene | NGEDIIS = MSET !MSET=mth iteration |
| 56 | 1 | equemene | MPLUS = MSET + 1 |
| 57 | 1 | equemene | MM = MPLUS * MPLUS |
| 58 | 1 | equemene | |
| 59 | 1 | equemene | ! CONSTRUCT THE GEDIIS MATRIX: |
| 60 | 1 | equemene | ! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)> |
| 61 | 1 | equemene | JJ=0 |
| 62 | 1 | equemene | INV=-NCoord |
| 63 | 1 | equemene | DO I=1,MSET |
| 64 | 1 | equemene | INV=INV+NCoord |
| 65 | 1 | equemene | JNV=-NCoord |
| 66 | 1 | equemene | DO J=1,MSET |
| 67 | 1 | equemene | JNV=JNV+NCoord |
| 68 | 1 | equemene | JJ = JJ + 1 |
| 69 | 1 | equemene | B(JJ)=0.D0 |
| 70 | 1 | equemene | DO K=1, NCoord |
| 71 | 1 | equemene | B(JJ) = B(JJ) + (((GradSet(INV+K)-GradSet(JNV+K))*(GeomSet(INV+K)-GeomSet(JNV+K)))/2.D0) |
| 72 | 1 | equemene | END DO |
| 73 | 1 | equemene | END DO |
| 74 | 1 | equemene | END DO |
| 75 | 1 | equemene | |
| 76 | 1 | equemene | ! The following shifting is required to correct indices of B_ij elements in the GEDIIS matrix. |
| 77 | 1 | equemene | ! The correction is needed because the last coloumn of the matrix contains all 1 and one zero. |
| 78 | 1 | equemene | DO I=MSET-1,1,-1 |
| 79 | 1 | equemene | DO J=MSET,1,-1 |
| 80 | 1 | equemene | B(I*MSET+J+I) = B(I*MSET+J) |
| 81 | 1 | equemene | END DO |
| 82 | 1 | equemene | END DO |
| 83 | 1 | equemene | |
| 84 | 1 | equemene | ! for last row and last column of GEDIIS matrix |
| 85 | 1 | equemene | DO I=1,MPLUS |
| 86 | 1 | equemene | B(MPLUS*I) = 1.D0 |
| 87 | 1 | equemene | B(MPLUS*MSET+I) = 1.D0 |
| 88 | 1 | equemene | END DO |
| 89 | 1 | equemene | B(MM) = 0.D0 |
| 90 | 1 | equemene | |
| 91 | 1 | equemene | DO I=1, MPLUS |
| 92 | 1 | equemene | !WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS)) |
| 93 | 1 | equemene | END DO |
| 94 | 1 | equemene | |
| 95 | 1 | equemene | ! ELIMINATE ERROR VECTORS WITH THE LARGEST NORM: |
| 96 | 1 | equemene | 80 CONTINUE |
| 97 | 1 | equemene | DO I=1,MM !MM = (MSET+1) * (MSET+1) |
| 98 | 1 | equemene | BS(I) = B(I) !just a copy of the original B (GEDIIS) matrix |
| 99 | 1 | equemene | END DO |
| 100 | 1 | equemene | |
| 101 | 1 | equemene | IF (NGEDIIS .NE. MSET) THEN |
| 102 | 1 | equemene | DO II=1,MSET-NGEDIIS |
| 103 | 1 | equemene | XMAX = -1.D10 |
| 104 | 1 | equemene | ITERA = 0 |
| 105 | 1 | equemene | DO I=1,MSET |
| 106 | 1 | equemene | XNORM = 0.D0 |
| 107 | 1 | equemene | INV = (I-1) * MPLUS |
| 108 | 1 | equemene | DO J=1,MSET |
| 109 | 1 | equemene | XNORM = XNORM + ABS(B(INV + J)) |
| 110 | 1 | equemene | END DO |
| 111 | 1 | equemene | IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN |
| 112 | 1 | equemene | XMAX = XNORM |
| 113 | 1 | equemene | ITERA = I |
| 114 | 1 | equemene | IONE = INV + I |
| 115 | 1 | equemene | ENDIF |
| 116 | 1 | equemene | END DO |
| 117 | 1 | equemene | |
| 118 | 1 | equemene | DO I=1,MPLUS |
| 119 | 1 | equemene | INV = (I-1) * MPLUS |
| 120 | 1 | equemene | DO J=1,MPLUS |
| 121 | 1 | equemene | JNV = (J-1) * MPLUS |
| 122 | 1 | equemene | IF (J.EQ.ITERA) B(INV + J) = 0.D0 |
| 123 | 1 | equemene | B(JNV + I) = B(INV + J) |
| 124 | 1 | equemene | END DO |
| 125 | 1 | equemene | END DO |
| 126 | 1 | equemene | B(IONE) = 1.0D0 |
| 127 | 1 | equemene | END DO |
| 128 | 1 | equemene | END IF ! matches IF (NGEDIIS .NE. MSET) THEN |
| 129 | 1 | equemene | |
| 130 | 1 | equemene | ! SCALE GEDIIS MATRIX BEFORE INVERSION: |
| 131 | 1 | equemene | DO I=1,MPLUS |
| 132 | 1 | equemene | II = MPLUS * (I-1) + I ! B(II)=diagonal elements of B matrix |
| 133 | 1 | equemene | GSAVE(I) = 1.D0 / DSQRT(1.D-20+DABS(B(II))) |
| 134 | 1 | equemene | !Print *, 'GSAVE(',I,')=', GSAVE(I)
|
| 135 | 1 | equemene | END DO |
| 136 | 1 | equemene | GSAVE(MPLUS) = 1.D0 |
| 137 | 1 | equemene | DO I=1,MPLUS |
| 138 | 1 | equemene | DO J=1,MPLUS |
| 139 | 1 | equemene | IJ = MPLUS * (I-1) + J |
| 140 | 1 | equemene | B(IJ) = B(IJ) * GSAVE(I) * GSAVE(J) |
| 141 | 1 | equemene | END DO |
| 142 | 1 | equemene | END DO |
| 143 | 1 | equemene | |
| 144 | 1 | equemene | ! INVERT THE GEDIIS MATRIX B: |
| 145 | 1 | equemene | DO I=1, MPLUS |
| 146 | 1 | equemene | !WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS)) |
| 147 | 1 | equemene | END DO |
| 148 | 1 | equemene | |
| 149 | 1 | equemene | CALL MINV(B,MPLUS,DET) ! matrix inversion. |
| 150 | 1 | equemene | |
| 151 | 1 | equemene | DO I=1, MPLUS |
| 152 | 1 | equemene | !WRITE(*,'(10(1X,F20.16))') B((I-1)*MPLUS+1:I*(MPLUS)) |
| 153 | 1 | equemene | END DO |
| 154 | 1 | equemene | |
| 155 | 1 | equemene | DO I=1,MPLUS |
| 156 | 1 | equemene | DO J=1,MPLUS |
| 157 | 1 | equemene | IJ = MPLUS * (I-1) + J |
| 158 | 1 | equemene | B(IJ) = B(IJ) * GSAVE(I) * GSAVE(J) |
| 159 | 1 | equemene | END DO |
| 160 | 1 | equemene | END DO |
| 161 | 1 | equemene | |
| 162 | 1 | equemene | ! COMPUTE THE NEW INTERPOLATED PARAMETER VECTOR (Geometry): |
| 163 | 1 | equemene | ci=0.d0 |
| 164 | 1 | equemene | ci_tmp=0.d0 |
| 165 | 1 | equemene | |
| 166 | 1 | equemene | ci_lt_zero= .FALSE. |
| 167 | 1 | equemene | DO I=1, MSET |
| 168 | 1 | equemene | DO J=1, MSET ! B matrix is read column-wise |
| 169 | 1 | equemene | ci(I)=ci(I)+B((J-1)*(MPLUS)+I)*ESET(J) !ESET is energy set, yet to be fixed. |
| 170 | 1 | equemene | END DO |
| 171 | 1 | equemene | ci(I)=ci(I)+B((MPLUS-1)*(MPLUS)+I) |
| 172 | 1 | equemene | !Print *, 'NO ci < 0 yet, c(',I,')=', ci(I)
|
| 173 | 1 | equemene | IF((ci(I) .LT. 0.0D0) .OR. (ci(I) .GT. 1.0D0)) THEN |
| 174 | 1 | equemene | ci_lt_zero=.TRUE. |
| 175 | 1 | equemene | EXIT |
| 176 | 1 | equemene | END IF |
| 177 | 1 | equemene | END DO !matches DO I=1, MSET |
| 178 | 1 | equemene | |
| 179 | 1 | equemene | IF (ci_lt_zero) Then |
| 180 | 1 | equemene | cis_zero = 0 |
| 181 | 1 | equemene | ER_star = 0.D0 |
| 182 | 1 | equemene | ER_star_tmp = 1e32 |
| 183 | 1 | equemene | |
| 184 | 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. |
| 185 | 1 | equemene | JJ=0 |
| 186 | 1 | equemene | INV=-NCoord |
| 187 | 1 | equemene | DO IX=1,MSET |
| 188 | 1 | equemene | INV=INV+NCoord |
| 189 | 1 | equemene | JNV=-NCoord |
| 190 | 1 | equemene | DO JX=1,MSET |
| 191 | 1 | equemene | JNV=JNV+NCoord |
| 192 | 1 | equemene | JJ = JJ + 1 |
| 193 | 1 | equemene | BST(JJ)=0.D0 |
| 194 | 1 | equemene | DO KX=1, NCoord |
| 195 | 1 | equemene | BST(JJ) = BST(JJ) + (((GradSet(INV+KX)-GradSet(JNV+KX))*(GeomSet(INV+KX)-GeomSet(JNV+KX)))/2.D0) |
| 196 | 1 | equemene | END DO |
| 197 | 1 | equemene | END DO |
| 198 | 1 | equemene | END DO |
| 199 | 1 | equemene | |
| 200 | 1 | equemene | DO I=1, (2**MSET)-2 ! all (2**MSET)-2 combinations of cis, except the one where all cis are .GT. 0 and .LT. 1 |
| 201 | 1 | equemene | ci(:)=1.D0 |
| 202 | 1 | equemene | ! find out which cis are zero in I: |
| 203 | 1 | equemene | DO IX=1, MSET |
| 204 | 1 | equemene | JJ=IAND(I, 2**(IX-1)) |
| 205 | 1 | equemene | IF(JJ .EQ. 0) Then |
| 206 | 1 | equemene | ci(IX)=0.D0 |
| 207 | 1 | equemene | END IF |
| 208 | 1 | equemene | END DO |
| 209 | 1 | equemene | |
| 210 | 1 | equemene | ci_lt_zero = .FALSE. |
| 211 | 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. |
| 212 | 1 | equemene | DO IX=1, MSET*MSET |
| 213 | 1 | equemene | B(IX) = BST(IX) !just a copy of the original B (GEDIIS) matrix |
| 214 | 1 | equemene | END DO |
| 215 | 1 | equemene | |
| 216 | 1 | equemene | ! Removal of KXth row and KXth column in order to accomodate cI to be zero: |
| 217 | 1 | equemene | current_size_B_mat=MSET |
| 218 | 1 | equemene | cis_zero = 0 |
| 219 | 1 | equemene | ! The bits of I (index of the upper loop 'DO I=1, (2**MSET)-2'), gives which cis are zero. |
| 220 | 1 | equemene | DO KX=1, MSET ! searching for each bit of I (index of the upper loop 'DO I=1, (2**MSET)-2') |
| 221 | 1 | equemene | IF (ci(KX) .EQ. 0.D0) Then !remove KXth row and KXth column |
| 222 | 1 | equemene | cis_zero = cis_zero + 1 |
| 223 | 1 | equemene | |
| 224 | 1 | equemene | ! First row removal: (B matrix is read column-wise) |
| 225 | 1 | equemene | JJ=0 |
| 226 | 1 | equemene | DO IX=1,current_size_B_mat ! columns reading |
| 227 | 1 | equemene | DO JX=1,current_size_B_mat ! rows reading |
| 228 | 1 | equemene | IF (JX .NE. KX) Then |
| 229 | 1 | equemene | JJ = JJ + 1 |
| 230 | 1 | equemene | B_tmp(JJ) = B((IX-1)*current_size_B_mat+JX) |
| 231 | 1 | equemene | END IF |
| 232 | 1 | equemene | END DO |
| 233 | 1 | equemene | END DO |
| 234 | 1 | equemene | |
| 235 | 1 | equemene | DO IX=1,current_size_B_mat*(current_size_B_mat-1) |
| 236 | 1 | equemene | B(IX) = B_tmp(IX) |
| 237 | 1 | equemene | END DO |
| 238 | 1 | equemene | |
| 239 | 1 | equemene | ! Now column removal: |
| 240 | 1 | equemene | JJ=0 |
| 241 | 1 | equemene | DO IX=1,KX-1 ! columns reading |
| 242 | 1 | equemene | DO JX=1,current_size_B_mat-1 ! rows reading |
| 243 | 1 | equemene | JJ = JJ + 1 |
| 244 | 1 | equemene | B_tmp(JJ) = B(JJ) |
| 245 | 1 | equemene | END DO |
| 246 | 1 | equemene | END DO |
| 247 | 1 | equemene | |
| 248 | 1 | equemene | DO IX=KX+1,current_size_B_mat |
| 249 | 1 | equemene | DO JX=1,current_size_B_mat-1 |
| 250 | 1 | equemene | JJ = JJ + 1 |
| 251 | 1 | equemene | B_tmp(JJ) = B(JJ+current_size_B_mat-1) |
| 252 | 1 | equemene | END DO |
| 253 | 1 | equemene | END DO |
| 254 | 1 | equemene | |
| 255 | 1 | equemene | DO IX=1,(current_size_B_mat-1)*(current_size_B_mat-1) |
| 256 | 1 | equemene | B(IX) = B_tmp(IX) |
| 257 | 1 | equemene | END DO |
| 258 | 1 | equemene | current_size_B_mat = current_size_B_mat - 1 |
| 259 | 1 | equemene | END IF ! matches IF (ci(KX) .EQ. 0.D0) Then !remove |
| 260 | 1 | equemene | END DO ! matches DO KX=1, MSET |
| 261 | 1 | equemene | |
| 262 | 1 | equemene | ! The following shifting is required to correct indices of B_ij elements in the GEDIIS matrix. |
| 263 | 1 | equemene | ! The correction is needed because the last coloumn and row of the matrix contains all 1 and one zero. |
| 264 | 1 | equemene | DO IX=MSET-cis_zero-1,1,-1 |
| 265 | 1 | equemene | DO JX=MSET-cis_zero,1,-1 |
| 266 | 1 | equemene | B(IX*(MSET-cis_zero)+JX+IX) = B(IX*(MSET-cis_zero)+JX) |
| 267 | 1 | equemene | END DO |
| 268 | 1 | equemene | END DO |
| 269 | 1 | equemene | |
| 270 | 1 | equemene | ! for last row and last column of GEDIIS matrix |
| 271 | 1 | equemene | DO IX=1,MPLUS-cis_zero |
| 272 | 1 | equemene | B((MPLUS-cis_zero)*IX) = 1.D0 |
| 273 | 1 | equemene | B((MPLUS-cis_zero)*(MSET-cis_zero)+IX) = 1.D0 |
| 274 | 1 | equemene | END DO |
| 275 | 1 | equemene | B((MPLUS-cis_zero) * (MPLUS-cis_zero)) = 0.D0 |
| 276 | 1 | equemene | |
| 277 | 1 | equemene | DO IX=1, MPLUS |
| 278 | 1 | equemene | !WRITE(*,'(10(1X,F20.4))') B((IX-1)*MPLUS+1:IX*(MPLUS)) |
| 279 | 1 | equemene | END DO |
| 280 | 1 | equemene | |
| 281 | 1 | equemene | ! ELIMINATE ERROR VECTORS WITH THE LARGEST NORM: |
| 282 | 1 | equemene | IF (NGEDIIS .NE. MSET) THEN |
| 283 | 1 | equemene | JX = min(MSET-NGEDIIS,MSET-cis_zero-1) |
| 284 | 1 | equemene | DO II=1,JX |
| 285 | 1 | equemene | XMAX = -1.D10 |
| 286 | 1 | equemene | ITERA = 0 |
| 287 | 1 | equemene | DO IX=1,MSET-cis_zero |
| 288 | 1 | equemene | XNORM = 0.D0 |
| 289 | 1 | equemene | INV = (IX-1) * (MPLUS-cis_zero) |
| 290 | 1 | equemene | DO J=1,MSET-cis_zero |
| 291 | 1 | equemene | XNORM = XNORM + ABS(B(INV + J)) |
| 292 | 1 | equemene | END DO |
| 293 | 1 | equemene | IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN |
| 294 | 1 | equemene | XMAX = XNORM |
| 295 | 1 | equemene | ITERA = IX |
| 296 | 1 | equemene | IONE = INV + IX |
| 297 | 1 | equemene | ENDIF |
| 298 | 1 | equemene | END DO |
| 299 | 1 | equemene | |
| 300 | 1 | equemene | DO IX=1,MPLUS-cis_zero |
| 301 | 1 | equemene | INV = (IX-1) * (MPLUS-cis_zero) |
| 302 | 1 | equemene | DO J=1,MPLUS-cis_zero |
| 303 | 1 | equemene | JNV = (J-1) * (MPLUS-cis_zero) |
| 304 | 1 | equemene | IF (J.EQ.ITERA) B(INV + J) = 0.D0 |
| 305 | 1 | equemene | B(JNV + IX) = B(INV + J) |
| 306 | 1 | equemene | END DO |
| 307 | 1 | equemene | END DO |
| 308 | 1 | equemene | B(IONE) = 1.0D0 |
| 309 | 1 | equemene | END DO |
| 310 | 1 | equemene | END IF ! matches IF (NGEDIIS .NE. MSET) THEN |
| 311 | 1 | equemene | |
| 312 | 1 | equemene | ! SCALE GEDIIS MATRIX BEFORE INVERSION: |
| 313 | 1 | equemene | DO IX=1,MPLUS-cis_zero |
| 314 | 1 | equemene | II = (MPLUS-cis_zero) * (IX-1) + IX ! B(II)=diagonal elements of B matrix |
| 315 | 1 | equemene | GSAVE(IX) = 1.D0 / DSQRT(1.D-20+DABS(B(II))) |
| 316 | 1 | equemene | END DO |
| 317 | 1 | equemene | GSAVE(MPLUS-cis_zero) = 1.D0 |
| 318 | 1 | equemene | DO IX=1,MPLUS-cis_zero |
| 319 | 1 | equemene | DO JX=1,MPLUS-cis_zero |
| 320 | 1 | equemene | IJ = (MPLUS-cis_zero) * (IX-1) + JX |
| 321 | 1 | equemene | B(IJ) = B(IJ) * GSAVE(IX) * GSAVE(JX) |
| 322 | 1 | equemene | END DO |
| 323 | 1 | equemene | END DO |
| 324 | 1 | equemene | |
| 325 | 1 | equemene | ! INVERT THE GEDIIS MATRIX B: |
| 326 | 1 | equemene | CALL MINV(B,MPLUS-cis_zero,DET) ! matrix inversion. |
| 327 | 1 | equemene | |
| 328 | 1 | equemene | DO IX=1,MPLUS-cis_zero |
| 329 | 1 | equemene | DO JX=1,MPLUS-cis_zero |
| 330 | 1 | equemene | IJ = (MPLUS-cis_zero) * (IX-1) + JX |
| 331 | 1 | equemene | B(IJ) = B(IJ) * GSAVE(IX) * GSAVE(JX) |
| 332 | 1 | equemene | END DO |
| 333 | 1 | equemene | END DO |
| 334 | 1 | equemene | |
| 335 | 1 | equemene | DO IX=1, MPLUS |
| 336 | 1 | equemene | !WRITE(*,'(10(1X,F20.4))') B((IX-1)*MPLUS+1:IX*(MPLUS)) |
| 337 | 1 | equemene | END DO |
| 338 | 1 | equemene | |
| 339 | 1 | equemene | ! ESET is rearranged to handle zero cis and stored in ESET_tmp: |
| 340 | 1 | equemene | JJ=0 |
| 341 | 1 | equemene | DO IX=1, MSET |
| 342 | 1 | equemene | IF (ci(IX) .NE. 0) Then |
| 343 | 1 | equemene | JJ=JJ+1 |
| 344 | 1 | equemene | ESET_tmp(JJ) = ESET(IX) |
| 345 | 1 | equemene | END IF |
| 346 | 1 | equemene | END DO |
| 347 | 1 | equemene | |
| 348 | 1 | equemene | ! DETERMINATION OF nonzero cis: |
| 349 | 1 | equemene | MyPointer=1 |
| 350 | 1 | equemene | DO IX=1, MSET-cis_zero |
| 351 | 1 | equemene | tmp = 0.D0 |
| 352 | 1 | equemene | DO J=1, MSET-cis_zero ! B matrix is read column-wise |
| 353 | 1 | equemene | tmp=tmp+B((J-1)*(MPLUS-cis_zero)+IX)*ESET_tmp(J) |
| 354 | 1 | equemene | END DO |
| 355 | 1 | equemene | tmp=tmp+B((MPLUS-cis_zero-1)*(MPLUS-cis_zero)+IX) |
| 356 | 1 | equemene | IF((tmp .LT. 0.0D0) .OR. (tmp .GT. 1.0D0)) THEN |
| 357 | 1 | equemene | ci_lt_zero=.TRUE. |
| 358 | 1 | equemene | EXIT |
| 359 | 1 | equemene | ELSE |
| 360 | 1 | equemene | DO JX=MyPointer,MSET |
| 361 | 1 | equemene | IF (ci(JX) .NE. 0) Then |
| 362 | 1 | equemene | ci(JX) = tmp |
| 363 | 1 | equemene | MyPointer=JX+1 |
| 364 | 1 | equemene | EXIT |
| 365 | 1 | equemene | END IF |
| 366 | 1 | equemene | END DO |
| 367 | 1 | equemene | END IF |
| 368 | 1 | equemene | END DO !matches DO I=1, MSET-cis_zero |
| 369 | 1 | equemene | !Print *, 'Local set of cis, first 10:, MSET=', MSET, ', I of (2**MSET)-2=', I |
| 370 | 1 | equemene | !WRITE(*,'(10(1X,F20.4))') ci(1:MSET) |
| 371 | 1 | equemene | !Print *, 'Local set of cis ends:****************************************' |
| 372 | 1 | equemene | |
| 373 | 1 | equemene | ! new set of cis determined based on the lower energy (ER_star): |
| 374 | 1 | equemene | IF (.NOT. ci_lt_zero) Then |
| 375 | 1 | equemene | Call Energy_GEDIIS(MRESET,MSET,ci,GeomSet,GradSet,ESET,NCoord,ER_star) |
| 376 | 1 | equemene | IF (ER_star .LT. ER_star_tmp) Then |
| 377 | 1 | equemene | ci_tmp=ci |
| 378 | 1 | equemene | ER_star_tmp = ER_star |
| 379 | 1 | equemene | END IF |
| 380 | 1 | equemene | END IF ! matches IF (.NOT. ci_lt_zero) Then |
| 381 | 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 |
| 382 | 1 | equemene | ci = ci_tmp |
| 383 | 1 | equemene | END IF! matches IF (ci_lt_zero) Then |
| 384 | 1 | equemene | |
| 385 | 1 | equemene | Print *, 'Final set of cis, first 10:***********************************' |
| 386 | 1 | equemene | WRITE(*,'(10(1X,F20.4))') ci(1:MSET) |
| 387 | 1 | equemene | Print *, 'Final set of cis ends:****************************************' |
| 388 | 1 | equemene | Geom_new(:) = 0.D0 |
| 389 | 1 | equemene | DO I=1, MSET |
| 390 | 1 | equemene | Geom_new(:) = Geom_new(:) + (ci(I)*GeomSet((I-1)*NCoord+1:I*NCoord)) !MPLUS=MSET+1 |
| 391 | 1 | equemene | ! R_(N+1)=R*+DeltaR: |
| 392 | 1 | equemene | DO J=1, NCoord |
| 393 | 1 | equemene | tmp=0.D0 |
| 394 | 1 | equemene | DO K=1,NCoord |
| 395 | 1 | equemene | !tmp=tmp+Hess((J-1)*NCoord+K)*GradSet((I-1)*NCoord+K) ! If Hinv=.False., then we need to invert Hess |
| 396 | 1 | equemene | END DO |
| 397 | 1 | equemene | Geom_new(J) = Geom_new(J) - (ci(I)*tmp) |
| 398 | 1 | equemene | END DO |
| 399 | 1 | equemene | END DO |
| 400 | 1 | equemene | |
| 401 | 1 | equemene | DX(:) = Geom(:) - Geom_new(:) |
| 402 | 1 | equemene | |
| 403 | 1 | equemene | XNORM = SQRT(DOT_PRODUCT(DX,DX)) |
| 404 | 1 | equemene | IF (PRINT) THEN |
| 405 | 1 | equemene | WRITE (6,'(/10X,''DEVIATION IN X '',F10.4,8X,''DETERMINANT '',G9.3)') XNORM, DET |
| 406 | 1 | equemene | !WRITE(*,'(10X,''GEDIIS COEFFICIENTS'')') |
| 407 | 1 | equemene | !WRITE(*,'(10X,5F12.5)') (B(MPLUS*MSET+I),I=1,MSET) |
| 408 | 1 | equemene | ENDIF |
| 409 | 1 | equemene | |
| 410 | 1 | equemene | ! THE FOLLOWING TOLERENCES FOR XNORM AND DET ARE SOMEWHAT ARBITRARY! |
| 411 | 1 | equemene | THRES = MAX(10.D0**(-NCoord), 1.D-25) |
| 412 | 1 | equemene | IF (XNORM.GT.2.D0 .OR. DABS(DET) .LT. THRES) THEN |
| 413 | 1 | equemene | IF (PRINT)THEN |
| 414 | 1 | equemene | WRITE(*,*) "THE GEDIIS MATRIX IS ILL CONDITIONED" |
| 415 | 1 | equemene | WRITE(*,*) " - PROBABLY, VECTORS ARE LINEARLY DEPENDENT - " |
| 416 | 1 | equemene | WRITE(*,*) "THE GEDIIS STEP WILL BE REPEATED WITH A SMALLER SPACE" |
| 417 | 1 | equemene | END IF |
| 418 | 1 | equemene | DO K=1,MM |
| 419 | 1 | equemene | B(K) = BS(K) ! why this is reverted? Because "IF (NGEDIIS .GT. 0) GO TO 80", see below |
| 420 | 1 | equemene | END DO |
| 421 | 1 | equemene | NGEDIIS = NGEDIIS - 1 |
| 422 | 1 | equemene | IF (NGEDIIS .GT. 0) GO TO 80 |
| 423 | 1 | equemene | IF (PRINT) WRITE(*,'(10X,''NEWTON-RAPHSON STEP TAKEN'')') |
| 424 | 1 | equemene | Geom_new(:) = Geom(:) ! Geom_new is set to original Geom, thus DX = Geom(:) - Geom_new(:)=zero |
| 425 | 1 | equemene | END IF ! matches IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN |
| 426 | 1 | equemene | |
| 427 | 1 | equemene | !******************************************************************************************************************* |
| 428 | 1 | equemene | Geom_new(:) = 0.D0 |
| 429 | 1 | equemene | DO I=1, MSET |
| 430 | 1 | equemene | Geom_new(:) = Geom_new(:) + (ci(I)*GeomSet((I-1)*NCoord+1:I*NCoord)) !MPLUS=MSET+1 |
| 431 | 1 | equemene | ! R_(N+1)=R*+DeltaR: |
| 432 | 1 | equemene | DO J=1, NCoord |
| 433 | 1 | equemene | tmp=0.D0 |
| 434 | 1 | equemene | DO K=1,NCoord |
| 435 | 1 | equemene | tmp=tmp+Hess((J-1)*NCoord+K)*GradSet((I-1)*NCoord+K) ! If Hinv=.False., then we need to invert Hess |
| 436 | 1 | equemene | END DO |
| 437 | 1 | equemene | Geom_new(J) = Geom_new(J) - (ci(I)*tmp) |
| 438 | 1 | equemene | END DO |
| 439 | 1 | equemene | END DO |
| 440 | 1 | equemene | !******************************************************************************************************************* |
| 441 | 1 | equemene | |
| 442 | 1 | equemene | IF (PRINT) WRITE(*,'(/,'' END GEDIIS '',/)') |
| 443 | 1 | equemene | |
| 444 | 1 | equemene | END SUBROUTINE Step_GEDIIS |