Statistics
| Revision:

root / src / Step_RFO_all.f90 @ 7

History | View | Annotate | Download (8.4 kB)

1 1 pfleura2
SUBROUTINE Step_RFO_all(NCoord,Step,IGeom,Geom,Grad,Hess,Tangent)
2 1 pfleura2
3 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4 1 pfleura2
!
5 1 pfleura2
!This subroutine calculates new coordinates given a Hessian, the forces
6 1 pfleura2
! and the old coordinates
7 1 pfleura2
!
8 1 pfleura2
!!!!!!!!!!
9 1 pfleura2
! v1.0
10 1 pfleura2
! Uses only basic RFO step
11 1 pfleura2
!
12 1 pfleura2
!!!!!!!!!!
13 1 pfleura2
! v 2.0
14 1 pfleura2
! We add the test of ADF here, so that it should be more efficient.
15 1 pfleura2
! This version uses HInv to decide wether it should work with the Hessian
16 1 pfleura2
! or its inverse.
17 1 pfleura2
!
18 1 pfleura2
! v 3.0
19 1 pfleura2
! Uses DIIS. Done by P. Dayal.
20 1 pfleura2
!
21 1 pfleura2
! v3.1
22 1 pfleura2
! Back to origins ! We remove the ADF test here to see what we loose...
23 1 pfleura2
! untill of course its replacement by something more efficient... and free :)
24 1 pfleura2
!!!!!!!!!!!!!!!!!
25 1 pfleura2
!
26 1 pfleura2
! Input:
27 1 pfleura2
!   - NCoord: INTEGER(KINT), Number of degrees of freedom
28 1 pfleura2
!   - IGeom : INTEGER(KINT), Index of the geometry along the path
29 1 pfleura2
!   - Geom  : REAL(KREAL)(NCOORD), Geometry expressed in full coordinates
30 1 pfleura2
!   - Hess  : REAL(KREAL)(Ncoord,NCoord), Hessian at the current geometry
31 1 pfleura2
!
32 1 pfleura2
!
33 1 pfleura2
! Output:
34 1 pfleura2
!   - Step  : REAL(KREAL)(Ncoord), Calculated step
35 1 pfleura2
!
36 1 pfleura2
! In/Out:
37 1 pfleura2
!   - Tangent: REAL(KREAL)(Ncoord), Tangent to the path, at the current point
38 1 pfleura2
!             if Tangent=0, then we are optimizing a minimum,  tangent is not considered.
39 1 pfleura2
!               Input: Expressed in full coordinates.
40 1 pfleura2
!              Output: Replaced by the tangent expressed in Free vectors
41 1 pfleura2
! (Changes by PFL 3 Jan 2008)
42 1 pfleura2
!   -  Grad: READ(KREAL) (NCoord), Gradient at the current point
43 1 pfleura2
!               Input: Expressed in full coordinates.
44 1 pfleura2
! PFL 18 Jan 2008: Grad is no longer changed on output ->
45 1 pfleura2
!              Output: Replaced by the gradient expressed in Free vector, in the subspace orthogonal to the tangent.
46 1 pfleura2
! <- PFL 18 Jan 2008: Grad is no longer changed on output
47 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
48 1 pfleura2
49 1 pfleura2
  use Io_module, only : IoOut
50 1 pfleura2
51 1 pfleura2
  use Path_module, only : NGeomF, Hinv,Coord,Vfree, NGintMax
52 1 pfleura2
  ! VFree(Ncoord,Ncoord)
53 1 pfleura2
54 1 pfleura2
  IMPLICIT NONE
55 1 pfleura2
56 1 pfleura2
  INTEGER, PARAMETER :: KINT=KIND(1)
57 1 pfleura2
  INTEGER, PARAMETER :: KREAL=KIND(1.D0)
58 1 pfleura2
59 1 pfleura2
  INTEGER(KINT), INTENT(IN) :: NCoord,IGeom
60 1 pfleura2
  REAL(KREAL), INTENT(IN) :: Hess(NCoord,NCoord)
61 1 pfleura2
  REAL(KREAL), INTENT(IN) :: Grad(NCoord),Geom(NCoord)
62 1 pfleura2
  REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord) ! Destroyed on output
63 1 pfleura2
  REAL(KREAL), INTENT(OUT) :: Step(NCoord)
64 1 pfleura2
65 1 pfleura2
  INTEGER(KINT) :: NFree
