root / src / Step_GEDIIS_All.f90 @ 8
History  View  Annotate  Download (21.6 kB)
1 
! Geom = input parameter vector (Geometry), Grad = input gradient vector, HEAT is Energy(Geom) 

2 
SUBROUTINE Step_GEDIIS_All(NGeomF,IGeom,Step,Geom,Grad,HEAT,Hess,NCoord,allocation_flag,Tangent) 
3 
!SUBROUTINE Step_GEDIIS(Geom_new,Geom,Grad,HEAT,Hess,NCoord,FRST) 
4 
use Io_module 
5 
use Path_module, only : Nom, Atome, OrderInv, indzmat, Pi, Nat, Vfree 
6 
IMPLICIT NONE 
7  
8 
INTEGER(KINT) :: NGeomF,IGeom 
9 
INTEGER(KINT), INTENT(IN) :: NCoord 
10 
REAL(KREAL) :: Geom(NCoord), Grad(NCoord), Hess(NCoord*NCoord), Step(NCoord) 
11 
REAL(KREAL) :: HEAT ! HEAT= Energy 
12 
LOGICAL :: allocation_flag 
13 
REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord) 
14 

15 
! MRESET = maximum number of iterations. 
16 
INTEGER(KINT), PARAMETER :: MRESET=15, M2=(MRESET+1)*(MRESET+1) !M2 = 256 
17 
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:,:), GradSet(:,:) ! NGeomF,MRESET*NCoord 
18 
REAL(KREAL), ALLOCATABLE, SAVE :: GSAVE(:,:) !NGeomF,NCoord 
19 
REAL(KREAL), ALLOCATABLE, SAVE :: ESET(:,:) 
20 
REAL(KREAL) :: ESET_tmp(MRESET), B(M2), BS(M2), BST(M2), B_tmp(M2) ! M2=256 
21 
LOGICAL :: DEBUG, PRINT, ci_lt_zero 
22 
INTEGER(KINT), ALLOCATABLE, SAVE :: MSET(:) ! mth Iteration 
23 
LOGICAL, ALLOCATABLE, SAVE :: FRST(:) 
24 
REAL(KREAL) :: ci(MRESET), ci_tmp(MRESET) ! MRESET = maximum number of iterations. 
25 
INTEGER(KINT) :: NGEDIIS, MPLUS, INV, ITERA, MM, cis_zero 
26 
INTEGER(KINT) :: I, J, K, JJ, JNV, II, IONE, IJ, IX, JX, KX 
27 
INTEGER(KINT) :: current_size_B_mat, MyPointer, Isch, NFree, Idx 
28 
REAL(KREAL) :: XMax, XNorm, DET, THRES, tmp, ER_star, ER_star_tmp, Norm 
29 
REAL(KREAL), PARAMETER :: eps=1e12 
30 
REAL(KREAL), PARAMETER :: crit=1e8 
31 
REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord 
32 
REAL(KREAL), ALLOCATABLE :: HFree(:) ! NFree*NFree 
33 
REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree 
34 
REAL(KREAL), ALLOCATABLE :: Grad_free(:), Step_free(:) ! NFree 
35 
REAL(KREAL), ALLOCATABLE :: Geom_free(:), Geom_new_free(:) ! NFree 
36 
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet_free(:,:), GradSet_free(:,:) 
37  
38 
DEBUG=.TRUE. 
39 
PRINT=.FALSE. 
40 

41 
IF (PRINT) WRITE(*,'(/,'' BEGIN Step_GEDIIS_ALL '')') 
42 

