Latest updates

src/Makefile

@@ -33,7 +33,7 @@ $(OBJDIR)/%.o: %.f90
 
 
 deps:
-	@fortdepend -o Makefile.dep -i mpi -b obj -w
+	 @makedepf90 -b obj *.f90 > Makefile.dep
 
 #----------------- CLEAN UP -----------------#
 

src/Makefile.dep

@@ -1,26 +1,190 @@
-obj/atoms.o : atoms.f90 obj/functions.o obj/parameters.o 
+# This file is generated automatically. DO NOT EDIT!
-obj/box.o : box.f90 obj/functions.o obj/parameters.o 
+
-obj/caller.o : caller.f90 obj/box.o obj/opt_bubble.o obj/opt_slip_plane.o obj/opt_redef_box.o obj/opt_delete.o obj/opt_deform.o obj/opt_orient.o obj/opt_group.o obj/opt_disl.o obj/parameters.o obj/mode_calc.o obj/mode_metric.o obj/mode_merge.o obj/mode_convert.o obj/mode_create.o obj/mode_da.o 
+obj/main :  \
-obj/elements.o : elements.f90 obj/box.o obj/sorts.o obj/subroutines.o obj/functions.o obj/parameters.o 
+	obj/atoms.o \
-obj/functions.o : functions.f90 obj/parameters.o 
+	obj/box.o \
-obj/io.o : io.f90 obj/str.o obj/box.o obj/atoms.o obj/parameters.o obj/elements.o 
+	obj/caller.o \
-obj/main.o : main.f90 obj/caller.o obj/io.o obj/elements.o obj/parameters.o 
+	obj/elements.o \
-obj/mode_calc.o : mode_calc.f90 obj/box.o obj/elements.o obj/subroutines.o obj/io.o obj/parameters.o 
+	obj/functions.o \
-obj/mode_convert.o : mode_convert.f90 obj/io.o obj/elements.o obj/box.o obj/parameters.o 
+	obj/io.o \
-obj/mode_create.o : mode_create.f90 obj/box.o obj/elements.o obj/subroutines.o obj/io.o obj/atoms.o obj/parameters.o 
+	obj/main.o \
-obj/mode_da.o : mode_da.f90 obj/neighbors.o obj/elements.o obj/io.o obj/parameters.o 
+	obj/mode_calc.o \
-obj/mode_merge.o : mode_merge.f90 obj/neighbors.o obj/elements.o obj/subroutines.o obj/io.o obj/atoms.o obj/parameters.o 
+	obj/mode_convert.o \
-obj/mode_metric.o : mode_metric.f90 obj/neighbors.o obj/elements.o obj/io.o obj/parameters.o 
+	obj/mode_create.o \
-obj/neighbors.o : neighbors.f90 obj/functions.o obj/subroutines.o obj/elements.o obj/parameters.o 
+	obj/mode_da.o \
-obj/opt_bubble.o : opt_bubble.f90 obj/opt_group.o obj/box.o obj/elements.o obj/parameters.o 
+	obj/mode_merge.o \
-obj/opt_deform.o : opt_deform.f90 obj/box.o obj/elements.o obj/subroutines.o obj/parameters.o 
+	obj/mode_metric.o \
-obj/opt_delete.o : opt_delete.f90 obj/neighbors.o obj/elements.o obj/subroutines.o obj/parameters.o 
+	obj/neighbors.o \
-obj/opt_disl.o : opt_disl.f90 obj/box.o obj/subroutines.o obj/elements.o obj/parameters.o 
+	obj/opt_bubble.o \
-obj/opt_group.o : opt_group.f90 obj/sorts.o obj/box.o obj/subroutines.o obj/elements.o obj/parameters.o 
+	obj/opt_deform.o \
-obj/opt_orient.o : opt_orient.f90 obj/box.o obj/elements.o obj/subroutines.o obj/parameters.o 
+	obj/opt_delete.o \
-obj/opt_redef_box.o : opt_redef_box.f90 obj/subroutines.o obj/elements.o obj/box.o 
+	obj/opt_disl.o \
-obj/opt_slip_plane.o : opt_slip_plane.f90 obj/subroutines.o obj/functions.o obj/elements.o obj/parameters.o 
+	obj/opt_group.o \
-obj/parameters.o : parameters.f90 
+	obj/opt_orient.o \
-obj/sorts.o : sorts.f90 obj/parameters.o 
+	obj/opt_redef_box.o \
-obj/str.o : str.f90 
+	obj/opt_slip_plane.o \
-obj/subroutines.o : subroutines.f90 obj/str.o obj/box.o obj/functions.o obj/parameters.o 
+	obj/parameters.o \
+	obj/sorts.o \
+	obj/str.o \
+	obj/subroutines.o
+
+obj/atoms.o :  \
+	obj/functions.o \
+	obj/parameters.o
+
+obj/box.o :  \
+	obj/functions.o \
+	obj/parameters.o
+
+obj/caller.o :  \
+	obj/box.o \
+	obj/mode_calc.o \
+	obj/mode_convert.o \
+	obj/mode_create.o \
+	obj/mode_da.o \
+	obj/mode_merge.o \
+	obj/mode_metric.o \
+	obj/opt_bubble.o \
+	obj/opt_deform.o \
+	obj/opt_delete.o \
+	obj/opt_disl.o \
+	obj/opt_group.o \
+	obj/opt_orient.o \
+	obj/opt_redef_box.o \
+	obj/opt_slip_plane.o \
+	obj/parameters.o
+
+obj/elements.o :  \
+	obj/atoms.o \
+	obj/box.o \
+	obj/functions.o \
+	obj/parameters.o \
+	obj/sorts.o \
+	obj/subroutines.o
+
+obj/functions.o :  \
+	obj/parameters.o
+
+obj/io.o :  \
+	obj/atoms.o \
+	obj/box.o \
+	obj/elements.o \
+	obj/parameters.o \
+	obj/str.o
+
+obj/main.o :  \
+	obj/caller.o \
+	obj/elements.o \
+	obj/io.o \
+	obj/parameters.o
+
+obj/mode_calc.o :  \
+	obj/box.o \
+	obj/elements.o \
+	obj/io.o \
+	obj/parameters.o \
+	obj/subroutines.o
+
+obj/mode_convert.o :  \
+	obj/box.o \
+	obj/elements.o \
+	obj/io.o \
+	obj/parameters.o
+
+obj/mode_create.o :  \
+	obj/atoms.o \
+	obj/box.o \
+	obj/elements.o \
+	obj/io.o \
+	obj/parameters.o \
+	obj/subroutines.o
+
+obj/mode_da.o :  \
+	obj/elements.o \
+	obj/io.o \
+	obj/neighbors.o \
+	obj/parameters.o
+
+obj/mode_merge.o :  \
+	obj/atoms.o \
+	obj/elements.o \
+	obj/io.o \
+	obj/neighbors.o \
+	obj/parameters.o \
+	obj/subroutines.o
+
+obj/mode_metric.o :  \
+	obj/elements.o \
+	obj/io.o \
+	obj/neighbors.o \
+	obj/parameters.o
+
+obj/neighbors.o :  \
+	obj/box.o \
+	obj/elements.o \
+	obj/functions.o \
+	obj/parameters.o \
+	obj/subroutines.o
+
+obj/opt_bubble.o :  \
+	obj/atoms.o \
+	obj/box.o \
+	obj/elements.o \
+	obj/opt_group.o \
+	obj/parameters.o
+
+obj/opt_deform.o :  \
+	obj/box.o \
+	obj/elements.o \
+	obj/parameters.o \
+	obj/subroutines.o
+
+obj/opt_delete.o :  \
+	obj/elements.o \
+	obj/neighbors.o \
+	obj/parameters.o \
+	obj/subroutines.o
+
+obj/opt_disl.o :  \
+	obj/box.o \
+	obj/elements.o \
+	obj/parameters.o \
+	obj/subroutines.o
+
+obj/opt_group.o :  \
+	obj/atoms.o \
+	obj/box.o \
+	obj/elements.o \
+	obj/parameters.o \
+	obj/sorts.o \
+	obj/subroutines.o
+
+obj/opt_orient.o :  \
+	obj/box.o \
+	obj/elements.o \
+	obj/parameters.o \
+	obj/subroutines.o
+
+obj/opt_redef_box.o :  \
+	obj/box.o \
+	obj/elements.o \
+	obj/subroutines.o
+
+obj/opt_slip_plane.o :  \
+	obj/elements.o \
+	obj/functions.o \
+	obj/parameters.o \
+	obj/subroutines.o
+
+obj/parameters.o : 
+
+obj/sorts.o :  \
+	obj/parameters.o
+
+obj/str.o : 
+
+obj/subroutines.o :  \
+	obj/box.o \
+	obj/functions.o \
+	obj/parameters.o \
+	obj/str.o

src/caller.f90

@@ -60,6 +60,8 @@ module caller
             call group(arg_pos)
         case('-ow')
             arg_pos = arg_pos + 1
+        case('-novel')
+            arg_pos = arg_pos+1
         case('-wrap')
             arg_pos = arg_pos + 1
         case('-norefine')

src/elements.f90