66 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: eigval(:), eigvec(:,:), eigvli(:)
67 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: Tanf(:)
68 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: Grad_f(:),Step_f(:) ! NFree
69 1 pfleura2
70 1 pfleura2
  REAL(KREAL), PARAMETER :: tiny=1e-20, zero=0., dele3=1e-6, eps=1e-12
71 1 pfleura2
  REAL(KREAL), PARAMETER :: crit=1e-8
72 2 pfleura2
  INTEGER(KINT) :: i, j, idx, k, isch
73 2 pfleura2
  REAL(KREAL) :: grd, eval, norm
74 2 pfleura2
  REAL(KREAL) :: dx
75 2 pfleura2
  CHARACTER(120) :: fmt2
76 1 pfleura2
77 1 pfleura2
  LOGICAL :: Debug
78 1 pfleura2
79 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: HFree(:,:) ! NFree,NFree
80 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree
81 1 pfleura2
82 2 pfleura2
!  LOGICAL, SAVE :: LFirst=.TRUE.
83 1 pfleura2
84 1 pfleura2
  INTERFACE
85 1 pfleura2
     function valid(string) result (isValid)
86 1 pfleura2
       CHARACTER(*), intent(in) :: string
87 1 pfleura2
       logical                  :: isValid
88 1 pfleura2
     END function VALID
89 1 pfleura2
  END INTERFACE
90 1 pfleura2
91 1 pfleura2
  debug=valid("step_rfo_all")
92 1 pfleura2
  if (debug) WRITE(*,*) "=========================== Entering Step_RFO_All ===================="
93 1 pfleura2
  if (debug) WRITE(*,*) "=========================== IGeom=",IGeom," ===================="
94 1 pfleura2
95 1 pfleura2
! PFL 22.Nov.2007 ->
96 1 pfleura2
! Nfree is equal to Ncoord-1 only when optimizing in the plane
97 1 pfleura2
! orthogonal to the tangent.
98 1 pfleura2
! In the case of a regular opt, NFree=Ncoord.
99 1 pfleura2
! Next line has been moved later
100 1 pfleura2
!  NFree=Ncoord-1
101 1 pfleura2
! <- PFL 22.Nov.2007
102 1 pfleura2
103 1 pfleura2
  IF (debug) THEN
104 1 pfleura2
     WRITE(*,*) 'DBG Step_RFO_All Hess'
105 1 pfleura2
     WRITE(*,*) 'DBG Step_RFO_All NCoord=',NCoord
106 1 pfleura2
107 1 pfleura2
     Idx=min(12,NCoord)
108 1 pfleura2
     DO I=1,NCoord
109 1 pfleura2
        WRITE(*,'(12(1X,F10.4))') Hess(I,1:Idx)
110 1 pfleura2
     END DO
111 1 pfleura2
     WRITE(*,*) 'DBG Step_RFO_All Grad'
112 1 pfleura2
     WRITE(*,'(12(1X,F10.4))') Grad(1:Idx)
113 1 pfleura2
     WRITE(*,*) 'DBG Step_RFO_All Coord'
114 1 pfleura2
     WRITE(*,'(12(1X,F10.4))') Geom(1:Idx)
115 1 pfleura2
     WRITE(*,*) 'DBG Step_RFO_All Tangent'
116 1 pfleura2
     WRITE(*,'(12(1X,F10.4))') tangent(1:Idx)
117 1 pfleura2
118 1 pfleura2
  END IF
119 1 pfleura2
120 1 pfleura2
  Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord)
121 1 pfleura2
! we orthogonalize Vfree to the tangent vector of this geom
122 1 pfleura2
! 2007/05/29 Only if Tangent/=0.d0
123 1 pfleura2
  Norm=sqrt(dot_product(Tangent,Tangent))
124 1 pfleura2
  IF (Norm.GT.eps) THEN
125 1 pfleura2
     ALLOCATE(Tanf(NCoord))
126 1 pfleura2
127 1 pfleura2
     ! We normalize Tangent
128 1 pfleura2
     Tangent=Tangent/Norm
129 1 pfleura2
130 1 pfleura2
     ! We convert Tangent into vfree only displacements
131 1 pfleura2
     ! This is useless for now (2007.Apr.23) as vfree=Id matrix
132 1 pfleura2
     ! but it will be usefull as soon as we introduce Constraints
133 1 pfleura2
     DO I=1,NCoord
134 1 pfleura2
        Tanf(I)=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent)
135 1 pfleura2
     END DO
136 1 pfleura2
     Tangent=0.d0
137 1 pfleura2
     DO I=1,NCoord
138 1 pfleura2
        Tangent=Tangent+Tanf(i)*vfree(:,I)
