Statistics
| Revision:

root / src / Extrapol_cart.f90 @ 5

 1 SUBROUTINE Extrapol_cart(iopt,s,dist,x0,y0,z0,xgeom,Coef)   ! This subroutine constructs the path, and if dist<>Infinity, it samples   ! the path to obtain geometries.   ! Basically, you call it twice: i) dist=infinity, it will calculate the length of the path   ! ii) dist finite, it will give you the images you want along the path.   !   ! For now, it gives equidistant geometries   !   use Path_module   use Io_module   IMPLICIT NONE   REAL(KREAL), INTENT(OUT) :: s   ! X0(Nat),Y0(Nat),Z0(Nat): reference geometry.   REAL(KREAL), INTENT(IN) :: dist,X0(Nat),Y0(Nat),Z0(Nat)   REAL(KREAL), INTENT(IN) :: Xgeom(NGeomI),Coef(NGeomI,NCoord)  ! Iopt: Number of the cycles for the optimization   INTEGER(KINT), INTENT(IN) :: Iopt   INTEGER(KINT) :: IdxGeom, I, J, K, Iat  ! NSpline is the number of points along the interpolating path   INTEGER(KINT) :: NSpline  ! FileSpline: Filename to save the interpolating path coordinates   CHARACTER(LCHARS) :: FileSpline,TmpChar   REAL(KREAL) :: Rmsd, MRot(3, 3), ds, u   REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)   REAL(KREAL), ALLOCATABLE :: DerXyz(:,:) ! Nat,3   LOGICAL :: debug, print, printspline   LOGICAL, EXTERNAL :: valid   !We will calculate the length of the path, in MW coordinates...   ! this is done is a stupid way: we interpolate the zmatrix values,   ! convert them into cartesian, weight the cartesian   ! and calculate the evolution of the distance !   ! We have to follow the same procedure for every geometry   ! so even for the first one, we have to convert it from zmat   ! to cartesian !   debug=valid("pathcreate")   print=valid("printgeom")   printspline=(valid("printspline").AND.(dist<=1e30))   if (debug) Call Header("Entering Extrapol_cart")  ! We want 100 points along the interpolating path   NSpline=int(NMaxPtPath/100)   if (printspline) THEN   WRITE(TmpChar,'(I5)') Iopt   FileSpline=Trim(adjustL(PathName)) // '_spline.' // AdjustL(TRIM(TmpChar))   OPEN(IOTMP,FILE=FileSpline)   END IF   ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),DerXyz(Nat,3))  ! PFL 2011 June 22  ! There was an alignment procedure there.  ! It has been moved to PathCreate before the  ! computation of the Spline coefficient.   if (debug) THEN   WRITE(*,*) "Extrapol_cart- L72 - geometries "   DO I=1,NGeomI   WRITE(*,*) "Geom ",I   DO J=1,Nat   If (renum) THEN   Iat=Order(J)   WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(J)),XyzGeomI(I,1:3,Iat)   ELSE   Iat=OrderInv(J)   WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),XyzGeomI(I,1:3,J)   END IF   END DO   END DO   END IF   XyzGeomF(1,:,:)=XyzGeomI(1,:,:)   ! We initialize the first geometry   XyzTmp=Reshape(XyzGeomI(1,:,:),(/Nat,3/),ORDER=(/2,1/))   s=0.d0   SGeom(1)=0.d0   if (printspline) THEN   u=0.d0   DO Iat=1,Nat   ! We generate the interpolated coordinates   if (Linear) THEN   call LinearInt(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat))   call LinearInt(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat))   call LinearInt(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat))   ELSE   call splintder(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat),Coef(1,3*Iat-2))   call splintder(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat),Coef(1,3*Iat-1))   call splintder(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat),Coef(1,3*Iat))   END IF   END DO   WRITE(IOTMP,*) Nat   WRITE(IOTMP,'(1X,A,1X,F15.6)') "Debug Spline - Coord=Cart,s=",s   DO Iat=1,Nat   WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)), &   (XyzTmp2(Iat,1:3))   END DO   END IF   IdxGeom=1   DO K=1,NMaxPtPath   u=real(K)/NMaxPtPath*(NGeomI-1.)   DO Iat=1,Nat   ! We generate the interpolated coordinates   if (Linear) THEN   call LinearInt(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat))   call LinearInt(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat))   call LinearInt(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat))   ELSE   call splintder(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat),Coef(1,3*Iat-2))   call splintder(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat),Coef(1,3*Iat-1))   call splintder(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat),Coef(1,3*Iat))   END IF   END DO  ! PFL 2011 June 22: We now align on the previous one  ! as we do for internal coord interpolation  ! We align this geometry with the previous one   ! PFL 17/July/2006: only if we have more than 4 atoms.  ! IF (Nat.GT.4) THEN  ! PFL 24 Nov 2008 ->  ! If we have frozen atoms we align only those ones.  ! PFL 8 Feb 2011 ->  ! I add a flag to see if we should align or not.  ! For small systems, it might be better to let the user align himself   IF (Align) THEN   if (NFroz.GT.0) THEN   Call AlignPartial(Nat, &   xyzTmp(1,1),xyzTmp(1,2),xyzTMP(1,3), &   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &   FrozAtoms,MRot)   ELSE   Call CalcRmsd(Nat, &   xyzTmp(1,1),xyzTmp(1,2),xyzTMP(1,3), &   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &   MRot,rmsd,.TRUE.,.TRUE.)   END IF  ! <- PFL 24 Nov 2008   END IF  ! -> PFL 8 Feb 2011  ! END IF   ds=0.   DO I=1,Nat   DO J=1,3   ds=ds+MassAt(I)*(XYZTMp2(I,J)-XYZTmp(I,J))**2   XYZTmp(I,J)=XyzTMP2(I,J)   ENDDO   ENDDO   s=s+sqrt(ds)   ! if (debug) WRITE(*,*) "Debug u,s,dist",u,s,dist   if ((printspline).AND.(MOD(K,NSpline).EQ.0)) THEN   WRITE(IOTMP,*) Nat   WRITE(IOTMP,'(1X,A,1X,F15.6)') "Debug Spline - Coord=Cart,s=",s   DO Iat=1,Nat   WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)), &   (XyzTmp2(Iat,1:3))   END DO   END IF   if (s>=dist) THEN   if (debug) THEN   WRITE(*,*) "DBG Interpol_cart",s   DO i=1,Nat   WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), &   (XyzTmp2(I,J),J=1,3)   END DO   END IF   s=s-dist   IdxGeom=IdxGeom+1   SGeom(IdxGeom)=s+dist*(IdxGeom-1)   ! We align this geometry with the original one   ! PFL 17/July/2006: only if we have more than 4 atoms.  ! IF (Nat.GT.4) THEN  ! PFL 24 Nov 2008 ->  ! If we have frozen atoms we align only those ones.  ! PFL 8 Feb 2011 ->  ! I add a flag to see if we should align or not.  ! For small systems, it might be better to let the user align himself   IF (Align) THEN   if (NFroz.GT.0) THEN   Call AlignPartial(Nat,x0,y0,z0, &   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &   FrozAtoms,MRot)   ELSE   Call CalcRmsd(Nat, x0,y0,z0, &   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &   MRot,rmsd,.TRUE.,.TRUE.)   END IF  ! <- PFL 24 Nov 2008   END IF  ! -> PFL 8 Feb 2011  ! END IF   XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))   XyzTangent(IdxGeom,:)=Reshape(DerXyz(:,:),(/3*Nat/))   if (print) THEN   WRITE(IOOUT,'(1X,I5)') Nat   WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K   DO i=1,Nat   WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), &   (XyzTmp2(I,J),J=1,3)   END DO   END IF   END IF ! s>= dist   ENDDO ! K   if (s>=0.9*dist) THEN   s=s-dist   XyzTmp2=Reshape(XyzGeomI(NGeomI,:,:),(/Nat,3/),ORDER=(/2,1/))   if (printspline) THEN   WRITE(IOTMP,*) Nat   WRITE(IOTMP,*) "Debug Spline - Coord=Cart,s=",s   DO Iat=1,Nat   WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)), &   (XyzTmp2(Iat,1:3))   END DO   END IF   ! We align this geometry with the original one   ! PFL 17/July/2006: only if we have more than 4 atoms.  ! IF (Nat.GT.4) THEN  ! PFL 24 Nov 2008 ->  ! If we have frozen atoms we align only those ones.  ! PFL 8 Feb 2011 ->  ! I add a flag to see if we should align or not.  ! For small systems, it might be better to let the user align himself   IF (Align) THEN   if (NFroz.GT.0) THEN   Call AlignPartial(Nat,x0,y0,z0, &   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &   FrozAtoms,MRot)   ELSE   Call CalcRmsd(Nat, x0,y0,z0, &   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &   MRot,rmsd,.TRUE.,.TRUE.)   END IF  ! <- PFL 24 Nov 2008   END IF  ! -> PFL 8 Feb 2011  ! END IF   IdxGeom=IdxGeom+1   SGeom(IdxGeom)=s+dist*(IdxGeom-1)   IF (IdxGeom.GT.NGeomF) THEN   WRITE(IOOUT,*) "!!! ERROR in Extrapol_cart !!!!"   WRITE(IOOUT,*) "Too many structures. Increase NMaxPath"   WRITE(*,*) "** PathCreate ***"   WRITE(*,*) "Distribution of points along the path is wrong."   WRITE(*,*) "Increase value of NMaxPtPath in the input file"   WRITE(*,*) "Present value is:", NMaxPtPath   STOP   END IF   ! We align this geometry with the original one   ! PFL 17/July/2006: only if we have more than 4 atoms.  ! IF (Nat.GT.4) THEN  ! PFL 24 Nov 2008 ->  ! If we have frozen atoms we align only those ones.  ! PFL 8 Feb 2011 ->  ! I add a flag to see if we should align or not.  ! For small systems, it might be better to let the user align himself   IF (Align) THEN   if (NFroz.GT.0) THEN   Call AlignPartial(Nat,x0,y0,z0, &   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &   FrozAtoms,MRot)   ELSE   Call CalcRmsd(Nat, x0,y0,z0, &   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &   MRot,rmsd,.TRUE.,.TRUE.)   END IF  ! <- PFL 24 Nov 2008   END IF  ! -> PFL 8 Feb 2011  ! END IF   XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))   XyzTangent(IdxGeom,:)=Reshape(DerXyz(:,:),(/3*Nat/))   if (print) THEN   WRITE(IOOUT,'(1X,I5)') Nat   WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K   DO i=1,Nat   WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), &   (XyzTmp2(I,J),J=1,3)   END DO   END IF   END IF   DEALLOCATE(XyzTmp,XyzTmp2)   if (debug) WRITE(*,*) 's final =',s   if (printspline) CLOSE(IOTMP)   if (debug) Call Header("Extrapol_cart over")  END SUBROUTINE EXTRAPOL_CART