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/opt_group.f90

200 lines
7.4 KiB

module opt_group
!This module contains all code associated with dislocations
use parameters
use elements
use subroutines
use box
implicit none
integer :: group_ele_num, group_atom_num
character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape
real(kind=dp) :: block_bd(6), disp_vec(3)
logical :: displace, wrap
integer, allocatable :: element_index(:), atom_index(:)
public
contains
subroutine group(arg_pos)
!Main calling function for the group option
integer, intent(inout) :: arg_pos
group_ele_num = 0
group_atom_num = 0
if(allocated(element_index)) deallocate(element_index)
if(allocated(atom_index)) deallocate(atom_index)
call parse_group(arg_pos)
call get_group
!Now call the transformation functions for the group
if(displace) call displace_group
end subroutine group
subroutine parse_group(arg_pos)
!Parse the group command
integer, intent(inout) :: arg_pos
integer :: i,arglen
character(len=100) :: textholder
!Parse type and shape command
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, type, arglen)
if (arglen==0) STOP "Missing select_type in group command"
select case(trim(adjustl(type)))
case('atoms', 'elements', 'both')
continue
case default
print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", &
"Please select from atoms, nodes, or both."
end select
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, shape, arglen)
if (arglen==0) STOP "Missing group_shape in group command"
!Now parse the arguments required by the user selected shape
select case(trim(adjustl(shape)))
case('block')
do i= 1, 6
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing block boundary in dislgen command"
call parse_pos(int((i+1)/2), textholder, block_bd(i))
end do
case default
print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", &
"for accepted group shapes."
end select
!Now parse the additional options which may be present
do while(.true.)
if(arg_pos > command_argument_count()) exit
!Pull out the next argument which should either be a keyword or an option
arg_pos=arg_pos+1
call get_command_argument(arg_pos, textholder)
textholder=adjustl(textholder)
!Choose what to based on what the option string is
select case(trim(textholder))
case('displace')
displace = .true.
do i = 1,3
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) stop "Missing vector component for shift command"
read(textholder, *) disp_vec(i)
end do
case('wrap')
wrap = .true.
case default
!If it isn't an available option to opt_disl then we just exit
exit
end select
end do
end subroutine parse_group
subroutine get_group
!This subroutine finds all elements and/or atoms within the group boundaries
!specified by the user.
integer :: i, inod, ibasis
integer, allocatable :: resize_array(:)
real(kind=dp) :: r_center(3)
select case(trim(adjustl(shape)))
case('block')
!Allocate variables to arbitrary size
allocate(element_index(1024), atom_index(1024))
!Check the type to see whether we need to find the elements within the group
select case(trim(adjustl(type)))
case('elements', 'both')
do i = 1, ele_num
r_center(:) = 0.0_dp
do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
r_center = r_center + r_node(:,ibasis,inod,i)/(basisnum(lat_ele(i))*ng_node(lat_ele(i)))
end do
end do
if (in_block_bd(r_center, block_bd)) then
group_ele_num = group_ele_num + 1
if(group_ele_num > size(element_index)) then
allocate(resize_array(size(element_index) + 1024))
resize_array(1:group_ele_num-1) = element_index
resize_array(group_ele_num:) = 0
call move_alloc(resize_array, element_index)
end if
element_index(group_ele_num) = i
end if
end do
end select
!Check the type to see if we need to find the atoms within the group
select case(trim(adjustl(type)))
case('atoms','both')
do i = 1, atom_num
if(in_block_bd(r_atom(:,i),block_bd)) then
group_atom_num = group_atom_num + 1
if (group_atom_num > size(atom_index)) then
allocate(resize_array(size(atom_index) + 1024))
resize_array(1:group_atom_num -1) = atom_index
resize_array(group_atom_num:) = 0
call move_alloc(resize_array, atom_index)
end if
atom_index(group_atom_num) = i
end if
end do
end select
end select
end subroutine get_group
subroutine displace_group
!This subroutine applies a displacement to elements/atoms in the groups
integer :: i, inod, ibasis
!Displace atoms
do i = 1, group_atom_num
r_atom(:,atom_index(i)) = r_atom(:,atom_index(i)) + disp_vec
end do
!Displace elements
do i = 1, group_ele_num
do inod = 1, ng_node(lat_ele(element_index(i)))
do ibasis = 1, basisnum(lat_ele(element_index(i)))
r_node(:,ibasis,inod,element_index(i)) = r_node(:,ibasis,inod,element_index(i)) + disp_vec
end do
end do
end do
!Now either apply periodic boundaries if wrap command was passed or adjust box dimensions
!Now we check if we have to wrap the atoms, nodes are not wrapped. For elements the periodic
!boundary conditions are applied in the actual CAC codes
if(wrap) then
do i = 1, atom_num
call apply_periodic(r_atom(:,i))
end do
!If we don't include the wrap command then we have to increase the size of the box
else
do i = 1,3
if (disp_vec(i) < -lim_zero) then
box_bd(2*i-1) = box_bd(2*i-1) - disp_vec(i)
else if (disp_vec(i) > lim_zero) then
box_bd(2*i) = box_bd(2*i) + disp_vec(i)
end if
end do
end if
end subroutine displace_group
end module opt_group