SUBROUTINE Calc_baker(atome,r_cov,fact,frozen,IGeomRef)
  !     
  ! This subroutine analyses a geometry to construct the baker
  ! delocalized internal coordinates
  ! v1.0
  ! We use only one geometry
  !     
  Use Path_module, only : max_Z,NMaxL,Nom,MaxFroz,NFroz,Pi,Massat,a0,BMat_BakerT, &
       Nat,NCoord,XyzGeomI,NGeomI,NGeomF,UMat,NPrim,BTransInv, &
       IntCoordI,Coordinate,CurrentCoord,ScanCoord,BprimT,BBT, &
       BBT_inv,Xprimitive,XPrimitiveF,Symmetry_elimination,UMat_local, &
       BTransInv_local, UMatF,BTransInvF,      &
       indzmat, dzdc,renum,order,OrderInv
  ! BMat_BakerT(3*Nat,NCoord), NCoord=3*Nat or NFree=3*Nat-6-Symmetry_elimination
  ! depending upon the coordinate choice. IntCoordI(NGeomI,NCoord) where 
  ! UMat(NPrim,NCoord), NCoord number of vectors in UMat matrix, i.e. NCoord
  ! Baker coordinates. NPrim is the number of primitive internal coordinates.

  Use Io_module
  IMPLICIT NONE

  ! IGeom: Geometry index, IGeomRef: Reference Geometry.
  INTEGER(KINT), INTENT(IN) :: atome(Nat),IGeomRef
  REAL(KREAL), INTENT(IN) :: fact
  REAL(KREAL), INTENT(IN)  :: r_cov(0:Max_Z)
  ! Frozen contains the indices of frozen atoms
  INTEGER(KINT), INTENT(IN) :: Frozen(*)

  INTEGER(KINT), ALLOCATABLE :: LIAISONS(:,:) ! (Nat,0:NMaxL)
  REAL(KREAL), ALLOCATABLE :: Geom(:,:) !(3,Nat)
  ! NPrim is the number of primitive coordinates and NCoord is the number
  ! of internal coordinates. BMat is actually (NPrim,3*Nat).
  REAL(KREAL), ALLOCATABLE :: GMat(:,:) !(NPrim,NPrim)
  ! EigVec(..) contains ALL eigevectors of BMat times BprimT, NOT only Baker Coordinate vectors.
  REAL(KREAL), ALLOCATABLE :: EigVec(:,:), EigVal(:) ! EigVec(NPrim,NPrim)
  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)
  REAL(KREAL), ALLOCATABLE :: x_refGeom(:), y_refGeom(:), z_refGeom(:)
  !REAL(KREAL), ALLOCATABLE :: V(:,:) ! (3*Nat,3*Nat)
  !REAL(KREAL), ALLOCATABLE :: P(:) ! Used in Baker Coordinates: Matrix Inversion
  REAL(KREAL), ALLOCATABLE :: UMat_tmp(:,:)
  INTEGER(KINT) :: NbBonds, NbAngles, NbDihedrals, IGeom

  real(KREAL) :: vx, vy, vz, Norm
  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

  INTEGER(KINT) :: I, J, IAt, K
  INTEGER(KINT) :: Kat, Lat, L

  REAL(KREAL) :: sAngleIatIKat, sAngleIIatLat

  LOGICAL :: debug, bond, AddPrimitiveCoord

  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)
       use Path_module, only :  Pi,KINT, KREAL
       real(KREAL) ::  v1x,v1y,v1z,norm1
       real(KREAL) ::  v2x,v2y,v2z,norm2
       real(KREAL) ::  angle
     END FUNCTION ANGLE



     FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3)
       use Path_module, only :  Pi,KINT, KREAL
       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

 


     SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef)
       !     
       ! This subroutine uses the description of a list of Coordinates
       ! to compute the values of the coordinates for a given geometry.
       !
!!!!!!!!!!
       ! Input:
       !   Na: INTEGER, Number of atoms
       !  x,y,z(Na): REAL, cartesian coordinates of the considered geometry
       ! Coordinate (Pointer(ListCoord)): description of the wanted coordiantes
       ! NPrim, INTEGER: Number of coordinates to compute
       !
       ! Optional: XPrimRef(NPrim) REAL: array that contains coordinates values for
       ! a former geometry. Useful for Dihedral angles evolution...

!!!!!!!!!!!
       ! Output:
       ! XPrimimite(NPrim) REAL: array that will contain the values of the coordinates
       !