139 1 pfleura2
     END DO
140 1 pfleura2
141 1 pfleura2
     ! first we substract Tangent from vfree
142 1 pfleura2
     DO I=1,NCoord
143 1 pfleura2
        Norm=dot_product(reshape(vfree(:,i),(/NCoord/)),Tangent)
144 1 pfleura2
        Vfree(:,I)=vfree(:,i)-Norm*Tangent
145 1 pfleura2
     END DO
146 1 pfleura2
147 1 pfleura2
     Idx=0.
148 1 pfleura2
     ! Schmidt orthogonalization of the vfree vectors
149 1 pfleura2
     DO I=1,NCoord
150 1 pfleura2
        ! We substract the first vectors
151 1 pfleura2
        ! we do it twice as the Schmidt procedure is not numerically stable
152 1 pfleura2
        DO isch=1,2
153 1 pfleura2
           DO J=1,Idx
154 1 pfleura2
              Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),reshape(vfree(:,J),(/NCoord/)))
155 1 pfleura2
              vfree(:,I)=vfree(:,I)-Norm*vfree(:,J)
156 1 pfleura2
           END DO
157 1 pfleura2
        END DO
158 1 pfleura2
        Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),reshape(vfree(:,I),(/NCoord/)))
159 1 pfleura2
        IF (Norm.GE.crit) THEN
160 1 pfleura2
           idx=idx+1
161 1 pfleura2
           vfree(:,idx)=vfree(:,i)/sqrt(Norm)
162 1 pfleura2
        END IF
163 1 pfleura2
     END DO
164 1 pfleura2
165 1 pfleura2
     If (Idx/= NCoord-1) THEN
166 1 pfleura2
        WRITE(*,*) "Pb in orthogonalizing vfree to tangent for geom",IGeom
167 1 pfleura2
        WRITE(Ioout,*) "Pb in orthogonalizing vfree to tangent for geom",IGeom
168 1 pfleura2
        STOP
169 1 pfleura2
     END IF
170 1 pfleura2
171 2 pfleura2
     if (debug) WRITE(*,*) "Deallocating Tanf"
172 1 pfleura2
     DEALLOCATE(Tanf)
173 1 pfleura2
174 1 pfleura2
     NFree=Idx
175 1 pfleura2
176 1 pfleura2
  ELSE ! Tangent =0, matches IF (Norm.GT.eps) THEN
177 1 pfleura2
     if (debug) WRITE(*,*) "Tangent=0, using full displacement"
178 1 pfleura2
     NFree=NCoord
179 1 pfleura2
  END IF !IF (Norm.GT.eps) THEN
180 1 pfleura2
181 1 pfleura2
  if (debug) WRITE(*,*) 'DBG Step_RFO_All NFree=',NFree
182 1 pfleura2
183 1 pfleura2
184 1 pfleura2
! We now calculate the new step
185 1 pfleura2
! we project the hessian onto the free vectors
186 1 pfleura2
     ALLOCATE(hfree(NFree,NFree),htmp(NCoord,NFree),Grad_f(NFree),Step_f(NFree))
187 1 pfleura2
     DO J=1,NFree
188 1 pfleura2
        DO I=1,NCoord
189 1 pfleura2
           Htmp(I,J)=0.d0
190 1 pfleura2
           DO K=1,NCoord
191 1 pfleura2
              Htmp(I,J)=htmp(I,J)+Hess(I,K)*vfree(K,J)
192 1 pfleura2
           END DO
193 1 pfleura2
        END DO
194 1 pfleura2
     END DO
195 1 pfleura2
     DO J=1,NFree
196 1 pfleura2
        DO I=1,NFree
197 1 pfleura2
           HFree(i,j)=0.d0
198 1 pfleura2
           DO K=1,NCoord
199 1 pfleura2
              HFree(i,j)=hfree(i,j)+vfree(k,i)*htmp(k,j)
200 1 pfleura2
           END DO
201 1 pfleura2
        END DO
202 1 pfleura2
     END DO
203 1 pfleura2
204 1 pfleura2
! We diagonalize the hessian or its inverse
205 1 pfleura2
  ALLOCATE(eigval(NFree), eigvec(NFree,NFree), eigvli(NFree))
206 1 pfleura2
207 1 pfleura2
  call Jacobi(Hfree,NFree,eigvli,eigvec,NFree)
208 1 pfleura2
209 1 pfleura2
  if (debug) THEN
210 1 pfleura2
     WRITE(*,*) 'Eigensystem, Step_RFO_all.f90, L225'
