Statistiques
| Révision :

root / src / Step_RFO_all.f90 @ 8

Historique | Voir | Annoter | Télécharger (8,38 ko)

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=1e-20, zero=0., dele3=1e-6, eps=1e-12
71
  REAL(KREAL), PARAMETER :: crit=1e-8
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 Ncoord-1 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=Ncoord-1
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/= NCoord-1) 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