Révision 8 src/Opt_Geom.f90

Opt_Geom.f90 (revision 8)
3 3
!so  It has been designed to be almost independent of the rest of the code.
4 4
! It is now an (almost) officiel option.
5 5

  
6
 SUBROUTINE Opt_geom(IGeom,Nat,NCoord,Coord,Geom,E,Flag_Opt_Geom,Step_method)
6
SUBROUTINE Opt_geom(IGeom,Nat,NCoord,Coord,Geom,E,Flag_Opt_Geom,Step_method)
7 7

  
8 8
  use VarTypes
9 9
  use Path_module, only : maxcyc,sthresh,gthresh,Hinv,Smax,Nom,Atome,OrderInv,indzmat,&
10
                                            XyzGeomF,IntCoordF,UMat,ScanCoord,Coordinate,NPrim,&
11
                                            BTransInv,BTransInv_local,XyzGeomI,Xprimitive_t,ScanCoord,&
12
                                            Coordinate,Pi,UMat_local,NCart, DynMaxStep,FPrintGeom,GeomFile, AtName,Renum,Order
13
  use Io_module, only : IoGeom
10
       XyzGeomF,IntCoordF,UMat,ScanCoord,Coordinate,NPrim,&
11
       BTransInv,BTransInv_local,XyzGeomI,Xprimitive_t,ScanCoord,&
12
       Coordinate,Pi,UMat_local,NCart, DynMaxStep,FPrintGeom,GeomFile, AtName,Renum,Order, &
13
       FormAna,NbVar, PrintGeomFactor,AnaGeom
14 14

  
15
  use Io_module
16

  
15 17
  IMPLICIT NONE
16 18

  
17 19
  INTEGER(KINT), INTENT(IN) :: Nat,NCoord,IGeom
......
30 32

  
31 33
     subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old)
32 34

  
33
  ! This routines calculates the energy E and the gradient Grad of 
34
  ! a molecule with Geometry Geom (may be in internal coordinates),
35
  ! using for now, either Gaussian or Ext, more general later.
35
       ! This routines calculates the energy E and the gradient Grad of 
36
       ! a molecule with Geometry Geom (may be in internal coordinates),
37
       ! using for now, either Gaussian or Ext, more general later.
36 38

  
37
  use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, &
38
       prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, &
39
       GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local &
40
       , BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &
41
      , Order,OrderInv, XPrimitiveF
39
       use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, &
40
            prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, &
41
            GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local &
42
            , BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &
43
            , Order,OrderInv, XPrimitiveF
42 44

  
43
  ! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)
44
  ! allocated in Path.f90
45
       ! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)
46
       ! allocated in Path.f90
45 47

  
46
  use Io_module
48
       use Io_module
47 49

  
48
  ! Energy (calculated if F300K=.F., else estimated)
49
  REAL(KREAL), INTENT (OUT) :: E
50
  ! NCoord: Number of the degrees of freedom
51
  ! IGeom: index of the geometry.
52
  INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt
53
  ! Geometry at which gradient is calculated (cf Factual also):
54
  REAL(KREAL), INTENT (INOUT) :: Geom(NCoord)
55
  ! Gradient calculated at Geom geometry:
56
  REAL(KREAL), INTENT (OUT) :: Grad(NCoord)
57
  ! Cartesian geometry corresponding to (Internal Geometry) Geom:
58
  REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3)
50
       ! Energy (calculated if F300K=.F., else estimated)
51
       REAL(KREAL), INTENT (OUT) :: E
52
       ! NCoord: Number of the degrees of freedom
53
       ! IGeom: index of the geometry.
54
       INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt
55
       ! Geometry at which gradient is calculated (cf Factual also):
56
       REAL(KREAL), INTENT (INOUT) :: Geom(NCoord)
57
       ! Gradient calculated at Geom geometry:
58
       REAL(KREAL), INTENT (OUT) :: Grad(NCoord)
59
       ! Cartesian geometry corresponding to (Internal Geometry) Geom:
60
       REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3)
59 61
!!! Optional, just for geometry optimization with Baker coordinates
60
  REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3)
61
  REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord)