43 
! Initialization 
44 
IF (allocation_flag) THEN 
45 
! allocation_flag will be set to False in SPACE_GEDIIS, so no need to modify it here 
46 
IF (ALLOCATED(GeomSet)) THEN 
47 
IF (PRINT) WRITE(*,'(/,'' In allocation_flag, GEDIIS_ALL Dealloc '')') 
48 
DEALLOCATE(GeomSet,GradSet,GSave,GeomSet_free,GradSet_free) 
49 
RETURN 
50 
ELSE 
51 
IF (PRINT) WRITE(*,'(/,'' In allocation_flag, GEDIIS_ALL Alloc '')') 
52 
ALLOCATE(GeomSet(NGeomF,MRESET*NCoord),GradSet(NGeomF,MRESET*NCoord),GSAVE(NGeomF,NCoord)) 
53 
ALLOCATE(GeomSet_free(NGeomF,MRESET*NCoord),GradSet_free(NGeomF,MRESET*NCoord)) 
54 
ALLOCATE(MSET(NGeomF),FRST(NGeomF),ESET(NGeomF,MRESET)) 
55 
DO I=1,NGeomF 
56 
FRST(I) = .TRUE. 
57 
GeomSet(I,:) = 0.d0 
58 
GradSet(I,:) = 0.d0 
59 
GSAVE(I,:)=0.d0 
60 
GeomSet_free(I,:) = 0.d0 
61 
GradSet_free(I,:) = 0.d0 
62 
END DO 
63 
MSET(:)=0 
64 
END IF 
65 
allocation_flag = .FALSE. 
66 
END IF ! IF (allocation_flag) THEN 
67 

68 
! ADDED FROM HERE: 
69 
Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord), as of now, an Identity matrix. 
70 
! we orthogonalize Vfree to the tangent vector of this geom only if Tangent/=0.d0 
71 
Norm=sqrt(dot_product(Tangent,Tangent)) 
72 
IF (Norm.GT.eps) THEN 
73 
ALLOCATE(Tanf(NCoord)) 
74  
75 
! We normalize Tangent 
76 
Tangent=Tangent/Norm 
77  
78 
! We convert Tangent into Vfree only displacements. This is useless for now (2007.Apr.23) 
79 
! as Vfree=Id matrix but it will be usefull as soon as we introduce constraints. 
80 
DO I=1,NCoord 
81 
Tanf(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Tangent) 
82 
END DO 
83 
Tangent=0.d0 
84 
DO I=1,NCoord 
85 
Tangent=Tangent+Tanf(I)*Vfree(:,I) 
86 
END DO 
87 
! first we subtract Tangent from vfree 
88 
DO I=1,NCoord 
89 
Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent) 
90 
Vfree(:,I)=Vfree(:,I)Norm*Tangent 
91 
END DO 
92  
93 
Idx=0 
94 
! Schmidt orthogonalization of the Vfree vectors 
95 
DO I=1,NCoord 
96 
! We subtract the first vectors, we do it twice as the Schmidt procedure is not numerically stable. 
97 
DO Isch=1,2 
98 
DO J=1,Idx 
99 
Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,J),(/NCoord/))) 
100 
Vfree(:,I)=Vfree(:,I)Norm*Vfree(:,J) 
101 
END DO 
102 
END DO 
103 
Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,I),(/NCoord/))) 
104 
IF (Norm.GE.crit) THEN 
105 
Idx=Idx+1 
106 
Vfree(:,Idx)=Vfree(:,I)/sqrt(Norm) 
107 
END IF 
108 
END DO 
109 

