Merge pull request #32 from aselimov/ft--add-square-vacancy-disloop

Added rectangular vacancy disloop

src/opt_disl.f90

@@ -530,8 +530,8 @@ module opt_disl
     subroutine vacancy_disloop
         !This subroutine actually creates the dislocation loop.
 
-        real(kind=dp) :: neighbor_dis, temp_box(6), dis, normal_dim
-        integer :: i, j, index, delete_num, delete_index(atom_num)
+        real(kind=dp) :: neighbor_dis, temp_box(6), dis
+        integer :: i, j, index, delete_num, delete_index(atom_num), normal_dim
 
         neighbor_dis = HUGE(1.0_dp)
         index= 0
@@ -560,10 +560,10 @@ module opt_disl
         centroid = r_atom(:,index)
         !Now reset the list for the scanning algorithm
 
-        !Now scan over all atoms again and find the closest vloop_size number of atoms to the initial atom 
-        !that reside on the same plane.
-
         delete_num = 0
+        !Now scan over all atoms again and find the closest vloop_size number of atoms to the initial atom 
+        !that reside on the same plane. If loop_radius is > 0 then we build a circular vacancy cluster
+        if(loop_radius > 0) then 
             do i = 1, atom_num
                 !Check to see if it is on the same plane
                 if (in_block_bd(r_atom(:,i), temp_box)) then 
@@ -576,6 +576,25 @@ module opt_disl
                 end if
             end do 
 
+        !If the loop radius is < 0 we build a square vacancy cluster 
+        else if (loop_radius < 0) then 
+            !Set the new box boundaries
+            do i = 1, 3
+                if (i /= normal_dim) then 
+                    temp_box(2*i) = centroid(i) + abs(loop_radius)
+                    temp_box(2*i-1) = centroid(i) - abs(loop_radius)
+                end if
+            end do 
+
+            !Now delete all atoms with box box boundary which should create a square planar vacancy cluster
+            do i = 1, atom_num
+                !Check to see if the atom is in the bounding cube 
+                if (in_block_bd(r_atom(:,i), temp_box)) then 
+                    delete_num = delete_num + 1
+                    delete_index(delete_num) = i
+                end if
+            end do  
+        end if
         !Now delete the atoms
         call delete_atoms(delete_num, delete_index(1:delete_num))