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.
867 lines
33 KiB
867 lines
33 KiB
module elements
|
|
|
|
!This module contains the elements data structures, structures needed for building regions
|
|
!and operations that are done on elements
|
|
use parameters
|
|
use functions
|
|
use subroutines
|
|
use sorts
|
|
use box
|
|
implicit none
|
|
|
|
!Data structures used to represent the CAC elements. Each index represents an element
|
|
character(len=100), allocatable :: type_ele(:) !Element type
|
|
integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:), tag_ele(:) !Element size
|
|
real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
|
|
!Element result data structures
|
|
real(kind=dp), allocatable :: force_node(:,:,:,:), virial_node(:,:,:,:), energy_node(:,:,:)
|
|
|
|
integer, save :: ele_num !Number of elements
|
|
integer, save :: node_num !Total number of nodes
|
|
integer, save :: node_atoms !Count of all basis atoms at nodes summed over all nodes
|
|
|
|
!Data structure used to represent atoms
|
|
integer, allocatable :: type_atom(:)!atom type
|
|
integer, allocatable :: sbox_atom(:), tag_atom(:)
|
|
real(kind =dp),allocatable :: r_atom(:,:) !atom position
|
|
integer :: atom_num=0 !Number of atoms
|
|
!Atom result data structures information
|
|
real(kind=8), allocatable :: force_atom(:,:), virial_atom(:,:), energy_atom(:)
|
|
|
|
!Mapping atom type to provided name
|
|
character(len=2), dimension(10) :: type_to_name
|
|
integer :: atom_types = 0
|
|
|
|
!Variables for creating elements based on primitive cells
|
|
real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3)
|
|
integer :: cubic_faces(4,6)
|
|
|
|
!Below are lattice type arrays which provide information on the general form of the elements.
|
|
!We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased.
|
|
integer :: lattice_types = 0
|
|
integer :: max_ng_node, ng_node(10) !Max number of nodes per element and number of nodes per element for each lattice type
|
|
integer :: max_esize=0 !Maximum number of atoms per side of element
|
|
real(kind=dp) :: lapa(10)
|
|
|
|
!These variables contain information on the basis, for simplicities sake we limit
|
|
!the user to the definition of 10 lattice types with 10 basis atoms at each lattice point.
|
|
!This can be easily increased with no change to efficiency
|
|
integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type
|
|
integer :: basis_type(10,10)
|
|
|
|
!Additional module level variables we need
|
|
logical :: wrap_flag
|
|
|
|
!flags for data variables
|
|
logical :: vflag
|
|
public
|
|
contains
|
|
|
|
subroutine lattice_init
|
|
!This subroutine just intializes variables needed for building the different finite
|
|
!element types
|
|
|
|
!First initialize the cubic cell
|
|
cubic_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
|
|
1.0_dp, 0.0_dp, 0.0_dp, &
|
|
1.0_dp, 1.0_dp, 0.0_dp, &
|
|
0.0_dp, 1.0_dp, 0.0_dp, &
|
|
0.0_dp, 0.0_dp, 1.0_dp, &
|
|
1.0_dp, 0.0_dp, 1.0_dp, &
|
|
1.0_dp, 1.0_dp, 1.0_dp, &
|
|
0.0_dp, 1.0_dp, 1.0_dp /), &
|
|
shape(fcc_cell))
|
|
|
|
!Now we create a list containing the list of vertices needed to describe the 6 cube faces
|
|
cubic_faces(:,1) = (/ 1, 2, 3, 4 /)
|
|
cubic_faces(:,2) = (/ 1, 2, 6, 5 /)
|
|
cubic_faces(:,3) = (/ 2, 3, 7, 6 /)
|
|
cubic_faces(:,4) = (/ 3, 4, 8, 7 /)
|
|
cubic_faces(:,5) = (/ 1, 4, 8, 5 /)
|
|
cubic_faces(:,6) = (/ 5, 6, 7, 8 /)
|
|
|
|
!!Now initialize the fcc primitive cell
|
|
fcc_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
|
|
0.5_dp, 0.5_dp, 0.0_dp, &
|
|
0.5_dp, 1.0_dp, 0.5_dp, &
|
|
0.0_dp, 0.5_dp, 0.5_dp, &
|
|
0.5_dp, 0.0_dp, 0.5_dp, &
|
|
1.0_dp, 0.5_dp, 0.5_dp, &
|
|
1.0_dp, 1.0_dp, 1.0_dp, &
|
|
0.5_dp, 0.5_dp, 1.0_dp /), &
|
|
shape(fcc_cell))
|
|
|
|
fcc_mat = reshape((/ 0.5_dp, 0.5_dp, 0.0_dp, &
|
|
0.0_dp, 0.5_dp, 0.5_dp, &
|
|
0.5_dp, 0.0_dp, 0.5_dp /), &
|
|
shape(fcc_mat))
|
|
|
|
!Initialize the bcc primitive cell
|
|
bcc_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
|
|
0.5_dp, -0.5_dp, 0.5_dp, &
|
|
1.0_dp, 0.0_dp, 1.0_dp, &
|
|
0.5_dp, 0.5_dp, 0.5_dp, &
|
|
-0.5_dp, 0.5_dp, 0.5_dp, &
|
|
0.0_dp, 0.0_dp, 1.0_dp, &
|
|
0.5_dp, 0.5_dp, 1.5_dp, &
|
|
0.0_dp, 1.0_dp, 1.0_dp /), &
|
|
shape(bcc_cell))
|
|
|
|
bcc_mat = reshape((/ 0.5_dp, -0.5_dp, 0.5_dp, &
|
|
0.5_dp, 0.5_dp, 0.5_dp, &
|
|
-0.5_dp, 0.5_dp, 0.5_dp /), &
|
|
shape(bcc_mat))
|
|
|
|
|
|
|
|
call matrix_inverse(fcc_mat,3,fcc_inv)
|
|
call matrix_inverse(bcc_mat,3,bcc_inv)
|
|
|
|
max_basisnum = 0
|
|
basisnum(:) = 0
|
|
ng_node(:) = 0
|
|
node_num = 0
|
|
node_atoms = 0
|
|
ele_num = 0
|
|
atom_num = 0
|
|
end subroutine lattice_init
|
|
|
|
subroutine cell_init(lapa,esize,ele_type, orient_mat, cell_mat)
|
|
!This subroutine uses the user provided information to transform the finite element cell to the correct
|
|
!size, orientation, and dimensions using the esize, lattice parameter, element_type, and orientation provided
|
|
!by the user
|
|
|
|
integer, intent(in) :: esize
|
|
real(kind=dp), intent(in) :: lapa, orient_mat(3,3)
|
|
character(len=100), intent(in) :: ele_type
|
|
|
|
real(kind=dp), dimension(3,max_ng_node), intent(out) :: cell_mat
|
|
|
|
integer :: inod, i
|
|
real(kind=dp), dimension(3,max_ng_node) :: adjustVar
|
|
|
|
adjustVar(:,:) = 0.0_dp
|
|
vflag = .false.
|
|
|
|
select case(trim(ele_type))
|
|
|
|
case('fcc')
|
|
if(lmpcac) then
|
|
do inod = 1, 8
|
|
do i = 1,3
|
|
if(is_equal(cubic_cell(i, inod),0.0_dp)) then
|
|
adjustVar(i,inod) = -0.5_dp
|
|
else
|
|
adjustVar(i, inod) = 0.5_dp
|
|
end if
|
|
end do
|
|
end do
|
|
adjustVar(:,1:8) = matmul(fcc_mat, adjustVar(:,1:8))
|
|
end if
|
|
cell_mat(:, 1:8) = fcc_cell + adjustVar(:,1:8)
|
|
cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, cell_mat(:,1:8)))
|
|
case('bcc')
|
|
cell_mat(:,1:8) = bcc_cell
|
|
cell_mat(:,1:8) = lapa* ((esize-1)*matmul(orient_mat, cell_mat(:,1:8)))
|
|
case default
|
|
print *, "Element type ", trim(ele_type), " currently not accepted"
|
|
stop
|
|
end select
|
|
|
|
end subroutine cell_init
|
|
|
|
subroutine alloc_ele_arrays(n,m)
|
|
!This subroutine used to provide initial allocation for the atom and element arrays
|
|
|
|
integer, intent(in) :: n, m !n-size of element arrays, m-size of atom arrays
|
|
|
|
integer :: allostat
|
|
|
|
!Allocate element arrays
|
|
if(n > 0) then
|
|
allocate(type_ele(n), tag_ele(n), size_ele(n), lat_ele(n), sbox_ele(n), r_node(3,max_basisnum, max_ng_node,n), &
|
|
stat=allostat)
|
|
if(allostat > 0) then
|
|
print *, "Error allocating element arrays in elements.f90 because of: ", allostat
|
|
stop
|
|
end if
|
|
end if
|
|
|
|
if(m > 0) then
|
|
!Allocate atom arrays
|
|
allocate(type_atom(m), sbox_atom(m), tag_atom(m), r_atom(3,m), stat=allostat)
|
|
if(allostat > 0) then
|
|
print *, "Error allocating atom arrays in elements.f90 because of: ", allostat
|
|
stop
|
|
end if
|
|
end if
|
|
end subroutine
|
|
|
|
subroutine grow_ele_arrays(n, m)
|
|
|
|
integer, intent(in) :: n, m
|
|
|
|
integer :: ele_size, atom_size, buffer_size
|
|
integer, allocatable :: temp_int(:)
|
|
real(kind=dp), allocatable :: temp_ele_real(:,:,:,:), temp_real(:,:)
|
|
character(len=100), allocatable :: char_temp(:)
|
|
|
|
!The default size we grow the
|
|
buffer_size = 1024
|
|
|
|
!First check to make sure if it is allocated
|
|
if (allocated(size_ele)) then
|
|
|
|
!Figure out the size of the atom and element arrays
|
|
ele_size = size(size_ele)
|
|
|
|
!Check if we need to grow the ele_size, if so grow all the variables
|
|
if ( n+ele_num > size(size_ele)) then
|
|
|
|
allocate(temp_int(n+ele_size+buffer_size))
|
|
temp_int(1:ele_size) = lat_ele(1:ele_size)
|
|
temp_int(ele_size+1:) = 0
|
|
call move_alloc(temp_int, lat_ele)
|
|
|
|
allocate(temp_int(n+ele_size+buffer_size))
|
|
temp_int(1:ele_size) = tag_ele(1:ele_size)
|
|
temp_int(ele_size+1:) = 0
|
|
call move_alloc(temp_int, tag_ele)
|
|
|
|
allocate(temp_int(n+ele_size+buffer_size))
|
|
temp_int(1:ele_size) = size_ele(1:ele_size)
|
|
temp_int(ele_size+1:) = 0
|
|
call move_alloc(temp_int, size_ele)
|
|
|
|
allocate(temp_int(n+ele_size+buffer_size))
|
|
temp_int(1:ele_size) = sbox_ele(1:ele_size)
|
|
temp_int(ele_size+1:) = 0
|
|
call move_alloc(temp_int, sbox_ele)
|
|
|
|
allocate(char_temp(n+ele_size+buffer_size))
|
|
char_temp(1:ele_size) = type_ele(1:ele_size)
|
|
call move_alloc(char_temp, type_ele)
|
|
|
|
allocate(temp_ele_real(3, max_basisnum, max_ng_node, n+ele_size+buffer_size))
|
|
temp_ele_real(:,:,:,1:ele_size) = r_node(:,:,:,1:ele_size)
|
|
temp_ele_real(:,:,:,ele_size+1:) = 0.0_dp
|
|
call move_alloc(temp_ele_real, r_node)
|
|
end if
|
|
else
|
|
call alloc_ele_arrays(n,0)
|
|
end if
|
|
|
|
!Now grow atom arrays if needed
|
|
if (allocated(type_atom)) then
|
|
atom_size = size(type_atom)
|
|
if (m+atom_num > atom_size) then
|
|
allocate(temp_int(m+atom_size+buffer_size))
|
|
temp_int(1:atom_size) = type_atom
|
|
temp_int(atom_size+1:) = 0
|
|
call move_alloc(temp_int, type_atom)
|
|
|
|
allocate(temp_int(m+atom_size+buffer_size))
|
|
temp_int(1:atom_size) = tag_atom
|
|
temp_int(atom_size+1:) = 0
|
|
call move_alloc(temp_int, tag_atom)
|
|
|
|
allocate(temp_int(m+atom_size+buffer_size))
|
|
temp_int(1:atom_size) = sbox_atom
|
|
temp_int(atom_size+1:) = 0
|
|
call move_alloc(temp_int, sbox_atom)
|
|
|
|
allocate(temp_real(3,m+atom_size+buffer_size))
|
|
temp_real(:,1:atom_size) = r_atom
|
|
temp_real(:, atom_size+1:) = 0.0_dp
|
|
call move_alloc(temp_real, r_atom)
|
|
end if
|
|
else
|
|
call alloc_ele_arrays(0,m)
|
|
end if
|
|
end subroutine
|
|
|
|
subroutine add_element(tag, type, size, lat, sbox, r)
|
|
!Subroutine which adds an element to the element arrays
|
|
integer, intent(in) :: size, lat, sbox, tag
|
|
character(len=100), intent(in) :: type
|
|
real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node)
|
|
|
|
integer :: newtag
|
|
|
|
ele_num = ele_num + 1
|
|
node_num = node_num + ng_node(lat)
|
|
node_atoms = node_atoms + ng_node(lat)*basisnum(lat)
|
|
|
|
if (tag==0) then
|
|
newtag = ele_num !If we don't assign a tag then pass the tag as the ele_num
|
|
else
|
|
newtag = tag
|
|
end if
|
|
|
|
!Check to see if we need to grow the arrays
|
|
call grow_ele_arrays(1,0)
|
|
tag_ele(ele_num) = newtag
|
|
type_ele(ele_num) = type
|
|
size_ele(ele_num) = size
|
|
lat_ele(ele_num) = lat
|
|
sbox_ele(ele_num) = sbox
|
|
r_node(:,:,:,ele_num) = r(:,:,:)
|
|
|
|
|
|
end subroutine add_element
|
|
|
|
subroutine add_atom(tag, type, sbox, r)
|
|
!Subroutine which adds an atom to the atom arrays
|
|
integer, intent(in) :: type, sbox, tag
|
|
real(kind=dp), intent(in), dimension(3) :: r
|
|
|
|
integer :: newtag
|
|
|
|
atom_num = atom_num+1
|
|
if(tag==0) then
|
|
newtag = atom_num !If we don't assign a tag then pass the tag as the atom_num
|
|
else
|
|
newtag = tag
|
|
end if
|
|
!Check to see if we need to grow the arrays
|
|
call grow_ele_arrays(0,1)
|
|
tag_atom(atom_num) = tag
|
|
type_atom(atom_num) = type
|
|
r_atom(:,atom_num) = r(:)
|
|
sbox_atom(atom_num) = sbox
|
|
|
|
end subroutine add_atom
|
|
|
|
subroutine add_atom_type(type, inttype)
|
|
!This subroutine adds a new atom type to the module list of atom types
|
|
character(len=2), intent(in) :: type
|
|
integer, intent(out) :: inttype
|
|
|
|
integer :: i
|
|
logical :: exists
|
|
|
|
exists = .false.
|
|
do i=1, 10
|
|
if(type == type_to_name(i)) then
|
|
exists = .true.
|
|
inttype = i
|
|
exit
|
|
end if
|
|
end do
|
|
|
|
if (exists.eqv..false.) then
|
|
atom_types = atom_types+1
|
|
if(atom_types > 10) then
|
|
print *, "Defined atom types are greater than 10 which is currently not supported."
|
|
stop 3
|
|
end if
|
|
type_to_name(atom_types) = type
|
|
inttype = atom_types
|
|
end if
|
|
return
|
|
end subroutine add_atom_type
|
|
|
|
subroutine def_ng_node(n, element_types)
|
|
!This subroutine defines the maximum node number among n element types
|
|
integer, intent(in) :: n !Number of element types
|
|
character(len=100), dimension(n) :: element_types !Array of element type strings
|
|
|
|
integer :: i
|
|
|
|
do i=1, n
|
|
select case(trim(adjustl(element_types(i))))
|
|
case('fcc')
|
|
ng_node(i) = 8
|
|
case('bcc')
|
|
ng_node(i) = 8
|
|
end select
|
|
|
|
if(ng_node(i) > max_ng_node) max_ng_node = ng_node(i)
|
|
end do
|
|
end subroutine def_ng_node
|
|
|
|
subroutine set_max_esize
|
|
!This subroutine sets the maximum esize
|
|
if(allocated(size_ele)) then
|
|
max_esize=maxval(size_ele)
|
|
else
|
|
max_esize = 2
|
|
end if
|
|
end subroutine
|
|
|
|
subroutine interpolate_atoms(type, esize, lat_type, r_in, type_interp, r_interp, eng, f, v, data_interp)
|
|
!This subroutine returns the interpolated atoms from the elements.
|
|
|
|
!Arguments
|
|
character(len=100), intent(in) :: type !The type of element that it is
|
|
integer, intent(in) :: esize !The number of atoms per side
|
|
integer, intent(in) :: lat_type !The integer lattice type of the element
|
|
real(kind=dp), dimension(3,max_basisnum, max_ng_node), intent(in) :: r_in !Nodal positions
|
|
integer, dimension(max_basisnum*max_esize**3), intent(out) :: type_interp !Interpolated atomic positions
|
|
real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: r_interp !Interpolated atomic positions
|
|
real(kind = dp), optional, intent(in) :: eng(max_basisnum, max_ng_node), f(3, max_basisnum, max_ng_node), &
|
|
v(6, max_basisnum, max_ng_node)
|
|
real(kind=dp), dimension(10, max_basisnum*max_esize**3), optional,intent(out) :: data_interp !Interpolated atomic positions
|
|
|
|
!Internal variables
|
|
integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp
|
|
real(kind=dp), allocatable :: a_shape(:)
|
|
real(kind=dp) :: t, s, r
|
|
|
|
!Initialize some variables
|
|
r_interp(:,:) = 0.0_dp
|
|
type_interp(:) = 0
|
|
if(present(data_interp)) data_interp = 0.0_dp
|
|
ia = 0
|
|
|
|
!Define bnum based on the input lattice type. If lat_type=0 then we are interpolating lattice points which means
|
|
!the basis is 0,0,0, and the type doesn't matter
|
|
|
|
select case(lat_type)
|
|
case(0)
|
|
bnum=1
|
|
lat_type_temp = 1
|
|
case default
|
|
bnum = basisnum(lat_type)
|
|
lat_type_temp = lat_type
|
|
end select
|
|
|
|
select case(trim(adjustl(type)))
|
|
case('fcc','bcc')
|
|
allocate(a_shape(8))
|
|
!Now loop over all the possible sites
|
|
do it = 1, esize
|
|
t = (1.0_dp*(it-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
|
|
do is =1, esize
|
|
s = (1.0_dp*(is-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
|
|
do ir = 1, esize
|
|
r = (1.0_dp*(ir-1) - (esize-1)/2)/(1.0_dp*(esize-1)/2)
|
|
call rhombshape(r,s,t,a_shape)
|
|
|
|
do ibasis = 1, bnum
|
|
ia = ia + 1
|
|
do inod = 1, 8
|
|
type_interp(ia) = basis_type(ibasis,lat_type_temp)
|
|
r_interp(:,ia) = r_interp(:,ia) + a_shape(inod) * r_in(:,ibasis,inod)
|
|
|
|
if(present(data_interp)) then
|
|
!If data is present then interpolate data arrays as well
|
|
data_interp(1,ia) = data_interp(1,ia) + eng(ibasis, inod)*a_shape(inod)
|
|
data_interp(2:4,ia) = data_interp(2:4,ia) + f(:, ibasis, inod)*a_shape(inod)
|
|
data_interp(5:10,ia) = data_interp(5:10,ia) + v(:, ibasis, inod)*a_shape(inod)
|
|
end if
|
|
end do
|
|
end do
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
end select
|
|
|
|
if (ia /= bnum*esize**3) then
|
|
print *, "Incorrect interpolation"
|
|
stop 3
|
|
end if
|
|
return
|
|
end subroutine interpolate_atoms
|
|
|
|
subroutine rhombshape(r,s,t, shape_fun)
|
|
!Shape function for rhombohedral elements
|
|
real(kind=8), intent(in) :: r, s, t
|
|
real(kind=8), intent(out) :: shape_fun(8)
|
|
|
|
shape_fun(1) = (1.0-r)*(1.0-s)*(1.0-t)/8.0
|
|
shape_fun(2) = (1.0+r)*(1.0-s)*(1.0-t)/8.0
|
|
shape_fun(3) = (1.0+r)*(1.0+s)*(1.0-t)/8.0
|
|
shape_fun(4) = (1.0-r)*(1.0+s)*(1.0-t)/8.0
|
|
shape_fun(5) = (1.0-r)*(1.0-s)*(1.0+t)/8.0
|
|
shape_fun(6) = (1.0+r)*(1.0-s)*(1.0+t)/8.0
|
|
shape_fun(7) = (1.0+r)*(1.0+s)*(1.0+t)/8.0
|
|
shape_fun(8) = (1.0-r)*(1.0+s)*(1.0+t)/8.0
|
|
|
|
|
|
return
|
|
end subroutine rhombshape
|
|
|
|
subroutine delete_atoms(num, in_index)
|
|
!This subroutine deletes atoms from the atom arrays
|
|
integer, intent(in) :: num
|
|
integer, intent(in), dimension(num) :: in_index
|
|
|
|
real(kind=dp), dimension(num) :: for_sort
|
|
integer, dimension(num) :: sorted_index
|
|
|
|
integer :: i
|
|
|
|
for_sort = in_index
|
|
|
|
call dpquicksort(for_sort)
|
|
do i = 1, num
|
|
sorted_index(i) = nint(for_sort(i))
|
|
end do
|
|
|
|
!We go from largest index to smallest index just to make sure that we don't miss
|
|
!accidentally overwrite values which need to be deleted
|
|
do i = num, 1, -1
|
|
if(sorted_index(i) == atom_num) then
|
|
r_atom(:,sorted_index(i)) = 0.0_dp
|
|
type_atom(sorted_index(i)) = 0
|
|
else
|
|
r_atom(:,sorted_index(i)) = r_atom(:, atom_num)
|
|
type_atom(sorted_index(i)) = type_atom(atom_num)
|
|
end if
|
|
atom_num = atom_num - 1
|
|
end do
|
|
end subroutine
|
|
|
|
subroutine delete_elements(num, in_index)
|
|
!This subroutine deletes elements from the element array
|
|
integer, intent(in) :: num
|
|
integer, intent(in), dimension(num) :: in_index
|
|
|
|
integer :: i
|
|
real(kind=dp), dimension(num) :: for_sort
|
|
integer, dimension(num) :: sorted_index
|
|
|
|
|
|
for_sort = in_index
|
|
call dpquicksort(for_sort)
|
|
do i = 1, num
|
|
sorted_index(i) = nint(for_sort(i))
|
|
end do
|
|
|
|
!We go from largest index to smallest index just to make sure that we don't miss
|
|
!accidentally overwrite values which need to be deleted
|
|
do i = num, 1, -1
|
|
node_num = node_num - ng_node(lat_ele(sorted_index(i)))
|
|
if(sorted_index(i) == ele_num) then
|
|
r_node(:,:,:,sorted_index(i)) = 0.0_dp
|
|
type_ele(sorted_index(i)) =''
|
|
size_ele(sorted_index(i)) = 0
|
|
lat_ele(sorted_index(i)) = 0
|
|
sbox_ele(sorted_index(i)) = 0
|
|
tag_ele(sorted_index(i)) = 0
|
|
else
|
|
r_node(:,:,:,sorted_index(i)) = r_node(:,:,:,ele_num)
|
|
type_ele(sorted_index(i)) = type_ele(ele_num)
|
|
size_ele(sorted_index(i)) = size_ele(ele_num)
|
|
lat_ele(sorted_index(i)) = lat_ele(ele_num)
|
|
sbox_ele(sorted_index(i)) = sbox_ele(ele_num)
|
|
tag_ele(sorted_index(i)) = tag_ele(ele_num)
|
|
end if
|
|
ele_num = ele_num - 1
|
|
end do
|
|
|
|
end subroutine delete_elements
|
|
|
|
subroutine wrap_atoms
|
|
!This subroutine wraps atoms back into the simulation cell if they have exited for any reason
|
|
integer :: i
|
|
do i = 1, atom_num
|
|
call apply_periodic(r_atom(:,i))
|
|
end do
|
|
end subroutine wrap_atoms
|
|
|
|
subroutine def_new_box
|
|
!This subroutine calculates new box boundaries based on minimum and maximum nodal/atomic positions
|
|
integer :: i, j, inod, ibasis
|
|
|
|
real(kind=dp) :: max_bd(3), min_bd(3)
|
|
|
|
max_bd(:) = -huge(1.0_dp)
|
|
min_bd(:) = huge(1.0_dp)
|
|
do i = 1, atom_num
|
|
do j = 1, 3
|
|
if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + tol
|
|
if (r_atom(j,i) < min_bd(j)) min_bd(j) = r_atom(j,i) - tol
|
|
end do
|
|
end do
|
|
|
|
do i = 1, ele_num
|
|
do inod = 1, ng_node(lat_ele(i))
|
|
do ibasis = 1, basisnum(lat_ele(i))
|
|
do j = 1, 3
|
|
if (r_node(j,ibasis,inod,i) > max_bd(j)) max_bd(j) = r_node(j,ibasis,inod,i) + tol
|
|
if (r_node(j,ibasis,inod,i) < min_bd(j)) min_bd(j) = r_node(j,ibasis,inod,i) - tol
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
do j = 1, 3
|
|
if(box_bc(j:j) == 's') then
|
|
box_bd(2*j) = max_bd(j)
|
|
box_bd(2*j-1) = min_bd(j)
|
|
end if
|
|
end do
|
|
end subroutine
|
|
|
|
recursive subroutine parse_pos(i, pos_string, pos)
|
|
!This subroutine parses the pos command allowing for command which include inf
|
|
integer, intent(in) :: i !The dimension of the position
|
|
character(len=100), intent(in) :: pos_string !The position string
|
|
real(kind=dp), intent(out) :: pos !The output parsed position value
|
|
|
|
integer :: iospara, face, ele, randsize
|
|
real(kind=dp) :: rand, rone, rtwo, rand_ele_pos(3)
|
|
character(len=100) :: cone, ctwo
|
|
|
|
iospara = 0
|
|
|
|
if(trim(adjustl(pos_string)) == 'inf') then
|
|
pos=box_bd(2*i)
|
|
|
|
else if(trim(adjustl(pos_string)) == '-inf') then
|
|
pos=box_bd(2*i-1)
|
|
|
|
else if (trim(adjustl(pos_string)) == 'rand') then
|
|
call random_number(rand)
|
|
pos = (box_bd(2*i)-box_bd(2*i-1))*rand + box_bd(2*i-1)
|
|
|
|
else if (index(pos_string,'rande')>0) then
|
|
!First select a random element
|
|
call random_number(rand)
|
|
ele= 1 + floor(ele_num*rand)
|
|
!Now read the rest of the command which specifies the face we need and check
|
|
!to make sure it's an accepted face number
|
|
cone = pos_string(index(pos_string, '[')+1:index(pos_string, '[')+1)
|
|
|
|
!Check to see if we also pass an element size, if it was passed then make sure our
|
|
!random element is the right size
|
|
if(index(pos_string,':') > 0) then
|
|
ctwo = pos_string(index(pos_string, ':')+1:index(pos_string,']')-1)
|
|
read(ctwo, *) randsize
|
|
do while(randsize /= size_ele(ele))
|
|
call random_number(rand)
|
|
ele= 1 + floor(ele_num*rand)
|
|
end do
|
|
end if
|
|
|
|
read(cone, *) face
|
|
if ((face < 1).or.(face > 6)) stop "Current face number must be 1,2,3,4,5,6. Please check documentation"
|
|
!Now get the position
|
|
call offset_pos(ele, face, rand_ele_pos)
|
|
pos = rand_ele_pos(i)
|
|
|
|
else if (index(pos_string,'rand')>0) then
|
|
call random_number(rand)
|
|
cone = pos_string(index(pos_string, '[')+1:index(pos_string,':')-1)
|
|
call parse_pos(i, cone, rone)
|
|
ctwo = pos_string(index(pos_string, ':')+1:index(pos_string,']')-1)
|
|
call parse_pos(i, ctwo, rtwo)
|
|
pos = (rtwo - rone)*rand + rone
|
|
|
|
else if ((index(pos_string,'-') > 0).and.(index(pos_string,'inf')>0)) then
|
|
!Now extract the number we are reducing from infinity
|
|
if(index(pos_string,'inf') < index(pos_string,'-')) then
|
|
read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
|
|
else
|
|
read(pos_string(1:index(pos_string,'-')-1), *, iostat=iospara) pos
|
|
end if
|
|
pos = box_bd(2*i) - pos
|
|
|
|
else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then
|
|
!Now extract the number we are reducing from infinity
|
|
if(index(pos_string,'inf') < index(pos_string,'+')) then
|
|
read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
|
|
else
|
|
read(pos_string(1:index(pos_string,'+')-1), *, iostat=iospara) pos
|
|
end if
|
|
pos = box_bd(2*i-1) + pos
|
|
|
|
else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then
|
|
!Now extract the number we are reducing from infinity
|
|
if(index(pos_string,'inf') < index(pos_string,'*')) then
|
|
read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos
|
|
else
|
|
read(pos_string(1:index(pos_string,'*')-1), *, iostat=iospara) pos
|
|
end if
|
|
pos = (box_bd(2*i)-box_bd(2*i-1))*pos + box_bd(2*i-1)
|
|
|
|
else
|
|
read(pos_string, *, iostat=iospara) pos
|
|
|
|
end if
|
|
|
|
if (iospara > 0) then
|
|
print *, "Error reading position argument ", trim(adjustl(pos_string)), ". Please reformat and try again."
|
|
end if
|
|
end subroutine parse_pos
|
|
|
|
subroutine offset_pos(ie, iface, pos)
|
|
!This returns a position slightly offset from the center of an element face
|
|
!This is used to return a random position within an element discontinuity
|
|
|
|
integer, intent(in) :: ie !Element index
|
|
integer, intent(in) :: iface !Face number between 1 and 6
|
|
real(kind=dp), dimension(3), intent(out) :: pos !Position vector
|
|
|
|
!Other variables we need
|
|
integer :: esize
|
|
real(kind=dp) :: orient(3,3), ori_inv(3,3), lp, r_cubic_node(3,8)
|
|
|
|
!First find the offset vector for the face in the untransformed cubic cell
|
|
esize = size_ele(ie)
|
|
select case(iface)
|
|
case(1)
|
|
pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp) /)
|
|
case(2)
|
|
pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /)
|
|
case(3)
|
|
pos = (/ (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
|
|
case(4)
|
|
pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /)
|
|
case(5)
|
|
pos = (/ -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
|
|
case(6)
|
|
pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp) /)
|
|
end select
|
|
|
|
!Now transform it to real space and adjust it to the position of the element in the first node.
|
|
select case (trim(adjustl(type_ele(ie))))
|
|
|
|
case('fcc')
|
|
!First we have to extract the element lattice parameter
|
|
call matrix_inverse(sub_box_ori(:,:,sbox_ele(ie)),3,ori_inv)
|
|
r_cubic_node = r_node(:,1,:,ie)
|
|
r_cubic_node = matmul(fcc_inv,matmul(ori_inv,r_cubic_node))
|
|
lp = (maxval(r_cubic_node(1,:))-minval(r_cubic_node(1,:)))/(esize-1)
|
|
pos = matmul(sub_box_ori(:,:,sbox_ele(ie)),matmul(fcc_mat,pos))*lp
|
|
pos = pos + r_node(:,1, 1, ie)
|
|
|
|
case default
|
|
print *, trim(adjustl(type_ele(ie))), " is not a currently accepted element type for random element positions"
|
|
stop 3
|
|
|
|
end select
|
|
end subroutine
|
|
|
|
subroutine lattice_map(in_bnum, in_btypes, in_ngnodes, in_lapa, lat_type)
|
|
!This subroutine maps an input lattice type to either a new lattice type or an existing one depending on basis_type and
|
|
!number of nodes at the atoms
|
|
|
|
integer, intent(in) :: in_ngnodes, in_bnum, in_btypes(10) !Input variables
|
|
real(kind=dp), intent(in) :: in_lapa
|
|
integer, intent(out) :: lat_type
|
|
|
|
integer j, ibasis
|
|
|
|
lat_type = 0
|
|
lat_loop:do j = 1, lattice_types
|
|
!Check all the lattice level variables
|
|
if ((basisnum(j) == in_bnum).and.(ng_node(j) == in_ngnodes).and.(is_equal(lapa(j),in_lapa))) then
|
|
!Now check lattice level variables
|
|
do ibasis = 1, basisnum(j)
|
|
if(basis_type(ibasis,j) /= in_btypes(ibasis)) cycle lat_loop
|
|
end do
|
|
lat_type = j
|
|
exit lat_loop
|
|
end if
|
|
end do lat_loop
|
|
|
|
!If it doesn't match an existing lattice type we add it
|
|
if( lat_type == 0) then
|
|
lattice_types = lattice_types + 1
|
|
basisnum(lattice_types) = in_bnum
|
|
basis_type(:,lattice_types) = in_btypes
|
|
ng_node(lattice_types) = in_ngnodes
|
|
lapa(lattice_types) = in_lapa
|
|
lat_type = lattice_types
|
|
end if
|
|
|
|
end subroutine lattice_map
|
|
|
|
subroutine get_interp_pos(i,j,k, ie, rout)
|
|
!This returns the position of an interpolated basis from an element ie.
|
|
!i, j, k should be in natural coordinates
|
|
|
|
integer, intent(in) :: i, j, k
|
|
real(kind=dp), dimension(3,max_basisnum), intent(out) :: rout
|
|
|
|
integer :: ie, ibasis, inod
|
|
real(kind=dp) :: a_shape(8), r, s, t
|
|
|
|
r = (1.0_dp*(i-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2)
|
|
s = (1.0_dp*(j-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2)
|
|
t = (1.0_dp*(k-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2)
|
|
rout(:,:) = 0
|
|
do ibasis = 1, basisnum(lat_ele(ie))
|
|
do inod = 1, ng_node(lat_ele(ie))
|
|
call rhombshape(r,s,t,a_shape)
|
|
rout(:,ibasis) = rout(:,ibasis) + a_shape(inod) * r_node(:,ibasis,inod,ie)
|
|
end do
|
|
end do
|
|
|
|
end subroutine get_interp_pos
|
|
|
|
subroutine alloc_dat_arrays(n,m)
|
|
!This subroutine used to provide initial allocation for the atom and element data arrays
|
|
integer, intent(in) :: n,m !n-size of element arrays, m-size of atom arrays
|
|
integer :: allostat
|
|
|
|
!Allocate element arrays
|
|
if (n > 0) then
|
|
if (allocated(force_node)) then
|
|
deallocate(force_node, virial_node, energy_node)
|
|
end if
|
|
allocate(force_node(3,max_basisnum, max_ng_node, n), &
|
|
virial_node(6,max_basisnum, max_ng_node, n), &
|
|
energy_node(max_basisnum,max_ng_node,n), &
|
|
stat=allostat)
|
|
if(allostat > 0) then
|
|
print *, "Error allocating element data arrays in mode_metric because of:", allostat
|
|
stop
|
|
end if
|
|
|
|
end if
|
|
|
|
if (m > 0) then
|
|
if (allocated(force_atom)) then
|
|
deallocate(force_atom, virial_atom, energy_atom)
|
|
end if
|
|
allocate(force_atom(3, m), &
|
|
virial_atom(6, m), &
|
|
energy_atom(m), &
|
|
stat=allostat)
|
|
if(allostat > 0) then
|
|
print *, "Error allocating atom data arrays in mode_metric because of:", allostat
|
|
stop
|
|
end if
|
|
end if
|
|
|
|
end subroutine
|
|
|
|
subroutine add_atom_data(ia, eng, force, virial)
|
|
!Function which sets the atom data arrays
|
|
integer, intent(in) :: ia
|
|
real(kind=dp), intent(in) :: eng, force(3), virial(6)
|
|
|
|
energy_atom(ia) = eng
|
|
force_atom(:,ia) = force(:)
|
|
virial_atom(:,ia) = virial(:)
|
|
vflag = .true.
|
|
return
|
|
end subroutine add_atom_data
|
|
|
|
subroutine add_element_data(ie, eng, force, virial)
|
|
!Function which sets the element data arrays
|
|
integer, intent(in) :: ie
|
|
real(kind=dp), intent(in) :: eng(max_basisnum, max_ng_node), &
|
|
force(3,max_basisnum, max_ng_node), &
|
|
virial(6,max_basisnum,max_ng_node)
|
|
vflag = .true.
|
|
energy_node(:,:,ie) = eng
|
|
force_node(:,:,:,ie) = force
|
|
virial_node(:,:,:,ie) = virial
|
|
return
|
|
end subroutine add_element_data
|
|
|
|
subroutine reset_data
|
|
!This function resets all of the data arrays for the elements and atoms
|
|
atom_num = 0
|
|
ele_num = 0
|
|
node_num = 0
|
|
end subroutine reset_data
|
|
|
|
end module elements
|