Statistiques
| Révision :

root / src / Step_DIIS_all.f90 @ 8

Historique | Voir | Annoter | Télécharger (16,32 ko)

1
!C  HEAT is never used, not even in call of Space(...)
2
!C  Geom = input parameter vector (Geometry).
3
!C  Grad = input gradient vector.
4
!C  Geom_new = New Geometry.
5

    
6
    SUBROUTINE Step_diis_all(NGeomF,IGeom,Step,Geom,Grad,HP,HEAT,Hess,NCoord,allocation_flag,Tangent)
7
!      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
8

    
9
    USE Io_module, only : IoOut 
10
    USE Path_module, only : Vfree
11
    
12
    IMPLICIT NONE
13
    INTEGER, parameter :: KINT = kind(1)
14
    INTEGER, parameter :: KREAL = kind(1.0d0)
15
    
16
    !      INCLUDE 'SIZES'
17
    
18
    INTEGER(KINT) :: NGeomF,IGeom
19
    INTEGER(KINT), INTENT(IN) :: NCoord
20
    
21
    REAL(KREAL) :: Geom(NCoord), Grad(NCoord)
22
    REAL(KREAL) :: Hess(NCoord*NCoord),Step(NCoord)
23
    REAL(KREAL) :: HEAT,HP
24
    LOGICAL :: allocation_flag
25
    REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord)
26

    
27
!************************************************************************
28
!*                                                                      *
29
!*     DIIS PERFORMS DIRECT INVERSION IN THE ITERATIVE SUBSPACE         *
30
!*                                                                      *
31
!*     THIS INVOLVES SOLVING FOR C IN Geom(NEW) = Geom' - HG'      *
32
!*                                                                      *
33
!*  WHERE Geom' = SUM(C(I)Geom(I), THE C COEFFICIENTES COMING FROM      *
34
!*                                                                      *
35
!*                   | B   1 | . | C | = | 0 |                          *
36
!*                   | 1   0 |   |-L |   | 1 |                          *
37
!*                                                                      *
38
!*  WHERE B(I,J) =GRAD(I)H(T)HGRAD(J)  GRAD(I) = GRADIENT ON CYCLE I    *
39
!*                              Hess    = INVERSE HESSIAN        *
40
!*                                                                      *
41
!*                          REFERENCE                                   *
42
!*                                                                      *
43
!*  P. CSASZAR, P. PULAY, J. MOL. STRUCT. (THEOCHEM), 114, 31 (1984)    *
44
!*                                                                      *
45
!************************************************************************
46
!************************************************************************
47
!*                                                                      *
48
!*     GEOMETRY OPTIMIZATION USING THE METHOD OF DIRECT INVERSION IN    *
49
!*     THE ITERATIVE SUBSPACE (GDIIS), COMBINED WITH THE BFGS OPTIMIZER *
50
!*     (A VARIABLE METRIC METHOD)                                       *
51
!*                                                                      *
52
!*     WRITTEN BY PETER L. CUMMINS, UNIVERSITY OF SYDNEY, AUSTRALIA     *
53
!*                                                                      *
54
!*                              REFERENCE                               *
55
!*                                                                      *
56
!*      "COMPUTATIONAL STRATEGIES FOR THE OPTIMIZATION OF EQUILIBRIUM   *
57
!*     GEOMETRIES AND TRANSITION-STATE STRUCTURES AT THE SEMIEMPIRICAL  *
58
!*     LEVEL", PETER L. CUMMINS, JILL E. GREADY, J. COMP. CHEM., 10,    *
59
!*     939-950 (1989).                                                  *
60
!*                                                                      *
61
!*     MODIFIED BY JJPS TO CONFORM TO EXISTING MOPAC CONVENTIONS        *
62
!*                                                                      *
63
!************************************************************************
64

    
65
      ! MRESET = maximum number of iterations.
66
      INTEGER(KINT), PARAMETER :: MRESET=15, M2=(MRESET+1)*(MRESET+1) !M2 = 256
67
      REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:,:),GradSet(:,:),ERR(:,:) ! MRESET*NCoord
68
      REAL(KREAL), SAVE :: ESET(MRESET)
69
    REAL(KREAL), ALLOCATABLE, SAVE :: DXTMP(:,:),GSAVE(:,:) !NCoord, why DXTMP has SAVE attribute??
70
      REAL(KREAL) :: B(M2), BS(M2)
71
      LOGICAL DEBUG, PRINT
72
      INTEGER(KINT), ALLOCATABLE, SAVE :: MSET(:)
73
    LOGICAL, ALLOCATABLE, SAVE :: FRST(:)
