Statistics
| Revision:

## root / src / Calc_zmat_frag.f90 @ 8

 1 SUBROUTINE Calc_Zmat_frag(na,atome,ind_zmat,val_zmat,x,y,z,r_cov,fact)  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !  ! Goal: to compute a viable Z-Matrix starting from the  ! cartesian coordinates of the atoms  !  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !  ! Input:  ! na : Number or atoms  ! atome : Mass number of the atoms (H=1, He=2,...)  ! x,y,z : cartesian coordinates of the atoms  ! r_cov : array containing the covalent radii of the atoms  ! fact : multiplicative factor used to determine if two atoms are linked.  ! see CalcCnct for more details.  !  ! Output:  ! ind_zmat : INTEGER(na,5) contains the indices of the Z-matrix  ! val_zmat : REAL(na,3) contains the values of the Z-matrix  !  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  ! History  !  ! v1.0 written for Cart a long time ago. Does not use fragment analysis, but was quite good !  !  ! v2.0 12/2007  ! Comes from Calc_zmat_constr_frag, based on a fragment analysis of the  ! system under study.  ! It should be more flexible and robust.  !  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   Use Path_module, only : max_Z, NMaxL, Nom, Pi   Use Io_module   IMPLICIT NONE   integer(KINT), INTENT(IN) :: na,atome(na)   real(KREAL), INTENT(IN) :: x(Na),y(Na),z(Na),fact   real(KREAL), INTENT(IN) :: r_cov(0:Max_Z)   INTEGER(KINT), INTENT(OUT) :: ind_zmat(Na,5)   real(KREAL), INTENT(OUT) :: val_zmat(Na,3)   INTEGER(KINT) :: idx_zmat(NA)  ! Frozen contains the indices of frozen atoms  ! INTEGER(KINT) Frozen(*),NFroz  ! Number of fragment, Index of the current fragment for loops   INTEGER(KINT) :: NbFrag  ! Fragment(I) contains the fragment index to which I belongs  ! NbAtFrag(j) contains the number of atoms of fragment j   INTEGER(KINT), ALLOCATABLE :: Fragment(:),NbAtFrag(:) !na  ! FragAt(i,:) lists the atoms of fragment i   INTEGER(KINT), ALLOCATABLE :: FragAt(:,:) !na,na  ! MaxLFrag(i,1) contains the maximum of links for an atom for fragment i  ! MaxLFrag(i,2) is the atom that has this number of linked atoms   INTEGER(KINT), ALLOCATABLE :: MaxLFrag(:,:) !na,2  ! INTEGER(KINT), ALLOCATABLE :: IdxFragAt(:) !na  ! INTEGER(KINT), ALLOCATABLE :: FrozBlock(:,:) !(na,0:na)  ! DistFrag contains the distance between a given atom and some other atoms   REAL(KREAL), ALLOCATABLE :: DistFrag(:) !na  ! FragDist(I) contains the index of the atoms corresponding to DistFrag(I)   INTEGER(KINT), ALLOCATABLE :: FragDist(:) ! na   INTEGER(KINT) :: IdxCaFaire, IAfaire   INTEGER(KINT), ALLOCATABLE :: LIAISONS(:,:) ! (Na,0:NMaxL)  ! INTEGER(KINT), ALLOCATABLE :: LiaisonsIni(:,:) ! (Na,0:NMaxL)   INTEGER(KINT), ALLOCATABLE :: CAFaire(:) ! (Na)   real(KREAL) :: vx1,vy1,vz1,norm1   real(KREAL) :: vx2,vy2,vz2,norm2   real(KREAL) :: vx3,vy3,vz3,norm3   real(KREAL) :: vx4,vy4,vz4,norm4   real(KREAL) :: vx5,vy5,vz5,norm5   real(KREAL) :: val,val_d, d12, d13,d23,d   Logical Debug, FirstAt, DebugGaussian   LOGICAL, ALLOCATABLE :: DejaFait(:), FCaf(:) !(na)  ! Logical, ALLOCATABLE :: FrozAt(:) !T if this atom is frozen   LOGICAL F1213, F1223,F1323   INTEGER(KINT) :: I,J,n0,n1,n2,n3,n4,IAt,IL,JL,IFrag,ITmp, K, KMax   INTEGER(KINT) :: I3, I1, Ip   INTEGER(KINT) :: I0, Izm, JAt   INTEGER(KINT) :: OrderZmat   REAL(KREAL) :: sDihe   INTERFACE   function valid(string) result (isValid)   CHARACTER(*), intent(in) :: string   logical :: isValid   END function VALID   FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)   INTEGER, PARAMETER :: KREAL=KIND(1.0D0)   real(KREAL) :: v1x,v1y,v1z,norm1   real(KREAL) :: v2x,v2y,v2z,norm2   real(KREAL) :: angle   END FUNCTION ANGLE   FUNCTION SinAngle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)   INTEGER, PARAMETER :: KREAL=KIND(1.0D0)   real(KREAL) :: v1x,v1y,v1z,norm1   real(KREAL) :: v2x,v2y,v2z,norm2   real(KREAL) :: SinAngle   END FUNCTION SINANGLE   FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3)   INTEGER, PARAMETER :: KREAL=KIND(1.0D0)   real(KREAL) :: v1x,v1y,v1z,norm1   real(KREAL) :: v2x,v2y,v2z,norm2   real(KREAL) :: v3x,v3y,v3z,norm3   real(KREAL) :: angle_d,ca,sa   END FUNCTION ANGLE_D   END INTERFACE   debug=valid("calc_zmat").OR.valid("calc_zmat_frag")   debugGaussian=valid("zmat_gaussian")     if (debug) WRITE(*,*) "=============================== Entering Calc_zmat_frag ========================"   if (na.le.2) THEN   WRITE(*,*) "I do not work for less than 2 atoms :-p"   RETURN   END IF   ALLOCATE(FragDist(Na),Fragment(na), NbAtFrag(na),FragAt(na,na))  ! ALLOCATE(FrozFrag(na,3), IdxFragAt(na), FrozBlock(na,0:na))  ! ALLOCATE(FrozFrag(na,3))   ALLOCATE(DistFrag(na),Liaisons(na,0:NMaxL))  ! ALLOCATE(Liaisons(na,0:NMaxL),LiaisonsIni(na,0:NMaxL))  ! ALLOCATE(CaFaire(na),DejaFait(Na),FCaf(na),FrozAt(na))   ALLOCATE(CaFaire(na+1),DejaFait(Na),FCaf(na))   if (debug) THEN   WRITE(*,*) "DBG Calc_zmat_frag - Cartesian coordinates"   DO I=1,na   WRITE(*,'(1X,A3,3(1X,F15.8))') Nom(atome(i)),x(i),y(i),z(i)   END DO   END if   DO I=1,na   DO J=1,5   Ind_Zmat(I,J)=0   END DO   END DO    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !  ! Easy case : 3 or less atoms  !  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   if (Na.eq.3) THEN   d12=sqrt((x(1)-x(2))**2+(y(1)-y(2))**2+(z(1)-z(2))**2)   d13=sqrt((x(1)-x(3))**2+(y(1)-y(3))**2+(z(1)-z(3))**2)   d23=sqrt((x(3)-x(2))**2+(y(3)-y(2))**2+(z(3)-z(2))**2)   F1213=(d12<=d13)   F1323=(d13<=d23)   F1223=(d12<=d23)   if (debug) THEN   WRITE(*,*) "DEBUG Calc_Zmat 3 atoms"   WRITE(*,*) "d12,d13,d23:",d12,d13,d23   WRITE(*,*) "F1213,F1323,F1223:",F1213,F1323,F1223   END IF   OrderZmat=0   if (F1213) orderZmat=OrderZmat+4   if (F1323) orderZmat=OrderZmat+2   if (F1223) orderZmat=OrderZmat+1   if (debug) WRITE(*,*) 'OrderZmat=',OrderZmat   SELECT CASE (OrderZmat)   CASE (0)  ! F F F ordre 2-3----1   ind_zmat(1,1)=3   ind_zmat(2,1)=2   ind_zmat(2,2)=3   ind_zmat(3,1)=1   ind_zmat(3,2)=3   ind_zmat(3,3)=2   CASE (2)  ! F T F ordre 1-3----2   ind_zmat(1,1)=3   ind_zmat(2,1)=1   ind_zmat(2,2)=3   ind_zmat(3,1)=2   ind_zmat(3,2)=3   ind_zmat(3,3)=1   CASE (3)  ! F T T ordre 2---1-3   ind_zmat(1,1)=1   ind_zmat(2,1)=3   ind_zmat(2,2)=1   ind_zmat(3,1)=2   ind_zmat(3,2)=1   ind_zmat(3,3)=3   CASE (5)  ! T F T ordre 1-2----3   ind_zmat(1,1)=2   ind_zmat(2,1)=1   ind_zmat(2,2)=2   ind_zmat(3,1)=3   ind_zmat(3,2)=2   ind_zmat(3,3)=1   CASE (7)  ! T T T ordre 3----1-2   ind_zmat(1,1)=1   ind_zmat(2,1)=2   ind_zmat(2,2)=1   ind_zmat(3,1)=3   ind_zmat(3,2)=1   ind_zmat(3,3)=2   END SELECT   IF (debug) THEN   WRITE(*,*) "DBG Calc_zmat_frag - Nat=3 -"   DO i=1,na   WRITE(*,'(1X,A3,5(1X,I5))') Nom(Atome(ind_zmat(i,1))),(ind_zmat(i,j),j=1,5)   END DO   END IF   ! We have ind_zmat, we fill val_zmat   val_zmat=0.d0   n2=ind_zmat(2,1)   n1=ind_zmat(2,2)   d=sqrt((x(n1)-x(n2))**2+(y(n1)-y(n2))**2+(z(n1)-z(n2))**2)   val_zmat(2,1)=d   n1=ind_zmat(3,1)   n2=ind_zmat(3,2)   n3=ind_zmat(3,3)   CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)   if (debug) write(*,*) n1,n2,norm1   val_zmat(3,1)=norm1   CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)   val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)   if (debug) write(*,*) n2,n3,norm2,val   val_zmat(3,2)=val   RETURN   END IF !matches if (Na.eq.3) THEN  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !  ! End of Easy case : 3 or less atoms  !  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !  ! General Case  !  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !  ! Initialization   DejaFait=.False.   Liaisons=0   ind_zmat=0   val_zmat=0.d0   if (debug) WRITE(*,*) "Coucou from Calc_zmat_frag.f90; L240"     if (debug) THEN   WRITE(*,*) "Bonds initialized"   WRITE(*,*) 'Covalent radii used'   DO iat=1,na   i=atome(iat)   WRITE(*,*) Nom(I),Iat,r_cov(i),r_cov(i)*fact   END DO   END IF  1003 FORMAT(1X,I4,' - ',25(I5))  ! First step : connectivity are calculated   Call CalcCnct(na,atome,x,y,z,LIAISONS,r_cov,fact)   if (debug) THEN   WRITE(*,*) "Connections calculated"   DO IL=1,na   WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)   END DO   END IF   FCaf=.TRUE.   Call Decomp_frag(na,liaisons,FCaf,nbfrag,Fragment,NbAtFrag,FragAt)     IF (debug) THEN   WRITE(*,*) 'I found ',NbFrag, 'fragments'   WRITE(*,*) (NbAtFrag(I),I=1,NbFrag)   DO I=1,NbFrag   WRITE(*,*) NbAtFrag(I)   WRITE(*,*) 'Fragment ', i   DO J=1,Na   IF (Fragment(J).EQ.I) WRITE(*,'(1X,I3,1X,A5,3(1X,F10.6))') J,Nom(Atome(J)) &   ,X(J),Y(J),Z(J)   END DO   WRITE(*,*) "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))   END DO   END IF   ALLOCATE(MaxLFrag(NbFrag,2))   MaxLFrag=0   DO I=1,NbFrag   MaxLFrag(I,1)=Liaisons(FragAt(I,1),0)   MaxLFrag(I,2)=FragAt(I,1)   DO J=1,NbAtFrag(I)   Iat=FragAt(I,J)   IF (Liaisons(IAt,0).GT.MaxLFrag(I,1)) THEN   MaxLFrag(I,1)=Liaisons(Iat,0)   MaxLFrag(I,2)=Iat   END IF   END DO   IF (debug) WRITE(*,*) 'Frag :',I,', atom ',MaxLFrag(I,2), ' has ',MaxLFrag(I,2),' links'   END DO   ! We will now build the molecule fragment by fragment   ! We choose the starting fragment with two criteria:   ! 1- Number of linked atoms:   ! * >=3 is good as it fully defines the coordinate space   ! * 2 is ok as we can either use a 3rd atom from the same fragment   ! or add a X atom somewhere but this complicates quite a lot the way   ! to treat the conversion from cartesian to zmat latter   ! * 1 is bad...   ! 2- Size of the fragment   ! this allows us to deal more easily with cases 1- when number of   ! directly linked atoms is less than 3   IFrag=0  ! I0 is the number of connections of the best fragment   I0=0  ! I1 is the number of atoms of the best fragment   I1=0   IAt=0   DO I=1,NbFrag   SELECT CASE(MaxLFrag(I,1)-I0)   CASE (1:)   IFrag=I   I0=MaxLFrag(I,1)   I1=NbAtFrag(I)   IAt=MaxLFrag(I,2)   CASE (0)   if (NbAtFrag(I).GT.I1) THEN   IFrag=I   I0=MaxLFrag(I,1)   I1=NbAtFrag(I)   IAt=MaxLFrag(I,2)   END IF   END SELECT   END DO   if (debug) WRITE(*,'(1X,A,I5,A,I5,A,I5,A,I5)') 'Starting with fragment:',IFrag,' with ',I0 &   ,' max links for atom',IAt,' fragment size',NbAtFrag(IFrag)   ! We will build the first fragment in a special way, as it will   ! set the coordinates system   if (debug) WRITE(*,*) 'Fragment 1, starting with atom:',IAt, &   'with ',I0,' connections'   DejaFait=.FALSE.   FCaf=.FALSE.     izm=0   SELECT CASE (I0)   CASE(3:)   if (debug) WRITE(*,*) 'DBG select case I0 3'   n0=Iat   ITmp=2   sDihe=0.   n2=IAt   n3=Liaisons(Iat,1)   ! We search for the third atom while making sure that it is not aligned with first two !   DO While ((ITmp.LE.Liaisons(Iat,0)).AND.(sDihe.LE.0.09d0))   n4=Liaisons(Iat,Itmp)   CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)   CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)   val_d=angle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)   sDiHe=abs(sin(val_d*pi/180.d0))   if (debug) Write(*,*) 'Trying n2,n3,n4,sdihe,val_d',n2,n3,n4,sdihe,val_d   Itmp=Itmp+1   END DO   If (debug) WRITE(*,*) 'Itmp,Liaisons',Itmp,Liaisons(Iat,1:NMaxL)   Liaisons(Iat,Itmp-1)=Liaisons(iat,2)   Liaisons(Iat,2)=n4   If (debug) WRITE(*,*) 'Itmp,Liaisons',Itmp,Liaisons(Iat,1:NMaxL)   if (sDihe.LE.0.09d0) THEN   WRITE(*,*) "PROBLEM !!! All atoms linked to ",n0," are aligned..."   WRITE(*,*) "Surprising as this atom has at least three bonds... For NOW STOP"   STOP   END IF   CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, vx5,vy5,vz5,norm5)  ! We search for the fourth atom while checking that it is not aligned with 1st and 2nd atoms.   Itmp=3   sDiHe=0.  ! PFL 28 Dec 2007 ->  ! I had a test on the dihedral angle, but I cannot see why it is important to have  ! a non planar fragment at the begining... ethylene works and is fully planar  ! I thus suppress this test  !  ! DO While ((ITmp.LE.Liaisons(Iat,0)).AND.(sDihe.LE.0.09d0))  ! ITmp=ITmp+1   n1=Liaisons(Iat,Itmp)   if (debug) WRITe(*,*) 'trying n1,n2,n3,n4',n1,n2,n3,n4   CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)   ! Is this atom aligned with n2-n3 ?   val_d=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)   sDiHe=abs(sin(val_d*pi/180.d0))   if (debug) Write(*,*) 'Angle n3-n2-n1',val_d   if (sDiHe.le.0.09d0) THEN   ! As atoms n2,n3 and n4 are not aligned, we interchange n3 and n4 so that n1,n2 and n3 are not aligned   if (debug) WRITE(*,*) "n3-n2-n1 aligned, we interchange n3 and n4"   n1=n3   n3=n4   n4=n1   n1=Liaisons(Iat,ITmp)   CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)   CALL vecteur(n2,n4,x,y,z,vx3,vy3,vz3,norm3)   val_d=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)   if (debug) Write(*,*) 'NEW Angle n3-n2-n1',val_d   END IF  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   ! avant de tester l'angle diedre, il faut verifier que ce 4e atome n'est pas   ! aligne avec les 2 premiers.   ! comme on a deja test? que les 3 premiers ne sont pas alignes,   ! s'il est align? avec les 2 premiers, on peut inverser le role de 2 et 3.   ! On pourrait tout simplier en faisant une bete recherche parmi tous les atomes geles   ! de ce bloc (au moins 4 ?) avec le critere 1) on les range par distance croissante   ! 2) on les scanne tant que l'angle valence n'est pas bon, puis tant que diehedre pas bon   ! on pourrait comme ca faire un tableau avec les atomes ranges d'abord pour le fragment s?lectionn?   ! puis les atomes des autres fragment par distance croissante.   ! les autres fragments ne seraient additonn?s que si l'on ne trouve pas notre bonheur dans le premier bloc  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)   val_d=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &   vx2,vy2,vz2,norm2)   sDihe=abs(sin(val_d*pi/180.d0))   if (debug) WRITE(*,*) 'n2,n3,n4,n1, angle_d',n2,n3,n4,n1,val_d  ! END DO   DejaFait(n2)=.TRUE.   DejaFait(n3)=.TRUE.   DejaFait(n4)=.TRUE.  ! if (sDihe.LE.0.09d0) THEN  ! ! None of the atoms linked to IAt are good to define the third  ! ! direction in space...  ! ! We will look at the other atoms  ! ! we might improve the search so as to take the atom closest to IAt  ! if (debug) WRITE(*,*) "All atoms linked to ",Iat," are in a plane. Looking for other atoms"  ! ITmp=0  ! DO I=1,NbAtFrag(IFrag)  ! JAt=FragAt(IFrag,I)  ! if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN  ! n1=JAt  ! CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)  ! CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)  ! val_d=angle_d(vx4,vy4,vz4,norm4, &  ! vx5,vy5,vz5,norm5, &  ! vx2,vy2,vz2,norm2)  ! if (abs(sin(val_d)).GE.0.09D0) THEN  ! ITmp=ITmp+1  ! DistFrag(ITmp)=Norm1  ! FragDist(ITmp)=JAt  ! END IF  ! END IF  ! END DO  ! IF (ITmp.EQ.0) THEN  ! ! All dihedral angles between frozen atoms are less than 5?  ! ! (or more than 175?). We have to look at other fragments :-(  ! DO I=1,NFroz  ! Jat=Frozen(I)  ! if (.NOT.DejaFait(JAt)) THEN  ! n1=JAt  ! CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)  ! CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)  ! val_d=angle_d(vx4,vy4,vz4,norm4, &  ! vx5,vy5,vz5,norm5, &  ! vx2,vy2,vz2,norm2)  ! if (abs(sin(val_d)).GE.0.09D0) THEN  ! ITmp=ITmp+1  ! DistFrag(ITmp)=Norm1  ! FragDist(ITmp)=JAt  ! END IF  ! END IF  ! END DO  ! IF (ITmp.EQ.0) THEN  ! ! All frozen atoms are in a plane... too bad  ! WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'  ! WRITE(*,*) 'For now, I do not treat this case'  ! STOP  ! END IF  ! END IF  ! I1=0  ! d=1e9  ! DO I=1,ITmp  ! IF (DistFrag(I).LE.d) THEN  ! d=DistFrag(I)  ! I1=FragDist(I)  ! END IF  ! END DO  ! ELSE !(sDihe.LE.0.09d0)  ! I1=FrozBlock(IAt,ITmp)  ! if (debug) WRITE(*,*) 'I1,n1:',I1,n1  ! END IF !(sDihe.LE.0.09d0)  ! ! we now have the atom that is closer to the first one and that  ! ! forms a non 0/Pi dihedral angle  !  ! <- PFL 28 Dec 2007  ! We construct the begining of the Z-Matrix   ind_zmat(1,1)=n2   ind_zmat(2,1)=n3   ind_zmat(2,2)=n2   ind_zmat(3,1)=n4   ind_zmat(3,2)=n2   ind_zmat(3,3)=n3   DejaFait(n2)=.TRUE.   DejaFait(n3)=.TRUE.   DejaFait(n4)=.TRUE.   CaFaire(1)=n3   CaFaire(2)=n4  ! PFL 28 Dec 2007  ! We have not selected a fourth atom, so that the following is not needed  ! ind_zmat(4,1)=I1  ! ind_zmat(4,2)=n2  ! ind_zmat(4,3)=n3  ! ind_zmat(4,4)=n4  ! DejaFait(I1)=.TRUE.  ! CaFaire(3)=I1  ! CaFaire(4)=0  ! IdxCaFaire=4  ! izm=4  ! FCaf(I1)=.TRUE.  !!!!!!!  ! and replaced by:   CaFaire(3)=0   IdxCaFaire=3   izm=3  !  ! <- PFL 28 Dec 2007     FCaf(n2)=.TRUE.   FCaf(n3)=.TRUE.   FirstAt=.TRUE.   DO I=3,Liaisons(Iat,0)   IF (.NOT.DejaFait(Liaisons(Iat,I))) THEN   izm=izm+1   ! ind_zmat(izm,1)=Liaisons(Iat,I)   ! ind_zmat(izm,2)=n2   ! ind_zmat(izm,3)=n3   ! ind_zmat(izm,4)=n4   Call add2indzmat(na,izm,Liaisons(Iat,I),n2,n3,n4,ind_zmat,x,y,z)   if (FirstAt) THEN   n4=Liaisons(Iat,I)   FirstAt=.FALSE.   END IF   IF (.NOT.FCaf(Liaisons(Iat,I))) THEN   CaFaire(IdxCaFaire)=Liaisons(Iat,I)   IdxCaFaire=IdxCaFaire+1   CaFaire(IdxCaFaire)=0   FCaf(Liaisons(Iat,I))=.TRUE.   END IF   DejaFait(Liaisons(Iat,I))=.TRUE.   END IF   END DO   if (debug) THEN   WRITE(*,*) "Ind_zmat 0 - SELECT CASE I0 3: -- izm=",izm   WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)   WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)   WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)   DO I=4,izm   WRITE(*,'(1X,4(1X,I5))') ind_zmat(I,1), ind_zmat(I,2), &   ind_zmat(I,3), ind_zmat(I,4)   END DO   END IF   ! First four atoms (at least) have been put we go on with next parts   ! of this fragment... later   CASE(2)   if (debug) WRITE(*,*) 'DBG select case I0 2'   WRITE(*,*) "PFL 28 Dec 2007: No test of alignment here :( TO DO TO DO TO DO"   ind_zmat(1,1)=IAt   ind_zmat(2,1)=Liaisons(IAt,1)   ind_zmat(2,2)=IAt   ind_zmat(3,1)=Liaisons(IAt,2)   ind_zmat(3,2)=IAt   ind_zmat(3,3)=Liaisons(IAt,1)   DejaFait(IAt)=.TRUE.   DejaFait(Liaisons(Iat,1))=.TRUE.   DejaFait(Liaisons(Iat,2))=.TRUE.   CaFaire(1)=Liaisons(Iat,1)   CaFaire(2)=Liaisons(Iat,2)   FCaf(Liaisons(Iat,1))=.TRUE.   FCaf(Liaisons(Iat,2))=.TRUE.  ! PFL 28 Dec 2007 ->  ! We do NOT need a fourth atom !!! The third direction in space is defined by the cross  ! product of the first two directions  !  ! the following is thus commented  !  ! ! We search for a fourth atom, first in the FrozBlock atoms  ! ITmp=2  ! sDihe=0.  ! n2=IAt  ! n3=Liaisons(Iat,1)  ! n4=Liaisons(Iat,2)  ! CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)  ! CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)  ! CALL produit_vect(vx3,vy3,vz3, vx2,vy2,vz2, vx5,vy5,vz5,norm5)  !  ! DO While ((ITmp.LE.FrozBlock(Iat,0)).AND.(sDihe.LE.0.09d0))  ! ITmp=ITmp+1  ! n1=FrozBlock(Iat,Itmp)  ! CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)  ! CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)  ! val_d=angle_d(vx4,vy4,vz4,norm4, &  ! vx5,vy5,vz5,norm5, &  ! vx2,vy2,vz2,norm2)  ! sDihe=abs(sin(val_d))  ! END DO  ! IF (debug) WRITE(*,*) 'DBG 4th atom, ITmp, sDihe',ITmp, sDihe  ! if (sDihe.LE.0.09d0) THEN  ! ! None of the frozen atoms linked to IAt are good to define the third  ! ! direction in space...  ! ! We will look at the other frozen atoms  ! ! we might improve the search so as to take the atom closest to IAt  ! ITmp=0  ! DO I=1,NbAtFrag(IFrag)  ! JAt=FragAt(IFrag,I)  ! if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN  ! n1=JAt  ! CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)  ! CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)  ! val_d=angle_d(vx4,vy4,vz4,norm4, &  ! vx5,vy5,vz5,norm5, &  ! vx2,vy2,vz2,norm2)  ! if (abs(sin(val_d)).GE.0.09D0) THEN  ! ITmp=ITmp+1  ! DistFrag(ITmp)=Norm1  ! FragDist(ITmp)=JAt  ! END IF  ! END IF  ! END DO  ! IF (ITmp.EQ.0) THEN  ! ! All dihedral angles between frozen atoms are less than 5?  ! ! (or more than 175?). We have to look at other fragments :-(  ! DO I=1,NFroz  ! Jat=Frozen(I)  ! if (.NOT.DejaFait(JAt)) THEN  ! n1=JAt  ! CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)  ! CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)  ! val_d=angle_d(vx4,vy4,vz4,norm4, &  ! vx5,vy5,vz5,norm5, &  ! vx2,vy2,vz2,norm2)  ! if (abs(sin(val_d)).GE.0.09D0) THEN  ! ITmp=ITmp+1  ! DistFrag(ITmp)=Norm1  ! FragDist(ITmp)=JAt  ! END IF  ! END IF  ! END DO  ! IF (ITmp.EQ.0) THEN  ! ! All frozen atoms are in a plane... too bad  ! WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'  ! WRITE(*,*) 'For now, I do not treat this case'  ! STOP  ! END IF  ! END IF  ! I1=0  ! d=1e9  ! DO I=1,ITmp  ! IF (DistFrag(I).LE.d) THEN  ! d=DistFrag(I)  ! I1=FragDist(I)  ! END IF  ! END DO  ! ELSE !(sDihe.LE.0.09d0)  ! I1=FrozBlock(IAt,ITmp)  ! END IF !(sDihe.LE.0.09d0)  ! ! we now have the atom that is closer to the first one and that  ! ! forms a non 0/Pi dihedral angle  ! ! ind_zmat(4,1)=I1  ! ! ind_zmat(4,2)=IAt  ! ! ind_zmat(4,3)=Liaisons(Iat,1)  ! ! ind_zmat(4,4)=Liaisons(Iat,2)  ! n3=Liaisons(Iat,1)  ! n4=Liaisons(Iat,2)  ! Call add2indzmat(na,4,I1,Iat,n3,n4,ind_zmat,x,y,z)  ! Liaisons(Iat,1)=n3  ! Liaisons(Iat,2)=n4  ! DejaFait(I1)=.TRUE.  ! CaFaire(3)=I1  ! CaFaire(4)=0  ! IdxCaFaire=4  ! izm=4  ! FCaf(I1)=.TRUE.  !  !!!!!! <- PFL 28 Dec 2007   CaFaire(3)=0   IdxCaFaire=3   izm=3   CASE(1)   if (debug) WRITE(*,*) 'DBG select case I0 1, NbAtFrag=',NbAtFrag(IFrag)   ind_zmat(1,1)=IAt   ind_zmat(2,1)=Liaisons(IAt,1)   ind_zmat(2,2)=IAt   DejaFait(IAt)=.TRUE.   DejaFait(Liaisons(Iat,1))=.TRUE.   IdxCaFaire=2   CaFaire(1)=Liaisons(Iat,1)   CaFaire(2)=0   FCaf(Liaisons(Iat,1))=.TRUE.  ! PFL 28 Dec 2007 ->  ! We do NOT need a fourth atom. So we will look only for a third atom  !  !!!!  !   ! We search for a third and fourth atoms, first in the FrozBlock atoms   ! It should not be possible to have (FrozBlock(Iat,0).GT.2) and   ! iat linked to only one atom !   ! we calculate the distances between Iat and all other frozen   ! atoms of this fragment, and store only those for which   ! valence angles are not too close to 0/Pi. (limit:5?)   ITmp=0   CALL vecteur(Liaisons(Iat,1),IAt,x,y,z,vx2,vy2,vz2,norm2)  ! PFL 28 Dec 2007: As MaxL=1 I think that there is at most 2 atoms in this fragment...  ! so that the following loop is useless... this should be tested more carefully   DO I=1,NbAtFrag(IFrag)   JAt=FragAt(IFrag,I)   if (.NOT.DejaFait(JAt)) THEN   CALL vecteur(JAt,IAt,x,y,z,vx1,vy1,vz1,norm1)   if (abs(cos(angle(vx1,vy1,vz1,norm1, &   vx2,vy2,vz2,norm2))).LE.0.996d0) THEN   ITmp=ITmp+1   DistFrag(ITmp)=Norm1   FragDist(ITmp)=JAt   END IF   END IF   END DO   IF (ITMP.EQ.0) THEN   ! We have no atoms correct in this fragment, we use   ! atoms from other fragments   DO Jat=1,Na  ! DejaFait(Iat)=.TRUE. so that we do not need to test Jat/=Iat   if (.NOT.DejaFait(JAt)) THEN   CALL vecteur(JAt,IAt,x,y,z,vx1,vy1,vz1,norm1)   if (abs(cos(angle(vx1,vy1,vz1,norm1, &   vx2,vy2,vz2,norm2))).LE.0.996d0) THEN   ITmp=ITmp+1   DistFrag(ITmp)=Norm1   FragDist(ITmp)=JAt   END IF   END IF   END DO   IF (ITMP.EQ.0) THEN   WRITE(*,*) 'It seems all atoms are aligned'   WRITE(*,*) 'Case non treated for now :-( '   STOP   END IF   END IF   I1=0   d=1e9  ! PFL 28 Dec 2007: There exists some F90 intrinsics to find the smallest element of an array.  ! The following loop should be replaced by it !   DO I=1,ITmp   IF (DistFrag(I).LE.d) THEN   I1=FragDist(I)   d=DistFrag(I)   END IF   END DO   ! we now have the atom that is closer to the first one and that   ! forms a non 0/Pi valence angle   ind_zmat(3,1)=I1   ind_zmat(3,2)=IAt   ind_zmat(3,3)=Liaisons(Iat,1)   DejaFait(I1)=.TRUE.   CaFaire(2)=I1   FCaf(I1)=.TRUE.  ! PFL 28 Dec 2007 ->  ! We do NOT need a fourth atom so that the following is suppressed  !  ! ! we search for a fourth atom !  ! ! We search for a fourth atom, first in the FrozBlock atoms  ! ITmp=2  ! sDihe=0.  ! n2=IAt  ! n3=Liaisons(Iat,1)  ! n4=I1  ! CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)  ! CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)  ! CALL produit_vect(vx3,vy3,vz3, vx2,vy2,vz2,vx5,vy5,vz5,norm5)  !  ! ! We will look at the other frozen atoms  ! ! we might improve the search so as to take the atom closest to IAt  ! ITmp=0  ! DO I=1,NbAtFrag(IFrag)  ! JAt=FragAt(IFrag,I)  ! if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN  ! n1=JAt  ! CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)  ! CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)  ! val_d=angle_d(vx4,vy4,vz4,norm4, &  ! vx5,vy5,vz5,norm5, &  ! vx2,vy2,vz2,norm2)  ! if (abs(sin(val_d)).GE.0.09D0) THEN  ! ITmp=ITmp+1  ! DistFrag(ITmp)=Norm1  ! FragDist(ITmp)=JAt  ! END IF  ! END IF  ! END DO  ! IF (ITmp.EQ.0) THEN  ! ! All dihedral angles between frozen atoms are less than 5?  ! ! (or more than 175?). We have to look at other fragments :-(  ! DO I=1,NFroz  ! Jat=Frozen(I)  ! if (.NOT.DejaFait(JAt)) THEN  ! n1=JAt  ! CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)  ! CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)  ! val_d=angle_d(vx4,vy4,vz4,norm4, &  ! vx5,vy5,vz5,norm5, &  ! vx2,vy2,vz2,norm2)  ! if (abs(sin(val_d)).GE.0.09D0) THEN  ! ITmp=ITmp+1  ! DistFrag(ITmp)=Norm1  ! FragDist(ITmp)=JAt  ! END IF  ! END IF  ! END DO  ! IF (ITmp.EQ.0) THEN  ! ! All frozen atoms are in a plane... too bad  ! WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'  ! WRITE(*,*) 'For now, I do not treat this case'  ! STOP  ! END IF  ! END IF ! ITmp.EQ.0 after scaning fragment  ! I1=0  ! d=1e9  ! DO I=1,ITmp  ! IF (DistFrag(I).LE.d) THEN  ! d=DistFrag(I)  ! I1=FragDist(I)  ! END IF  ! END DO  !  ! ! we now have the atom that is closer to the first one and that  ! ! forms a non 0/Pi dihedral angle  ! ! ind_zmat(4,1)=I1  ! ! ind_zmat(4,2)=IAt  ! ! ind_zmat(4,3)=ind_zmat(2,1)  ! ! ind_zmat(4,4)=ind_zmat(3,1)  ! n3=ind_zmat(2,1)  ! n4=ind_zmat(3,1)  ! Call add2indzmat(na,4,I1,IAt,n3,n4,ind_zmat,x,y,z)  ! ind_zmat(2,1)=n3  ! ind_zmat(3,1)=n4  ! DejaFait(I1)=.TRUE.  ! CaFaire(3)=I1  ! CaFaire(4)=0  ! IdxCaFaire=4  ! izm=4  ! FCaf(I1)=.TRUE.  !!!!!!!!!!!  !  ! <- PFL 28 Dec 2007   CaFaire(3)=0   IdxCaFaire=3     CASE(0)   WRITE(*,*) 'All atoms are separated .. '   WRITE(*,*) 'this case should be treated separately !'   STOP   END SELECT   if (debug) THEN   WRITE(*,*) 'ind_zmat 1 izm=',izm   WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)   WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)   WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)   DO I=4,izm   WRITE(*,'(1X,4(1X,I5))') ind_zmat(I,1), ind_zmat(I,2), &   ind_zmat(I,3), ind_zmat(I,4)   END DO   END IF   DO I=1,izm   Idx_zmat(ind_zmat(I,1))=i   END Do   ! at least first three atoms of this fragment done...   ! we empty the 'cafaire' array before going on   IAFaire=1   DO WHILE (CaFaire(IaFaire).NE.0)   n1=CaFaire(IaFaire)   n2=ind_zmat(idx_zmat(N1),2)   if (idx_zmat(N1).EQ.2) THEN   ! We have a (small) problem: we have to add atoms linked to the 2nd   ! atom of the zmat. This is a pb because we do not know   ! which atom to use to define the dihedral angle   ! we take the third atom of the zmat   n3=ind_zmat(3,1)   ELSE   n3=ind_zmat(idx_zmat(n1),3)   END IF     FirstAt=.TRUE.   DO I=1,Liaisons(n1,0)   IAt=Liaisons(n1,I)  ! PFL 29.Aug.2008 ->  ! We dissociate the test on 'DejaFait' that indicates that this atom  ! has already been put in the Zmat  ! from the test on FCaf that indicates that this atom has been put in the  ! 'CAFaire' list that deals with identifying its connectivity.  ! Those two test might differ in some cases.   IF (.NOT.DejaFait(Iat)) THEN   izm=izm+1   if (debug) WRITE(*,*) ">1< Adding atom ",Iat,"position izm=",izm   ! ind_zmat(izm,1)=iat   ! ind_zmat(izm,2)=n1   ! ind_zmat(izm,3)=n2   ! ind_zmat(izm,4)=n3   Call add2indzmat(na,izm,iat,n1,n2,n3,ind_zmat,x,y,z)   if (FirstAt) THEN   n3=Iat   FirstAt=.FALSE.   END IF   idx_zmat(iat)=izm   DejaFait(iat)=.TRUE.   END IF   IF (.NOT.FCaf(Iat)) THEN   CaFaire(IdxCaFaire)=iat   IdxCaFaire=IdxCaFaire+1   CaFaire(IdxCaFaire)=0   FCaf(Iat)=.TRUE.   END IF  ! <- PFL 29.Aug.2008   END DO   IaFaire=IaFaire+1   END Do ! DO WHILE CaFaire   if (debug) THEN   WRITE(*,*) 'ind_zmat 2, izm=',izm   WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)   WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)   WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)   DO I=4,izm   WRITE(*,'(1X,4(1X,I5))') ind_zmat(I,1), ind_zmat(I,2), &   ind_zmat(I,3), ind_zmat(I,4)   END DO   END IF   ! We have finished putting atoms linked to the first one   ! There should not be any atom left from this fragment. We check:   ! we will add other atoms of this fragment   DO I=1,NbAtFrag(IFrag)   Iat=FragAt(IFrag,I)   if (debug) WRITE(*,*) "DBG: I,Iat,dejafait",I,Iat,DejaFait(Iat)   IF (.NOT.DejaFait(Iat)) THEN   WRITE(*,*) 'Treating atom I,Iat',I,Iat     END IF   END DO   NbAtFrag(Ifrag)=0   MaxLFrag(IFrag,1)=0     ! we start again with the rest of the molecule...   ! v 1.01 We add the fragment in the order corresponding to NbAtFrag   KMax=NbFrag-1   IF (DEBUG) WRITE(*,*) "Adding the ",Kmax," remaining fragments"   DO K=1, KMax   IFrag=0   I0=0   IAt=0   I1=0   DO I=1,NbFrag   SELECT CASE(MaxLFrag(I,1)-I0)   CASE (1:)   IFrag=I   I0=MaxLFrag(I,1)   IAt=MaxLFrag(I,2)   I1=NbAtFrag(I)   CASE (0)   if (NbAtFrag(I).GT.I1) THEN   IFrag=I   I0=MaxLFrag(I,1)   IAt=MaxLFrag(I,2)   I1=NbAtFrag(I)   END IF   END SELECT     END DO   if (debug) WRITE(*,'(1X,A,I5,A,I5,A,I5,A,I5)') 'Adding fragment:',IFrag,' with ',I0 &   ,' max links for atom',IAt,' fragment size',NbAtFrag(IFrag)     MaxLFrag(IFrag,1)=0  ! We search for the closest atoms of the previous fragments to the atom with max links   d=1e9   DO J=1,izm   Call vecteur(iat,ind_zmat(j,1),x,y,z,vx1,vy1,vz1,norm1)   if (norm1.le.d) THEN   Jat=j   d=norm1   END IF   END DO   n2=ind_zmat(jat,1)   SELECT CASE(jat)   CASE (1)   n3=ind_zmat(2,1)   n4=ind_zmat(3,1)   CASE (2)   n3=ind_zmat(1,1)   n4=ind_zmat(3,1)   CASE DEFAULT   n3=ind_zmat(jAt,2)   n4=ind_zmat(jat,3)   END SELECT   izm=izm+1   Call add2indzmat(na,izm,iat,n2,n3,n4,ind_zmat,x,y,z)   idx_zmat(iat)=izm   DejaFait(iat)=.TRUE.   IdxCaFaire=2   CaFaire(1)=iat   CaFaire(2)=0   FCaf(Iat)=.TRUE.   IaFaire=1   DO WHILE (CaFaire(IaFaire).NE.0)   n1=CaFaire(IaFaire)   n2=ind_zmat(idx_zmat(N1),2)   if (idx_zmat(N1).EQ.2) THEN   ! We have a (small) problem: we have to add atoms linked to the 2nd   ! atom of the zmat. This is a pb because we do not know   ! which atom to use to define the dihedral angle   ! we take the third atom of the zmat   n3=ind_zmat(3,1)   ELSE   n3=ind_zmat(idx_zmat(n1),3)   END IF   DO I3=1,Liaisons(n1,0)   IAt=Liaisons(n1,I3)  ! PFL 29.Aug.2008 ->  ! We dissociate the test on 'DejaFait' that indicates that this atom  ! has already been put in the Zmat  ! from the test on FCaf that indicates that this atom has been put in the  ! 'CAFaire' list that deals with identifying its connectivity.  ! Those two test might differ for a frozen atom linked to non frozen atoms.   IF (.NOT.DejaFait(Iat)) THEN   izm=izm+1   Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)   idx_zmat(iat)=izm   n3=iat   DejaFait(Iat)=.TRUE.   END IF   IF (.NOT.FCaf(Iat)) THEN   CaFaire(IdxCaFaire)=iat   IdxCaFaire=IdxCaFaire+1   CaFaire(IdxCaFaire)=0   FCaf(Iat)=.TRUE.   END IF  ! <- PFL 29.Aug.2008   END DO   IaFaire=IaFaire+1   END Do ! DO WHILE CaFaire   if (debug) THEN   WRITE(*,*) 'ind_zmat 4'   WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)   WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)   WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)   DO Ip=4,izm   WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1), ind_zmat(Ip,2), &   ind_zmat(Ip,3), ind_zmat(Ip,4)   END DO   END IF   END DO ! Loop on all fragments   ! We have ind_zmat. We calculate val_zmat :-)   if (debug) WRITE(*,*) "Calculating val_zmat"   val_zmat(1,1)=0.d0   val_zmat(1,2)=0.d0   val_zmat(1,3)=0.d0   val_zmat(2,2)=0.d0   val_zmat(2,3)=0.d0   val_zmat(3,3)=0.d0   n1=ind_zmat(2,1)   n2=ind_zmat(2,2)   CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)   val_zmat(2,1)=norm1   n1=ind_zmat(3,1)   n2=ind_zmat(3,2)   n3=ind_zmat(3,3)   CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)   val_zmat(3,1)=norm1   CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)   val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)   val_zmat(3,2)=val   DO i=4,na   n1=ind_zmat(i,1)   n2=ind_zmat(i,2)   n3=ind_zmat(i,3)   n4=ind_zmat(i,4)   if (debug) WRITE(*,*) "Doing i,n1,n2,n3,n4",i,n1,n2,n3,n4   CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)   CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)   val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)   CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)   CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)   CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,vx5,vy5,vz5,norm5)   val_d=angle_d(vx4,vy4,vz4,norm4, vx5,vy5,vz5,norm5, &   vx2,vy2,vz2,norm2)   ! write(*,11) n1,n2,norm1,n3,val,n4,val_d   !11 format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))   val_zmat(i,1)=norm1   val_zmat(i,2)=val   val_zmat(i,3)=val_d   END DO   if (debug) THEN   WRITE(*,*) 'DBG Cre_Zmat_Frag: ind_zmat'   DO I=1,na   WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)   END DO   WRITE(*,*) 'DBG Cre_Zmat_Frag: Full zmat'   DO I=1,na   WRITE(*,'(1X,I5,1X,I5,F8.4,2(1X,I5,1X,F7.2))') ind_zmat(i,1),(ind_zmat(i,j+1),val_zmat(i,j),j=1,3)   END DO   END IF   if (debugGaussian) THEN   WRITE(*,*) 'DBG Cre_Zmat_Frag: Gaussian Zmat - START'   Call zmat_g92(na,atome,ind_zmat,val_zmat)   WRITE(*,*) 'DBG Cre_Zmat_Frag: Gaussian Zmat - END'   END IF   if (debug) WRITE(*,*) "Deallocate (FragDist,Fragment, NbAtFrag,FragAt)"   DEALLOCATE(FragDist,Fragment, NbAtFrag,FragAt)   if (debug) WRITE(*,*) "Deallocate (DistFrag,Liaisons)"   DEALLOCATE(DistFrag,Liaisons)   if (debug) WRITE(*,*) "Deallocate(CaFaire,DejaFait)"   DEALLOCATE(CaFaire,DejaFait,FCaf,MaxLFrag)     if (debug) WRITE(*,*) "=============================== Exiting Calc_zmat_frag ========================"  END SUBROUTINE Calc_Zmat_frag  SUBROUTINE zmat_g92(na,atome,ind_zmat,val_zmat)  ! This subroutine comes for Cart. Slightly modified to be f90   Use Path_module, only : Nom   Use Io_module   IMPLICIT NONE   integer(KINT), INTENT(IN) :: na,atome(na)   INTEGER(KINT), INTENT(IN) :: ind_zmat(Na,5)   real(KREAL), INTENT(IN) :: val_zmat(Na,3)   character(6) :: at1, at2, at3, at4, d, a, dh   character(SCHARS), ALLOCATABLE :: tab(:) ! 3*na   character(LCHARS) :: ligne     INTEGER(KINT) :: i,n1,n2,n3,n4   ALLOCATE(tab(3*na))   DO i=1,na   IF (i .GE. 1) THEN   n1=ind_zmat(i,1)   write(at1,11) nom(atome(n1)),n1  11 format(a2,i3)   Call ecris_sb(at1,at1)   write(ligne,4) at1   END IF   IF (i .GE. 2) THEN   n2=ind_zmat(i,2)   write(at2,11) nom(atome(n2)),n2   Call ecris_sb(At2,at2)   write(d,11) 'R',i-1   Call ecris_sb(D,d)   write(ligne,4) at1,at2,d   write(tab(i-1),12) d,val_zmat(i,1)  12 format(a8,f8.4)   END IF   IF (i .GE. 3) THEN   n3=ind_zmat(i,3)   write(at3,11) nom(atome(n3)),n3   Call ecris_sb(At3,at3)   write(a,11) 'A',na+i-3   Call ecris_sb(A,A)   write(ligne,4) at1,at2,d,at3,a   write(tab(na+i-3),13) a,val_zmat(i,2)  13 format(a8,f8.3)   END IF   IF (i .GE. 4) THEN   n4=ind_zmat(i,4)   write(at4,11) nom(atome(n4)),n4   Call ecris_sb(At4,at4)   write(dh,11) 'DH',na+na+i-6   Call ecris_sb(dh,dh)   write(ligne,4) at1,at2,d,at3,a,at4,dh  4 format(7a6)   write(tab(na+na+i-6),13) dh,val_zmat(i,3)   END IF   write(IOOUT,*) TRIM(ligne)   END DO     write(IOOUT,*)   IF (na .EQ. 2) THEN   write(IOOUT,*) TRIM(tab(1))   ELSE   DO i=1,3*na-6   write(IOOUT,*) TRIM(tab(i))   END DO   END IF   write(IOOUT,*)   DEALLOCATE(Tab)   END SUBROUTINE zmat_g92   SUBROUTINE ecris_sb(inter1,inter)   character(6) :: inter,inter1   k=1   DO j=1,len(inter)   IF (inter(j:j) .NE. ' ' ) THEN   inter1(k:k)=inter(j:j)   k=k+1   END IF   END DO   inter1(k:6)=' '   END SUBROUTINE ecris_sb