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