Statistiques
| Révision :

root / src / Opt_Geom.f90 @ 9

Historique | Voir | Annoter | Télécharger (18,94 ko)

1
! This subroutine optimizes a geometry. 
2
! Initially, it was mainly here for  debugging purposes..
3
!so  It has been designed to be almost independent of the rest of the code.
4
! It is now an (almost) officiel option.
5

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

    
8
  use VarTypes
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
       FormAna,NbVar, PrintGeomFactor,AnaGeom
14

    
15
  use Io_module
16

    
17
  IMPLICIT NONE
18

    
19
  INTEGER(KINT), INTENT(IN) :: Nat,NCoord,IGeom
20
  CHARACTER(32), INTENT(IN) :: Coord
21
  CHARACTER(32), INTENT(IN) :: Step_method
22
  REAL(KREAL), INTENT(INOUT) :: Geom(NCoord)
23
  REAL(KREAL), INTENT(OUT) :: E
24
  LOGICAL, INTENT(IN) :: Flag_Opt_Geom
25

    
26

    
27
  INTERFACE
28
     function valid(string) result (isValid)
29
       CHARACTER(*), intent(in) :: string
30
       logical                  :: isValid
31
     END function VALID
32

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

    
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.
38

    
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
44

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

    
48
       use Io_module
49

    
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)
61
!!! Optional, just for geometry optimization with Baker coordinates
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
69

    
70
     END subroutine Egrad
71

    
72
  END INTERFACE
73

    
74

    
75
  LOGICAL  :: debug
76
  LOGICAL :: FiniS,FiniG,Fini
77
  LOGICAL, SAVE :: FRST=.TRUE.
78

    
79
  ! Variables
80
  INTEGER(KINT) :: IOpt, I,J,Iat, IBEG
81
  REAL(KREAL), ALLOCATABLE :: GradTmp(:),ZeroVec(:) ! 3*Nat or 3*Nat-6
82
  REAL(KREAL), ALLOCATABLE :: Step(:) ! 3*Nat or 3*Nat-6
83
  REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat-6)
84
  REAL(KREAL), ALLOCATABLE :: GeomRef(:) ! (3*Nat or 3*Nat-6)
85
  REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! 3*Nat or 3*Nat-6
86
  REAL(KREAL), ALLOCATABLE :: Hess_local(:,:)
87
  REAL(KREAL), ALLOCATABLE :: GeomCart(:,:) !(Nat,3)
88
  REAL(KREAL), ALLOCATABLE :: GeomCart_old(:,:) !(Nat,3)
89
  REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim)
90
  REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat-6))
91
  REAL(KREAL), ALLOCATABLE :: NewGeom(:) !NCoord=3*Nat-6
92
  REAL(KREAL), ALLOCATABLE :: NewGradTmp(:)
93
  REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis
94
  REAL(KREAL) :: NormStep, FactStep, HP, MaxStep,NormGrad
95
  REAL(KREAL) :: Eold, Eini
96
  ! Values contains the values for the geometry analizes
97
  REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar
98
  CHARACTER(LCHARS) :: Line
99

    
100
  debug=valid('OptGeom')
101

    
102
  E=0.
103
  Eold=0.
104
  Eini=0.
105
  MaxStep=SMax
106

    
107
  if (debug) Call Header("Entering Geom Optimization")
108

    
109
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord),Step(NCoord))
110
  ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))
111
  ALLOCATE(GeomRef(NCoord))
112
  ALLOCATE(Hess_local(NCoord,NCoord))
113
  ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3))
114
  ALLOCATE(NewGeom(NCoord))
115
  ALLOCATE(NewGradTmp(NCoord))
116
  ALLOCATE(Hess_local_inv(NCoord,NCoord))
117

    
118
  if (NbVar>0) THEN
119
     ALLOCATE(Values(NbVar))
120
  END IF
121
  FormAna(5:8)=" I5 "
122
  IOpt=0
123
  ZeroVec=0.d0
124

    
125
  ! Initialize the Hessian
126
  Hess_local=0.
127

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

    
146
     IF (Step_method .EQ. "RFO") Then
147
        Hess_local=0.
148
        DO I=1,NCoord
