root / src / Step_RFO_all.f90 @ 8
History  View  Annotate  Download (8.4 kB)
1 
SUBROUTINE Step_RFO_all(NCoord,Step,IGeom,Geom,Grad,Hess,Tangent) 

2  
3 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
4 
! 
5 
!This subroutine calculates new coordinates given a Hessian, the forces 
6 
! and the old coordinates 
7 
! 
8 
!!!!!!!!!! 
9 
! v1.0 
10 
! Uses only basic RFO step 
11 
! 
12 
!!!!!!!!!! 
13 
! v 2.0 
14 
! We add the test of ADF here, so that it should be more efficient. 
15 
! This version uses HInv to decide wether it should work with the Hessian 
16 
! or its inverse. 
17 
! 
18 
! v 3.0 
19 
! Uses DIIS. Done by P. Dayal. 
20 
! 
21 
! v3.1 
22 
! Back to origins ! We remove the ADF test here to see what we loose... 
23 
! untill of course its replacement by something more efficient... and free :) 
24 
!!!!!!!!!!!!!!!!! 
25 
! 
26 
! Input: 
27 
!  NCoord: INTEGER(KINT), Number of degrees of freedom 
28 
!  IGeom : INTEGER(KINT), Index of the geometry along the path 
29 
!  Geom : REAL(KREAL)(NCOORD), Geometry expressed in full coordinates 
30 
!  Hess : REAL(KREAL)(Ncoord,NCoord), Hessian at the current geometry 
31 
! 
32 
! 
33 
! Output: 
34 
!  Step : REAL(KREAL)(Ncoord), Calculated step 
35 
! 
36 
! In/Out: 
37 
!  Tangent: REAL(KREAL)(Ncoord), Tangent to the path, at the current point 
38 
! if Tangent=0, then we are optimizing a minimum, tangent is not considered. 
39 
! Input: Expressed in full coordinates. 
40 
! Output: Replaced by the tangent expressed in Free vectors 
41 
! (Changes by PFL 3 Jan 2008) 
42 
!  Grad: READ(KREAL) (NCoord), Gradient at the current point 
43 
! Input: Expressed in full coordinates. 
44 
! PFL 18 Jan 2008: Grad is no longer changed on output > 
45 
! Output: Replaced by the gradient expressed in Free vector, in the subspace orthogonal to the tangent. 
46 
! < PFL 18 Jan 2008: Grad is no longer changed on output 
47 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
48  
49 
use Io_module, only : IoOut 
50 

