From 338587b7b74cddbcdeb94546cf304ec68ccc16c5 Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Thu, 4 Jun 2020 12:13:42 -0400
Subject: [PATCH] Added sphere as a possible group shape

---
 README.md         |  5 +++++
 src/opt_group.f90 | 40 ++++++++++++++++++++++++++++++++--------
 2 files changed, 37 insertions(+), 8 deletions(-)

diff --git a/README.md b/README.md
index 4c4e5ed..d285495 100644
--- a/README.md
+++ b/README.md
@@ -218,6 +218,11 @@ This shape is similar to a wedge shape except instead of becoming atomically sha
 `bw` - Base width
 `tr` - Tip radius
 
+*Sphere*
+
+`-group atoms sphere x y z r`
+This shape selects all atoms within a sphere centered at `(x,y,z)` with radius `r`.
+
 **Displace:**
 
 ```
diff --git a/src/opt_group.f90 b/src/opt_group.f90
index 080d0e6..281ab28 100644
--- a/src/opt_group.f90
+++ b/src/opt_group.f90
@@ -10,7 +10,7 @@ module opt_group
 
     integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num
     character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape
-    real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), tip_radius, bwidth
+    real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth
     logical :: displace, delete, max_remesh, refine, group_nodes, flip
 
     integer, allocatable :: element_index(:), atom_index(:) 
@@ -171,7 +171,7 @@ module opt_group
             arg_pos  = arg_pos + 1
             call get_command_argument(arg_pos, textholder, arglen)
             if (arglen==0) STOP "Missing tip radius in group notch command"
-            read(textholder,*)  tip_radius
+            read(textholder,*)  radius
 
             !Calculate the vertex positions
             vertices(:,1) = centroid
@@ -201,7 +201,7 @@ module opt_group
              
             !Now update the centroid so that the desire tip diameter matches the wedge with
             if (bwidth > 0) then 
-                centroid(dim1) = centroid(dim1) + 2*tip_radius*(H)/bwidth
+                centroid(dim1) = centroid(dim1) + 2*radius*(H)/bwidth
             end if
         !Read the ID type shape for group 
         case('id')
@@ -313,6 +313,21 @@ module opt_group
                     if(i ==1) arg_pos = arg_pos + 1
                 end do
             end select
+
+        case('sphere')
+            !First extract the sphere centroid
+            do i = 1, 3
+                arg_pos = arg_pos + 1
+                call get_command_argument(arg_pos, textholder, arglen)
+                if (arglen==0) STOP "Missing sphere centroid in group command"
+                call parse_pos(i, textholder, centroid(i))   
+            end do
+            !Now get the tip radius
+            arg_pos=arg_pos+1
+            call get_command_argument(arg_pos, textholder, arglen)
+            if (arglen==0) STOP "Missing sphere radius in group command"
+            read(textholder, *) radius
+
         case default
             print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", &
                      "for accepted group shapes."
@@ -376,11 +391,13 @@ module opt_group
         case ('wedge')
             print *, "Group has wedge shape with dim1", dim1, "and dim2", dim2, "and vertices ", vertices
         case ('notch')
-            print *, "Group has notch shape with dim1", dim1, "and dim2", dim2, " tip radius ", tip_radius,  "and vertices ", &
+            print *, "Group has notch shape with dim1", dim1, "and dim2", dim2, " tip radius ", radius,  "and vertices ", &
                       vertices
         case('id')
             print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms."
             return
+        case('sphere')
+            print *, "Group has sphere shape with centroid ", centroid(:), " and radius ", radius
         end select
 
         !Reset group if needed
@@ -837,7 +854,7 @@ module opt_group
     function in_group(r)
         !This subroutine determines if a point is within the group boundaries
         real(kind=dp), intent(in) :: r(3)
-        real(kind=dp) :: r_notch
+        real(kind=dp) :: rnorm
 
         integer :: dim3, i
         logical :: in_group
@@ -857,9 +874,16 @@ module opt_group
                 if (r(dim1) < centroid(dim1)) in_group = .false.
             end if
 
-            r_notch = sqrt((r(dim1) - centroid(dim1))**2 + (r(dim3)-centroid(dim3))**2)
-            in_group = in_group.or.(r_notch < tip_radius)
+            rnorm = sqrt((r(dim1) - centroid(dim1))**2 + (r(dim3)-centroid(dim3))**2)
+            in_group = in_group.or.(rnorm < radius)
 
-                   end select 
+        case('sphere') 
+            rnorm = norm2(r(:) - centroid(:))
+            if (rnorm <= radius) then 
+                in_group = .true. 
+            else
+                in_group = .false.
+            end if
+        end select 
     end function in_group
 end module opt_group