MODULE lattid
!
!**********************************************************************************
!*  LATTID                                                                        *
!**********************************************************************************
!* This module aims at identifying the lattice type each atom belongs to.         *
!* Lattice "identifiers" (latid) follow the numbering used in OVITO:              *
!*   0 = Other, unknown coordination structure                                    *
!*   1 = FCC, face-centered cubic                                                 *
!*   2 = HCP, hexagonal close-packed                                              *
!*   3 = BCC, body-centered cubic                                                 *
!*   4 = ICO, icosahedral coordination                                            *
!*   5 = SC, simple cubic                                                         *
!*   6 = Cubic diamond                                                            *
!*   7 = Hexagonal diamond                                                        *
!*   8 = Graphene                                                                 *
!* In addition, Atomsk also recognizes the following structures, for which        *
!* additional auxiliary property may also be computed (indicated in [brackets]):  *
!*  10 = rock-salt                                                                *
!*  11 = perovskite [Goldschmidt tolerance factor]                                *
!**********************************************************************************
!* (C) October 2022 - Pierre Hirel                                                *
!*     Université de Lille, Sciences et Technologies                              *
!*     UMR CNRS 8207, UMET - C6, F-59655 Villeneuve D'Ascq, France                *
!*     pierre.hirel@univ-lille.fr                                                 *
!* Last modification: P. Hirel - 23 Jan. 2023                                     *
!**********************************************************************************
!* This program is free software: you can redistribute it and/or modify           *
!* it under the terms of the GNU General Public License as published by           *
!* the Free Software Foundation, either version 3 of the License, or              *
!* (at your option) any later version.                                            *
!*                                                                                *
!* This program 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 General Public License for more details.                                   *
!*                                                                                *
!* You should have received a copy of the GNU General Public License              *
!* along with this program.  If not, see <http://www.gnu.org/licenses/>.          *
!**********************************************************************************
!
USE atoms
USE comv
USE constants
USE messages
USE files
USE neighbors
USE subroutines
USE readin
USE options
USE writeout
!
!
CONTAINS
!
!
SUBROUTINE ID_LATTICE(H,P,AUX,AUXNAMES)
!
!
IMPLICIT NONE
CHARACTER(LEN=*),INTENT(IN):: inputfile
CHARACTER(LEN=*),INTENT(IN):: prefix
!
CHARACTER(LEN=2):: species
CHARACTER(LEN=128):: msg, temp
CHARACTER(LEN=5),DIMENSION(:),ALLOCATABLE:: outfileformats !list of output file formats
CHARACTER(LEN=128),DIMENSION(:),ALLOCATABLE:: AUXNAMES, newAUXNAMES !names of auxiliary properties
CHARACTER(LEN=128),DIMENSION(:),ALLOCATABLE:: comment
CHARACTER(LEN=128),DIMENSION(:),ALLOCATABLE:: options_array !options and their parameters
LOGICAL,DIMENSION(:),ALLOCATABLE:: SELECT  !mask for atom list
INTEGER:: c     !column in AUX that will contain the central symmetry parameters c(i)
INTEGER:: latid !column in AUX that will contain the lattice identifier
INTEGER:: i, ipairs, j, k, ktemp
INTEGER:: Mdefault      !most common number of neighbors
INTEGER,PARAMETER:: Nmax=14 !maximum number of neighbors for any lattice
INTEGER:: Nneigh2           !number of neighbors of an atom
INTEGER,DIMENSION(3):: Nneigh  !number of 1st, 2nd, and 3rd neighbors of an atom
INTEGER,DIMENSION(3):: Neighsp !species of 1st, 2nd, and 3rd neighbors of an atom
INTEGER:: Nspecies      !number of species among an atom's neighbours
INTEGER,DIMENSION(40):: latscore !score of different lattice identifiers
INTEGER,DIMENSION(:,:),ALLOCATABLE:: NeighList !list of index of neighbors
INTEGER,DIMENSION(:),ALLOCATABLE:: newindex  !list of index after sorting
INTEGER,DIMENSION(:,:),ALLOCATABLE:: pairs  !indexes of pairs of atoms
REAL(dp):: distance, dmin
REAL(dp):: sum_dj       !sum of |d_j|²
REAL(dp),DIMENSION(3,3):: Huc   !Box vectors of unit cell (unknown, set to 0 here)
REAL(dp),DIMENSION(3,3):: H       !Base vectors of the supercell
REAL(dp),DIMENSION(3,3):: ORIENT  !crystallographic orientation of the system (mode create)
REAL(dp),DIMENSION(9,9):: C_tensor  !elastic tensor
REAL(dp),DIMENSION(:),ALLOCATABLE:: pairs_distances  !distances
REAL(dp),DIMENSION(:),ALLOCATABLE:: three_angles     !angles
REAL(dp),DIMENSION(:,:),ALLOCATABLE:: aentries    !array containing result
REAL(dp),DIMENSION(:,:),ALLOCATABLE:: AUX, newAUX !auxiliary properties
REAL(dp),DIMENSION(:,:),ALLOCATABLE:: P, S     !positions of atoms, shells
REAL(dp),DIMENSION(:,:),ALLOCATABLE:: PosList
!
!
!Initialize variables
Mdefault = 12
Nspecies = 0
Huc(:,:) = 0.d0
 C_tensor(:,:) = 0.d0
