module opt_group

    !This module contains all code associated with dislocations

    use atoms
    use parameters
    use elements
    use subroutines
    use box 
    use sorts

    implicit none

    integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize, insert_type, &
               insert_site, num_species, group_atom_types, exclude_num, exclude_types(10)
    character(len=15) :: type, gshape!Type indicates what element type is selected and shape is the group shape
    character(len=2) :: species_type(10)
    real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness, insert_conc, &
    insert_lattice, s_fractions(10), max_height

    logical :: displace, delete, max_remesh, refine, group_nodes, flip, efill, refinefill, alloy

    integer, allocatable :: element_index(:), atom_index(:) 

    public 
    contains

    subroutine group(arg_pos)
        !Main calling function for the group option
        integer, intent(inout) :: arg_pos

        print *, '------------------------------------------------------------'
        print *, 'Option Group'
        print *, '------------------------------------------------------------'

        group_ele_num = 0
        group_atom_num = 0
        remesh_size=0
        insert_type = 0
        random_num=0
        group_type=0
        group_atom_types=0
        notsize=0
        displace=.false.
        delete=.false.
        max_remesh=.false.
        refine = .false.
        flip = .false.
        group_nodes=.false.
        num_species = 0
        exclude_num = 0

        if(allocated(element_index)) deallocate(element_index)
        if(allocated(atom_index)) deallocate(atom_index)

        call parse_group(arg_pos)


        !Now call the transformation functions for the group
         if(refine) then 
            call get_group
            call refine_group
        end if

        if(refinefill) then 
            call get_group
            call refinefill_group
        end if

        if(delete)then 
            call get_group
            call delete_group
        end if

        if(remesh_size > 0)then 
            print *, "Remesh command has been dropped"
            stop 1
            call get_group
