You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
275 lines
8.6 KiB
275 lines
8.6 KiB
module functions
|
|
|
|
use parameters
|
|
|
|
implicit none
|
|
|
|
public
|
|
contains
|
|
|
|
|
|
! Functions below this comment are taken from the functions module of atomsk
|
|
!********************************************************
|
|
! STRUPCASE
|
|
! This function reads a string of any length
|
|
! and capitalizes all letters.
|
|
!********************************************************
|
|
FUNCTION StrUpCase (input_string) RESULT (UC_string)
|
|
!
|
|
IMPLICIT NONE
|
|
CHARACTER(*),PARAMETER:: lower_case = 'abcdefghijklmnopqrstuvwxyz'
|
|
CHARACTER(*),PARAMETER:: upper_case = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
|
|
CHARACTER(*),INTENT(IN):: input_string
|
|
CHARACTER(LEN(Input_String)):: UC_string !Upper-Case string that is produced
|
|
INTEGER:: i, n
|
|
!
|
|
IF(LEN(input_string)==0) RETURN
|
|
UC_string = input_string
|
|
! Loop over string elements
|
|
DO i=1,LEN(UC_string)
|
|
!Find location of letter in lower case constant string
|
|
n = INDEX( lower_case, UC_string(i:i) )
|
|
!If current substring is a lower case letter, make it upper case
|
|
IF(n>0) THEN
|
|
UC_string(i:i) = upper_case(n:n)
|
|
ENDIF
|
|
END DO
|
|
!
|
|
END FUNCTION StrUpCase
|
|
|
|
!********************************************************
|
|
! STRDNCASE
|
|
! This function reads a string of any length
|
|
! and transforms all letters to lower case.
|
|
!********************************************************
|
|
FUNCTION StrDnCase (input_string) RESULT (lc_string)
|
|
!
|
|
IMPLICIT NONE
|
|
CHARACTER(*),PARAMETER:: lower_case = 'abcdefghijklmnopqrstuvwxyz'
|
|
CHARACTER(*),PARAMETER:: upper_case = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
|
|
CHARACTER(*),INTENT(IN):: input_string
|
|
CHARACTER(LEN(Input_String)):: lc_string !Lower-Case string that is produced
|
|
INTEGER:: i, n
|
|
!
|
|
IF(LEN(input_string)==0) RETURN
|
|
lc_string = input_string
|
|
! Loop over string elements
|
|
DO i=1,LEN(lc_string)
|
|
!Find location of letter in upper case constant string
|
|
n = INDEX( upper_case, lc_string(i:i) )
|
|
!If current substring is an upper case letter, make it lower case
|
|
IF(n>0) THEN
|
|
lc_string(i:i) = lower_case(n:n)
|
|
ENDIF
|
|
END DO
|
|
!
|
|
END FUNCTION StrDnCase
|
|
|
|
pure function matrix_normal(a, n)
|
|
|
|
integer :: i
|
|
integer, intent(in) :: n
|
|
real(kind = dp), dimension(n) :: v
|
|
real(kind = dp), dimension(n, n),intent(in) :: a
|
|
real(kind = dp), dimension(n,n) :: matrix_normal
|
|
|
|
matrix_normal(:, :) = a(:, :)
|
|
|
|
do i = 1, n
|
|
|
|
v(:) = a(i,:)
|
|
|
|
matrix_normal(i, :) = v(:) / norm2(v)
|
|
|
|
end do
|
|
|
|
return
|
|
end function matrix_normal
|
|
|
|
pure function cross_product(a, b)
|
|
!Function which finds the cross product of two vectors
|
|
|
|
real(kind = dp), dimension(3), intent(in) :: a, b
|
|
real(kind = dp), dimension(3) :: cross_product
|
|
|
|
cross_product(1) = a(2) * b(3) - a(3) * b(2)
|
|
cross_product(2) = a(3) * b(1) - a(1) * b(3)
|
|
cross_product(3) = a(1) * b(2) - a(2) * b(1)
|
|
|
|
return
|
|
end function cross_product
|
|
|
|
pure function identity_mat(n)
|
|
!Returns the nxn identity matrix
|
|
|
|
integer :: i
|
|
integer, intent(in) :: n
|
|
real(kind = dp), dimension(n, n) :: identity_mat
|
|
|
|
identity_mat(:, :) = 0.0_dp
|
|
do i = 1, n
|
|
identity_mat(i, i) = 1.0_dp
|
|
end do
|
|
|
|
return
|
|
end function identity_mat
|
|
|
|
pure function triple_product(a, b, c)
|
|
!triple product between three 3*1 vectors
|
|
|
|
real(kind = dp) :: triple_product
|
|
real(kind = dp), dimension(3), intent(in) :: a, b, c
|
|
triple_product = dot_product(a, cross_product(b, c))
|
|
|
|
return
|
|
end function triple_product
|
|
|
|
function in_bd_lat(v, box_faces, box_norms)
|
|
!This function returns whether the point is within the transformed box boundaries. The transformed
|
|
!box being the transformed simulation cell in the lattice basis
|
|
|
|
!Input/output variables
|
|
real(kind=dp), dimension(3), intent(in) :: v !integer lattice position
|
|
real(kind=dp), dimension(3,6), intent(in) :: box_faces !Centroid of all the box faces
|
|
real(kind=dp), dimension(3,6), intent(in) :: box_norms !Box face normals
|
|
logical :: in_bd_lat
|
|
|
|
!Other variables
|
|
integer :: i
|
|
real(kind=dp) :: pt_to_face(3)
|
|
|
|
in_bd_lat = .true.
|
|
|
|
!Check if point is in box bounds, this works by comparing the dot product of the face normal and the
|
|
!vector from the point to the face. If the dot product is greater than 0 then the point is behind the face
|
|
!if it is equal to zero then the point is on the face, if is less than 0 the the point is in front of the face.
|
|
do i = 1, 6
|
|
pt_to_face(:) = box_faces(:, i) - v
|
|
if(dot_product(pt_to_face, box_norms(:,i)) <= 0) then
|
|
in_bd_lat = .false.
|
|
exit
|
|
end if
|
|
end do
|
|
|
|
return
|
|
end function in_bd_lat
|
|
|
|
function in_block_bd(v, box_bd)
|
|
!This function determines whether a point is within a block in 3d
|
|
!Input/output
|
|
real(kind=dp), dimension(3), intent(in) :: v
|
|
real(kind=dp), dimension(6), intent(in) :: box_bd
|
|
logical :: in_block_bd
|
|
|
|
!Other variables
|
|
integer :: i
|
|
|
|
in_block_bd = .true.
|
|
|
|
do i =1 ,3
|
|
!Check upper bound
|
|
if(v(i) > (box_bd(2*i)+10.0_dp**(-10)) ) then
|
|
in_block_bd =.false.
|
|
exit
|
|
!Check lower bound
|
|
else if (v(i) < (box_bd(2*i-1)-10.0_dp**(-10))) then
|
|
in_block_bd = .false.
|
|
exit
|
|
end if
|
|
end do
|
|
end function in_block_bd
|
|
|
|
function in_wedge_bd(r,vertex)
|
|
!This code determines whether the 2dimensional projection of a point lies within the 2 dimensional projection of a wedge.
|
|
real(kind=dp), intent(in) :: r(3) !This is the point position
|
|
real(kind=dp), intent(in) :: vertex(3,3) !These are the relevant vertex positions for the wedge
|
|
real(kind=dp) :: v1(3), v2(3), v3(3), c1(3), c2(3) !Vertex vector to point and cross products
|
|
integer :: i
|
|
logical :: in_wedge_bd
|
|
|
|
in_wedge_bd = .true.
|
|
do i = 1, 3
|
|
v1 = vertex(:,mod(i,3)+1) - vertex(:,i)
|
|
v2 = r - vertex(:,i)
|
|
v3 = vertex(:,mod(i+1,3)+1) - vertex(:,i)
|
|
c1 = cross_product(v1,v2)
|
|
c2 = cross_product(v1,v3)
|
|
if(dot_product(c1,c2) < 0) then
|
|
in_wedge_bd=.false.
|
|
exit
|
|
end if
|
|
end do
|
|
end function in_wedge_bd
|
|
|
|
function lcm(a,b)
|
|
!This function returns the smallest least common multiple of two numbers
|
|
|
|
real(kind=dp), intent(in) :: a, b
|
|
real(kind=dp) :: lcm
|
|
|
|
integer :: aint, bint, gcd, remainder, placeholder
|
|
|
|
!Cast the vector positions to ints. There will be some error associated with this calculation
|
|
aint = a*10**2
|
|
bint = b*10**2
|
|
|
|
!Calculate greated common divisor
|
|
gcd = aint
|
|
placeholder = bint
|
|
do while(placeholder /= 0)
|
|
remainder = modulo(gcd, placeholder)
|
|
gcd = placeholder
|
|
placeholder=remainder
|
|
end do
|
|
lcm = real((aint*bint),dp)/(real(gcd,dp))* 10.0_dp**(-2.0_dp)
|
|
end function lcm
|
|
|
|
function is_neighbor(rl, rk, r_in, r_out)
|
|
!This function checks to see if two atoms are within a shell with an inner radius r_in and outer radius
|
|
!r_out
|
|
real(kind=dp), intent(in) :: r_in, r_out
|
|
real(kind=dp), dimension(3), intent(in) :: rl, rk
|
|
logical :: is_neighbor
|
|
|
|
!Internal variable
|
|
real(kind=dp) :: rlk
|
|
|
|
rlk = norm2(rk - rl)
|
|
is_neighbor=.true.
|
|
if((rlk>r_out).or.(rlk < r_in)) is_neighbor = .false.
|
|
|
|
return
|
|
end function is_neighbor
|
|
|
|
function is_equal(A, B)
|
|
!This function checks if too numbers are equal within a tolerance
|
|
real(kind=dp), intent(in) :: A, B
|
|
logical :: is_equal
|
|
|
|
if((A>(B - 10.0_dp**(-10))).and.(A < (B+10.0_dp**(-10)))) then
|
|
is_equal = .true.
|
|
else
|
|
is_equal = .false.
|
|
end if
|
|
return
|
|
end function is_equal
|
|
|
|
pure function unitvec(n,vec)
|
|
integer, intent(in) :: n
|
|
real(kind=dp), intent(in) :: vec(n)
|
|
real(kind=dp) :: unitvec(n)
|
|
|
|
unitvec = vec/norm2(vec)
|
|
return
|
|
end function unitvec
|
|
|
|
pure function norm_dis(rl,rk)
|
|
!This just returns the magnitude of the vector between two points
|
|
real(kind=dp), dimension(3), intent(in) :: rl, rk
|
|
real(kind=dp) :: norm_dis(4)
|
|
|
|
norm_dis(1:3) = (rk - rl)
|
|
norm_dis(4) = norm2(rk-rl)
|
|
end function
|
|
end module functions
|