Merge pull request #23 from aselimov/ft--opt-delete

Ft  opt delete

README.md

@@ -245,6 +245,24 @@ This allows the user to specify the boundary conditions for the model being outp
 
 **Example:** `-boundary psp`
 
+### Option Delete
+
+```
+-delete keywords
+```
+
+Delete requires the usage of additional keywords to specify which delete action will be taken. These additional keywords are below:
+
+**overlap**
+
+```
+-delete overlap rc_off
+```
+
+This command will delete all overlapping atoms within a specific cutoff radius `rc_off`. This currently does not affect elements. 
+
+****
+
 ## Position Specification
 
 Specifying positions in cacmb can be done through a variety of ways. Examples of each format is shown below.

src/Makefile

@@ -1,8 +1,8 @@
 FC=ifort
-FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
+FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays
-#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin
+#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays
 MODES=mode_create.o mode_merge.o mode_convert.o
-OPTIONS=opt_disl.o opt_group.o opt_orient.o
+OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o 
 OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o  $(MODES) $(OPTIONS) call_option.o
 
 .SUFFIXES:

src/call_option.f90

@@ -3,6 +3,7 @@ subroutine call_option(option, arg_pos)
     use opt_disl
     use opt_group
     use opt_orient
+    use opt_delete
     use box
     implicit none
 
@@ -27,6 +28,8 @@ subroutine call_option(option, arg_pos)
         arg_pos=arg_pos+1
         call get_command_argument(arg_pos, box_bc)
         arg_pos=arg_pos+1
+    case('-delete')
+        call run_delete(arg_pos)
     case default
         print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.'
         stop 3 

src/mode_merge.f90

@@ -130,8 +130,9 @@ module mode_merge
 
         integer :: i, ibasis, inod
         real(kind=dp), dimension(3) :: current_shift
+        character(len=3) :: alldims
 
-
+        alldims = 'xyz'
         !Calculate the current shift which is the filenum-1 multiplied by the user specified shift
         current_shift = (filenum-1)*shift_vec
 
@@ -154,11 +155,15 @@ module mode_merge
         !If we don't include the wrap command then we have to increase the size of the box
         if(.not.(wrap_flag)) then 
             do i = 1,3
+                if (alldims(i:i) /= trim(adjustl(dim))) then 
                     if (current_shift(i) < -lim_zero) then 
                         box_bd(2*i-1) = box_bd(2*i-1) - current_shift(i)
                     else if (current_shift(i) > lim_zero) then 
                         box_bd(2*i) = box_bd(2*i) + current_shift(i)
                     end if
+                else if (alldims(i:i) == trim(adjustl(dim))) then 
+                    box_bd(2*i) = box_bd(2*i) + current_shift(i)
+                end if
             end do
         end if
 

src/opt_delete.f90

@@ -0,0 +1,131 @@
+module opt_delete
+    
+    use parameters
+    use subroutines
+    use elements
+
+    implicit none 
+
+    real(kind=dp) :: rc_off
+
+    public  
+    contains
+
+    subroutine run_delete(arg_pos)  
+
+        integer, intent(inout) :: arg_pos
+
+        rc_off = -lim_zero
+        !Main calling function for delete option.
+        print *, '-----------------------Option Delete------------------------'
+
+        call parse_delete(arg_pos)
+
+        print *, atom_num
+
+        if (rc_off > 0.0_dp) call delete_overlap
+    end subroutine run_delete
+
+    subroutine parse_delete(arg_pos)
+        !Parse the delete command
+
+        integer, intent(inout) :: arg_pos
+
+        integer :: arg_len 
+        character(len=100) :: textholder
+        arg_pos = arg_pos + 1
+
+        call get_command_argument(arg_pos, textholder, arg_len)
+        if(arg_len==0) stop "Missing argument to delete command"
+
+        select case(textholder)
+            case('overlap')
+                arg_pos=arg_pos + 1
+                call get_command_argument(arg_pos, textholder, arg_len)
+                if(arg_len==0) stop "Missing argument to delete overlap command"
+                read(textholder, *) rc_off
+            case default
+                print *, "Command ", trim(adjustl(textholder)), " is not accepted for option delete"
+                stop 3
+        end select 
+
+        arg_pos = arg_pos + 1
+    end subroutine parse_delete
+
+    subroutine delete_overlap
+        !This subroutine deletes all overlapping atoms, which is defined as atoms which are separated by a distance of 
+        !less then rc_off
+
+        integer :: i, c(3), ci, cj, ck, num_nei, nei, delete_num
+        integer, dimension(atom_num) :: for_delete
+
+        !These are the variables containing the cell list information
+        integer, dimension(3) :: cell_num
+        integer, allocatable :: num_in_cell(:,:,:), which_cell(:,:)
+        integer, allocatable :: cell_list(:,:,:,:)
+        
+        allocate(which_cell(3,atom_num))
+
+        print *, atom_num
+        print *, for_delete(atom_num)
+        print *, which_cell(1,1)
+
+        !First pass the atom list and atom num to the algorithm which builds the cell list
+        call build_cell_list(atom_num, r_atom, rc_off, cell_num, num_in_cell, cell_list, which_cell)
+
+        !Now loop over every atom and figure out if it has neighbors within the rc_off
+        delete_num = 0
+        atom_loop: do i = 1, atom_num
+            
+            !c is the position of the cell that the atom belongs to 
+            c = which_cell(:,i)
+
+            !Check to make sure it hasn't already been deleted
+            if(all(c /= 0)) then 
+                !Now loop over all neighboring cells
+                do ci = -1, 1, 1
+                    do cj = -1, 1, 1
+                        do ck = -1, 1, 1  
+
+                            if (any((c + (/ ck, cj, ci /)) == 0)) cycle 
+
+                            if( (c(1) + ck > cell_num(1)).or.(c(2) + cj > cell_num(2)).or. &
+                                (c(3) + ci > cell_num(3))) cycle
+
+
+                            do num_nei = 1, num_in_cell(c(1) + ck, c(2) + cj, c(3) + ci)
+                                nei = cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci)
+
+                                !Check to make sure the atom isn't the same index as the atom we are checking 
+                                !and that the neighbor hasn't already been deleted
+                                if((nei /= i).and.(nei/= 0)) then 
+
+                                    !Now check to see if it is in the cutoff radius, if it is add it to the delete code
+                                    if (norm2(r_atom(:,nei)-r_atom(:,i)) < rc_off) then 
+                                        
+                                        delete_num = delete_num + 1
+                                        for_delete(delete_num) = max(i,nei)
+
+                                        !Now zero out the larger index
+                                        if(i > nei) then 
+                                            which_cell(:,i) = 0
+                                            cycle atom_loop
+                                        else
+                                            which_cell(:,nei) = 0 
+                                            cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci) = 0
+                                        end if
+                                    end if
+                                end if
+                            end do 
+                        end do
+                    end do 
+                end do 
+            end if
+
+        end do atom_loop
+
+        print *, "Overlap command deletes ", delete_num, " atoms"
+        !Now delete all the atoms
+        call delete_atoms(delete_num, for_delete(1:delete_num))
+    end subroutine delete_overlap
+end module opt_delete

