From c698c31ede5a5d90ce9ccf62665a4c503a1a439e Mon Sep 17 00:00:00 2001
From: Alex Selimov <aselimov3@gatech.edu>
Date: Mon, 24 Feb 2020 14:23:32 -0500
Subject: [PATCH] First version that doesn't crash, incorrectly places elements
 and atoms

---
 src/opt_group.f90 | 51 ++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 39 insertions(+), 12 deletions(-)

diff --git a/src/opt_group.f90 b/src/opt_group.f90
index cfd6499..8ff789a 100644
--- a/src/opt_group.f90
+++ b/src/opt_group.f90
@@ -263,7 +263,7 @@ module opt_group
 
         integer :: i, j, ix, iy, iz, inod, ibasis, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num,   &
                    current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3)
-        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3),  ori_inv(3,3),  &
+        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3),  ori_inv(3,3),  r(3), &
                          r_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8)
         logical, allocatable :: lat_points(:,:,:)
         character(len=100) :: remesh_ele_type
@@ -295,8 +295,8 @@ module opt_group
 
                 group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/remesh_lat_pam))
                 do i = 1, 3
-                    bd_in_lat(2*i-1) = minval(group_in_lat(i,:))
-                    bd_in_lat(2*i) = maxval(group_in_lat(i,:))
+                    bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:)))
+                    bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:)))
                 end do
 
         end select 
@@ -306,8 +306,11 @@ module opt_group
 
         !Now place all group atoms and group interpolated atoms into lat_points
         do i = 1, group_atom_num
-            r_lat = r_atom(:,atom_index(i))/remesh_lat_pam
-            r_lat = matmul(fcc_inv,matmul(ori_inv,r_lat)) 
+            r = r_atom(:,atom_index(i))/remesh_lat_pam
+            r = matmul(fcc_inv,matmul(ori_inv,r))
+            do j = 1, 3
+                r_lat(j) = nint(r(j))
+            end do 
             !Do a check to make sure the code is working and that lattice points aren't being written on top of each other. 
             !This is primarily a debugging statement
             if((r_lat(1)==13).and.(r_lat(2)==13).and.(r_lat(3)==-13)) then 
@@ -349,10 +352,20 @@ module opt_group
         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
 
-        do iz = 1, bd_in_array(3)
-            do iy = 1, bd_in_array(2)
-                do ix = 1, bd_in_array(1)
+        zloop: do iz = 1, bd_in_array(3)
+            yloop: do iy = 1, bd_in_array(2)
+                xloop: do ix = 1, bd_in_array(1)
                     if (lat_points(ix, iy,iz)) then 
+                        
+                        !Check to see if the element overshoots the bound
+                        if (iz+remesh_size-1 > bd_in_array(3)) then
+                            exit zloop
+                        else if (iy+remesh_size-1 > bd_in_array(2)) then
+                            cycle zloop
+                        else if (ix+remesh_size-1 > bd_in_array(1)) then 
+                            cycle yloop
+                        end if  
+
                         if (all(lat_points(ix:ix+remesh_size-1,iy:iy+remesh_size-1,iz:iz+remesh_size-1))) then 
                             do inod = 1, 8
                                 vlat = ele(:,inod) + (/ix, iy, iz /)
@@ -369,11 +382,25 @@ module opt_group
 
                         end if
                     end if
-                end do 
-            end do 
-        end do 
-
+                end do xloop
+            end do yloop
+        end do zloop
 
+        !Now we have to add any leftover lattice points as atoms
+        do iz = 1, bd_in_array(3)
+            do iy=1, bd_in_array(2)
+                do ix = 1, bd_in_array(1)
+                    if(lat_points(ix,iy,iz)) then 
+                        vlat = (/ ix, iy, iz /)
+                        do i = 1, 3
+                            vlat = vlat + bd_in_lat(2*i-1)-5
+                        end do 
+                       r = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam
+                       call add_atom(remesh_type, sbox_atom(atom_index(1)), r)
+                    end if
+                end do
+            end do 
+        end do
         
     end subroutine remesh_group