74
      INTEGER(KINT) :: NDIIS, MPLUS, INV, ITERA, MM, NFree, I, J, K
75
      INTEGER(KINT) :: JJ, KJ, JNV, II, IONE, IJ, INK,ITmp, Isch, Idx
76
      REAL(KREAL) :: XMax, XNorm, S, DET, THRES, Norm
77
    REAL(KREAL), PARAMETER :: eps=1e-12
78
    REAL(KREAL), PARAMETER :: crit=1e-8
79
    REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord
80
    REAL(KREAL), ALLOCATABLE :: HFree(:) ! NFree*NFree
81
    REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree
82
    REAL(KREAL), ALLOCATABLE :: Grad_free(:),Grad_new_free_inter(:),Step_free(:) ! NFree
83
    REAL(KREAL), ALLOCATABLE :: Geom_free(:),Geom_new_free_inter(:),Geom_new_free(:) ! NFree
84
    REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet_free(:,:),GradSet_free(:,:)
85
    
86
      DEBUG=.TRUE.
87
      PRINT=.TRUE.
88

    
89
      IF (PRINT)  WRITE(*,'(/,''      BEGIN GDIIS   '')')
90
    
91
      ! Initialization
92
      IF (allocation_flag) THEN ! allocation_flag = .TRUE. at the begining and effective for all geometries in path.
93
      ! FRST(IGeom) will be set to False in Space, so no need to modify it here
94
         IF (ALLOCATED(GeomSet)) THEN
95
            IF (PRINT)  WRITE(*,'(/,''    In FRST, GDIIS Dealloc  '')')
96
            DEALLOCATE(GeomSet,GradSet,ERR,DXTMP,GSave,GeomSet_free,GradSet_free)
97
            RETURN
98
         ELSE
99
        ! these allocated arrays need to be properly deallocated.
100
            IF (PRINT)  WRITE(*,'(/,''     In FRST,  GDIIS Alloc  '')')
101
            ALLOCATE(GeomSet(NGeomF,MRESET*NCoord),GradSet(NGeomF,MRESET*NCoord),ERR(NGeomF,MRESET*NCoord))
102
      ALLOCATE(GeomSet_free(NGeomF,MRESET*NCoord),GradSet_free(NGeomF,MRESET*NCoord))
103
            ALLOCATE(DXTMP(NGeomF,NCoord),GSAVE(NGeomF,NCoord),MSET(NGeomF),FRST(NGeomF))
104
      DO I=1,NGeomF
105
         FRST(I) = .TRUE.
106
         GeomSet(I,:) = 0.d0
107
         GradSet(I,:) = 0.d0
108
         ERR(I,:) = 0.d0
109
         GeomSet_free(I,:) = 0.d0
110
         GradSet_free(I,:) = 0.d0
111
         DXTMP(I,:)=0.d0
112
         GSAVE(I,:)=0.d0
113
      END DO
114
      MSET(:)=0
115
         END IF
116
     allocation_flag = .FALSE.
117
      END IF
118
    
119
      ! Addded from here:
120
    Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord)
121
      ! we orthogonalize Vfree to the tangent vector of this geom only if Tangent/=0.d0
122
      Norm=sqrt(dot_product(Tangent,Tangent))
123
      IF (Norm.GT.eps) THEN
124
         ALLOCATE(Tanf(NCoord))
125

    
126
         ! We normalize Tangent
127
         Tangent=Tangent/Norm
128

    
129
         ! We convert Tangent into Vfree only displacements. This is useless for now (2007.Apr.23)
130
     ! as Vfree=Id matrix but it will be usefull as soon as we introduce constraints.
131
         DO I=1,NCoord
132
            Tanf(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Tangent)
133
         END DO
134
         Tangent=0.d0
135
         DO I=1,NCoord
136
            Tangent=Tangent+Tanf(I)*Vfree(:,I)
137
         END DO
138
         ! first we subtract Tangent from vfree
139
         DO I=1,NCoord
140
            Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent)
141
            Vfree(:,I)=Vfree(:,I)-Norm*Tangent
142
     END DO
143

    
144
         Idx=0
145
         ! Schmidt orthogonalization of the Vfree vectors
146
         DO I=1,NCoord
147
            ! We subtract the first vectors, we do it twice as the Schmidt procedure is not numerically stable.
148
            DO Isch=1,2
149
               DO J=1,Idx
150
                  Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,J),(/NCoord/)))
151
                  Vfree(:,I)=Vfree(:,I)-Norm*Vfree(:,J)
152
               END DO
153
            END DO
154
            Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,I),(/NCoord/)))
155
            IF (Norm.GE.crit) THEN
