Working planar vacancy cluster algorithm

README.md

@@ -135,4 +135,18 @@ This options adds an arbitrarily oriented dislocation into your model based on u
 
 `char_angle` - Character angle of the dislocation (0 is screw and 90 is edge)
 
-`poisson` - Poisson's ratio used for the displacement field.
+`poisson` - Poisson's ratio used for the displacement field.
+
+### Option disloop
+
+````
+-disloop loop_normal loop_size x y z
+````
+
+This option deletes vacancies on a plane which when minimized should result in a dislocation loop structure. The arguments are below:
+
+`dim` - The box dimension which defines the normal to the loop plane. As of now this dimension must be a closed back direction, meaning that for fcc a box dimension has to be of the (111) family of planes
+
+`n` - The number of atoms to delete on the loop plane
+
+`x y z` - The centroid of the loop. As a note, 

src/Makefile

@@ -1,6 +1,6 @@
 FC=ifort
-#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
-FFLAGS=-mcmodel=large -Ofast -no-wrap-margin
+FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
+#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin
 MODES=mode_create.o mode_merge.o mode_convert.o
 OPTIONS=opt_disl.o
 OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o call_option.o $(MODES) $(OPTIONS)
@@ -19,12 +19,16 @@ clean:
 	$(RM) cacmb  *.o 
 
 testfuncs: testfuncs.o functions.o subroutines.o
-	$(FC) testfuncs.o functions.o subroutines.o elements.o -o $@
+	$(FC) testfuncs.o functions.o subroutines.o box.o elements.o -o $@
 
 .PHONY: cleantest
 cleantest:
 	$(RM) testfuncs testfuncs.o
 
+.PHONY: test
+test: testfuncs
+	./testfuncs
+
 $(OBJECTS) : parameters.o
 atoms.o subroutines.o testfuncs.o box.o : functions.o
 main.o io.o $(MODES) $(OPTIONS) : elements.o 

src/box.f90

@@ -85,4 +85,21 @@ module box
         return
     end subroutine grow_box
 
+
+    subroutine in_sub_box(r, which_sub_box)
+        !This returns which sub_box a point is in. It returns the first sub_box with boundaries which 
+        !contain the point.
+        real(kind=dp), dimension(3), intent(in) :: r
+        integer, intent(out) :: which_sub_box
+
+        integer :: i 
+
+        do i = 1, sub_box_num
+            if( in_block_bd(r, sub_box_bd(:,i))) then 
+                which_sub_box = i
+                exit
+            end if
+        end do 
+    return 
+    end subroutine in_sub_box
 end module box

src/elements.f90

@@ -3,6 +3,7 @@ module elements
     !This module contains the elements data structures, structures needed for building regions
     !and operations that are done on elements
     use parameters
+    use functions
     use subroutines
     implicit none
 
@@ -355,4 +356,26 @@ module elements
         return
     end subroutine rhombshape
 
+    subroutine delete_atoms(num, index)
+        !This subroutine deletes atoms from the atom arrays 
+        integer, intent(in) :: num
+        integer, intent(inout), dimension(num) :: index
+
+        integer :: i, j
+
+        call heapsort(index)
+
+        !We go from largest index to smallest index just to make sure that we don't miss 
+        !accidentally overwrite values which need to be deleted
+        do i = num, 1, -1
+            if(index(i) == atom_num) then 
+                r_atom(:,index(i)) = 0.0_dp
+                type_atom(index(i)) = 0
+            else
+                r_atom(:,index(i)) = r_atom(:, atom_num)
+                type_atom(index(i)) = type_atom(atom_num)
+            end if            
+            atom_num = atom_num - 1
+        end do 
+    end subroutine
 end module elements

src/functions.f90

@@ -242,4 +242,12 @@ END FUNCTION StrDnCase
         return
     end function unitvec
 
+    pure function norm_dis(rl,rk)
+        !This just returns the magnitude of the vector between two points
+        real(kind=dp), dimension(3), intent(in) :: rl, rk
+        real(kind=dp) :: norm_dis(4)
+
+        norm_dis(1:3) = (rk - rl)
+        norm_dis(4) = norm2(rk-rl)
+    end function
 end module functions

src/subroutines.f90

@@ -203,4 +203,50 @@ module subroutines
         end if
     end subroutine parse_pos
 
+    subroutine heapsort(a)
+     
+       integer, intent(in out) :: a(0:)
+       integer :: start, n, bottom
+       real :: temp
+     
+       n = size(a)
+       do start = (n - 2) / 2, 0, -1
+         call siftdown(a, start, n);
+       end do
+     
+       do bottom = n - 1, 1, -1
+         temp = a(0)
+         a(0) = a(bottom)
+         a(bottom) = temp;
+         call siftdown(a, 0, bottom)
+       end do
+     
+    end subroutine heapsort
+     
+    subroutine siftdown(a, start, bottom)
+     
+      integer, intent(in out) :: a(0:)
+      integer, intent(in) :: start, bottom
+      integer :: child, root
+      real :: temp
+     
+      root = start
+      do while(root*2 + 1 < bottom)
+        child = root * 2 + 1
+     
+        if (child + 1 < bottom) then
+          if (a(child) < a(child+1)) child = child + 1
+        end if
+     
+        if (a(root) < a(child)) then
+          temp = a(child)
+          a(child) = a (root)
+          a(root) = temp
+          root = child
+        else
+          return
+        end if  
+      end do      
+     
+    end subroutine siftdown
 end module subroutines