Statistics
| Revision:

## root / src / Extrapol_int.f90 @ 7

 1 SUBROUTINE Extrapol_int(iopt,s,dist,x0,y0,z0,xgeom,Coef,XgeomF)   ! 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   !   ! Default value of FRot=.TRUE.   ! Default value of NMaxPtPath = 1000   ! IntCoordI(:,:) ! (NGeomI,3*Nat-6)   ! XyzGeomF(:,:,:) ! (NGeomF,3,Nat)   ! Default value of FAlign=.TRUE.   use Path_module, only : NMaxPtPath, IntCoordI, Pi, IndZmat, XyzGeomF, &   IntCoordF, IntTangent, Renum, Nom, Order, MassAt, SGeom, XyzGeomI, &   Atome, Nat, NGeomI, NCoord, NGeomF, OrderInv,NFroz,FrozAtoms,Align   ! IndZmat(Nat,5)   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   REAL(KREAL), INTENT(OUT) :: XGeomF(NGeomF)   INTEGER(KINT) :: IdxGeom, I, J, K, Idx   REAL(KREAL) :: Rmsd,MRot(3,3), ds, u, v   REAL(KREAL) :: a_val, d   REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)   REAL(KREAL), ALLOCATABLE :: ValZmat(:,:),DerInt(:) ! (Nat,3)   LOGICAL :: debug, print,printspline   LOGICAL, EXTERNAL :: valid   INTEGER(KINT) :: NSpline   CHARACTER(LCHARS) :: FileSpline,TmpChar   ! LOGICAL :: FAlign=.FALSE., FRot=.FALSE.   ! 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_int")  ! 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),valzmat(Nat,3),DerInt(NCoord))   IdxGeom=1   XYZTmp2(1,1)=0.   XYZTmp2(1,2)=0.   XYZTmp2(1,3)=0.   XYZTmp2(2,2)=0.   XYZTmp2(2,3)=0.   XYZTmp2(3,3)=0.   valzmat=0.   valzmat(2,1)=IntCoordI(1,1)   valzmat(3,1)=IntCoordI(1,2)   valzmat(3,2)=IntCoordI(1,3)*180./Pi   Idx=4   DO I=4,Nat   ValZmat(I,1)=IntCoordI(1,Idx)   Idx=Idx+1   DO J=2,3   ValZmat(I,J)=IntCoordI(1,Idx)*180./Pi   Idx=Idx+1   END DO   END DO   XyzTmp2(2,1)=valzmat(2,1)   d=valzmat(3,1)   a_val=valzmat(3,2)/180.*Pi   ! write(*,*) "aval,pi",a_val,valzmat(3,2),pi   if (Nat.GE.3) THEN   if (IndZmat(3,2).EQ.1) THEN   XyzTmp2(3,1)=XYzTmp2(1,1)+d*cos(a_val)   ELSE   XyzTmp2(3,1)=XyzTmp2(2,1)-d*cos(a_val)   ENDIF   XyzTmp2(3,2)=d*sin(a_val)   ENDIF   ! i=1   ! WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)   ! WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)   ! i=2   ! WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)   ! WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)   ! i=3   ! WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)   ! WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)   DO i=4,Nat   call ConvertZmat_cart(i,IndZmat,valzmat, &   XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))   ! WRITE(*,*) 'TOTOZma:',i,IndZmat(I,1), &   ! (IndZmat(I,J+1),valzmat(I,J),J=1,3)   ! WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)   END DO   ! We align this geometry with the original one   ! PFL 17/July/2006: only if we have more than 4 atoms.  ! PFL 2013 Feb 27 ... or if the user asks for it   IF ((Nat.GE.4).OR.Align) THEN   ! The rotation matrix MRot has INTENT(OUT) in subroutine rotation_matrix(...),   ! which is called in the CalcRmsd(...).   ! PFL 24 Nov 2008 ->   ! If we have frozen atoms we align only those ones.   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   XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))   IntCoordF(IdxGeom,:)=IntCoordI(1,:)   ! We calculate the first derivatives   u=0.d0   DO Idx=1,NCoord   call splintder(u,v,DerInt(iDX),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))   END DO   IntTangent(IdxGeom,:)=DerInt   if (print.AND.(Dist.LE.1e20)) THEN   WRITE(IOOUT,'(1X,I5)') Nat   WRITE(IOOUT,*) "# Cartesian Coordinates for geom",IdxGeom   DO i=1,Nat   If (Renum) THEN   WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), &   (XyzTmp2(Order(I),J),J=1,3)   ELSE   WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))), &   (XyzTmp2(I,J),J=1,3)   END IF   END DO   END IF   ! We initialize the first geometry   ! PFL 26 Nov 2008 ->   ! Now, I copy XyzTmp2 into XyzTmp because XyzTmp2 has been aligned   ! with the reference geometry   XyzTmp=XyzTmp2   ! XyzTmp(1,1)=0.   ! XyzTmp(1,2)=0.   ! XyzTmp(1,3)=0.   ! XyzTmp(2,2)=0.   ! XyzTmp(2,3)=0.   ! XyzTmp(3,3)=0.   ! XyzTmp(2,1)=valzmat(2,1)   ! d=valzmat(3,1)   ! a_val=valzmat(3,2)/180.*Pi   ! if (Nat.GE.3) THEN   ! if (IndZmat(3,2).EQ.1) THEN   ! XyzTmp(3,1)=XYzTmp(1,1)+d*cos(a_val)   ! ELSE   ! XyzTmp(3,1)=XyzTmp(2,1)-d*cos(a_val)   ! ENDIF   ! XyzTmp(3,2)=d*sin(a_val)   ! ENDIF   ! DO i=4,Nat   ! call ConvertZmat_cart(i,IndZmat,valzmat, &   ! XYzTMP(1,1), XYzTMP(1,2),XYzTMP(1,3))   ! END DO   ! <- PFL 26 Nov 2008   s=0.   if (printspline) THEN   SELECT CASE (Nat)   CASE(2)   WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)   CASE (3)   WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)   CASE(4:)   WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)   END SELECT   END IF   valzmat=0.d0   DO K=1,NMaxPtPath   u=real(K)/NMaxPtPath*(NGeomI-1.)   XYZTmp2(1,1)=0.   XYZTmp2(1,2)=0.   XYZTmp2(1,3)=0.   XYZTmp2(2,2)=0.   XYZTmp2(2,3)=0.   XYZTmp2(3,3)=0.   ! We generate the interpolated Zmatrix   call splintder(u,v,DerInt(1),NGeomI,xgeom(1),IntCoordI(1,1),Coef(1,1))   valzmat(2,1)=v   call splintder(u,v,DerInt(2),NGeomI,xgeom(1),IntCoordI(1,2),Coef(1,2))   valzmat(3,1)=v   call splintder(u,v,DerInt(3),NGeomI,xgeom(1),IntCoordI(1,3),Coef(1,3))   valzmat(3,2)=v*180./Pi   Idx=4   DO I=4,Nat   call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))   valzmat(I,1)=v   Idx=Idx+1   DO J=2,3   call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))   valzmat(I,J)=v*180./pi   Idx=Idx+1   END DO   END DO ! matches DO I=4,Nat   ! We convert it into Cartesian coordinates   XyzTmp2(2,1)=valzmat(2,1)   d=valzmat(3,1)   a_val=valzmat(3,2)/180.*Pi   if (Nat.GE.3) THEN   if (IndZmat(3,2).EQ.1) THEN   XyzTmp2(3,1)=XYzTmp2(1,1)+d*cos(a_val)   ELSE   XyzTmp2(3,1)=XyzTmp2(2,1)-d*cos(a_val)   ENDIF   XyzTmp2(3,2)=d*sin(a_val)   ENDIF   DO I=4,Nat   call ConvertZmat_cart(I,IndZmat,valzmat, &   XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))   END DO   IF ((Nat.GE.4).OR.Align) THEN   ! PFL 24 Nov 2008 ->   ! If we have frozen atoms we align only those ones.   if (NFroz.GT.0) THEN   Call AlignPartial(Nat,x0,y0,z0, &   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &   FrozAtoms,MRot)   ELSE   ! PFL 26 Nov 2008!   ! Note to myself: in Extrapol_cart we always align on x0,y0,z0 but here   ! we align on the previous geom... what is the best ? Is there a difference ?   Call CalcRmsd(Nat, XyzTmp(1:Nat,1),XyzTmp(1:Nat,2),XyzTmp(1:Nat,3), &   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &   MRot,rmsd,.TRUE.,.TRUE.)   END IF   ! <- PFL 24 Nov 2008   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)   END DO   END DO     s=s+sqrt(ds)   if ((printspline).AND.(MOD(K,NSpline).EQ.0)) THEN   SELECT CASE (Nat)   CASE(2)   WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)   CASE (3)   WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)   CASE(4:)   WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)   END SELECT   END IF   ! if (debug) WRITE(*,*) "Debug u,s,dist",u,s,dist   if (s>=dist) THEN   if (debug) THEN   WRITE(*,*) "DBG Zmat",s   WRITE(*,'(1X,I5)') IndZmat(1,1)   WRITE(*,'(1X,I5,I5,1X,F10.4)') IndZmat(2,1),IndZmat(2,2),valzmat(2,1)   WRITE(*,'(1X,I5,2(I5,1X,F10.4))') IndZmat(3,1),IndZmat(3,2),valzmat(3,1) &   ,IndZmat(3,3),valzmat(3,2)   DO I=4,Nat   WRITE(*,'(1X,I5,3(I5,1X,F10.4))') IndZmat(I,1),IndZmat(I,2),valzmat(I,1) &   ,IndZmat(I,3),valzmat(I,2),IndZmat(I,4),valzmat(I,3)   END DO   END IF ! matches if (debug) THEN   s=s-dist   IdxGeom=IdxGeom+1   SGeom(IdxGeom)=s+IdxGeom*dist   XgeomF(IdxGeom)=u   XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))   IntCoordF(IdxGeom,1)=valzmat(2,1)   IntCoordF(IdxGeom,2)=valzmat(3,1)   IntCoordF(IdxGeom,3)=valzmat(3,2)/180.*Pi   Idx=4   DO I=4,Nat   IntCoordF(IdxGeom,Idx)=valzmat(I,1)   IntCoordF(IdxGeom,Idx+1)=valzmat(I,2)/180.*Pi   IntCoordF(IdxGeom,Idx+2)=valzmat(I,3)/180.*Pi   Idx=Idx+3   END DO   IntTangent(IdxGeom,:)=DerInt   if (print) THEN   WRITE(IOOUT,'(1X,I5)') Nat   WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K   ! PFL 17/July/2006: only if we have more than 4 atoms.   IF ((Nat.GE.4).OR.Align) THEN   ! PFL 24 Nov 2008 ->   ! If we have frozen atoms we align only those ones.   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   DO I=1,Nat   If (Renum) THEN   WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), &   (XyzTmp2(Order(I),J),J=1,3)   ELSE   WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))), &   (XyzTmp2(I,J),J=1,3)   END IF   END DO   END IF ! matches if (print) THEN   END IF ! matches if (s>=dist) THEN   END DO ! matches DO K=1,NMaxPtPath, Is it correct??   if (printspline) THEN   SELECT CASE (Nat)   CASE(2)   WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)   CASE (3)   WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)   CASE(4:)   WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)   END SELECT   END IF   if (s>=0.9*dist) THEN   if (debug) WRITE(*,*) "DBG Extrapol_int L383: Adding last geom"   if (debug) write(*,*) "Extrapol_int u,xgeom(NGeomI),s,dist,s-dist",u,xgeom(NGeomI),s,dist,s-dist  ! u=xgeom(NGeomI)   s=s-dist  ! ! We generate the interpolated Zmatrix  ! call splintder(u,v,DerInt(1),NGeomI,xgeom(1),IntCoordI(1,1),Coef(1,1))  ! valzmat(2,1)=v  ! call splintder(u,v,DerInt(2),NGeomI,xgeom(1),IntCoordI(1,2),Coef(1,2))  ! valzmat(3,1)=v  ! call splintder(u,v,DerInt(3),NGeomI,xgeom(1),IntCoordI(1,3),Coef(1,3))  ! valzmat(3,2)=v*180./Pi  ! Idx=4  ! DO I=4,Nat  ! call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))  ! valzmat(I,1)=v  ! Idx=Idx+1  ! DO J=2,3  ! call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))  ! valzmat(I,J)=v*180./pi  ! Idx=Idx+1  ! END DO  ! END DO ! matches DO I=4,Nat  ! ! We convert it into Cartesian coordinates  ! XyzTmp2(2,1)=valzmat(2,1)  ! d=valzmat(3,1)  ! a_val=valzmat(3,2)/180.*Pi  ! if (Nat.GE.3) THEN  ! if (IndZmat(3,2).EQ.1) THEN  ! XyzTmp2(3,1)=XYzTmp2(1,1)+d*cos(a_val)  ! ELSE  ! XyzTmp2(3,1)=XyzTmp2(2,1)-d*cos(a_val)  ! ENDIF  ! XyzTmp2(3,2)=d*sin(a_val)  ! ENDIF  ! DO I=4,Nat  ! call ConvertZmat_cart(I,IndZmat,valzmat, &  ! XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))  ! END DO  ! IF (Nat.GE.4) THEN  ! ! PFL 24 Nov 2008 ->  ! ! If we have frozen atoms we align only those ones.  ! if (NFroz.GT.0) THEN  ! Call AlignPartial(Nat,x0,y0,z0, &  ! xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &  ! FrozAtoms,MRot)  ! ELSE  ! ! PFL 26 Nov 2008!  ! ! Note to myself: in Extrapol_cart we always align on x0,y0,z0 but here  ! ! we align on the previous geom... what is the best ? Is there a difference ?  ! Call CalcRmsd(Nat, XyzTmp(1:Nat,1),XyzTmp(1:Nat,2),XyzTmp(1:Nat,3), &  ! xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), &  ! MRot,rmsd,.TRUE.,.TRUE.)  ! END IF  ! ! <- PFL 24 Nov 2008  ! END IF   IdxGeom=IdxGeom+1   IF (IdxGeom.GT.NGeomF) THEN   WRITE(IOOUT,*) "!!! ERROR in Extrapol_int !!!!"   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   SGeom(IdxGeom)=s+IdxGeom*dist   XgeomF(IdxGeom)=min(u,NGeomI-1.d0)   XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))   ! XyzGeomF(IdxGeom,:,:)=XyzTmp2(:,:)   IntCoordF(IdxGeom,1)=valzmat(2,1)   IntCoordF(IdxGeom,2)=valzmat(3,1)   IntCoordF(IdxGeom,3)=valzmat(3,2)/180.*Pi   Idx=4   DO I=4,Nat   IntCoordF(IdxGeom,Idx)=valzmat(I,1)   IntCoordF(IdxGeom,Idx+1)=valzmat(I,2)/180.*Pi   IntCoordF(IdxGeom,Idx+2)=valzmat(I,3)/180.*Pi   Idx=Idx+3   END DO   IntTangent(IdxGeom,:)=DerInt   if (print) THEN   WRITE(IOOUT,'(1X,I5)') Nat   WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K   ! PFL 17/July/2006: only if we have more than 4 atoms.   IF ((Nat.GE.4).OR.Align) THEN   ! PFL 24 Nov 2008 ->   ! If we have frozen atoms we align only those ones.   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   DO I=1,Nat   If (Renum) THEN   WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), &   (XyzTmp2(Order(I),J),J=1,3)   ELSE   WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))), &   (XyzTmp2(I,J),J=1,3)   END IF   END DO ! matches DO I=1,Nat   END IF ! matches if (print) THEN   END IF ! matches if (s>=0.9*dist) THEN   if (debug) WRITE(*,*) 's final =',s   DEALLOCATE(XyzTmp,XyzTmp2,valzmat)   if (printspline) CLOSE(IOTMP)   if (debug) Call Header("Extrapol_int Over")  END SUBROUTINE EXTRAPOL_INT