156
               Idx=Idx+1
157
               Vfree(:,Idx)=Vfree(:,I)/sqrt(Norm)
158
            END IF
159
         END DO
160
       
161
         Print *, 'Idx=', Idx
162
         IF (Idx/= NCoord-1) THEN
163
            WRITE(*,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom
164
            WRITE(IoOut,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom
165
            STOP
166
         END IF
167
     
168
         DEALLOCATE(Tanf)
169
         NFree=Idx
170
      ELSE ! Tangent =0, matches IF (Norm.GT.eps) THEN
171
         if (debug) WRITE(*,*) "Tangent=0, using full displacement"
172
         NFree=NCoord
173
      END IF !IF (Norm.GT.eps) THEN
174
  
175
      if (debug) WRITE(*,*) 'DBG Step_DIIS_All, IGeom, NFree=', IGeom, NFree
176

    
177
      ! We now calculate the new step
178
      ! we project the hessian onto the free vectors
179
      ALLOCATE(HFree(NFree*NFree),Htmp(NCoord,NFree),Grad_free(NFree),Grad_new_free_inter(NFree))
180
    ALLOCATE(Geom_free(NFree),Step_free(NFree),Geom_new_free_inter(NFree),Geom_new_free(NFree))
181
      DO J=1,NFree
182
        DO I=1,NCoord
183
           Htmp(I,J)=0.d0
184
           DO K=1,NCoord
185
              Htmp(I,J)=Htmp(I,J)+Hess(((I-1)*NCoord)+K)*Vfree(K,J)
186
           END DO
187
        END DO
188
      END DO
189
      DO J=1,NFree
190
        DO I=1,NFree
191
       HFree(I+((J-1)*NFree))=0.d0
192
           DO K=1,NCoord
193
        HFree(I+((J-1)*NFree))=HFree(I+((J-1)*NFree))+Vfree(K,I)*Htmp(K,J)
194
           END DO
195
        END DO
196
      END DO
197

    
198
      DO I=1,NFree
199
         Grad_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Grad)
200
     Geom_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Geom)
201
      END DO
202
      !Added Ends here.***********************************************
203
    
204
!C  SPACE SIMPLY LOADS THE CURRENT VALUES OF Geom AND GRAD INTO
205
!C  THE ARRAYS GeomSet AND GradSet
206
!C  HEAT is never used, not even in Space_all(...)
207

    
208
      CALL Space_all(NGeomF,IGeom,MRESET,MSET,Geom,Grad,HEAT,NCoord,GeomSet,GradSet,ESET,FRST)
209
    
210
      IF (PRINT)  WRITE(*,'(/,''       GDIIS after Space  '')')
211
    
212
    DO J=1,MSet(IGeom)
213
         DO K=1,NFree
214
            GradSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&
215
               GradSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
216
        GeomSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&
217
                GeomSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
218
         END DO
219
    END DO    
220
!C
221
!C     INITIALIZE SOME VARIABLES AND CONSTANTS
222
!C
223
      NDIIS = MSET(IGeom)
224
      MPLUS = MSET(IGeom) + 1
225
      MM = MPLUS * MPLUS            
226
!C
227
!C     COMPUTE THE APPROXIMATE ERROR VECTORS
228
!C
229
      INV=-NFree
230
      DO 30 I=1,MSET(IGeom)
231
         INV = INV + NFree
232
         DO 30 J=1,NFree
233
            S = 0.D0
234
            KJ=(J*(J-1))/2
235
            DO 10 K=1,J
236
               KJ = KJ+1
237
   10       S = S - HFree(KJ) * GradSet_free(IGeom,INV+K)
238
            DO 20 K=J+1,NFree
239
               KJ = (K*(K-1))/2+J
240
   20       S = S - HFree(KJ) * GradSet_free(IGeom,INV+K)
241
   30 ERR(IGeom,INV+J) = S
242
   
243
!C
244
!C     CONSTRUCT THE GDIIS MATRIX
245
!C
246
      DO 40 I=1,MM
247
   40 B(I) = 1.D0
248
   
249
      JJ=0
250
      INV=-NFree
251
      DO 50 I=1,MSET(IGeom)
252
         INV=INV+NFree
253
         JNV=-NFree
254
         DO 50 J=1,MSET(IGeom)
255
            JNV=JNV+NFree
256
            JJ = JJ + 1
257
            B(JJ)=0.D0
258
            DO 50 K=1,NFree
259
      !Print *, 'B(',JJ,')=', B(JJ) + ERR(IGeom,INV+K) * ERR(IGeom,JNV+K)
260
   50 B(JJ) = B(JJ) + ERR(IGeom,INV+K) * ERR(IGeom,JNV+K)
