First working version of model builder with several output file types and mode_create working

src/Makefile

@@ -1,24 +1,24 @@
 FC=ifort
-#FFLAGS=-c -mcmodel=large -debug -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone
+FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone
-FFLAGS=-c -mcmodel=large -Ofast
+#FFLAGS=-c -mcmodel=large -Ofast
 MODES=mode_create.o
-OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o build_subroutines.o $(MODES)
+OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o  $(MODES)
 
 .SUFFIXES:
 .SUFFIXES: .c .f .f90 .F90 .o
 
 cacmb: $(OBJECTS) 
-	$(FC) $(OBJECTS) -o $@
+	$(FC) $(FFLAGS) $(OBJECTS) -o $@
 
 .f90.o: 
-	$(FC) $(FFLAGS) $< 
+	$(FC) $(FFLAGS) -c $< 
 
 .PHONY: clean
 clean: 
 	$(RM) cacmb  *.o 
 
 testfuncs: testfuncs.o functions.o subroutines.o
-	$(FC) testfuncs.o functions.o build_subroutines.o subroutines.o elements.o -o $@
+	$(FC) testfuncs.o functions.o subroutines.o elements.o -o $@
 
 .PHONY: cleantest
 cleantest:
@@ -31,4 +31,3 @@ call_mode.o : $(MODES)
 $(MODES) io.o: atoms.o 
 $(MODES) main.o : io.o
 testfuncs.o elements.o mode_create.o: subroutines.o
-testfuncs.o : build_subroutines.o

src/call_mode.f90

@@ -12,7 +12,7 @@ subroutine call_mode(arg_num,mode)
 
     select case(mode)
     case('--create')
-        call create(arg_num)
+        call create()
 
     case default
         print *, "Mode ", mode, " currently not accepted. Please check documentation for ", &

src/elements.f90

@@ -13,13 +13,17 @@ module elements
     real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
 
     integer :: ele_num=0 !Number of elements
+    integer :: node_num=0 !Total number of nodes
 
     !Data structure used to represent atoms 
-    integer, allocatable :: tag_atom(:) !atom id 
+    integer, allocatable :: tag_atom(:), type_atom(:)!atom id 
-    character(len=100), allocatable:: type_atom(:) !atom type
     real(kind =dp),allocatable :: r_atom(:,:) !atom position
     integer :: atom_num=0 !Number of atoms
 
+    !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)
     integer :: cubic_faces(4,6)
@@ -28,12 +32,12 @@ module elements
     !We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased.
     
     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
 
     !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 :: max_basisnum, basisnum(10), basis_type(10,10)!Max basis atom number, number of basis atoms in each lattice type 
-    character(len=2) :: basis_type(10,10) !Atom type of each basis
     real(kind=dp) :: basis_pos(3,10,10) !Basis atom positions
 
     !Simulation cell parameters
@@ -150,13 +154,14 @@ module elements
         size_ele(ele_num) = size
         lat_ele(ele_num) = lat
        r_node(:,:,:,ele_num) = r(:,:,:)
+       node_num = node_num + ng_node(lat)
 
 
     end subroutine add_element
 
     subroutine add_atom(type, r)
         !Subroutine which adds an atom to the atom arrays
-        character(len=2), intent(in) :: type
+        integer, intent(in) :: type
         real(kind=dp), intent(in), dimension(3) :: r
 
         atom_num = atom_num+1
@@ -166,6 +171,32 @@ module elements
 
     end subroutine add_atom
 
+    subroutine add_atom_type(type, inttype)
+        !This subroutine adds a new atom type to the module list of atom types
+        character(len=2), intent(in) :: type
+        integer, intent(out) :: inttype
+
+        integer :: i
+        logical :: exists
+        
+        exists = .false.
+        do i=1, 10
+            if(type == type_to_name(i)) exists = .true.
+            inttype = i
+        end do 
+
+        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
@@ -184,4 +215,93 @@ module elements
             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