@@ -2,6 +2,7 @@ module elements
     
     !This module contains the elements data structures, structures needed for building regions
     !and operations that are done on elements
+    use atoms
     use parameters
     use functions
     use subroutines
@@ -31,20 +32,21 @@ module elements
     !Mapping atom type to provided name
     character(len=2), dimension(10) :: type_to_name
     integer :: atom_types = 0
+    real(kind=dp) :: masses(10)
 
     !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)
+    integer :: cubic_faces(4,6), oneonetwopairs(2,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)
+    real(kind=dp) :: lapa(10)
 
-    !These variables contain information on the basis, for simplicities sake we limit 
+    !These variables contain informatione 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 
@@ -104,6 +106,14 @@ module elements
         cubic_faces(:,5) = (/ 4, 3, 7, 8 /)
         cubic_faces(:,6) = (/ 5, 6, 7, 8 /)
         
+        !Now mark which node pairs represent the 112 directions. This is used for the dislocation analysis.
+        oneonetwopairs(:,1) = (/ 1, 3 /)
+        oneonetwopairs(:,2) = (/ 1, 6 /)
+        oneonetwopairs(:,3) = (/ 2, 7 /)
+        oneonetwopairs(:,4) = (/ 1, 8 /)
+        oneonetwopairs(:,5) = (/ 4, 7 /)
+        oneonetwopairs(:,6) = (/ 5, 7 /)
+
         !!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, &
@@ -345,9 +355,9 @@ module elements
 
     end subroutine add_atom
 
-    subroutine add_atom_type(type, inttype, force_new)
+    subroutine add_atom_type(mass, inttype, force_new)
         !This subroutine adds a new atom type to the module list of atom types
-        character(len=2), intent(in) :: type
+        real(kind=dp), intent(in) :: mass
         integer, intent(out) :: inttype
 
         logical, optional, intent(in) :: force_new
@@ -364,7 +374,7 @@ module elements
         exists = .false.
         if(.not.force) then 
             do i=1, 10
-                if(type == type_to_name(i)) then 
+                if(is_equal(mass, masses(i))) then 
                     exists = .true.
                     inttype = i
                     exit
@@ -378,7 +388,8 @@ module elements
                 print *, "Defined atom types are greater than 10 which is currently not supported."
                 stop 3
             end if
-            type_to_name(atom_types) = type
+            call atommassspecies(mass, type_to_name(atom_types))
+            masses(atom_types) = mass
             inttype = atom_types
         end if
         return
@@ -426,10 +437,10 @@ module elements
         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
+        real(kind=dp), dimension(:,:), optional,intent(out) :: data_interp !Interpolated atomic positions
 
         !Internal variables
-        integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp
+        integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp, adjust
         real(kind=dp) :: a_shape(max_ng_node)
         real(kind=dp) :: t, s, r
         
@@ -451,16 +462,15 @@ module elements
             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)
+                t=-1.0_dp +(it-1)*(2.0_dp/(esize-1))   
                 do is =1, esize
-                    s = (1.0_dp*(is-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
+                    s=-1.0_dp +(is-1)*(2.0_dp/(esize-1))   
                     do ir = 1, esize
-                        r = (1.0_dp*(ir-1) - (esize-1)/2)/(1.0_dp*(esize-1)/2)
+                        r=-1.0_dp +(ir-1)*(2.0_dp/(esize-1))   
                         call rhombshape(r,s,t,a_shape)
 
                         do ibasis = 1, bnum
@@ -518,6 +528,68 @@ module elements
         return
     end subroutine interpolate_atoms
 
+    subroutine interpolate_vel(type, esize, lat_type, vel_in, vel_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) :: vel_in !Nodal positions
+        real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: vel_interp !Interpolated atomic positions
+
+        !Internal variables
+        integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp, adjust
+        real(kind=dp) :: a_shape(max_ng_node)
+        real(kind=dp) :: t, s, r
+        
+        !Initialize some variables
+        vel_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)*(2.0_dp/(esize-1))   
+                do is =1, esize
+                    s=-1.0_dp +(is-1)*(2.0_dp/(esize-1))   
+                    do ir = 1, esize
+                        r=-1.0_dp +(ir-1)*(2.0_dp/(esize-1))   
+                        call rhombshape(r,s,t,a_shape)
+
+                        do ibasis = 1, bnum
+                            ia = ia + 1
+                            do inod = 1, 8
+                                vel_interp(:,ia) = vel_interp(:,ia) + a_shape(inod) * vel_in(:,ibasis,inod)
+
+                            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_vel
+
     subroutine rhombshape(r,s,t, shape_fun)
         !Shape function for rhombohedral elements
         real(kind=8), intent(in) :: r, s, t
@@ -651,6 +723,7 @@ module elements
         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)
@@ -1008,11 +1081,36 @@ do i = 1, atom_num
     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)
+        real(kind=dp), intent(in) :: eng, force(:), virial(:)
+
+        integer :: size_atom_array
+        real(kind=dp), allocatable :: tmp_eng(:), tmp_force(:,:), tmp_virial(:,:)
+
+        size_atom_array=size(tag_atom)
+        if(ia > atom_num) then 
+            stop "ia in add_atom_dat shouldn't be larger than atom_num"
+        else if(ia > size(energy_atom)) then
+            allocate(tmp_eng(size(energy_atom)+1024))
+            allocate(tmp_force(3, size(energy_atom)+1024)) 
+            allocate(tmp_virial(6, size(energy_atom)+1024))
+
+            tmp_eng=0.0_dp
+            tmp_eng(1:size(energy_atom)) = energy_atom
+            call move_alloc(tmp_eng, energy_atom)
+
+            tmp_force=0.0_dp
+            tmp_force(:, 1:size(force_atom, 2)) = force_atom
+            call move_alloc(tmp_force, force_atom)
+
+            tmp_virial=0.0_dp
+            tmp_virial(:, 1:size(virial_atom,2)) = virial_atom
+            call move_alloc(tmp_virial, virial_atom)
+        end if
 
         energy_atom(ia) = eng
         force_atom(:,ia) = force(:) 
         virial_atom(:,ia) = virial(:) 
+
         vflag = .true.
         return
     end subroutine add_atom_data

src/functions.f90

@@ -168,11 +168,11 @@ END FUNCTION StrDnCase
 
         do i =1 ,3
             !Check upper bound
-            if(v(i) > (box_bd(2*i)+10.0_dp**(-10)) ) then 
+            if(v(i) > (box_bd(2*i))) then 
                 in_block_bd =.false.
                 exit
             !Check lower bound
-            else if (v(i) < (box_bd(2*i-1)-10.0_dp**(-10))) then 
+            else if (v(i) < (box_bd(2*i-1))+1d-6) then 
                 in_block_bd = .false.
                 exit
             end if

src/io.f90

@@ -139,7 +139,7 @@ module io
         do i = 1, ele_num
             do inod = 1, ng_node(lat_ele(i))
                 do ibasis = 1, basisnum(lat_ele(i))
-                        write(11, '(2i16, 3f23.15)') basis_type(ibasis,lat_ele(i)), 1, r_node(:,ibasis,inod,i)
+                        write(11, *) basis_type(ibasis,lat_ele(i)), 1, r_node(:,ibasis,inod,i)
                         outn = outn + 1
                 end do 
             end do 
@@ -151,7 +151,7 @@ module io
 
         !Write atom positions
         do i = 1, atom_num
-            write(11, '(2i16, 3f23.15)') type_atom(i), 0, r_atom(:,i) 
+            write(11, *) type_atom(i), 0, r_atom(:,i) 
             outn = outn + 1
         end do
 
@@ -207,7 +207,7 @@ module io
         write(11, '(a)') 'Atoms'
         write(11, '(a)') ' '
         do i = 1, atom_num
-            write(11, '(2i16, 3f23.15)') tag_atom(i), type_atom(i), r_atom(:,i) 
+            write(11, *) tag_atom(i), type_atom(i), r_atom(:,i) 
         end do
 
         !Write refined element atomic positions
@@ -219,7 +219,7 @@ module io
                 do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
                     interp_num = interp_num+1
                     call apply_periodic(r_interp(:,iatom))
