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),  basis_pos(3,10)
    integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
               esize_nums, esize_index(10)
    logical :: dup_flag, dim_flag, efill, crossb(3)

    real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
    integer, allocatable :: elat(:)
    public 
    contains

    subroutine create(arg_pos)
        ! Main subroutine which controls execution

        integer, intent(out) :: arg_pos

        integer :: i, ibasis, inod, ei, curr_esize
        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
        efill = .false.

        !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 and get the new maxlen
            select case(trim(adjustl(element_type)))
            case('fcc', '20fcc')
                box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
                maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2)
            case('bcc')
                box_lat_vert = int(matmul(bcc_inv, matmul(orient_inv, box_vert)))+1
                maxbd = maxval(matmul(orient,matmul(bcc_mat,box_lat_vert)),2)
            end select
                
            !Find the new maxlen
            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)
                box_len(i) = box_bd(2*i) - box_bd(2*i-1)
            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
            select case(trim(adjustl(element_type)))
            case('fcc', '20fcc')
                box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
            case('bcc')
                box_lat_vert = int(matmul(bcc_inv, matmul(orient_inv, box_vert)))+1
            end select

            !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,:,:))+10.0_dp**(-6.0_dp)
                box_bd(2*i-1) = minval(r_node_temp(i,:,:)) - 10.0_dp**(-6.0_dp)
            end do
            call add_element(0,element_type, esize, 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, 8, fcc_inv)
            case('bcc')
                call build_with_rhomb(box_lat_vert, bcc_mat, 8, bcc_inv)
            case('20fcc')
                call build_with_rhomb(box_lat_vert, fcc_mat, 20, fcc_inv)
            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(0,basis_type(ibasis, 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(0,element_type, elat(i), 1, r_node_temp)
                end do 
            end if
        end if

        !If any elements are fully outside the box then wrap them back in
        if (any(crossb)) then 
            call wrap_elements
        end if
        box_ori = orient
    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=100) :: orient_string
        character(len=2) :: btype
        logical :: isortho, isrighthanded, bool
        real(kind=dp) :: mass


        !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
        max_esize = 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

                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('crossb')
                do i = 1, 3
                    call get_command_argument(arg_pos, textholder)
                    read(textholder, *) crossb(i) 
                    arg_pos = arg_pos + 1
                end do
            case('basis')
                call get_command_argument(arg_pos, textholder)
                read(textholder, *) basisnum(1)
                arg_pos = arg_pos + 1
                max_basisnum=basisnum(1)
                do i = 1, max_basisnum
                    call get_command_argument(arg_pos, btype)
                    arg_pos = arg_pos+1
                    call atommass(btype, mass)
                    call add_atom_type(mass, basis_type(i,1))
                    do j = 1, 3
                        call get_command_argument(arg_pos, textholder)
                        read(textholder, *) basis_pos(j,i)
                        arg_pos=arg_pos+1
                    end do
                end do

            case('efill')
                efill=.true.
            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
        !Now normalize the orientation matrix
        orient =  matrix_normal(orient,3)
        !Check special periodicity relations
        select case(trim(adjustl(element_type)))
        case('fcc', '20fcc')
            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

        case('bcc')
            do i = 1, 3
                !Check if the direction is 111
                if((is_equal(abs(orient(i,1)),abs(orient(i,2)))).and.(is_equal(abs(orient(i,2)),abs(orient(i,3))))) then 
                    lattice_space(i) = 0.5_dp * lattice_space(i)
                end if
            end do
        end select 
        !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 atommass(name, mass)
            call add_atom_type(mass, 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, nn, transform_inverse)
        !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
        real(kind=dp), dimension(3,3), intent(in) :: transform_inverse !The inverse transform
        integer, intent(in) :: nn
        !Internal variables
        integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz,  numlatpoints, ele(3,nn), rzero(3), efill_size, &
                   vlat(3), temp_lat(3,nn), m, n, o, j, k, nump_ele, efill_temp_lat(3,nn), filzero(3), &
                   bd_ele_lat(6), efill_ele(3,nn), ebd(6), shift_vec(3), type_interp(max_basisnum*max_esize**3)
        real(kind=dp) :: v(3), temp_nodes(3,1,nn), r(3), centroid_bd(6), vreal(3), r_interp(3, max_basisnum*max_esize**3)
        logical, allocatable :: lat_points(:,:,:)
        logical :: node_in_bd(nn), add, lat_points_ele(esize,esize,esize), intersect_bd(3)

        !Do some value initialization
        max_esize = esize
        do i = 1,3
            centroid_bd(2*i) = -huge(1.0_dp)
            centroid_bd(2*i-1) = huge(1.0_dp)
        end do

        !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

        if(numlatpoints < 0) stop "number of lattice points is negative, this can occur if the model is too big"

        !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 
            select case(element_type)
            case('fcc','bcc')
                ele(:,:) = (esize-1) * cubic_cell(:,:)
            case('20fcc')
                ele(:,:) = (esize-1)*cube20
            end select
            
            !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(:) = .false.
                        do inod = 1, nn
                            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,nn,numlatpoints/esize), elat(numlatpoints/esize))

            select case(element_type)
            case('fcc','bcc')
                ele(:,:) = (esize-1) * cubic_cell(:,:)
            case('20fcc')
                ele(:,:) = (esize-1)*cube20
            end select
 
            !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, nn
                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, nn
                            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.
                            else if(any(crossb)) then 
                                vreal=0
                                do i = 1, 3
                                    if(crossb(i)) then 
                                        if(temp_nodes(i,1,inod) < box_bd(2*i-1)) then 
                                            vreal(i) = temp_nodes(i,1,inod)+box_len(i)
                                        else if(temp_nodes(i,1,inod) > box_bd(2*i)) then
                                            vreal(i) = temp_nodes(i,1,inod)-box_len(i)
                                        else 
                                            vreal(i) = temp_nodes(i,1,inod)
                                        end if
                                    else
                                        vreal(i) = temp_nodes(i,1,inod)
                                    end if
                                end do
                                v = matmul(transform_inverse, matmul(orient_inv, vreal))
                                do i = 1, 3
                                    vlat(i) = nint(v(i) - bd_in_lat(2*i-1)+5)
                                end do
                                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 if
                        end do

                        !If we are on the first round of element building then we can just add the element if all(node_in_bd) is
                        !true
                        if(all(node_in_bd)) then 
                            lat_ele_num = lat_ele_num+1
                            r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
                            elat(lat_ele_num) = esize
                            if(any(crossb)) then 
                                call interpolate_atoms('fcc', esize, 0, temp_nodes, type_interp, r_interp)
                                j= 0
                                 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,:))
                                            j=j+1
                                            do i = 1, 3
                                                if(crossb(i)) then 
                                                    if(r_interp(i,j) < box_bd(2*i-1)) then 
                                                        vreal(i) = r_interp(i,j)+box_len(i)
                                                    else if(r_interp(i,j) > box_bd(2*i))  then
                                                        vreal(i) = r_interp(i,j)-box_len(i)
                                                    else 
                                                        vreal(i) = r_interp(i,j)
                                                    end if
                                                else
                                                    vreal(i) = r_interp(i,j)
                                                end if
                                            end do
                                            v = matmul(transform_inverse, matmul(orient_inv, vreal))
                                            do i = 1, 3
                                                vlat(i) = nint(v(i) - bd_in_lat(2*i-1)+5)
                                            end do

                                            if(lat_points(vlat(1), vlat(2), vlat(3))) then  
                                                lat_points(vlat(1), vlat(2), vlat(3)) = .false.
                                            else
                                                print *, "Lat points should be true not false"
                                                stop 2
                                            end if
                                        end do 
                                    end do
                                end do                               
                            else
                                !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

                        !If any nodes are in the boundary and we want to efill then run the efill code
                        else if(any(node_in_bd).and.efill) then

                            !Pull out the section of the lat points array 
                            lat_points_ele(:,:,:)=.false.
                            do i = 1,3
                                if (minval(temp_lat(i,:)) <lim_zero) then 
                                    bd_ele_lat(2*i-1) = 1
                                else 
                                    bd_ele_lat(2*i-1) = minval(temp_lat(i,:))
                                end if
                                if(maxval(temp_lat(i,:))>size(lat_points,i)) then 
                                    bd_ele_lat(2*i) = size(temp_lat(i,:))
                                else 
                                    bd_ele_lat(2*i) = maxval(temp_lat(i,:))
                                end if
                            end do

                            lat_points_ele(1:(bd_ele_lat(2)-bd_ele_lat(1)+1),1:(bd_ele_lat(4)-bd_ele_lat(3)+1),&
                            1:(bd_ele_lat(6)-bd_ele_lat(5))+1)= lat_points(bd_ele_lat(1):bd_ele_lat(2), &
                                                                         bd_ele_lat(3):bd_ele_lat(4), &
                                                                         bd_ele_lat(5):bd_ele_lat(6)) 
                            !Now start looping through elements and try to fit as many as you can
                            efill_size = esize-2
                            i=1
                            j=1
                            k=1
                            nump_ele = count(lat_points_ele)
                            do i = 1, 3
                                filzero(i) = bd_ele_lat(2*i-1) -1
                            end do
                            do while(efill_size>min_efillsize)
                                !First check whether there are enough lattice points to house the current element size
                                efill_ele=cubic_cell*(efill_size-1)
                                if  (nump_ele < efill_size**3) then 
                                    efill_size = efill_size - 2
                                else 
                                    ze: do k = 1, (esize-efill_size)
                                        ye: do j = 1, (esize-efill_size)
                                            xe: do i = 1, (esize-efill_size)
                                                do inod = 1,nn
                                                    vlat = efill_ele(:,inod) + (/ i, j, k /)
                                                    if (.not.lat_points_ele(vlat(1),vlat(2),vlat(3))) cycle xe
                                                     do o = 1,3
                                                        v(o) = real(vlat(o) + filzero(o) + bd_in_lat(2*o-1) -5)
                                                    end do
                                                    temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v))
                                                    efill_temp_lat(:,inod) = vlat 
                                                end do

                                                do o = 1,3 
                                                    ebd(2*o-1) = minval(efill_temp_lat(o,:))
                                                    ebd(2*o) = maxval(efill_temp_lat(o,:))
                                                end do
                                                lat_ele_num = lat_ele_num+1
                                                r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
                                                elat(lat_ele_num) = efill_size
                                                nump_ele = nump_ele - efill_size**3
                                                !Now set all the lattice points contained within an element to false
                                                do o = ebd(5), ebd(6)
                                                    do n = ebd(3), ebd(4)
                                                        do m = ebd(1), ebd(2)
                                                            lat_points(m+filzero(1),n+filzero(2),o+filzero(3)) = .false.
                                                            lat_points_ele(m,n,o) = .false.
                                                        end do 
                                                    end do
                                                end do
                                            end do xe
                                        end do ye
                                    end do ze
                                    efill_size = efill_size-2
                                end if
                            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
