forked from lijiext/lammps
206 lines
11 KiB
Fortran
206 lines
11 KiB
Fortran
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 !*********************************************************************************
|
|
|