261

    
262
     ! The following shifting is required to correct indices of B_ij elements in the GDIIS matrix.
263
   ! The correction is needed because the last coloumn of the matrix contains all 1 and one zero.
264
      DO 60 I=MSET(IGeom)-1,1,-1
265
         DO 60 J=MSET(IGeom),1,-1
266
   60 B(I*MSET(IGeom)+J+I) = B(I*MSET(IGeom)+J)
267
   
268
      ! For the last row and last column of GEDIIS matrix:
269
      DO 70 I=1,MPLUS
270
         B(MPLUS*I) = 1.D0
271
   70 B(MPLUS*MSET(IGeom)+I) = 1.D0
272
      B(MM) = 0.D0
273
!C
274
!C     ELIMINATE ERROR VECTORS WITH THE LARGEST NORM
275
!C
276
   80 CONTINUE
277
      DO 90 I=1,MM
278
   90 BS(I) = B(I)
279
   
280
      IF (NDIIS .EQ. MSET(IGeom)) GO TO 140
281
      DO 130 II=1,MSET(IGeom)-NDIIS
282
         XMAX = -1.D10
283
         ITERA = 0
284
         DO 110 I=1,MSET(IGeom)
285
            XNORM = 0.D0
286
            INV = (I-1) * MPLUS
287
            DO 100 J=1,MSET(IGeom)
288
  100       XNORM = XNORM + ABS(B(INV + J))
289
            IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN
290
               XMAX = XNORM
291
               ITERA = I
292
               IONE = INV + I
293
            ENDIF
294
  110    CONTINUE
295
  
296
         DO 120 I=1,MPLUS
297
            INV = (I-1) * MPLUS
298
            DO 120 J=1,MPLUS
299
               JNV = (J-1) * MPLUS
300
               IF (J.EQ.ITERA) B(INV + J) = 0.D0
301
               B(JNV + I) = B(INV + J)
302
         !Print *,'B(JNV + I)=',B(JNV + I)
303
  120    CONTINUE
304
         B(IONE) = 1.0D0
305
  130 CONTINUE
306
  140 CONTINUE
307
!C
308
!C     OUTPUT THE GDIIS MATRIX
309
!C
310
      IF (DEBUG) THEN
311
         WRITE(*,'(/5X,'' GDIIS MATRIX'')')
312
         ITmp=min(12,MPLUS)
313
         DO IJ=1,MPLUS
314
            WRITE(*,'(12(F12.4,1X))') B((IJ-1)*MPLUS+1:(IJ-1)*MPLUS+ITmp)
315
         END DO
316
      ENDIF
317
!C
318
!C     SCALE DIIS MATRIX BEFORE INVERSION
319
!C
320
      DO 160 I=1,MPLUS
321
         II = MPLUS * (I-1) + I
322
     !Print *, 'B(',II,')=', B(II)
323
     !Print *, 'GSave(',IGeom,',',I,')=', 1.D0 / DSQRT(1.D-20+DABS(B(II)))
324
  160 GSAVE(IGeom,I) = 1.D0 / DSQRT(1.D-20+DABS(B(II)))
325
  
326
      GSAVE(IGeom,MPLUS) = 1.D0
327
      !Print *, 'GSave(',IGeom,',',MPlus,')=1.D0'  
328
        
329
      DO 170 I=1,MPLUS
330
         DO 170 J=1,MPLUS
331
            IJ = MPLUS * (I-1) + J
332
  170 B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J)  
333
!C
334
!C     OUTPUT SCALED GDIIS MATRIX
335
!C
336
      IF (DEBUG) THEN
337
         WRITE(*,'(/5X,'' GDIIS MATRIX (SCALED)'')')
338
         ITmp=min(12,MPLUS)
339
         DO IJ=1,MPLUS
340
            WRITE(*,'(12(F12.4,1X))') B((IJ-1)*MPLUS+1:(IJ-1)*MPLUS+ITmp)
341
         END DO            
342
      ENDIF
343
!C
344
!C     INVERT THE GDIIS MATRIX B
345
!C
346
      CALL MINV(B,MPLUS,DET) ! matrix inversion.
347

    
348
      DO 190 I=1,MPLUS
349
         DO 190 J=1,MPLUS
350
            IJ = MPLUS * (I-1) + J
351
      !Print *, 'B(',IJ,')=', B(IJ)
352
      !Print *, 'GSAVE(',IGeom,',',I,')=', GSAVE(IGeom,I)
353
      !Print *, 'GSAVE(',IGeom,',',J,')=', GSAVE(IGeom,J)
