Revision 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*Nat6 
82 
REAL(KREAL), ALLOCATABLE :: Step(:) ! 3*Nat or 3*Nat6 

82  83 
REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat6) 
83  84 
REAL(KREAL), ALLOCATABLE :: GeomRef(:) ! (3*Nat or 3*Nat6) 
84  85 
REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! 3*Nat or 3*Nat6 
...  ...  
90  91 
REAL(KREAL), ALLOCATABLE :: NewGeom(:) !NCoord=3*Nat6 
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(I1,1))*(GeomCart(I,1)GeomCart(I1,1)))& 
471  516 
! + ((GeomCart(I,2)GeomCart(I1,2))*(GeomCart(I,2)GeomCart(I1,2)))& 
472  517 
! + ((GeomCart(I,3)GeomCart(I1,3))*(GeomCart(I,3)GeomCart(I1,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 
Also available in: Unified diff