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*Nat6 

28 
REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat6) 

29 
REAL(KREAL), ALLOCATABLE :: GeomRef(:) ! (3*Nat or 3*Nat6) 

30 
REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! 3*Nat or 3*Nat6 

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*Nat6)) 

36 
REAL(KREAL), ALLOCATABLE :: NewGeom(:) !NCoord=3*Nat6 

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*Nat6 

80 
REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat6) 

81 
REAL(KREAL), ALLOCATABLE :: GeomRef(:) ! (3*Nat or 3*Nat6) 

82 
REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! 3*Nat or 3*Nat6 

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*Nat6)) 

88 
REAL(KREAL), ALLOCATABLE :: NewGeom(:) !NCoord=3*Nat6 

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*Iat2: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*I2: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=GradTmpGradOld 
...  ...  
313  342 
! If (IOpt.LT.5) GradTmp=(IOpt*GradTmp+(5IOpt)*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