!!!!!!!!!

       Use VarTypes
       Use Io_module
       Use Path_module, only : pi

       IMPLICIT NONE

       Type (ListCoord), POINTER :: Coordinate
       INTEGER(KINT), INTENT(IN) :: Nat,NPrim
       REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat)
       REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) 
       REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim)

     END SUBROUTINE CALC_XPRIM
  END INTERFACE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!! PFL Debug
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  INTEGER(KINT) :: I1,I2,I3,I4
  LOGICAL :: debugPFL
  REAL(KREAL), ALLOCATABLE :: val_zmat(:,:)
  INTEGER(KINT), ALLOCATABLE :: indzmat0(:,:)


  debug=valid("Calc_baker")
  debugPFL=valid("BakerPFL")
  if (debug) WRITE(*,*) "============= Entering Cal_baker =============="

  IF (Nat.le.2) THEN
     WRITE(*,*) "I do not work for less than 2 atoms :-p"
     RETURN
  END IF

  IF (debug) THEN
     WRITE(*,*) "DBG Calc_baker, geom"
     DO I=1,Nat
        WRITE(*,'(1X,I5,":",I3,3(1X,F15.3))') I,atome(I),XyzGeomI(IGeomRef,1:3,I)
     END DO
  END IF




  ! we construct the list of bonds, valence angles and dihedral angles using this connectivity
  ALLOCATE(Coordinate)
  NULLIFY(Coordinate%next)
  NbBonds=0
  NbAngles=0
  NbDihedrals=0

  CurrentCoord => Coordinate

  ALLOCATE(Liaisons(Nat,0:NMaxL),Geom(3,Nat),x(Nat),y(Nat),z(Nat))
  ALLOCATE(x_refGeom(Nat),y_refGeom(Nat),z_refGeom(Nat))

  x_refGeom(1:Nat) = XyzGeomI(IGeomRef,1,1:Nat)
  y_refGeom(1:Nat) = XyzGeomI(IGeomRef,2,1:Nat)
  z_refGeom(1:Nat) = XyzGeomI(IGeomRef,3,1:Nat) 

!!!!!!!!!!!!!!!!!!!!
     !
     !   PFL Debug
     !
!!!!!!!!!!!!!!!!!!!!
     if (debugPFL) THEN
        WRITE(*,*) "DBG PFL Calc_BAKER - Calculating Zmat"
     if (.not.ALLOCATED(indzmat)) THEN
        ALLOCATE(indzmat(Nat,5))
     END IF
     IF (.NOT.ALLOCATED(DzDc)) THEN
        ALLOCATE(DzDc(3,nat,3,Nat))
     END IF
     IF (.NOT.ALLOCATED(Val_zmat)) THEN
        ALLOCATE(val_zmat(nat,3))
     END IF
     IF (.NOT.ALLOCATED(indzmat0)) THEN
        ALLOCATE(indzmat0(nat,5))
     END IF
! We construct a Zmat
          IF (Frozen(1).NE.0) THEN 
           Renum=.TRUE.
!!! remplcaer indzmat par indzmat0 puis renumerotation
           call Calc_zmat_const_frag(Nat,atome,IndZmat0,val_zmat, &
             x_refGeom,y_refGeom,z_refGeom, &
                r_cov,fact,frozen)
        ELSE
           !no zmat... we create our own 
           call Calc_zmat_frag(Nat,atome,IndZmat0,val_zmat, &
             x_refGeom,y_refGeom,z_refGeom, &
                r_cov,fact)
        END IF

         DO I=1,Nat
           Order(IndZmat0(I,1))=I
           OrderInv(I)=IndZmat0(I,1)
        END DO
        IndZmat(1,1)=Order(IndZmat0(1,1))
        IndZmat(2,1)=Order(IndZmat0(2,1))
        IndZmat(2,2)=Order(IndZmat0(2,2))
        IndZmat(3,1)=Order(IndZmat0(3,1))
        IndZmat(3,2)=Order(IndZmat0(3,2))
        IndZmat(3,3)=Order(IndZmat0(3,3))
        DO I=4,Nat
           DO J=1,4
              IndZmat(I,J)=Order(IndZmat0(I,J))
           END DO
        end do


        WRITE(*,*) "DBG PFL CalcBaker, IndZmat0"
        DO I=1,Nat
            WRITE(*,*) I, indZmat0(I,:)
        end do

        WRITE(*,*) "DBG PFL CalcBaker, IndZmat"
        DO I=1,Nat
            WRITE(*,*) I, indZmat(I,:)
        end do


     DO I=1,Nat
        I1=indzmat0(I,1)
        I2=indzmat0(I,2)
        I3=indzmat0(I,3)
        I4=indzmat0(I,4)