62
! FOptGeom is a flag indicating if we are doing a geom optimization
63
! it can be omitted so that we use a local flag: Flag_Opt_Geom
64
  LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom
65
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F
66
  LOGICAL  :: Flag_Opt_Geom
62
       REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3)
63
       REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord)
64
       ! FOptGeom is a flag indicating if we are doing a geom optimization
65
       ! it can be omitted so that we use a local flag: Flag_Opt_Geom
66
       LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom
67
       ! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F
68
       LOGICAL  :: Flag_Opt_Geom
67 69

  
68
 END subroutine Egrad
70
     END subroutine Egrad
69 71

  
70 72
  END INTERFACE
71 73

  
......
74 76
  LOGICAL :: Fini
75 77
  LOGICAL, SAVE :: FRST=.TRUE.
76 78

  
77
! Variables
79
  ! Variables
78 80
  INTEGER(KINT) :: IOpt, I,J,Iat, IBEG
79 81
  REAL(KREAL), ALLOCATABLE :: GradTmp(:),ZeroVec(:) ! 3*Nat or 3*Nat-6
80 82
  REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat-6)
......
89 91
  REAL(KREAL), ALLOCATABLE :: NewGradTmp(:)
90 92
  REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis
91 93
  REAL(KREAL) :: NormStep, FactStep, HP, MaxStep
92
  REAL(KREAL) :: Eold  
94
  REAL(KREAL) :: Eold, Eini
95
  ! Values contains the values for the geometry analizes
96
  REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar
93 97

  
94 98
  debug=valid('OptGeom')
95 99

  
96
    E=0.
97
    Eold=0.
98
    MaxStep=SMax
100
  E=0.
101
  Eold=0.
102
  Eini=0.
103
  MaxStep=SMax
99 104

  
100
    if (debug) Call Header("Entering Geom Optimization")
105
  if (debug) Call Header("Entering Geom Optimization")
101 106

  
102
     ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord))
103
     ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))
104
     ALLOCATE(GeomRef(NCoord))
105
     ALLOCATE(Hess_local(NCoord,NCoord))
106
     ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3))
107
     ALLOCATE(NewGeom(NCoord))
108
     ALLOCATE(NewGradTmp(NCoord))
109
     ALLOCATE(Hess_local_inv(NCoord,NCoord))
110
        
111
     IOpt=0
112
     ZeroVec=0.d0
113
     
114
     ! Initialize the Hessian
115
     Hess_local=0.
107
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord))
108
  ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))
109
  ALLOCATE(GeomRef(NCoord))
110
  ALLOCATE(Hess_local(NCoord,NCoord))
111
  ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3))
112
  ALLOCATE(NewGeom(NCoord))
113
  ALLOCATE(NewGradTmp(NCoord))
114
  ALLOCATE(Hess_local_inv(NCoord,NCoord))
116 115

  
117
     SELECT CASE (COORD)
118
        CASE ('ZMAT')
119
        ! We use the fact that we know that approximate good hessian values
120
        ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
121
        Hess_local(1,1)=0.5d0
122
        Hess_local(2,2)=0.5d0
123
        Hess_local(3,3)=0.1d0
124
        DO J=1,Nat-3
125
           Hess_local(3*J+1,3*J+1)=0.5d0
126
           Hess_local(3*J+2,3*J+2)=0.1d0
127
           Hess_local(3*J+3,3*J+3)=0.02d0
116
  if (NbVar>0) THEN
117
     ALLOCATE(Values(NbVar))
118
  END IF
119
  FormAna(5:8)=" I5 "
120
  IOpt=0
121
  ZeroVec=0.d0
122

  
123
  ! Initialize the Hessian
124
  Hess_local=0.
125

  
126
  SELECT CASE (COORD)
127
 CASE ('ZMAT')
128
     ! We use the fact that we know that approximate good hessian values
129
     ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
130
     Hess_local(1,1)=0.5d0
131
     Hess_local(2,2)=0.5d0
132
     Hess_local(3,3)=0.1d0
133
     DO J=1,Nat-3
134
        Hess_local(3*J+1,3*J+1)=0.5d0
135
        Hess_local(3*J+2,3*J+2)=0.1d0
136
        Hess_local(3*J+3,3*J+3)=0.02d0
137
     END DO