-                    write(11, '(2i16, 3f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom)
+                    write(11, *) atom_num+interp_num, type_interp(iatom), r_interp(:,iatom)
                 end do
             end select
         end do
@@ -277,12 +277,12 @@ module io
         
         if (write_dat) then 
             do i = 1, atom_num
-                write(11, '(2i16, 13f23.15)') tag_atom(i), type_atom(i), r_atom(:,i), energy_atom(i), &
+                write(11, *) tag_atom(i), type_atom(i), r_atom(:,i), energy_atom(i), &
                                               force_atom(:,i), virial_atom(:,i)
             end do
         else
             do i = 1, atom_num
-                write(11, '(2i16, 3f23.15)') tag_atom(i), type_atom(i), r_atom(:,i)
+                write(11, *) tag_atom(i), type_atom(i), r_atom(:,i)
             end do
         end if
 
@@ -294,12 +294,12 @@ module io
                 do inod =1, ng_node(lat_ele(i))
                     do ibasis=1, basisnum(lat_ele(i))
                         if(write_dat) then 
-                            write(11, '(2i16, 13f23.15)') atom_num+interp_num, basis_type(ibasis,lat_ele(i)), &
+                            write(11, *) atom_num+interp_num, basis_type(ibasis,lat_ele(i)), &
                                         r_node(:, ibasis, inod, i), energy_node(ibasis,inod,i), force_node(:, ibasis, inod, i), &
                                         virial_node(:, ibasis, inod, i)
                         else
     
-                            write(11, '(2i16, 3f23.15)') atom_num+interp_num, basis_type(ibasis,lat_ele(i)), &
+                            write(11, *) atom_num+interp_num, basis_type(ibasis,lat_ele(i)), &
                                                          r_node(:, ibasis, inod, i) 
                         end if
                         interp_num = interp_num +1
@@ -321,14 +321,14 @@ module io
                         do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
                             interp_num = interp_num+1
                             call apply_periodic(r_interp(:,iatom))
-                            write(11, '(2i16, 13f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom), &
+                            write(11, *) atom_num+interp_num, type_interp(iatom), r_interp(:,iatom), &
                                                           data_interp(:,iatom)
                         end do
                     else
                         do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
                             interp_num = interp_num+1
                             call apply_periodic(r_interp(:,iatom))
-                            write(11, '(2i16, 3f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom)
+                            write(11, *) atom_num+interp_num, type_interp(iatom), r_interp(:,iatom)
                         end do
                     end if
                 end select
@@ -436,6 +436,7 @@ module io
 16      format(/'POINT_DATA', i16)
 17      format('SCALARS weight float', / &
                'LOOKUP_TABLE default')
+
 18      format('SCALARS atom_type float', / &
                'LOOKUP_TABLE default')
 
@@ -445,6 +446,9 @@ module io
 21      format('SCALARS esize float', /&
                'LOOKUP_TABLE default')
 
+22      format('SCALARS energy 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')
 
@@ -467,6 +471,14 @@ module io
         do i = 1, atom_num
             write(11, '(i16)') type_atom(i)
         end do
+
+        !Write the atom data 
+        if(allocated(force_atom)) then
+            write(11,22)
+            do i =1, atom_num
+                write(11, "(f23.15)") energy_atom(i)
+            end do
+        end if
         close(11)
 
         !Now we write the vtk file for the elements
@@ -499,6 +511,18 @@ module io
         do i = 1, ele_num
             write(11, '(i16)') size_ele(i)
         end do
+        
+        if(allocated(force_node)) then 
+            write(11,16) node_num
+            write(11,22)
+            do i = 1, ele_num
+                do inod=1, ng_node(lat_ele(i))
+                    do ibasis=1, basisnum(lat_ele(i))
+                        write(11,'(f23.15)') energy_node(ibasis, inod, i)
+                    end do
+                end do
+            end do
+        end if
         close(11)
     end subroutine
 
@@ -509,7 +533,7 @@ module io
         character(len=100), intent(in) :: file
         integer :: interp_max, i, j, inod, ibasis, ip, unique_index(50), unique_size(50), unique_type(50), unique_num, &
                    etype
-        real(kind=dp) :: box_vec(3), masses(10)
+        real(kind=dp) :: box_vec(3)
 
 1  format('time' / i16, f23.15)
 2  format('number of elements' / i16)
@@ -548,9 +572,6 @@ module io
         write(11,3) node_num
         write(11,19) max_ng_node, max_basisnum
         write(11,5) atom_num
-        do i = 1, atom_types
-            call atommass(type_to_name(i),masses(i))
-        end do
         write(11,*) "masses: " 
         write(11, *) (masses(i), i = 1, atom_types)
         write(11,7) box_bc(1:1), box_bc(2:2), box_bc(3:3)
@@ -603,12 +624,12 @@ module io
             if(allocated(vel_atom)) then 
                 write(11, 16)
                 do i = 1, atom_num
-                    write(11, '(2i16, 6f23.15)') i, type_atom(i), r_atom(:,i), vel_atom(:, i)
+                    write(11, '(2i16, 6f23.15)') tag_atom(i), type_atom(i), r_atom(:,i), vel_atom(:, i)
                 end do
             else
                 write(11,14)
                 do i = 1, atom_num
-                    write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i)
+                    write(11, '(2i16, 3f23.15)') tag_atom(i), type_atom(i), r_atom(:,i)
                 end do
             end if
         end if
@@ -631,7 +652,7 @@ module io
         write(11,*) box_bd(:)
 
         !Write the number of atom types in the current model and all of their names
-        write(11,*) atom_types, (type_to_name(i)//' ', i=1, atom_types)
+        write(11,*) atom_types, (masses(i), i=1, atom_types)
         !Write the number of lattice_types, basisnum and number of nodes for each lattice type
         write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types)
         !Now for every lattice type write the basis atom types
@@ -754,9 +775,8 @@ module io
                    new_type_to_type(10), new_lattice_types, new_lattice_map(10), btype_map(10,10), bnum_map(10), &
                    node_map(10)
         character(len=100) :: etype
-        real(kind=dp) ::  r(3), newdisplace(3), lapa_map(10)
+        real(kind=dp) ::  r(3), newdisplace(3), lapa_map(10), new_masses(10)
         real(kind=dp), allocatable :: r_innode(:,:,:)
-        character(len = 2) :: new_type_to_name(10)
         !First open the file
         open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
 
@@ -777,11 +797,11 @@ module io
         call grow_box(temp_box_bd)
 
         !Read in the number of atom types and all their names
-        read(11, *) new_atom_types, (new_type_to_name(i), i = 1, new_atom_types)
+        read(11, *) new_atom_types, (new_masses(i), i = 1, new_atom_types)
         !Now fit these into the global list of atom types, after this new_type_to_type is the actual global 
         !type of the atoms within this file
         do i = 1, new_atom_types
-            call add_atom_type(new_type_to_name(i), new_type_to_type(i), all_new)
+            call add_atom_type(new_masses(i), new_type_to_type(i), all_new)
         end do
         !Read the number of lattice types, basisnum, and number of nodes for each lattice type
         read(11,*) new_lattice_types, (bnum_map(i), i = 1, new_lattice_types),(node_map(i), i =1, new_lattice_types)
@@ -893,8 +913,7 @@ module io
         !Read define atom_types by mass
         print *, j
         do i = 1, j
-            call realatomspecies(atomic_masses(i), atomic_element)
+            call add_atom_type(atomic_masses(i), atom_type_map(i), all_new)
-            call add_atom_type(atomic_element, atom_type_map(i), all_new)
             print *, i, atom_type_map(i)
         end do 
 
@@ -979,7 +998,7 @@ module io
 
         !Internal Variables
         integer :: i, j, in_eles, in_atoms, inbtypes(10), lat_type, ia, ie, inod, &
-            id, type_node, ilat, esize, tag, type, bnum, n, ibasis, ip, atom_type_map(100)
+            id, type_node, ilat, esize, tag, atype, bnum, n, ibasis, ip, atom_type_map(100)
         real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), vela(3),&
                          ee(10,20), fe(3,10,20), ve(6,10,20), re(3,10,20), vele(3,10,20), atomic_masses(10)
         character(len=100) :: textholder, fcc
@@ -1029,8 +1048,7 @@ module io
 
         !Read define atom_types by mass
         do i = 1, j
-            call realatomspecies(atomic_masses(i), atomic_element)
+            call add_atom_type(atomic_masses(i), atom_type_map(i), all_new)
-            call add_atom_type(atomic_element, atom_type_map(i), all_new)
         end do 
 
         if(in_atoms > 0 ) then
@@ -1046,11 +1064,11 @@ module io
             do ia = 1, in_atoms
                 read(11,'(a)')  line(:)
                 if(read_vel) then 
-                    read(line,*) tag, type, ra(:), ea, fa(:), va(:), vela(:)         
+                    read(line,*) tag, atype, ra(:), ea, fa(:), va(:), vela(:)         
                 else
-                    read(line,*) tag, type, ra(:), ea, fa(:), va(:)         
+                    read(line,*) tag, atype, ra(:), ea, fa(:), va(:)         
                 end if
-                call add_atom(tag, atom_type_map(type), ra)
+                call add_atom(tag, atom_type_map(atype), ra)
                 call add_atom_data(atom_num, ea, fa, va)
                 if(read_vel) vel_atom(:,atom_num) = vela 
             end do
@@ -1152,8 +1170,7 @@ module io
         !Read atomic masses
         do i = 1, type_in
             read(11,*) j, mass
-            call realatomspecies(mass, atom_species) 
+            call add_atom_type(mass, type_map(i), all_new)
-            call add_atom_type(atom_species, type_map(i), all_new)
         end do
 
         !Read useless info
@@ -1226,8 +1243,7 @@ module io
         !Read atomic masses
         do i = 1, type_in
             read(11,*) j, mass, textholder
-            call realatomspecies(mass, atom_species) 
+            call add_atom_type(mass, type_map(i), all_new)
-            call add_atom_type(atom_species, type_map(i), all_new)
         end do
 
         !Read useless info
@@ -1311,6 +1327,7 @@ module io
         integer :: i, j,arglen, arg_pos, ntypes
 
         character(len=100) :: textholder
+        real(kind=dp) :: mass
 
         arg_pos = apos + 1
         call get_command_argument(arg_pos, textholder, arglen)
@@ -1320,7 +1337,8 @@ module io
         do i=1,ntypes
             arg_pos = arg_pos + 1
             call get_command_argument(arg_pos, textholder, arglen)
-            call add_atom_type(textholder, j)
+            call atommass(textholder, mass)
+            call add_atom_type(mass, j)
         end do
 
         return

src/mode_calc.f90

@@ -25,6 +25,9 @@ module mode_calc
         case('tot_virial')
             allocate(calculated(6))
             call calc_tot_virial
+        case('tot_energy')
+            allocate(calculated(1))
+            call calc_tot_energy
         case default
             print *, trim(adjustl(calc_opt)), " is not accepted as a calc option in mode_calc"
             stop 3
@@ -87,5 +90,26 @@ module mode_calc
 
         print *, "Total virial is calculated as : (v11, v22, v33, v32, v31, v21)"
         print *, calculated
-    end subroutine 
+    end subroutine calc_tot_virial
+
+    subroutine calc_tot_energy
+        integer :: i, inod, ibasis, j
+        calculated=0
+        !Sum atom energies
+        do i = 1, atom_num
+            calculated(1) = calculated(1) + energy_atom(i)
+        end do
+
+        !Sum element energies
+        do i =1, ele_num
+            do inod=1, ng_node(lat_ele(i))
+                do ibasis=1, basisnum(lat_ele(i))
+                    calculated(1) = calculated(1) + energy_node(ibasis, inod, i)*(size_ele(i)**3)/ng_node(lat_ele(i))
+                end do
+            end do
+        end do
+
+        print *, 'Total energy in eV is:'
+        print *, calculated(1)
+    end subroutine calc_tot_energy
 end module mode_calc

src/mode_create.f90

@@ -12,9 +12,9 @@ module mode_create
 
     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)