src/subroutines.f90

@@ -289,4 +289,63 @@ module subroutines
         end do
     end subroutine
 
+    subroutine build_cell_list(numinlist, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell)
+        !This subroutine builds a cell list based on rc_off
+        
+        !----------------------------------------Input/output variables-------------------------------------------
+
+        integer, intent(in) :: numinlist !The number of points within r_list
+
+        real(kind=dp), dimension(3,numinlist), intent(in) :: r_list !List of points to be used for the construction of 
+                                                                    !the cell list.
+
+        real(kind=dp), intent(in) :: rc_off ! Cutoff radius which dictates the size of the cells
+
+        integer, dimension(3), intent(inout) :: cell_num !Number of cells in each dimension.
+
+        integer, allocatable, intent(inout) :: num_in_cell(:,:,:) !Number of points within each cell
+
+        integer, allocatable, intent(inout) :: cell_list(:,:,:,:) !Index of points from r_list within each cell.
+
+        integer, dimension(3,numinlist), intent(out) :: which_cell !The cell index for each point in r_list
+
+        !----------------------------------------Begin Subroutine -------------------------------------------
+
+        integer :: i, j, cell_lim, c(3)
+        real(kind=dp) :: box_len(3)
+        integer, allocatable :: resize_cell_list(:,:,:,:)
+
+        !First calculate the number of cells that we need in each dimension
+        do i = 1,3
+            box_len(i) = box_bd(2*i) - box_bd(2*i-1)
+            cell_num(i) = int(box_len(i)/(rc_off/2))+1
+        end do 
+
+        !Initialize/allocate variables 
+        cell_lim = 10
+        allocate(num_in_cell(cell_num(1),cell_num(2),cell_num(3)), cell_list(cell_lim, cell_num(1), cell_num(2), cell_num(3)))
+
+        !Now place points within cell
+        do i = 1, numinlist
+            !c is the position of the cell that the point belongs to 
+            do j = 1, 3
+                c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off/2)) + 1
+            end do
+
+            !Place the index in the correct position, growing if necessary
+            num_in_cell(c(1),c(2),c(3)) = num_in_cell(c(1),c(2),c(3)) + 1
+            if (num_in_cell(c(1),c(2),c(3)) > cell_lim) then 
+                allocate(resize_cell_list(cell_lim+10,cell_num(1),cell_num(2),cell_num(3)))
+                resize_cell_list(1:cell_lim, :, :, :) = cell_list
+                resize_cell_list(cell_lim+1:, :, :, :) = 0
+                call move_alloc(resize_cell_list, cell_list)
+            end if
+
+            cell_list(num_in_cell(c(1),c(2),c(3)),c(1),c(2),c(3)) = i
+            which_cell(:,i) = c
+        end do 
+
+        return
+    end subroutine build_cell_list
+
 end module subroutines