Statistiques
| Révision :

root / src / Opt_Geom.f90 @ 3

Historique | Voir | Annoter | Télécharger (13,38 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, PARAMETER :: debug=.TRUE.
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
    E=0.
42

    
43
    if (debug) WRITE(*,*) "=================== Geom Optimization =================="
44

    
45
     ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord))
46
     ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))
47
     ALLOCATE(GeomRef(NCoord))
48
     ALLOCATE(Hess_local(NCoord,NCoord))
49
	 ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3))
50
	 ALLOCATE(NewGeom(NCoord))
51
 	 ALLOCATE(NewGradTmp(NCoord))
52
	 ALLOCATE(Hess_local_inv(NCoord,NCoord))
53
	 
54
     IOpt=0
55
     ZeroVec=0.d0
56
     
57
     ! Initialize the Hessian
58
     Hess_local=0.
59

    
60
     SELECT CASE (COORD)
61
	 CASE ('ZMAT')
62
        ! We use the fact that we know that approximate good hessian values
63
        ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
64
        Hess_local(1,1)=0.5d0
65
        Hess_local(2,2)=0.5d0
66
        Hess_local(3,3)=0.1d0
67
        DO J=1,Nat-3
68
           Hess_local(3*J+1,3*J+1)=0.5d0
69
           Hess_local(3*J+2,3*J+2)=0.1d0
70
           Hess_local(3*J+3,3*J+3)=0.02d0
71
        END DO
72
        IF (HInv) THEN
73
           DO I=1,NCoord
74
              Hess_local(I,I)=1.d0/Hess_local(I,I)
75
           END DO
76
        END IF
77

    
78
        IF (Step_method .EQ. "RFO") Then
79
                   Hess_local=0.
80
           DO I=1,NCoord
81
              Hess_local(I,I)=0.5d0
82
           END DO
83
        END IF
84

    
85
	 CASE ('BAKER')
86
	     ! UMat(NPrim,3*Nat-6)
87
		 BTransInv_local = BTransInv
88
		 UMat_local = UMat
89
         ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,NCoord))
90
         Hprim=0.d0
91
         ScanCoord=>Coordinate
92
         I=0
93
         DO WHILE (Associated(ScanCoord%next))
94
            I=I+1
95
            SELECT CASE (ScanCoord%Type)
96
			CASE ('BOND')
97
			     Hprim(I,I) = 0.5d0
98
            CASE ('ANGLE')
99
			     Hprim(I,I) = 0.2d0
100
            CASE ('DIHEDRAL')
101
			     Hprim(I,I) = 0.1d0
102
	        END SELECT
103
            ScanCoord => ScanCoord%next
104
         END DO
105
		 
106
         ! Hprim U:
107
         HprimU=0.d0
108
         DO I=1, NCoord
109
		    DO J=1,NPrim
110
               HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
111
            END DO
112
         END DO
113

    
114
         ! Hess = U^T Hprim U:
115
         Hess_local=0.d0
116
         DO I=1, NCoord
117
            DO J=1,NPrim
118
		       ! UMat^T is needed below.
119
			   Hess_local(:,I) = Hess_local(:,I) + UMat(J,:)*HprimU(J,I)
120
            END DO
121
         END DO
122
		 
123
		 !Print *, 'Hprim='
124
		 DO I=1,NPrim
125
		 !   WRITE(*,'(10(1X,F6.3))') Hprim(I,:)
126
	     END DO
127
		 !Print *, 'UMat='
128
		 DO I=1,NPrim
129
		  !  WRITE(*,'(3(1X,F6.3))') UMat(I,1:3)
130
	     END DO
131
		 !Print *, 'HprimU='
132
		 DO I=1,NPrim
133
		  !  WRITE(*,'(3(1X,F6.3))') HprimU(I,1:3)
134
	     END DO
135
		 !Print *, 'Hess_local='
136
		 DO I=1,NCoord
137
		  !  WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
138
	     END DO
139
		 
140
	     DEALLOCATE(Hprim,HprimU)
141
		 
142
         IF (HInv) THEN
143
			Call GenInv(NCoord,Hess_local,Hess_local_inv,NCoord)
144
		    Hess_local = Hess_local_inv
145
        END IF
146
		
147
		!Print *, 'Hess_local after inversion='
148
		 DO I=1,NCoord
149
		 !   WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
150
	     END DO
151
				 
152
        IF (Step_method .EQ. "RFO") Then
153
		   Hess_local=0.
154
           DO I=1,NCoord
155
              Hess_local(I,I)=0.5d0
156
           END DO
157
        END IF
158
		   
159
	 CASE ('MIXED','CART','HYBRID')
160
        DO J=1,NCoord
161
           Hess_local(J,J)=1.
162
        END DO
163
	 CASE DEFAULT
164
        WRITE (*,*) "Coord=",TRIM(COORD)," not recognized, opt_Geom.f90, L164. Stop"
