module TPMLib !*************************************************************************************
!
! Common constants, types, and functions for TPM force field
!
!---------------------------------------------------------------------------------------------------
!
! Intel Fortran
!
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, 2020, Version 13.00
!
!***************************************************************************************************

implicit none

!---------------------------------------------------------------------------------------------------
! Mathematical constants
!---------------------------------------------------------------------------------------------------

        real*8, parameter  :: M_PI_2            = 1.57079632679489661923
        real*8, parameter  :: M_PI              = 3.14159265358979323846
        real*8, parameter  :: M_3PI_2           = 4.71238898038468985769
        real*8, parameter  :: M_2PI             = 6.28318530717958647692
        real*8, parameter  :: M_PI_180          = 0.017453292519943295769

!---------------------------------------------------------------------------------------------------
! Physical unit constants
!---------------------------------------------------------------------------------------------------

        real*8, parameter  :: K_AMU             = 1.66056E-27                           ! a.m.u. (atomic mass unit, Dalton)
        real*8, parameter  :: K_EV              = 1.60217646e-19                        ! eV (electron-volt)

        real*8, parameter  :: K_MDLU            = 1.0E-10                               ! MD length unit (m)
        real*8, parameter  :: K_MDEU            = K_EV                                  ! MD energy unit (J)
        real*8, parameter  :: K_MDMU            = K_AMU                                 ! MD mass unit (kg)
        real*8, parameter  :: K_MDFU            = K_MDEU / K_MDLU                       ! MD force unit (N)
        real*8, parameter  :: K_MDCU            = K_MDEU / K_MDMU                       ! MD specific heat unit (J/(kg*K))

!---------------------------------------------------------------------------------------------------
! Global variables
!---------------------------------------------------------------------------------------------------

        integer*4          :: StdUID            = 31

contains !******************************************************************************************

!---------------------------------------------------------------------------------------------------
! Simple mathematical functions
!---------------------------------------------------------------------------------------------------
        
        real*8 function rad ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        real*8, intent(in)      :: X
        !-------------------------------------------------------------------------------------------
                rad = X * M_PI_180
        end function rad !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        real*8 function sqr ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        real*8, intent(in)      :: X
        !-------------------------------------------------------------------------------------------
                sqr = X * X
        end function sqr !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        integer*4 function signum ( X )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        real*8, intent(in)      :: X
        !-------------------------------------------------------------------------------------------   
                if ( X > 0 ) then
                        signum = 1
                else if ( X < 0 ) then
                        signum = -1
                else
                        signum = 0
                end if
        end function signum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!---------------------------------------------------------------------------------------------------