149
           Hess_local(I,I)=1.d0
150
        END DO
151
     END IF
152

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

    
174
     ! Hprim U:
175
     HprimU=0.d0
176
     DO I=1, NCoord
177
        DO J=1,NPrim
178
           HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
179
        END DO
180
     END DO
181

    
182
     ! Hess = U^T Hprim U:
183
     Hess_local=0.d0
184
     DO I=1, NCoord
185
        DO J=1,NPrim
186
           ! UMat^T is needed below.
187
           Hess_local(:,I) = Hess_local(:,I) + UMat(J,:)*HprimU(J,I)
188
        END DO
189
     END DO
190

    
191
     !Print *, 'Hprim='
192
     DO I=1,NPrim
193
        !   WRITE(*,'(10(1X,F6.3))') Hprim(I,:)
194
     END DO
195
     !Print *, 'UMat='
196
     DO I=1,NPrim
197
        !  WRITE(*,'(3(1X,F6.3))') UMat(I,1:3)
198
     END DO
199
     !Print *, 'HprimU='
200
     DO I=1,NPrim
201
        !  WRITE(*,'(3(1X,F6.3))') HprimU(I,1:3)
202
     END DO
203
     !Print *, 'Hess_local='
204
     DO I=1,NCoord
205
        !  WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
206
     END DO
207

    
208
     DEALLOCATE(Hprim,HprimU)
209

    
210
     IF (HInv) THEN
211
        Call GenInv(NCoord,Hess_local,Hess_local_inv,NCoord)
212
        Hess_local = Hess_local_inv
213
     END IF
214

    
215
     !Print *, 'Hess_local after inversion='
216
!     DO I=1,NCoord
217
        !   WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
218
!     END DO
219

    
220
     IF (Step_method .EQ. "RFO") Then
221
        Hess_local=0.
222
        DO I=1,NCoord
223
           Hess_local(I,I)=0.5d0
224
        END DO
225
     END IF
226

    
227
  CASE ('MIXED','CART','HYBRID')
228
     DO J=1,NCoord
229
        Hess_local(J,J)=1.
230
     END DO
231
  CASE DEFAULT
232
     WRITE (*,*) "Coord=",TRIM(COORD)," not recognized, opt_Geom.f90, L164. Stop"
233
     STOP
234
  END SELECT
235

    
236
  ! Go to optimization
237
  GeomOld=0.d0 ! Internal coordinates
238
  GeomCart_old=0.d0 ! Cartesian coordinates
239

    
240
  IF (FPrintGeom) THEN
241
     OPEN(IOGeom,File=GeomFile)
242
  END IF
243

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

    
248
     Write(Line,'(1X,A,I5)') "Iteration ",Iopt
249
     Call Header(TRIM(Line))
250
     WRITE(IoOut,*) "Current Geometry"
251
     Line=""
252
     Call PrintGeom(Line,Nat,NCoord,Geom,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)
253
     if (COORD/="CART") THEN
254
        WRITE(IoOut,*) "Current Geometry in Cartesian coordinates"
255
        Call PrintGeom(Line,Nat,NCoord,GeomCart,"CART",IOOUT,Atome,Order,OrderInv,IndZmat)
256
     END IF
257

    
258
     WRITE(IoOut,*) "Computing energy and gradient"
259
     ! Calculate the  energy and gradient
260
     IF (debug) THEN 
261
        WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
262
        WRITE(*,*) "Energy and Coordinates:"
263
        WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
264
        WRITE(*,*) "Geom Old:"
265
        WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13))
266
        WRITE(*,*) "Grad Old:"
267
        WRITE(*,'(F12.6,12(1X,F8.3))') GradOld(1:min(NCoord,13))
268
     END IF
269

    
270
     !Call EGrad(E,Geom,GradTmp,NCoord)
271
     IF (COORD.EQ.'BAKER') THEN
272
        ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
273
        ! GeomCart has INTENT(OUT)
274
        ! GeomRef is modified in Egrag_baker, so we should not use GeomOld variable
275
        GeomRef=GeomOld
276
        Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, &
277
             GeomRef,GeomCart_old)
278
        GeomCart_old=GeomCart
279
     ELSE
280
        Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart)
281
     END IF
