1! ------------ ----------------------------------------------------------
2!   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3!   https://www.lammps.org/ Sandia National Laboratories
4!   Steve Plimpton, sjplimp@sandia.gov
5!
6!   Copyright (2003) Sandia Corporation.  Under the terms of Contract
7!   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8!   certain rights in this software.  This software is distributed under
9!   the GNU General Public License.
10!
11!   See the README file in the top-level LAMMPS directory.
12!
13!   Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
14!-------------------------------------------------------------------------
15
16module TPMLib !*************************************************************************************
17!
18! Basic constants, types, and mathematical functions.
19!
20!---------------------------------------------------------------------------------------------------
21!
22! Intel Fortran
23!
24! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
25!
26!***************************************************************************************************
27use iso_c_binding, only : c_int, c_double, c_char
28implicit none
29
30!---------------------------------------------------------------------------------------------------
31! Mathematical constants
32!---------------------------------------------------------------------------------------------------
33
34        real(c_double), parameter  :: M_PI_2            = 1.57079632679489661923
35        real(c_double), parameter  :: M_PI              = 3.14159265358979323846
36        real(c_double), parameter  :: M_3PI_2           = 4.71238898038468985769
37        real(c_double), parameter  :: M_2PI             = 6.28318530717958647692
38        real(c_double), parameter  :: M_PI_180          = 0.017453292519943295769
39
40!---------------------------------------------------------------------------------------------------
41! Physical unit constants
42!---------------------------------------------------------------------------------------------------
43
44        real(c_double), parameter  :: K_AMU             = 1.66056E-27           ! a.m.u. (atomic mass unit, Dalton)
45        real(c_double), parameter  :: K_EV              = 1.60217646e-19        ! eV (electron-volt)
46
47        real(c_double), parameter  :: K_MDLU            = 1.0E-10               ! MD length unit (m)
48        real(c_double), parameter  :: K_MDEU            = K_EV                  ! MD energy unit (J)
49        real(c_double), parameter  :: K_MDMU            = K_AMU                 ! MD mass unit (kg)
50        real(c_double), parameter  :: K_MDFU            = K_MDEU / K_MDLU       ! MD force unit (N)
51        real(c_double), parameter  :: K_MDCU            = K_MDEU / K_MDMU       ! MD specific heat unit (J/(kg*K))
52
53!---------------------------------------------------------------------------------------------------
54! Global variables
55!---------------------------------------------------------------------------------------------------
56
57        integer(c_int)             :: StdUID            = 31
58
59contains !******************************************************************************************
60
61!---------------------------------------------------------------------------------------------------
62! Simple mathematical functions
63!---------------------------------------------------------------------------------------------------
64
65        real(c_double) function rad ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
66        real(c_double), intent(in)      :: X
67        !-------------------------------------------------------------------------------------------
68                rad = X * M_PI_180
69        end function rad !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
70
71        real(c_double) function sqr ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
72        real(c_double), intent(in)      :: X
73        !-------------------------------------------------------------------------------------------
74                sqr = X * X
75        end function sqr !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
76
77        integer(c_int) function signum ( X )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
78        real(c_double), intent(in)      :: X
79        !-------------------------------------------------------------------------------------------
80                if ( X > 0 ) then
81                        signum = 1
82                else if ( X < 0 ) then
83                        signum = -1
84                else
85                        signum = 0
86                end if
87        end function signum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
88
89!---------------------------------------------------------------------------------------------------
90! Vector & matrix functions
91!---------------------------------------------------------------------------------------------------
92
93        real(c_double) function S_V3xx ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
94        real(c_double), dimension(0:2), intent(in)      :: V
95        !-------------------------------------------------------------------------------------------
96                S_V3xx = V(0) * V(0) + V(1) * V(1) + V(2) * V(2)
97        end function S_V3xx !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
98
99        real(c_double) function S_V3xV3 ( V1, V2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
100        real(c_double), dimension(0:2), intent(in)      :: V1, V2
101        !-------------------------------------------------------------------------------------------
102                S_V3xV3 = V1(0) * V2(0) + V1(1) * V2(1) + V1(2) * V2(2)
103        end function S_V3xV3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
104
105        real(c_double) function S_V3norm3 ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
106        real(c_double), dimension(0:2), intent(in)      :: V
107        !-------------------------------------------------------------------------------------------
108                S_V3norm3 = dsqrt ( V(0) * V(0) + V(1) * V(1) + V(2) * V(2) )
109        end function S_V3norm3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
110
111        subroutine V3_ort ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
112        real(c_double), dimension(0:2), intent(inout)   :: V
113        !-------------------------------------------------------------------------------------------
114        real(c_double)                                  :: Vabs
115        !-------------------------------------------------------------------------------------------
116                Vabs = S_V3norm3 ( V )
117                V(0) = V(0) / Vabs
118                V(1) = V(1) / Vabs
119                V(2) = V(2) / Vabs
120        end subroutine V3_ort !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
121
122        subroutine V3_V3xxV3 ( V, V1, V2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
123        real(c_double), dimension(0:2), intent(out)     :: V
124        real(c_double), dimension(0:2), intent(in)      :: V1, V2
125        !-------------------------------------------------------------------------------------------
126                V(0) = V1(1) * V2(2) - V1(2) * V2(1)
127                V(1) = V1(2) * V2(0) - V1(0) * V2(2)
128                V(2) = V1(0) * V2(1) - V1(1) * V2(0)
129        end subroutine V3_V3xxV3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
130
131!---------------------------------------------------------------------------------------------------
132! Handling of spherical and Euler angles
133!---------------------------------------------------------------------------------------------------
134
135        subroutine RotationMatrix3  ( M, Psi, Tet, Phi ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
136        ! Ksi, Tet and Phi are Euler angles
137        !-------------------------------------------------------------------------------------------
138        real(c_double), dimension(0:2,0:2), intent(out) :: M
139        real(c_double), intent(in)                      :: Psi, Tet, Phi
140        !-------------------------------------------------------------------------------------------
141        real(c_double)                                  :: cK, sK, cT, sT, cP, sP
142        !-------------------------------------------------------------------------------------------
143                cK = dcos ( Psi )
144                sK = dsin ( Psi )
145                cT = dcos ( Tet )
146                sT = dsin ( Tet )
147                cP = dcos ( Phi )
148                sP = dsin ( Phi )
149                M(0,0) = cP * cK - sK * sP * cT
150                M(0,1) = cP * sK + sP * cT * cK
151                M(0,2) = sP * sT
152                M(1,0) = - sP * cK - cP * cT * sK
153                M(1,1) = - sP * sK + cP * cT * cK
154                M(1,2) = cP * sT
155                M(2,0) = sT * sK
156                M(2,1) = - sT * cK
157                M(2,2) = cT
158        end subroutine RotationMatrix3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
159
160        subroutine EulerAngles ( Psi, Tet, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
161        real(c_double), intent(out)                     :: Tet, Psi
162        real(c_double), dimension(0:2), intent(in)      :: L
163        !-------------------------------------------------------------------------------------------
164                Tet = acos ( L(2) )
165                Psi = atan2 ( L(1), L(0) )
166                if ( Psi > M_3PI_2 ) then
167                        Psi = Psi - M_3PI_2
168                else
169                        Psi = Psi + M_PI_2
170                end if
171        end subroutine EulerAngles !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
172
173!---------------------------------------------------------------------------------------------------
174! File input and output
175!---------------------------------------------------------------------------------------------------
176
177        integer(c_int) function OpenFile ( Name, Params, Path ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
178        character*(*), intent(in)       :: Name, Params, Path
179        !-------------------------------------------------------------------------------------------
180        integer(c_int)                       :: Fuid
181        character*512                   :: FullName, Msg, Name1, Action1, Status1, Form1, Position1
182        !-------------------------------------------------------------------------------------------
183                OpenFile = StdUID + 1
184                if ( Params(1:1) == 'r' ) then
185                        Action1   = 'read'
186                        Status1   = 'old'
187                        Position1 = 'rewind'
188                else if ( Params(1:1) == 'w' ) then
189                        Action1 = 'write'
190                        Status1 = 'replace'
191                        Position1 = 'rewind'
192                else if ( Params(1:1) == 'a' ) then
193                        Action1 = 'write'
194                        Status1 = 'old'
195                        Position1 = 'append'
196                endif
197                if ( Params(2:2) == 'b' ) then
198                        Form1 = 'binary'
199                else
200                        Form1 = 'formatted'
201                endif
202                open ( unit = OpenFile, file = Name, form = Form1, action = Action1, status = Status1, position = Position1 )
203                StdUID = StdUID + 1
204                return
205        end function OpenFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
206
207        subroutine CloseFile ( Fuid ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
208        integer(c_int), intent(inout)        :: Fuid
209        !-------------------------------------------------------------------------------------------
210                if ( Fuid < 0 ) return
211                close ( unit = Fuid )
212                Fuid = -1
213        end subroutine CloseFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
214
215end module TPMLib !*********************************************************************************
216