CACmb/src/elements.f90

module elements
    
    !This module contains the elements data structures, structures needed for building regions
    !and operations that are done on elements
    use parameters
    use functions
    use subroutines
    use sorts
    use box
    implicit none

    !Data structures used to represent the CAC elements. Each index represents an element
    character(len=100), allocatable :: type_ele(:) !Element type 
    integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:), tag_ele(:) !Element size
    real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
    !Element result data structures
    real(kind=dp), allocatable :: force_node(:,:,:,:), virial_node(:,:,:,:), energy_node(:,:,:), vel_node(:,:,:,:)

    integer, save :: ele_num !Number of elements
    integer, save :: node_num !Total number of nodes
    integer, save :: node_atoms !Count of all basis atoms at nodes summed over all nodes

    !Data structure used to represent atoms 
    integer, allocatable ::  type_atom(:)!atom type 
    integer, allocatable :: sbox_atom(:), tag_atom(:)
    real(kind =dp),allocatable :: r_atom(:,:) !atom position
    integer :: atom_num=0 !Number of atoms
    !Atom result data structures information
    real(kind=8), allocatable :: force_atom(:,:), virial_atom(:,:), energy_atom(:), vel_atom(:,:)

    !Mapping atom type to provided name
    character(len=2), dimension(10) :: type_to_name
    integer :: atom_types = 0

    !Variables for creating elements based on primitive cells
    real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3), &
                     cube20(3,20)   
    integer :: cubic_faces(4,6)

    !Below are lattice type arrays which provide information on the general form of the elements.
    !We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased.
    integer :: lattice_types = 0
    integer :: max_ng_node, ng_node(10) !Max number of nodes per element and number of nodes per element for each lattice type
    integer :: max_esize=0 !Maximum number of atoms per side of element
    real(kind=dp) :: lapa(10), masses(10)

    !These variables contain information on the basis, for simplicities sake we limit 
    !the user to the definition of 10 lattice types with 10 basis atoms at each lattice point.
    !This can be easily increased with no change to efficiency
    integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type 
    integer :: basis_type(10,10)

    !Additional module level variables we need 
    logical :: wrap_flag

    !flags for data variables
    logical :: vflag
    public 
    contains

    subroutine lattice_init
        !This subroutine just intializes variables needed for building the different finite 
        !element types

        !First initialize the cubic cell
        cubic_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
                              1.0_dp, 0.0_dp, 0.0_dp, &
                              1.0_dp, 1.0_dp, 0.0_dp, &
                              0.0_dp, 1.0_dp, 0.0_dp, &
                              0.0_dp, 0.0_dp, 1.0_dp, &
                              1.0_dp, 0.0_dp, 1.0_dp, &
                              1.0_dp, 1.0_dp, 1.0_dp, &
                              0.0_dp, 1.0_dp, 1.0_dp /), &
                             shape(fcc_cell))

        !Now initialize the 20 node cube element
        cube20 = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
                             1.0_dp, 0.0_dp, 0.0_dp, &
                             1.0_dp, 1.0_dp, 0.0_dp, &
                             0.0_dp, 1.0_dp, 0.0_dp, &
                             0.0_dp, 0.0_dp, 1.0_dp, &
                             1.0_dp, 0.0_dp, 1.0_dp, &
                             1.0_dp, 1.0_dp, 1.0_dp, &
                             0.0_dp, 1.0_dp, 1.0_dp, &
                             0.5_dp, 0.0_dp, 0.0_dp, &
                             1.0_dp, 0.5_dp, 0.0_dp, &
                             0.5_dp, 1.0_dp, 0.0_dp, &
                             0.0_dp, 0.5_dp, 0.0_dp, &
                             0.0_dp, 0.0_dp, 0.5_dp, &
                             1.0_dp, 0.0_dp, 0.5_dp, &
                             1.0_dp, 1.0_dp, 0.5_dp, &
                             0.0_dp, 1.0_dp, 0.5_dp, &
                             0.5_dp, 0.0_dp, 1.0_dp, &
                             1.0_dp, 0.5_dp, 1.0_dp, &
                             0.5_dp, 1.0_dp, 1.0_dp, &
                             0.0_dp, 0.5_dp, 1.0_dp /), &
                             shape(cube20))

        !Now we create a list containing the list of vertices needed to describe the 6 cube faces
        cubic_faces(:,1) = (/ 1, 2, 3, 4 /) 
        cubic_faces(:,2) = (/ 1, 2, 6, 5 /)
        cubic_faces(:,3) = (/ 2, 3, 7, 6 /)
        cubic_faces(:,4) = (/ 3, 4, 8, 7 /)
        cubic_faces(:,5) = (/ 1, 4, 8, 5 /)
        cubic_faces(:,6) = (/ 5, 6, 7, 8 /)

        !!Now initialize the fcc primitive cell
        fcc_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
                              0.5_dp, 0.5_dp, 0.0_dp, &
                             0.5_dp, 1.0_dp, 0.5_dp, &
                             0.0_dp, 0.5_dp, 0.5_dp, &
                             0.5_dp, 0.0_dp, 0.5_dp, &
                             1.0_dp, 0.5_dp, 0.5_dp, &
                             1.0_dp, 1.0_dp, 1.0_dp, &
                             0.5_dp, 0.5_dp, 1.0_dp /), &
                             shape(fcc_cell))

        fcc_mat = reshape((/ 0.5_dp, 0.5_dp, 0.0_dp, &
                             0.0_dp, 0.5_dp, 0.5_dp, &
                             0.5_dp, 0.0_dp, 0.5_dp /), &
                             shape(fcc_mat))

        !Initialize the bcc primitive cell
        bcc_cell = reshape((/ 0.0_dp,  0.0_dp,  0.0_dp, &
                              0.5_dp, -0.5_dp,  0.5_dp, &
                              1.0_dp,  0.0_dp,  1.0_dp, &
                              0.5_dp,  0.5_dp,  0.5_dp, &
                             -0.5_dp,  0.5_dp,  0.5_dp, &
                              0.0_dp,  0.0_dp,  1.0_dp, &
                              0.5_dp,  0.5_dp,  1.5_dp, &
                              0.0_dp,  1.0_dp,  1.0_dp /),  &
                              shape(bcc_cell))

        bcc_mat = reshape((/ 0.5_dp, -0.5_dp,  0.5_dp, &
                             0.5_dp,  0.5_dp,  0.5_dp, &
                            -0.5_dp,  0.5_dp,  0.5_dp /), &
                             shape(bcc_mat))
                             
                              
        
        call matrix_inverse(fcc_mat,3,fcc_inv)
        call matrix_inverse(bcc_mat,3,bcc_inv)

        max_basisnum = 0
        basisnum(:) = 0
        basis_type(:,:) = 0
        ng_node(:) = 0
        node_num = 0
        node_atoms = 0
        ele_num = 0
        atom_num = 0
    end subroutine lattice_init

    subroutine cell_init(lapa,esize,ele_type, orient_mat, cell_mat)
        !This subroutine uses the user provided information to transform the finite element cell to the correct
        !size, orientation, and dimensions using the esize, lattice parameter, element_type, and orientation provided 
        !by the user 

        integer, intent(in) :: esize
        real(kind=dp), intent(in) :: lapa, orient_mat(3,3)
        character(len=100), intent(in) :: ele_type

        real(kind=dp), dimension(3,max_ng_node), intent(out) :: cell_mat

        integer :: inod, i
        real(kind=dp), dimension(3,max_ng_node) :: adjustVar

        adjustVar(:,:) = 0.0_dp
        vflag = .false.

        select case(trim(ele_type))

        case('fcc')
            if(lmpcac) then 
                do inod = 1, 8
                    do i = 1,3
                        if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
                            adjustVar(i,inod) = -0.5_dp
                        else 
                            adjustVar(i, inod) = 0.5_dp
                        end if
                    end do
                end do
                adjustVar(:,1:8) = matmul(fcc_mat, adjustVar(:,1:8))
            end if
            cell_mat(:, 1:8) = fcc_cell + adjustVar(:,1:8)
            cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, cell_mat(:,1:8)))
        case('bcc')
            cell_mat(:,1:8) = bcc_cell
            cell_mat(:,1:8) = lapa* ((esize-1)*matmul(orient_mat, cell_mat(:,1:8)))
        case default
            print *, "Element type ", trim(ele_type), " currently not accepted"
            stop
        end select

    end subroutine cell_init

    subroutine alloc_ele_arrays(n,m)
        !This subroutine used to provide initial allocation for the atom and element arrays

        integer, intent(in) :: n, m !n-size of element arrays, m-size of atom arrays

        integer :: allostat
        
        !Allocate element arrays
        if(n > 0) then 
            allocate(type_ele(n), tag_ele(n), size_ele(n), lat_ele(n), sbox_ele(n), r_node(3,max_basisnum, max_ng_node,n), &
                       stat=allostat)
            if(allostat > 0) then 
                print *, "Error allocating element arrays in elements.f90 because of: ", allostat
                stop
            end if
        end if

        if(m > 0) then 
            !Allocate atom arrays
            allocate(type_atom(m), sbox_atom(m), tag_atom(m), r_atom(3,m), stat=allostat)
            if(allostat > 0) then 
                print *, "Error allocating atom arrays in elements.f90 because of: ", allostat
                stop
            end if
        end if    
    end subroutine

    subroutine grow_ele_arrays(n, m)

        integer, intent(in) :: n, m

        integer :: ele_size, atom_size, buffer_size
        integer, allocatable :: temp_int(:)
        real(kind=dp), allocatable :: temp_ele_real(:,:,:,:), temp_real(:,:)
        character(len=100), allocatable :: char_temp(:)

        !The default size we grow the 
        buffer_size = 1024
        
        !First check to make sure if it is allocated
        if (allocated(size_ele)) then 

            !Figure out the size of the atom and element arrays
            ele_size = size(size_ele)

            !Check if we need to grow the ele_size, if so grow all the variables
            if ( n+ele_num > size(size_ele)) then 
                
                allocate(temp_int(n+ele_size+buffer_size))
                temp_int(1:ele_size) = lat_ele(1:ele_size)
                temp_int(ele_size+1:) = 0
                call move_alloc(temp_int, lat_ele)

                allocate(temp_int(n+ele_size+buffer_size))
                temp_int(1:ele_size) = tag_ele(1:ele_size)
                temp_int(ele_size+1:) = 0
                call move_alloc(temp_int, tag_ele)

                allocate(temp_int(n+ele_size+buffer_size))
                temp_int(1:ele_size) = size_ele(1:ele_size)
                temp_int(ele_size+1:) = 0
                call move_alloc(temp_int, size_ele)

                allocate(temp_int(n+ele_size+buffer_size))
                temp_int(1:ele_size) = sbox_ele(1:ele_size)
                temp_int(ele_size+1:) = 0
                call move_alloc(temp_int, sbox_ele)

                allocate(char_temp(n+ele_size+buffer_size))
                char_temp(1:ele_size) = type_ele(1:ele_size)
                call move_alloc(char_temp, type_ele)

                allocate(temp_ele_real(3, max_basisnum, max_ng_node, n+ele_size+buffer_size))
                temp_ele_real(:,:,:,1:ele_size) = r_node(:,:,:,1:ele_size)
                temp_ele_real(:,:,:,ele_size+1:) = 0.0_dp
                call move_alloc(temp_ele_real, r_node)
            end if
        else 
            call alloc_ele_arrays(n,0)
        end if

        !Now grow atom arrays if needed
        if (allocated(type_atom)) then 
            atom_size = size(type_atom)
            if (m+atom_num > atom_size) then 
                allocate(temp_int(m+atom_size+buffer_size))
                temp_int(1:atom_size) = type_atom
                temp_int(atom_size+1:) = 0
                call move_alloc(temp_int, type_atom)

                allocate(temp_int(m+atom_size+buffer_size))
                temp_int(1:atom_size) = tag_atom
                temp_int(atom_size+1:) = 0
                call move_alloc(temp_int, tag_atom)

                allocate(temp_int(m+atom_size+buffer_size))
                temp_int(1:atom_size) = sbox_atom
                temp_int(atom_size+1:) = 0
                call move_alloc(temp_int, sbox_atom)

                allocate(temp_real(3,m+atom_size+buffer_size))
                temp_real(:,1:atom_size) = r_atom
                temp_real(:, atom_size+1:) = 0.0_dp
                call move_alloc(temp_real, r_atom)
            end if
        else
            call alloc_ele_arrays(0,m)
        end if
    end subroutine

    subroutine add_element(tag, type, size, lat, sbox, r)
        !Subroutine which adds an element to the element arrays
        integer, intent(in) :: size, lat, sbox, tag
        character(len=100), intent(in) :: type
        real(kind=dp), intent(in) :: r(:,:,:)

        integer :: newtag

        ele_num = ele_num + 1
        node_num = node_num + ng_node(lat)
        node_atoms = node_atoms + ng_node(lat)*basisnum(lat)

        if (tag==0) then
            newtag = ele_num !If we don't assign a tag then pass the tag as the ele_num
        else
            newtag = tag
        end if

        !Check to see if we need to grow the arrays
        call grow_ele_arrays(1,0)
        tag_ele(ele_num) = newtag
        type_ele(ele_num) = type
        size_ele(ele_num) = size
        lat_ele(ele_num) = lat
        sbox_ele(ele_num) = sbox
        r_node(:,:,:,ele_num) = r(:,:,:)


    end subroutine add_element

    subroutine add_atom(tag, type, sbox, r)
        !Subroutine which adds an atom to the atom arrays
        integer, intent(in) :: type, sbox, tag
        real(kind=dp), intent(in), dimension(3) :: r

        integer :: newtag

        atom_num = atom_num+1
        if(tag==0) then
            newtag = atom_num !If we don't assign a tag then pass the tag as the atom_num
        else
            newtag = tag
        end if
        !Check to see if we need to grow the arrays
        call grow_ele_arrays(0,1)
        tag_atom(atom_num) = newtag
        type_atom(atom_num) = type
        r_atom(:,atom_num) = r(:)
        sbox_atom(atom_num) = sbox

    end subroutine add_atom

    subroutine add_atom_type(type, inttype, force_new)
        !This subroutine adds a new atom type to the module list of atom types
        character(len=2), intent(in) :: type
        integer, intent(out) :: inttype

        logical, optional, intent(in) :: force_new

        integer :: i
        logical :: exists, force
        
        if(present(force_new)) then 
            force = force_new
        else 
            force = .false.
        end if

        exists = .false.
        if(.not.force) then 
            do i=1, 10
                if(type == type_to_name(i)) then 
                    exists = .true.
                    inttype = i
                    exit
                end if
            end do 
        end if

        if (exists.eqv..false.) then 
            atom_types = atom_types+1
            if(atom_types > 10) then 
                print *, "Defined atom types are greater than 10 which is currently not supported."
                stop 3
            end if
            type_to_name(atom_types) = type
            inttype = atom_types
        end if
        return
    end subroutine add_atom_type

    subroutine def_ng_node(n, element_types)
        !This subroutine defines the maximum node number among n element types
        integer, intent(in) :: n !Number of element types
        character(len=100), dimension(n) :: element_types !Array of element type strings

        integer :: i

        do i=1, n
            select case(trim(adjustl(element_types(i))))
            case('fcc')
                ng_node(i) = 8
            case('bcc')
                ng_node(i) = 8
            case('20fcc')
                ng_node(i) = 20
            end select

            if(ng_node(i) > max_ng_node) max_ng_node = ng_node(i)
        end do 
    end subroutine def_ng_node

    subroutine set_max_esize
        !This subroutine sets the maximum esize
        if(allocated(size_ele)) then 
            max_esize=maxval(size_ele)
        else 
            max_esize = 2
        end if
    end subroutine

    subroutine interpolate_atoms(type, esize, lat_type, r_in, type_interp, r_interp, eng, f, v, data_interp)
        !This subroutine returns the interpolated atoms from the elements. 

        !Arguments
        character(len=*), intent(in) :: type !The type of element that it is
        integer, intent(in) :: esize !The number of atoms per side
        integer, intent(in) :: lat_type !The integer lattice type of the element
        real(kind=dp), dimension(:,:,:), intent(in) :: r_in !Nodal positions
        integer, dimension(:), intent(out) :: type_interp !Interpolated atomic positions
        real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: r_interp !Interpolated atomic positions
        real(kind = dp), optional, intent(in) :: eng(:,:), f(:,:,:), &
                                                 v(:,:,:)
        real(kind=dp), dimension(10, max_basisnum*max_esize**3), optional,intent(out) :: data_interp !Interpolated atomic positions

        !Internal variables
        integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp
        real(kind=dp) :: a_shape(max_ng_node)
        real(kind=dp) :: t, s, r
        
        !Initialize some variables
        r_interp(:,:) = 0.0_dp
        type_interp(:) = 0
        if(present(data_interp)) data_interp = 0.0_dp
        ia = 0

        !Define bnum based on the input lattice type. If lat_type=0 then we are interpolating lattice points which means
        !the basis is 0,0,0, and the type doesn't matter

        select case(lat_type)
        case(0)
            bnum=1
            lat_type_temp = 1
        case default
            bnum = basisnum(lat_type)
            lat_type_temp = lat_type
        end select


        select case(type)
        case('fcc','bcc')
            !Now loop over all the possible sites
            do it = 1, esize
                t = (1.0_dp*(it-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
                do is =1, esize
                    s = (1.0_dp*(is-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
                    do ir = 1, esize
                        r = (1.0_dp*(ir-1) - (esize-1)/2)/(1.0_dp*(esize-1)/2)
                        call rhombshape(r,s,t,a_shape)

                        do ibasis = 1, bnum
                            ia = ia + 1
                            do inod = 1, 8
                                type_interp(ia) = basis_type(ibasis,lat_type_temp)
                                r_interp(:,ia) = r_interp(:,ia) + a_shape(inod) * r_in(:,ibasis,inod)

                                if(present(data_interp)) then 
                                    !If data is present then interpolate data arrays as well
                                    data_interp(1,ia) = data_interp(1,ia) + eng(ibasis, inod)*a_shape(inod)
                                    data_interp(2:4,ia) = data_interp(2:4,ia) + f(:, ibasis, inod)*a_shape(inod)
                                    data_interp(5:10,ia) = data_interp(5:10,ia) + v(:, ibasis, inod)*a_shape(inod)
                                end if
                            end do
                        end do

                    end do
                end do
            end do 
        case('20fcc')
            !Now loop over all the possible sites
            do it = 1, esize
                t = (1.0_dp*(it-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
                do is =1, esize
                    s = (1.0_dp*(is-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
                    do ir = 1, esize
                        r = (1.0_dp*(ir-1) - (esize-1)/2)/(1.0_dp*(esize-1)/2)
                        call quad_rhombshape(r,s,t,a_shape)

                        do ibasis = 1, bnum
                            ia = ia + 1
                            do inod = 1, 20
                                type_interp(ia) = basis_type(ibasis,lat_type_temp)
                                r_interp(:,ia) = r_interp(:,ia) + a_shape(inod) * r_in(:,ibasis,inod)

                                if(present(data_interp)) then 
                                    !If data is present then interpolate data arrays as well
                                    data_interp(1,ia) = data_interp(1,ia) + eng(ibasis, inod)*a_shape(inod)
                                    data_interp(2:4,ia) = data_interp(2:4,ia) + f(:, ibasis, inod)*a_shape(inod)
                                    data_interp(5:10,ia) = data_interp(5:10,ia) + v(:, ibasis, inod)*a_shape(inod)
                                end if
                            end do
                        end do

                    end do
                end do
            end do 
        end select

        if (ia /= bnum*esize**3) then 
            print *, "Incorrect interpolation"
            stop 3
        end if
        return
    end subroutine interpolate_atoms

    subroutine rhombshape(r,s,t, shape_fun)
        !Shape function for rhombohedral elements
        real(kind=8), intent(in) :: r, s, t
        real(kind=8), intent(out) :: shape_fun(:)

        shape_fun(1) =  (1.0-r)*(1.0-s)*(1.0-t)/8.0
        shape_fun(2) =  (1.0+r)*(1.0-s)*(1.0-t)/8.0
        shape_fun(3) =  (1.0+r)*(1.0+s)*(1.0-t)/8.0
        shape_fun(4) =  (1.0-r)*(1.0+s)*(1.0-t)/8.0
        shape_fun(5) =  (1.0-r)*(1.0-s)*(1.0+t)/8.0
        shape_fun(6) =  (1.0+r)*(1.0-s)*(1.0+t)/8.0
        shape_fun(7) =  (1.0+r)*(1.0+s)*(1.0+t)/8.0
        shape_fun(8) =  (1.0-r)*(1.0+s)*(1.0+t)/8.0


        return
    end subroutine rhombshape

    subroutine quad_rhombshape(r,s,t, shape_fun)
        !Shape function for 20 node quadratic rhombohedral elements
        real(kind=8), intent(in) :: r, s, t
        real(kind=8), intent(out) :: shape_fun(:)

        !Corner nodes
        shape_fun(1) = (1.0_dp-r)*(1.0_dp-s)*(1.0_dp-t)*(-r-s-t-2)/8.0_dp
        shape_fun(2) = (1.0_dp+r)*(1.0_dp-s)*(1.0_dp-t)*(r-s-t-2)/8.0_dp
        shape_fun(3) = (1.0_dp+r)*(1.0_dp+s)*(1.0_dp-t)*(r+s-t-2)/8.0_dp
        shape_fun(4) = (1.0_dp-r)*(1.0_dp+s)*(1.0_dp-t)*(-r+s-t-2)/8.0_dp
        shape_fun(5) = (1.0_dp-r)*(1.0_dp-s)*(1.0_dp+t)*(-r-s+t-2)/8.0_dp
        shape_fun(6) = (1.0_dp+r)*(1.0_dp-s)*(1.0_dp+t)*(r-s+t-2)/8.0_dp
        shape_fun(7) = (1.0_dp+r)*(1.0_dp+s)*(1.0_dp+t)*(r+s+t-2)/8.0_dp
        shape_fun(8) = (1.0_dp-r)*(1.0_dp+s)*(1.0_dp+t)*(-r+s+t-2)/8.0_dp

        !Side nodes, first node r is zero
        shape_fun(9) = (1-r*r)*(1-s)*(1-t)/4.0_dp
        shape_fun(11) = (1-r*r)*(1+s)*(1-t)/4.0_dp
        shape_fun(17) = (1-r*r)*(1-s)*(1+t)/4.0_dp
        shape_fun(19) = (1-r*r)*(1+s)*(1+t)/4.0_dp

        !node s is zero
        shape_fun(10) = (1+r)*(1-s*s)*(1-t)/4.0_dp
        shape_fun(12) = (1-r)*(1-s*s)*(1-t)/4.0_dp
        shape_fun(18) = (1+r)*(1-s*s)*(1+t)/4.0_dp
        shape_fun(20) = (1-r)*(1-s*s)*(1+t)/4.0_dp

        !node t is zero
        shape_fun(13) = (1-r)*(1-s)*(1-t*t)/4.0_dp
        shape_fun(14) = (1+r)*(1-s)*(1-t*t)/4.0_dp
        shape_fun(15) = (1+r)*(1+s)*(1-t*t)/4.0_dp
        shape_fun(16) = (1-r)*(1+s)*(1-t*t)/4.0_dp


    end subroutine quad_rhombshape

    subroutine delete_atoms(num, in_index)
        !This subroutine deletes atoms from the atom arrays 
        integer, intent(in) :: num
        integer, intent(in), dimension(num) :: in_index
        
        real(kind=dp), dimension(num) :: for_sort
        integer, dimension(num) :: sorted_index

        integer :: i

        for_sort = in_index 

        call dpquicksort(for_sort)
        do i = 1, num
            sorted_index(i) = nint(for_sort(i))
        end do

        !We go from largest index to smallest index just to make sure that we don't miss 
        !accidentally overwrite values which need to be deleted
        do i = num, 1, -1
            if(sorted_index(i) == atom_num) then 
                r_atom(:,sorted_index(i)) = 0.0_dp
                type_atom(sorted_index(i)) = 0
            else
                r_atom(:,sorted_index(i)) = r_atom(:, atom_num)
                type_atom(sorted_index(i)) = type_atom(atom_num)
            end if            
            atom_num = atom_num - 1
        end do 
    end subroutine

    subroutine delete_elements(num, in_index)
        !This subroutine deletes elements from the element array 
        integer, intent(in) :: num
        integer, intent(in), dimension(num) :: in_index

        integer :: i
        real(kind=dp), dimension(num) :: for_sort
        integer, dimension(num) :: sorted_index

        
        for_sort = in_index
        call dpquicksort(for_sort)
        do i = 1, num
            sorted_index(i) = nint(for_sort(i))
        end do

        !We go from largest index to smallest index just to make sure that we don't miss 
        !accidentally overwrite values which need to be deleted
        do i = num, 1, -1
            node_num = node_num - ng_node(lat_ele(sorted_index(i)))
            node_atoms = node_atoms - ng_node(lat_ele(sorted_index(i)))*basisnum(lat_ele(sorted_index(i)))
            if(sorted_index(i) == ele_num) then 
                r_node(:,:,:,sorted_index(i)) = 0.0_dp
                type_ele(sorted_index(i)) =''
                size_ele(sorted_index(i)) = 0
                lat_ele(sorted_index(i)) = 0
                sbox_ele(sorted_index(i)) = 0
                tag_ele(sorted_index(i)) = 0
            else
                r_node(:,:,:,sorted_index(i)) = r_node(:,:,:,ele_num)
                type_ele(sorted_index(i)) = type_ele(ele_num)
                size_ele(sorted_index(i)) = size_ele(ele_num)
                lat_ele(sorted_index(i)) = lat_ele(ele_num)
                sbox_ele(sorted_index(i)) = sbox_ele(ele_num)
                tag_ele(sorted_index(i)) = tag_ele(ele_num)
            end if            
            ele_num = ele_num - 1
        end do 

    end subroutine delete_elements

    subroutine wrap_atoms
        !This subroutine wraps atoms back into the simulation cell if they have exited for any reason

        integer :: i
        do i = 1, atom_num
            call apply_periodic(r_atom(:,i))
        end do 

    end subroutine wrap_atoms
    subroutine wrap_elements
        integer :: i, inod, ibasis, j, node_in_bd(3,max_ng_node)
        real(kind =dp) :: box_len(3)

        do j = 1, 3
            box_len(j) = box_bd(2*j) - box_bd(2*j-1)
        end do
        print *, box_bd
        do i = 1, ele_num
            node_in_bd = 0
            do inod = 1, ng_node(lat_ele(i))
                do ibasis = 1, basisnum(lat_ele(i))
                    do j = 1, 3
                        if(r_node(j,ibasis,inod,i) < box_bd(2*j-1)) then 
                            node_in_bd(j,inod) = 1
                        else if(r_node(j,ibasis,inod,i) > box_bd(2*j)) then 
                            node_in_bd(j,inod) = -1
                        end if
                    end do
                end do
            end do

            do j = 1, 3
                if(all(node_in_bd(j,:) == 1)) then 
                    r_node(j,:,:,i) = r_node(j,:,:,i) + box_len(j)
                else if(all(node_in_bd(j,:) == -1)) then 
                    r_node(j,:,:,i) = r_node(j,:,:,i) - box_len(j)
                end if
            end do
        end do
    end subroutine wrap_elements


    subroutine def_new_box
        !This subroutine calculates new box boundaries based on minimum and maximum nodal/atomic positions
        integer :: i, j, inod, ibasis

        real(kind=dp) :: max_bd(3), min_bd(3)

        max_bd(:) = -huge(1.0_dp)
        min_bd(:) = huge(1.0_dp)
do i = 1, atom_num
            do j = 1, 3
                if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + tol 
                if (r_atom(j,i) < min_bd(j)) min_bd(j) = r_atom(j,i) - tol 
            end do 
        end do 

        do i = 1, ele_num
            do inod = 1, ng_node(lat_ele(i))
                do ibasis = 1, basisnum(lat_ele(i))
                    do j = 1, 3
                        if (r_node(j,ibasis,inod,i) > max_bd(j)) max_bd(j) = r_node(j,ibasis,inod,i) + tol
                        if (r_node(j,ibasis,inod,i) < min_bd(j)) min_bd(j) = r_node(j,ibasis,inod,i) - tol
                    end do 
                end do 
            end do 
        end do 

        do j = 1, 3
            if(box_bc(j:j) == 's') then 
                box_bd(2*j) = max_bd(j)
                box_bd(2*j-1) = min_bd(j)
            end if
        end do 
    end subroutine

    recursive subroutine parse_pos(i, pos_string, pos)
        !This subroutine parses the pos command allowing for command which include inf
        integer, intent(in) ::  i !The dimension of the position
        character(len=100), intent(in) :: pos_string !The position string
        real(kind=dp), intent(out) :: pos !The output parsed position value

        integer :: iospara, face, ele, randsize
        real(kind=dp) :: rand, rone, rtwo, rand_ele_pos(3)
        character(len=100) :: cone, ctwo

        iospara = 0
        
        if(trim(adjustl(pos_string)) == 'inf') then 
            pos=box_bd(2*i)

        else if(trim(adjustl(pos_string)) == '-inf') then
            pos=box_bd(2*i-1)

        else if(trim(adjustl(pos_string)) == 'mid') then 
            pos=(box_bd(2*i)+box_bd(2*i-1))/2

        else if (trim(adjustl(pos_string)) == 'rand') then 
            call random_number(rand)
            pos = (box_bd(2*i)-box_bd(2*i-1))*rand + box_bd(2*i-1)

        else if (index(pos_string,'rande')>0) then 
            !First select a random element 
            call random_number(rand)
            ele= 1 + floor(ele_num*rand)
            !Now read the rest of the command which specifies the face we need and check 
            !to make sure it's an accepted face number
            cone = pos_string(index(pos_string, '[')+1:index(pos_string, '[')+1)
            
            !Check to see if we also pass an element size, if it was passed then make sure our 
            !random element is the right size 
            if(index(pos_string,':') > 0) then 
                ctwo = pos_string(index(pos_string, ':')+1:index(pos_string,']')-1) 
                read(ctwo, *) randsize
                do while(randsize /= size_ele(ele))
                    call random_number(rand)
                    ele= 1 + floor(ele_num*rand)
                end do 
            end if

            if ((face < 1).or.(face > 6)) stop "Current face number must be 1,2,3,4,5,6. Please check documentation" 
            !Now get the position
            call offset_pos(ele, face, rand_ele_pos)
            pos = rand_ele_pos(i)

        else if (index(pos_string,'rand')>0) then 
            call random_number(rand)
            cone = pos_string(index(pos_string, '[')+1:index(pos_string,':')-1)
            call parse_pos(i, cone, rone)
            ctwo = pos_string(index(pos_string, ':')+1:index(pos_string,']')-1)
            call parse_pos(i, ctwo, rtwo)
            pos = (rtwo - rone)*rand + rone            

        else if ((index(pos_string,'-') > 0).and.(index(pos_string,'inf')>0)) then 
            !Now extract the number we are reducing from infinity
            if(index(pos_string,'inf') < index(pos_string,'-')) then 
                read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
            else
                read(pos_string(1:index(pos_string,'-')-1), *, iostat=iospara) pos
            end if
            pos = box_bd(2*i) - pos

        else if ((index(pos_string,'-') > 0).and.(index(pos_string,'mid')>0)) then 
            !Now extract the number we are reducing from midinity
            if(index(pos_string,'mid') < index(pos_string,'-')) then 
                read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
            else
                read(pos_string(1:index(pos_string,'-')-1), *, iostat=iospara) pos
            end if
            pos = (box_bd(2*i)+box_bd(2*i-1))/2.0_dp - pos

        else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then 
            !Now extract the number we are reducing from infinity
           if(index(pos_string,'inf') < index(pos_string,'+')) then 
                read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
            else
                read(pos_string(1:index(pos_string,'+')-1), *, iostat=iospara) pos
            end if
            pos = box_bd(2*i-1) + pos

        else if ((index(pos_string,'+') > 0).and.(index(pos_string,'mid')>0)) then 
            !Now extract the number we are reducing from midinity
           if(index(pos_string,'mid') < index(pos_string,'+')) then 
                read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
            else
                read(pos_string(1:index(pos_string,'+')-1), *, iostat=iospara) pos
            end if
            pos = (box_bd(2*i)+box_bd(2*i-1))/2.0_dp + pos

        else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then 
            !Now extract the number we are reducing from infinity
           if(index(pos_string,'inf') < index(pos_string,'*')) then 
                read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos
            else
                read(pos_string(1:index(pos_string,'*')-1), *, iostat=iospara) pos
            end if
            pos = (box_bd(2*i)-box_bd(2*i-1))*pos + box_bd(2*i-1)

        else
            read(pos_string, *, iostat=iospara) pos

        end if

        if (iospara > 0) then 
            print *, "Error reading position argument ", trim(adjustl(pos_string)), ". Please reformat and try again."
        end if
    end subroutine parse_pos

    subroutine offset_pos(ie, iface, pos)
    !This returns a position slightly offset from the center of an element face
    !This is used to return a random position within an element discontinuity

    integer, intent(in) :: ie !Element index
    integer, intent(in) :: iface !Face number between 1 and 6 
    real(kind=dp), dimension(3), intent(out) :: pos !Position vector

    !Other variables we need
    integer :: esize
    real(kind=dp) :: orient(3,3), ori_inv(3,3), lp, r_cubic_node(3,8)

    !First find the offset vector for the face in the untransformed cubic cell
    esize = size_ele(ie)
    select case(iface)
        case(1)
            pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp) /)
        case(2)
            pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /)
        case(3)
            pos = (/ (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
        case(4)
            pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp /)
        case(5)
            pos = (/ -10.0_dp**(-2.0_dp), real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
        case(6)
            pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**(-2.0_dp) /)
    end select

    !Now transform it to real space and adjust it to the position of the element in the first node.
    select case (trim(adjustl(type_ele(ie))))

    case('fcc')
        !First we have to extract the element lattice parameter
        call matrix_inverse(sub_box_ori(:,:,sbox_ele(ie)),3,ori_inv)
        r_cubic_node = r_node(:,1,:,ie)
        r_cubic_node = matmul(fcc_inv,matmul(ori_inv,r_cubic_node))
        lp = (maxval(r_cubic_node(1,:))-minval(r_cubic_node(1,:)))/(esize-1)
        pos = matmul(sub_box_ori(:,:,sbox_ele(ie)),matmul(fcc_mat,pos))*lp
        pos = pos + r_node(:,1, 1, ie)

    case default
        print *, trim(adjustl(type_ele(ie))), " is not a currently accepted element type for random element positions"
        stop 3

    end select
    end subroutine

    subroutine lattice_map(in_bnum, in_btypes, in_ngnodes, in_lapa, lat_type)
        !This subroutine maps an input lattice type to either a new lattice type or an existing one depending on basis_type and 
        !number of nodes at the atoms

        integer, intent(in) :: in_ngnodes, in_bnum, in_btypes(10) !Input variables
        real(kind=dp), intent(in) :: in_lapa
        integer, intent(out) :: lat_type

        integer j, ibasis

        lat_type = 0
        lat_loop:do j = 1, lattice_types
            !Check all the lattice level variables
            if ((basisnum(j) == in_bnum).and.(ng_node(j) == in_ngnodes).and.(is_equal(lapa(j),in_lapa))) then 
                !Now check lattice level variables 
                do ibasis = 1, basisnum(j)
                    if(basis_type(ibasis,j) /= in_btypes(ibasis)) cycle lat_loop
                end do
                lat_type = j
                exit lat_loop
            end if
        end do lat_loop

        !If it doesn't match an existing lattice type we add it 
        if( lat_type == 0) then 
            lattice_types = lattice_types + 1
            basisnum(lattice_types) = in_bnum
            basis_type(:,lattice_types) = in_btypes
            ng_node(lattice_types) = in_ngnodes
            lapa(lattice_types) = in_lapa
            lat_type = lattice_types
        end if

    end subroutine lattice_map

    subroutine get_interp_pos(i,j,k, ie, rout)
        !This returns the position of an interpolated basis from an element ie.
        !i, j, k should be in natural coordinates

        integer, intent(in) :: i, j, k
        real(kind=dp), dimension(3,max_basisnum), intent(out) :: rout

        integer :: ie, ibasis, inod
        real(kind=dp) :: a_shape(8), r, s, t

        r = (1.0_dp*(i-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2)
        s = (1.0_dp*(j-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2)
        t = (1.0_dp*(k-1)-(size_ele(ie)-1)/2)/(1.0_dp*(size_ele(ie)-1)/2)
        rout(:,:) = 0
        do ibasis = 1, basisnum(lat_ele(ie))
            do inod = 1, ng_node(lat_ele(ie)) 
                call rhombshape(r,s,t,a_shape)
                rout(:,ibasis) = rout(:,ibasis) + a_shape(inod) * r_node(:,ibasis,inod,ie)
            end do
        end do

    end subroutine get_interp_pos

    subroutine alloc_dat_arrays(n,m)
        !This subroutine used to provide initial allocation for the atom and element data arrays
        integer, intent(in) :: n,m !n-size of element arrays, m-size of atom arrays
        integer :: allostat

        !Allocate element arrays
        if (n > 0) then 
            if (allocated(force_node)) then 
                deallocate(force_node, virial_node, energy_node)
            end if 
            allocate(force_node(3,max_basisnum, max_ng_node, n), &
                     virial_node(6,max_basisnum, max_ng_node, n), &
                     energy_node(max_basisnum,max_ng_node,n),    &
                     stat=allostat)
            if(allostat > 0) then 
                print *, "Error allocating element data arrays in mode_metric because of:", allostat
                stop
            end if

        end if

        if (m > 0) then
            if (allocated(force_atom)) then 
                deallocate(force_atom, virial_atom, energy_atom)
            end if 
            allocate(force_atom(3, m), &
                     virial_atom(6, m), &
                     energy_atom(m),    &
                     stat=allostat)
            if(allostat > 0) then 
                print *, "Error allocating atom data arrays in mode_metric because of:", allostat
                stop
            end if
        end if

    end subroutine 

    subroutine alloc_vel_arrays(n,m)
        !This subroutine used to provide initial allocation for the atom and element data arrays
        integer, intent(in) :: n,m !n-size of element arrays, m-size of atom arrays
        integer :: allostat

        !Allocate element arrays
        if (n > 0) then 
            if (allocated(vel_node)) then 
                deallocate(vel_node)
            end if 
            allocate(vel_node(3,max_basisnum,max_ng_node,n), stat=allostat)
            if(allostat > 0) then 
                print *, "Error allocating element data arrays in mode_metric because of:", allostat
                stop
            end if

        end if

        if (m > 0) then
            if (allocated(vel_atom)) then 
                deallocate(vel_atom)
            end if 
            allocate(vel_atom(3,m), stat=allostat)
            if(allostat > 0) then 
                print *, "Error allocating atom data arrays in mode_metric because of:", allostat
                stop
            end if
        end if

    end subroutine alloc_vel_arrays


    subroutine add_atom_data(ia, eng, force, virial)
        !Function which sets the atom data arrays 
        integer, intent(in) :: ia
        real(kind=dp), intent(in) :: eng, force(3), virial(6)

        energy_atom(ia) = eng
        force_atom(:,ia) = force(:) 
        virial_atom(:,ia) = virial(:) 
        vflag = .true.
        return
    end subroutine add_atom_data

    subroutine add_element_data(ie, eng, force, virial)
        !Function which sets the element data arrays
        integer, intent(in) :: ie
        real(kind=dp), intent(in) :: eng(max_basisnum, max_ng_node), &
                                     force(3,max_basisnum, max_ng_node), &
                                     virial(6,max_basisnum,max_ng_node)
        integer :: inod
        vflag = .true.
        energy_node(:,:,ie)   = eng
        force_node(:,:,:,ie)  = force
        virial_node(:,:,:,ie) = virial
        
        return
    end subroutine add_element_data

    subroutine reset_data
        !This function resets all of the data arrays for the elements and atoms
        atom_num = 0
        ele_num = 0
        node_num = 0
    end subroutine reset_data

    function ele_in_bounds(bd, etype, esize, lat_type, r_in )
        real(kind=dp), intent(in) :: bd(6), r_in(3, max_basisnum, max_ng_node)
        integer, intent(in) ::  lat_type, esize
        character(len=*), intent(in) :: etype

        integer :: i
        real(kind=dp) :: rinterp(3, max_basisnum*max_esize**3) 
        integer :: tinterp(max_basisnum*max_esize**3) 
        logical :: ele_in_bounds

        ele_in_bounds=.false.
        call interpolate_atoms(etype, esize, lat_type, r_in, tinterp, rinterp) 

        do i=1, esize**3
            if(in_block_bd(rinterp(:,i), bd))then 
                ele_in_bounds = .true.
                exit
            end if
        end do

        return
    end function ele_in_bounds

end module elements