Statistics
| Revision:

## root / src / Step_RFO_all.f90 @ 7

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