 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.


  Use Path_module
  Use Io_module

  IMPLICIT NONE


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

  REAL(KREAL) :: Latrloc(3,3),Bloc(3),Xt,Yt,Zt
  LOGICAL :: Debug

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

  debug=valid('checkperiodicbound')
if (debug) WRITE(*,*) "========================= Entering CheckPeriodicBound ==========="

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

  V_direct=adjustl(Vasp_direct)
  Call upcase(V_direct)

! The test is easy in DIRECT space, so we will convert the geometries
! into this space. 
! We might optimize this as Read_geom converts geometries into NORMAL
! coordinates... so it might be that we are converting them back to direct here...

     Latrloc(1:3,1)=Lat_a
     Latrloc(1:3,2)=Lat_b
     Latrloc(1:3,3)=Lat_c
     Bloc=1.
     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

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


 


   
