!================================================================
! Calculation of the connectivity.
!
! Input:
! Na: (INTEGER) Number of atoms in the system
! Atome: (INTEGER) Number of mass of the atoms.
!                  Needed to find their covalent raii.
!  x(Nat), y(nat), z(nat) : (REAL) Cartesian coordinates of the system.
!  LIAISONS: (INTEGER Na,0:NMaxL) Array containing the connectivity.
!            That is Liaisons(i,:) contains the atoms linked to i.
!            Liaisons(i,0) contains the number of atoms linked to i.
! r_cov(0:Max_Z) (REAL) Covalent radii
! Fact (REAL) Factor used to determined if two atoms are linked :
!            i and j are linked if (dist(i,j)*fact <= rcov(i)+rcov(j)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Taken from Path_Module
! NMaxL (INTEGER): Maximum number of bonds for one atome
! Max_Z (INTEGER): Last atom for which I have the atomic of mass and covalent radii.
! Nom(Nat) (STRING): Name of the atoms.
! Prog (SCHARS): Name of the program used to calculate energies and gradients
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! v1.0 (c) PFL
! v1.1 (c) PFL 11.2007
! A test has been added for periodic system, ie the linked test
! is done between i in the central cell, and j in all surrounding cells.
!
!================================================================

SUBROUTINE CalcCnct(na,atome,x,y,z,LIAISONS,r_cov,fact)

  use Path_module, only : NMaxL, max_Z,Nom,Prog, KINT, KREAL

  character(2) :: ATOM
  integer(KINT) :: na,atome(Na),at,ind_zmat(Na,5),long
  real(KREAL) ::  x(Na),y(Na),z(Na),fact
  real(KREAL) :: vx,vy,vz,dist
  REAL(KREAL) :: r_cov(0:Max_Z)
  INTEGER(KINT) :: LIAISONS(Na,0:NMaxL),Nbli,Nblj

  ! For periodic systems
  INTEGER(KINT) :: Ka,Kb,Kc
  LOGICAL   :: Bound

  ! Internals
  logical(KREAL) ::debug

  INTERFACE
     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID
  END INTERFACE

  debug=valid('calccnct')
  if (debug) WRITE(*,*) "======================= Entering CalcCnct ======================="

  DO i=1,na
     Liaisons(i,0)=0
  END DO

  if (debug) THEN
     WRITE(*,*) 'CalcCnct'
     DO iat=1,na
        i=atome(iat)
        WRITE(*,*) Nom(I),I,r_cov(i),r_cov(i)*fact
     END DO
  END IF

  IF (PROG/="VASP") THEN
     DO i=1,na
        NbLi=LIAISONS(i,0)
        DO j=i+1,na
           CALL vecteur(j,i,x,y,z,vx,vy,vz,dist)
           dist=dist/fact
           distth=(r_cov(atome(i))+r_cov(atome(j)))/100.
!           if (debug) WRITE(*,*) atome(i),atome(j),dist,distth
           if (dist.le.distth) THEN
              NbLi=NbLi+1
              LIAISONS(i,NbLi)=j;
              NBlj=LIAISONS(j,0)+1
              LIAISONS(j,Nblj)=i;
              LIAISONS(j,0)=Nblj
           END IF
        END DO
        LIAISONS(i,0)=Nbli
     END DO
  ELSE
     DO i=1,na
        NbLi=LIAISONS(i,0)
        DO j=i+1,na
           Bound=.FALSE.
           DO ka=-1,1
              DO Kb=-1,1
                 DO Kc=-1,1
                    CALL VectorPer(j,i,ka,kb,kc,x,y,z,vx,vy,vz,dist)
                    dist=dist/fact
                    distth=(r_cov(atome(i))+r_cov(atome(j)))/100.
!                    if (debug) WRITE(*,*) atome(i),atome(j),dist,distth
                    if (dist.le.distth) Bound=.TRUE.
                 END DO
              END DO
           END DO
           IF (Bound) THEN
              if (debug) WRITE(*,*) "Adding a bond between:",i,j
              NbLi=NbLi+1
              LIAISONS(i,NbLi)=j;
              NBlj=LIAISONS(j,0)+1
              LIAISONS(j,Nblj)=i;
              LIAISONS(j,0)=Nblj
           END IF
        END DO
        LIAISONS(i,0)=Nbli
     END DO
  END IF
  if (debug) WRITE(*,*) "======================= CalcCnct  over ======================="
END SUBROUTINE CalcCnct