+        max_esize=maxval(size_ele)
+    end subroutine
+
+    subroutine interpolate_atoms(type, esize, lat_type, r_in, type_interp, r_interp)
+        !This subroutine returns the interpolated atoms from the elements. 
+
+        !Arguments
+        character(len=100), 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(3,max_basisnum, max_ng_node), intent(in) :: r_in !Nodal positions
+        integer, dimension(max_basisnum*max_esize**3), intent(out) :: type_interp !Interpolated atomic positions
+        real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: r_interp !Interpolated atomic positions
+
+        !Internal variables
+        integer :: i, it, is, ir, ibasis, inod, ia, bnum, lat_type_temp
+        real(kind=dp), allocatable :: a_shape(:)
+        real(kind=dp) :: t, s, r
+        
+        !Initialize some variables
+        r_interp(:,:) = 0.0_dp
+        type_interp(:) = 0
+        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(trim(adjustl(type)))
+        case('fcc')
+            allocate(a_shape(8))
+            !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)
+
+                            end do
+                        end do
+
+                    end do
+                end do
+            end do 
+        end select
+
+        if (ia /= 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(8)
+
+        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
 end module elements

src/functions.f90

@@ -225,7 +225,7 @@ END FUNCTION StrDnCase
         real(kind=dp), intent(in) :: A, B
         logical :: is_equal
 
-        if((A>(B - 10.0_dp**-10)).and.(A < (B+10.0_dp**-10))) then 
+        if((A>(B - 10.0_dp**(-10))).and.(A < (B+10.0_dp**(-10)))) then 
             is_equal = .true.
         else
             is_equal = .false.

src/io.f90

@@ -64,10 +64,13 @@ module io
                 outfilenum=outfilenum+1
                 outfiles(outfilenum) = temp_outfile
                 exit
-
+            case('vtk')
+                outfilenum=outfilenum+1
+                outfiles(outfilenum)=temp_outfile
+                exit
             case default    
                 print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), "not currently accepted. ", &
-                         "please input a filename with extension from following list: xyz."
+                         "please input a filename with extension from following list: xyz, lmp, vtk."
                 read(*,*) temp_outfile
 
             end select
@@ -81,6 +84,9 @@ module io
 
         integer :: i 
 
+        !Find max esize which will be needed later
+        call set_max_esize
+
         do i = 1, outfilenum
             !Pull out the extension of the file and call the correct write subroutine
             select case(trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))))
@@ -88,6 +94,8 @@ module io
                 call write_xyz(outfiles(i))
             case('lmp')
                 call write_lmp(outfiles(i))
+            case('vtk')
+                call write_vtk(outfiles(i))
             case default
                 print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
                          " is not accepted for writing. Please select from: xyz and try again"
@@ -138,10 +146,10 @@ module io
     end subroutine write_xyz
 
     subroutine write_lmp(file)
-
-        integer :: write_num, i
-        character(len=100), intent(in) :: file
         !This subroutine writes out a .lmp style dump file
+        character(len=100), intent(in) :: file
+        integer :: write_num, i, iatom, type_interp(max_basisnum*max_esize**3)
+        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), mass
 
         open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
 
@@ -150,10 +158,13 @@ module io
         write(11, '(a)')
         !Calculate total atom number
         write_num = atom_num
+        do i = 1,ele_num
+            if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3
+        end do
         !Write total number of atoms + elements
         write(11, '(i16, a)') write_num, ' atoms'
         !Write number of atom types 
-        write(11, '(i16, a)') 1, ' atom types'
+        write(11, '(i16, a)') atom_types, ' atom types'
 
         write(11,'(a)') ' '
         !Write box bd
@@ -163,15 +174,106 @@ module io
 
         !Masses
         write(11, '(a)') 'Masses'
+
         write(11, '(a)') ' '
-        write(11, '(i16, f23.15)') 1, 63.546
+        do i =1, atom_types
+            call atommass(type_to_name(i),mass)
+            write(11, '(i16, f23.15)') i, mass
+        end do
         write(11, '(a)') ' '
 
         !Write atom positions
         write(11, '(a)') 'Atoms'
         write(11, '(a)') ' '
         do i = 1, atom_num
-            write(11, '(2i16, 3f23.15)') i, 1, r_atom(:,i) 
+            write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i) 
+        end do
+
+        !Write refined element atomic positions
+        do i = 1, ele_num
+            call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp)
+            select case(trim(adjustl(type_ele(i))))
+            case('fcc')
+                do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
+                    write(11, '(2i16, 3f23.15)') atom_num+iatom, type_interp(iatom), r_interp(:,iatom)
+                end do
+            end select
         end do
     end subroutine write_lmp