latscore(:) = 0
!
!
msg = 'ENTERING LOCAL_SYM...'
CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
!
CALL ATOMSK_MSG(4069,(/inputfile/),(/0.d0/))
!
!
100 CONTINUE
! Read the input file
CALL READ_AFF(inputfile,H,P,S,comment,AUXNAMES,AUX)
!
! Apply options if any
CALL OPTIONS_AFF(options_array,Huc,H,P,S,AUXNAMES,AUX,ORIENT,SELECT,C_tensor)
IF(nerr>0) GOTO 1000
!
!If auxiliary properties already exist, add a column to save the central symmetry parameter
IF( ALLOCATED(AUX) .AND. SIZE(AUX,2)==SIZE(AUXNAMES) ) THEN
  !If "local_symmetry" and/or "latid" already exist in AUX, use these columns
  c = 0
  latid = 0
  DO i=1,SIZE(AUXNAMES)
    IF( AUXNAMES(i)=="local_symmetry" ) c=i
    IF( AUXNAMES(i)=="latid" ) latid=i
  ENDDO
  !Otherwise, create a new column to store local symmetry parameter
  k = SIZE(AUXNAMES)
  IF( c==0 ) THEN
    k = k+1
    c = k
  ENDIF
  IF( latid==0 ) THEN
    k = k+1
    latid = k
  ENDIF
  !
  IF( k>SIZE(AUXNAMES) ) THEN
    !Resize arrays
    ALLOCATE(newAUXNAMES(k))
    DO i=1,SIZE(AUXNAMES)
      newAUXNAMES(i) = AUXNAMES(i)
    ENDDO
    DEALLOCATE(AUXNAMES)
    ALLOCATE(AUXNAMES(SIZE(newAUXNAMES)))
    AUXNAMES(:) = newAUXNAMES(:)
    DEALLOCATE(newAUXNAMES)
    ALLOCATE(newAUX(SIZE(P,1),k))
    newAUX(:,:) = 0.d0
    DO i=1,SIZE(AUX,1)
      newAUX(i,:) = AUX(i,:)
    ENDDO
    DEALLOCATE(AUX)
    ALLOCATE(AUX(SIZE(newAUX,1),SIZE(newAUX,2)))
    AUX(:,:) = newAUX(:,:)
    DEALLOCATE(newAUX)
  ENDIF
ELSE
  !No auxiliary prop. existed: allocate new arrays
  IF(ALLOCATED(AUX)) DEALLOCATE(AUX)
  IF(ALLOCATED(AUXNAMES)) DEALLOCATE(AUXNAMES)
  c = 1
  latid = 2
  ALLOCATE(AUXNAMES(2))
  ALLOCATE( AUX(SIZE(P,1),2) )
  AUX(:,:) = 0.d0
ENDIF
!
AUXNAMES(c) = "local_symmetry"
AUXNAMES(latid) = "latid"
!
!
!
200 CONTINUE
! Determine how many different atom species are present, and their number
CALL FIND_NSP(P(:,4),aentries)
WRITE(msg,*) "Number of different species:", SIZE(aentries,1)
CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
IF( SIZE(aentries,1)>1 ) THEN
  !Sort them by increasing number of atoms
  CALL BUBBLESORT(aentries(:,:),2,'up  ',newindex)