110 
IF (Idx/= NCoord1) THEN 
111 
WRITE(*,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom 
112 
WRITE(IOOut,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom 
113 
STOP 
114 
END IF 
115 

116 
DEALLOCATE(Tanf) 
117 
NFree=Idx 
118 
ELSE ! Tangent =0, matches IF (Norm.GT.eps) THEN 
119 
if (debug) WRITE(*,*) "Tangent=0, using full displacement" 
120 
NFree=NCoord 
121 
END IF !IF (Norm.GT.eps) THEN 
122 

123 
if (debug) WRITE(*,*) 'DBG Step_GEDIIS_All, IGeom, NFree=', IGeom, NFree 
124  
125 
! We now calculate the new step 
126 
! we project the hessian onto the free vectors 
127 
ALLOCATE(HFree(NFree*NFree),Htmp(NCoord,NFree),Grad_free(NFree)) 
128 
ALLOCATE(Geom_free(NFree),Step_free(NFree),Geom_new_free(NFree)) 
129 
DO J=1,NFree 
130 
DO I=1,NCoord 
131 
Htmp(I,J)=0.d0 
132 
DO K=1,NCoord 
133 
Htmp(I,J)=Htmp(I,J)+Hess(((I1)*NCoord)+K)*Vfree(K,J) 
134 
END DO 
135 
END DO 
136 
END DO 
137 
DO J=1,NFree 
138 
DO I=1,NFree 
139 
HFree(I+((J1)*NFree))=0.d0 
140 
DO K=1,NCoord 
141 
HFree(I+((J1)*NFree))=HFree(I+((J1)*NFree))+Vfree(K,I)*Htmp(K,J) 
142 
END DO 
143 
END DO 
144 
END DO 
145  
146 
DO I=1,NFree 
147 
Grad_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Grad) 
148 
Geom_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Geom) 
149 
END DO 
150 
!ADDED ENDS HERE.*********************************************** 
151 

152 
! SPACE_GEDIIS SIMPLY LOADS THE CURRENT VALUES OF Geom AND Grad INTO THE ARRAYS GeomSet 
153 
! AND GradSet, MSET is set to zero and then 1 in SPACE_GEDIIS_All at first iteration. 
154 
CALL SPACE_GEDIIS_All(NGeomF,IGeom,MRESET,MSET,Geom,Grad,HEAT,NCoord,GeomSet,GradSet,ESET,FRST) 
155  
156 
IF (PRINT) WRITE(*,'(/,'' GEDIIS after SPACE_GEDIIS_ALL '')') 
157 