! Adding BOND between at1 and at2
      WRITE(*,*)  "DBG Indzmat",I,I1,I2,I3,I4
         if (I2.NE.0) THEN
           NbBonds=NbBonds+1
    ! values of the coordinate (bond) is calculated for the reference geometry
    ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.)
    ! are determined using all geometries. vecteur is defined in bib_oper.f
                  Call vecteur(I1,I2,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
!                  CurrentCoord%value=Norm2
                  WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Bond",NbBonds," Norm2,val_zmat",Norm2,val_zmat(I,1)
                  CurrentCoord%value=val_zmat(I,1)
                  CurrentCoord%SignDihedral = 0
                  CurrentCoord%At1=I1
                  CurrentCoord%At2=I2
                  CurrentCoord%Type="BOND"
                  IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
                     ALLOCATE(CurrentCoord%next)
                     NULLIFY(CurrentCoord%next%next)
                  END IF
                  CurrentCoord => CurrentCoord%next
        END IF !IE.NE.0
       if (I3.NE.0) THEN

                       NbAngles=NbAngles+1
                       ! values of the coordinate (bond) is calculated for the reference geometry
                       ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are 
                       ! determined using all geometries.
                       Call vecteur(i2,I3,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
                       Call vecteur(i2,I1,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
!                       CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
                  WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Angle",NbAngles," angle,val_zmat", &
                       angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2),val_zmat(I,2)
                       CurrentCoord%value=val_zmat(I,2)*Pi/180.
                       CurrentCoord%SignDihedral = 0
                       CurrentCoord%At1=I1
                       CurrentCoord%At2=I2
                       CurrentCoord%At3=I3
                       CurrentCoord%Type="ANGLE"
                       IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
                          ALLOCATE(CurrentCoord%next)
                          NULLIFY(CurrentCoord%next%next)
                       END IF
                       CurrentCoord => CurrentCoord%next
                    END IF ! matches IF (I3.NE.0)
         if (I4.NE.0) THEN
                       NbDihedrals=NbDihedrals+1
              Call vecteur(I2,I1,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
              Call vecteur(I2,I3,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
              Call vecteur(I3,I4,x_refGeom,y_refGeom,z_refGeom,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)
!                       CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
!                            vx2,vy2,vz2,norm2)*Pi/180.
                  WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Dihedral",NbDihedrals,    &
                       " angle_d,val_zmat", &
                       angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
                       vx2,vy2,vz2,norm2) &
                       ,val_zmat(I,3)

                       CurrentCoord%value = val_zmat(I,3)*Pi/180.
                       CurrentCoord%SignDihedral = int(sign(1.D0,val_zmat(I,3)))
                       CurrentCoord%At1=I1
                       CurrentCoord%At2=I2
                       CurrentCoord%At3=I3
                       CurrentCoord%At4=I4
                       CurrentCoord%Type="DIHEDRAL"
                       IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
                          ALLOCATE(CurrentCoord%next)
                          NULLIFY(CurrentCoord%next%next)
                       END IF
                       CurrentCoord => CurrentCoord%next
                    END IF ! matches IF (I4.NE.0)
              END DO ! matches DO I=1,Nat
      deallocate(indzmat0)

  ELSE !! if (debugPFL)

  DO IGeom=1, NGeomI
     !IGeom=1 ! need when using only first geometry.
     x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
     y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
     z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) 

     ! First step: we calculate the connectivity
     Liaisons=0

     Call CalcCnct(Nat,atome,x,y,z,LIAISONS,r_cov,fact)

     IF (debug) THEN 
        WRITE(*,*) 'Connectivity done'
        DO I=1,Nat
           WRITE(*,'(1X,I5,":",I3,"=",15I4)') I,(Liaisons(I,J),J=0,NMaxL)
           WRITE(*,'(1X,I5,":",I3,"=",15I4)') I,(Liaisons(I,:))
        END DO
        DO I=1,Nat
           WRITE(*,'(1X,I5,":",I3," r_cov=",F15.6)') I,atome(I),r_cov(atome(I))
        END DO
     END IF

     DO I=1,Nat
        IF (Liaisons(I,0).NE.0) THEN
           DO J=1,Liaisons(I,0)
              ! We add the bond
              Iat=Liaisons(I,J)
              IF (Iat.GT.I) THEN
                 ! Checking the uniqueness of bond.
                 ! Duplicate bonds should not be added.
                 AddPrimitiveCoord=.True.
                 ScanCoord=>Coordinate
                 DO WHILE (Associated(ScanCoord%next))
                    IF (ScanCoord%Type .EQ. 'BOND') THEN
                       IF (ScanCoord%At1 .EQ. Iat) THEN
                          IF (ScanCoord%At2 .EQ. I) THEN
                             AddPrimitiveCoord=.False.
                          END IF
                       END IF
                    END IF
                    ScanCoord => ScanCoord%next
                 END DO

                 IF (AddPrimitiveCoord) THEN
                    NbBonds=NbBonds+1
                    ! values of the coordinate (bond) is calculated for the reference geometry
                    ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.)
                    ! are determined using all geometries. vecteur is defined in bib_oper.f
                    Call vecteur(I,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
                    CurrentCoord%value=Norm2
                    CurrentCoord%SignDihedral = 0
                    CurrentCoord%At1=Iat
                    CurrentCoord%At2=I
                    CurrentCoord%Type="BOND"
                    IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
                       ALLOCATE(CurrentCoord%next)
                       NULLIFY(CurrentCoord%next%next)
                    END IF
                    CurrentCoord => CurrentCoord%next
                 END IF ! matches IF (AddPrimitiveCoord) THEN
              END IF ! matches IF (Iat.GT.I) THEN

              !           Call Annul(Liaisons,IAt,I,Nat,NMaxL)
              !           Liaisons(Iat,0)=Liaisons(Iat,0)-1

              ! We add the valence angles
              DO K=J+1,Liaisons(I,0)
                 Kat=Liaisons(I,K)
                 IF (Kat.GT.I) THEN
                    ! Checking the uniqueness of valence angle.
                    ! Duplicate bonds should not be added.
                    AddPrimitiveCoord=.True.
                    ScanCoord=>Coordinate
                    DO WHILE (Associated(ScanCoord%next))
                       IF (ScanCoord%Type .EQ. 'ANGLE') THEN
                          IF (ScanCoord%At1 .EQ. Iat) THEN
                             IF (ScanCoord%At2 .EQ. I) THEN
                                IF (ScanCoord%At3 .EQ. Kat) THEN
                                   AddPrimitiveCoord=.False.
                                END IF
                             END IF
                          END IF
                       END IF
                       ScanCoord => ScanCoord%next
                    END DO ! matches DO WHILE (Associated(ScanCoord%next))

                    IF (AddPrimitiveCoord) THEN      
                       IF (debug) WRITE(*,*) "Adding angle:",Iat,I,Kat
                       NbAngles=NbAngles+1
                       ! values of the coordinate (bond) is calculated for the reference geometry
                       ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are 
                       ! determined using all geometries.
                       Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
                       Call vecteur(i,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
                       CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
                       CurrentCoord%SignDihedral = 0
                       sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.))
                       CurrentCoord%At1=Iat
                       CurrentCoord%At2=I
                       CurrentCoord%At3=KAt
                       CurrentCoord%Type="ANGLE"
                       IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
                          ALLOCATE(CurrentCoord%next)
                          NULLIFY(CurrentCoord%next%next)
                       END IF
                       CurrentCoord => CurrentCoord%next
                    END IF ! matches IF (AddPrimitiveCoord) THEN
                 ELSE ! matches IF (Kat.GT.I) THEN
                    ! Kat is less than I. If there is a bond between I and Kat, then this angle
                    ! has already been taken into account.
                    if (debug) WRITE(*,*) "Testing angle:",Iat,I,Kat
                    Bond=.FALSE.
                    DO L=1,Liaisons(Iat,0)
                       IF (Liaisons(Iat,L).EQ.Kat) Bond=.TRUE.
                    END DO
                    IF (.NOT.Bond) THEN
                       IF (debug) WRITE(*,*) "Not Bond Adding angle:",Iat,I,Kat

                       ! Checking the uniqueness of valence angle.
                       ! Duplicate bonds should not be added.        
                       AddPrimitiveCoord=.True.
                       ScanCoord=>Coordinate
                       DO WHILE (Associated(ScanCoord%next))
                          IF (ScanCoord%Type .EQ. 'ANGLE') THEN
                             IF (ScanCoord%At1 .EQ. Iat) THEN
                                IF (ScanCoord%At2 .EQ. I) THEN
                                   IF (ScanCoord%At3 .EQ. Kat) THEN
                                      AddPrimitiveCoord=.False.
                                   END IF
                                END IF
                             END IF
                          END IF
                          ScanCoord => ScanCoord%next
                       END DO

                       IF (AddPrimitiveCoord) THEN  
                          NbAngles=NbAngles+1
                          ! values of the coordinate (bond) is calculated for the reference geometry
                          ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) 
                          ! are determined using all geometries.
                          Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
                          Call vecteur(i,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
                          CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
                           CurrentCoord%SignDihedral = 0
                          CurrentCoord%At1=Iat
                          CurrentCoord%At2=I
                          CurrentCoord%At3=KAt
                          CurrentCoord%Type="ANGLE"
                          IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
                             ALLOCATE(CurrentCoord%next)
                             NULLIFY(CurrentCoord%next%next)
                          END IF
                          CurrentCoord => CurrentCoord%next
                       END IF  ! matches IF (AddPrimitiveCoord) THEN
                    END IF ! matches IF (.NOT.Bond) THEN
                 END IF ! matches IF (Kat.GT.I) THEN, after else
              END DO ! matches DO K=J+1,Liaisons(I,0)
           END DO ! matches DO J=1,Liaisons(I,0)
        END IF ! matches IF (Liaisons(I,0).NE.0) THEN

        ! We add dihedral angles. We do not concentrate on necessary angles, ie those that are
        ! needed to build the molecule. Instead we collect all of them, and the choice of the best
        ! ones will be done when diagonalizing the Baker matrix.
        IF (Liaisons(I,0).GE.3) Then
           DO J=1,Liaisons(I,0)
              Iat=Liaisons(I,J)
              ! values of the coordinate (bond) is calculated for the reference geometry
              ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are 
              ! determined using all geometries.
              Call vecteur(Iat,i,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)

              DO K=J+1,Liaisons(I,0)
                 Kat=Liaisons(I,K)
                 Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
                 sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,-vx2,-vy2,-vz2,Norm2)*Pi/180.))

                 DO L=K+1,Liaisons(I,0)
                    Lat=Liaisons(I,L)
                    if (debug) WRITE(*,*) "0-Dihe Kat,I,Iat:",Kat,I,Iat,angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)
                    Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3)
                    sAngleIIatLat=abs(sin(angle(vx3,vy3,vz3,Norm3,vx2,vy2,vz2,Norm2)*Pi/180.))
                    if (debug) WRITE(*,*) "0-Dihe I,Iat,Lat:",angle(vx3,vy3,vz3,Norm3,vx2,vy2,vz2,Norm2)
                    IF ((sAngleIatIKat.ge.0.01d0).AND.(sAngleIIatLat.ge.0.01d0)) THEN
                       if (debug) WRITE(*,*) "Adding dihedral:",Kat,I,Iat,Lat

                       ! Checking for the uniqueness of valence angle.
                       ! Duplicate bonds should not be added.
                       AddPrimitiveCoord=.True.
                       ScanCoord=>Coordinate
                       DO WHILE (Associated(ScanCoord%next))
                          IF (ScanCoord%Type .EQ. 'DIHEDRAL') THEN
                             IF (ScanCoord%At1 .EQ. Kat) THEN
                                IF (ScanCoord%At2 .EQ. I) THEN
                                   IF (ScanCoord%At3 .EQ. Iat) THEN
                                      IF (ScanCoord%At4 .EQ. Lat) THEN
                                         AddPrimitiveCoord=.False.
                                      END IF
                                   END IF
                                END IF
                             END IF
                          END IF
                          ScanCoord => ScanCoord%next
                       END DO ! matches DO WHILE (Associated(ScanCoord%next))

                       ! IF (AddPrimitiveCoord) THEN
                       !   NbDihedrals=NbDihedrals+1
                       !  CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
                            !      vx5,vy5,vz5,norm5)
                       !CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
                            !    vx4,vy4,vz4,norm4)
                       !         CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
                       !             vx2,vy2,vz2,norm2)*Pi/180.
                       !          CurrentCoord%SignDihedral = 0
                       !       CurrentCoord%At1=Kat
                       !      CurrentCoord%At2=I
                       !     CurrentCoord%At3=IAt
                       !    CurrentCoord%At4=LAt
                       !WRITE(*,'(1X,":(dihed 1)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') &
                       !        ScanCoord%At1,ScanCoord%At2,ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi
                       !         CurrentCoord%Type="DIHEDRAL"
                       !        IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
                       !          ALLOCATE(CurrentCoord%next)
                       !         NULLIFY(CurrentCoord%next%next)
                       !     END IF
                       !    CurrentCoord => CurrentCoord%next
                       !END IF ! matches IF (AddPrimitiveCoord) THEN

                    END IF ! matches IF ((sAngleIatIKat.ge.0.01d0).AND.(sAngleIIatLat.ge.0.01d0)) THEN
                 END DO ! matches DO L=K+1,Liaisons(I,0)
              END DO ! matches DO K=J+1,Liaisons(I,0)
           END DO ! matches DO J=1,Liaisons(I,0)
        END IF !matches  IF (Liaisons(I,0).GE.3) Then

        ! we now add other dihedrals
        DO J=1,Liaisons(I,0)
           Iat=Liaisons(I,J)
           If (Iat.LE.I) Cycle
           ! values of the coordinate (bond) is calculated for the reference geometry
           ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are 
           ! determined using all geometries.
           Call vecteur(Iat,i,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)

           DO K=1,Liaisons(I,0)
              IF (Liaisons(I,K).EQ.Iat) Cycle 
              Kat=Liaisons(I,K)
              Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
              sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,-vx2,-vy2,-vz2,Norm2)*Pi/180.))

              DO L=1,Liaisons(Iat,0)
                 Lat=Liaisons(Iat,L)
                 IF (Lat.eq.I) Cycle

                 if (debug) WRITE(*,*) "Dihe Kat,I,Iat:",angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2),sAngleIatIKat
                 Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx,vy,vz,Norm)

                 sAngleIIatLat=abs(sin(angle(vx,vy,vz,Norm,vx2,vy2,vz2,Norm2)*Pi/180.))
                 if (debug) WRITE(*,*) "Dihe I,Iat,Lat:",angle(vx,vy,vz,Norm,vx2,vy2,vz2,Norm2)
                 ! we add the dihedral angle Kat-I-Iat-Lat
                 if (debug) WRITE(*,*) "Trying dihedral:",Kat,I,Iat,Lat,sAngleIatIKat,sAngleIIatLat
                 IF ((Lat.NE.Kat).AND.(Lat.Ne.I).AND.(sAngleIatIKat.ge.0.01d0) &
                      .AND.(sAngleIIatLat.ge.0.01d0)) THEN
                    if (debug) WRITE(*,*) "Adding dihedral:",Kat,I,Iat,Lat

                    ! Checking for the uniqueness of valence angle.
                    ! Duplicate bonds should not be added.
                    AddPrimitiveCoord=.True.
                    ScanCoord=>Coordinate
                    DO WHILE (Associated(ScanCoord%next))
                       IF (ScanCoord%Type .EQ. 'DIHEDRAL') THEN
                          IF (ScanCoord%At1 .EQ. Kat) THEN
                             IF (ScanCoord%At2 .EQ. I) THEN
                                IF (ScanCoord%At3 .EQ. Iat) THEN
                                   IF (ScanCoord%At4 .EQ. Lat) THEN
                                      AddPrimitiveCoord=.False.
                                   END IF
                                END IF
                             END IF
                          END IF
                       END IF
                       ScanCoord => ScanCoord%next
                    END DO

                    IF (AddPrimitiveCoord) THEN
                       NbDihedrals=NbDihedrals+1
                       Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3)
                       CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
                            vx5,vy5,vz5,norm5)
                       CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
                            vx4,vy4,vz4,norm4)
                       CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
                            vx2,vy2,vz2,norm2)*Pi/180.
                       CurrentCoord%At1=Kat
                       CurrentCoord%At2=I
                       CurrentCoord%At3=IAt
                       CurrentCoord%At4=LAt
                       CurrentCoord%Type="DIHEDRAL"
                       CurrentCoord%SignDihedral=int(sign(1.d0,angle_d(vx4,vy4,vz4,norm4,   &
                            vx5,vy5,vz5,norm5, vx2,vy2,vz2,norm2)))
                       WRITE(*,'(1X,":(dihed 2)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') ScanCoord%At1,&
                            ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi
                       IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
                          ALLOCATE(CurrentCoord%next)
                          NULLIFY(CurrentCoord%next%next)
                       END IF
                       CurrentCoord => CurrentCoord%next
                    END IF ! matches IF (AddPrimitiveCoord) THEN
                 END IF ! matches IF ((Lat.NE.Kat).AND.(Lat.Ne.I).AND.(sAngleIatIKat.ge.0.01d0) .....
              END DO ! matches DO L=1,Liaisons(Iat,0)
           END DO ! matches DO K=1,Liaisons(I,0)
        END DO ! matches DO J=1,Liaisons(I,0)
     END DO ! end of loop over all atoms, DO I=1, Nat
  END DO ! matches DO IGeom=1, NGeomI



 END IF ! if (debugPFL)
  DEALLOCATE(Liaisons)

  WRITE(*,*) "I have found"
  WRITE(*,*) NbBonds, " bonds"
  WRITE(*,*) NbAngles, " angles"
  WRITE(*,*) NbDihedrals, " dihedrals"

  WRITE(*,*) 'All in all...'
  ScanCoord=>Coordinate
  I=0
  DO WHILE (Associated(ScanCoord%next))
     I=I+1
     SELECT CASE (ScanCoord%Type)
     CASE ('BOND')
        !WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,ScanCoord%Value
        WRITE(*,'(1X,I3,":(bond )",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,ScanCoord%Value
     CASE ('ANGLE')
        !WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, &
        !    ScanCoord%At2, ScanCoord%At3,ScanCoord%Value*180./Pi
        WRITE(*,'(1X,I3,":(angle)",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, &
             ScanCoord%At2, ScanCoord%At3,ScanCoord%Value*180./Pi   
     CASE ('DIHEDRAL')
        !WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,&
        !    ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi
        WRITE(*,'(1X,I3,":(dihed)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,&
             ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi
     END SELECT
     ScanCoord => ScanCoord%next
  END DO

  ! NPrim is the number of primitive coordinates.
  NPrim=NbBonds+NbAngles+NbDihedrals

  ! Now calculating values of all primitive bonds for all geometries:
  ALLOCATE(Xprimitive(NgeomI,NPrim),XPrimitiveF(NGeomF,NPrim))
  ALLOCATE(UMatF(NGeomF,NPrim,NCoord))
  ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat))

! We calculate the Primitive values for the first geometry, this plays a special role
! as it will serve as a reference for other geometries !

  x(1:Nat) = XyzGeomI(1,1,1:Nat)
  y(1:Nat) = XyzGeomI(1,2,1:Nat)
  z(1:Nat) = XyzGeomI(1,3,1:Nat)

  Call Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive(1,:))

  DO IGeom=2, NGeomI

        x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
        y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
        z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)

        Call Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive(IGeom,:),Xprimitive(1,:))

