Revision 5 src/Opt_Geom.f90

Opt_Geom.f90 (revision 5)
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.
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.
3 5

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

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

  
11 15
  IMPLICIT NONE
12 16

  
13 17
  INTEGER(KINT), INTENT(IN) :: Nat,NCoord,IGeom
......
17 21
  REAL(KREAL), INTENT(OUT) :: E
18 22
  LOGICAL, INTENT(IN) :: Flag_Opt_Geom
19 23

  
20
! As this subroutine is here for debugging, debug=T !
21
  LOGICAL  :: debug
22
  LOGICAL :: Fini
23
  LOGICAL, SAVE :: FRST=.TRUE.
24 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 25
  INTERFACE
42 26
     function valid(string) result (isValid)
43 27
       CHARACTER(*), intent(in) :: string
......
85 69

  
86 70
  END INTERFACE
87 71

  
72

  
73
  LOGICAL  :: debug
74
  LOGICAL :: Fini
75
  LOGICAL, SAVE :: FRST=.TRUE.
76

  
77
! Variables
78
  INTEGER(KINT) :: IOpt, I,J,Iat, IBEG
79
  REAL(KREAL), ALLOCATABLE :: GradTmp(:),ZeroVec(:) ! 3*Nat or 3*Nat-6
80
  REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat-6)
81
  REAL(KREAL), ALLOCATABLE :: GeomRef(:) ! (3*Nat or 3*Nat-6)
82
  REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! 3*Nat or 3*Nat-6
83
  REAL(KREAL), ALLOCATABLE :: Hess_local(:,:)
84
  REAL(KREAL), ALLOCATABLE :: GeomCart(:,:) !(Nat,3)
85
  REAL(KREAL), ALLOCATABLE :: GeomCart_old(:,:) !(Nat,3)
86
  REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim)
87
  REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat-6))
88
  REAL(KREAL), ALLOCATABLE :: NewGeom(:) !NCoord=3*Nat-6
89
  REAL(KREAL), ALLOCATABLE :: NewGradTmp(:)
90
  REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis
91
  REAL(KREAL) :: NormStep, FactStep, HP, MaxStep
92
  REAL(KREAL) :: Eold  
93

  
88 94
  debug=valid('OptGeom')
89 95

  
90 96
    E=0.
97
    Eold=0.
98
    MaxStep=SMax
91 99

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

  
......
95 103
     ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))
96 104
     ALLOCATE(GeomRef(NCoord))
97 105
     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))
106
     ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3))
107
     ALLOCATE(NewGeom(NCoord))
108
     ALLOCATE(NewGradTmp(NCoord))
109
     ALLOCATE(Hess_local_inv(NCoord,NCoord))
102 110
        
103 111
     IOpt=0
104 112
     ZeroVec=0.d0
......
212 220
        CASE DEFAULT
213 221
        WRITE (*,*) "Coord=",TRIM(COORD)," not recognized, opt_Geom.f90, L164. Stop"
214 222
        STOP
215
     END SELECT       
223
     END SELECT
216 224

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

  
229
     IF (FPrintGeom) THEN
230
        OPEN(IOGeom,File=GeomFile)
231
     END IF
232

  
221 233
     Fini=.FALSE.
222 234
     DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
223 235
        IOpt=IOpt+1
......
251 263
           WRITE(*,*) "Energy and Coordinates:"
252 264
           WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
253 265
           WRITE(*,*) "Geom Old:"
254
           WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13))
255
           END IF
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
256 270

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

  
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
258 287
        ! We have enough geometries and gradient to update the hessian (or its 
259 288
        ! inverse)
260 289
        GradTmp2=GradTmp-GradOld
......
313 342
        !   If (IOpt.LT.5) GradTmp=(IOpt*GradTmp+(5-IOpt)*0.01*Grad(IGeom,:))/5.d0
314 343
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
315 344
              WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh
316
        FactStep=min(1.d0,SMax/NormStep)
345
        FactStep=min(1.d0,MaxStep/NormStep)
317 346
        Fini=(NormStep.LE.SThresh)
318 347
        NormStep=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients.
319 348
              WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormStep (Gradient)=', NormStep, 'Gradient Threshold=', GThresh
......
321 350
        if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep
322 351
        GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all.
323 352

  
353

  
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

  
362
        END IF
363

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

  
368
        EOld=E
369

  
328 370
        IF (debug) THEN
329 371
           WRITE(*,*) "L235, DBG OptGeom, IGeom=",IGeom," IOpt=", IOpt
330 372
           SELECT CASE (COORD)
......
390 432
     IF (Fini) THEN
391 433
        WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles, Energy =",E
392 434
     ELSE
393
        WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom,"*NOT* converged in ",Iopt," cycles, Last Energy = ",E
435
        WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom," *NOT* converged in ",Iopt," cycles, Last Energy = ",E
394 436
     END IF
395 437
                
396 438
        WRITE(*,*) "Initial Geometry:"

Also available in: Unified diff