158 
DO J=1,MSet(IGeom) 
159 
DO K=1,NFree 
160 
GradSet_free(IGeom,((J1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),& 
161 
GradSet(IGeom,((J1)*NCoord)+1:((J1)*NCoord)+NCoord)) 
162 
GeomSet_free(IGeom,((J1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),& 
163 
GeomSet(IGeom,((J1)*NCoord)+1:((J1)*NCoord)+NCoord)) 
164 
END DO 
165 
END DO 
166  
167 
! INITIALIZE SOME VARIABLES AND CONSTANTS: 
168 
NGEDIIS = MSET(IGeom) !MSET=mth iteration 
169 
MPLUS = MSET(IGeom) + 1 
170 
MM = MPLUS * MPLUS 
171  
172 
! CONSTRUCT THE GEDIIS MATRIX: 
173 
! B_ij calculations from <B_ij=(g_ig_j)(R_iR_j)> 
174 
JJ=0 
175 
INV=NFree 
176 
DO I=1,MSET(IGeom) 
177 
INV=INV+NFree 
178 
JNV=NFree 
179 
DO J=1,MSET(IGeom) 
180 
JNV=JNV+NFree 
181 
JJ = JJ + 1 
182 
B(JJ)=0.D0 
183 
DO K=1, NFree 
184 
B(JJ) = B(JJ) + (((GradSet_free(IGeom,INV+K)GradSet_free(IGeom,JNV+K))* & 
185 
(GeomSet_free(IGeom,INV+K)GeomSet_free(IGeom,JNV+K)))/2.D0) 
186 
END DO 
187 
END DO 
188 
END DO 
189  
190 
! The following shifting is required to correct indices of B_ij elements in the GEDIIS matrix. 
191 
! The correction is needed because the last coloumn of the matrix contains all 1 and one zero. 
192 
DO I=MSET(IGeom)1,1,1 
193 
DO J=MSET(IGeom),1,1 
194 
B(I*MSET(IGeom)+J+I) = B(I*MSET(IGeom)+J) 
195 
END DO 
196 
END DO 
197 

198 
! For the last row and last column of GEDIIS matrix: 
199 
DO I=1,MPLUS 
200 
B(MPLUS*I) = 1.D0 
201 
B(MPLUS*MSET(IGeom)+I) = 1.D0 
202 
END DO 
203 
B(MM) = 0.D0 
204 

205 
DO I=1, MPLUS 
206 
!WRITE(*,'(10(1X,F20.4))') B((I1)*MPLUS+1:I*(MPLUS)) 
207 
END DO 
208 

209 
! ELIMINATE ERROR VECTORS WITH THE LARGEST NORM: 
210 
80 CONTINUE 
211 
DO I=1,MM !MM = (MSET(IGeom)+1) * (MSET(IGeom)+1) 
212 
BS(I) = B(I) !just a copy of the original B (GEDIIS) matrix 
213 
END DO 
214 

215 
IF (NGEDIIS .NE. MSET(IGeom)) THEN 
216 
DO II=1,MSET(IGeom)NGEDIIS 
217 
XMAX = 1.D10 
218 
ITERA = 0 
219 
DO I=1,MSET(IGeom) 
220 
XNORM = 0.D0 
221 
INV = (I1) * MPLUS 
222 
DO J=1,MSET(IGeom) 
223 
XNORM = XNORM + ABS(B(INV + J)) 
224 
END DO 
225 
IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN 
226 
XMAX = XNORM 
227 
ITERA = I 
228 
IONE = INV + I 
229 
ENDIF 
230 
END DO 
231 

232 
DO I=1,MPLUS 
233 
INV = (I1) * MPLUS 
234 
DO J=1,MPLUS 
235 
JNV = (J1) * MPLUS 
236 
IF (J.EQ.ITERA) B(INV + J) = 0.D0 
237 
B(JNV + I) = B(INV + J) 
238 
END DO 
239 
END DO 
240 
B(IONE) = 1.0D0 
241 
END DO 
242 
END IF ! matches IF (NGEDIIS .NE. MSET(IGeom)) THEN 
243 

244 
! SCALE GEDIIS MATRIX BEFORE INVERSION: 
245 
DO I=1,MPLUS 
246 
II = MPLUS * (I1) + I ! B(II)=diagonal elements of B matrix 
247 
GSAVE(IGeom,I) = 1.D0 / DSQRT(1.D20+DABS(B(II))) 
248 
!Print *, 'GSAVE(',IGeom,',',I,')=', GSAVE(IGeom,I) 
249 
END DO 
250 
GSAVE(IGeom,MPLUS) = 1.D0 
251 
DO I=1,MPLUS 
252 
DO J=1,MPLUS 
253 
IJ = MPLUS * (I1) + J 
254 
B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J) 
255 
END DO 
256 
END DO 
257 

258 
! INVERT THE GEDIIS MATRIX B: 
259 
DO I=1, MPLUS 
260 
!WRITE(*,'(10(1X,F20.4))') B((I1)*MPLUS+1:I*(MPLUS)) 
261 
END DO 
262 

263 
CALL MINV(B,MPLUS,DET) ! matrix inversion. 
264 

265 
DO I=1, MPLUS 
266 
!WRITE(*,'(10(1X,F20.16))') B((I1)*MPLUS+1:I*(MPLUS)) 
267 
END DO 
268 

269 
DO I=1,MPLUS 
270 
DO J=1,MPLUS 
271 
IJ = MPLUS * (I1) + J 
272 
B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J) 
273 
END DO 
274 
END DO 
275 

276 
! COMPUTE THE NEW INTERPOLATED PARAMETER VECTOR (Geometry): 
277 
ci=0.d0 
278 
ci_tmp=0.d0 
279 

280 
ci_lt_zero= .FALSE. 
281 
DO I=1, MSET(IGeom) 
282 
DO J=1, MSET(IGeom) ! B matrix is read columnwise 
283 
ci(I)=ci(I)+B((J1)*(MPLUS)+I)*ESET(IGeom,J) !ESET is energy set. 
284 
END DO 
285 
ci(I)=ci(I)+B((MPLUS1)*(MPLUS)+I) 
286 
!Print *, 'NO ci < 0 yet, c(',I,')=', ci(I) 
287 
IF((ci(I) .LT. 0.0D0) .OR. (ci(I) .GT. 1.0D0)) THEN 
288 
ci_lt_zero=.TRUE. 
289 
EXIT 
290 
END IF 
291 
END DO !matches DO I=1, MSET(IGeom) 
292 

