Chimie Théorique » Opt'n Path

!================================================================

! 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)

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

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

!  Centre National de la Recherche Scientifique,

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

!

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

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

!

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

!  Contact: optnpath@gmail.com

!

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

!

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

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

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

!  or (at your option) any later version.

!

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

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

!

!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

!  GNU Affero General Public License for more details.

!

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

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

!

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

! for commercial licensing opportunities.

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

  use Path_module, only : NMaxL, max_Z,Nom, KINT, KREAL, FPBC, &

       kaBeg,kaEnd,kbBeg,kbEnd,kcBeg,kcEnd

  IMPLICIT NONE

  integer(KINT) :: na, atome(Na)

  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

  REAL(KREAL) :: DistTh

  INTEGER(KINT) :: I,j,iat

  INTERFACE

     function valid(string) result (isValid)

       CHARACTER(*), intent(in) :: string

       logical                  :: isValid

     END function VALID

  END INTERFACE

  debug=valid('calccnct')

  if (debug)  Call Header (" Entering CalcCnct ")

  DO i=1,na

     Liaisons(i,0)=0

  END DO

  if (debug) THEN

     WRITE(*,*) 'CalcCnct : covalent radii used'

     DO iat=1,na

        i=atome(iat)

        WRITE(*,*) Nom(I),I,r_cov(i)*fact

     END DO

     WRITE(*,*) 'Coordinates'

     DO iat=1,na

        i=atome(iat)

        WRITE(*,*) Nom(I),x(iat),y(iat),z(iat)

     END DO

  END IF

  IF (FPBC) THEN

     DO i=1,na

        NbLi=LIAISONS(i,0)

        DO j=i+1,na

           Bound=.FALSE.

           DO ka=kaBeg,kaEnd

              DO Kb=kbBeg,kbEnd

                 DO Kc=kcBeg,kcEnd

                    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

  ELSE

     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

  END IF

  if (debug) Call Header (" CalcCnct over ")

END SUBROUTINE CalcCnct