Statistiques
| Révision :

root / src / Opt_Geom.f90 @ 8

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

1
! This subroutine optimizes a geometry. It is mainly here for  debugging purposes..
2
! It has been designed to be almost independent of the rest of the code.
3

    
4
 SUBROUTINE Opt_geom(IGeom,Nat,NCoord,Coord,Geom,E,Flag_Opt_Geom,Step_method)
5

    
6
  use Path_module, only : maxcyc,sthresh,gthresh,Hinv,Smax,Nom,Atome,OrderInv,indzmat,&
7
                                            XyzGeomF,IntCoordF,UMat,ScanCoord,Coordinate,NPrim,&
8
                                            BTransInv,BTransInv_local,XyzGeomI,Xprimitive_t,ScanCoord,&
9
                                            Coordinate,Pi,UMat_local,NCart
10
  use VarTypes
11
  IMPLICIT NONE
12

    
13
  INTEGER(KINT), INTENT(IN) :: Nat,NCoord,IGeom
14
  CHARACTER(32), INTENT(IN) :: Coord
15
  CHARACTER(32), INTENT(IN) :: Step_method
16
  REAL(KREAL), INTENT(INOUT) :: Geom(NCoord)
17
  REAL(KREAL), INTENT(OUT) :: E
18
  LOGICAL, INTENT(IN) :: Flag_Opt_Geom
19

    
20
! As this subroutine is here for debugging, debug=T !
21
  LOGICAL  :: debug
22
  LOGICAL :: Fini
23
  LOGICAL, SAVE :: FRST=.TRUE.
24

    
25
! Variables
26
  INTEGER(KINT) :: IOpt, I,J,Iat, IBEG
27
  REAL(KREAL), ALLOCATABLE :: GradTmp(:),ZeroVec(:) ! 3*Nat or 3*Nat-6
28
  REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat-6)
29
  REAL(KREAL), ALLOCATABLE :: GeomRef(:) ! (3*Nat or 3*Nat-6)
30
  REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! 3*Nat or 3*Nat-6
31
  REAL(KREAL), ALLOCATABLE :: Hess_local(:,:)
32
  REAL(KREAL), ALLOCATABLE :: GeomCart(:,:) !(Nat,3)
33
  REAL(KREAL), ALLOCATABLE :: GeomCart_old(:,:) !(Nat,3)
34
  REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim)
35
  REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat-6))
36
  REAL(KREAL), ALLOCATABLE :: NewGeom(:) !NCoord=3*Nat-6
37
  REAL(KREAL), ALLOCATABLE :: NewGradTmp(:)
38
  REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis
39
  REAL(KREAL) :: NormStep, FactStep, HP
40

    
41
  INTERFACE
42
     function valid(string) result (isValid)
43
       CHARACTER(*), intent(in) :: string
44
       logical                  :: isValid
45
     END function VALID
46

    
47
     subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old)
48

    
49
  ! This routines calculates the energy E and the gradient Grad of 
50
  ! a molecule with Geometry Geom (may be in internal coordinates),
51
  ! using for now, either Gaussian or Ext, more general later.
52

    
53
  use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, &
54
       prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, &
55
       GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local &
56
       , BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &
57
      , Order,OrderInv, XPrimitiveF
58

    
59
  ! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)
60
  ! allocated in Path.f90
61

    
62
  use Io_module
63

    
64
  ! Energy (calculated if F300K=.F., else estimated)
65
  REAL(KREAL), INTENT (OUT) :: E
66
  ! NCoord: Number of the degrees of freedom
67
  ! IGeom: index of the geometry.
68
  INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt
69
  ! Geometry at which gradient is calculated (cf Factual also):
70
  REAL(KREAL), INTENT (INOUT) :: Geom(NCoord)
71
  ! Gradient calculated at Geom geometry:
72
  REAL(KREAL), INTENT (OUT) :: Grad(NCoord)
73
  ! Cartesian geometry corresponding to (Internal Geometry) Geom:
74
  REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3)
75
!!! Optional, just for geometry optimization with Baker coordinates
76
  REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3)
77
  REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord)
78
! FOptGeom is a flag indicating if we are doing a geom optimization
79
! it can be omitted so that we use a local flag: Flag_Opt_Geom
80
  LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom
81
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F
82
  LOGICAL  :: Flag_Opt_Geom
83

    
84
 END subroutine Egrad
85

    
86
  END INTERFACE
87

    
88
  debug=valid('OptGeom')
89

    
90
    E=0.
91

    
92
    if (debug) Call Header("Entering Geom Optimization")
93

    
94
     ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord))
95
     ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))
96
     ALLOCATE(GeomRef(NCoord))