293 
IF (ci_lt_zero) Then 
294 
cis_zero = 0 
295 
ER_star = 0.D0 
296 
ER_star_tmp = 1e32 
297 

298 
! B_ij calculations from <B_ij=(g_ig_j)(R_iR_j)>, Full B matrix created first and then rows and columns are removed. 
299 
JJ=0 
300 
INV=NFree 
301 
DO IX=1,MSET(IGeom) 
302 
INV=INV+NFree 
303 
JNV=NFree 
304 
DO JX=1,MSET(IGeom) 
305 
JNV=JNV+NFree 
306 
JJ = JJ + 1 
307 
BST(JJ)=0.D0 
308 
DO KX=1, NFree 
309 
BST(JJ) = BST(JJ) + (((GradSet_free(IGeom,INV+KX)GradSet_free(IGeom,JNV+KX))* & 
310 
(GeomSet_free(IGeom,INV+KX)GeomSet_free(IGeom,JNV+KX)))/2.D0) 
311 
END DO 
312 
END DO 
313 
END DO 
314 

315 
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 
!Print *, 'Entering into DO I=1, (2**MSET(IGeom))2 loop, MSET(IGeom)=', MSET(IGeom), ', I=', I 
317 
ci(:)=1.D0 
318 
! find out which cis are zero in I: 
319 
DO IX=1, MSET(IGeom) 
320 
JJ=IAND(I, 2**(IX1)) 
321 
IF(JJ .EQ. 0) Then 
322 
ci(IX)=0.D0 
323 
END IF 
324 
END DO 
325 

326 
ci_lt_zero = .FALSE. 
327 
! B_ij calculations from <B_ij=(g_ig_j)(R_iR_j)>, Full B matrix created first and then rows and columns are removed. 
328 
DO IX=1, MSET(IGeom)*MSET(IGeom) 
329 
B(IX) = BST(IX) !just a copy of the original B (GEDIIS) matrix 
330 
END DO 
331 

332 
! Removal of KXth row and KXth column in order to accomodate cI to be zero: 
333 
current_size_B_mat=MSET(IGeom) 
334 
cis_zero = 0 
335 
! The bits of I (index of the upper loop 'DO I=1, (2**MSET(IGeom))2'), gives which cis are zero. 
336 
DO KX=1, MSET(IGeom) ! searching for each bit of I (index of the upper loop 'DO I=1, (2**MSET(IGeom))2') 
337 
IF (ci(KX) .EQ. 0.D0) Then !remove KXth row and KXth column 
338 
cis_zero = cis_zero + 1 
339 

340 
! First row removal: (B matrix is read columnwise) 
341 
JJ=0 
342 
DO IX=1,current_size_B_mat ! columns reading 
343 
DO JX=1,current_size_B_mat ! rows reading 
344 
IF (JX .NE. KX) Then 
345 
JJ = JJ + 1 
346 
B_tmp(JJ) = B((IX1)*current_size_B_mat+JX) 
347 
END IF 
348 
END DO 
349 
END DO 
350 

351 
DO IX=1,current_size_B_mat*(current_size_B_mat1) 
352 
B(IX) = B_tmp(IX) 
353 
END DO 
354 

355 
! Now column removal: 
356 
JJ=0 
357 
DO IX=1,KX1 ! columns reading 
358 
DO JX=1,current_size_B_mat1 ! rows reading 
359 
JJ = JJ + 1 
360 
B_tmp(JJ) = B(JJ) 
361 
END DO 
362 
END DO 
363 

