Révision 9 src/Opt_Geom.f90

Opt_Geom.f90 (revision 9)
73 73

  
74 74

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

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

  
98 100
  debug=valid('OptGeom')
99 101

  
......
104 106

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

  
107
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord))
109
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord),Step(NCoord))
108 110
  ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))
109 111
  ALLOCATE(GeomRef(NCoord))
110 112
  ALLOCATE(Hess_local(NCoord,NCoord))
......
144 146
     IF (Step_method .EQ. "RFO") Then
145 147
        Hess_local=0.
146 148
        DO I=1,NCoord
147
           Hess_local(I,I)=0.5d0
149
           Hess_local(I,I)=1.d0
148 150
        END DO
149 151
     END IF
150 152

  
......
211 213
     END IF
212 214

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

  
218 220
     IF (Step_method .EQ. "RFO") Then
219 221
        Hess_local=0.
......
243 245
  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
244 246
     IOpt=IOpt+1
245 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"
246 259
     ! Calculate the  energy and gradient
247 260
     IF (debug) THEN 
248 261
        WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
......
266 279
     ELSE
267 280
        Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart)
268 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 !!!!!!!!!!!!!!!!!!!!
269 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

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

  
272 298
     IF (debug) THEN 
......
329 355
     ! GradTmp becomes Step in Step_RFO_all.
330 356
     SELECT CASE (Step_method)
331 357
     CASE ('RFO')
332
        Call Step_RFO_all(NCoord,GradTmp,IGeom,Geom,GradTmp,Hess_local,ZeroVec)
358
        Call Step_RFO_all(NCoord,Step,IGeom,Geom,GradTmp,Hess_local,ZeroVec)
359
        GradTmp=Step
333 360
     CASE ('GDIIS')
334 361
        Call Step_DIIS(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)
335 362
        ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
......
367 394
     NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
368 395
     if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh
369 396
     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))
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
374 402
     if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep
375
     GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all.
403
     GradTmp=GradTmp*FactStep ! GradTmp is step  here, from Step_RFO_all.
376 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
377 417

  
378 418
     if (DynMaxStep.AND.(IOpt>1)) THEN
379 419
        If (E<EOld) Then
......
387 427

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

  
390
     Geom=GeomOld+GradTmp ! GradTmp is step size here, from Step_RFO_all.
430
! We add the step to the old geom
431
     Geom=GeomOld+GradTmp
391 432

  
392 433
     EOld=E
393 434

  
......
446 487

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

  
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 490
  IF (Fini) THEN
457
     WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles, Energy =",E
491
     WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles"
458 492
  ELSE
459
     WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom," *NOT* converged in ",Iopt," cycles, Last Energy = ",E
493
     WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom," *NOT* converged in ",Iopt," cycles"
460 494
  END IF
495
  Call Header(Line)
496
  WRITE(IoOut,*) "Last Energy E=",E
461 497

  
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,:)
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,:)
469 514
     !IF (I .GT. 1) Then
470 515
     !       WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)-GeomCart(I-1,1))*(GeomCart(I,1)-GeomCart(I-1,1)))&
471 516
     !        + ((GeomCart(I,2)-GeomCart(I-1,2))*(GeomCart(I,2)-GeomCart(I-1,2)))&
472 517
     !        + ((GeomCart(I,3)-GeomCart(I-1,3))*(GeomCart(I,3)-GeomCart(I-1,3))))
473 518
     !END IF
474
  END DO
519
!  END DO
475 520

  
476 521
  IF (COORD .EQ. "BAKER") Then
477 522
     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
......
493 538
  END IF
494 539

  
495 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)
496 547

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

  

Formats disponibles : Unified diff