Statistics
| Revision:

## root / src / Step_RFO_all.f90 @ 8

 1 SUBROUTINE Step_RFO_all(NCoord,Step,IGeom,Geom,Grad,Hess,Tangent)  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !  !This subroutine calculates new coordinates given a Hessian, the forces  ! and the old coordinates  !  !!!!!!!!!!  ! v1.0  ! Uses only basic RFO step  !  !!!!!!!!!!  ! v 2.0  ! We add the test of ADF here, so that it should be more efficient.  ! This version uses HInv to decide wether it should work with the Hessian  ! or its inverse.  !  ! v 3.0  ! Uses DIIS. Done by P. Dayal.  !  ! v3.1  ! Back to origins ! We remove the ADF test here to see what we loose...  ! untill of course its replacement by something more efficient... and free :)  !!!!!!!!!!!!!!!!!  !  ! Input:  ! - NCoord: INTEGER(KINT), Number of degrees of freedom  ! - IGeom : INTEGER(KINT), Index of the geometry along the path  ! - Geom : REAL(KREAL)(NCOORD), Geometry expressed in full coordinates  ! - Hess : REAL(KREAL)(Ncoord,NCoord), Hessian at the current geometry  !  !  ! Output:  ! - Step : REAL(KREAL)(Ncoord), Calculated step  !  ! In/Out:  ! - Tangent: REAL(KREAL)(Ncoord), Tangent to the path, at the current point  ! if Tangent=0, then we are optimizing a minimum, tangent is not considered.  ! Input: Expressed in full coordinates.  ! Output: Replaced by the tangent expressed in Free vectors  ! (Changes by PFL 3 Jan 2008)  ! - Grad: READ(KREAL) (NCoord), Gradient at the current point  ! Input: Expressed in full coordinates.  ! PFL 18 Jan 2008: Grad is no longer changed on output ->  ! Output: Replaced by the gradient expressed in Free vector, in the subspace orthogonal to the tangent.  ! <- PFL 18 Jan 2008: Grad is no longer changed on output  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   use Io_module, only : IoOut     use Path_module, only : NGeomF, Hinv,Coord,Vfree, NGintMax   ! VFree(Ncoord,Ncoord)   IMPLICIT NONE   INTEGER, PARAMETER :: KINT=KIND(1)   INTEGER, PARAMETER :: KREAL=KIND(1.D0)   INTEGER(KINT), INTENT(IN) :: NCoord,IGeom   REAL(KREAL), INTENT(IN) :: Hess(NCoord,NCoord)   REAL(KREAL), INTENT(IN) :: Grad(NCoord),Geom(NCoord)   REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord) ! Destroyed on output   REAL(KREAL), INTENT(OUT) :: Step(NCoord)   INTEGER(KINT) :: NFree   REAL(KREAL), ALLOCATABLE :: eigval(:), eigvec(:,:), eigvli(:)   REAL(KREAL), ALLOCATABLE :: Tanf(:)   REAL(KREAL), ALLOCATABLE :: Grad_f(:),Step_f(:) ! NFree   REAL(KREAL), PARAMETER :: tiny=1e-20, zero=0., dele3=1e-6, eps=1e-12   REAL(KREAL), PARAMETER :: crit=1e-8   INTEGER(KINT) :: i, j, idx, k, isch   REAL(KREAL) :: grd, eval, norm   REAL(KREAL) :: dx   CHARACTER(120) :: fmt2   LOGICAL :: Debug   REAL(KREAL), ALLOCATABLE :: HFree(:,:) ! NFree,NFree   REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree  ! LOGICAL, SAVE :: LFirst=.TRUE.   INTERFACE   function valid(string) result (isValid)   CHARACTER(*), intent(in) :: string   logical :: isValid   END function VALID   END INTERFACE   debug=valid("step_rfo_all")   if (debug) WRITE(*,*) "=========================== Entering Step_RFO_All ===================="   if (debug) WRITE(*,*) "=========================== IGeom=",IGeom," ===================="  ! PFL 22.Nov.2007 ->  ! Nfree is equal to Ncoord-1 only when optimizing in the plane  ! orthogonal to the tangent.  ! In the case of a regular opt, NFree=Ncoord.  ! Next line has been moved later  ! NFree=Ncoord-1  ! <- PFL 22.Nov.2007   IF (debug) THEN   WRITE(*,*) 'DBG Step_RFO_All Hess'   WRITE(*,*) 'DBG Step_RFO_All NCoord=',NCoord   Idx=min(12,NCoord)   DO I=1,NCoord   WRITE(*,'(12(1X,F10.4))') Hess(I,1:Idx)   END DO   WRITE(*,*) 'DBG Step_RFO_All Grad'   WRITE(*,'(12(1X,F10.4))') Grad(1:Idx)   WRITE(*,*) 'DBG Step_RFO_All Coord'   WRITE(*,'(12(1X,F10.4))') Geom(1:Idx)   WRITE(*,*) 'DBG Step_RFO_All Tangent'   WRITE(*,'(12(1X,F10.4))') tangent(1:Idx)   END IF   Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord)  ! we orthogonalize Vfree to the tangent vector of this geom  ! 2007/05/29 Only if Tangent/=0.d0   Norm=sqrt(dot_product(Tangent,Tangent))   IF (Norm.GT.eps) THEN   ALLOCATE(Tanf(NCoord))   ! We normalize Tangent   Tangent=Tangent/Norm   ! We convert Tangent into vfree only displacements   ! This is useless for now (2007.Apr.23) as vfree=Id matrix   ! but it will be usefull as soon as we introduce Constraints   DO I=1,NCoord   Tanf(I)=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent)   END DO   Tangent=0.d0   DO I=1,NCoord   Tangent=Tangent+Tanf(i)*vfree(:,I)   END DO   ! first we substract Tangent from vfree   DO I=1,NCoord   Norm=dot_product(reshape(vfree(:,i),(/NCoord/)),Tangent)   Vfree(:,I)=vfree(:,i)-Norm*Tangent   END DO   Idx=0   ! Schmidt orthogonalization of the vfree vectors   DO I=1,NCoord   ! We substract the first vectors   ! we do it twice as the Schmidt procedure is not numerically stable   DO isch=1,2   DO J=1,Idx   Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),reshape(vfree(:,J),(/NCoord/)))   vfree(:,I)=vfree(:,I)-Norm*vfree(:,J)   END DO   END DO   Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),reshape(vfree(:,I),(/NCoord/)))   IF (Norm.GE.crit) THEN   idx=idx+1   vfree(:,idx)=vfree(:,i)/sqrt(Norm)   END IF   END DO     If (Idx/= NCoord-1) THEN   WRITE(*,*) "Pb in orthogonalizing vfree to tangent for geom",IGeom   WRITE(Ioout,*) "Pb in orthogonalizing vfree to tangent for geom",IGeom   STOP   END IF     if (debug) WRITE(*,*) "Deallocating Tanf"   DEALLOCATE(Tanf)   NFree=Idx   ELSE ! Tangent =0, matches IF (Norm.GT.eps) THEN   if (debug) WRITE(*,*) "Tangent=0, using full displacement"   NFree=NCoord   END IF !IF (Norm.GT.eps) THEN     if (debug) WRITE(*,*) 'DBG Step_RFO_All NFree=',NFree    ! We now calculate the new step  ! we project the hessian onto the free vectors   ALLOCATE(hfree(NFree,NFree),htmp(NCoord,NFree),Grad_f(NFree),Step_f(NFree))   DO J=1,NFree   DO I=1,NCoord   Htmp(I,J)=0.d0   DO K=1,NCoord   Htmp(I,J)=htmp(I,J)+Hess(I,K)*vfree(K,J)   END DO   END DO   END DO   DO J=1,NFree   DO I=1,NFree   HFree(i,j)=0.d0   DO K=1,NCoord   HFree(i,j)=hfree(i,j)+vfree(k,i)*htmp(k,j)   END DO   END DO   END DO  ! We diagonalize the hessian or its inverse   ALLOCATE(eigval(NFree), eigvec(NFree,NFree), eigvli(NFree))   call Jacobi(Hfree,NFree,eigvli,eigvec,NFree)   if (debug) THEN   WRITE(*,*) 'Eigensystem, Step_RFO_all.f90, L225'   fmt2='(I5,":",F8.3," -", (1X,F10.4))'   WRITE(fmt2(19:20),'(I2)') min(NFree,99)   DO I=1,NFree   WRITE(*,fmt2) I,Eigvli(I),eigvec(:,I)   END DO   END IF   ! We inverse the eigenvalues if Hess corresponds to Hessian inverse   IF (Hinv) THEN   do J=1,NFree   if (abs(eigvli(j)) <= tiny) eigvli(j) = zero   eigval(j) = sign (eigvli(j)**2 / (dele3 + abs(eigvli(j))**3), eigvli(j))   if (abs(eigval(j)) <= tiny) eigval(j) = zero   end do   ELSE   eigval=eigvli   END IF   call trie(NFree,eigval,eigvec,NFree)   DO I=1,NFree   Grad_f(I)=dot_product(reshape(vfree(:,I),(/NCoord/)),grad)   END DO   ! We now construct the step   Step_f=0.d0   DO I=1,NFree   grd= dot_product(grad_f,reshape(eigvec(:,i),(/NFree/)))   eval = 0.5 * (abs(eigval(i)) + sqrt (eigval(i)**2 + 4.*grd**2))   dx = (-grd) / eval   step_f = step_f + dx * eigvec(:,i)   if (debug) write (*,*) i,eigval(i),eval,grd,dx   END DO   Step=0.d0  ! PFL 04 Apr 2008 ->  ! Grad is once again modified  ! Grad=0.d0   DO I=1,NFree   Step=Step+ step_f(i)*vfree(:,i)  ! Grad=Grad+ Grad_f(i)*vfree(:,i)   END DO  ! <- PFL 04 Apr 2008     if (debug) WRITE(*,*) "Deallocataing Eigval"   DEALLOCATE(eigval)   if (debug) WRITE(*,*) "Deallocataing Eigvec"   DEALLOCATE(eigvec)   if (debug) WRITE(*,*) "Deallocataing Eigvli"   DEALLOCATE(eigvli)   if (debug) WRITE(*,*) "Deallocataing hfree"   DEALLOCATE(hfree)   if (debug) WRITE(*,*) "Deallocataing htmp"   DEALLOCATE(htmp)   if (debug) WRITE(*,*) "Deallocataing grad_f"   DEALLOCATE(grad_f)   if (debug) WRITE(*,*) "Deallocataing step_f"   DEALLOCATE(step_f)   if (debug) THEN   WRITE(*,*) 'DBG Step_RFO_All Step'   WRITE(*,'(12(1X,F8.4))') Step(1:Idx)   END IF   if (debug) WRITE(*,*) "=========================== Step_RFO_All Over ===================="  END SUBROUTINE STEP_RFO_ALL