364 
DO IX=KX+1,current_size_B_mat 
365 
DO JX=1,current_size_B_mat1 
366 
JJ = JJ + 1 
367 
B_tmp(JJ) = B(JJ+current_size_B_mat1) 
368 
END DO 
369 
END DO 
370 

371 
DO IX=1,(current_size_B_mat1)*(current_size_B_mat1) 
372 
B(IX) = B_tmp(IX) 
373 
END DO 
374 
current_size_B_mat = current_size_B_mat  1 
375 
END IF ! matches IF (ci(KX) .EQ. 0.D0) Then !remove 
376 
END DO ! matches DO KX=1, MSET(IGeom) 
377  
378 
! The following shifting is required to correct indices of B_ij elements in the GEDIIS matrix. 
379 
! The correction is needed because the last coloumn and row of the matrix contains all 1 and one zero. 
380 
DO IX=MSET(IGeom)cis_zero1,1,1 
381 
DO JX=MSET(IGeom)cis_zero,1,1 
382 
B(IX*(MSET(IGeom)cis_zero)+JX+IX) = B(IX*(MSET(IGeom)cis_zero)+JX) 
383 
END DO 
384 
END DO 
385 

386 
! for last row and last column of GEDIIS matrix 
387 
DO IX=1,MPLUScis_zero 
388 
B((MPLUScis_zero)*IX) = 1.D0 
389 
B((MPLUScis_zero)*(MSET(IGeom)cis_zero)+IX) = 1.D0 
390 
END DO 
391 
B((MPLUScis_zero) * (MPLUScis_zero)) = 0.D0 
392 

393 
DO IX=1, MPLUS 
394 
!WRITE(*,'(10(1X,F20.4))') B((IX1)*MPLUS+1:IX*(MPLUS)) 
395 
END DO 
396 

397 
! ELIMINATE ERROR VECTORS WITH THE LARGEST NORM: 
398 
IF (NGEDIIS .NE. MSET(IGeom)) THEN 
399 
JX = min(MSET(IGeom)NGEDIIS,MSET(IGeom)cis_zero1) 
400 
DO II=1,JX 
401 
XMAX = 1.D10 
402 
ITERA = 0 
403 
DO IX=1,MSET(IGeom)cis_zero 
404 
XNORM = 0.D0 
405 
INV = (IX1) * (MPLUScis_zero) 
406 
DO J=1,MSET(IGeom)cis_zero 
407 
XNORM = XNORM + ABS(B(INV + J)) 
408 
END DO 
409 
IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN 
410 
XMAX = XNORM 
411 
ITERA = IX 
412 
IONE = INV + IX 
413 
ENDIF 
414 
END DO 
415 

416 
DO IX=1,MPLUScis_zero 
417 
INV = (IX1) * (MPLUScis_zero) 
418 
DO J=1,MPLUScis_zero 
419 
JNV = (J1) * (MPLUScis_zero) 
420 
IF (J.EQ.ITERA) B(INV + J) = 0.D0 
421 
B(JNV + IX) = B(INV + J) 
422 
END DO 
423 
END DO 
424 
B(IONE) = 1.0D0 
425 
END DO 
426 
END IF ! matches IF (NGEDIIS .NE. MSET(IGeom)) THEN 
427  
428 
! SCALE GEDIIS MATRIX BEFORE INVERSION: 
429 
DO IX=1,MPLUScis_zero 
430 
II = (MPLUScis_zero) * (IX1) + IX ! B(II)=diagonal elements of B matrix 
431 
GSAVE(IGeom,IX) = 1.D0 / DSQRT(1.D20+DABS(B(II))) 
432 
END DO 
433 
GSAVE(IGeom,MPLUScis_zero) = 1.D0 
434 
DO IX=1,MPLUScis_zero 
435 
DO JX=1,MPLUScis_zero 
436 
IJ = (MPLUScis_zero) * (IX1) + JX 
437 
B(IJ) = B(IJ) * GSAVE(IGeom,IX) * GSAVE(IGeom,JX) 
438 
END DO 
439 
END DO 
440  
441 
! INVERT THE GEDIIS MATRIX B: 
442 
CALL MINV(B,MPLUScis_zero,DET) ! matrix inversion. 
443 

