Final changes fixing remesh group

src/opt_group.f90

@@ -302,12 +302,12 @@ module opt_group
         !This command is used to remesh the group to a desired element size
 
         integer :: i, j, k, 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), old_ele, old_atom, &
+                   current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3), new_ele, new_atom, &
-                   max_loops, working_esize, group_lat_num, lat_list(10), group_sbox_num, sbox_list(100), is, ilat
+                   max_loops, working_esize, group_lat_num, lat_list(10), group_sbox_num, sbox_list(100), is, ilat, tot_dof
 
         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), group_bd(6)
-        logical, allocatable :: lat_points(:,:,:)
+        logical, allocatable :: lat_points(:,:,:), new_sbox, new_lat
         character(len=100) :: remesh_ele_type
 
 
@@ -318,35 +318,67 @@ module opt_group
         group_sbox_num = 0
         group_lat_num = 0
         do i = 1, group_atom_num 
+            new_sbox=.true.
+            new_lat=.true.
             do j = 1, group_sbox_num
-                if (sbox_list(j) == sbox_atom(atom_index(i))) exit
+                if (sbox_list(j) == sbox_atom(atom_index(i))) then 
+                    new_sbox=.false.
+                    exit
+                end if
+            end do
+
+            if(new_sbox) then 
                 group_sbox_num = group_sbox_num + 1 
                 sbox_list(group_sbox_num) = sbox_atom(atom_index(i))
-            end do
+            end if
 
             do j = 1, group_lat_num
-                if (basis_type(1,lat_list(j)) == type_atom(atom_index(i))) exit
+                if (basis_type(1,lat_list(j)) == type_atom(atom_index(i))) then 
+                    new_lat = .false.
+                    exit
+                end if
+            end do 
+
+            if (new_lat) then 
                 group_lat_num = group_lat_num + 1 
                 do k = 1, lattice_types
-                    if (basis_type(1,k) == type_atom(atom_index(i))) lat_list(group_lat_num) = k
+                    if (basis_type(1,k) == type_atom(atom_index(i)))  lat_list(group_lat_num) = k
                 end do
-            end do 
+            end if
+
         end do 
 
         do i = 1, group_ele_num
+            new_sbox=.true.
+            new_lat = .true.
             do j = 1, group_sbox_num
-                if (sbox_list(j) == sbox_ele(element_index(i))) exit
+                if (sbox_list(j) == sbox_ele(element_index(i))) then 
-                group_sbox_num = group_sbox_num + 1
+                    new_sbox = .false.
-                sbox_list(group_sbox_num) = sbox_ele(element_index(i))
+                    exit
+                end if
             end do
             
+            if (new_sbox) then
+                 group_sbox_num = group_sbox_num + 1
+                sbox_list(group_sbox_num) = sbox_ele(element_index(i))
+            end if
+
             do j = 1, group_lat_num
-                if (lat_list(group_lat_num) == lat_ele(element_index(i))) exit
+                if (lat_list(group_lat_num) == lat_ele(element_index(i))) then 
+                    new_lat=.false.
+                    exit
+                end if
+            end do
+
+            if (new_lat) then
                 group_lat_num = group_lat_num + 1
                 lat_list(group_lat_num) = lat_ele(element_index(i))
-            end do
+            end if
         end do
 
+        new_atom = 0
+        new_ele=0
+        tot_dof=0
         do is = 1, group_sbox_num 
 
             orient = sub_box_ori(:, :, sbox_list(is))
@@ -388,6 +420,7 @@ module opt_group
 
                 !If for some reason there are no dof in this loop then cycle out
                 if(dof == 0) cycle
+                tot_dof = tot_dof+dof
 
                 group_in_lat = reshape((/ group_bd(1),group_bd(3),group_bd(5), &
                                           group_bd(2),group_bd(3),group_bd(5), &
@@ -410,6 +443,7 @@ module opt_group
 
                 !Now place all group atoms and group interpolated atoms into lat_points
                 do i = 1, group_atom_num
+                    if (.not.((sbox_atom(atom_index(i)) == is).and.(basis_type(1,ilat) == type_atom(atom_index(i))))) cycle
                     r = r_atom(:,atom_index(i))/lapa(ilat)
                     r = matmul(fcc_inv,matmul(ori_inv,r))
                     do j = 1, 3
@@ -427,6 +461,7 @@ module opt_group
 
                 !Now place interpolated atoms within lat_points array
                 do i =1, group_ele_num
+                    if (.not.((sbox_ele(element_index(i)) == is).and.( lat_ele(element_index(i)) == ilat))) cycle
                     ie = element_index(i)
                     call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
                     do j = 1, size_ele(ie)**3 * basisnum(lat_ele(ie))
@@ -446,17 +481,6 @@ module opt_group
                     end do 
                 end do 
 
-                print *, "Group has ", dof, " degrees of freedom to remesh"
-
-                !Delete all elements and atoms to make space for new elements and atoms
-                call delete_atoms(group_atom_num, atom_index)
-                call delete_elements(group_ele_num, element_index)
-
-                old_atom = atom_num
-                old_ele = ele_num
-
-
-
                 !Now run remeshing algorithm, not the most optimized or efficient but gets the job done
                 !Figure out new looping boundaries
                 bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10
@@ -499,6 +523,7 @@ module opt_group
                                         lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1) = .false.
                                         !Add the element, for the sbox we just set it to the same sbox that we get the orientation 
                                         !from. In this case it is from the sbox of the first atom in the group.
+                                        new_ele = new_ele+1
                                         call add_element(remesh_ele_type, working_esize, ilat, sbox_atom(atom_index(1)),r_new_node)
 
                                     end if
@@ -519,15 +544,22 @@ module opt_group
                                 end do 
                                 lat_points(ix,iy,iz) = .false.
                                r = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat)
-                               call add_atom(basis_type(1,ilat), sbox_atom(atom_index(1)), r)
+                               new_atom=new_atom+1
+                               call add_atom(basis_type(1,ilat), is, r)
                             end if
                         end do
                     end do 
                 end do
+                deallocate(lat_points)
             end do
         end do
 
-        print *, "Remeshing has created ", ele_num-old_ele, " elements and ", atom_num-old_atom, " atoms."
+        !Delete all elements and atoms to make space for new elements and atoms
+        call delete_atoms(group_atom_num, atom_index)
+        call delete_elements(group_ele_num, element_index)
+    
+        print *, tot_dof, " degrees of freedom in group"
+        print *, "remesh_group creates ", new_atom, " atoms and ", new_ele, " elements."
     end subroutine remesh_group
 
     subroutine delete_group