97
     ALLOCATE(Hess_local(NCoord,NCoord))
98
        ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3))
99
        ALLOCATE(NewGeom(NCoord))
100
         ALLOCATE(NewGradTmp(NCoord))
101
        ALLOCATE(Hess_local_inv(NCoord,NCoord))
102
        
103
     IOpt=0
104
     ZeroVec=0.d0
105
     
106
     ! Initialize the Hessian
107
     Hess_local=0.
108

    
109
     SELECT CASE (COORD)
110
        CASE ('ZMAT')
111
        ! We use the fact that we know that approximate good hessian values
112
        ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
113
        Hess_local(1,1)=0.5d0
114
        Hess_local(2,2)=0.5d0
115
        Hess_local(3,3)=0.1d0
116
        DO J=1,Nat-3
117
           Hess_local(3*J+1,3*J+1)=0.5d0
118
           Hess_local(3*J+2,3*J+2)=0.1d0
119
           Hess_local(3*J+3,3*J+3)=0.02d0
120
        END DO
121
        IF (HInv) THEN
122
           DO I=1,NCoord
123
              Hess_local(I,I)=1.d0/Hess_local(I,I)
124
           END DO
125
        END IF
126

    
127
        IF (Step_method .EQ. "RFO") Then
128
                   Hess_local=0.
129
           DO I=1,NCoord
130
              Hess_local(I,I)=0.5d0
131
           END DO
132
        END IF
133

    
134
     CASE ('BAKER')
135
            ! UMat(NPrim,3*Nat-6)
136
        BTransInv_local = BTransInv
137
        UMat_local = UMat
138
         ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,NCoord))
139
         Hprim=0.d0
140
         ScanCoord=>Coordinate
141
         I=0
142
         DO WHILE (Associated(ScanCoord%next))
143
            I=I+1
144
            SELECT CASE (ScanCoord%Type)
145
                     CASE ('BOND')
146
                          Hprim(I,I) = 0.5d0
147
            CASE ('ANGLE')
148
                          Hprim(I,I) = 0.2d0
149
            CASE ('DIHEDRAL')
150
                          Hprim(I,I) = 0.1d0
151
               END SELECT
152
            ScanCoord => ScanCoord%next
153
         END DO
154
               
155
         ! Hprim U:
156
         HprimU=0.d0
157
         DO I=1, NCoord
158
                  DO J=1,NPrim
159
               HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
160
            END DO
161
         END DO
162

    
163
         ! Hess = U^T Hprim U:
164
         Hess_local=0.d0
165
         DO I=1, NCoord
166
            DO J=1,NPrim
167
                     ! UMat^T is needed below.
168
                        Hess_local(:,I) = Hess_local(:,I) + UMat(J,:)*HprimU(J,I)
169
            END DO
170
         END DO
171
               
172
               !Print *, 'Hprim='
173
               DO I=1,NPrim
174
               !   WRITE(*,'(10(1X,F6.3))') Hprim(I,:)
175
            END DO
176
               !Print *, 'UMat='
177
               DO I=1,NPrim
178
                !  WRITE(*,'(3(1X,F6.3))') UMat(I,1:3)
179
            END DO
180
               !Print *, 'HprimU='
181
               DO I=1,NPrim
182
                !  WRITE(*,'(3(1X,F6.3))') HprimU(I,1:3)
183
            END DO
184
               !Print *, 'Hess_local='
185
               DO I=1,NCoord
186
                !  WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
187
            END DO
188
               
189
            DEALLOCATE(Hprim,HprimU)
190
               
191
         IF (HInv) THEN
192
                     Call GenInv(NCoord,Hess_local,Hess_local_inv,NCoord)
193
                  Hess_local = Hess_local_inv
194
        END IF
195
              
196
              !Print *, 'Hess_local after inversion='
197
               DO I=1,NCoord
198
               !   WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
199
            END DO
200
                             
201
        IF (Step_method .EQ. "RFO") Then
202
                 Hess_local=0.
203
           DO I=1,NCoord
204
              Hess_local(I,I)=0.5d0
205
           END DO
206
        END IF
207
                 
208
        CASE ('MIXED','CART','HYBRID')
209
        DO J=1,NCoord
210
           Hess_local(J,J)=1.
211
        END DO
212
        CASE DEFAULT
213
        WRITE (*,*) "Coord=",TRIM(COORD)," not recognized, opt_Geom.f90, L164. Stop"
214
        STOP
215
     END SELECT       
216

    
217
     ! Go to optimization
218
        GeomOld=0.d0 ! Internal coordinates
219
        GeomCart_old=0.d0 ! Cartesian coordinates
220

    
221
     Fini=.FALSE.
222
     DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
223
        IOpt=IOpt+1