!      ScanCoord=>Coordinate
!      I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
!      DO WHILE (Associated(ScanCoord%next))
!         I=I+1
!         SELECT CASE (ScanCoord%Type)
!         CASE ('BOND')
!            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)
!            Xprimitive(IGeom,I) = Norm2  
!         CASE ('ANGLE')
!            Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx1,vy1,vz1,Norm1)
!            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)
!            Xprimitive(IGeom,I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
!         CASE ('DIHEDRAL')


!            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1)
!            Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2)
!            Call vecteur(ScanCoord%At3,ScanCoord%At4,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)

!                  DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
!                          vx2,vy2,vz2,norm2)
!                  Xprimitive(IGeom,I) = DiheTmp*Pi/180.
! ! We treat large dihedral angles differently as +180=-180 mathematically and physically
! ! but this causes lots of troubles when converting baker to cart.
! ! So we ensure that large dihedral angles always have the same sign
!                     if (abs(ScanCoord%SignDihedral).NE.1) THEN
!                         ScanCoord%SignDihedral=Int(Sign(1.D0,DiheTmp))
!                     ELSE
!                         If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN
!                           Xprimitive(IGeom,I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi
!                        END IF
!                   END IF
! ! Another test... will all this be consistent ??? 
! ! We want the shortest path, so we check that the change in dihedral angles is less than Pi:
!                   IF (IGeom.GT.1) THEN
!                      IF (abs(Xprimitive(IGeom,I)-XPrimitive(IGeom-1,I)).GE.Pi) THEN
!                         if ((Xprimitive(IGeom,I)-XPrimitive(IGeom-1,I)).GE.Pi) THEN
!                            Xprimitive(IGeom,I)=Xprimitive(IGeom,I)-2.d0*Pi
!                         ELSE
!                            Xprimitive(IGeom,I)=Xprimitive(IGeom,I)+2.d0*Pi
!                         END IF
!                      END IF
!                   END IF
!         END SELECT
!         ScanCoord => ScanCoord%next
!      END DO ! matches DO WHILE (Associated(ScanCoord%next))
  END DO ! matches DO IGeom=1, NGeomI

  IF (debug) THEN
     WRITE(*,*) "DBG Calc_Baker Xprimitives for all geometries"
  DO IGeom=1, NGeomI
     WRITE(*,*) "Geometry ",IGeom
     ScanCoord=>Coordinate
     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
     DO WHILE (Associated(ScanCoord%next))
        I=I+1
        SELECT CASE (ScanCoord%Type)
        CASE ('BOND')
           WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive(IGeom,I)
        CASE ('ANGLE')
           WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, &
            ScanCoord%At2, ScanCoord%At3,Xprimitive(IGeom,I)*180./Pi
        CASE ('DIHEDRAL')
           WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,&
           ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive(IGeom,I)*180./Pi
        END SELECT
        ScanCoord => ScanCoord%next
     END DO ! matches DO WHILE (Associated(ScanCoord%next))
  END DO ! matches DO IGeom=1, NGeomI
  END IF
  ! Here reference geometry is assigned to Geom. Earlier x,y,z(Nat) were assigned from XyzGeomI(...)
  ! before the call of this (Calc_baker) subroutine and passed to this subroutine as arguments of the
  ! subroutine.
! PFL 19 Feb 2008 ->
! In order to reproduce the B matrix in the Baker article, one has
! to divide Geom by a0. 
! However, this is inconsitent with our way of storing geometries
! so we do not do it !

  Geom(1,:)=XyzGeomI(IGeomRef,1,1:Nat) ! XyzGeomI(NGeomI,3,Nat)
  Geom(2,:)=XyzGeomI(IGeomRef,2,1:Nat)
  Geom(3,:)=XyzGeomI(IGeomRef,3,1:Nat)

! <- PFL 19 Feb 2008

  ! We now have to compute dq/dx...
  ALLOCATE(BprimT(3*Nat,NPrim))
  BprimT=0.d0

  ScanCoord=>Coordinate
  I=0
  DO WHILE (Associated(ScanCoord%next))
     I=I+1
     SELECT CASE (ScanCoord%Type)
     CASE ('BOND') 
        CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2, &
             Geom,BprimT(1,I))
     CASE ('ANGLE')
        CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2, &
             ScanCoord%At3,Geom,BprimT(1,I))
     CASE ('DIHEDRAL')
        CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2, &
             ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I))