ENDIF
!
! Determine what type of system it is: unary, binary, ternary, other
IF( SIZE(aentries,1)==1 ) THEN
  !No ambiguity, it is a unary system
  Nspecies = 1
ELSEIF( SIZE(aentries,1)==2 ) THEN
  Nspecies = 2
  !It may be a unary system with just a few solute atoms
  !Check if a species is a minority
  IF( aentries(1,2) > 5.d0*aentries(2,2) .OR. aentries(2,2) > 5.d0*aentries(1,2) ) THEN
    !One of the species is a minority => consider that it is a unary system
    Nspecies = 1
  ELSE
    !Consider that it is a binary system
    Nspecies = 2
  ENDIF
ELSEIF( SIZE(aentries,1)==3 ) THEN
  Nspecies = 3
  !It may be a unary system with just a few solute atoms
  !Check which species is (or are) a minority
  IF( aentries(1,2) > 3.9d0*aentries(2,2) .AND. aentries(1,2) > 3.9d0*aentries(3,2) ) THEN
    !Species #1 is a majority compared to #2 and #3 => consider that it is a unary system
    Nspecies = 1
  ELSEIF( aentries(2,2) > 3.9d0*aentries(1,2) .AND. aentries(2,2) > 3.9d0*aentries(3,2) ) THEN
    !Species #2 is a majority compared to #1 and #3 => consider that it is a unary system
    Nspecies = 1
  ELSEIF( aentries(3,2) > 3.9d0*aentries(1,2) .AND. aentries(3,2) > 3.9d0*aentries(2,2) ) THEN
    !Species #3 is a majority compared to #1 and #2 => consider that it is a binary system
    Nspecies = 1
  ELSEIF( aentries(1,2) > 3.9d0*aentries(3,2) .AND. aentries(2,2) > 3.9d0*aentries(3,2) ) THEN
    !Species #3 is a minority compared to #1 and #2 => consider that it is a binary system
    Nspecies = 2
  ELSEIF( aentries(1,2) > 3.9d0*aentries(2,2) .AND. aentries(3,2) > 3.9d0*aentries(2,2) ) THEN
    !Species #2 is a minority compared to #1 and #3 => consider that it is a binary system
    Nspecies = 2
  ELSEIF( aentries(2,2) > 3.9d0*aentries(1,2) .AND. aentries(3,2) > 3.9d0*aentries(1,2) ) THEN
    !Species #1 is a minority compared to #2 and #3 => consider that it is a binary system
    Nspecies = 2
  ELSE
    !No minority species => it is a ternary material
    Nspecies = 3
  ENDIF
ELSE
  Nspecies = 0
  !PRINT*, "X!X ERROR Nspecies is zero"
  GOTO 1000
ENDIF
WRITE(msg,*) "Nspecies:", Nspecies
CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
!
!
!
300 CONTINUE
! Compute the central symmetry parameter for each atom
!
!Construct neighbor list
CALL ATOMSK_MSG(11,(/""/),(/0.d0/))
CALL VERLET_LIST(H,P,6.d0,NeighList)
CALL ATOMSK_MSG(15,(/""/),(/0.d0/))
!PRINT*, "SIZE NeighList = ", SIZE(NeighList,1), SIZE(NeighList,2)
IF( .NOT.ALLOCATED(NeighList) .OR. SIZE(NeighList)<1 ) THEN
  nerr=nerr+1
  CALL ATOMSK_MSG(1815,(/""/),(/0.d0/))
  GOTO 1000
ENDIF
!
IF( verbosity==4 ) THEN
  !Some debug messages
  IF( ALLOCATED(NeighList) .AND. SIZE(NeighList,1)>0 ) THEN
    WRITE(msg,*) "Sample of neighbor list:"
    CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
    DO i=1,MIN(10,SIZE(NeighList,1))
      WRITE(msg,'(i5,a3,20i5)') i, " | ", NeighList(i,1:MIN(SIZE(NeighList,2),16))
      msg = TRIM(ADJUSTL(msg))//' (...)'
      CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
    ENDDO
    IF( i>=10 ) THEN
      WRITE(msg,*) '      (...discontinued...)'
      CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
    ENDIF
  ENDIF
