CACmb/src/mode_create.f90

module mode_create
    !This mode is intended for creating element/atom regions and writing them to specific files.

    use parameters
    use atoms   
    use io
    use subroutines
    use elements
    use box

    implicit none

    character(len=100) :: name, element_type
    real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), &
                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3)
    integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
               basis_pos(3,10)
    logical :: dup_flag, dim_flag

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

    subroutine create(arg_pos)
        ! Main subroutine which controls execution

        character(len=100) :: textholder
        integer, intent(out) :: arg_pos

        integer :: i, ibasis, inod
        real(kind=dp), allocatable :: r_node_temp(:,:,:)

        !Initialize default parameters
        orient = reshape((/ 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp /), shape(orient))
        cell_mat(:,:)=0.0_dp
        name =''
        element_type = ''
        esize=0
        lattice_parameter=0.0_dp
        duplicate(:) = 0
        box_len(:) = 0.0_dp
        dup_flag = .false.
        dim_flag = .false.
        basisnum = 0
        lat_ele_num = 0
        lat_atom_num = 0

        !First we parse the command
        call parse_command(arg_pos)

        !Now we setup the unit element and call other init subroutines
        call def_ng_node(1, element_type)

        allocate(r_node_temp(3,max_basisnum,max_ng_node))

        !Get the inverse orientation matrix we will need later 
        call matrix_inverse(orient,3,orient_inv)

        if(dup_flag) then 
            
            !Define box vertices
            do i = 1, 8
                box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + (origin(:)/lattice_parameter) 
            end do
            !Now get the rotated box vertex positions in lattice space. Should be integer units
            box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
            !Find the new maxlen
            maxbd = maxval(matmul(orient,matmul(fcc_mat,box_lat_vert)),2)
            do i = 1, 3
                box_bd(2*i) = maxval(box_vert(i,:)) - 0.25_dp*lattice_space(i)
                box_bd(2*i-1) = origin(i)-0.25_dp*lattice_space(i)
            end do



        else if(dim_flag) then 
            !As a note everything is defined so that the lattice parameter is multiplied in at the end 
            !so we have to divide all the real Angstroms units by the lattice parameter

            !Define box_vertices
            do i = 1, 8
                box_vert(:,i) = (cubic_cell(:,i)*box_len(:) + origin(:))/lattice_parameter
            end do
            !Now get the rotated box vertex positions in lattice space. Should be integer units
            box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1

            !Now get the box_bd in lattice units
            do i = 1, 3
                box_bd(2*i) = (box_len(i)+origin(i))/lattice_parameter
                box_bd(2*i-1) = origin(i)/lattice_parameter
            end do 
        else

            call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) 
            !If the user doesn't pass any build instructions than we just put the cell mat into the element_array 
            call alloc_ele_arrays(1,0)

            !Add the basis atoms to the unit cell
            do inod = 1, max_ng_node
                do ibasis = 1, basisnum(1)
                    r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis) + origin(:)
                end do
            end do
            do i = 1,3
                box_bd(2*i) = maxval(r_node_temp(i,:,:))
                box_bd(2*i-1) = origin(i)
            end do
            call add_element(element_type, esize, 1, r_node_temp)
        end if
        
        !If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays
        if(dup_flag.or.dim_flag) then 
            !Call the build function with the correct transformation matrix
            select case(trim(adjustl(element_type)))
            case('fcc')

                call build_with_rhomb(box_lat_vert, fcc_mat)
            case default
                print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", &
                         "element type"
                stop 3
            end select
            !Now that it is built multiply by the lattice parameter
            box_bd = box_bd*lattice_parameter

            !Allocate variables 
            call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1))
            if(lat_atom_num > 0) then 
                do i = 1, lat_atom_num
                    do ibasis = 1, basisnum(1)
                        call add_atom(basis_type(ibasis, 1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
                    end do
                end do
                deallocate(r_atom_lat)
            end if

            if(lat_ele_num > 0) then 
                do i = 1, lat_ele_num
                    do inod= 1, ng_node(1)
                        do ibasis = 1, basisnum(1)
                            r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis)
                        end do
                    end do 
                    call add_element(element_type, esize, 1, r_node_temp)
                end do 
            end if
        end if

        !The last thing we do is setup the sub_box_boundaries
        call alloc_sub_box(1)
        sub_box_num = 1
        sub_box_ori(:,:,1) = orient
        sub_box_bd(:,1) = box_bd
        sub_box_array_bd(1,:,1) = 1
        sub_box_array_bd(2,1,1) = atom_num
        sub_box_array_bd(2,2,1) = ele_num
    end subroutine create
    !This subroutine parses the command and pulls out information needed for mode_create
    subroutine parse_command(arg_pos)

        integer, intent(out) :: arg_pos
        integer ::  ori_pos, i, j, arglen, stat
        character(len=100) :: textholder
        character(len=8) :: orient_string


        !Pull out all required positional arguments
        call get_command_argument(2, name, arglen)
        if(arglen==0) STOP "Name is missing in mode create"

        call get_command_argument(3,element_type, arglen)
        if(arglen==0) STOP "Element_type is missing in mode create"

        call get_command_argument(4,textholder, arglen)
        if(arglen==0) STOP "Lattice Parameter is missing in mode create"
        read(textholder, *, iostat=stat) lattice_parameter
        if(stat > 0) STOP "Error reading lattice parameter"

        call get_command_argument(5, textholder, arglen)
        if(arglen==0) STOP "Esize missing in mode create"
        read(textholder, *, iostat=stat) esize
        if(stat > 0) STOP "Error reading esize"

        arg_pos = 6 
        !Check for optional keywords
        do while(.true.)
            if(arg_pos > command_argument_count()) exit
            !Pull out the next argument which should either be a keyword or an option
            call get_command_argument(arg_pos, textholder)
            textholder=adjustl(textholder)
            arg_pos=arg_pos+1

            !Choose what to based on what the option string is 
            select case(trim(textholder))

            !If orient command is passed extract the orientation to numeric array format
            case('orient')
                do i = 1, 3
                    call get_command_argument(arg_pos, orient_string, arglen)
                    if (arglen==0) STOP "Missing orientation in orient command of mode create"
                    arg_pos = arg_pos+1
                    call parse_ori_vec(orient_string, orient(i,:))
                end do



            !If the duplicate command is passed then we extract the information on the new bounds.
            case('duplicate')
                if(dim_flag) STOP "Both duplicate and dim options cannot be used in mode_create"
                dup_flag = .true.
                do i = 1, 3
                    call get_command_argument(arg_pos, textholder)
                    read(textholder, *) duplicate(i)
                    arg_pos = arg_pos + 1
                end do
            case('dim')
                if(dup_flag) STOP "Both duplicate and dim options cannot be used in mode_create"
                dim_flag = .true.
                do i = 1, 3
                    call get_command_argument(arg_pos, textholder)
                    read(textholder, *) box_len(i)
                    arg_pos = arg_pos + 1
                end do
            case('origin')
                do i = 1, 3
                    call get_command_argument(arg_pos, textholder)
                    read(textholder, *) origin(i)
                    arg_pos = arg_pos + 1
                end do

            case default
                !If it isn't an option then you have to exit
                exit
            end select
        end do

        !Calculate the lattice periodicity length in lattice units
        do i = 1, 3
            lattice_space(i) = norm2(orient(i,:))
        end do

        !Check special periodicity relations
        select case(trim(adjustl(element_type)))
        case('fcc')
            do i = 1,3
                !Check if one of the directions is 110
                if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(orient(i,3),0.0_dp))).or.&
                    (is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(orient(i,1),0.0_dp))).or.&
                    (is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(orient(i,2),0.0_dp)))) then 
                    
                    lattice_space(i) = 0.5_dp * lattice_space(i)

                !Check if one direction is 112
                else if ((is_equal(abs(orient(i,1)), abs(orient(i,2))).and.(is_equal(abs(orient(i,3)),2.0_dp*abs(orient(i,1))))).or.&
                         (is_equal(abs(orient(i,2)), abs(orient(i,3))).and.(is_equal(abs(orient(i,1)),2.0_dp*abs(orient(i,2))))).or.&
                         (is_equal(abs(orient(i,3)), abs(orient(i,1))).and.(is_equal(abs(orient(i,2)),2.0_dp*abs(orient(i,3))))))&
                          then 

                    lattice_space(i) = 0.5_dp * lattice_space(i)

                end if
            end do
        end select 
        !Now normalize the orientation matrix
        orient =  matrix_normal(orient,3)

        !Set lattice_num to 1 
        lattice_types = 1

        !If we haven't defined a basis then define the basis and add the default basis atom type and position
        if (basisnum(1) == 0) then
            max_basisnum = 1
            basisnum(1) = 1
            call add_atom_type(name, basis_type(1,1)) !If basis command not defined then we use name as the atom_name
            basis_pos(:,1) = 0.0_dp
        end if
        !
    end subroutine
    
    subroutine build_with_rhomb(box_in_lat, transform_matrix)
        !This subroutine returns all the lattice points in the box in r_lat

        !Inputs
        integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space
        real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space
        !Internal variables
        integer :: i,  inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz,  numlatpoints, ele(3,8), rzero(3), &
                   vlat(3), temp_lat(3,8), m, n, o
        real(kind=dp) :: v(3), temp_nodes(3,1,8)
        real(kind=dp), allocatable :: resize_lat_array(:,:)
        logical, allocatable :: lat_points(:,:,:)
        logical :: node_in_bd(8)

        !Do some value initialization
        max_esize = esize

        !First find the bounding lattice points (min and max points for the box in each dimension)
        numlatpoints = 1
        do i = 1, 3
            bd_in_lat(2*i-1) = minval(box_in_lat(i,:))
            bd_in_lat(2*i) = maxval(box_in_lat(i,:))
            numlatpoints = numlatpoints*(bd_in_lat(2*i)-bd_in_lat(2*i-1))
        end do

        !Allocate the correct lat variable
        select case(esize)
        !Atomistics
        case(2)
            allocate(r_atom_lat(3,numlatpoints))
        case default
            continue
        end select


        !Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case 
        !and one for the regular case
        if (esize==2) then  
            !atomistics
            do iz = bd_in_lat(5)-5, bd_in_lat(6)+5
                do iy = bd_in_lat(3)-5, bd_in_lat(4)+5
                    do ix = bd_in_lat(1)-5, bd_in_lat(2)+5
                        v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /)

                        !Transform point back to real space for easier checking
                        v = matmul(orient, matmul(transform_matrix,v))
                        !If within the boundaries
                        if(in_block_bd(v, box_bd)) then 
                            lat_atom_num = lat_atom_num + 1
                            r_atom_lat(:,lat_atom_num) =  v
                        end if
                    end do
                end do
            end do
        else
            !If we are working with elements we have to use more complex code

            !Initialize finite element 
            ele(:,:) = (esize-1) * cubic_cell(:,:)
            
            !Make a 3 dimensional array of lattice points. This array is indexed by the integer lattice position. 
            !A value of true means that the coordinate is a lattice point which is within the box_bd
            allocate(lat_points(bd_in_lat(2)-bd_in_lat(1)+10,bd_in_lat(4)-bd_in_lat(3)+10,bd_in_lat(6)-bd_in_lat(5)+10))
            lat_points(:,:,:) = .false.
            do iz = bd_in_lat(5)-5, bd_in_lat(6)+5
                do iy = bd_in_lat(3)-5, bd_in_lat(4)+5
                    do ix = bd_in_lat(1)-5, bd_in_lat(2)+5
                        v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /)

                        !Transform point back to real space for easier checking
                        v = matmul(orient, matmul(transform_matrix,v))
                        !If within the boundaries
                        if(in_block_bd(v, box_bd)) then 
                            lat_points(ix-bd_in_lat(1)+5,iy-bd_in_lat(3)+5,iz-bd_in_lat(5) + 5) = .true.
                        end if
                    end do
                end do
            end do

            !Now we redefine bd_in_lat The first 3 indices contains limits for the lat_points array
            bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10
            bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10
            bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10
            !Figure out where the starting point is. This is the first piont which fully contains the finite element
            outerloop: do iz = 1, bd_in_array(3)
                do iy = 1, bd_in_array(2)
                    do ix = 1, bd_in_array(1)
                        node_in_bd(8) = .false.
                        do inod = 1, 8
                            vlat = ele(:,inod) + (/ ix, iy, iz /)

                            !Check to see if the node resides at a position containing a lattice point within the box
                            if(any(vlat > shape(lat_points))) then 
                                continue
                            else if(lat_points(vlat(1),vlat(2),vlat(3))) then 
                                node_in_bd(inod) = .true.
                            end if
                        end do

                        if(all(node_in_bd)) then 
                            rzero = (/ ix, iy, iz /)
                            exit outerloop
                        end if
                    end do
                end do
            end do outerloop
            
            !Now build the finite element region
            lat_ele_num = 0
            lat_atom_num = 0
            allocate(r_lat(3,8,numlatpoints/esize))
            
            !Redefined the second 3 indices as the number of elements that fit in the bounds
            do i = 1, 3
                bd_in_array(3+i) = bd_in_array(i)/esize
            end do

            !Now start the element at rzero
            do inod=1, 8
                ele(:,inod) = ele(:,inod) + rzero
            end do
            do iz = -bd_in_array(6), bd_in_array(6)
                do iy = -bd_in_array(5), bd_in_array(5)
                    do ix = -bd_in_array(4), bd_in_array(4)
                        node_in_bd(:) = .false.
                        temp_nodes(:,:,:) = 0.0_dp
                        temp_lat(:,:) = 0
                        do inod = 1, 8
                            vlat= ele(:,inod) + (/ ix*(esize), iy*(esize), iz*(esize) /)
                            !Transform point back to real space for easier checking
                            ! v = matmul(orient, matmul(transform_matrix,v))
                            do i = 1,3
                                v(i) = real(vlat(i) + bd_in_lat(2*i-1) - 5)
                            end do
                            temp_nodes(:,1, inod) = matmul(orient, matmul(transform_matrix, v))
                            temp_lat(:,inod) = vlat

                            !Check to see if the lattice point values are greater then the array limits
                            if(any(vlat > shape(lat_points)).or.any(vlat < 1)) then 
                                continue
                            !If within array boundaries check to see if it is a lattice point
                            else if(lat_points(vlat(1),vlat(2),vlat(3))) then 
                                node_in_bd(inod) = .true.
                            end if
                        end do

                        if(all(node_in_bd)) then 
                            lat_ele_num = lat_ele_num+1
                            r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)

                            !Now set all the lattice points contained within an element to false
                            do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
                                do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))
                                    do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:))
                                        lat_points(m,n,o) = .false.
                                    end do 
                                end do
                            end do

                        end if
                    end do
                end do
            end do

            !Now figure out how many lattice points could not be contained in elements
            allocate(r_atom_lat(3,count(lat_points)))
            lat_atom_num = 0
            do iz  = 1, bd_in_array(3)
                do iy = 1, bd_in_array(2)
                    do ix = 1, bd_in_array(1)
                        !If this point is a lattice point then save the lattice point as an atom
                        if (lat_points(ix,iy,iz)) then 
                            v= (/ real(ix,dp), real(iy, dp), real(iz,dp) /)
                            do i = 1,3
                                v(i) = v(i) + real(bd_in_lat(2*i-1) - 5)
                            end do
                            !Transform point back to real space 
                            v = matmul(orient, matmul(transform_matrix,v))
                            lat_atom_num = lat_atom_num + 1
                            r_atom_lat(:,lat_atom_num) = v
                        end if
                    end do
                end do
            end do

        end if 

    end subroutine build_with_rhomb


    subroutine error_message(errorid)

        integer, intent(in) :: errorid

        select case(errorid)
        case(1) 
            STOP "Name is missing."
        case(2)
            print *, "Element_type is missing"
        case(3)
            print *, "Lattice_parameter is missing"
        case(4)
            print *, "Lattice parameter is not in float form"
        case(5)
            print *, "Esize is missing"
        case(6) 
            print *, "Esize is not in integer form"
        case(7)
            print *, "Cx(1) should equal Cx(2) in plane centroid finding algorithm. Please double check implementation"
        end select

        STOP 3
    end subroutine error_message


end module mode_create