+                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3),  basis_pos(3,10)
     integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
-               basis_pos(3,10), esize_nums, esize_index(10)
+               esize_nums, esize_index(10)
     logical :: dup_flag, dim_flag, efill, crossb(3)
 
     real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
@@ -186,6 +186,7 @@ module mode_create
         character(len=100) :: orient_string
         character(len=2) :: btype
         logical :: isortho, isrighthanded, bool
+        real(kind=dp) :: mass
 
 
         !Pull out all required positional arguments
@@ -265,7 +266,8 @@ module mode_create
                 do i = 1, max_basisnum
                     call get_command_argument(arg_pos, btype)
                     arg_pos = arg_pos+1
-                    call add_atom_type(btype, basis_type(i,1))
+                    call atommass(btype, mass)
+                    call add_atom_type(mass, basis_type(i,1))
                     do j = 1, 3
                         call get_command_argument(arg_pos, textholder)
                         read(textholder, *) basis_pos(j,i)
@@ -334,7 +336,8 @@ module mode_create
         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
+            call atommass(name, mass)
+            call add_atom_type(mass, basis_type(1,1)) !If basis command not defined then we use name as the atom_name
             basis_pos(:,1) = 0.0_dp
         end if
         
@@ -371,6 +374,8 @@ module mode_create
             numlatpoints = numlatpoints*(bd_in_lat(2*i)-bd_in_lat(2*i-1))
         end do
 
+        if(numlatpoints < 0) stop "number of lattice points is negative, this can occur if the model is too big"
+
         !Allocate the correct lat variable
         select case(esize)
         !Atomistics

src/mode_da.f90

@@ -8,7 +8,9 @@ module mode_da
     
     implicit none 
 
+    integer :: num_d_points
     integer, allocatable :: is_slipped(:,:,:) 
+    real(kind=dp), allocatable :: disl_line(:,:), disl_data(:)
 
     private :: write_xyz
     public 
@@ -31,16 +33,28 @@ module mode_da
         if(allocated(is_slipped)) deallocate(is_slipped)
         allocate(is_slipped(max_basisnum, max_ng_node, ele_num))
 
+        if(allocated(disl_line)) deallocate(disl_line)
+        num_d_points=0
+        allocate(disl_line(3,1024))
+        allocate(disl_data(1024))
+
         call cga_da
 
         !Now create the output file num and write out to xyz format
+        ppos = scan(trim(infiles(1)),".", BACK= .true.)
+        if ( ppos > 0 ) then 
+            outfile = 'da_'//infiles(1)(1:ppos)//'vtk'
+        else 
+            outfile = 'da_'//infiles(1)//'.vtk'
+        end if
+        call write_da_vtk(outfile)        
+
         ppos = scan(trim(infiles(1)),".", BACK= .true.)
         if ( ppos > 0 ) then 
             outfile = 'da_'//infiles(1)(1:ppos)//'xyz'
         else 
             outfile = 'da_'//infiles(1)//'.xyz'
         end if
-
         call write_da_xyz(outfile)
         
         return
@@ -65,16 +79,39 @@ module mode_da
 
     subroutine cga_da
 
-        integer :: i, j, k, l, ibasis, nei
+        integer :: i, j, k, l, ibasis, nei, close_neighbors(2,4,6), far_neighbor(4,6), minindex, minindices(2), linevec(3)
 
-        integer :: face_types(ele_num*6), face_ele(6*ele_num)
+        integer :: face_types(ele_num*6), face_ele(6*ele_num), diff_pair(2,6)
-        real(kind = dp) :: face_centroids(3, ele_num*6), r(3), rc, vnode(3, max_basisnum, 4), vnorm(max_basisnum, 4), &
+        real(kind = dp) :: face_centroids(3, ele_num*6), r(3), rc, vnode(3, max_basisnum, 4), vnorm, &
-                           max_node_dist, ndiff, rmax(3), rmin(3), rnorm
+                           max_node_dist, ndiff, rmax(3), rmin(3), rnorm, v1(3), v2(3), v1norm, v2norm, theta1, theta2, &
+                            diff_mag(6), diffvec(3)
+        real(kind=dp), allocatable :: disl_line_temp(:,:), disl_data_temp(:)
+        logical :: is_edge
 
         !Initialize variables
         l = 0
         max_node_dist = 0
 
+        !Now save the close and far neighbors of the nodes. This is done to attempt to figure out the character of the dislocation
+        !by comparing how the burgers vector is distributed over the face.
+        do j = 1, 6
+            close_neighbors(1, 1, j) = cubic_faces(4, j)
+            close_neighbors(2, 1, j) = cubic_faces(2, j)
+            far_neighbor(1,j) = cubic_faces(3,j)
+
+            close_neighbors(1, 2, j) = cubic_faces(1, j)
+            close_neighbors(2, 2, j) = cubic_faces(3, j)
+            far_neighbor(2,j) = cubic_faces(4,j)
+
+            close_neighbors(1, 3, j) = cubic_faces(2, j)
+            close_neighbors(2, 3, j) = cubic_faces(4, j)
+            far_neighbor(3,j) = cubic_faces(1,j)
+
+            close_neighbors(1, 4, j) = cubic_faces(3, j)
+            close_neighbors(2, 4, j) = cubic_faces(1, j)
+            far_neighbor(4,j) = cubic_faces(2,j)
+        end do
+
         !First calculate all of the face centroids
         do i = 1, ele_num
             do j = 1, 6
@@ -101,8 +138,6 @@ module mode_da
 
 
             end do
-            !Calculate the largest difference between nodes for each face 
-
         end do
 
         !Now set the cut off distance as half the distance from opposite nodes in the largest element
@@ -131,31 +166,121 @@ module mode_da
                     end do
                 end do
 
-                do j = 1, 3
+                !Now calculate the maximum difference between nodes
-                    vnode(j,1:basisnum(lat_ele(face_ele(i))),:) = vnode(j,1:basisnum(lat_ele(face_ele(i))),:) - &
+                vnorm = 0
-                                                                  minval(vnode(j,1:basisnum(lat_ele(face_ele(i))),:))
-                end do
                 do j = 1, 4
-                    do ibasis = 1, basisnum(lat_ele(face_ele(i)))
+                    do k=j,4
-                        vnorm(ibasis,j) = norm2(vnode(:, ibasis, j))
+                        v1=vnode(:,1,j) - vnode(:,1,k)
+                        v1norm=norm2(v1)
+                        if (vnorm < v1norm) then
+                            vnorm = v1norm
+                            diffvec= v1
+                        end if 
                     end do
                 end do
                 !Now calculate the difference between the largest norm and the smallest norm, if it's larger than 0.5 then mark it
                 !as slipped. This value probably can be converted to a variable value that depends on the current lattice parameter
                 !I think 0.5 works ok though.
-                if (any(vnorm > 0.5_dp)) then 
+                if (vnorm>0.5_dp) then 
                     print *, "Element number ", face_ele(i), " is dislocated along face ", face_types(i), &
-                             " with neighbor ", face_ele(nei), " with max displacement of ", maxval(vnorm) 
+                             " with neighbor ", face_ele(nei), " with max displacement of ", vnorm
+                    l=0
                     do j = 1, 4
                         is_slipped(:, cubic_faces(j,face_types(i)), face_ele(i)) = 1
