Chimie Théorique » Opt'n Path

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.

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!----------------------------------------------------------------------

!  Copyright 2003-2014 Ecole Normale Sup?rieure de Lyon,

!  Centre National de la Recherche Scientifique,

!  Universit? Claude Bernard Lyon 1. All rights reserved.

!

!  This work is registered with the Agency for the Protection of Programs

!  as IDDN.FR.001.100009.000.S.P.2014.000.30625

!

!  Authors: P. Fleurat-Lessard, P. Dayal

!  Contact: optnpath@gmail.com

!

! This file is part of "Opt'n Path".

!

!  "Opt'n Path" is free software: you can redistribute it and/or modify

!  it under the terms of the GNU Affero General Public License as

!  published by the Free Software Foundation, either version 3 of the License,

!  or (at your option) any later version.

!

!  "Opt'n Path" is distributed in the hope that it will be useful,

!  but WITHOUT ANY WARRANTY; without even the implied warranty of

!

!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

!  GNU Affero General Public License for more details.

!

!  You should have received a copy of the GNU Affero General Public License

!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.

!

! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,

! for commercial licensing opportunities.

!----------------------------------------------------------------------

  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

  REAL(KREAL), PARAMETER :: LocalNCart=0.d0

  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

        ! 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