138
     IF (HInv) THEN
139
        DO I=1,NCoord
140
           Hess_local(I,I)=1.d0/Hess_local(I,I)
128 141
        END DO
129
        IF (HInv) THEN
130
           DO I=1,NCoord
131
              Hess_local(I,I)=1.d0/Hess_local(I,I)
132
           END DO
133
        END IF
142
     END IF
134 143

  
135
        IF (Step_method .EQ. "RFO") Then
136
                   Hess_local=0.
137
           DO I=1,NCoord
138
              Hess_local(I,I)=0.5d0
139
           END DO
140
        END IF
144
     IF (Step_method .EQ. "RFO") Then
145
        Hess_local=0.
146
        DO I=1,NCoord
147
           Hess_local(I,I)=0.5d0
148
        END DO
149
     END IF
141 150

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

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

  
225
     ! Go to optimization
226
     GeomOld=0.d0 ! Internal coordinates
227
     GeomCart_old=0.d0 ! Cartesian coordinates
180
     ! Hess = U^T Hprim U:
181
     Hess_local=0.d0
182
     DO I=1, NCoord
183
        DO J=1,NPrim
184
           ! UMat^T is needed below.
185
           Hess_local(:,I) = Hess_local(:,I) + UMat(J,:)*HprimU(J,I)
186
        END DO
187
     END DO
228 188

  
229
     IF (FPrintGeom) THEN
230
        OPEN(IOGeom,File=GeomFile)
189
     !Print *, 'Hprim='
190
     DO I=1,NPrim
191
        !   WRITE(*,'(10(1X,F6.3))') Hprim(I,:)
192
     END DO
193
     !Print *, 'UMat='
194
     DO I=1,NPrim
195
        !  WRITE(*,'(3(1X,F6.3))') UMat(I,1:3)
196
     END DO
197
     !Print *, 'HprimU='
198
     DO I=1,NPrim
199
        !  WRITE(*,'(3(1X,F6.3))') HprimU(I,1:3)
200
     END DO
201
     !Print *, 'Hess_local='
202
     DO I=1,NCoord
203
        !  WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
204
     END DO
205

  
206
     DEALLOCATE(Hprim,HprimU)
207

  
208
     IF (HInv) THEN
209
        Call GenInv(NCoord,Hess_local,Hess_local_inv,NCoord)
210
        Hess_local = Hess_local_inv
231 211
     END IF
232 212

  
233
     Fini=.FALSE.
234
     DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
235
        IOpt=IOpt+1
213
     !Print *, 'Hess_local after inversion='
214
     DO I=1,NCoord
215
        !   WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
216
     END DO
236 217

  
237
        ! Calculate the  energy and gradient
238
        IF (debug) THEN 
239
           WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
240
           WRITE(*,*) "Energy and Coordinates:"
241
           WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
242
           WRITE(*,*) "Geom Old:"
243
           WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13))
244
           WRITE(*,*) "Grad Old:"
245
           WRITE(*,'(F12.6,12(1X,F8.3))') GradOld(1:min(NCoord,13))
246
        END IF
218
     IF (Step_method .EQ. "RFO") Then
219
        Hess_local=0.
220
        DO I=1,NCoord
221
           Hess_local(I,I)=0.5d0
222
        END DO
223
     END IF
247 224

  
248
        !Call EGrad(E,Geom,GradTmp,NCoord)
249
        IF (COORD.EQ.'BAKER') THEN
250
           ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
251
           ! GeomCart has INTENT(OUT)
252
           ! GeomRef is modified in Egrag_baker, so we should not use GeomOld variable
253
           GeomRef=GeomOld
254
           Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, &
255
                GeomRef,GeomCart_old)
256
           GeomCart_old=GeomCart
225
  CASE ('MIXED','CART','HYBRID')
226
     DO J=1,NCoord
227
        Hess_local(J,J)=1.
228
     END DO
229
  CASE DEFAULT
230
     WRITE (*,*) "Coord=",TRIM(COORD)," not recognized, opt_Geom.f90, L164. Stop"
231
     STOP
232
  END SELECT
233

  
234
  ! Go to optimization
235
  GeomOld=0.d0 ! Internal coordinates