-                        print *, j, vnode(:,1,j)
+
+                        !This portion of the code is used to determine what the character of the dislocation most likely is
+                        !effectively how this works is we look for the two nodes with the most similar burgers vector.
+                        !The vector between these two nodes is close to the line direction and as a result we can probably estimate
+                        !if it's close to edge, screw, or if it's pretty fairly mixed.
+                        do k = j+1, 4
+                            l=l+1
+                            diff_mag(l) = norm2(vnode(:, 1, j)-vnode(:,1,k))
+                            diff_pair(1,l)=cubic_faces(j, face_types(i))
+                            diff_pair(2,l)=cubic_faces(k, face_types(i))
+                        end do
+                    end do
+                    minindex = minloc(diff_mag,1)
+                    !Now figure out if the min difference between nodes is associated with a 112 direction    
+                    is_edge = .false.
+                    do j = 1,6
+                        if((diff_pair(1,minindex) == oneonetwopairs(1,j)).and.(diff_pair(2,minindex)==oneonetwopairs(2,j))) then
+                            is_edge=.true.
+                            exit
+                        else if((diff_pair(2,minindex)==oneonetwopairs(1,j)).and.(diff_pair(1,minindex)==oneonetwopairs(2,j)))then
+                            is_edge=.true.
+                            exit
+                        end if
                     end do
+
+                    if(is_edge) then 
+                        print *, 'Dislocation has primarily edge character'
+                    else
+                        print *, 'Dislocation has primarily screw character'
+                    end if
+
+                    if( i < nei) then 
+                        num_d_points=num_d_points+2
+                        if(num_d_points > size(disl_line,2)) then 
+                            allocate(disl_line_temp(3,size(disl_line,2)+1024))
+                            disl_line_temp=0.0_dp
+                            disl_line_temp(:, 1:size(disl_line,2))=disl_line
+                            call move_alloc(disl_line_temp, disl_line)
+                        end if
+                        if(num_d_points/2 > size(disl_data)) then 
+                            allocate(disl_data_temp(size(disl_data)+1024))
+                            disl_data_temp=0
+                            disl_data_temp(1:size(disl_data)) = disl_data
+                            call move_alloc(disl_data_temp, disl_data)
+                        end if
+                        print *, num_d_points/2, vnorm
+                        disl_data(num_d_points/2)=vnorm
+                        disl_line(:,num_d_points-1) =r_node(:, 1, diff_pair(1,minindex), face_ele(i))
+                        disl_line(:,num_d_points) =r_node(:, 1, diff_pair(2,minindex), face_ele(i))
+                    end if
                 end if
             end if
         end do
 
     end subroutine cga_da
 
+    subroutine write_da_vtk(outfile)
+        character (len=*), intent(in) :: outfile
+
+        integer :: i 
+        open(unit=11, file=outfile, action='write', status='replace', position='rewind')
+
+        write(11, '(a)') '# vtk DataFile Version 4.0.1'
+        write(11, '(a)') 'Dislocation line analyzed via cacmb dislocation analysis'
+        write(11, '(a)') 'ASCII'
+        write(11, '(a)') 'DATASET UNSTRUCTURED_GRID'
+        write(11, *) 'POINTS', num_d_points, 'float'
+       
+        !Write the points to the file
+        do i = 1, num_d_points
+            write(11, '(3f23.15)') disl_line(:, i)
+        end do
+
+        write(11,*) 'CELLS', num_d_points/2, num_d_points/2*3
+        do i=1, num_d_points, 2
+            write(11,'(3i16)') 2, i-1, i
+        end do
+        
+        write(11,*) 'CELL_TYPES', num_d_points/2
+        do i=1, num_d_points, 2
+            write(11,'(i16)') 3
+        end do
+
+        write(11,*) 'CELL_DATA', num_d_points/2
+        write(11, '(a)') 'SCALARS burgers_vec, float'
+        write(11, '(a)') 'LOOKUP_TABLE default'
+        do i=1, num_d_points, 2
+            write(11, '(f23.15)') disl_data((i+1)/2)
+        end do
+
+        close(11)
+
+    end subroutine write_da_vtk
+
     subroutine write_da_xyz(outfile)
         !This subroutine write the element positions to a .xyz file and marks whether they are slipped or not 
         character(len=*), intent(in) :: outfile
@@ -184,4 +309,4 @@ module mode_da
         
         close(11)
     end subroutine write_da_xyz
-end module
+end module mode_da

src/mode_metric.f90

@@ -10,7 +10,7 @@ module mode_metric
     
     integer :: nfiles
     character(len=100) :: metric_type
-    real(kind=dp) :: rc_off
+    real(kind=dp) :: rc_off, old_box_len(3)
 
     !Save reference positions
     integer :: np, npreal, nmet
@@ -57,6 +57,10 @@ module mode_metric
         print *, "Getting initial neighbor list"
         call calc_neighbor(rc_off, r_zero, np)
 
+        !Save box lengths for applying periodic boundaries 
+        do i=1,3
+            old_box_len(i) = box_bd(2*i) - box_bd(2*i-1)
+        end do
         !Reset element and box
         call reset_data
         call reset_box 
@@ -128,8 +132,11 @@ module mode_metric
         !This subroutine calculates the continuum metric that we require
         integer :: i, j, k, nei, ip, jp
         real(kind=dp) :: def_grad(3,3), omega(3,3), eta(3,3), rij(3), eta_inv(3,3), ftf(3,3), &
-                         U(3,3), R(3,3), Rskew(3,3), oldrij(3)
+                         U(3,3), R(3,3), Rskew(3,3), oldrij(3), box_len(3)
 
+        do i = 1, 3
+            box_len(i) = box_bd(2*i) - box_bd(2*i-1)
+        end do
         !Loop over all points
         do ip = 1, np
             eta(:,:) = 0.0_dp
@@ -138,8 +145,17 @@ module mode_metric
             do jp = 1, nei_num(ip)
                 !Calculate the neighbor vec in current configuration
                 nei = nei_list(jp, ip)
-                rij = r_curr(:,nei) - r_curr(:,ip)
+                rij = r_curr(:,nei) - r_curr(:,ip)  +periodvec(:,jp,ip)*box_len 
-                oldrij = r_zero(:,nei) - r_zero(:,ip)
+                if (norm2(rij) > 5*rc_off) then 
+                    do j = 1, 3
+                        if (r_curr(j,nei) > r_curr(j,ip)) then 
+                            rij(j) = rij(j) - box_len(j)
+                        else if(r_curr(j,nei) < r_curr(j,ip)) then 
+                            rij(j) = rij(j) + box_len(j)
+                        end if
+                    end do
+                end if
+                oldrij = r_zero(:,nei) - r_zero(:,ip) + periodvec(:,jp,ip)*old_box_len 
 
                 !Calculate  eta and omega
                 do i = 1,3

src/neighbors.f90

@@ -4,10 +4,12 @@ module neighbors
     use elements
     use subroutines
     use functions
+    use box
 
-    integer, allocatable :: nei_list(:,:), nei_num(:), nn(:)
+    integer, allocatable :: nei_list(:,:), nei_num(:), nn(:), periodvec(:,:,:)
     real(kind=dp), allocatable :: init_vec(:,:,:), output(:), microrotation(:,:)
     public 
+
     contains
  
     subroutine build_cell_list(numinlist, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell)
@@ -80,44 +82,70 @@ module neighbors
         
         real(kind=dp), intent(in) :: rc_off
         integer, intent(in) :: n
-        real(kind=dp), dimension(3,n) :: r_list
+        real(kind=dp), dimension(:, :), intent(in) :: r_list
 
-        integer :: i, c(3), ci, cj, ck, num_nei, nei
+        integer :: i, j, c(3),cn(3), ci, cj, ck, num_nei, nei, v(3), period_dir(3)
         !Variables for cell list code
         integer, dimension(3) ::cell_num
         integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:)
         integer :: which_cell(3,n)
+        real(kind=dp) :: r(3), box_len(3)
+        logical :: period_bd(3)
         
-        !First reallocate the neighbor list codes
+        !First reallocate the neighbor list variables
         if (allocated(nei_list)) then
-            deallocate(nei_list,nei_num)
+            deallocate(nei_list,nei_num, periodvec)
         end if
 
-        allocate(nei_list(100,n),nei_num(n))
+        allocate(nei_list(100,n),nei_num(n), periodvec(3,100,n))
+        nei_list=0
+        periodvec=0
+        nei_num=0
 
         !Now first pass the position list and and point num to the cell algorithm
         call build_cell_list(n, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell)
 
+        do i=1, 3
+            if (box_bc(i:i) == 'p') then 
+                period_bd(i) = .true.
+            else
+                period_bd(i)=.false.
+            end if
+            box_len(i) = box_bd(2*i) - box_bd(2*i-1)
+        end do
         !Now loop over every point and find it's neighbors
         pointloop: do i = 1, n
             
             !First check to see if the point is a filler point, if so then skip it 
             if(r_list(1,i) < box_bd(1)) cycle
 