165
        STOP
166
     END SELECT	
167

    
168
     ! Go to optimization
169
	 GeomOld=0.d0 ! Internal coordinates
170
	 GeomCart_old=0.d0 ! Cartesian coordinates
171

    
172
     Fini=.FALSE.
173
     DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
174
        IOpt=IOpt+1
175

    
176
        ! Calculate the  energy and gradient
177
        IF (debug) THEN 
178
           WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
179
           WRITE(*,*) "Energy and Coordinates:"
180
           WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
181
           WRITE(*,*) "Geom Old:"
182
           WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13))
183
           WRITE(*,*) "Grad Old:"
184
           WRITE(*,'(F12.6,12(1X,F8.3))') GradOld(1:min(NCoord,13))
185
        END IF
186

    
187
        !Call EGrad(E,Geom,GradTmp,NCoord)
188
        IF (COORD.EQ.'BAKER') THEN
189
           ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
190
           ! GeomCart has INTENT(OUT)
191
           ! GeomRef is modified in Egrag_baker, so we should not use GeomOld variable
192
           GeomRef=GeomOld
193
           Call EGrad_baker(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, &
194
                GeomRef,GeomCart_old)
195
           GeomCart_old=GeomCart
196
        ELSE
197
           Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart)
198
        END IF
199
		
200
        IF (debug) THEN 
201
           WRITE(*,*) 'L198, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
202
           WRITE(*,*) "Energy and Coordinates:"
203
           WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
204
           WRITE(*,*) "Geom Old:"
205
           WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13))
206
	    END IF
207

    
208
     IF (IOpt.GE.2) THEN
209
        ! We have enough geometries and gradient to update the hessian (or its 
210
        ! inverse)
211
        GradTmp2=GradTmp-GradOld
212
        GeomTmp2=Geom-GeomOld
213

    
214
        if (debug) Write(*,*) "Call Hupdate_all, Geom"
215
        IF (debug) THEN
216
		   WRITE(*,*) "dx before calling",SHAPE(GeomTmp2)
217
		   WRITE(*,'(12(1X,F8.4))') GeomTmp2(1:NCoord)
218
		   WRITE(*,*) "dgrad before calling",SHAPE(GradTmp2)
219
		   WRITE(*,'(12(1X,F8.4))') GradTmp2(1:NCoord)
220
		END IF
221
           Call Hupdate_all(NCoord,GeomTmp2,GradTmp2,Hess_local)
222
        END IF ! matches IF (IOpt.GE.2) THEN
223
		
224
        GradOld=GradTmp
225
        GeomOld=Geom
226

    
227
        ! GradTmp becomes Step in Step_RFO_all.
228
		SELECT CASE (Step_method)
229
        CASE ('RFO')
230
          Call Step_RFO_all(NCoord,GradTmp,IGeom,Geom,GradTmp,Hess_local,ZeroVec)
231
		CASE ('GDIIS')
232
          Call Step_DIIS(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)
233
		 ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
234
			Geom=0.d0
235
            DO I=1, NCoord
236
			   ! If Hinv=.False., then we need to invert Hess_local
237
			   Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)
238
			END DO
239
			Geom(:) = NewGeom(:) - Geom(:)
240
			! GradTmp now becomes "step" below:
241
		    GradTmp = Geom - GeomOld
242
		CASE ('GDIIST')
243
          Call Step_GDIIS_Simple_Err(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)
244
		 ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
245
			Geom=0.d0
246
            DO I=1, NCoord
247
			   ! If Hinv=.False., then we need to invert Hess_local
248
			   Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)
249
			END DO
250
			Geom(:) = NewGeom(:) - Geom(:)
251
			! GradTmp now becomes "step" below:
252
		    GradTmp = Geom - GeomOld
253
		CASE ('GEDIIS')
254
          Call Step_GEDIIS(NewGeom,Geom,GradTmp,E,Hess_local,NCoord,FRST)
255
			! GradTmp is actually "step" below:
256
		     GradTmp = NewGeom - GeomOld
257
		  !Call Step_GEDIIS_All(3,1,GradTmp,Geom,GradTmp,E,Hess_local,NCoord,FRST,ZeroVec)
258
 	    CASE DEFAULT
259
          WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, opt_Geom.f90, L225. Stop"
260
          STOP
261
        END SELECT	
262

    
263
        Fini=.TRUE.
264
        !   If (IOpt.LT.5) GradTmp=(IOpt*GradTmp+(5-IOpt)*0.01*Grad(IGeom,:))/5.d0
265
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
266
		WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh
267
        FactStep=min(1.d0,SMax/NormStep)
268
        Fini=(NormStep.LE.SThresh)
269
        NormStep=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients.
270
		WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormStep (Gradient)=', NormStep, 'Gradient Threshold=', GThresh
271
        Fini=(Fini.AND.(NormStep.LE.GThresh))
272
        if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep
273
        GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all.