236
  GeomCart_old=0.d0 ! Cartesian coordinates
237

  
238
  IF (FPrintGeom) THEN
239
     OPEN(IOGeom,File=GeomFile)
240
  END IF
241

  
242
  Fini=.FALSE.
243
  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
244
     IOpt=IOpt+1
245

  
246
     ! Calculate the  energy and gradient
247
     IF (debug) THEN 
248
        WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
249
        WRITE(*,*) "Energy and Coordinates:"
250
        WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
251
        WRITE(*,*) "Geom Old:"
252
        WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13))
253
        WRITE(*,*) "Grad Old:"
254
        WRITE(*,'(F12.6,12(1X,F8.3))') GradOld(1:min(NCoord,13))
255
     END IF
256

  
257
     !Call EGrad(E,Geom,GradTmp,NCoord)
258
     IF (COORD.EQ.'BAKER') THEN
259
        ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
260
        ! GeomCart has INTENT(OUT)
261
        ! GeomRef is modified in Egrag_baker, so we should not use GeomOld variable
262
        GeomRef=GeomOld
263
        Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, &
264
             GeomRef,GeomCart_old)
265
        GeomCart_old=GeomCart
266
     ELSE
267
        Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart)
268
     END IF
269

  
270
     If (Iopt==1) EIni=E
271

  
272
     IF (debug) THEN 
273
        WRITE(*,*) 'L198, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
274
        WRITE(*,*) "Energy and Coordinates:"
275
        WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
276
        WRITE(*,*) "Geom Old:"
277
        WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,12))
278
        WRITE(*,*) "Grad:"
279
        WRITE(*,'(F12.6,12(1X,F8.3))') GradTmp(1:min(NCoord,12))
280
     END IF
281

  
282
     IF (FPrintGeom) THEN
283
        WRITE(IoGeom,*) Nat
284
        WRITE(IoGeom,"(' Geom. ',I5,' Energy= ',F15.6)") IOpt, E
285

  
286
        DO I=1,Nat
287
           If (renum) THEN
288
              Iat=Order(I)
289
              WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(I)),GeomCart(Iat,:),GradTmp(3*Iat-2:3*Iat)
290
           ELSE
291
              Iat=OrderInv(I)
292
              WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:),GradTmp(3*I-2:3*I)
293
           END IF
294
        END DO
295
     END IF
296

  
297
     ! do we have to analyze geometries ?
298
     If (AnaGeom) THEN
299
        If (NbVar>0) THEN
300
           Call AnalyzeGeom(GeomCart,Values)
301
           WRITE(IoDat,FormAna) Iopt,Values*PrintGeomFactor,(E-Eini)*ConvE,E
302
           if (debug) WRITE(*,FormAna) Iopt,Values*PrintGeomFactor,(E-Eini)*ConvE,E
257 303
        ELSE
258
           Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart)
304
           WRITE(IoDat,FormAna) Iopt,(E-Eini)*ConvE,E
305
           if (debug) WRITE(*,FormAna) Iopt,(E-Eini)*ConvE,E
259 306
        END IF
260
              
261
        IF (debug) THEN 
262
           WRITE(*,*) 'L198, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
263
           WRITE(*,*) "Energy and Coordinates:"
264
           WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
265
           WRITE(*,*) "Geom Old:"
266
           WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,12))
267
           WRITE(*,*) "Grad:"
268
           WRITE(*,'(F12.6,12(1X,F8.3))') GradTmp(1:min(NCoord,12))
269
        END IF
307
     END IF
270 308

  
271
        IF (FPrintGeom) THEN
272
           WRITE(IoGeom,*) Nat
273
           WRITE(IoGeom,"(' Geom. ',I5,' Energy= ',F15.6)") IOpt, E
274 309

  
275
           DO I=1,Nat
276
              If (renum) THEN
277
                 Iat=Order(I)
278
                 WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(I)),GeomCart(Iat,:),GradTmp(3*Iat-2:3*Iat)
279
              ELSE
280
                 Iat=OrderInv(I)
281
                 WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:),GradTmp(3*I-2:3*I)
282
              END IF
283
           END DO
284
        END IF
285

  
286
        IF (IOpt.GE.2) THEN
310
     IF (IOpt.GE.2) THEN