-            !c is the position of the cell that the point
+            !c is the position of the cell that the point belongs to 
             c = which_cell(:,i)
 
             !Loop over all neighboring cells
             do ci = -1, 1, 1
                 do cj = -1, 1, 1
                     do ck = -1, 1, 1
-                        !First check to make sure that the neighboring cell exists
+                        !First try to apply periodic boundaries 
+                        v=(/ ck, cj, ci /)
+                        cn=0
+                        period_dir=0
+                        do j=1, 3
+                            if ((c(j) + v(j) == 0).and.period_bd(j)) then 
+                                cn(j) = cell_num(j)
+                                period_dir(j) = 1
+                            else if ((c(j) + v(j) > cell_num(j)).and.period_bd(j)) then 
+                                cn(j) = 1
+                                period_dir(j) = -1
+
+                            else if ((c(j)+v(j) >= 1) .and. (c(j)+v(j) <= cell_num(j))) then
+                                cn(j) = c(j) + v(j)
+                            end if
+                        end do
                         if(any((c + (/ ck, cj, ci /)) == 0)) cycle
-                        if( (c(1) + ck > cell_num(1)).or.(c(2) + cj > cell_num(2)).or. &
-                            (c(3) + ci > cell_num(3))) cycle
 
-                        do num_nei = 1, num_in_cell(c(1) + ck, c(2) + cj, c(3) + ci)
+                        do num_nei = 1, num_in_cell(cn(1),cn(2),cn(3))
-                            nei = cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci)
+                            nei = cell_list(num_nei,cn(1),cn(2),cn(3))
 
                             !Check to make sure the atom isn't the same index as the atom we are checking 
                             !and that the neighbor hasn't already been deleted
@@ -125,10 +153,13 @@ module neighbors
 
                                 !Now check to see if it is in the cutoff radius, if it is add it to the neighbor list for that
                                 !atom and calculate the initial neighbor vector
-                                if (norm2(r_list(:,nei)-r_list(:,i)) < rc_off) then 
+                        
+                                r = r_list(:,nei) + period_dir*box_len
+                                if (norm2(r-r_list(:,i)) < rc_off) then 
                                     
                                     nei_num(i) = nei_num(i) + 1
                                     nei_list(nei_num(i), i) = nei
+                                    periodvec(:,nei_num(i),i) = period_dir
 
                                 end if
                             end if

src/opt_bubble.f90

@@ -1,6 +1,7 @@
 module opt_bubble
     !This module contains the bubble option which can be used to create voids with specific pressures of certain atoms
 
+    use atoms
     use parameters
     use elements 
     use box
@@ -8,7 +9,7 @@ module opt_bubble
 
     implicit none
 
-    real(kind=dp), private :: br, bt, bp, c(3)
+    real(kind=dp), private :: br, brat, c(3)
     character(len=2), private :: species
 
     public 
@@ -16,8 +17,8 @@ module opt_bubble
     subroutine bubble(arg_pos)
         integer, intent(inout) :: arg_pos
 
-        integer :: new_type, n, j, i
+        integer :: new_type, n, j, i, atom_num_pre
-        real(kind = dp)  :: p(3), rand, factor, per, vol
+        real(kind = dp)  :: p(3), rand, factor, per, vol, mass
 
         print *, '------------------------------------------------------------'
         print *, 'Option Bubble'
@@ -28,34 +29,36 @@ module opt_bubble
         !Now we use the existing group code to delete a sphere which represents the bubble
         centroid=c
         radius = br
-        type='both'
+        type='all'
         gshape='sphere'
         group_nodes = .true.
-
+        group_atom_types=0
         call get_group
         call refine_group
         call get_group
+        atom_num_pre= atom_num
         call delete_group
 
         !Now we create a new atom type with the desired species
-        call add_atom_type(species, new_type)
+        call atommass(species, mass)
+        call add_atom_type(mass, new_type)
 
         !Now we calculate the number of atoms we need to add for the desired pressure
-        print *, "Creating a bubble with pressure", bp, " at temperature ", bt, " with radius ", br
+        !print *, "Creating a bubble with pressure", bp, " at temperature ", bt, " with radius ", br
-        
+        !
-        factor = 1.0e24/6.02214e23
+        !factor = 1.0e24/6.02214e23
-        if (bp <= 10.0) then
+        !if (bp <= 10.0) then
-            per=factor*(3.29674113+4.51777872*bp**(-0.80473167))
+        !    per=factor*(3.29674113+4.51777872*bp**(-0.80473167))
-        else if (bp .le. 20.0) then
+        !else if (bp .le. 20.0) then
-            per=factor*(4.73419689-0.072919483*bp)
+        !    per=factor*(4.73419689-0.072919483*bp)
-        else
+        !else
-            per=factor*(4.73419689-0.072919483*bp)
+        !    per=factor*(4.73419689-0.072919483*bp)
-            print *, 'warning: pressure is too high'
+        !    print *, 'warning: pressure is too high'
-            print *, 'equation of state is only valid for < 20 GPa'
+        !    print *, 'equation of state is only valid for < 20 GPa'
-        endif 
+        !endif 
-        vol = 4.0*pi/3.0*br**3.0
+        !vol = 4.0*pi/3.0*br**3.0
-
+
-        n = vol/per
+        n = brat*(atom_num_pre-atom_num)
         print *, "Adding ", n, " atoms of species ", species
 
         !Now add n atoms randomly within the sphere 
@@ -77,39 +80,66 @@ module opt_bubble
         integer, intent(inout) :: arg_pos
 
         integer :: i, arglen
+        real(kind=dp) :: mass
         character(len=100) :: tmptxt
 
+        
         !First read in the bubble centroid
         do i = 1, 3
             arg_pos = arg_pos + 1 
             call get_command_argument(arg_pos, tmptxt, arglen)
+            print *,  trim(adjustl(tmptxt))
             call parse_pos(i, tmptxt, c(i))
         end do
 
         !Now the bubble radius
         arg_pos = arg_pos + 1
         call get_command_argument(arg_pos, tmptxt, arglen)
+        print *,  trim(adjustl(tmptxt))
         if(arglen == 0) stop "Missing bubble radius"
         read(tmptxt, *) br
 
-        !Now bubble pressure
+        !Now bubble ratio
-        arg_pos = arg_pos + 1
-        call get_command_argument(arg_pos, tmptxt, arglen)
-        if(arglen == 0) stop "Missing bubble pressure"
-        read(tmptxt, *) bp
-        
-        !Now bubble temperature
         arg_pos = arg_pos + 1
         call get_command_argument(arg_pos, tmptxt, arglen)
-        if(arglen == 0) stop "Missing bubble temperature"
+        print *,  trim(adjustl(tmptxt))
-        read(tmptxt, *) bt
+        if(arglen == 0) stop "Missing bubble ratio"
+        read(tmptxt, *) brat
         
         !Now the bubble species
         arg_pos = arg_pos + 1
         call get_command_argument(arg_pos, species, arglen)
+        print *,  species
         if(arglen == 0) stop "Missing bubble species"
 
+
+        !OPtional arguments
+        do while(.true.)
+            if(arg_pos > command_argument_count()) exit
             arg_pos=arg_pos+1
+
+            call get_command_argument(arg_pos, tmptxt)
+            tmptxt=adjustl(tmptxt)
+            select case(trim(tmptxt))
+
+            case('excludetypes')
+                arg_pos=arg_pos + 1
+                call get_command_argument(arg_pos, tmptxt, arglen)
+                if(arglen==0) stop "Missing number of atoms to exclude"
+                read(tmptxt, *) exclude_num
+                do i=1,exclude_num
+                    arg_pos = arg_pos + 1
+                    call get_command_argument(arg_pos, tmptxt, arglen)
+                    if(arglen==0) stop "Missing exclude atom types"
+                    call atommass(tmptxt, mass)
+                    call add_atom_type(mass, exclude_types(i))
+                end do
+
+            case default
+                exit
+            end select
+        end do
+        
         return
     end subroutine parse_bubble
 end module opt_bubble

src/opt_delete.f90

@@ -87,7 +87,7 @@ module opt_delete
         allocate(which_cell(3,atom_num))
 
         !First pass the atom list and atom num to the algorithm which builds the cell list
-        call build_cell_list(atom_num, r_atom, rc_off, cell_num, num_in_cell, cell_list, which_cell)
+        call build_cell_list(atom_num, r_atom, 5*rc_off, cell_num, num_in_cell, cell_list, which_cell)
 
         !Now loop over every atom and figure out if it has neighbors within the rc_off
         delete_num = 0

src/opt_disl.f90

@@ -520,7 +520,7 @@ module opt_disl
         end do
         return
 
-    end subroutine
+    end subroutine disloop
 
     !********************************************************
     ! SOLIDANGLE

src/opt_group.f90

@@ -2,6 +2,7 @@ module opt_group
 
     !This module contains all code associated with dislocations
 
+    use atoms
     use parameters
     use elements
     use subroutines
@@ -11,7 +12,7 @@ module opt_group
     implicit none
 
     integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize, insert_type, &
-               insert_site, num_species
+               insert_site, num_species, group_atom_types, exclude_num, exclude_types(10)
     character(len=15) :: type, gshape!Type indicates what element type is selected and shape is the group shape
     character(len=2) :: species_type(10)
     real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness, insert_conc, &