+
+    subroutine write_vtk(file)
+        !This subroutine writes out a vtk style dump file
+        integer :: i, j, inod, ibasis
+        character(len=100), intent(in) :: file
+
+1       format('# vtk DataFile Version 4.0.1', / &
+              'CAC output -- cg', / &
+              'ASCII')
+11      format('# vtk DataFile Version 4.0.1', / &
+               'CACmb output -- atoms', / &
+               'ASCII')
+2       format('DATASET UNSTRUCTURED_GRID')
+3       format('POINTS', i16, ' float')
+4       format(/'CELLS', 2i16)
+5       format(/'CELL_TYPES', i16)
+12      format(/'CELL_DATA', i16)
+16      format(/'POINT_DATA', i16)
+17      format('SCALARS weight float', / &
+               'LOOKUP_TABLE default')
+18      format('SCALARS atom_type float', / &
+               'LOOKUP_TABLE default')
+
+20      format('SCALARS lattice_type float', /&
+               'LOOKUP_TABLE default')
+
+        !First we write the vtk file containing the atoms
+        open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
+
+        write(11, 11)
+        write(11, 2)
+        write(11, 3) atom_num
+        do i = 1, atom_num
+            write(11, '(3f23.15)') r_atom(:,i)
+        end do 
+        write(11,4) atom_num, atom_num*2
+        do i = 1, atom_num
+            write(11, '(2i16)') 1, i-1
+        end do
+        write(11, 5) atom_num
+        do i = 1, atom_num
+            write(11, '(i16)') 1 
+        end do
+        write(11, 16) atom_num
+        write(11, 18)
+        do i = 1, atom_num
+            write(11, '(i16)') type_atom(i)
+        end do
+        close(11)
+
+        open(unit=11, file='cg_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
+        write(11,1)
+        write(11,2)
+        write(11,3) node_num
+        do i = 1, ele_num
+            do inod=1, ng_node(lat_ele(i))  
+                do ibasis = 1, basisnum(lat_ele(i))
+                    write(11, '(3f23.1)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i))
+                end do 
+            end do 
+        end do
+        write(11, 4) ele_num, ele_num + node_num
+        do i =1, ele_num
+            write(11, '(9i16)') ng_node(lat_ele(i)), (j, j = (i-1)*ng_node(lat_ele(i)), i*ng_node(lat_ele(i))-1)
+        end do
+        write(11,5) ele_num
+        do i = 1, ele_num
+            if(trim(adjustl(type_ele(i))) == 'fcc') write(11, '(i16)') 12
+        end do
+        write(11,12) ele_num
+        write(11,20)
+        do i = 1, ele_num
+            write(11, '(i16)') lat_ele(i)
+        end do
+        close(11)
+    end subroutine
 end module io

src/mode_create.f90

@@ -5,28 +5,27 @@ module mode_create
     use atoms   
     use io
     use subroutines
+    use elements
 
     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), ox_bd(6), maxbd(3), lattice_space(3)
+                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3)
-    integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_num, lat_atom_num, bd_in_lat(6)
+    integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6)
     logical :: dup_flag, dim_flag
 
     real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
     public 
     contains
 
-    subroutine create(arg_num)
+    subroutine create()
         ! Main subroutine which controls execution
 
-        integer, intent(in) :: arg_num
         character(len=100) :: textholder
     
         integer :: i, ibasis, inod
-        real(kind=dp) :: r(3), periodvone(3), periodvtwo(3)
+        real(kind=dp), allocatable :: r_node_temp(:,:,:)
-        real(kind=dp), allocatable :: r_node(:,:,:)
 
         !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))
@@ -40,11 +39,11 @@ module mode_create
         dup_flag = .false.
         dim_flag = .false.
         basisnum = 0
-        lat_num = 0
+        lat_ele_num = 0
         lat_atom_num = 0
 
         !First we parse the command
-        call parse_command(arg_num)
+        call parse_command()
 
         ! Before building do a check on the file
         if (outfilenum == 0) then
@@ -55,7 +54,7 @@ module mode_create
         !Now we setup the unit element and call other init subroutines
         call def_ng_node(1, element_type)
 
-        allocate(r_node(3,max_basisnum,max_ng_node))
+        allocate(r_node_temp(3,max_basisnum,max_ng_node))
 
         if(dup_flag) then 
 
@@ -64,7 +63,7 @@ module mode_create
             
 
             do i = 1, 8
-                box_vert(:,i) = duplicate(:)*lattice_space(:)*cubic_cell(:,i) + origin(:) 
+                box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + origin(:) 
             end do
             call matrix_inverse(orient,3,orient_inv)
             !Now get the rotated box vertex positions in lattice space. Should be integer units
@@ -72,31 +71,14 @@ module mode_create
             !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.1_dp*lattice_space(i)