51 
use Path_module, only : NGeomF, Hinv,Coord,Vfree, NGintMax 
52 
! VFree(Ncoord,Ncoord) 
53  
54 
IMPLICIT NONE 
55  
56 
INTEGER, PARAMETER :: KINT=KIND(1) 
57 
INTEGER, PARAMETER :: KREAL=KIND(1.D0) 
58  
59 
INTEGER(KINT), INTENT(IN) :: NCoord,IGeom 
60 
REAL(KREAL), INTENT(IN) :: Hess(NCoord,NCoord) 
61 
REAL(KREAL), INTENT(IN) :: Grad(NCoord),Geom(NCoord) 
62 
REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord) ! Destroyed on output 
63 
REAL(KREAL), INTENT(OUT) :: Step(NCoord) 
64  
65 
INTEGER(KINT) :: NFree 
66 
REAL(KREAL), ALLOCATABLE :: eigval(:), eigvec(:,:), eigvli(:) 
67 
REAL(KREAL), ALLOCATABLE :: Tanf(:) 
68 
REAL(KREAL), ALLOCATABLE :: Grad_f(:),Step_f(:) ! NFree 
69  
70 
REAL(KREAL), PARAMETER :: tiny=1e20, zero=0., dele3=1e6, eps=1e12 
71 
REAL(KREAL), PARAMETER :: crit=1e8 
72 
INTEGER(KINT) :: i, j, idx, k, isch 
73 
REAL(KREAL) :: grd, eval, norm 
74 
REAL(KREAL) :: dx 
75 
CHARACTER(120) :: fmt2 
76  
77 
LOGICAL :: Debug 
78  
79 
REAL(KREAL), ALLOCATABLE :: HFree(:,:) ! NFree,NFree 
80 
REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree 
81  
82 
! LOGICAL, SAVE :: LFirst=.TRUE. 
83  
84 
INTERFACE 
85 
function valid(string) result (isValid) 
86 
CHARACTER(*), intent(in) :: string 
87 
logical :: isValid 
88 
END function VALID 
89 
END INTERFACE 
90  
91 
debug=valid("step_rfo_all") 
92 
if (debug) WRITE(*,*) "=========================== Entering Step_RFO_All ====================" 
93 
if (debug) WRITE(*,*) "=========================== IGeom=",IGeom," ====================" 
94  
95 
! PFL 22.Nov.2007 > 
96 
! Nfree is equal to Ncoord1 only when optimizing in the plane 
97 
! orthogonal to the tangent. 
98 
! In the case of a regular opt, NFree=Ncoord. 
99 
! Next line has been moved later 
100 
! NFree=Ncoord1 
101 
! < PFL 22.Nov.2007 
102  
103 
IF (debug) THEN 
104 
WRITE(*,*) 'DBG Step_RFO_All Hess' 
105 
WRITE(*,*) 'DBG Step_RFO_All NCoord=',NCoord 
106  
107 
Idx=min(12,NCoord) 
108 
DO I=1,NCoord 
109 
WRITE(*,'(12(1X,F10.4))') Hess(I,1:Idx) 
110 
END DO 
111 
WRITE(*,*) 'DBG Step_RFO_All Grad' 
112 
WRITE(*,'(12(1X,F10.4))') Grad(1:Idx) 
113 
WRITE(*,*) 'DBG Step_RFO_All Coord' 
114 
WRITE(*,'(12(1X,F10.4))') Geom(1:Idx) 
115 
WRITE(*,*) 'DBG Step_RFO_All Tangent' 
116 
WRITE(*,'(12(1X,F10.4))') tangent(1:Idx) 
117  
118 
END IF 
119  
120 
Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord) 
121 
! we orthogonalize Vfree to the tangent vector of this geom 
122 
! 2007/05/29 Only if Tangent/=0.d0 
123 
Norm=sqrt(dot_product(Tangent,Tangent)) 
124 
IF (Norm.GT.eps) THEN 
125 
ALLOCATE(Tanf(NCoord)) 
126  
127 
! We normalize Tangent 
128 
Tangent=Tangent/Norm 
129  
130 
! We convert Tangent into vfree only displacements 
131 
! This is useless for now (2007.Apr.23) as vfree=Id matrix 
132 
! but it will be usefull as soon as we introduce Constraints 
133 
DO I=1,NCoord 
134 
Tanf(I)=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent) 
135 
END DO 
136 
Tangent=0.d0 
137 
DO I=1,NCoord 
138 
Tangent=Tangent+Tanf(i)*vfree(:,I) 
139 
END DO 
140  
141 
! first we substract Tangent from vfree 
142 
DO I=1,NCoord 
143 
Norm=dot_product(reshape(vfree(:,i),(/NCoord/)),Tangent) 
144 
Vfree(:,I)=vfree(:,i)Norm*Tangent 
145 
END DO 
146  
147 
Idx=0 
148 
! Schmidt orthogonalization of the vfree vectors 
149 
DO I=1,NCoord 
150 
! We substract the first vectors 
151 
! we do it twice as the Schmidt procedure is not numerically stable 
152 
DO isch=1,2 
153 
DO J=1,Idx 
154 
Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),reshape(vfree(:,J),(/NCoord/))) 
155 
vfree(:,I)=vfree(:,I)Norm*vfree(:,J) 
156 
END DO 
157 
END DO 
158 
Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),reshape(vfree(:,I),(/NCoord/))) 
159 
IF (Norm.GE.crit) THEN 
160 
idx=idx+1 
161 
vfree(:,idx)=vfree(:,i)/sqrt(Norm) 
162 
END IF 
163 
END DO 
164 

165 
If (Idx/= NCoord1) THEN 
166 
WRITE(*,*) "Pb in orthogonalizing vfree to tangent for geom",IGeom 
167 
WRITE(Ioout,*) "Pb in orthogonalizing vfree to tangent for geom",IGeom 
168 
STOP 
169 
END IF 
170 

171 
if (debug) WRITE(*,*) "Deallocating Tanf" 
172 
DEALLOCATE(Tanf) 
173  
174 
NFree=Idx 
175  
176 
ELSE ! Tangent =0, matches IF (Norm.GT.eps) THEN 
177 
if (debug) WRITE(*,*) "Tangent=0, using full displacement" 
178 
NFree=NCoord 
179 
END IF !IF (Norm.GT.eps) THEN 
180 

181 
if (debug) WRITE(*,*) 'DBG Step_RFO_All NFree=',NFree 
182 