287 311
        ! We have enough geometries and gradient to update the hessian (or its 
288 312
        ! inverse)
289 313
        GradTmp2=GradTmp-GradOld
......
291 315

  
292 316
        if (debug) Write(*,*) "Call Hupdate_all, Geom"
293 317
        IF (debug) THEN
294
                 WRITE(*,*) "dx before calling",SHAPE(GeomTmp2)
295
                 WRITE(*,'(12(1X,F8.4))') GeomTmp2(1:NCoord)
296
                 WRITE(*,*) "dgrad before calling",SHAPE(GradTmp2)
297
                 WRITE(*,'(12(1X,F8.4))') GradTmp2(1:NCoord)
298
              END IF
299
           Call Hupdate_all(NCoord,GeomTmp2,GradTmp2,Hess_local)
300
        END IF ! matches IF (IOpt.GE.2) THEN
301
              
302
        GradOld=GradTmp
303
        GeomOld=Geom
318
           WRITE(*,*) "dx before calling",SHAPE(GeomTmp2)
319
           WRITE(*,'(12(1X,F8.4))') GeomTmp2(1:NCoord)
320
           WRITE(*,*) "dgrad before calling",SHAPE(GradTmp2)
321
           WRITE(*,'(12(1X,F8.4))') GradTmp2(1:NCoord)
322
        END IF
323
        Call Hupdate_all(NCoord,GeomTmp2,GradTmp2,Hess_local)
324
     END IF ! matches IF (IOpt.GE.2) THEN
304 325

  
305
        ! GradTmp becomes Step in Step_RFO_all.
306
              SELECT CASE (Step_method)
307
        CASE ('RFO')
308
          Call Step_RFO_all(NCoord,GradTmp,IGeom,Geom,GradTmp,Hess_local,ZeroVec)
309
              CASE ('GDIIS')
310
          Call Step_DIIS(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)
311
               ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
312
                     Geom=0.d0
313
            DO I=1, NCoord
314
                        ! If Hinv=.False., then we need to invert Hess_local
315
                        Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)
316
                     END DO
317
                     Geom(:) = NewGeom(:) - Geom(:)
318
                     ! GradTmp now becomes "step" below:
319
                  GradTmp = Geom - GeomOld
320
              CASE ('GDIIST')
321
          Call Step_GDIIS_Simple_Err(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,NCoord,FRST)
322
               ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
323
                     Geom=0.d0
324
            DO I=1, NCoord
325
                        ! If Hinv=.False., then we need to invert Hess_local
326
                        Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)
327
                     END DO
328
                     Geom(:) = NewGeom(:) - Geom(:)
329
                     ! GradTmp now becomes "step" below:
330
                  GradTmp = Geom - GeomOld
331
              CASE ('GEDIIS')
332
          Call Step_GEDIIS(NewGeom,Geom,GradTmp,E,Hess_local,NCoord,FRST)
333
                     ! GradTmp is actually "step" below:
334
                   GradTmp = NewGeom - GeomOld
335
                !Call Step_GEDIIS_All(3,1,GradTmp,Geom,GradTmp,E,Hess_local,NCoord,FRST,ZeroVec)
336
            CASE DEFAULT
337
          WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, opt_Geom.f90, L225. Stop"
338
          STOP
339
        END SELECT       
326
     GradOld=GradTmp
327
     GeomOld=Geom
340 328

  
341
        Fini=.TRUE.
342
        !   If (IOpt.LT.5) GradTmp=(IOpt*GradTmp+(5-IOpt)*0.01*Grad(IGeom,:))/5.d0
343
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
344
              WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh
345
        FactStep=min(1.d0,MaxStep/NormStep)
346
        Fini=(NormStep.LE.SThresh)
347
        NormStep=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients.
348
              WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormStep (Gradient)=', NormStep, 'Gradient Threshold=', GThresh
349
        Fini=(Fini.AND.(NormStep.LE.GThresh))
350
        if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep
351
        GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all.
329
     ! GradTmp becomes Step in Step_RFO_all.
330
     SELECT CASE (Step_method)
331
     CASE ('RFO')
332
        Call Step_RFO_all(NCoord,GradTmp,IGeom,Geom,GradTmp,Hess_local,ZeroVec)