! Vector & matrix functions
!---------------------------------------------------------------------------------------------------

        real*8 function S_V3xx ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        real*8, dimension(0:2), intent(in)      :: V
        !-------------------------------------------------------------------------------------------
                S_V3xx = V(0) * V(0) + V(1) * V(1) + V(2) * V(2)
        end function S_V3xx !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        real*8 function S_V3xV3 ( V1, V2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        real*8, dimension(0:2), intent(in)      :: V1, V2
        !-------------------------------------------------------------------------------------------
                S_V3xV3 = V1(0) * V2(0) + V1(1) * V2(1) + V1(2) * V2(2)
        end function S_V3xV3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        real*8 function S_V3norm3 ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        real*8, dimension(0:2), intent(in)      :: V
        !-------------------------------------------------------------------------------------------
                S_V3norm3 = dsqrt ( V(0) * V(0) + V(1) * V(1) + V(2) * V(2) )
        end function S_V3norm3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        subroutine V3_ort ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Vector production
        !-------------------------------------------------------------------------------------------
        real*8, dimension(0:2), intent(inout)   :: V
        !-------------------------------------------------------------------------------------------
        real*8                                  :: Vabs 
        !-------------------------------------------------------------------------------------------
                Vabs = S_V3norm3 ( V )
                V(0) = V(0) / Vabs
                V(1) = V(1) / Vabs
                V(2) = V(2) / Vabs
        end subroutine V3_ort !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        subroutine V3_V3xxV3 ( V, V1, V2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Vector production
        !-------------------------------------------------------------------------------------------
        real*8, dimension(0:2), intent(out)     :: V
        real*8, dimension(0:2), intent(in)      :: V1, V2
        !-------------------------------------------------------------------------------------------
                V(0) = V1(1) * V2(2) - V1(2) * V2(1)
                V(1) = V1(2) * V2(0) - V1(0) * V2(2)
                V(2) = V1(0) * V2(1) - V1(1) * V2(0)
        end subroutine V3_V3xxV3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!---------------------------------------------------------------------------------------------------
! Handling of spherical and Euler angles
!---------------------------------------------------------------------------------------------------
        
        subroutine RotationMatrix3  ( M, Psi, Tet, Phi ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Ksi, Tet and Phi are Euler angles
        !-------------------------------------------------------------------------------------------
        real*8, dimension(0:2,0:2), intent(out) :: M
        real*8, intent(in)                      :: Psi, Tet, Phi
        !-------------------------------------------------------------------------------------------
        real*8                                  :: cK, sK, cT, sT, cP, sP
        !-------------------------------------------------------------------------------------------
                cK = dcos ( Psi )
                sK = dsin ( Psi )
                cT = dcos ( Tet )
                sT = dsin ( Tet )
                cP = dcos ( Phi )
                sP = dsin ( Phi )
                M(0,0) = cP * cK - sK * sP * cT
                M(0,1) = cP * sK + sP * cT * cK
                M(0,2) = sP * sT
                M(1,0) = - sP * cK - cP * cT * sK
                M(1,1) = - sP * sK + cP * cT * cK
                M(1,2) = cP * sT
                M(2,0) = sT * sK
                M(2,1) = - sT * cK
                M(2,2) = cT
        end subroutine RotationMatrix3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        
        subroutine EulerAngles ( Psi, Tet, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        real*8, intent(out)                     :: Tet, Psi
        real*8, dimension(0:2), intent(in)      :: L
        !-------------------------------------------------------------------------------------------
                Tet = acos ( L(2) )
                Psi = atan2 ( L(1), L(0) )
                if ( Psi > M_3PI_2 ) then
                        Psi = Psi - M_3PI_2
                else 
                        Psi = Psi + M_PI_2
                end if
        end subroutine EulerAngles !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!---------------------------------------------------------------------------------------------------
! File inout and output
!---------------------------------------------------------------------------------------------------

        integer*4 function OpenFile ( Name, Params, Path ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        character*(*), intent(in)       :: Name, Params, Path
        !-------------------------------------------------------------------------------------------
        integer*4                       :: Fuid
        character*512                   :: FullName, Msg, Name1, Action1, Status1, Form1, Position1
        !-------------------------------------------------------------------------------------------
                OpenFile = StdUID + 1
                if ( Params(1:1) == 'r' ) then
                        Action1   = 'read'
                        Status1   = 'old'
                        Position1 = 'rewind'
                else if ( Params(1:1) == 'w' ) then
                        Action1 = 'write'
                        Status1 = 'replace'
                        Position1 = 'rewind'
                else if ( Params(1:1) == 'a' ) then
                        Action1 = 'write'
                        Status1 = 'old'
                        Position1 = 'append'
                endif
                if ( Params(2:2) == 'b' ) then
                        Form1 = 'binary'
                else
                        Form1 = 'formatted'
                endif
                open ( unit = OpenFile, file = Name, form = Form1, action = Action1, status = Status1, position = Position1 )
                StdUID = StdUID + 1
                return
        end function OpenFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        subroutine CloseFile ( Fuid ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        integer*4, intent(inout)        :: Fuid
        !-------------------------------------------------------------------------------------------
                if ( Fuid < 0 ) return
                close ( unit = Fuid )
                Fuid = -1
        end subroutine CloseFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        
end module TPMLib !*********************************************************************************