183  
184 
! We now calculate the new step 
185 
! we project the hessian onto the free vectors 
186 
ALLOCATE(hfree(NFree,NFree),htmp(NCoord,NFree),Grad_f(NFree),Step_f(NFree)) 
187 
DO J=1,NFree 
188 
DO I=1,NCoord 
189 
Htmp(I,J)=0.d0 
190 
DO K=1,NCoord 
191 
Htmp(I,J)=htmp(I,J)+Hess(I,K)*vfree(K,J) 
192 
END DO 
193 
END DO 
194 
END DO 
195 
DO J=1,NFree 
196 
DO I=1,NFree 
197 
HFree(i,j)=0.d0 
198 
DO K=1,NCoord 
199 
HFree(i,j)=hfree(i,j)+vfree(k,i)*htmp(k,j) 
200 
END DO 
201 
END DO 
202 
END DO 
203  
204 
! We diagonalize the hessian or its inverse 
205 
ALLOCATE(eigval(NFree), eigvec(NFree,NFree), eigvli(NFree)) 
206  
207 
call Jacobi(Hfree,NFree,eigvli,eigvec,NFree) 
208  
209 
if (debug) THEN 
210 
WRITE(*,*) 'Eigensystem, Step_RFO_all.f90, L225' 
211 
fmt2='(I5,":",F8.3," ", (1X,F10.4))' 
212 
WRITE(fmt2(19:20),'(I2)') min(NFree,99) 
213 
DO I=1,NFree 
214 
WRITE(*,fmt2) I,Eigvli(I),eigvec(:,I) 
215 
END DO 
216 
END IF 
217  
218 
! We inverse the eigenvalues if Hess corresponds to Hessian inverse 
219 
IF (Hinv) THEN 
220 
do J=1,NFree 
221 
if (abs(eigvli(j)) <= tiny) eigvli(j) = zero 
222 
eigval(j) = sign (eigvli(j)**2 / (dele3 + abs(eigvli(j))**3), eigvli(j)) 
223 
if (abs(eigval(j)) <= tiny) eigval(j) = zero 
224 
end do 
225 
ELSE 
226 
eigval=eigvli 
227 
END IF 
228  
229 
call trie(NFree,eigval,eigvec,NFree) 
230  
231 
DO I=1,NFree 
232 
Grad_f(I)=dot_product(reshape(vfree(:,I),(/NCoord/)),grad) 
233 
END DO 
234  
235 
! We now construct the step 
236  
237 
Step_f=0.d0 
238 
DO I=1,NFree 
239 
grd= dot_product(grad_f,reshape(eigvec(:,i),(/NFree/))) 
240 
eval = 0.5 * (abs(eigval(i)) + sqrt (eigval(i)**2 + 4.*grd**2)) 
241 
dx = (grd) / eval 
242 
step_f = step_f + dx * eigvec(:,i) 
243 
if (debug) write (*,*) i,eigval(i),eval,grd,dx 
244 
END DO 
245  
246 
Step=0.d0 
247 
! PFL 04 Apr 2008 > 
248 
! Grad is once again modified 
249 
! Grad=0.d0 
250 
DO I=1,NFree 
251 
Step=Step+ step_f(i)*vfree(:,i) 
252 
! Grad=Grad+ Grad_f(i)*vfree(:,i) 
253 
END DO 
254 
! < PFL 04 Apr 2008 
255 

256 
if (debug) WRITE(*,*) "Deallocataing Eigval" 
257 
DEALLOCATE(eigval) 
258 
if (debug) WRITE(*,*) "Deallocataing Eigvec" 
259 
DEALLOCATE(eigvec) 
260 
if (debug) WRITE(*,*) "Deallocataing Eigvli" 
261 
DEALLOCATE(eigvli) 
262 
if (debug) WRITE(*,*) "Deallocataing hfree" 
263 
DEALLOCATE(hfree) 
264 
if (debug) WRITE(*,*) "Deallocataing htmp" 
265 
DEALLOCATE(htmp) 
266 
if (debug) WRITE(*,*) "Deallocataing grad_f" 
267 
DEALLOCATE(grad_f) 
268 
if (debug) WRITE(*,*) "Deallocataing step_f" 
269 
DEALLOCATE(step_f) 
270  
271 
if (debug) THEN 
272 
WRITE(*,*) 'DBG Step_RFO_All Step' 
273 
WRITE(*,'(12(1X,F8.4))') Step(1:Idx) 
274 
END IF 
275  
276 
if (debug) WRITE(*,*) "=========================== Step_RFO_All Over ====================" 
277  
278 
END SUBROUTINE STEP_RFO_ALL 
279  
280 