282
!!!!!!!!!!! PFL March 2013 !!!!!!!!!!!!!
283
! We have a pb with the format of GradTmp(Nat,3) whereas it is 3,Nat in Egrad
284
!
285
! This is a path for CART coordinates !!!
286
     IF (COORD=="CART") THEN
287
        Step=reshape(reshape(GradTmp,(/3,Nat/),ORDER=(/2,1/)),(/3*Nat/))
288
        GradTmp=Step
289
     END IF
290
!!!!!!!!!!!!!!!!!! PFL March 2013 !!!!!!!!!!!!!!!!!!!!
291

    
292
     WRITE(IoOut,'(1X," Energy E=",F15.6,A)') E,TRIM(UnitProg)
293
     WRITE(IoOut,*) "Gradient for current geometry"
294
     Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)
295

    
296
     If (Iopt==1) EIni=E
297

    
298
     IF (debug) THEN 
299
        WRITE(*,*) 'L198, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
300
        WRITE(*,*) "Energy and Coordinates:"
301
        WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
302
        WRITE(*,*) "Geom Old:"
303
        WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,12))
304
        WRITE(*,*) "Grad:"
305
        WRITE(*,'(F12.6,12(1X,F8.3))') GradTmp(1:min(NCoord,12))
306
     END IF
307

    
308
     IF (FPrintGeom) THEN
309
        WRITE(IoGeom,*) Nat
310
        WRITE(IoGeom,"(' Geom. ',I5,' Energy= ',F15.6)") IOpt, E
311

    
312
        DO I=1,Nat
313
           If (renum) THEN
314
              Iat=Order(I)
315
              WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(I)),GeomCart(Iat,:),GradTmp(3*Iat-2:3*Iat)
316
           ELSE
317
              Iat=OrderInv(I)
318
              WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:),GradTmp(3*I-2:3*I)
319
           END IF
320
        END DO
321
     END IF
322

    
323
     ! do we have to analyze geometries ?
324
     If (AnaGeom) THEN
325
        If (NbVar>0) THEN
326
           Call AnalyzeGeom(GeomCart,Values)
327
           WRITE(IoDat,FormAna) Iopt,Values*PrintGeomFactor,(E-Eini)*ConvE,E
328
           if (debug) WRITE(*,FormAna) Iopt,Values*PrintGeomFactor,(E-Eini)*ConvE,E
329
        ELSE
330
           WRITE(IoDat,FormAna) Iopt,(E-Eini)*ConvE,E
331
           if (debug) WRITE(*,FormAna) Iopt,(E-Eini)*ConvE,E
332
        END IF
333
     END IF
334

    
335

    
336
     IF (IOpt.GE.2) THEN
337
        ! We have enough geometries and gradient to update the hessian (or its 
338
        ! inverse)
339
        GradTmp2=GradTmp-GradOld
340
        GeomTmp2=Geom-GeomOld
341

    
342
        if (debug) Write(*,*) "Call Hupdate_all, Geom"
343
        IF (debug) THEN
344
           WRITE(*,*) "dx before calling",SHAPE(GeomTmp2)
345
           WRITE(*,'(12(1X,F8.4))') GeomTmp2(1:NCoord)
346
           WRITE(*,*) "dgrad before calling",SHAPE(GradTmp2)
347
           WRITE(*,'(12(1X,F8.4))') GradTmp2(1:NCoord)
348
        END IF
349
        Call Hupdate_all(NCoord,GeomTmp2,GradTmp2,Hess_local)
350
     END IF ! matches IF (IOpt.GE.2) THEN
351

    
352
     GradOld=GradTmp
353
     GeomOld=Geom
354

    
355
     ! GradTmp becomes Step in Step_RFO_all.
356
     SELECT CASE (Step_method)
357
     CASE ('RFO')
358
        Call Step_RFO_all(NCoord,Step,IGeom,Geom,GradTmp,Hess_local,ZeroVec)
359
        GradTmp=Step
360
     CASE ('GDIIS')
361
        Call Step_DIIS(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)
362
        ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
363
        Geom=0.d0
364
        DO I=1, NCoord
365
           ! If Hinv=.False., then we need to invert Hess_local
366
           Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)
367
        END DO
