Chimie Théorique » Opt'n Path

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