444 
DO IX=1,MPLUScis_zero 
445 
DO JX=1,MPLUScis_zero 
446 
IJ = (MPLUScis_zero) * (IX1) + JX 
447 
B(IJ) = B(IJ) * GSAVE(IGeom,IX) * GSAVE(IGeom,JX) 
448 
END DO 
449 
END DO 
450 

451 
DO IX=1, MPLUS 
452 
!WRITE(*,'(10(1X,F20.4))') B((IX1)*MPLUS+1:IX*(MPLUS)) 
453 
END DO 
454 

455 
! ESET is rearranged to handle zero cis and stored in ESET_tmp: 
456 
JJ=0 
457 
DO IX=1, MSET(IGeom) 
458 
IF (ci(IX) .NE. 0) Then 
459 
JJ=JJ+1 
460 
ESET_tmp(JJ) = ESET(IGeom,IX) 
461 
END IF 
462 
END DO 
463 

464 
! DETERMINATION OF nonzero cis: 
465 
MyPointer=1 
466 
DO IX=1, MSET(IGeom)cis_zero 
467 
tmp = 0.D0 
468 
DO J=1, MSET(IGeom)cis_zero ! B matrix is read columnwise 
469 
tmp=tmp+B((J1)*(MPLUScis_zero)+IX)*ESET_tmp(J) 
470 
END DO 
471 
tmp=tmp+B((MPLUScis_zero1)*(MPLUScis_zero)+IX) 
472 
IF((tmp .LT. 0.0D0) .OR. (tmp .GT. 1.0D0)) THEN 
473 
ci_lt_zero=.TRUE. 
474 
EXIT 
475 
ELSE 
476 
DO JX=MyPointer,MSET(IGeom) 
477 
IF (ci(JX) .NE. 0) Then 
478 
ci(JX) = tmp 
479 
MyPointer=JX+1 
480 
EXIT 
481 
END IF 
482 
END DO 
483 
END IF 
484 
END DO !matches DO I=1, MSET(IGeom)cis_zero 
485 
!Print *, 'Local set of cis, first 10:, MSET(IGeom)=', MSET(IGeom), ', I of (2**MSET(IGeom))2=', I 
486 
!WRITE(*,'(10(1X,F20.4))') ci(1:MSET(IGeom)) 
487 
!Print *, 'Local set of cis ends:****************************************' 
488  
489 
! new set of cis determined based on the lower energy (ER_star): 
490 
IF (.NOT. ci_lt_zero) Then 
491 
Call Energy_GEDIIS(MRESET,MSET(IGeom),ci,GeomSet_free(IGeom,:),GradSet_free(IGeom,:),ESET(IGeom,:),NFree,ER_star) 
492 
IF (ER_star .LT. ER_star_tmp) Then 
493 
ci_tmp=ci 
494 
ER_star_tmp = ER_star 
495 
END IF 
496 
END IF ! matches IF (.NOT. ci_lt_zero) Then 
497 
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 
ci = ci_tmp 
499 
END IF! matches IF (ci_lt_zero) Then 
500 

501 
Print *, 'Final set of cis, first 10:***********************************' 
502 
WRITE(*,'(10(1X,F20.4))') ci(1:MSET(IGeom)) 
503 
Print *, 'Final set of cis ends:****************************************' 
504 
Geom_new_free(:) = 0.D0 
505 
DO I=1, MSET(IGeom) 
506 
Geom_new_free(:) = Geom_new_free(:) + (ci(I)*GeomSet_free(IGeom,(I1)*NFree+1:I*NFree)) !MPLUS=MSET(IGeom)+1 
507 
! R_(N+1)=R*+DeltaR: 
508 
DO J=1, NFree 
509 
tmp=0.D0 
510 
DO K=1,NFree 
511 
! this can be commented: 
512 
!tmp=tmp+HFree((J1)*NFree+K)*GradSet_free(IGeom,(I1)*NFree+K) ! If Hinv=.False., then we need to invert Hess 
513 
END DO 
514 
Geom_new_free(J) = Geom_new_free(J)  (ci(I)*tmp) 
515 
END DO 
516 
END DO 
517 

