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,Pi,Massat,BMat_BakerT, &

       Nat,NCoord,XyzGeomI,NGeomI,NGeomF,UMat,NPrim,BTransInv, &

       IntCoordI,Coordinate,CurrentCoord,ScanCoord,BprimT,BBT, &

       BBT_inv,Xprimitive,XPrimitiveF,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) Call Header(" 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