+                box_bd(2*i) = maxval(box_vert(i,:)) - 0.25_dp*lattice_space(i)
-                box_bd(2*i-1) = origin(i)
+                box_bd(2*i-1) = origin(i)-0.25_dp*lattice_space(i)
             end do
-        
             !and then call the build function with the correct transformation matrix
             select case(trim(adjustl(element_type)))
             case('fcc')
-                ! periodvone(:) = matmul(orient, reshape((/ 1, 1, 0 /),(/ 3 /)))
+
-                ! periodvtwo(:) = matmul(orient, reshape((/ 1, 1, 2 /),(/ 3 /)))
+                call build_with_rhomb(box_lat_vert, fcc_mat)
-
-                ! do i = 1, 3
-                !     if (periodvone(i) < lim_zero) then 
-                !         ! box_bd(2*i) =floor(box_bd(2*i)/periodvtwo(i))*periodvtwo(i)
-                !         box_bd(2*i) = box_bd(2*i) - 0.5*periodvtwo(i)
-                !     else if(periodvtwo(i) < lim_zero) then 
-                !         ! box_bd(2*i) =floor(box_bd(2*i)/periodvone(i))*periodvone(i)
-                !         box_bd(2*i) = box_bd(2*i) - 0.5*periodvone(i)
-                !     else    
-                !         ! box_bd(2*i) = floor(box_bd(2*i)/lcm(periodvone(i),periodvtwo(i)))*lcm(periodvone(i),periodvtwo(i))
-                !         box_bd(2*i) = box_bd(2*i) - 0.5*lcm(periodvone(i),periodvtwo(i))
-
-                !     end if
-                ! end do
-
-                call lattice_in_box(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"
@@ -116,19 +98,21 @@ module mode_create
             !Add the basis atoms to the unit cell
             do inod = 1, max_ng_node
                 do ibasis = 1, basisnum(1)
-                    r_node(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis,1) + origin(:)
+                    r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis,1) + origin(:)
                 end do
             end do
-
+            do i = 1,3
-            call add_element(element_type, esize, 1, r_node)
+                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 
             !Allocate variables 
-            call alloc_ele_arrays(lat_num, lat_atom_num*basisnum(1))
+            call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1))
             if(lat_atom_num > 0) then 
-                !Check for periodicity
                 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,1))
@@ -136,13 +120,23 @@ module mode_create
                 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,1)
+                        end do
+                    end do 
+                    call add_element(element_type, esize, 1, r_node_temp)
+                end do 
+            end if
         end if
 
     end subroutine create
     !This subroutine parses the command and pulls out information needed for mode_create
-    subroutine parse_command(arg_num)
+    subroutine parse_command()
 
-        integer, intent(in) :: arg_num
 
         integer :: arg_pos, ori_pos, i, j, arglen, stat
         character(len=100) :: textholder
@@ -217,7 +211,6 @@ module mode_create
                     read(textholder, *) origin(i)
                     arg_pos = arg_pos + 1
                 end do
-                print *, origin
             !If a filetype is passed then we add name.ext to the outfiles list
             case('xyz')
                 textholder = trim(adjustl(name)) //'.xyz'
@@ -244,13 +237,11 @@ module mode_create
         select case(trim(adjustl(element_type)))
         case('fcc')
             do i = 1,3
-                print *, orient(i,:) 
                 !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 
                     
-                    print *, '110', i
                     lattice_space(i) = 0.5_dp * lattice_space(i)
 
                 !Check if one direction is 112
@@ -259,7 +250,6 @@ module mode_create
                          (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 
 
-                    print *, '112 ', i
                     lattice_space(i) = 0.5_dp * lattice_space(i)
 
                 end if
@@ -272,21 +262,27 @@ module mode_create
         if (basisnum(1) == 0) then
             max_basisnum = 1
             basisnum(1) = 1
-            basis_type(1,1) = name !If basis command not defined then we use name as the atom_name
+            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,1) = 0.0_dp
         end if
     end subroutine
     
-    subroutine lattice_in_box(box_in_lat, transform_matrix)
+    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, j, bd_in_lat(6), ix, iy, iz,  numlatpoints
+        integer :: i,  inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz,  numlatpoints, templatpoints, ele(3,8), rzero(3), ilat, &
-        real(kind=dp) :: box_face_center(3,6), face_normals(3,6), Cx(2), Cy, Cz, A(2), v(3)
+                   type_interp(basisnum(1)*esize**3), vlat(3), temp_lat(3,8), m, n, o
+        real(kind=dp) :: v(3), temp_nodes(3,1,8), ele_atoms(3,esize**3), r_interp(3,basisnum(1)*esize**3)
         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