333
     CASE ('GDIIS')
334
        Call Step_DIIS(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)
335
        ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
336
        Geom=0.d0
337
        DO I=1, NCoord
338
           ! If Hinv=.False., then we need to invert Hess_local
339
           Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)
340
        END DO
341
        Geom(:) = NewGeom(:) - Geom(:)
342
        ! GradTmp now becomes "step" below:
343
        GradTmp = Geom - GeomOld
344
     CASE ('GDIIST')
345
        Call Step_GDIIS_Simple_Err(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,NCoord,FRST)
346
        ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
347
        Geom=0.d0
348
        DO I=1, NCoord
349
           ! If Hinv=.False., then we need to invert Hess_local
350
           Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)
351
        END DO
352
        Geom(:) = NewGeom(:) - Geom(:)
353
        ! GradTmp now becomes "step" below:
354
        GradTmp = Geom - GeomOld
355
     CASE ('GEDIIS')
356
        Call Step_GEDIIS(NewGeom,Geom,GradTmp,E,Hess_local,NCoord,FRST)
357
        ! GradTmp is actually "step" below:
358
        GradTmp = NewGeom - GeomOld
359
        !Call Step_GEDIIS_All(3,1,GradTmp,Geom,GradTmp,E,Hess_local,NCoord,FRST,ZeroVec)
360
     CASE DEFAULT
361
        WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, opt_Geom.f90, L225. Stop"
362
        STOP
363
     END SELECT
352 364

  
365
     Fini=.TRUE.
366
     !   If (IOpt.LT.5) GradTmp=(IOpt*GradTmp+(5-IOpt)*0.01*Grad(IGeom,:))/5.d0
367
     NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
368
     if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh
369
     FactStep=min(1.d0,MaxStep/NormStep)
370
     Fini=(NormStep.LE.SThresh)
371
     NormStep=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients.
372
     if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormStep (Gradient)=', NormStep, 'Gradient Threshold=', GThresh
373
     Fini=(Fini.AND.(NormStep.LE.GThresh))
374
     if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep
375
     GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all.
353 376

  
354
        if (DynMaxStep.AND.(IOpt>1)) THEN
355
           If (E<EOld) Then
356
              MaxStep=min(1.1*MaxStep,2.*SMax)
357
           ELSE
358
              MaxStep=max(0.8*MaxStep,SMax/2.)
359
           END IF
360
           if (debug) WRITE(*,*) "Dynamically updated  Maximum step in Opt_Geom:",MaxStep,MaxStep/SMax
361 377

  
378
     if (DynMaxStep.AND.(IOpt>1)) THEN
379
        If (E<EOld) Then
380
           MaxStep=min(1.1*MaxStep,2.*SMax)
381
        ELSE
382
           MaxStep=max(0.8*MaxStep,SMax/2.)
362 383
        END IF
384
        if (debug) WRITE(*,*) "Dynamically updated  Maximum step in Opt_Geom:",MaxStep,MaxStep/SMax
363 385

  
364
        Call Check_step(Coord,Nat,NCoord,GeomOld,GradTmp)
365
       
366
        Geom=GeomOld+GradTmp ! GradTmp is step size here, from Step_RFO_all.
386
     END IF
367 387

  
368
        EOld=E
388
     Call Check_step(Coord,Nat,NCoord,GeomOld,GradTmp)
369 389

  
370
        IF (debug) THEN
371
           WRITE(*,*) "L235, DBG OptGeom, IGeom=",IGeom," IOpt=", IOpt
372
           SELECT CASE (COORD)
373
           CASE ('ZMAT','BAKER')
374
              !WRITE(*,'(12(1X,F8.3))') Geom(1:min(12,NCoord))
375
           CASE('CART','HYBRID')
376
              DO Iat=1,Nat
377
! PFL 30 Apr 2008 ->
378
! Deleted ref to orderinv that is created only for coord in zmat, baker, hybrid.
379
!                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
380
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), &
381
                      Geom(3*Iat-2:3*Iat)
382
! <- PFL 30 Apr 2008
383
              END DO
384
           CASE('MIXED')
385
             DO Iat=1,NCart
386
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
387
                      Geom(3*Iat-2:3*Iat)