211 1 pfleura2
     fmt2='(I5,":",F8.3," -",  (1X,F10.4))'
212 1 pfleura2
     WRITE(fmt2(19:20),'(I2)') min(NFree,99)
213 1 pfleura2
     DO I=1,NFree
214 1 pfleura2
        WRITE(*,fmt2) I,Eigvli(I),eigvec(:,I)
215 1 pfleura2
     END DO
216 1 pfleura2
  END IF
217 1 pfleura2
218 1 pfleura2
  ! We inverse the eigenvalues if Hess corresponds to Hessian inverse
219 1 pfleura2
  IF (Hinv) THEN
220 1 pfleura2
     do J=1,NFree
221 1 pfleura2
        if (abs(eigvli(j)) <= tiny) eigvli(j) = zero
222 1 pfleura2
        eigval(j) = sign (eigvli(j)**2 / (dele3 + abs(eigvli(j))**3), eigvli(j))
223 1 pfleura2
        if (abs(eigval(j)) <= tiny) eigval(j) = zero
224 1 pfleura2
     end do
225 1 pfleura2
  ELSE
226 1 pfleura2
     eigval=eigvli
227 1 pfleura2
  END IF
228 1 pfleura2
229 1 pfleura2
  call trie(NFree,eigval,eigvec,NFree)
230 1 pfleura2
231 1 pfleura2
  DO I=1,NFree
232 1 pfleura2
     Grad_f(I)=dot_product(reshape(vfree(:,I),(/NCoord/)),grad)
233 1 pfleura2
  END DO
234 1 pfleura2
235 1 pfleura2
  ! We now construct the step
236 1 pfleura2
237 1 pfleura2
  Step_f=0.d0
238 1 pfleura2
  DO I=1,NFree
239 1 pfleura2
     grd= dot_product(grad_f,reshape(eigvec(:,i),(/NFree/)))
240 1 pfleura2
     eval = 0.5 * (abs(eigval(i)) + sqrt (eigval(i)**2 + 4.*grd**2))
241 1 pfleura2
     dx   = (-grd) / eval
242 1 pfleura2
     step_f = step_f + dx * eigvec(:,i)
243 1 pfleura2
     if (debug) write (*,*) i,eigval(i),eval,grd,dx
244 1 pfleura2
  END DO
245 1 pfleura2
246 1 pfleura2
  Step=0.d0
247 1 pfleura2
! PFL 04 Apr 2008 ->
248 1 pfleura2
! Grad is once again modified
249 1 pfleura2
!  Grad=0.d0
250 1 pfleura2
  DO I=1,NFree
251 1 pfleura2
     Step=Step+ step_f(i)*vfree(:,i)
252 1 pfleura2
!     Grad=Grad+ Grad_f(i)*vfree(:,i)
253 1 pfleura2
  END DO
254 1 pfleura2
! <- PFL 04 Apr 2008
255 1 pfleura2
256 2 pfleura2
  if (debug) WRITE(*,*) "Deallocataing Eigval"
257 2 pfleura2
  DEALLOCATE(eigval)
258 2 pfleura2
  if (debug) WRITE(*,*) "Deallocataing Eigvec"
259 2 pfleura2
  DEALLOCATE(eigvec)
260 2 pfleura2
  if (debug) WRITE(*,*) "Deallocataing Eigvli"
261 2 pfleura2
  DEALLOCATE(eigvli)
262 2 pfleura2
  if (debug) WRITE(*,*) "Deallocataing hfree"
263 2 pfleura2
  DEALLOCATE(hfree)
264 2 pfleura2
  if (debug) WRITE(*,*) "Deallocataing htmp"
265 2 pfleura2
  DEALLOCATE(htmp)
266 2 pfleura2
  if (debug) WRITE(*,*) "Deallocataing grad_f"
267 2 pfleura2
  DEALLOCATE(grad_f)
268 2 pfleura2
  if (debug) WRITE(*,*) "Deallocataing step_f"
269 2 pfleura2
  DEALLOCATE(step_f)
270 1 pfleura2
271 1 pfleura2
  if (debug) THEN
272 1 pfleura2
     WRITE(*,*) 'DBG Step_RFO_All Step'
273 1 pfleura2
     WRITE(*,'(12(1X,F8.4))') Step(1:Idx)
274 1 pfleura2
  END IF
275 1 pfleura2
276 1 pfleura2
  if (debug) WRITE(*,*) "=========================== Step_RFO_All Over ===================="
277 1 pfleura2
278 1 pfleura2
END SUBROUTINE STEP_RFO_ALL
279 1 pfleura2