224

    
225
        ! Calculate the  energy and gradient
226
        IF (debug) THEN 
227
           WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
228
           WRITE(*,*) "Energy and Coordinates:"
229
           WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
230
           WRITE(*,*) "Geom Old:"
231
           WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13))
232
           WRITE(*,*) "Grad Old:"
233
           WRITE(*,'(F12.6,12(1X,F8.3))') GradOld(1:min(NCoord,13))
234
        END IF
235

    
236
        !Call EGrad(E,Geom,GradTmp,NCoord)
237
        IF (COORD.EQ.'BAKER') THEN
238
           ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
239
           ! GeomCart has INTENT(OUT)
240
           ! GeomRef is modified in Egrag_baker, so we should not use GeomOld variable
241
           GeomRef=GeomOld
242
           Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, &
243
                GeomRef,GeomCart_old)
244
           GeomCart_old=GeomCart
245
        ELSE
246
           Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart)
247
        END IF
248
              
249
        IF (debug) THEN 
250
           WRITE(*,*) 'L198, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
251
           WRITE(*,*) "Energy and Coordinates:"
252
           WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
253
           WRITE(*,*) "Geom Old:"
254
           WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13))
255
           END IF
256

    
257
     IF (IOpt.GE.2) THEN
258
        ! We have enough geometries and gradient to update the hessian (or its 
259
        ! inverse)
260
        GradTmp2=GradTmp-GradOld
261
        GeomTmp2=Geom-GeomOld
262

    
263
        if (debug) Write(*,*) "Call Hupdate_all, Geom"
264
        IF (debug) THEN
265
                 WRITE(*,*) "dx before calling",SHAPE(GeomTmp2)
266
                 WRITE(*,'(12(1X,F8.4))') GeomTmp2(1:NCoord)
267
                 WRITE(*,*) "dgrad before calling",SHAPE(GradTmp2)
268
                 WRITE(*,'(12(1X,F8.4))') GradTmp2(1:NCoord)
269
              END IF
270
           Call Hupdate_all(NCoord,GeomTmp2,GradTmp2,Hess_local)
271
        END IF ! matches IF (IOpt.GE.2) THEN
272
              
273
        GradOld=GradTmp
274
        GeomOld=Geom
275

    
276
        ! GradTmp becomes Step in Step_RFO_all.
277
              SELECT CASE (Step_method)
278
        CASE ('RFO')
279
          Call Step_RFO_all(NCoord,GradTmp,IGeom,Geom,GradTmp,Hess_local,ZeroVec)
280
              CASE ('GDIIS')
281
          Call Step_DIIS(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)
282
               ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
283
                     Geom=0.d0
284
            DO I=1, NCoord
285
                        ! If Hinv=.False., then we need to invert Hess_local
286
                        Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)
287
                     END DO
288
                     Geom(:) = NewGeom(:) - Geom(:)
289
                     ! GradTmp now becomes "step" below:
290
                  GradTmp = Geom - GeomOld
291
              CASE ('GDIIST')
292
          Call Step_GDIIS_Simple_Err(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)
293
               ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
294
                     Geom=0.d0
295
            DO I=1, NCoord
296
                        ! If Hinv=.False., then we need to invert Hess_local
297
                        Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)
298
                     END DO
299
                     Geom(:) = NewGeom(:) - Geom(:)
300
                     ! GradTmp now becomes "step" below:
301
                  GradTmp = Geom - GeomOld
302
              CASE ('GEDIIS')
303
          Call Step_GEDIIS(NewGeom,Geom,GradTmp,E,Hess_local,NCoord,FRST)
304
                     ! GradTmp is actually "step" below:
305
                   GradTmp = NewGeom - GeomOld
306
                !Call Step_GEDIIS_All(3,1,GradTmp,Geom,GradTmp,E,Hess_local,NCoord,FRST,ZeroVec)
307
            CASE DEFAULT
308
          WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, opt_Geom.f90, L225. Stop"
309
          STOP
310
        END SELECT       
311

    
312
        Fini=.TRUE.
313
        !   If (IOpt.LT.5) GradTmp=(IOpt*GradTmp+(5-IOpt)*0.01*Grad(IGeom,:))/5.d0
314
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
315
              WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh
316
        FactStep=min(1.d0,SMax/NormStep)
317
        Fini=(NormStep.LE.SThresh)
318
        NormStep=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients.
319
              WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormStep (Gradient)=', NormStep, 'Gradient Threshold=', GThresh
320
        Fini=(Fini.AND.(NormStep.LE.GThresh))
321
        if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep
322
        GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all.
323

    
324
        Call Check_step(IGeom,Coord,Nat,NCoord,GeomOld,GradTmp)
325
       