388
             END DO
389
            
390
             SELECT CASE (NCart)
391
             CASE(1)
392
               if (NAt.GE.2) THEN
393
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
390
     Geom=GeomOld+GradTmp ! GradTmp is step size here, from Step_RFO_all.
391

  
392
     EOld=E
393

  
394
     IF (debug) THEN
395
        WRITE(*,*) "L235, DBG OptGeom, IGeom=",IGeom," IOpt=", IOpt
396
        SELECT CASE (COORD)
397
        CASE ('ZMAT','BAKER')
398
           !WRITE(*,'(12(1X,F8.3))') Geom(1:min(12,NCoord))
399
        CASE('CART','HYBRID')
400
           DO Iat=1,Nat
401
              ! PFL 30 Apr 2008 ->
402
              ! Deleted ref to orderinv that is created only for coord in zmat, baker, hybrid.
403
              !                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
404
              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), &
405
                   Geom(3*Iat-2:3*Iat)
406
              ! <- PFL 30 Apr 2008
407
           END DO
408
        CASE('MIXED')
409
           DO Iat=1,NCart
410
              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
411
                   Geom(3*Iat-2:3*Iat)
412
           END DO
413

  
414
           SELECT CASE (NCart)
415
           CASE(1)
416
              if (NAt.GE.2) THEN
417
                 WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
394 418
                      IndZmat(2,2),Geom(4)
395
                IBEG=5
396
               END IF
397
                IF (NAT.GE.3) THEN
398
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
419
                 IBEG=5
420
              END IF
421
              IF (NAT.GE.3) THEN
422
                 WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
399 423
                      IndZmat(3,2),Geom(5),IndZmat(3,3),Geom(6)
400
                IBeg=7
401
               END IF
402
              CASE(2)
403
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))')  Nom(Atome(OrderInv(indzmat(3,1)))), &
404
                      IndZmat(3,2),Geom(7),IndZmat(3,3),Geom(8)
405
                IBeg=9
406
              CASE DEFAULT
407
                 IBeg=3*NCart+1
408
              END SELECT
409

  
410
              DO Iat=max(4,NCart),Nat
411
                WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
412
                      IndZmat(iat,1),Geom(IBeg),IndZmat(3,2),Geom(IBeg+1), &
413
                      IndZmat(3,3),Geom(IBeg+2)*180./pi
414
                IBeg=IBeg+3
415
              END DO
416

  
424
                 IBeg=7
425
              END IF
426
           CASE(2)
427
              WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))')  Nom(Atome(OrderInv(indzmat(3,1)))), &
428
                   IndZmat(3,2),Geom(7),IndZmat(3,3),Geom(8)
429
              IBeg=9
417 430
           CASE DEFAULT
418
               WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in OptGeom. L294."
419
               STOP
431
              IBeg=3*NCart+1
420 432
           END SELECT
421
        END IF ! matches IF (debug) THEN
422
              
423
     END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
424
        
425
     DEALLOCATE(GradTmp2,GeomTmp2,GradTmp)
426
     DEALLOCATE(GradOld,GeomOld)
427
     DEALLOCATE(Hess_local)
428
        DEALLOCATE(GeomCart_old)
429
        DEALLOCATE(NewGeom,NewGradTmp)
430
     DEALLOCATE(Hess_local_inv)
431 433

  
432
     IF (Fini) THEN
433
        WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles, Energy =",E
434
     ELSE
435
        WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom," *NOT* converged in ",Iopt," cycles, Last Energy = ",E
436
     END IF
437
                
438
        WRITE(*,*) "Initial Geometry:"
439
        DO I=1, Nat
440
           WRITE(*,'(3(1X,F8.4))')  XyzGeomI(IGeom,:,I)
441
        END DO
442
         WRITE(*,*) "Final Geometry:"
443
        DO I=1, Nat
444
           WRITE(*,'(3(1X,F8.4))') GeomCart(I,:)
445
              !IF (I .GT. 1) Then
446
        !       WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)-GeomCart(I-1,1))*(GeomCart(I,1)-GeomCart(I-1,1)))&
447
              !        + ((GeomCart(I,2)-GeomCart(I-1,2))*(GeomCart(I,2)-GeomCart(I-1,2)))&