368
        Geom(:) = NewGeom(:) - Geom(:)
369
        ! GradTmp now becomes "step" below:
370
        GradTmp = Geom - GeomOld
371
     CASE ('GDIIST')
372
        Call Step_GDIIS_Simple_Err(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,NCoord,FRST)
373
        ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
374
        Geom=0.d0
375
        DO I=1, NCoord
376
           ! If Hinv=.False., then we need to invert Hess_local
377
           Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)
378
        END DO
379
        Geom(:) = NewGeom(:) - Geom(:)
380
        ! GradTmp now becomes "step" below:
381
        GradTmp = Geom - GeomOld
382
     CASE ('GEDIIS')
383
        Call Step_GEDIIS(NewGeom,Geom,GradTmp,E,Hess_local,NCoord,FRST)
384
        ! GradTmp is actually "step" below:
385
        GradTmp = NewGeom - GeomOld
386
        !Call Step_GEDIIS_All(3,1,GradTmp,Geom,GradTmp,E,Hess_local,NCoord,FRST,ZeroVec)
387
     CASE DEFAULT
388
        WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, opt_Geom.f90, L225. Stop"
389
        STOP
390
     END SELECT
391

    
392
     Fini=.TRUE.
393
     !   If (IOpt.LT.5) GradTmp=(IOpt*GradTmp+(5-IOpt)*0.01*Grad(IGeom,:))/5.d0
394
     NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
395
     if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh
396
     FactStep=min(1.d0,MaxStep/NormStep)
397
     FiniS=((NormStep*FactStep).LE.SThresh)
398
     NormGrad=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients.
399
     if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormGrad (Gradient)=', NormGrad, 'Gradient Threshold=', GThresh
400
     FiniG=(NormGrad.LE.GThresh)
401
     Fini=FiniS.AND.FiniG
402
     if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep
403
     GradTmp=GradTmp*FactStep ! GradTmp is step  here, from Step_RFO_all.
404

    
405
     WRITE(IoOut,*) " Converged ?"
406
     WRITE(IoOut,'(1X,A18,1X,A15,1X,A13,5X,A3)') " Criterion","Current Value","Threshold","Ok?"
407
     IF (FiniS) THEN
408
        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Stepsize :",NormStep,SThresh,"YES"
409
     ELSE
410
        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Stepsize :",NormStep,SThresh,"NO"
411
     END IF
412
     IF (FiniG) THEN
413
        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Grad norm :",NormGrad,GThresh,"YES"
414
     ELSE
415
        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Grad norm :",NormGrad,GThresh,"NO"
416
     END IF
417

    
418
     if (DynMaxStep.AND.(IOpt>1)) THEN
419
        If (E<EOld) Then
420
           MaxStep=min(1.1*MaxStep,2.*SMax)
421
        ELSE
422
           MaxStep=max(0.8*MaxStep,SMax/2.)
423
        END IF
424
        if (debug) WRITE(*,*) "Dynamically updated  Maximum step in Opt_Geom:",MaxStep,MaxStep/SMax
425

    
426
     END IF
427

    
428
     Call Check_step(Coord,Nat,NCoord,GeomOld,GradTmp)
429

    
430
! We add the step to the old geom
431
     Geom=GeomOld+GradTmp
432

    
433
     EOld=E
434

    
435
     IF (debug) THEN
436
        WRITE(*,*) "L235, DBG OptGeom, IGeom=",IGeom," IOpt=", IOpt
437
        SELECT CASE (COORD)
438
        CASE ('ZMAT','BAKER')
439
           !WRITE(*,'(12(1X,F8.3))') Geom(1:min(12,NCoord))
440
        CASE('CART','HYBRID')
441
           DO Iat=1,Nat
442
              ! PFL 30 Apr 2008 ->
443
              ! Deleted ref to orderinv that is created only for coord in zmat, baker, hybrid.
444
              !                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
445
              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), &
446
                   Geom(3*Iat-2:3*Iat)
447
              ! <- PFL 30 Apr 2008
448
           END DO
449
        CASE('MIXED')
450
           DO Iat=1,NCart
451
              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
452
                   Geom(3*Iat-2:3*Iat)
453
           END DO
454

    
455
           SELECT CASE (NCart)