!        CALL CONSTRAINTS_TORSION_DER(Nat,ScanCoord%At1,ScanCoord%AT2, &
!             ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I))

     END SELECT
     ScanCoord => ScanCoord%next
  END DO


!!!!!!!!!!!!!!!!!!!!
     !
     !   PFL Debug
     !
!!!!!!!!!!!!!!!!!!!!

     if (debugPFL) THEN
        WRITe(*,*) "DBG PFL Calc_Baker ====="
     DzDc=0.d0
     WRITE(*,*) "PFL DBG, DzdC"
  do I=1,Nat
  Geom(1,i)=XyzGeomI(IGeomRef,1,orderInv(i)) ! XyzGeomI(NGeomI,3,Nat)
  Geom(2,i)=XyzGeomI(IGeomRef,2,OrderInv(i))
  Geom(3,i)=XyzGeomI(IGeomRef,3,OrderInv(i))
  WRITE(*,*) I,OrderInv(I),Geom(:,I)
    end do
     CALL CalcBMat_int (nat, Geom, indzmat, dzdc, massat,atome)     

     WRITE(*,*) "PFL DBG, BPrimT"
     DO I=1,NPrim
        WRITE(*,'(50(1X,F12.6))') BPrimT(:,I)
     END DO
        WRITe(*,*) "DBG PFL Calc_Baker ====="
     END IF