!            call remesh_group
        end if

        if(group_type > 0) then 
            call get_group
            call change_group_type
        end if

        if(displace)then 
            call get_group
            call displace_group
        end if
        if(insert_type > 0) then 
            !print *, "Insert command has been dropped"
            !stop 1
            call get_group
            call insert_group
        end if

        if(num_species > 0) then 
            call get_group
            call alloy_group
        end if


    end subroutine group

    subroutine parse_group(arg_pos)
        !Parse the group command
        integer, intent(inout) :: arg_pos

        integer :: i, j, arglen, in_num
        character(len=100) :: textholder, type_spec
        real(kind=dp)  H, mass

        !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', 'all')
            continue
        case default
            print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", &
                     "Please select from atoms, elements, or all."
        end select

        arg_pos = arg_pos + 1
        call get_command_argument(arg_pos, gshape, 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(gshape)))
        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 group command"
                call parse_pos(int((i+1)/2), textholder, block_bd(i))   
            end do 

        case('wedge')
            arg_pos  = arg_pos + 1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing normal dim in group wedge command"
            read(textholder,*) dim1
            
            arg_pos  = arg_pos + 1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing normal dim in group wedge command"
            read(textholder,*) dim2

            do i = 1, 3
                arg_pos = arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if (arglen==0) STOP "Missing centroid in group wedge command"
                call parse_pos(i, textholder, centroid(i))           
            end do

            arg_pos  = arg_pos + 1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing base width in group wedge command"
            read(textholder,*)  bwidth

            !Calculate the vertex positions
            vertices(:,1) = centroid
            vertices(dim2,1) = 0.0_dp
            do i = 1, 3
                if (i == dim1) then 
                    if (bwidth > 0) then 
                        vertices(i,2) = box_bd(2*i)
                        vertices(i,3) = box_bd(2*i)
                    else if (bwidth < 0) then
                        vertices(i,2) = box_bd(2*i-1)
                        vertices(i,3) = box_bd(2*i-1)
                    else
                        print *, "bwidth cannot be 0 in wedge gshaped group"
                        stop 3
                    end if
                else if (i == dim2) then 
                    vertices(i,2) = 0.0_dp
                    vertices(i,3) = 0.0_dp
                else 
                    vertices(i,2) = centroid(i) + bwidth
                    vertices(i,3) = centroid(i) - bwidth
                end if
            end do 

        case('wedge_cut')
            arg_pos  = arg_pos + 1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing normal dim in group wedge_cut command"
            read(textholder,*) dim1
            
            arg_pos  = arg_pos + 1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing normal dim in group wedge_cut command"
            read(textholder,*) dim2

            do i = 1, 3
                arg_pos = arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if (arglen==0) STOP "Missing centroid in group wedge_cut command"
                call parse_pos(i, textholder, centroid(i))           
            end do

            arg_pos  = arg_pos + 1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing base width in group wedge_cut command"
            read(textholder,*)  bwidth

            arg_pos  = arg_pos + 1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing base width in group wedge_cut command"
            read(textholder,*)  max_height 

            !Calculate the vertex positions
            vertices(:,1) = centroid
            vertices(dim2,1) = 0.0_dp
            do i = 1, 3
                if (i == dim1) then 
                    if (bwidth > 0) then 
                        vertices(i,2) = box_bd(2*i)
                        vertices(i,3) = box_bd(2*i)
                    else if (bwidth < 0) then
                        vertices(i,2) = box_bd(2*i-1)
                        vertices(i,3) = box_bd(2*i-1)
                    else
                        print *, "bwidth cannot be 0 in wedge gshaped group"
                        stop 3
                    end if
                else if (i == dim2) then 
                    vertices(i,2) = 0.0_dp
                    vertices(i,3) = 0.0_dp
                else 
                    vertices(i,2) = centroid(i) + bwidth
                    vertices(i,3) = centroid(i) - bwidth
                end if
            end do 

        case('notch')
            arg_pos  = arg_pos + 1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing normal dim in group notch command"
            read(textholder,*) dim1
            
            arg_pos  = arg_pos + 1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing normal dim in group notch command"
            read(textholder,*) dim2

            do i = 1, 3
                arg_pos = arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if (arglen==0) STOP "Missing centroid in group notch command"
                call parse_pos(i, textholder, centroid(i))           
            end do

            arg_pos  = arg_pos + 1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing base width in group notch command"
            read(textholder,*)  bwidth

            arg_pos  = arg_pos + 1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing tip radius in group notch command"
            read(textholder,*)  radius

            !Calculate the vertex positions
            vertices(:,1) = centroid
            vertices(dim2,1) = 0.0_dp
            do i = 1, 3
                if (i == dim1) then 
                    if (bwidth > 0) then 
                        vertices(i,2) = box_bd(2*i)
                        vertices(i,3) = box_bd(2*i)
                        H= (box_bd(2*i) - centroid(i)) !Calculate the height of the wedge
                    else if (bwidth < 0) then
                        vertices(i,2) = box_bd(2*i-1)
                        vertices(i,3) = box_bd(2*i-1)
                        H = (centroid(i) - box_bd(2*i-1))
                    else
                        print *, "bwidth cannot be 0 in wedge shaped group"
                        stop 3
                    end if
                else if (i == dim2) then 
                    vertices(i,2) = 0.0_dp
                    vertices(i,3) = 0.0_dp
                else 
                    vertices(i,2) = centroid(i) + 0.5_dp*bwidth
                    vertices(i,3) = centroid(i) - 0.5_dp*bwidth
                end if
            end do 
             
            !Now update the centroid so that the desire tip diameter matches the wedge with
            if (bwidth > 0) then 
                centroid(dim1) = centroid(dim1) + 2*radius*(H)/bwidth
            end if
        !Read the ID type shape for group 
        case('id')
            arg_pos = arg_pos + 1
            !For this type we have to call different options if type is atoms, elements, or all. all is the most complex.
            select case(trim(adjustl(type)))
            case('atoms')
                !Read number of ids
                call get_command_argument(arg_pos, textholder, arglen)
                if(arglen == 0) then 
                    print *, "Missing number of input atoms for group id"
                end if
                read(textholder,*) in_num
                
                !allocate arrays
                allocate(atom_index(in_num))
                !Read ids
                do i = 1, in_num
                    arg_pos = arg_pos + 1
                    call get_command_argument(arg_pos, textholder, arglen)
                    if(arglen == 0) then 
                        print *, "Missing atom id in group atom id"
                        stop 3
                    end if 
                    read(textholder,*) atom_index(i)
                end do 
                group_atom_num = in_num

            case('elements')
                !Read number of ids
                 call get_command_argument(arg_pos, textholder, arglen)
                if(arglen == 0) then 
                    print *, "Missing number of input elements for group id"
                end if
                read(textholder, *) in_num

                !allocate arrays
                allocate(element_index(in_num))
                !Read ids
                do i = 1, in_num
                    arg_pos = arg_pos + 1
                    call get_command_argument(arg_pos, textholder, arglen)
                    if(arglen == 0) then 
                        print *, "Missing element id in group element id"
                        stop 3
                    end if 
                    read(textholder,*) element_index(i)
                end do 
                group_ele_num = in_num

            case('all')
                !We repeat this code twice, once for the atoms and once for the elements
                allocate(element_index(1024),atom_index(1024))
                do i = 1, 2
                    !Read the first id type (either atoms or elements)
                    call get_command_argument(arg_pos, type_spec, arglen)
                    if (arglen == 0) then
                        print *, "Missing type specifier in group id command"
                        stop 3
                    end if
                    
                    !Now read the first in_num
                    arg_pos = arg_pos + 1
                    call get_command_argument(arg_pos, textholder, arglen)
                    if(arglen == 0) then 
                        print *, "Missing input number in group_id"
                        stop 3 
                    end if
                    read(textholder, *) in_num

                    select case(trim(adjustl(type_spec)))
                    case('atoms','atom')
                        if (group_atom_num > 0) then
                            print *, "Atoms specifier used more than once in group id command with type all, either use type ", &
                                     "atoms or include elements identifier"
                            stop 3

                            do j = 1, in_num
                                arg_pos = arg_pos + 1
                                call get_command_argument(arg_pos, textholder, arglen)
                                if (arglen == 0) then 
                                    print *, "Missing atom id in group atom id"
                                    stop 3
                                end if
                                read(textholder, *) atom_index(j)
                            end do
                            group_atom_num = in_num
                        end if

                    case('elements','element')
                        if (group_ele_num > 0) then
                            print *, "Elements specifier used more than once in group id command with type all, either use type ",&
                                     "elements or include atoms identifier"
                            stop 3

                            do j = 1, in_num
                                arg_pos = arg_pos + 1
                                call get_command_argument(arg_pos, textholder, arglen)
                                if (arglen == 0) then 
                                    print *, "Missing element id in group element id"
                                    stop 3
                                end if
                                read(textholder, *) element_index(j)
                            end do
                            group_ele_num = in_num
                        end if

                    end select
                    if(i ==1) arg_pos = arg_pos + 1
                end do
            end select

        case('sphere')
            !First extract the sphere centroid
            do i = 1, 3
                arg_pos = arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if (arglen==0) STOP "Missing sphere centroid in group command"
                call parse_pos(i, textholder, centroid(i))   
            end do
            !Now get the tip radius
            arg_pos=arg_pos+1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing sphere radius in group command"
            read(textholder, *) radius

        case('shell')

            !First extract the shell centroid
            do i = 1, 3
                arg_pos = arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if (arglen==0) STOP "Missing shell centroid in group command"
                call parse_pos(i, textholder, centroid(i))   
            end do

            !Now get the radius
            arg_pos=arg_pos+1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing shell radius in group command"
            read(textholder, *) radius

            !Now get the shell thickness
            arg_pos=arg_pos+1
            call get_command_argument(arg_pos, textholder, arglen)
            if (arglen==0) STOP "Missing shell thickness in group command"
            read(textholder, *) shell_thickness

        case('all')
            !Do nothing if the shape is all
            continue

        case default
            print *, "Group shape ", trim(adjustl(gshape)), " 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('shift')
                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('refine')
                refine=.true.
            case('refinefill')
                refinefill=.true.
            case('selecttype')
                arg_pos = arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if (arglen==0) stop "Missing select type"
                call atommass(textholder, mass)
                call add_atom_type(mass, group_atom_types)
            case('excludetypes')
                arg_pos=arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if(arglen==0) stop "Missing number of atoms to exclude"
                read(textholder, *) exclude_num
                do i=1,exclude_num
                    arg_pos = arg_pos + 1
                    call get_command_argument(arg_pos, textholder, arglen)
                    if(arglen==0) stop "Missing exclude atom types"
                    call atommass(textholder, mass)
                    call add_atom_type(mass, exclude_types(i))
                end do
            case('remesh')
                arg_pos = arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if (arglen==0) stop "Missing remesh element size in group command"
                read(textholder, *) remesh_size
            case('max')
                max_remesh =.true.
            case('delete')
                delete=.true.
            case('nodes')
                group_nodes=.true.
            case('random')
                arg_pos = arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if (arglen==0) stop "Missing number of random atoms in group command"
                read(textholder, *) random_num
            case('flip')
                flip=.true.
            case('efill')
                efill=.true.
            case('type')
                arg_pos = arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if (arglen==0) stop "Missing atom type for group"
                call atommass(textholder, mass)
                call add_atom_type(mass, group_type)
            case('notsize')
                arg_pos = arg_pos + 1
                call get_command_argument(arg_pos, textholder, arglen)
                if(arglen ==0) stop "Missing notsize size"
                read(textholder, *) notsize
                print *, "Ignoring elements with size ", notsize
            case('alloy')
                alloy=.true.
                arg_pos=arg_pos+1
                call get_command_argument(arg_pos, textholder, arglen)
                print *, textholder
                if(arglen == 0) stop "Missing alloy num"
                read(textholder, *) num_species
                do i = 1, num_species
                    arg_pos=arg_pos + 1
                    call get_command_argument(arg_pos, textholder, arglen)
                    print *, textholder
                    read(textholder, *) species_type(i)
                    arg_pos=arg_pos + 1
                    call get_command_argument(arg_pos, textholder, arglen)
                    print *, textholder
                    read(textholder, *) s_fractions(i)
                end do
            case('insert')
                arg_pos=arg_pos+1
                call get_command_argument(arg_pos, textholder, arglen)
                if (arglen==0) stop "Missing element type for insert command"
                call atommass(textholder, mass)
                call add_atom_type(mass, insert_type)
                arg_pos=arg_pos + 1 
                call get_command_argument(arg_pos, textholder, arglen)
                select case(trim(adjustl(textholder)))
                case('tetra')
                    insert_site=1
                case('octa')
                    insert_site=2
                case('mixed')
                    insert_site=3
                case default
                    print *, "site value must be tetra, octa, or mixed and not ", trim(adjustl(textholder))
                    stop 3
                end select
                arg_pos=arg_pos+1
                call get_command_argument(arg_pos, textholder, arglen)
                if(arglen ==0) stop "missing lattice_type in insert command"
                read(textholder, *) insert_lattice
                arg_pos=arg_pos+1
                call get_command_argument(arg_pos, textholder, arglen)
                if(arglen ==0) stop "missing concentration in insert command"
                read(textholder, *) insert_conc
            case default
                !if it isn't an available option to opt_group 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, j, inod, ibasis, temp, node_in_out(max_ng_node)
        integer, allocatable :: resize_array(:) 
        real(kind=dp) :: r_center(3), rand

        select case(trim(adjustl(gshape)))
        case('block')
            print *, "Group has block shape with boundaries: ", block_bd
        case ('wedge')
            print *, "Group has wedge shape with dim1", dim1, "and dim2", dim2, "and vertices ", vertices
        case ('notch')
            print *, "Group has notch shape with dim1", dim1, "and dim2", dim2, " tip radius ", radius,  "and vertices ", &
                      vertices
        case('id')
            print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms."
            return
        case('sphere')
            print *, "Group has sphere shape with centroid ", centroid(:), " and radius ", radius
        case('all')
            print *, "Group has all of type ", type
        end select

        !Reset group if needed
        if(allocated(element_index)) deallocate(element_index,atom_index)

        group_ele_num = 0
        group_atom_num = 0
        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', 'all')
            do i = 1, ele_num
                if(in_group_ele(size_ele(i), lat_ele(i), r_node(:,:,:,i))) 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

            if(random_num > 0) then 
                !If we have the random option enabled then we select random_num number of elements from the group and overwrite 
                !the group with those elements
                do i = 1, random_num
                    call random_number(rand)
                    j = i + floor((group_ele_num+1-i)*rand)
                    temp = element_index(j)
                    element_index(j) = element_index(i)
                    element_index(i) = temp
                end do

                group_ele_num = random_num
            end if

        end select
                !Check the type to see if we need to find the atoms within the group
        select case(trim(adjustl(type)))
        case('atoms','all')
            do i = 1, atom_num
                if(in_group(type_atom(i), r_atom(:,i)).neqv.flip) 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
            if(random_num > 0) then 
                !If we have the random option enabled then we select random_num number of atom from the group and overwrite 
                !the group with those atom
                do i = 1, random_num
                    call random_number(rand)
                    j = i + floor((group_atom_num+1-i)*rand)
                    temp = atom_index(j)
                    atom_index(j) = atom_index(i)
                    atom_index(i) = temp
                end do

                group_atom_num = random_num
            end if
        end select
        
        print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms."
    end subroutine get_group 

    subroutine displace_group
        !This subroutine applies a displacement to elements/atoms in the groups

        integer :: i, inod, ibasis

        print *, "Elements/atoms in group displaced by ", disp_vec

        !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

        !If we don't include the wrap command then we have to increase the size of the box
        if (.not.(wrap_flag)) then 
            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

    subroutine refine_group
        !This command is used to refine the group to full atomistics

        integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num
        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), data_interp(10, max_basisnum*max_esize**3),  &
                         vel_interp(3, max_basisnum*max_esize**3) 
        real(kind=dp), allocatable :: vel_tmp(:,:)
        logical :: refine_dat


        if (allocated(force_node)) then
            refine_dat=.true.
        else 
            refine_dat=.false.
        end if

        print *, size(data_interp)
        !Refining to atoms
        if(group_ele_num > 0) then 
            orig_atom_num = atom_num
            !Estimate number of atoms we are adding, this doesn't have to be exact
            add_atom_num = group_ele_num*max_basisnum*(max_esize)**3
            call grow_ele_arrays(0,add_atom_num)
            do i = 1, group_ele_num
                ie = element_index(i)
                !Get the interpolated atom positions
                if (refine_dat) then 
                    call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp, &
                                           energy_node(:,:,i), force_node(:,:,:,i), virial_node(:,:,:,i), data_interp) 
                else
                    call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
                end if
                if(allocated(vel_atom)) then 
                    call interpolate_vel(type_ele(i), size_ele(i), lat_ele(i), vel_node(:,:,:,i), vel_interp)
                end if
                !Loop over all interpolated atoms and add them to the system, we apply periodic boundaries 
                !here as well to make sure they are in the box
                do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3
                    call apply_periodic(r_interp(:,j))
                    call add_atom(0,type_interp(j), r_interp(:,j))
                    if(refine_dat) then 
                        call add_atom_data(atom_num, data_interp(1, j), data_interp(2:4, j), data_interp(5:10, j))
                    end if
                    if(allocated(vel_atom)) then
                        if(size(vel_atom, 2) < size(type_atom)) then 
                            allocate(vel_tmp(3, size(type_atom)))
                            vel_tmp=0.0_dp
                            vel_tmp(:, 1:size(vel_atom,2)) = vel_atom
                            call move_alloc(vel_tmp, vel_atom)
                        end if
                        vel_atom(:, atom_num) = vel_interp(:,j)
                    end if
                end do 
            end do
            !Once all atoms are added we delete all of the elements
            call delete_elements(group_ele_num, element_index)
            print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms."

        end if                

    end subroutine refine_group

    subroutine refinefill_group
        !This command is used to refine the group to full atomistics

        integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num, m, n, o, esize, &
                   ele(3,8), new_ele_num, ibasis, inod, vlat(3), nump_ele, added_points
        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum,max_ng_node), ravg(3), ratom(3,max_basisnum)
        logical :: lat_points(max_esize, max_esize, max_esize)

        !Refining to atoms
        if(group_ele_num > 0) then 
            orig_atom_num = atom_num
            new_ele_num = 0
            !Estimate number of atoms we are adding, this doesn't have to be exact
            add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3
            call grow_ele_arrays(0,add_atom_num)
            do i = 1, group_ele_num
                ie = element_index(i)
                !Find all possible elements that we can make while making sure they aren't in the group
                lat_points(1:size_ele(ie),1:size_ele(ie),1:size_ele(ie)) = .true.

                !Now calculate the number of elemets which are available for remeshing
                nump_ele = size_ele(ie)**3
                do o =1, size_ele(ie)
                    do n = 1, size_ele(ie)
                        do m =1, size_ele(ie)
                            call get_interp_pos(m,n,o,i,rfill(:,:,1))

                            ravg(:) = 0
                            do ibasis = 1, basisnum(lat_ele(ie))
                                ravg = ravg + rfill(:,ibasis, 1)/basisnum(lat_ele(ie))
                            end do 

                            if( in_group(0,ravg)) then 
                                nump_ele = nump_ele - 1
                            end if
                        end do
                    end do
                end do

                !Now start the remeshing loop for the element
                esize = size_ele(ie) - 2
                added_points=0
                do while(esize > min_efillsize)
                    if(nump_ele < min_efillsize**3) then 
                        exit
                    else if (nump_ele < esize**3) then 
                        esize = esize - 2
                    else
                        ele = cubic_cell*(esize-1)
                        do o = 1, size_ele(ie) - esize
                            do n = 1, size_ele(ie) - esize
                                latloop:do m = 1, size_ele(ie) - esize
                                    do inod = 1, ng_node(lat_ele(ie))
                                        vlat = ele(:,inod) + (/ m, n, o /)
                                        if (.not.lat_points(vlat(1), vlat(2),vlat(3))) cycle latloop
                                        call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod))
                                    end do

                                    !Check to make sure all lattice points exist for the current element
                                    if(any(.not.lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1))) cycle latloop 
                                    if (.not.in_group_ele(esize, lat_ele(ie), rfill)) then 
                                        nump_ele=nump_ele - esize**3
                                        lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) = .false.
                                        call add_element(0,type_ele(ie), esize, lat_ele(ie),  rfill)
                                        new_ele_num = new_ele_num + 1
                                        added_points = added_points + esize**3
                                    end if

                                end do latloop
                            end do
                        end do
                        esize=esize-2
                    end if
                end do
                !Now add the leftover lattice points as atoms
                do o = 1, size_ele(ie)
                    do n = 1, size_ele(ie)
                        do m = 1, size_ele(ie)
                            if(lat_points(m,n,o)) then 
                                call get_interp_pos(m,n,o, ie, ratom(:,:))
                                do ibasis = 1, basisnum(lat_ele(ie))
                                    call apply_periodic(ratom(:,ibasis))
                                    call add_atom(0, basis_type(ibasis,lat_ele(ie)), ratom(:,ibasis))
                                    added_points=added_points + 1
                                end do
                            end if
                        end do
                    end do
                end do

                if (added_points /= (size_ele(ie)**3)) then 

                    print *, "Element ", ie, " is refined incorrectly in refinefill"
                end if
            end do
            !Once all atoms are added we delete all of the elements
            call delete_elements(group_ele_num, element_index)
            print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms and ", new_ele_num, &
                     " elements."
        end if                

    end subroutine refinefill_group