326
        Geom=GeomOld+GradTmp ! GradTmp is step size here, from Step_RFO_all.
327

    
328
        IF (debug) THEN
329
           WRITE(*,*) "L235, DBG OptGeom, IGeom=",IGeom," IOpt=", IOpt
330
           SELECT CASE (COORD)
331
           CASE ('ZMAT','BAKER')
332
              !WRITE(*,'(12(1X,F8.3))') Geom(1:min(12,NCoord))
333
           CASE('CART','HYBRID')
334
              DO Iat=1,Nat
335
! PFL 30 Apr 2008 ->
336
! Deleted ref to orderinv that is created only for coord in zmat, baker, hybrid.
337
!                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
338
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), &
339
                      Geom(3*Iat-2:3*Iat)
340
! <- PFL 30 Apr 2008
341
              END DO
342
           CASE('MIXED')
343
             DO Iat=1,NCart
344
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
345
                      Geom(3*Iat-2:3*Iat)
346
             END DO
347
            
348
             SELECT CASE (NCart)
349
             CASE(1)
350
               if (NAt.GE.2) THEN
351
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
352
                      IndZmat(2,2),Geom(4)
353
                IBEG=5
354
               END IF
355
                IF (NAT.GE.3) THEN
356
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
357
                      IndZmat(3,2),Geom(5),IndZmat(3,3),Geom(6)
358
                IBeg=7
359
               END IF
360
              CASE(2)
361
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))')  Nom(Atome(OrderInv(indzmat(3,1)))), &
362
                      IndZmat(3,2),Geom(7),IndZmat(3,3),Geom(8)
363
                IBeg=9
364
              CASE DEFAULT
365
                 IBeg=3*NCart+1
366
              END SELECT
367

    
368
              DO Iat=max(4,NCart),Nat
369
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
370
                      IndZmat(iat,1),Geom(IBeg),IndZmat(3,2),Geom(IBeg+1), &
371
                      IndZmat(3,3),Geom(IBeg+2)*180./pi
372
                IBeg=IBeg+3
373
              END DO
374

    
375
           CASE DEFAULT
376
               WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in OptGeom. L294."
377
               STOP
378
           END SELECT
379
        END IF ! matches IF (debug) THEN
380
              
381
     END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
382
        
383
     DEALLOCATE(GradTmp2,GeomTmp2,GradTmp)
384
     DEALLOCATE(GradOld,GeomOld)
385
     DEALLOCATE(Hess_local)
386
        DEALLOCATE(GeomCart_old)
387
        DEALLOCATE(NewGeom,NewGradTmp)
388
     DEALLOCATE(Hess_local_inv)
389

    
390
     IF (Fini) THEN
391
        WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles, Energy =",E
392
     ELSE
393
        WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom,"*NOT* converged in ",Iopt," cycles, Last Energy = ",E
394
     END IF
395
                
396
        WRITE(*,*) "Initial Geometry:"
397
        DO I=1, Nat
398
           WRITE(*,'(3(1X,F8.4))')  XyzGeomI(IGeom,:,I)
399
        END DO
400
         WRITE(*,*) "Final Geometry:"
401
        DO I=1, Nat
402
           WRITE(*,'(3(1X,F8.4))') GeomCart(I,:)
403
              !IF (I .GT. 1) Then
404
        !       WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)-GeomCart(I-1,1))*(GeomCart(I,1)-GeomCart(I-1,1)))&
405
              !        + ((GeomCart(I,2)-GeomCart(I-1,2))*(GeomCart(I,2)-GeomCart(I-1,2)))&
406
              !        + ((GeomCart(I,3)-GeomCart(I-1,3))*(GeomCart(I,3)-GeomCart(I-1,3))))
407
              !END IF
408
        END DO
409
        
410
        IF (COORD .EQ. "BAKER") Then
411
            I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
412
            ScanCoord=>Coordinate
413
         DO WHILE (Associated(ScanCoord%next))
414
                  I=I+1
415
                     SELECT CASE (ScanCoord%Type)
416
                     CASE ('BOND')
417
                            WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive_t(I)
418
                     CASE ('ANGLE')
419
                            WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1, &
420
                             ScanCoord%At2, ScanCoord%At3,Xprimitive_t(I)*180./Pi              
421
                     CASE ('DIHEDRAL')
422
                            WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,&
423
                            ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive_t(I)*180./Pi
424
                     END SELECT
425
                     ScanCoord => ScanCoord%next
426
               END DO ! matches DO WHILE (Associated(ScanCoord%next))
427
        END IF
428
        
429
     DEALLOCATE(GeomCart)
430

    
431
    if (debug) Call Header("Geom Optimization Over")
432

    
433
   END SUBROUTINE Opt_geom