448
              !        + ((GeomCart(I,3)-GeomCart(I-1,3))*(GeomCart(I,3)-GeomCart(I-1,3))))
449
              !END IF
450
        END DO
451
        
452
        IF (COORD .EQ. "BAKER") Then
453
            I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
454
            ScanCoord=>Coordinate
455
         DO WHILE (Associated(ScanCoord%next))
456
                  I=I+1
457
                     SELECT CASE (ScanCoord%Type)
458
                     CASE ('BOND')
459
                            WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive_t(I)
460
                     CASE ('ANGLE')
461
                            WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1, &
462
                             ScanCoord%At2, ScanCoord%At3,Xprimitive_t(I)*180./Pi              
463
                     CASE ('DIHEDRAL')
464
                            WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,&
465
                            ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive_t(I)*180./Pi
466
                     END SELECT
467
                     ScanCoord => ScanCoord%next
468
               END DO ! matches DO WHILE (Associated(ScanCoord%next))
469
        END IF
470
        
471
     DEALLOCATE(GeomCart)
434
           DO Iat=max(4,NCart),Nat
435
              WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
436
                   IndZmat(iat,1),Geom(IBeg),IndZmat(3,2),Geom(IBeg+1), &
437
                   IndZmat(3,3),Geom(IBeg+2)*180./pi
438
              IBeg=IBeg+3
439
           END DO
472 440

  
473
    if (debug) Call Header("Geom Optimization Over")
441
        CASE DEFAULT
442
           WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in OptGeom. L294."
443
           STOP
444
        END SELECT
445
     END IF ! matches IF (debug) THEN
474 446

  
475
   END SUBROUTINE Opt_geom
447
  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
448

  
449
  DEALLOCATE(GradTmp2,GeomTmp2,GradTmp)
450
  DEALLOCATE(GradOld,GeomOld)
451
  DEALLOCATE(Hess_local)
452
  DEALLOCATE(GeomCart_old)
453
  DEALLOCATE(NewGeom,NewGradTmp)
454
  DEALLOCATE(Hess_local_inv)
455

  
456
  IF (Fini) THEN
457
     WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles, Energy =",E
458
  ELSE
459
     WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom," *NOT* converged in ",Iopt," cycles, Last Energy = ",E
460
  END IF
461

  
462
  WRITE(*,*) "Initial Geometry:"
463
  DO I=1, Nat
464
     WRITE(*,'(3(1X,F8.4))')  XyzGeomI(IGeom,:,I)
465
  END DO
466
  WRITE(*,*) "Final Geometry:"
467
  DO I=1, Nat
468
     WRITE(*,'(3(1X,F8.4))') GeomCart(I,:)
469
     !IF (I .GT. 1) Then
470
     !       WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)-GeomCart(I-1,1))*(GeomCart(I,1)-GeomCart(I-1,1)))&
471
     !        + ((GeomCart(I,2)-GeomCart(I-1,2))*(GeomCart(I,2)-GeomCart(I-1,2)))&
472
     !        + ((GeomCart(I,3)-GeomCart(I-1,3))*(GeomCart(I,3)-GeomCart(I-1,3))))
473
     !END IF
474
  END DO
475

  
476
  IF (COORD .EQ. "BAKER") Then
477
     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
478
     ScanCoord=>Coordinate
479
     DO WHILE (Associated(ScanCoord%next))
480
        I=I+1
481
        SELECT CASE (ScanCoord%Type)
482
        CASE ('BOND')
483
           WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive_t(I)
484
        CASE ('ANGLE')
485
           WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1, &
486
                ScanCoord%At2, ScanCoord%At3,Xprimitive_t(I)*180./Pi              
487
        CASE ('DIHEDRAL')
488
           WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,&
489
                ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive_t(I)*180./Pi
490
        END SELECT
491
        ScanCoord => ScanCoord%next
492
     END DO ! matches DO WHILE (Associated(ScanCoord%next))
493
  END IF
494

  
495
  DEALLOCATE(GeomCart)
496

  
497
  if (debug) Call Header("Geom Optimization Over")
498

  
499
END SUBROUTINE Opt_geom

Formats disponibles : Unified diff