@@ -38,6 +39,7 @@ module opt_group
         insert_type = 0
         random_num=0
         group_type=0
+        group_atom_types=0
         notsize=0
         displace=.false.
         delete=.false.
@@ -46,6 +48,7 @@ module opt_group
         flip = .false.
         group_nodes=.false.
         num_species = 0
+        exclude_num = 0
 
         if(allocated(element_index)) deallocate(element_index)
         if(allocated(atom_index)) deallocate(atom_index)
@@ -86,10 +89,10 @@ module opt_group
             call displace_group
         end if
         if(insert_type > 0) then 
-            print *, "Insert command has been dropped"
+            !print *, "Insert command has been dropped"
-            stop 1
+            !stop 1
             call get_group
-!            call insert_group
+            call insert_group
         end if
 
         if(num_species > 0) then 
@@ -106,7 +109,7 @@ module opt_group
 
         integer :: i, j, arglen, in_num
         character(len=100) :: textholder, type_spec
-        real(kind=dp)  H
+        real(kind=dp)  H, mass
 
         !Parse type and shape command
         arg_pos = arg_pos + 1
@@ -469,6 +472,24 @@ module opt_group
                 refine=.true.
             case('refinefill')
                 refinefill=.true.
+            case('selecttype')
+                arg_pos = arg_pos + 1
+                call get_command_argument(arg_pos, textholder, arglen)
+                if (arglen==0) stop "Missing select type"
+                call atommass(textholder, mass)
+                call add_atom_type(mass, group_atom_types)
+            case('excludetypes')
+                arg_pos=arg_pos + 1
+                call get_command_argument(arg_pos, textholder, arglen)
+                if(arglen==0) stop "Missing number of atoms to exclude"
+                read(textholder, *) exclude_num
+                do i=1,exclude_num
+                    arg_pos = arg_pos + 1
+                    call get_command_argument(arg_pos, textholder, arglen)
+                    if(arglen==0) stop "Missing exclude atom types"
+                    call atommass(textholder, mass)
+                    call add_atom_type(mass, exclude_types(i))
+                end do
             case('remesh')
                 arg_pos = arg_pos + 1
                 call get_command_argument(arg_pos, textholder, arglen)
@@ -493,7 +514,8 @@ module opt_group
                 arg_pos = arg_pos + 1
                 call get_command_argument(arg_pos, textholder, arglen)
                 if (arglen==0) stop "Missing atom type for group"
-                call add_atom_type(textholder, group_type)
+                call atommass(textholder, mass)
+                call add_atom_type(mass, group_type)
             case('notsize')
                 arg_pos = arg_pos + 1
                 call get_command_argument(arg_pos, textholder, arglen)
@@ -521,7 +543,8 @@ module opt_group
                 arg_pos=arg_pos+1
                 call get_command_argument(arg_pos, textholder, arglen)
                 if (arglen==0) stop "Missing element type for insert command"
-                call add_atom_type(textholder, insert_type)
+                call atommass(textholder, mass)
+                call add_atom_type(mass, insert_type)
                 arg_pos=arg_pos + 1 
                 call get_command_argument(arg_pos, textholder, arglen)
                 select case(trim(adjustl(textholder)))
@@ -616,7 +639,7 @@ module opt_group
         select case(trim(adjustl(type)))
         case('atoms','all')
             do i = 1, atom_num
-                if(in_group(r_atom(:,i)).neqv.flip) then 
+                if(in_group(type_atom(i), r_atom(:,i)).neqv.flip) then 
                     group_atom_num = group_atom_num + 1
                     if (group_atom_num > size(atom_index)) then 
                         allocate(resize_array(size(atom_index) + 1024))
@@ -643,15 +666,6 @@ module opt_group
             end if
         end select
         
-        j = 0
-        do i = 1, group_atom_num
-            if (atom_index(i) == 23318348) then
-                j = j + 1
-            end if
-        end do  
-
-        if (j > 1) stop "Code broken"
-
         print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms."
     end subroutine get_group 
 
@@ -693,8 +707,19 @@ module opt_group
         !This command is used to refine the group to full atomistics
 
         integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num
-        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3)
+        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), data_interp(10, max_basisnum*max_esize**3),  &
+                         vel_interp(3, max_basisnum*max_esize**3) 
+        real(kind=dp), allocatable :: vel_tmp(:,:)
+        logical :: refine_dat
 
+
+        if (allocated(force_node)) then
+            refine_dat=.true.
+        else 
+            refine_dat=.false.
+        end if
+
+        print *, size(data_interp)
         !Refining to atoms
         if(group_ele_num > 0) then 
             orig_atom_num = atom_num
@@ -704,17 +729,38 @@ module opt_group
             do i = 1, group_ele_num
                 ie = element_index(i)
                 !Get the interpolated atom positions
+                if (refine_dat) then 
+                    call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp, &
+                                           energy_node(:,:,i), force_node(:,:,:,i), virial_node(:,:,:,i), data_interp) 
+                else
                     call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
+                end if
+                if(allocated(vel_atom)) then 
+                    call interpolate_vel(type_ele(i), size_ele(i), lat_ele(i), vel_node(:,:,:,i), vel_interp)
+                end if
                 !Loop over all interpolated atoms and add them to the system, we apply periodic boundaries 
                 !here as well to make sure they are in the box
                 do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3
                     call apply_periodic(r_interp(:,j))
                     call add_atom(0,type_interp(j), r_interp(:,j))
+                    if(refine_dat) then 
+                        call add_atom_data(atom_num, data_interp(1, j), data_interp(2:4, j), data_interp(5:10, j))
+                    end if
+                    if(allocated(vel_atom)) then
+                        if(size(vel_atom, 2) < size(type_atom)) then 
+                            allocate(vel_tmp(3, size(type_atom)))
+                            vel_tmp=0.0_dp
+                            vel_tmp(:, 1:size(vel_atom,2)) = vel_atom
+                            call move_alloc(vel_tmp, vel_atom)
+                        end if
+                        vel_atom(:, atom_num) = vel_interp(:,j)
+                    end if
                 end do 
             end do
             !Once all atoms are added we delete all of the elements
             call delete_elements(group_ele_num, element_index)
             print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms."
+
         end if                
 
     end subroutine refine_group
@@ -751,7 +797,7 @@ module opt_group
                                 ravg = ravg + rfill(:,ibasis, 1)/basisnum(lat_ele(ie))
                             end do 
 
-                            if( in_group(ravg)) then 
+                            if( in_group(0,ravg)) then 
                                 nump_ele = nump_ele - 1
                             end if
                         end do
@@ -1124,86 +1170,102 @@ module opt_group
 
     end subroutine change_group_type
     
