Chimie Théorique » Opt'n Path

subroutine egrad_siesta(e,geomcart,gradcart)

  ! This routines calculates the energy and the gradient of 

  ! a molecule, using Siesta 

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

!  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 : Nat, renum,Order,OrderInv, Coord, ProgExe,v_direct,latr

  use Io_module

  IMPLICIT NONE

  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

!--------------------------------------------------------------- Input Variables

      character(len=*), intent(in)           :: routine, msg

      character(len=*), intent(in), optional :: file

      integer(KINT), intent(in), optional      :: line, unit

    END SUBROUTINE die

    SUBROUTINE Warning(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 Warning

 SUBROUTINE WriteList(Input,Unit)

! This routine reads an input template for Siesta

  use VarTypes

  use Io_module

  IMPLICIT NONE

! Input

      TYPE (Input_line), POINTER, INTENT(IN) :: Input

      INTEGER(KINT), OPTIONAL, INTENT(IN) :: Unit

    END SUBROUTINE WriteList

  END INTERFACE

  ! Energy 

  REAL(KREAL), INTENT (OUT) :: e

  ! Geometry in cartesian coordinates

  REAL(KREAL), INTENT (IN) :: geomcart(Nat,3)

  ! Gradient calculated at Geom geometry

  REAL(KREAL), INTENT (OUT) :: gradcart(3*NAt)

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

  character(LCHARS) :: LINE

  logical           :: debug, TChk

  LOGICAL, SAVE :: first=.true.

  REAL(KREAL), ALLOCATABLE :: X(:),Y(:),Z(:) ! Nat

  REAL(KREAL), ALLOCATABLE :: GeomCartLoc(:,:) !Nat,3

  INTEGER(KINT) :: iat, i, J,Idx

  CHARACTER(LCHARS) :: FileIn, FileOut, FileForces,FileE

  CHARACTER(3) :: SepIn=' < ', SepOut=' > '

  CHARACTER(VLCHARS), SAVE :: RunCommand

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

  debug=valid('EGRAD').OR.valid('egrad_siesta')

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

  if (first) THEN

     First=.FALSE.

     call system("uname > Path.tmp")

     OPEN(IOTMP,FILE='Path.tmp')

     READ(IOTMP,'(A)') Line

     CLOSE(IOTMP)

     IF (index(Line,'CYGWIN').NE.0) THEN

        SepIn='   '

        SepOut='   '

     END IF

  END IF

  ALLOCATE(GeomCartLoc(Nat,3))

  gradcart=0.

  E=0.

  FileIn=Trim(CalcName) // "." // Trim(ISuffix)

  FileOut=Trim(CalcName) // "." // Trim(OSuffix)

  FileForces=Trim(Siesta_Label) // ".FA"

  FileE='BASIS_ENTHALPY'

  RunCommand=Trim(ProgExe) // SepIn // Trim(FileIn) // SepOut // Trim(FileOut)

  IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(RunCommand)

  ! we create the input file

  OPEN(IOTMP,File=FileIn)

  Call WriteList(Siesta_input,Unit=IOTMP)

  WRITE(IOTMP,*) 

  IF (V_Direct=='DIRECT') THEN

     ALLOCATE(X(Nat),Y(Nat),Z(Nat))

     X=0.

     Y=0.

     Z=0.

     DO I=1,Nat

        DO J=1,3

           X(I)=X(I)+GeomCart(I,J)*Latr(1,J)

           Y(I)=Y(I)+GeomCart(I,J)*Latr(2,J)

           Z(I)=Z(I)+GeomCart(I,J)*Latr(3,J)

        END DO

     END DO

     GeomCartLoc(1:Nat,1)=X(1:Nat)

     GeomCartLoc(1:Nat,2)=Y(1:Nat)

     GeomCartLoc(1:Nat,3)=Z(1:Nat)

     DEALLOCATE(X,Y,Z)

  ELSE

     GeomCartLoc=GeomCart

  END IF

  WRITE(IOTMP,'(1X,A)')  '%block AtomicCoordinatesAndAtomicSpecies' 

  DO I=1,Nat

     If (renum) THEN

        Iat=Order(I)

        WRITE(IOTMP,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(Iat,:)/Siesta_Unit_Read, IdxSpecies(Iat),TRIM(Siesta_Paste(I))

     ELSE

        Iat=OrderInv(I)

        WRITE(IOTMP,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(I,:)/Siesta_Unit_Read, IdxSpecies(Iat), TRIM(Siesta_Paste(Iat))

     END IF

  END DO

  WRITE(IOTMP,'(1X,A)')  '%endblock AtomicCoordinatesAndAtomicSpecies'

  WRITE(IOTMP,*) 

  CLOSE(IOTMP)

  IF (debug) THEN

     WRITE(*,*) 'DBG EGRAD ' //  Trim(FileIn),' COORD=',TRIM(COORD)

     call system('cat ' // Trim(FileIn))

  END IF

  call system(RunCommand)

  if (debug) WRITE(*,*) 'DBG EGRAD, back from calculation'

  ! Whatever the coordinate system, we read the cartesian gradient

  ! and we convert it to Zmat or other if needed

! in Siesta, the forces are in the file 'Label'.FA so this is easy to read :)

  OPEN(IOTMP,FILE=FileForces, STATUS='old')

! firt line is the atom number

  READ(IOTMP,'(A)') LINE

  if (debug) WRITE(*,*) "GradCart reading - renum=",Renum

  DO I=1,Nat

     Iat=I

     IF (renum) Iat=Order(I)

     READ(IOTMP,'(6X,3(F12.6))')  GradCart(3*Iat-2:3*Iat)

     if (debug) WRITE(*,*) "I,Iat,GradCart:",I,Iat,GradCart(3*Iat-2:3*Iat)

  END DO

! Siesta gives the Forces in eV/Ang, so we convert it into the 

! gradient in au/Angstrom 

  gradCart=-1.0d0*ev2au*gradCart

  CLOSE(IOTMP)

! We now have to read the energy

  INQUIRE(File=FileE,EXIST=TChk)

  If (Tchk) THEN

! We have a BASIS_ENTHALPY file. Great. 

! Energy is the last number on the second line

     OPEN(IOTMP,File=FileE,Status='Old')

     READ(IOTMP,*) Line

     READ(IOTMP,'(A)') Line

     if (debug) WRITE(*,*) "Line E:",TRIM(LINE)

     Idx=INDEX(TRIM(LINE),' ',BACK=.TRUE.)

     Line=Line(Idx+1:)

     if (debug) WRITE(*,*) "Line E:",TRIM(LINE)

     Read(Line,*) E

  ELSE

! We do not have this file... we have to look for it in the .out file

     FileE=FileOut

     OPEN(IOTMP,File=FileE,Status='Old')

! We search for  "siesta:         Total ="

     Line=""

     DO WHILE (INDEX(Line,"siesta:         Total =")/=1)

        READ(IOTMP,*) Line

     END DO

     Idx=INDEX(TRIM(LINE),' ',BACK=.TRUE.)

     Line=Line(Idx+1:)

     Read(Line,*) E

  END IF

  CLOSE (IOTMP)

  WRITE(*,*) 'E=',E

  WRITE(*,*) 'Gradcart='

  DO I=1,Nat

     WRITE(*,'(1X,3(1X,F15.8))') GradCart(3*I-2:3*I)

  END DO

  DeAllocate(GeomCartLoc)

!  Call Die('Egrad_Siesta','Ah ah.',UNIT=IOOUT)

  if (debug) Call header("Egrad_siesta Over")

  RETURN

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

end subroutine egrad_siesta