!!!!!!!!!!!!!!!!!!!!
     !
     !   PFL Debug
     !
!!!!!!!!!!!!!!!!!!!!

  DEALLOCATE(Geom)

  ! BprimT(3*Nat,NPrim)
  ! We now compute G=B(BT) matrix
  ALLOCATE(Gmat(NPrim,NPrim))
  GMat=0.d0
  DO I=1,NPrim
     DO J=1,3*Nat
        GMat(:,I)=Gmat(:,I)+BprimT(J,:)*BprimT(J,I) !*1.d0/mass(atome(int(K/3.d0)))
     END DO
  END DO

!!!!!!!!!!!!!!!!!!!!
     !
     !   PFL Debug --->
     !
!!!!!!!!!!!!!!!!!!!!
    if (debugPFL) THEN
  GMat=0.d0
  DO I=1,NPrim
        GMat(I,I)=1.
  END DO
   END IF

!!!!!!!!!!!!!!!!!!!!
     !
     !  <--- PFL Debug
     !
!!!!!!!!!!!!!!!!!!!!

  !DO I=1,NPrim
  !  WRITE(*,'(20(1X,F6.3))') GMat(1:min(20,NPrim),I)
  !END DO

  ! We diagonalize G
  ALLOCATE(EigVal(NPrim),EigVec(NPrim,NPrim))
  Call Jacobi(GMat,NPrim,EigVal,EigVec,NPrim)
  Call Trie(NPrim,EigVal,EigVec,NPrim)
  DO I=1,NPrim
     WRITE(*,'(1X,"Vector ",I3,": e=",F8.3)') I,EigVal(i)
     WRITE(*,'(20(1X,F8.4))') EigVec(1:min(20,NPrim),I)
  END DO

  !DO I=1,NPrim
  ! WRITE(*,'(1X,"Vector ",I3,(20(F8.3)))') I,EigVec(1:min(20,NPrim),I)
  ! END DO

  !DO I=1,NPrim
  !  WRITE(*,'(1X,"Vector ",I3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,&
  ! F8.3,F8.3,F8.3,F8.3,F8.3)') I,EigVec(1,I),EigVec(4,I),EigVec(6,I),EigVec(14,I),EigVec(16,I),&
  !EigVec(2,I),EigVec(3,I),EigVec(5,I),EigVec(12,I),EigVec(13,I),EigVec(15,I),EigVec(7,I),&
       !EigVec(8,I),EigVec(9,I),EigVec(10,I),EigVec(11,I),EigVec(17,I)
  !END DO

  ! We construct the new DzDc(b=baker) matrix
  ! UMat is nonredundant vector set, i.e. set of eigenvectors of BB^T corresponding
  ! to eigenvalues > zero.
  ALLOCATE(UMat(NPrim,NCoord))
  ALLOCATE(UMat_local(NPrim,NCoord))
  ALLOCATE(UMat_tmp(NPrim,3*Nat-6))
  ! BMat_BakerT(3*Nat,NCoord), allocated in Path.f90, 
  ! NCoord=3*Nat-6
  BMat_BakerT = 0.d0
  J=0
  UMat=0.d0
  UMat_tmp=0.d0
  DO I=1,NPrim
     IF (abs(eigval(I))>=1e-6) THEN
        J=J+1
        DO K=1,NPrim
           !DzDcb(:,J)=DzDcb(:,J)+Eigvec(K,I)*BprimT(:,K) ! original
           ! BprimT is transpose of B^prim.
           ! B = UMat^T B^prim, B^T = (B^prim)^T UMat
           BMat_BakerT(:,J)=BMat_BakerT(:,J)+BprimT(:,K)*Eigvec(K,I)
           !Dzdcb(:,J)=DzDcb(:,J)+Eigvec(K,:)*BprimT(I,K)       
        END DO
        IF(J .GT. 3*Nat-6) THEN
           WRITE(*,*) 'Number of vectors in Eigvec with eigval .GT. 1e-6(=UMat) (=' &
                ,J,') exceeded 3*Nat-6=',3*Nat-6, &
                'Stopping the calculation.'
           STOP
        END IF
        UMat_tmp(:,J) =  Eigvec(:,I)
     END IF
  END DO

  if (debug) THEN
     WRITE(*,*) "DBG Calc_Baker: UMat"
     DO I=1,3*Nat-6
        WRITE(*,'(50(1X,F12.8))') UMat_tmp(:,I)
     END DO
  END IF

  ! THIS IS TO BE CORRECTED:
  UMat = UMat_tmp
  !J=0
  !DO I=1, 3*Nat-6
  !  IF ((I .NE. 6) .AND. (I .NE. 7) .AND. (I .NE. 9)) Then
  !   J=J+1
  !  Print *, 'J=', J, ' I=', I
  ! UMat(:,J) = UMat_tmp(:,I)
  !END IF
  !END DO
  ! BMat_BakerT=0.d0
  !DO I=1,NCoord
  !  DO J=1,NPrim
  ! B = UMat^T B^prim, B^T = (B^prim)^T UMat
  !   BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat(J,I)
  !    END DO
  !END DO
  ! TO BE CORRECTED, ENDS HERE.

  IntCoordI=0.d0
  DO IGeom=1, NGeomI
     DO I=1, NPrim
        ! Transpose of UMat is needed below, that is why UMat(I,:).
        IntCoordI(IGeom,:) = IntCoordI(IGeom,:) + UMat(I,:)*Xprimitive(IGeom,I)
     END DO
  END DO

  if (debug) THEN
     WRITE(*,*) "DBG Calc_Baker IntCoordI"
     DO IGeom=1, NGeomI
        WRITE(*,'(1X,I5,":",30(1X,F10.6))') IGeom,IntCoordI(IGeom,:)
     END DO
  END IF


  ! the following might be used to get rid of some coordinates that 
  ! do not respect the symmetry of the problem.

  ! ! We look at the group symmetry of the molecule. Nothing complicate here
  ! ! we just look for planar symmetry.
  ! ! We first place it into the standard orientation
  !   Call Std_ori(Nat,x,y,z,a,b,c,massat)

  ! ! Is the molecule planar
  !   DO I=1,Nat
  !      WRITE(*,*) nom(atome(i)),x(i),y(i),z(i)
  !      END DO
  DEALLOCATE(BprimT,GMat,EigVal,EigVec) ! BprimT is B^prim^T
  ! Calculation of BTransInv starts here:
  ! Calculation of BBT(3*Nat-6,3*Nat-6)=BB^T:
  ! BMat_BakerT(3*Nat,NCoord) is Transpose of B = UMat^TB^prim
  ALLOCATE(BBT(NCoord,NCoord))
  BBT = 0.d0
  DO I=1, NCoord
     DO J=1, 3*Nat
        ! BBT(:,I) forms BB^T
        BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I)
     END DO
  END DO

  !ALLOCATE(V(3*Nat,3*Nat))
  !Call GINVSE(BBT,3*Nat,3*Nat,BBT_inv,3*Nat,3*Nat,3*Nat,1e-6,V,3*Nat,FAIL,NOUT)
  !DEALLOCATE(V)
  ALLOCATE(BBT_inv(NCoord,NCoord))
  ! GenInv is in Mat_util.f90
  Call GenInv(NCoord,BBT,BBT_inv,NCoord)

  ! Calculation of (B^T)^-1 = (BB^T)^-1B:
  !ALLOCATE(BTransInv(3*Nat-6,3*Nat))    ! allocated in Path.f90
  BTransInv = 0.d0
  DO I=1, 3*Nat
     DO J=1, NCoord
        BTransInv(:,I) = BTransInv(:,I) + BBT_inv(:,J)*BMat_BakerT(I,J)
     END DO
  END DO
  BTransInv_local = BTransInv
  UMat_local = UMat
  DEALLOCATE(BBT,BBT_inv,UMat_tmp)
  DEALLOCATE(x_refGeom,y_refGeom,z_refGeom,x,y,z)
  IF (debug) WRITE(*,*) "DBG Calc_baker over."
END SUBROUTINE Calc_baker