274

    
275
        Call Check_step(IGeom,Coord,Nat,NCoord,GeomOld,GradTmp)
276
       
277
        Geom=GeomOld+GradTmp ! GradTmp is step size here, from Step_RFO_all.
278

    
279
        IF (debug) THEN
280
           WRITE(*,*) "L235, DBG OptGeom, IGeom=",IGeom," IOpt=", IOpt
281
           SELECT CASE (COORD)
282
           CASE ('ZMAT','BAKER')
283
              !WRITE(*,'(12(1X,F8.3))') Geom(1:min(12,NCoord))
284
           CASE('CART','HYBRID')
285
              DO Iat=1,Nat
286
! PFL 30 Apr 2008 ->
287
! Deleted ref to orderinv that is created only for coord in zmat, baker, hybrid.
288
!                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
289
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), &
290
                      Geom(3*Iat-2:3*Iat)
291
! <- PFL 30 Apr 2008
292
              END DO
293
           CASE('MIXED')
294
             DO Iat=1,NCart
295
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
296
                      Geom(3*Iat-2:3*Iat)
297
             END DO
298
            
299
             SELECT CASE (NCart)
300
             CASE(1)
301
               if (NAt.GE.2) THEN
302
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
303
                      IndZmat(2,2),Geom(4)
304
                IBEG=5
305
               END IF
306
                IF (NAT.GE.3) THEN
307
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
308
                      IndZmat(3,2),Geom(5),IndZmat(3,3),Geom(6)
309
                IBeg=7
310
               END IF
311
              CASE(2)
312
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))')  Nom(Atome(OrderInv(indzmat(3,1)))), &
313
                      IndZmat(3,2),Geom(7),IndZmat(3,3),Geom(8)
314
                IBeg=9
315
              CASE DEFAULT
316
                 IBeg=3*NCart+1
317
              END SELECT
318

    
319
              DO Iat=max(4,NCart),Nat
320
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
321
                      IndZmat(iat,1),Geom(IBeg),IndZmat(3,2),Geom(IBeg+1), &
322
                      IndZmat(3,3),Geom(IBeg+2)*180./pi
323
                IBeg=IBeg+3
324
              END DO
325

    
326
           CASE DEFAULT
327
               WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in OptGeom. L294."
328
               STOP
329
           END SELECT
330
        END IF ! matches IF (debug) THEN
331
		
332
     END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
333
	 
334
     DEALLOCATE(GradTmp2,GeomTmp2,GradTmp)
335
     DEALLOCATE(GradOld,GeomOld)
336
     DEALLOCATE(Hess_local)
337
	 DEALLOCATE(GeomCart_old)
338
	 DEALLOCATE(NewGeom,NewGradTmp)
339
     DEALLOCATE(Hess_local_inv)
340

    
341
     IF (Fini) THEN
342
        WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles, Energy =",E
343
     ELSE
344
        WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom,"*NOT* converged in ",Iopt," cycles, Last Energy = ",E
345
     END IF
346
	 	 
347
	 WRITE(*,*) "Initial Geometry:"
348
	 DO I=1, Nat
349
	    WRITE(*,'(3(1X,F8.4))')  XyzGeomI(IGeom,:,I)
350
	 END DO
351
	  WRITE(*,*) "Final Geometry:"
352
	 DO I=1, Nat
353
	    WRITE(*,'(3(1X,F8.4))') GeomCart(I,:)
354
		!IF (I .GT. 1) Then
355
        !	WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)-GeomCart(I-1,1))*(GeomCart(I,1)-GeomCart(I-1,1)))&
356
		!	 + ((GeomCart(I,2)-GeomCart(I-1,2))*(GeomCart(I,2)-GeomCart(I-1,2)))&
357
		!	 + ((GeomCart(I,3)-GeomCart(I-1,3))*(GeomCart(I,3)-GeomCart(I-1,3))))
358
		!END IF
359
	 END DO
360
	 
361
	 IF (COORD .EQ. "BAKER") Then
362
    	 I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
363
	     ScanCoord=>Coordinate
364
         DO WHILE (Associated(ScanCoord%next))
365
		    I=I+1
366
			SELECT CASE (ScanCoord%Type)
367
			CASE ('BOND')
368
				WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive_t(I)
369
			CASE ('ANGLE')
370
				WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1, &
371
				 ScanCoord%At2, ScanCoord%At3,Xprimitive_t(I)*180./Pi		
372
			CASE ('DIHEDRAL')
373
				WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,&
374
				ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive_t(I)*180./Pi
375
			END SELECT
376
			ScanCoord => ScanCoord%next
377
		 END DO ! matches DO WHILE (Associated(ScanCoord%next))
378
	 END IF
379
	 
380
     DEALLOCATE(GeomCart)
381
	 if (debug) WRITE(*,*) "================== Geom Optimization Over================="
382
   END SUBROUTINE Opt_geom