SUBROUTINE CheckPeriodicBound
  ! This subroutine is only for periodic systems
  ! it checks that the atoms do not change their
  ! positions from one side to the other during the simulations.
  ! This would ruin the interpolation.


!----------------------------------------------------------------------
!  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
  Use Io_module

  IMPLICIT NONE


  INTEGER(KINT) :: IGeom, I,J,K
  REAL(KREAL), ALLOCATABLE :: GeomRef(:,:), GeomTmp(:,:) ! (3,Nat)

  REAL(KREAL) :: Latrloc(3,3),Bloc(3),Xt,Yt,Zt
  REAL(KREAL) :: NormA, NormB,Xrel,Yrel
  REAL(KREAL) :: Prod,LatrA(3),LatrB(3)
  REAL(KREAL) :: V(3)

  LOGICAL :: Debug

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


    SUBROUTINE die(routine, msg, file, line, unit)

      Use VarTypes
      Use io_module

      implicit none
      character(len=*), intent(in)           :: routine, msg
      character(len=*), intent(in), optional :: file
      integer(KINT), intent(in), optional      :: line, unit

    END SUBROUTINE die

  END INTERFACE

  debug=valid('checkperiodicbound')
  if (debug) call header("Entering CheckPeriodicBound")

  ALLOCATE(GeomRef(3,Nat), GeomTmp(3,Nat))

  ! The test is easy in DIRECT space, so we will convert the geometries
  ! into this space. 

  SELECT CASE (Iper)
  CASE(1)
     NormA=DOT_PRODUCT(Lat_a,Lat_a)
     DO IGeom=2,NGeomI
        DO I=1,NAt
           DO J=1,3
              v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I)
           END DO
           XRel=DOT_PRODUCT(Lat_a,v)/NormA
           if (abs(Xrel).GE.BoxTol) THEN
              If (debug) THEN
                 WRITE(*,'("Atom ",I5," moved because XRel=",F10.6)') I,XRel
                 WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)
                 WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)
                 WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Lat_a*Sign(1.0d0,XRel)
              END IF
              XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Lat_a*Sign(1.0d0,XRel)
           END IF
        END DO
     END DO
  CASE (2)
     NormA=DOT_PRODUCT(Lat_a,Lat_a)
     NormB=DOT_PRODUCT(Lat_b,Lat_b)
     Prod=DOT_PRODUCT(Lat_a,Lat_b)
     LatrA=(NormB*Lat_a-Prod*Lat_b)/(NormA*NormB-Prod**2)
     LatrB=(NormA*Lat_b-Prod*Lat_a)/(NormA*NormB-Prod**2)

     DO IGeom=2,NGeomI
        DO I=1,NAt
           DO J=1,3
              v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I)
           END DO
           XRel=DOT_PRODUCT(Latra,v)
           YRel=DOT_PRODUCT(Latrb,v)
           if (abs(Xrel).GE.BoxTol) THEN
              If (debug) THEN
                 WRITE(*,'("Atom ",I5," moved because XRel=",F10.6)') I,XRel
                 WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)
                 WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)
                 WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Latra*Sign(1.0d0,XRel)
              END IF
              XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Latra*Sign(1.0d0,XRel)
           END IF
           if (abs(Yrel).GE.BoxTol) THEN
              If (debug) THEN
                 WRITE(*,'("Atom ",I5," moved because YRel=",F10.6)') I,YRel
                 WRITE(*,'("Ref: ",3(1X,F10.6))') XyzGeomI(1,1:3,I)
                 WRITE(*,'("Old: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)
                 WRITE(*,'("New: ",3(1X,F10.6))') XyzGeomI(IGeom,1:3,I)-Latrb*Sign(1.0d0,YRel)
              END IF
              XyzGeomI(IGeom,1:3,I)=XyzGeomI(IGeom,1:3,I)-Latrb*Sign(1.0d0,YRel)
           END IF
        END DO
     END DO

  CASE(3)



!      if (debug) THEN
!         WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_a:",Lat_a
!         WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_b:",Lat_b
!         WRITE(*,'(1X,A,3(1X,F10.6))') "Lat_c:",Lat_c
!         WRITE(*,'(1X,A,3(1X,F10.6))') "BocTol:",BoxTol

!      NormA=DOT_PRODUCT(Lat_a,Lat_a)
!      NormB=DOT_PRODUCT(Lat_b,Lat_b)
!      Prod=DOT_PRODUCT(Lat_a,Lat_b)
!      LatrA=(NormB*Lat_a-Prod*Lat_b)/(NormA*NormB-Prod**2)
!      LatrB=(NormA*Lat_b-Prod*Lat_a)/(NormA*NormB-Prod**2)
! !     NormA=DOT_PRODUCT(Latra,Latra)
! !     NormB=DOT_PRODUCT(Latrb,Latrb)

!         DO IGeom=2,NGeomI
!            DO I=1,NAt
!               DO J=1,3
!                  v(J)=XyzGeomI(IGeom,J,I)-XyzGeomI(1,J,I)
!               END DO
!               XRel=DOT_PRODUCT(Latra,v)
!               YRel=DOT_PRODUCT(Latrb,v)
!               If (abs(XRel).GT.BoxTol)  WRITE(*,*) "Xrel too big:",IGeom,I,Xrel
!               If (abs(YRel).GT.BoxTol)  WRITE(*,*) "Yrel too big:",IGeom,I,Yrel
! !               WRITE(*,*) "Xrel :",IGeom,I,Xrel
! !               WRITE(*,*) "Yrel :",IGeom,I,Yrel
! !               WRITE(*,*) "Zrel :",IGeom,I,Zrel

!            END DO
!         END DO
!      END IF

     Latrloc(1:3,1)=Lat_a
     Latrloc(1:3,2)=Lat_b
     Latrloc(1:3,3)=Lat_c
     Bloc=1.
     ! TODO: Replace by LAPACK subroutine
     CALL Gaussj(Latrloc,3,3,Bloc,1,1)

     ! We put first geometry as the reference
     GeomRef=0.
     DO I=1,Nat
        DO K=1,3
           DO J=1,3
              GeomRef(K,I)=GeomRef(K,I)+XyzGeomI(1,J,I)*LatrLoc(K,J)
           END DO
        END DO
     END DO

     DO IGeom=2,NGeomI
        GeomTmp=0.
        DO I=1,Nat
           DO K=1,3
              DO J=1,3
                 GeomTmp(K,I)=GeomTmp(K,I)+XyzGeomI(IGeom,J,I)*LatrLoc(K,J)
              END DO
              if (abs(GeomTmp(K,I)-GeomRef(K,I)).GE.BoxTol) THEN
                 If (debug) & 
                      WRITE(*,'("Coord ",I1," has been changed for atom ",I5," by ",F4.0," ref:",F9.6," old:",F9.6," new:",F9.6)') &
                      K,I,Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I)),& 
                      GeomRef(K,I),GeomTmp(K,I), &
                      GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I))

                 GeomTmp(K,I)=GeomTmp(K,I)+Sign(1.d0,GeomRef(K,I)-GeomTmp(K,I))
              END IF
           END DO
           Xt=GeomTmp(1,I)
           Yt=GeomTmp(2,I)
           Zt=GeomTmp(3,I)
           XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1)
           XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2)
           XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3)
        END DO
     END DO

  CASE DEFAULT
     Call Die("CheckPeriodicBound"," Iper unknown")
  END SELECT


  if (debug) WRITE(*,*) "========================= Exiting CheckPeriodicBound ==========="
END SUBROUTINE CheckPeriodicBound


 


   