456
           CASE(1)
457
              if (NAt.GE.2) THEN
458
                 WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
459
                      IndZmat(2,2),Geom(4)
460
                 IBEG=5
461
              END IF
462
              IF (NAT.GE.3) THEN
463
                 WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
464
                      IndZmat(3,2),Geom(5),IndZmat(3,3),Geom(6)
465
                 IBeg=7
466
              END IF
467
           CASE(2)
468
              WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))')  Nom(Atome(OrderInv(indzmat(3,1)))), &
469
                   IndZmat(3,2),Geom(7),IndZmat(3,3),Geom(8)
470
              IBeg=9
471
           CASE DEFAULT
472
              IBeg=3*NCart+1
473
           END SELECT
474

    
475
           DO Iat=max(4,NCart),Nat
476
              WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
477
                   IndZmat(iat,1),Geom(IBeg),IndZmat(3,2),Geom(IBeg+1), &
478
                   IndZmat(3,3),Geom(IBeg+2)*180./pi
479
              IBeg=IBeg+3
480
           END DO
481

    
482
        CASE DEFAULT
483
           WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in OptGeom. L294."
484
           STOP
485
        END SELECT
486
     END IF ! matches IF (debug) THEN
487

    
488
  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
489

    
490
  IF (Fini) THEN
491
     WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles"
492
  ELSE
493
     WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom," *NOT* converged in ",Iopt," cycles"
494
  END IF
495
  Call Header(Line)
496
  WRITE(IoOut,*) "Last Energy E=",E
497

    
498
!  WRITE(*,*) "Initial Geometry:"
499
!  DO I=1, Nat
500
!     WRITE(*,'(3(1X,F8.4))')  XyzGeomI(IGeom,:,I)
501
!  END DO
502
  WRITE(IoOut,*) "Final Geometry:"
503

    
504
     Line=""
505
     Call PrintGeom(Line,Nat,NCoord,Geom,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)
506
     if (COORD/="CART") THEN
507
        WRITE(IoOut,*) "Last Geometry in Cartesian coordinates"
508
        Call PrintGeom(Line,Nat,NCoord,GeomCart,"CART",IOOUT,Atome,Order,OrderInv,IndZmat)
509
     END IF
510

    
511

    
512
!  DO I=1, Nat
513
!     WRITE(*,'(3(1X,F8.4))') GeomCart(I,:)
514
     !IF (I .GT. 1) Then
515
     !       WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)-GeomCart(I-1,1))*(GeomCart(I,1)-GeomCart(I-1,1)))&
516
     !        + ((GeomCart(I,2)-GeomCart(I-1,2))*(GeomCart(I,2)-GeomCart(I-1,2)))&
517
     !        + ((GeomCart(I,3)-GeomCart(I-1,3))*(GeomCart(I,3)-GeomCart(I-1,3))))
518
     !END IF
519
!  END DO
520

    
521
  IF (COORD .EQ. "BAKER") Then
522
     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
523
     ScanCoord=>Coordinate
524
     DO WHILE (Associated(ScanCoord%next))
525
        I=I+1
526
        SELECT CASE (ScanCoord%Type)
527
        CASE ('BOND')
528
           WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive_t(I)
529
        CASE ('ANGLE')
530
           WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1, &
531
                ScanCoord%At2, ScanCoord%At3,Xprimitive_t(I)*180./Pi              
532
        CASE ('DIHEDRAL')
533
           WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,&
534
                ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive_t(I)*180./Pi
535
        END SELECT
536
        ScanCoord => ScanCoord%next
537
     END DO ! matches DO WHILE (Associated(ScanCoord%next))
538
  END IF
539

    
540
  DEALLOCATE(GeomCart)
541
  DEALLOCATE(GradTmp2,GeomTmp2,GradTmp,Step)
542
  DEALLOCATE(GradOld,GeomOld)
543
  DEALLOCATE(Hess_local)
544
  DEALLOCATE(GeomCart_old)
545
  DEALLOCATE(NewGeom,NewGradTmp)
546
  DEALLOCATE(Hess_local_inv)
547

    
548
  if (debug) Call Header("Geom Optimization Over")
549

    
550
END SUBROUTINE Opt_geom