ENDIF
!
!
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP& PRIVATE(i,j,k,ktemp,ipairs,latscore,Nneigh,Neighsp,sum_dj,PosList,newindex) &
!$OMP& PRIVATE(pairs,pairs_distances,three_angles,dmin,distance,msg,temp)
DO i=1,SIZE(P,1)
  !i is the index of the central atom
  !Initialize variables and arrays for atom #i
  ipairs = 0
  ktemp = 0
  Nspecies = 0
  Nneigh(:) = 0
  Neighsp(:) = 0
  latscore(:) = 0
  sum_dj = 0.d0
  IF(ALLOCATED(pairs)) DEALLOCATE(pairs)
  IF(ALLOCATED(pairs_distances)) DEALLOCATE(pairs_distances)
  IF(ALLOCATED(three_angles)) DEALLOCATE(three_angles)
  !
  !Get the positions of neighbors of atom #i
  CALL NEIGHBOR_POS(H,P,P(i,1:3),NeighList(i,:),ALLOCATED(NeighList),6.d0,PosList)
  !
  IF( ALLOCATED(PosList) .AND. SIZE(PosList,1)>0 ) THEN
    !
    WRITE(msg,*) "Atom # ", i, " : size of PosList = ", SIZE(PosList,1)
    CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
    !Sort neighbors by increasing distance
    CALL BUBBLESORT(PosList(:,:),4,'up  ',newindex)
    !
    !
    !Count number of 1st, 2nd, and 3rd neighbours
    !Compute distance between central atom #i and its third neighbour in the list
    k = MIN(3,SIZE(PosList,1))
    distance = SUM(PosList(1:k,4)) / DBLE(k)
    !Get species of 1st neighbours: it is the most frequent species among 1st neighbours
    k = NINT(PosList(1,5))
    Neighsp(1) = P(k,4)
    !Loop on all neighbours. All neighbours at a distance smaller or comparable to 3rd neighbour,
    !*or* with the same species as the 1st-neighbours, are saved as first neighbours
    DO j=1,MIN(Nmax,SIZE(PosList,1))
      IF( PosList(j,4) < 1.3d0*distance .OR. NINT(P(k,4))==Neighsp(1) ) THEN
        !Atom #j's distance to central atom #i is comparable to that of 3rd neighbour
        !=> count atom #j as a first neighbour of atom #i
        Nneigh(1) = Nneigh(1)+1
        sum_dj = sum_dj + ( VECLENGTH(PosList(j,1:3)-P(i,1:3)) )**2
      ENDIF
    ENDDO
    !Now we know that there are Nneigh(1) first neighbours
    !Count 2nd neighbours
    DO j=Nneigh(1)+1,SIZE(PosList,1)
      !Compute distance between central atom #i and the Nneigh(1)+1 neighbour
      k = MIN(Nneigh(1)+1,SIZE(PosList,1))
      distance = PosList(k,4)
      IF( PosList(j,4) < 1.1d0*distance ) THEN
        !Atom #j is a second neighbour to atom #i
        Nneigh(2) = Nneigh(2)+1
      ENDIF
      !Save species of 2nd neighbour, assume it is the species of all 2nd neighbours
      k = NINT(PosList(Nneigh(1)+1,5))
      Neighsp(2) = P(k,4)
    ENDDO
    !Now we know that there are Nneigh(1) first neighbours and Nneigh(2) second neighbours
    !Count 3rd neighbours
    DO j=Nneigh(1)+Nneigh(2)+1,SIZE(PosList,1)
      !Compute distance between central atom #i and the Nneigh(1)+1 neighbour
      k = MIN(Nneigh(1)+Nneigh(2)+1,SIZE(PosList,1))
      distance = PosList(k,4)
      IF( PosList(j,4) < 1.1d0*distance ) THEN
        !Atom #j is a third neighbour to atom #i
        Nneigh(3) = Nneigh(3)+1
      ENDIF
      !Save species of 3rd neighbour, assume it is the species of all 3rd neighbours
      k = NINT(PosList(Nneigh(1)+Nneigh(2)+1,5))
      Neighsp(3) = P(k,4)
    ENDDO
    !
    !For the following, keep only first neighbors to compute local symmetry parameter
    !PosList(Nneigh(1)+1:,:) = 0.d0
    !
    !Based on number of 1st, 2nd and 3rd neighbours, increase score of matching lattice type(s)
    ! 0=unknown; 1=FCC; 2=HCP; 3=BCC; 4=ICO; 5=CUBIC; 6=c.DIAMOND; 7=h.DIAMOND; 8=GRAPHENE;
    ! 21=ROCK-SALT; 31=PEROVSKITE
    SELECT CASE(Nneigh(1))
    CASE(3)
      !3 first neighbours: lattice may be GRAPHENE
      latscore(8) = latscore(8)+1
      IF( Nneigh(3)==6 ) THEN
        latscore(8) = latscore(8)+1
      ENDIF
    CASE(4)
      !4 first neighbours: lattice may be DIAMOND
      latscore(6) = latscore(6)+1
      IF( Nneigh(2)==12 ) THEN
        latscore(6) = latscore(6)+1
      ENDIF
    CASE(6)
      !6 first neighbours: lattice may be ROCK-SALT or PEROVSKITE
      latscore(21) = latscore(21)+1  !rock-salt
      latscore(31) = latscore(31)+1  !perovskite
      !Check if species of central atom and 1st neighbours are different
      IF( NINT(P(i,4)).NE.Neighsp(1) ) THEN
        latscore(21) = latscore(21)+1  !rock-salt
        latscore(31) = latscore(31)+1  !perovskite
        !Check if species of central atom and 2nd neighbours are the same
        IF( NINT(P(i,4))==Neighsp(2) ) THEN
          latscore(21) = latscore(21)+1  !rock-salt
          latscore(31) = latscore(31)-1  !perovskite
        ELSE
          latscore(21) = latscore(21)-1  !rock-salt
          latscore(31) = latscore(31)+1  !perovskite
        ENDIF
        !Check if species of 1st and 3rd neighbours are the same
        IF( Neighsp(1)==Neighsp(3) ) THEN
          latscore(21) = latscore(21)+1  !rock-salt
          latscore(31) = latscore(31)-1  !perovskite
        ELSE
          latscore(21) = latscore(21)-1  !rock-salt
          latscore(31) = latscore(31)+1  !perovskite
        ENDIF
      ENDIF
    CASE(8)
      !8 first neighbours: lattice may be BCC
      latscore(3) = latscore(3)+1
    CASE(12)
      !12 first neighbours: lattice may be FCC or HCP
      latscore(1) = latscore(1)+1    !fcc
      latscore(2) = latscore(2)+1    !hcp
      !Compare 3rd neighbours: should be 24 for FCC, 2 for HCP
      IF( Nneigh(3)>4 ) THEN
        latscore(1) = latscore(1)+1    !fcc
        latscore(2) = latscore(2)-1    !hcp
      ELSE
        latscore(1) = latscore(1)-1  !fcc
        latscore(2) = latscore(2)+1  !hcp
      ENDIF
    END SELECT
    !
    !Make the PosList more compact by removing intermediate zeros
    DO j=2,SIZE(PosList,1)
      IF( NINT(PosList(j,5))==0 ) THEN
        DO k=j,SIZE(PosList,1)-1
          PosList(k,:) = PosList(k+1,:)
        ENDDO
      ENDIF
    ENDDO
    !
    IF( verbosity==4 ) THEN
      !Some debug messages
      CALL ATOMSPECIES(P(i,4),species)
      WRITE(msg,*) i
      WRITE(temp,*) Nneigh(1)
      WRITE(msg,*) "Number of 1st neighbors of atom # "//TRIM(ADJUSTL(msg))//&
                & "("//TRIM(species)//"): "//TRIM(ADJUSTL(temp))
      CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
      WRITE(msg,*) "   j |      x        y        z    |     d_ij    at.number"
      CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
      DO j=1,MIN(20,Nneigh(1))
        WRITE(msg,'(i5,a3,3f9.3,a3,2f9.3)') j, " | ", PosList(j,1:3), " | ", PosList(j,4), P(NINT(PosList(j,5)),4)
        CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
      ENDDO
      IF( j>=20 ) THEN
        WRITE(msg,*) '      (...discontinued...)'
        CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
      ENDIF
    ENDIF
    !
    IF( Nneigh(1) >= 2 ) THEN
      !Atom #i has many neighbors, we can work with that
      Nneigh2 = Nneigh(1)
      !
      !We must make pairs of atoms => Nneigh has to be an even number
      IF( MOD(Nneigh2,2) .NE. 0 ) THEN
        Nneigh2 = Nneigh2-1
      ENDIF
      !
      !We will need a table containing pairs of atoms
      ALLOCATE( pairs(Nneigh2/2,2) )
      pairs(:,:) = 0
      ALLOCATE( pairs_distances(Nneigh2/2) )
      pairs_distances(:) = 0.d0
      !
      !Find atoms j and k that are opposite
      ipairs=0
      DO j=1,MIN(Nneigh2,Mdefault)
        IF( .NOT. ANY(pairs(:,:)==j) ) THEN
          !We have the first member of a new pair
          ipairs = ipairs+1
          pairs(ipairs,1) = j
          !
          !Among the remaining atoms, find the atom #k that minimizes D = |dj + dk|²
          ! *AND* that is of the same species as atom #j
          dmin = 1.d12
          ktemp = 0
          DO k=j+1,MIN(Nneigh2,Mdefault)
            IF( NINT(PosList(k,4)) == NINT(PosList(j,4)) ) THEN
              !Atom #k is of same species as atom #j
              IF( .NOT. ANY(pairs(:,:)==k) ) THEN
                !Atom #k was not paired with another atom yet
                distance = ( VECLENGTH( PosList(j,1:3) + PosList(k,1:3) - 2.d0*P(i,1:3) ) )**2
                IF( distance < dmin ) THEN
                  dmin = distance
                  ktemp = k
                ENDIF
              ENDIF
            ELSE
              !Species of atom #k is different from atom #j
            ENDIF
          ENDDO
          !Now atom #ktemp is the best match to atom #j
          IF( ktemp>0 ) THEN
            pairs(ipairs,2) = ktemp
            pairs_distances(ipairs) = dmin
          ELSE
            !No suitable atom was found to pair atom #j => remove this entry
            pairs(ipairs,:) = 0
            ipairs = ipairs-1
          ENDIF
        ENDIF
      ENDDO
      !
      IF( verbosity==4 ) THEN
        WRITE(msg,*) i
        WRITE(msg,*) "Pairs of atoms around atom # "//TRIM(ADJUSTL(msg))//" :"
        CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
        WRITE(msg,*) "   j    k      |dj+dk|² "
        CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
        DO j=1,SIZE(pairs,1)
          WRITE(msg,'(2i5,f12.3)') pairs(j,1), pairs(j,2), pairs_distances(j)
          CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
        ENDDO
      ENDIF
      !
      !Compute the central symmetry parameter c of atom #i
      AUX(i,c) = SUM(pairs_distances(:)) / (2.d0*sum_dj)
      !
      DEALLOCATE(pairs)
      DEALLOCATE(pairs_distances)
      !
      !Update scores of lattice types
      
      !
      !
      !Check if atom #i has 4 neighbours: maybe it has a tetrahedral environment
      !e.g. in silicon or diamond lattice, or in minerals containing SiO4 tetrahedra
      !Actually do this calculation if Nneigh is 3, 4 or 5
      IF( Nneigh(1)>=3 .AND. Nneigh(1)<=5 ) THEN
        !Compute angles between pairs of neighbours, using atom #i as central atom
        !In a perfect tetrahedron we expect 6 unique pairs, but here Nneigh may be 3, 4 or 5
        !therefore there are  Nneigh*(Nneigh-1)/2  unique pairs
        ipairs = Nneigh(1)*(Nneigh(1)-1)/2
        ALLOCATE( three_angles(ipairs) )
        three_angles(:) = 0.d0
        ipairs=0
        DO j=1,Nneigh(1)-1
          DO k=j+1,Nneigh(1)
            ipairs=ipairs+1
            !Compute angle between position vectors r_ij and r_ik
            dmin = ANGVEC( PosList(j,1:3)-P(i,1:3) , PosList(k,1:3)-P(i,1:3) )
            !Compute the squared difference of cosine, save it into three_angles(:)
            !NOTE: in a perfect tetrahedron the angle is a=acos(-1/3)=109.471°, hence -cos(a)=1/3
            three_angles(ipairs) = ( DCOS(dmin) + 1.d0/3.d0 )**2
          ENDDO
        ENDDO
        !Compute the sum of all angles and normalize
        dmin = SUM(three_angles(:)) / SIZE(three_angles)
        WRITE(msg,*) dmin
        WRITE(msg,*) "Normalized sum of squared dihedral angles: "//TRIM(ADJUSTL(msg))
        CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
        !If this criterion is smaller than centro-symmetry criterion, use it instead
        IF( dmin<AUX(i,c) ) THEN
          AUX(i,c) = dmin
          !Increase score of diamond lattice
          latscore(6) = latscore(6)+1
        ENDIF
        DEALLOCATE(three_angles)
      ENDIF
      !
      !
      !Check if atom #i has 3 neighbours: maybe it is a sp2 ("graphene") environment
      !Actually do this calculation if Nneigh is 3 or 4
      IF( Nneigh(1)>=3 .AND. Nneigh(1)<=4 ) THEN
        !Compute angles between pairs of neighbours, using atom #i as central atom
        !In perfect sp2 we expect 3 unique pairs, but here Nneigh may be 3 or 4
        !therefore there are  Nneigh*(Nneigh-1)/2  unique pairs
        ipairs = Nneigh(1)*(Nneigh(1)-1)/2
        ALLOCATE( three_angles(ipairs) )
        three_angles(:) = 0.d0
        ipairs=0
        DO j=1,Nneigh(1)-1
          DO k=j+1,Nneigh(1)
            ipairs=ipairs+1
            !Compute angle between position vectors r_ij and r_ik
            dmin = ANGVEC( PosList(j,1:3)-P(i,1:3) , PosList(k,1:3)-P(i,1:3) )
            !Compute the squared difference of cosine, save it into three_angles(:)
            !NOTE: in a perfect sp2 the angle is 120°, hence cos(120°)=-1/2
            three_angles(ipairs) = ( DCOS(dmin) + 0.5d0 )**2
          ENDDO
        ENDDO
        !Compute the sum of all angles and normalize
        dmin = SUM(three_angles(:)) / SIZE(three_angles)
        WRITE(msg,*) dmin
        WRITE(msg,*) "Normalized sum of squared dihedral angles: "//TRIM(ADJUSTL(msg))
        CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
        !If this criterion is smaller than centro-symmetry criterion, use it instead
        IF( dmin<AUX(i,c) ) THEN
          AUX(i,c) = dmin
          !Increase score of graphene lattice
          latscore(8) = latscore(8)+1
        ENDIF
        DEALLOCATE(three_angles)
      ENDIF
      !
      !
      !Save lattice identifier
      AUX(i,latid) = MAXLOC(latscore,DIM=1)
      !
      !
    ELSEIF( Nneigh(1)==1 ) THEN
      !Atom #i has only one neighbor
      AUX(i,c) = 1.d0
      AUX(i,latid) = 0.d0
    ELSE
      !Isolated atom with no neighbor: set parameter to zero
      AUX(i,c) = 0.d0
      AUX(i,latid) = 0.d0
    ENDIF
    !
    WRITE(msg,*) i
    WRITE(msg,*) "c("//TRIM(ADJUSTL(msg))//") = ", AUX(i,c)
    CALL ATOMSK_MSG(999,(/TRIM(msg)/),(/0.d0/))
    !
    IF( AUX(i,c) > 1.d3 .OR. AUX(i,c) < 0.d0 ) THEN
      !Something went wrong in the calculation!
      WRITE(*,*) "/!\ WARNING: auxiliary property out-of-bound for atom: ", i, AUX(i,c)
      AUX(i,c) = 0.d0
    ENDIF
    !
  ENDIF
  !
ENDDO
!$OMP END PARALLEL DO
!
!
!
400 CONTINUE
!Write atom positions and central symmetry parameters into output file(s)
IF( .NOT.ALLOCATED(outfileformats) ) THEN
  ALLOCATE(outfileformats(1))
  outfileformats(1) = "cfg"
ENDIF
CALL WRITE_AFF(prefix,outfileformats,H,P,S,comment,AUXNAMES,AUX)
!
!
GOTO 1000
!
!
!
850 CONTINUE
nerr = nerr+1
GOTO 1000
!
!
!
1000 CONTINUE
IF( ALLOCATED(pairs) ) DEALLOCATE(pairs)
IF( ALLOCATED(pairs_distances) ) DEALLOCATE(pairs_distances)
IF( ALLOCATED(three_angles) ) DEALLOCATE(three_angles)
!
!
END SUBROUTINE ID_LATTICE
!
!
END MODULE lattid