354
      !Print *, 'B(',IJ,')=', B(IJ) * GSAVE(I) * GSAVE(J)
355
  190 B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J)
356
!C
357
!C     COMPUTE THE INTERMEDIATE INTERPOLATED PARAMETER AND GRADIENT VECTORS
358
!C
359
      !Print *, 'MSET(',IGeom,')=', MSET(IGeom), ' MPLUS=', MPLUS
360
      DO 200 K=1,NFree
361
         Geom_new_free_inter(K) = 0.D0
362
         Grad_new_free_inter(K) = 0.D0
363
         DO 200 I=1,MSET(IGeom)
364
            INK = (I-1) * NFree + K
365
            Geom_new_free_inter(K) = Geom_new_free_inter(K) + B(MPLUS*MSET(IGeom)+I) * GeomSet_free(IGeom,INK)
366
      !Print *, 'Geom_new_free_inter(',K,')=', Geom_new_free_inter(K)
367
      !Print *, 'B(MPLUS*MSET(',IGeom,')+',I,')=', B(MPLUS*MSET(IGeom)+I)
368
      !Print *, 'GeomSet_free(',IGeom,',',INK,')=', GeomSet_free(IGeom,INK)
369
  200 Grad_new_free_inter(K) = Grad_new_free_inter(K) + B(MPLUS*MSET(IGeom)+I) * GradSet_free(IGeom,INK)
370
      HP=0.D0
371
      DO 210 I=1,MSET(IGeom)
372
  210 HP=HP+B(MPLUS*MSET(IGeom)+I)*ESET(I)
373
      DO 220 K=1,NFree
374
  220 DXTMP(IGeom,K) = Geom_free(K) - Geom_new_free_inter(K)
375
      XNORM = SQRT(DOT_PRODUCT(DXTMP(IGeom,1:NFree),DXTMP(IGeom,1:NFree)))
376
      IF (PRINT) THEN
377
         WRITE (6,'(/10X,''DEVIATION IN X '',F10.6, 8X,''DETERMINANT '',G9.3)') XNORM,DET
378
         WRITE(*,'(10X,''GDIIS COEFFICIENTS'')')
379
         WRITE(*,'(10X,5F12.5)') (B(MPLUS*MSET(IGeom)+I),I=1,MSET(IGeom))
380
      ENDIF
381

    
382
!C     THE FOLLOWING TOLERENCES FOR XNORM AND DET ARE SOMEWHAT ARBITRARY!
383
      THRES = MAX(10.D0**(-NFree), 1.D-25)
384
      IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN
385
         IF (PRINT)THEN
386
            WRITE(*,*) "THE DIIS MATRIX IS ILL CONDITIONED"     
387
            WRITE(*,*) " - PROBABLY, VECTORS ARE LINEARLY DEPENDENT - "
388
            WRITE(*,*) "THE DIIS STEP WILL BE REPEATED WITH A SMALLER SPACE"
389
         END IF
390
         DO 230 K=1,MM
391
  230    B(K) = BS(K)
392
         NDIIS = NDIIS - 1
393
         IF (NDIIS .GT. 0) GO TO 80
394
         IF (PRINT) WRITE(*,'(10X,''NEWTON-RAPHSON STEP TAKEN'')')
395
         DO 240 K=1,NFree
396
            Geom_new_free_inter(K) = Geom_free(K)
397
  240       Grad_new_free_inter(K) = Grad_free(K)
398
      ENDIF ! matches IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN, L378
399
    
400
    ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
401
    Geom_new_free=0.d0
402
    DO I = 1, NFree
403
       DO J = 1, NFree
404
          ! If Hinv=.False., then we need to invert Hess
405
        !Geom_new_free(:) = Geom_new_free(:) + HFree(:,I)*Grad_new_free_inter(I)
406
       Geom_new_free(J) = Geom_new_free(J) + HFree(I+((J-1)*NFree))*Grad_new_free_inter(I)
407
     END DO
408
    END DO
409
    Geom_new_free(:) = Geom_new_free_inter(:) - Geom_new_free(:)
410
    
411
    Step_free = Geom_new_free - Geom_free
412
      
413
    Step = 0.d0
414
      DO I=1,NFree
415
     Step = Step + Step_free(I)*Vfree(:,I)
416
      END DO
417

    
418
      DEALLOCATE(Hfree,Htmp,Grad_free,Grad_new_free_inter,Step_free,Geom_free)
419
      DEALLOCATE(Geom_new_free_inter,Geom_new_free)
420
    
421
    IF (PRINT)  WRITE(*,'(/,''       END GDIIS  '',/)')
422
    
423
      END SUBROUTINE Step_diis_all