!    subroutine remesh_group
!        !This command is used to remesh the group to a desired element size
!
!        integer :: i, j, k, ix, iy, iz, inod, ibasis, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num,   &
!                   current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3), new_ele, new_atom, &
!                   max_loops, working_esize, group_lat_num, lat_list(10), is, ilat, tot_dof
!
!        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3),  ori_inv(3,3),  r(3), &
!                         r_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8), group_bd(6)
!        logical, allocatable :: lat_points(:,:,:), new_sbox, new_lat
!        character(len=100) :: remesh_ele_type
!
!
!        !Right now we just hardcode only remeshing to elements
!        remesh_ele_type = 'fcc'
!
!        ! Determine which sub_boxes and lattices types are within in the group
!        group_sbox_num = 0
!        group_lat_num = 0
!        do i = 1, group_atom_num 
!            new_sbox=.true.
!            new_lat=.true.
!            do j = 1, group_sbox_num
!                if (sbox_list(j) == sbox_atom(atom_index(i))) then 
!                    new_sbox=.false.
!                    exit
!                end if
!            end do
!
!            if(new_sbox) then 
!                group_sbox_num = group_sbox_num + 1 
!                sbox_list(group_sbox_num) = sbox_atom(atom_index(i))
!            end if
!
!            do j = 1, group_lat_num
!                if (basis_type(1,lat_list(j)) == type_atom(atom_index(i))) then 
!                    new_lat = .false.
!                    exit
!                end if
!            end do 
!
!            if (new_lat) then 
!                group_lat_num = group_lat_num + 1 
!                do k = 1, lattice_types
!                    if (basis_type(1,k) == type_atom(atom_index(i)))  lat_list(group_lat_num) = k
!                end do
!            end if
!
!        end do 
!
!        do i = 1, group_ele_num
!            new_sbox=.true.
!            new_lat = .true.
!            do j = 1, group_sbox_num
!                if (sbox_list(j) == sbox_ele(element_index(i))) then 
!                    new_sbox = .false.
!                    exit
!                end if
!            end do
!            
!            if (new_sbox) then
!                 group_sbox_num = group_sbox_num + 1
!                sbox_list(group_sbox_num) = sbox_ele(element_index(i))
!            end if
!
!            do j = 1, group_lat_num
!                if (lat_list(group_lat_num) == lat_ele(element_index(i))) then 
!                    new_lat=.false.
!                    exit
!                end if
!            end do
!
!            if (new_lat) then
!                group_lat_num = group_lat_num + 1
!                lat_list(group_lat_num) = lat_ele(element_index(i))
!            end if
!        end do
!
!        new_atom = 0
!        new_ele=0
!        tot_dof=0
!        do is = 1, group_sbox_num 
!
!            orient = sub_box_ori(:, :, sbox_list(is))
!            call matrix_inverse(orient,3,ori_inv)
!
!            do ilat = 1, group_lat_num
!
!                !First calculate max position in lattice space to be able to allocate lat_points array, also sum the total number
!                !of degrees of freedom which are added
!                dof = 0 
!                do j = 1, 3
!                    group_bd(2*j) = -huge(1.0_dp)
!                    group_bd(2*j-1) = huge(1.0_dp)
!                end do
!                do i = 1, group_atom_num
!                    if ((type_atom(atom_index(i)) == basis_type(1,ilat)).and.(sbox_atom(atom_index(i)) == is)) then
!                        do j =1 ,3
!                            if (r_atom(j,atom_index(i)) > group_bd(2*j)) group_bd(2*j) = r_atom(j,atom_index(i))
!                            if (r_atom(j,atom_index(i)) < group_bd(2*j-1)) group_bd(2*j-1) = r_atom(j,atom_index(i))
!                        end do
!                        dof = dof + 1
!                    end if
!                end do 
!                
!                do i = 1, group_ele_num
!                    if ((lat_ele(element_index(i)) == ilat).and.(sbox_ele(element_index(i)) == is)) then 
!                        do inod =1, ng_node(ilat)
!                            do ibasis = 1, basisnum(ilat)
!                                do j = 1, 3
!                                    r =r_node(j,ibasis,inod,element_index(i)) 
!                                    if (r(j) > group_bd(2*j)) group_bd(2*j) = r(j)
!                                    if (r(j) < group_bd(2*j-1)) group_bd(2*j-1) = r(j)
!                                end do
!                            end do
!                        end do
!                        dof = dof + size_ele(element_index(i))**3
!                    end if
!                end do
!
!                !If for some reason there are no dof in this loop then cycle out
!                if(dof == 0) cycle
!                tot_dof = tot_dof+dof
!
!                group_in_lat = reshape((/ group_bd(1),group_bd(3),group_bd(5), &
!                                          group_bd(2),group_bd(3),group_bd(5), &
!                                          group_bd(2),group_bd(4),group_bd(5), &
!                                          group_bd(1),group_bd(4),group_bd(5), &
!                                          group_bd(1),group_bd(3),group_bd(6), &
!                                          group_bd(2),group_bd(3),group_bd(6), &
!                                          group_bd(2),group_bd(4),group_bd(6), &
!                                          group_bd(1),group_bd(4),group_bd(6) /), [3,8])
!
!                group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/lapa(ilat)))
!                do i = 1, 3
!                    bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:)))
!                    bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:)))
!                end do
!
!                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.
!                dof = 0
!
!                !Now place all group atoms and group interpolated atoms into lat_points
!                do i = 1, group_atom_num
!                    if (.not.((sbox_atom(atom_index(i)) == is).and.(basis_type(1,ilat) == type_atom(atom_index(i))))) cycle
!                    r = r_atom(:,atom_index(i))/lapa(ilat)
!                    r = matmul(fcc_inv,matmul(ori_inv,r))
!                    do j = 1, 3
!                        r_lat(j) = nint(r(j))
!                    end do 
!                    !Do a check to make sure the code is working and that lattice points aren't being written on top of each other. 
!                    !This is primarily a debugging statement
!                    if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then 
!                        stop "Multiple atoms share same position in lat point array, this shouldn't happen"
!                    else 
!                        lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true.
!                        dof = dof + 1
!                    end if
!                end do 
!
!                !Now place interpolated atoms within lat_points array
!                do i =1, group_ele_num
!                    if (.not.((sbox_ele(element_index(i)) == is).and.( lat_ele(element_index(i)) == ilat))) cycle
!                    ie = element_index(i)
!                    call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
!                    do j = 1, size_ele(ie)**3 * basisnum(lat_ele(ie))
!                        r = r_interp(:,j)/lapa(ilat)
!                        r = matmul(fcc_inv,matmul(ori_inv,r))
!                        do k = 1, 3
!                            r_lat(k) = nint(r(k))
!                        end do 
!                        !Do a check to make sure the code is working and that lattice points aren't being written on top of each  
!                        !other. This is primarily a debugging statement
!                        if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then 
!                            stop "Multiple atoms/interpolated atoms share same position in lat point array, this shouldn't happen"
!                        else 
!                            lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true.
!                            dof = dof + 1
!                        end if
!                    end do 
!                end do 
!
!                !Now run remeshing algorithm, not the most optimized or efficient but gets the job done
!                !Figure out new looping boundaries
!                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
!
!
!                if (max_remesh) then 
!                    max_loops = (remesh_size-3)/2
!                else 
!                    max_loops = 1
!                end if
!                do j = 1, max_loops
!                    working_esize = remesh_size - 2*(j-1)
!                    ele = (working_esize-1)*cubic_cell
!                    zloop: do iz = 1, bd_in_array(3)
!                        yloop: do iy = 1, bd_in_array(2)
!                            xloop: do ix = 1, bd_in_array(1)
!                                if (lat_points(ix, iy,iz)) then 
!                                    r_new_node(:,:,:) = 0.0_dp
!                                    
!                                    !Check to see if the element overshoots the bound
!                                    if (iz+working_esize-1 > bd_in_array(3)) then
!                                        exit zloop
!                                    else if (iy+working_esize-1 > bd_in_array(2)) then
!                                        cycle zloop
!                                    else if (ix+working_esize-1 > bd_in_array(1)) then 
!                                        cycle yloop
!                                    end if  
!
!                                    if (all(lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1))) then 
!                                        do inod = 1, ng_node(ilat)
!                                            vlat = ele(:,inod) + (/ix, iy, iz /)
!                                            do i = 1, 3
!                                                vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5
!                                            end do 
!                                           r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat)
!                                        end do 
!
!                                        lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1) = .false.
!                                        !Add the element, for the sbox we just set it to the same sbox that we get the orientation 
!                                        !from. In this case it is from the sbox of the first atom in the group.
!                                        new_ele = new_ele+1
!                                        call add_element(0,remesh_ele_type, working_esize, ilat, &
!                                                         sbox_atom(atom_index(1)),r_new_node)
!
!                                    end if
!                                end if
!                            end do xloop
!                        end do yloop
!                    end do zloop
!                end do 
!
!                !Now we have to add any leftover lattice points as atoms
!                do iz = 1, bd_in_array(3)
!                    do iy=1, bd_in_array(2)
!                        do ix = 1, bd_in_array(1)
!                            if(lat_points(ix,iy,iz)) then 
!                                vlat = (/ ix, iy, iz /)
!                                do i = 1, 3
!                                    vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5
!                                end do 
!                                lat_points(ix,iy,iz) = .false.
!                               r = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat)
!                               new_atom=new_atom+1
!                               call add_atom(0,basis_type(1,ilat), is, r)
!                            end if
!                        end do
!                    end do 
!                end do
!                deallocate(lat_points)
!            end do
!        end do
!
!        !Delete all elements and atoms to make space for new elements and atoms
!        call delete_atoms(group_atom_num, atom_index)
!        call delete_elements(group_ele_num, element_index)
!    
!        print *, tot_dof, " degrees of freedom in group"
!        print *, "remesh_group creates ", new_atom, " atoms and ", new_ele, " elements."
!    end subroutine remesh_group

    subroutine delete_group
        !This subroutine deletes all atoms/elements within a group

        print *, "Deleting group containing ", group_atom_num, " atoms and ", group_ele_num, " elements."

        !Delete atoms
        call delete_atoms(group_atom_num, atom_index)

        !Delete elements
        call delete_elements(group_ele_num, element_index)

    end subroutine delete_group

    subroutine change_group_type
        !This subroutine changes all atoms and nodes at atoms within a group to a specific type
        integer :: i, j, ltype,ibasis, inod, basis_type(10)
        
        print *, "Changing ", group_atom_num, " atoms and ", group_ele_num, " elements to atom type ", group_type

        !Change all atom group types
        do i = 1, group_atom_num
            j = atom_index(i)
            type_atom(j) = group_type
        end do

        !Map to a new lattice type for all element
        do i =1, group_ele_num
            j = element_index(i)
            ltype = lat_ele(j)
            do ibasis = 1, basisnum(ltype)
               basis_type(ibasis) = group_type 
            end do
            call lattice_map(basisnum(ltype), basis_type, ng_node(ltype), lapa(ltype), lat_ele(j))
        end do

    end subroutine change_group_type
    
    subroutine insert_group
        !This code inserts atoms into interstitial sites. This only works on atoms within the group due to the limitations with the
        !Coarse-graining methodology which doesn't allow for concentration fluctuations.
        real(kind=dp)  interstitial_sites(3,14), rand, rinsert(3)
        integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi
        integer, allocatable :: used_sites(:,:)

        !First save all of the displacement vectors from a lattice site to interstitial site
        !The first 6 are the octohedral sites and the next 8 are the tetrahedral sites
        interstitial_sites= reshape( (/ -0.5_dp, 0.0_dp, 0.0_dp, &
                                         0.5_dp, 0.0_dp, 0.0_dp, &
                                         0.0_dp,-0.5_dp, 0.0_dp, &
                                         0.0_dp, 0.5_dp, 0.0_dp, &
                                         0.0_dp, 0.0_dp,-0.5_dp, &
                                         0.0_dp, 0.0_dp, 0.5_dp, &
                                        -0.25_dp,-0.25_dp,-0.25_dp, &
                                        -0.25_dp,-0.25_dp, 0.25_dp, &
                                        -0.25_dp, 0.25_dp,-0.25_dp, &
                                        -0.25_dp, 0.25_dp, 0.25_dp, &
                                         0.25_dp,-0.25_dp,-0.25_dp, &
                                         0.25_dp,-0.25_dp, 0.25_dp, &
                                         0.25_dp, 0.25_dp,-0.25_dp, &
                                         0.25_dp, 0.25_dp, 0.25_dp /), &
                                         shape(interstitial_sites))

        !First we calculate the number of atoms needed based on the number of atoms in the group and the concentration
        interstitial_sites=interstitial_sites*insert_lattice


        !Now set up the random number generator for the desired interstitial type
        select case(insert_site)
        case(1)
            rlo=1
            rhi=6
        case(2)
            rlo=7 
            rhi = 14
        case(3)
            rlo=1
            rhi=14
        end select

        if ((insert_conc > 1).or.(insert_conc < 0)) stop "Insert concentration should be between 0 and 1"
        add_num = (insert_conc*group_atom_num)
        allocate(used_sites(2,add_num))

        print *, "Inserting ", add_num, " atoms as atom type ", insert_type

        !Now add the atoms
        i = 1
        addloop:do while ( i < add_num)
            call random_number(rand)
            rindex = int(1+rand*(group_atom_num-1))
            ia=atom_index(rindex)
            call random_number(rand)
            sindex = int(rlo+rand*(rhi-rlo))
            do j = 1, i
                if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop
            end do
            rinsert = r_atom(:,ia) + matmul(box_ori,interstitial_sites(:,sindex))
            call add_atom(0, insert_type, rinsert)
            used_sites(1,i) = ia
            used_sites(2,i) = sindex
            i = i + 1
        end do addloop
    end subroutine insert_group

    subroutine alloy_group
        !This subroutine randomizes the atom types to reach desired concentrations, this only operates on atoms
        integer :: i, j, ia, type_map(10), added_types(num_species)
        real(kind=dp) ::  rand, mass, old_fractions(num_species)
        character(len=2) :: sorted_species_type(10)
        logical :: species_mapped(10)


        print *, "Alloying group with desired fractions", s_fractions(1:num_species)
        old_fractions=s_fractions(1:num_species)
        !First we sort the fractions
        call quicksort(s_fractions(1:num_species))
        
        species_mapped = .false.
        do i = 1,num_species
            do j=1, num_species
                if (is_equal(s_fractions(i), old_fractions(j)) .and. .not. species_mapped(j)) then 
                    sorted_species_type(i)=species_type(j)
                    species_mapped(j) = .true.
                    exit
                end if
            end do
        end do


        !Now get the atom type maps for all the atoms and make the fractions a running sum
        do i = 1, num_species
            call atommass(sorted_species_type(i), mass)
            call add_atom_type(mass, type_map(i))
            if(i > 1) s_fractions(i) = s_fractions(i) + s_fractions(i-1)
        end do
        !Now randomize the atom types
        added_types = 0
        do i = 1, group_atom_num
            ia = atom_index(i)
            call random_number(rand)
            do j = 1, num_species
                if(rand < s_fractions(j)) then 
                    type_atom(ia) = type_map(j)
                    added_types(j) = added_types(j) +1 
                    exit
                end if
            end do
        end do

        print *, "Converted ", added_types(1:num_species), " atoms"

        return
    end subroutine alloy_group

    function in_group(atype, r)
        !This subroutine determines if a point is within the group boundaries
        integer, intent(in) :: atype
        real(kind=dp), intent(in) :: r(3)
        real(kind=dp) :: rnorm

        integer :: dim3, i
        logical :: in_group
        select case(trim(adjustl(gshape)))
        case('block')
            in_group=in_block_bd(r,block_bd)
        case('wedge')
            in_group = in_wedge_bd(r,vertices)
        case('wedge_cut')
            in_group = in_wedge_bd(r,vertices)
            in_group = in_group.and.((abs(r(dim1) - vertices(dim1,2)) < max_height))
        case('notch')
            do i = 1, 3
                if (.not.((dim1==i).or.(dim2==i))) dim3 = i
            end do 

            in_group = in_wedge_bd(r,vertices)
            !Do a check to make sure the wedge isn't used if it should be the tip radius
            if (bwidth>0) then 
                if (r(dim1) < centroid(dim1)) in_group = .false.
            end if

            rnorm = sqrt((r(dim1) - centroid(dim1))**2 + (r(dim3)-centroid(dim3))**2)
            in_group = in_group.or.(rnorm < radius)

        case('sphere') 
            rnorm = norm2(r(:) - centroid(:))
            if (rnorm <= radius) then 
                in_group = .true. 
            else
                in_group = .false.
            end if
        case('shell') 
            rnorm = norm2(r(:) - centroid(:))
            if ((rnorm >= radius).and.(rnorm<=(radius + shell_thickness))) then 
                in_group = .true. 
            else
                in_group = .false.
            end if
        case('all')
            in_group = .true.
        end select 

        if(group_atom_types> 0) then 
            in_group = in_group.and.(atype== group_atom_types)
        elseif(exclude_num > 0) then 
            if(in_group) then
                do i=1,exclude_num 
                    if (atype == exclude_types(i)) then 
                        in_group=.false.
                        exit
                    end if
                end do
            end if
        end if
    end function in_group

    function in_group_ele(esize, elat, rn)
        !This figures out if the elements are in the group boundaries
        real(kind=dp), intent(in) :: rn(3,max_basisnum, max_ng_node)
        integer, intent(in) :: esize, elat
        logical :: in_group_ele
        
        integer :: i, inod, ibasis, node_in_out(max_ng_node)
        real(kind=dp) :: r_center(3) 

        in_group_ele=.false.

        if(trim(adjustl(gshape)) == 'shell') then 
            node_in_out(:) = -1
            !First calculate whether each element node is within the shell region, inside the shell sphere, or outside the
            !shell region
            nodeloop:do inod = 1, ng_node(elat)
                r_center(:)=0.0_dp
                do ibasis = 1, basisnum(elat)
                    r_center(:)= r_center(:) + rn(:,ibasis,inod)/basisnum(elat) 
                end do

                if((in_group(0,rn(:, ibasis, inod)).neqv.flip).and.(size_ele(i)/=notsize)) then 
                    node_in_out(inod) = 2
                    exit nodeloop 
                end if
                
                gshape ='sphere'
                if((in_group(0,rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then 
                    node_in_out(inod) = 1
                else
                    node_in_out(inod) = 0
                end if
                gshape='shell'
            end do nodeloop

            !If any nodes are within the shell region, or if the shell region interescts an element then add it to the group
            if(any(node_in_out == 2).or.(any(node_in_out==1).and.(any(node_in_out==0)))) then 
                in_group_ele=.true.
                return
            end if

        else if(.not.(group_nodes)) then 
            r_center(:) = 0.0_dp
            do inod = 1, ng_node(elat)
                do ibasis = 1, basisnum(elat)
                    r_center = r_center + rn(:,ibasis,inod)/(basisnum(elat)*ng_node(elat))
                end do
            end do

            if ((in_group(0,r_center).neqv.flip).and.(esize/= notsize)) then 
                in_group_ele=.true.
                return
            end if 

        else if(group_nodes) then 
            r_center(:) = 0.0_dp
            do inod = 1, ng_node(elat)
                do ibasis = 1, basisnum(elat)
                     if ((in_group(0,rn(:,ibasis,inod)).neqv.flip).and.(esize/=notsize)) then 
                         in_group_ele=.true.
                         return
                    end if 
                end do
            end do
        end if

    end function in_group_ele

end module opt_group