518 
Step_free(:) = Geom_new_free(:)  Geom_free(:) 
519 

520 
XNORM = SQRT(DOT_PRODUCT(Step_free,Step_free)) 
521 
IF (PRINT) THEN 
522 
WRITE (6,'(/10X,''DEVIATION IN X '',F10.4,8X,''DETERMINANT '',G9.3)') XNORM, DET 
523 
!WRITE(*,'(10X,''GEDIIS COEFFICIENTS'')') 
524 
!WRITE(*,'(10X,5F12.5)') (B(MPLUS*MSET(IGeom)+I),I=1,MSET(IGeom)) 
525 
ENDIF 
526  
527 
! THE FOLLOWING TOLERENCES FOR XNORM AND DET ARE SOMEWHAT ARBITRARY! 
528 
THRES = MAX(10.D0**(NFree), 1.D25) 
529 
IF (XNORM.GT.2.D0 .OR. DABS(DET) .LT. THRES) THEN 
530 
IF (PRINT)THEN 
531 
WRITE(*,*) "THE GEDIIS MATRIX IS ILL CONDITIONED" 
532 
WRITE(*,*) "  PROBABLY, VECTORS ARE LINEARLY DEPENDENT  " 
533 
WRITE(*,*) "THE GEDIIS STEP WILL BE REPEATED WITH A SMALLER SPACE" 
534 
END IF 
535 
DO K=1,MM 
536 
B(K) = BS(K) ! why this is reverted? Because "IF (NGEDIIS .GT. 0) GO TO 80", see below 
537 
END DO 
538 
NGEDIIS = NGEDIIS  1 
539 
IF (NGEDIIS .GT. 0) GO TO 80 
540 
IF (PRINT) WRITE(*,'(10X,''NEWTONRAPHSON STEP TAKEN'')') 
541 
Geom_new_free(:) = Geom_free(:) ! Geom_new is set to original Geom, thus Step = Geom(:)  Geom_new(:)=zero, the whole 
542 
! new update to Geom_new is discarded, since XNORM.GT.2.D0 .OR. DABS(DET) .LT. THRES 
543 
END IF ! matches IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN 
544 

545 
!****************************************************************************************************************** 
546 
Geom_new_free(:) = 0.D0 
547 
DO I=1, MSET(IGeom) 
548 
Geom_new_free(:) = Geom_new_free(:) + (ci(I)*GeomSet_free(IGeom,(I1)*NFree+1:I*NFree)) !MPLUS=MSET(IGeom)+1 
549 
! R_(N+1)=R*+DeltaR: 
550 
DO J=1, NFree 
551 
tmp=0.D0 
552 
DO K=1,NFree 
553 
tmp=tmp+HFree((J1)*NFree+K)*GradSet_free(IGeom,(I1)*NFree+K) ! If Hinv=.False., then we need to invert Hess 
554 
END DO 
555 
Geom_new_free(J) = Geom_new_free(J)  (ci(I)*tmp) 
556 
END DO 
557 
END DO 
558 

559 
Step_free(:) = Geom_new_free(:)  Geom_free(:) 
560 
!****************************************************************************************************************** 
561 
Step = 0.d0 
562 
DO I=1,NFree 
563 
Step = Step + Step_free(I)*Vfree(:,I) 
564 
END DO 
565 

566 
DEALLOCATE(Hfree,Htmp,Grad_free,Step_free,Geom_free,Geom_new_free) 
567 

568 
IF (PRINT) WRITE(*,'(/,'' END Step_GEDIIS_ALL '',/)') 
569  
570 
END SUBROUTINE Step_GEDIIS_All 
571 