@@ -304,82 +300,10 @@ module mode_create
         case default
             continue
         end select
-        !Calculate the box_face centroids and box face normals. This is used in the centroid code.
-        box_face_center(:,:) = 0.0_dp
-        face_normals = reshape((/ -1.0_dp,  0.0_dp,  0.0_dp, &
-                                   1.0_dp,  0.0_dp,  0.0_dp, &
-                                   0.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, &
-                                   0.0_dp,  0.0_dp,  1.0_dp /),&
-                                   shape(face_normals))
-        !Face normals
-        select case(trim(adjustl(element_type)))
-        case('fcc')
-            do i = 1, 6
-                !Box face normal
-                face_normals(:,i) = matmul(fcc_inv, matmul(orient_inv, face_normals(:,i)))
-            end do
-        end select
-
-        !Face centroids
-        do i =1, 6
-
-            !Initialize variables
-            A(:) = 0.0_dp
-            Cx(:) = 0.0_dp
-            Cy = 0.0_dp
-            Cz = 0.0_dp
-
-            !Calculate all terms using a projection onto the xy and xz planes and then using the 2d equation 
-            !for centroid of a plane. This is why we calculate the x centroid twice.
-            do j = 1, 4
-                ! A(1) = A(1) + 0.5*(box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) &
-                !                     - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i)))
-
-                ! !Centroid in x from xy plane
-                ! Cx(1) = Cx(1) + (box_in_lat(1,cubic_faces(j,i))+box_in_lat(1,cubic_faces(j+1,i)))* &
-                !                     (box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) &
-                !                     - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i)))
-                
-                ! !Centroid in y from xy plane
-                ! Cy    = Cy    + (box_in_lat(2,cubic_faces(j,i))+box_in_lat(2,cubic_faces(j+1,i)))* &
-                !                     (box_in_lat(1,cubic_faces(j,i))*box_in_lat(2,cubic_faces(j+1,i)) &
-                !                     - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(2,cubic_faces(j,i)))
-
-                ! A(2) = A(2) + 0.5*(box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) &
-                !                     - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i)))
-
-                ! !Centroid in x from xz plane
-                ! Cx(2) = Cx(2) + (box_in_lat(1,cubic_faces(j,i))+box_in_lat(1,cubic_faces(j+1,i)))* &
-                !                     (box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) &
-                !                     - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i)))
-
-                ! !Centroid in z from xz plane
-                ! Cz    = Cz    + (box_in_lat(3,cubic_faces(j,i))+box_in_lat(3,cubic_faces(j+1,i)))* &
-                !                     (box_in_lat(1,cubic_faces(j,i))*box_in_lat(3,cubic_faces(j+1,i)) &
-                !                     - box_in_lat(1,cubic_faces(j+1,i))*box_in_lat(3,cubic_faces(j,i)))
-
-                ! print *, j, i, Cx, Cy, Cz, A
-                Cx(1) = Cx(1) + box_in_lat(1,cubic_faces(j,i))*0.25
-                Cy = Cy + box_in_lat(2,cubic_faces(j,i))*0.25
-                Cz = Cz + box_in_lat(3,cubic_faces(j,i))*0.25
-
-            end do
 
-            ! Cx = Cx * 1/(6*A)
-            ! if(Cx(1) /= Cx(2)) then 
-            !     call error_message(7)
-            ! end if
-            ! Cy = Cy* 1/(6*A(1))
-            ! Cz = Cz*1/(6*A(2))
-
-            box_face_center(:,i) = (/ Cx(1), Cy, Cz /)
-        end do
 
         !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
-        print *, box_bd
         if (esize==2) then  
             !atomistics
             do iz = bd_in_lat(5)-5, bd_in_lat(6)+5
@@ -398,10 +322,141 @@ module mode_create
                 end do
             end do
         else
-            continue 
+            !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
+            print *, count(lat_points)
+            allocate(r_atom_lat(3,count(lat_points)))
+            lat_atom_num = 0
+            do ix  = 1, bd_in_array(3)
+                do iy = 1, bd_in_array(2)
+                    do iz = 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
+
+            print *, lat_atom_num
         end if 
 
-    end subroutine lattice_in_box
+    end subroutine build_with_rhomb
 
 
     subroutine error_message(errorid)