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)
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=1e-12
30
      REAL(KREAL), PARAMETER :: crit=1e-8
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/= NCoord-1) 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(((I-1)*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+((J-1)*NFree))=0.d0
140
           DO K=1,NCoord
141
        HFree(I+((J-1)*NFree))=HFree(I+((J-1)*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,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&
161
               GradSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
162
        GeomSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&
163
                GeomSet(IGeom,((J-1)*NCoord)+1:((J-1)*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_i-g_j)(R_i-R_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((I-1)*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 = (I-1) * 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 = (I-1) * MPLUS
234
              DO J=1,MPLUS
235
                 JNV = (J-1) * 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 * (I-1) + I ! B(II)=diagonal elements of B matrix
247
         GSAVE(IGeom,I) = 1.D0 / DSQRT(1.D-20+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 * (I-1) + 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((I-1)*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((I-1)*MPLUS+1:I*(MPLUS))
267
    END DO
268
    
269
      DO I=1,MPLUS
270
         DO J=1,MPLUS
271
            IJ = MPLUS * (I-1) + 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 column-wise
283
        ci(I)=ci(I)+B((J-1)*(MPLUS)+I)*ESET(IGeom,J) !ESET is energy set.                            
284
     END DO
285
     ci(I)=ci(I)+B((MPLUS-1)*(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_i-g_j)(R_i-R_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**(IX-1))
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_i-g_j)(R_i-R_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 column-wise)
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((IX-1)*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_mat-1)
352
                B(IX) = B_tmp(IX)
353
             END DO
354
       
355
             ! Now column removal:
356
             JJ=0   
357
                   DO IX=1,KX-1 ! columns reading
358
                      DO JX=1,current_size_B_mat-1 ! 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_mat-1
366
                 JJ = JJ + 1
367
                    B_tmp(JJ) = B(JJ+current_size_B_mat-1)
368
              END DO
369
             END DO
370
              
371
             DO IX=1,(current_size_B_mat-1)*(current_size_B_mat-1)
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_zero-1,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,MPLUS-cis_zero
388
        B((MPLUS-cis_zero)*IX) = 1.D0
389
        B((MPLUS-cis_zero)*(MSET(IGeom)-cis_zero)+IX) = 1.D0
390
       END DO
391
       B((MPLUS-cis_zero) * (MPLUS-cis_zero)) = 0.D0
392
         
393
           DO IX=1, MPLUS
394
              !WRITE(*,'(10(1X,F20.4))') B((IX-1)*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_zero-1)
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 = (IX-1) * (MPLUS-cis_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,MPLUS-cis_zero
417
                      INV = (IX-1) * (MPLUS-cis_zero)
418
                      DO J=1,MPLUS-cis_zero
419
                         JNV = (J-1) * (MPLUS-cis_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,MPLUS-cis_zero
430
        II = (MPLUS-cis_zero) * (IX-1) + IX ! B(II)=diagonal elements of B matrix
431
        GSAVE(IGeom,IX) = 1.D0 / DSQRT(1.D-20+DABS(B(II)))
432
       END DO
433
       GSAVE(IGeom,MPLUS-cis_zero) = 1.D0
434
       DO IX=1,MPLUS-cis_zero
435
        DO JX=1,MPLUS-cis_zero
436
           IJ = (MPLUS-cis_zero) * (IX-1) + 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,MPLUS-cis_zero,DET) ! matrix inversion.
443
    
444
       DO IX=1,MPLUS-cis_zero
445
        DO JX=1,MPLUS-cis_zero
446
           IJ = (MPLUS-cis_zero) * (IX-1) + 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((IX-1)*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 column-wise
469
           tmp=tmp+B((J-1)*(MPLUS-cis_zero)+IX)*ESET_tmp(J)
470
        END DO
471
            tmp=tmp+B((MPLUS-cis_zero-1)*(MPLUS-cis_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,(I-1)*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((J-1)*NFree+K)*GradSet_free(IGeom,(I-1)*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.D-25)
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,''NEWTON-RAPHSON 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,(I-1)*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((J-1)*NFree+K)*GradSet_free(IGeom,(I-1)*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