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.
CACmb/src/mode_create.f90

524 lines
22 KiB

module mode_create
!This mode is intended for creating element/atom regions and writing them to specific files.
use parameters
use atoms
use io
use subroutines
use elements
use box
implicit none
character(len=100) :: name, element_type
real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), &
orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3)
integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
basis_pos(3,10)
logical :: dup_flag, dim_flag
real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
public
contains
subroutine create(arg_pos)
! Main subroutine which controls execution
integer, intent(out) :: arg_pos
integer :: i, ibasis, inod
real(kind=dp), allocatable :: r_node_temp(:,:,:)
print *, '-----------------------Mode Create---------------------------'
!Initialize default parameters
orient = reshape((/ 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp /), shape(orient))
cell_mat(:,:)=0.0_dp
name =''
element_type = ''
esize=0
lattice_parameter=0.0_dp
duplicate(:) = 0
box_len(:) = 0.0_dp
dup_flag = .false.
dim_flag = .false.
basisnum = 0
lat_ele_num = 0
lat_atom_num = 0
!First we parse the command
call parse_command(arg_pos)
!Now we setup the unit element and call other init subroutines
call def_ng_node(1, element_type)
allocate(r_node_temp(3,max_basisnum,max_ng_node))
!Get the inverse orientation matrix we will need later
call matrix_inverse(orient,3,orient_inv)
if(dup_flag) then
!Define box vertices
do i = 1, 8
box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + (origin(:)/lattice_parameter)
end do
!Now get the rotated box vertex positions in lattice space. Should be integer units
box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
!Find the new maxlen
maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2)
do i = 1, 3
box_bd(2*i) = maxval(box_vert(i,:)) - 0.25_dp*lattice_space(i)
box_bd(2*i-1) = origin(i)-0.25_dp*lattice_space(i)
end do
else if(dim_flag) then
!As a note everything is defined so that the lattice parameter is multiplied in at the end
!so we have to divide all the real Angstroms units by the lattice parameter
!Define box_vertices
do i = 1, 8
box_vert(:,i) = (cubic_cell(:,i)*box_len(:) + origin(:))/lattice_parameter
end do
!Now get the rotated box vertex positions in lattice space. Should be integer units
box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
!Now get the box_bd in lattice units
do i = 1, 3
box_bd(2*i) = (box_len(i)+origin(i))/lattice_parameter
box_bd(2*i-1) = origin(i)/lattice_parameter
end do
else
print *, "Creating 1 element"
call cell_init(lattice_parameter, esize, element_type, orient, cell_mat)
!If the user doesn't pass any build instructions than we just put the cell mat into the element_array
call alloc_ele_arrays(1,0)
!Add the basis atoms to the unit cell
do inod = 1, max_ng_node
do ibasis = 1, basisnum(1)
r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis) + origin(:)
end do
end do
do i = 1,3
box_bd(2*i) = maxval(r_node_temp(i,:,:))
box_bd(2*i-1) = origin(i)
end do
call add_element(element_type, esize, 1, 1, r_node_temp)
end if
!If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays
if(dup_flag.or.dim_flag) then
!Call the build function with the correct transformation matrix
select case(trim(adjustl(element_type)))
case('fcc')
call build_with_rhomb(box_lat_vert, fcc_mat)
case default
print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", &
"element type"
stop 3
end select
!Now that it is built multiply by the lattice parameter
box_bd = box_bd*lattice_parameter
print *, "Using mode create, ", lat_ele_num, " elements are created and ", lat_atom_num*basisnum(1), &
" atoms are created."
!Allocate variables
call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1))
if(lat_atom_num > 0) then
do i = 1, lat_atom_num
do ibasis = 1, basisnum(1)
call add_atom(basis_type(ibasis, 1), 1, (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
end do
end do
deallocate(r_atom_lat)
end if
if(lat_ele_num > 0) then
do i = 1, lat_ele_num
do inod= 1, ng_node(1)
do ibasis = 1, basisnum(1)
r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis)
end do
end do
call add_element(element_type, esize, 1, 1, r_node_temp)
end do
end if
end if
!The last thing we do is setup the sub_box_boundaries
call alloc_sub_box(1)
sub_box_num = 1
sub_box_ori(:,:,1) = orient
sub_box_bd(:,1) = box_bd
sub_box_array_bd(1,:,1) = 1
sub_box_array_bd(2,1,1) = atom_num
sub_box_array_bd(2,2,1) = ele_num
end subroutine create
!This subroutine parses the command and pulls out information needed for mode_create
subroutine parse_command(arg_pos)
integer, intent(out) :: arg_pos
integer :: ori_pos, i, j, arglen, stat
character(len=100) :: textholder
character(len=8) :: orient_string
logical :: isortho, isrighthanded
!Pull out all required positional arguments
call get_command_argument(2, name, arglen)
if(arglen==0) STOP "Name is missing in mode create"
call get_command_argument(3,element_type, arglen)
if(arglen==0) STOP "Element_type is missing in mode create"
call get_command_argument(4,textholder, arglen)
if(arglen==0) STOP "Lattice Parameter is missing in mode create"
read(textholder, *, iostat=stat) lattice_parameter
if(stat > 0) STOP "Error reading lattice parameter"
call get_command_argument(5, textholder, arglen)
if(arglen==0) STOP "Esize missing in mode create"
read(textholder, *, iostat=stat) esize
if(stat > 0) STOP "Error reading esize"
arg_pos = 6
!Check for optional keywords
do while(.true.)
if(arg_pos > command_argument_count()) exit
!Pull out the next argument which should either be a keyword or an option
call get_command_argument(arg_pos, textholder)
textholder=adjustl(textholder)
arg_pos=arg_pos+1
!Choose what to based on what the option string is
select case(trim(textholder))
!If orient command is passed extract the orientation to numeric array format
case('orient')
do i = 1, 3
call get_command_argument(arg_pos, orient_string, arglen)
if (arglen==0) STOP "Missing orientation in orient command of mode create"
call parse_ori_vec(orient_string, orient(i,:))
arg_pos = arg_pos+1
! ori_pos=2
! do j = 1,3
! if (orient_string(ori_pos:ori_pos) == '-') then
! ori_pos = ori_pos + 1
! read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
! if (stat>0) STOP "Error reading orient value"
! orient(i,j) = -orient(i,j)
! ori_pos = ori_pos + 1
! else
! read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
! if(stat>0) STOP "Error reading orient value"
! ori_pos=ori_pos + 1
! end if
! end do
end do
!If the duplicate command is passed then we extract the information on the new bounds.
case('duplicate')
if(dim_flag) STOP "Both duplicate and dim options cannot be used in mode_create"
dup_flag = .true.
do i = 1, 3
call get_command_argument(arg_pos, textholder)
read(textholder, *) duplicate(i)
arg_pos = arg_pos + 1
end do
case('dim')
if(dup_flag) STOP "Both duplicate and dim options cannot be used in mode_create"
dim_flag = .true.
do i = 1, 3
call get_command_argument(arg_pos, textholder)
read(textholder, *) box_len(i)
arg_pos = arg_pos + 1
end do
case('origin')
do i = 1, 3
call get_command_argument(arg_pos, textholder)
read(textholder, *) origin(i)
arg_pos = arg_pos + 1
end do
case default
!If it isn't an option then you have to exit
arg_pos = arg_pos -1
exit
end select
end do
!Calculate the lattice periodicity length in lattice units
do i = 1, 3
lattice_space(i) = norm2(orient(i,:))
end do
!Check special periodicity relations
select case(trim(adjustl(element_type)))
case('fcc')
do i = 1,3
!Check if one of the directions is 110
if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(orient(i,3),0.0_dp))).or.&
(is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(orient(i,1),0.0_dp))).or.&
(is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(orient(i,2),0.0_dp)))) then
lattice_space(i) = 0.5_dp * lattice_space(i)
!Check if one direction is 112
else if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(abs(orient(i,3)),2.0_dp*abs(orient(i,1))))).or.&
(is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(abs(orient(i,1)),2.0_dp*abs(orient(i,2))))).or.&
(is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(abs(orient(i,2)),2.0_dp*abs(orient(i,3))))))&
then
lattice_space(i) = 0.5_dp * lattice_space(i)
end if
end do
end select
!Now normalize the orientation matrix
orient = matrix_normal(orient,3)
!Now check these to make sure they are right handed and orthogonal
call check_right_ortho(orient, isortho, isrighthanded)
if (.not.isortho) then
stop "Directions in orient are not orthogonal"
else if (.not.isrighthanded) then
stop "Directions in orient are not righthanded"
end if
!Set lattice_num to 1 and add the lattice_parameter to the elements module lattice paramter variable
lattice_types = 1
lapa(1) = lattice_parameter
!If we haven't defined a basis then define the basis and add the default basis atom type and position
if (basisnum(1) == 0) then
max_basisnum = 1
basisnum(1) = 1
call add_atom_type(name, basis_type(1,1)) !If basis command not defined then we use name as the atom_name
basis_pos(:,1) = 0.0_dp
end if
!
end subroutine
subroutine build_with_rhomb(box_in_lat, transform_matrix)
!This subroutine returns all the lattice points in the box in r_lat
!Inputs
integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space
real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space
!Internal variables
integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,8), rzero(3), &
vlat(3), temp_lat(3,8), m, n, o
real(kind=dp) :: v(3), temp_nodes(3,1,8)
logical, allocatable :: lat_points(:,:,:)
logical :: node_in_bd(8)
!Do some value initialization
max_esize = esize
!First find the bounding lattice points (min and max points for the box in each dimension)
numlatpoints = 1
do i = 1, 3
bd_in_lat(2*i-1) = minval(box_in_lat(i,:))
bd_in_lat(2*i) = maxval(box_in_lat(i,:))
numlatpoints = numlatpoints*(bd_in_lat(2*i)-bd_in_lat(2*i-1))
end do
!Allocate the correct lat variable
select case(esize)
!Atomistics
case(2)
allocate(r_atom_lat(3,numlatpoints))
case default
continue
end select
!Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case
!and one for the regular case
if (esize==2) then
!atomistics
do iz = bd_in_lat(5)-5, bd_in_lat(6)+5
do iy = bd_in_lat(3)-5, bd_in_lat(4)+5
do ix = bd_in_lat(1)-5, bd_in_lat(2)+5
v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /)
!Transform point back to real space for easier checking
v = matmul(orient, matmul(transform_matrix,v))
!If within the boundaries
if(in_block_bd(v, box_bd)) then
lat_atom_num = lat_atom_num + 1
r_atom_lat(:,lat_atom_num) = v
end if
end do
end do
end do
else
!If we are working with elements we have to use more complex code
!Initialize finite element
ele(:,:) = (esize-1) * cubic_cell(:,:)
!Make a 3 dimensional array of lattice points. This array is indexed by the integer lattice position.
!A value of true means that the coordinate is a lattice point which is within the box_bd
allocate(lat_points(bd_in_lat(2)-bd_in_lat(1)+10,bd_in_lat(4)-bd_in_lat(3)+10,bd_in_lat(6)-bd_in_lat(5)+10))
lat_points(:,:,:) = .false.
do iz = bd_in_lat(5)-5, bd_in_lat(6)+5
do iy = bd_in_lat(3)-5, bd_in_lat(4)+5
do ix = bd_in_lat(1)-5, bd_in_lat(2)+5
v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /)
!Transform point back to real space for easier checking
v = matmul(orient, matmul(transform_matrix,v))
!If within the boundaries
if(in_block_bd(v, box_bd)) then
lat_points(ix-bd_in_lat(1)+5,iy-bd_in_lat(3)+5,iz-bd_in_lat(5) + 5) = .true.
end if
end do
end do
end do
!Now we redefine bd_in_lat The first 3 indices contains limits for the lat_points array
bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10
bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10
bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10
!Figure out where the starting point is. This is the first piont which fully contains the finite element
outerloop: do iz = 1, bd_in_array(3)
do iy = 1, bd_in_array(2)
do ix = 1, bd_in_array(1)
node_in_bd(8) = .false.
do inod = 1, 8
vlat = ele(:,inod) + (/ ix, iy, iz /)
!Check to see if the node resides at a position containing a lattice point within the box
if(any(vlat > shape(lat_points))) then
continue
else if(lat_points(vlat(1),vlat(2),vlat(3))) then
node_in_bd(inod) = .true.
end if
end do
if(all(node_in_bd)) then
rzero = (/ ix, iy, iz /)
exit outerloop
end if
end do
end do
end do outerloop
!Now build the finite element region
lat_ele_num = 0
lat_atom_num = 0
allocate(r_lat(3,8,numlatpoints/esize))
!Redefined the second 3 indices as the number of elements that fit in the bounds
do i = 1, 3
bd_in_array(3+i) = bd_in_array(i)/esize
end do
!Now start the element at rzero
do inod=1, 8
ele(:,inod) = ele(:,inod) + rzero
end do
do iz = -bd_in_array(6), bd_in_array(6)
do iy = -bd_in_array(5), bd_in_array(5)
do ix = -bd_in_array(4), bd_in_array(4)
node_in_bd(:) = .false.
temp_nodes(:,:,:) = 0.0_dp
temp_lat(:,:) = 0
do inod = 1, 8
vlat= ele(:,inod) + (/ ix*(esize), iy*(esize), iz*(esize) /)
!Transform point back to real space for easier checking
! v = matmul(orient, matmul(transform_matrix,v))
do i = 1,3
v(i) = real(vlat(i) + bd_in_lat(2*i-1) - 5)
end do
temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v))
temp_lat(:,inod) = vlat
!Check to see if the lattice point values are greater then the array limits
if(any(vlat > shape(lat_points)).or.any(vlat < 1)) then
continue
!If within array boundaries check to see if it is a lattice point
else if(lat_points(vlat(1),vlat(2),vlat(3))) then
node_in_bd(inod) = .true.
end if
end do
if(all(node_in_bd)) then
lat_ele_num = lat_ele_num+1
r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
!Now set all the lattice points contained within an element to false
do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))
do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:))
lat_points(m,n,o) = .false.
end do
end do
end do
end if
end do
end do
end do
!Now figure out how many lattice points could not be contained in elements
allocate(r_atom_lat(3,count(lat_points)))
lat_atom_num = 0
do iz = 1, bd_in_array(3)
do iy = 1, bd_in_array(2)
do ix = 1, bd_in_array(1)
!If this point is a lattice point then save the lattice point as an atom
if (lat_points(ix,iy,iz)) then
v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /)
do i = 1,3
v(i) = v(i) + real(bd_in_lat(2*i-1) - 5)
end do
!Transform point back to real space
v = matmul(orient, matmul(transform_matrix,v))
lat_atom_num = lat_atom_num + 1
r_atom_lat(:,lat_atom_num) = v
end if
end do
end do
end do
end if
end subroutine build_with_rhomb
subroutine error_message(errorid)
integer, intent(in) :: errorid
select case(errorid)
case(1)
STOP "Name is missing."
case(2)
print *, "Element_type is missing"
case(3)
print *, "Lattice_parameter is missing"
case(4)
print *, "Lattice parameter is not in float form"
case(5)
print *, "Esize is missing"
case(6)
print *, "Esize is not in integer form"
case(7)
print *, "Cx(1) should equal Cx(2) in plane centroid finding algorithm. Please double check implementation"
end select
STOP 3
end subroutine error_message
end module mode_create