Chimie Théorique » Opt'n Path

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