-!    subroutine insert_group
+    subroutine insert_group
-!        !This code inserts atoms into interstitial sites. This only works on atoms within the group due to the limitations with the
+        !This code inserts atoms into interstitial sites. This only works on atoms within the group due to the limitations with the
-!        !Coarse-graining methodology which doesn't allow for concentration fluctuations.
+        !Coarse-graining methodology which doesn't allow for concentration fluctuations.
-!        real(kind=dp)  interstitial_sites(3,14), rand, rinsert(3)
+        real(kind=dp)  interstitial_sites(3,14), rand, rinsert(3)
-!        integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi
+        integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi
-!        integer, allocatable :: used_sites(:,:)
+        integer, allocatable :: used_sites(:,:)
-!
+
-!        !First save all of the displacement vectors from a lattice site to interstitial site
+        !First save all of the displacement vectors from a lattice site to interstitial site
-!        !The first 6 are the octohedral sites and the next 8 are the tetrahedral sites
+        !The first 6 are the octohedral sites and the next 8 are the tetrahedral sites
-!        interstitial_sites= reshape( (/ -0.5_dp, 0.0_dp, 0.0_dp, &
+        interstitial_sites= reshape( (/ -0.5_dp, 0.0_dp, 0.0_dp, &
-!                                         0.5_dp, 0.0_dp, 0.0_dp, &
+                                         0.5_dp, 0.0_dp, 0.0_dp, &
-!                                         0.0_dp,-0.5_dp, 0.0_dp, &
+                                         0.0_dp,-0.5_dp, 0.0_dp, &
-!                                         0.0_dp, 0.5_dp, 0.0_dp, &
+                                         0.0_dp, 0.5_dp, 0.0_dp, &
-!                                         0.0_dp, 0.0_dp,-0.5_dp, &
+                                         0.0_dp, 0.0_dp,-0.5_dp, &
-!                                         0.0_dp, 0.0_dp, 0.5_dp, &
+                                         0.0_dp, 0.0_dp, 0.5_dp, &
-!                                        -0.25_dp,-0.25_dp,-0.25_dp, &
+                                        -0.25_dp,-0.25_dp,-0.25_dp, &
-!                                        -0.25_dp,-0.25_dp, 0.25_dp, &
+                                        -0.25_dp,-0.25_dp, 0.25_dp, &
-!                                        -0.25_dp, 0.25_dp,-0.25_dp, &
+                                        -0.25_dp, 0.25_dp,-0.25_dp, &
-!                                        -0.25_dp, 0.25_dp, 0.25_dp, &
+                                        -0.25_dp, 0.25_dp, 0.25_dp, &
-!                                         0.25_dp,-0.25_dp,-0.25_dp, &
+                                         0.25_dp,-0.25_dp,-0.25_dp, &
-!                                         0.25_dp,-0.25_dp, 0.25_dp, &
+                                         0.25_dp,-0.25_dp, 0.25_dp, &
-!                                         0.25_dp, 0.25_dp,-0.25_dp, &
+                                         0.25_dp, 0.25_dp,-0.25_dp, &
-!                                         0.25_dp, 0.25_dp, 0.25_dp /), &
+                                         0.25_dp, 0.25_dp, 0.25_dp /), &
-!                                         shape(interstitial_sites))
+                                         shape(interstitial_sites))
-!
+
-!        !First we calculate the number of atoms needed based on the number of atoms in the group and the concentration
+        !First we calculate the number of atoms needed based on the number of atoms in the group and the concentration
-!        interstitial_sites=interstitial_sites*insert_lattice
+        interstitial_sites=interstitial_sites*insert_lattice
-!
+
-!        add_num = (insert_conc*group_atom_num)/(1-insert_conc)
+
-!        allocate(used_sites(2,add_num))
+        !Now set up the random number generator for the desired interstitial type
-!
+        select case(insert_site)
-!        print *, "Inserting ", add_num, " atoms as atom type ", insert_type
+        case(1)
-!
+            rlo=1
-!        !Now set up the random number generator for the desired interstitial type
+            rhi=6
-!        select case(insert_site)
+        case(2)
-!        case(1)
+            rlo=7 
-!            rlo=1
+            rhi = 14
-!            rhi=6
+        case(3)
-!        case(2)
+            rlo=1
-!            rlo=7 
+            rhi=14
-!            rhi = 14
+        end select
-!        case(3)
+
-!            rlo=1
+        if ((insert_conc > 1).or.(insert_conc < 0)) stop "Insert concentration should be between 0 and 1"
-!            rhi=14
+        add_num = (insert_conc*group_atom_num)
-!        end select
+        allocate(used_sites(2,add_num))
-!
+
-!        !Now add the atoms
+        print *, "Inserting ", add_num, " atoms as atom type ", insert_type
-!        i = 1
+
-!        addloop:do while ( i < add_num)
+        !Now add the atoms
-!            call random_number(rand)
+        i = 1
-!            rindex = int(1+rand*(group_atom_num-1))
+        addloop:do while ( i < add_num)
-!            ia=atom_index(rindex)
+            call random_number(rand)
-!            call random_number(rand)
+            rindex = int(1+rand*(group_atom_num-1))
-!            sindex = int(rlo+rand*(rhi-rlo))
+            ia=atom_index(rindex)
-!            do j = 1, i
+            call random_number(rand)
-!                if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop
+            sindex = int(rlo+rand*(rhi-rlo))
-!            end do
+            do j = 1, i
-!            rinsert = r_atom(:,ia) + matmul(sub_box_ori(:,:,sbox_atom(ia)),interstitial_sites(:,sindex))
+                if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop
-!            sbox = sbox_atom(ia)
+            end do
-!            call add_atom(0, insert_type, sbox, rinsert)
+            rinsert = r_atom(:,ia) + matmul(box_ori,interstitial_sites(:,sindex))
-!            used_sites(1,i) = ia
+            call add_atom(0, insert_type, rinsert)
-!            used_sites(2,i) = sindex
+            used_sites(1,i) = ia
-!            i = i + 1
+            used_sites(2,i) = sindex
-!        end do addloop
+            i = i + 1
-!    end subroutine insert_group
+        end do addloop
+    end subroutine insert_group
 
     subroutine alloy_group
         !This subroutine randomizes the atom types to reach desired concentrations, this only operates on atoms
         integer :: i, j, ia, type_map(10), added_types(num_species)
-        real(kind=dp) ::  rand
+        real(kind=dp) ::  rand, mass, old_fractions(num_species)
+        character(len=2) :: sorted_species_type(10)
+        logical :: species_mapped(10)
 
 
         print *, "Alloying group with desired fractions", s_fractions(1:num_species)
+        old_fractions=s_fractions(1:num_species)
         !First we sort the fractions
         call quicksort(s_fractions(1:num_species))
         
+        species_mapped = .false.
+        do i = 1,num_species
+            do j=1, num_species
+                if (is_equal(s_fractions(i), old_fractions(j)) .and. .not. species_mapped(j)) then 
+                    sorted_species_type(i)=species_type(j)
+                    species_mapped(j) = .true.
+                    exit
+                end if
+            end do
+        end do
+
 
         !Now get the atom type maps for all the atoms and make the fractions a running sum
         do i = 1, num_species
-            call add_atom_type(species_type(i), type_map(i))
+            call atommass(sorted_species_type(i), mass)
+            call add_atom_type(mass, type_map(i))
             if(i > 1) s_fractions(i) = s_fractions(i) + s_fractions(i-1)
         end do
         !Now randomize the atom types
@@ -1225,8 +1287,9 @@ module opt_group
         return
     end subroutine alloy_group
 
-    function in_group(r)
+    function in_group(atype, r)
         !This subroutine determines if a point is within the group boundaries
+        integer, intent(in) :: atype
         real(kind=dp), intent(in) :: r(3)
         real(kind=dp) :: rnorm
 
@@ -1271,6 +1334,19 @@ module opt_group
         case('all')
             in_group = .true.
         end select 
+
+        if(group_atom_types> 0) then 
+            in_group = in_group.and.(atype== group_atom_types)
+        elseif(exclude_num > 0) then 
+            if(in_group) then
+                do i=1,exclude_num 
+                    if (atype == exclude_types(i)) then 
+                        in_group=.false.
+                        exit
+                    end if
+                end do
+            end if
+        end if
     end function in_group
 
     function in_group_ele(esize, elat, rn)
@@ -1294,13 +1370,13 @@ module opt_group
                     r_center(:)= r_center(:) + rn(:,ibasis,inod)/basisnum(elat) 
                 end do
 
-                if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(size_ele(i)/=notsize)) then 
+                if((in_group(0,rn(:, ibasis, inod)).neqv.flip).and.(size_ele(i)/=notsize)) then 
                     node_in_out(inod) = 2
                     exit nodeloop 
                 end if
                 
                 gshape ='sphere'
-                if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then 
+                if((in_group(0,rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then 
                     node_in_out(inod) = 1
                 else
                     node_in_out(inod) = 0
@@ -1322,7 +1398,7 @@ module opt_group
                 end do
             end do
 
-            if ((in_group(r_center).neqv.flip).and.(esize/= notsize)) then 
+            if ((in_group(0,r_center).neqv.flip).and.(esize/= notsize)) then 
                 in_group_ele=.true.
                 return
             end if 
@@ -1331,7 +1407,7 @@ module opt_group
             r_center(:) = 0.0_dp
             do inod = 1, ng_node(elat)
                 do ibasis = 1, basisnum(elat)
-                     if ((in_group(rn(:,ibasis,inod)).neqv.flip).and.(esize/=notsize)) then 
+                     if ((in_group(0,rn(:,ibasis,inod)).neqv.flip).and.(esize/=notsize)) then 
                          in_group_ele=.true.
                          return
                     end if 

src/opt_slip_plane.f90

@@ -17,11 +17,13 @@ module opt_slip_plane
         !Main calling function for the slip_plane option
         integer, intent(inout) :: arg_pos 
 
-        integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), ibasis
+        integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), ibasis,&
+                    starting_anum
                 
         integer, allocatable :: slip_eles(:), temp_int(:)
         real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node), ratom(3,max_basisnum), &
-                         maxp, minp
+                         maxp, minp, vel_interp(3, max_basisnum*max_esize**3)
+        real(kind=dp), allocatable :: vel_tmp(:,:)
 
         integer :: type_interp(max_basisnum*max_esize**3)
         logical :: lat_points(max_esize,max_esize, max_esize)
@@ -58,11 +60,25 @@ module opt_slip_plane
                 
                 !If we aren't efilling then just refine the element
                 if(.not.efill) then 
+                    starting_anum=atom_num
                     call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
                     do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3           
                         call apply_periodic(r_interp(:,ia))
                         call add_atom(0, type_interp(ia), r_interp(:,ia))
                     end do
+
+                    if(allocated(vel_atom)) then  
+                        call interpolate_vel(type_ele(ie), size_ele(ie), lat_ele(ie), vel_node(:,:,:,ie), vel_interp)
+                        if(size(vel_atom,2) < size(type_atom)) then 
+                            allocate(vel_tmp(3, size(type_atom)))
+                            vel_tmp=0.0_dp
+                            vel_tmp(:, 1:size(vel_atom,2)) = vel_atom
+                            call move_alloc(vel_tmp, vel_atom)
+                        end if
+                        do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3           
+                            vel_atom(:, starting_anum+ia) = vel_interp(:, ia)
+                        end do
+                    end if
                 !If we are efilling then the code is slightly more complex
                 else
                     !First populate the lat points array