Manual merge of control box code

README.md

@@ -129,6 +129,11 @@ wrap
 
 This will wrap atomic positions back inside the box. Effectively as if periodic boundary conditions are applied so that atoms which exit from one side of the simulation cell enter back in through the other.
 
+###Mode Metric
+```
+--metric cmetric nfiles {filepaths} 
+```
+
 ## Options
 
 Options are additional portions of code which have additional functionality. Options are performed in the order that they appear in the argument list and can be added to any mode. If wanting to use strictly options use `--convert` to specify input and output files.

src/Makefile

@@ -8,9 +8,9 @@ FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -w
 #FFLAGS=-mcmodel=large -O3 -g
 #FFLAGS=-mcmodel=large -O0 -g -fbacktrace -fcheck=all
 
-MODES=mode_create.o mode_merge.o mode_convert.o
+MODES=mode_create.o mode_merge.o mode_convert.o mode_metric.o
 OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o opt_deform.o opt_redef_box.o opt_slip_plane.o
-OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o  $(MODES) $(OPTIONS) call_option.o sorts.o
+OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o neighbors.o $(MODES) $(OPTIONS) call_option.o sorts.o
 
 .SUFFIXES:
 .SUFFIXES: .c .f .f90 .F90 .o
@@ -49,3 +49,4 @@ elements.o : sorts.o
 $(MODES) $(OPTIONS) subroutines.o io.o : atoms.o box.o
 $(MODES) main.o : io.o
 testfuncs.o elements.o mode_create.o $(OPTIONS) $(MODES): subroutines.o
+opt_delete.o mode_metric.o : neighbors.o

src/box.f90

@@ -97,4 +97,10 @@ module box
         end do 
     return 
     end subroutine in_sub_box
+
+    subroutine reset_box
+        !This subroutine just resets the box boundary and position
+        box_bc = "ppp"
+        box_bd(:) = 0
+    end subroutine reset_box
 end module box

src/call_mode.f90

@@ -1,16 +1,16 @@
-subroutine call_mode(arg_pos,mode)
+subroutine call_mode(arg_pos)
     !This code is used to parse the command line argument for the mode information and calls the required 
     !mode module.
 
     use mode_create
     use mode_convert
     use mode_merge
+    use mode_metric
     use parameters
 
     implicit none    
 
     integer, intent(out) :: arg_pos
-    character(len=100), intent(in) :: mode
 
     select case(mode)
     case('--create')
@@ -19,6 +19,8 @@ subroutine call_mode(arg_pos,mode)
         call convert(arg_pos)  
     case('--merge')
         call merge(arg_pos) 
+    case('--metric')
+        call metric(arg_pos)
     case default
         print *, "Mode ", trim(adjustl(mode)), " currently not accepted. Please check documentation for ", &
                  "accepted modes and rerun."

src/call_option.f90

@@ -40,6 +40,8 @@ subroutine call_option(option, arg_pos)
         call run_delete(arg_pos)
     case('-set_cac')
         arg_pos=arg_pos +3
+    case('-set_types')
+        arg_pos = arg_pos + 3 + atom_types
     case('-redef_box')
         call redef_box(arg_pos)
     case('-slip_plane')

src/elements.f90

@@ -13,6 +13,8 @@ module elements
     character(len=100), allocatable :: type_ele(:) !Element type 
     integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:), tag_ele(:) !Element size
     real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
+    !Element result data structures
+    real(kind=8), allocatable :: force_node(:,:,:,:), virial_node(:,:,:,:), energy_node(:,:,:)
 
     integer, save :: ele_num !Number of elements
     integer, save :: node_num !Total number of nodes
@@ -22,6 +24,8 @@ module elements
     integer, allocatable :: sbox_atom(:), tag_atom(:)
     real(kind =dp),allocatable :: r_atom(:,:) !atom position
     integer :: atom_num=0 !Number of atoms
+    !Atom result data structures information
+    real(kind=8), allocatable :: force_atom(:,:), virial_atom(:,:), energy_atom(:)
 
     !Mapping atom type to provided name
     character(len=2), dimension(10) :: type_to_name
@@ -512,6 +516,7 @@ module elements
         !We go from largest index to smallest index just to make sure that we don't miss 
         !accidentally overwrite values which need to be deleted
         do i = num, 1, -1
+            node_num = node_num - ng_node(lat_ele(sorted_index(i)))
             if(sorted_index(i) == ele_num) then 
                 r_node(:,:,:,sorted_index(i)) = 0.0_dp
                 type_ele(sorted_index(i)) =''
@@ -520,7 +525,6 @@ module elements
                 sbox_ele(sorted_index(i)) = 0
                 tag_ele(sorted_index(i)) = 0
             else
-                node_num = node_num - ng_node(lat_ele(sorted_index(i)))
                 r_node(:,:,:,sorted_index(i)) = r_node(:,:,:,ele_num)
                 type_ele(sorted_index(i)) = type_ele(ele_num)
                 size_ele(sorted_index(i)) = size_ele(ele_num)
@@ -549,7 +553,6 @@ module elements
 
         max_bd(:) = -huge(1.0_dp)
         min_bd(:) = huge(1.0_dp)
-
 do i = 1, atom_num
             do j = 1, 3
                 if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + tol 
@@ -775,4 +778,71 @@ module elements
 
     end subroutine 
 
+    subroutine alloc_dat_arrays(n,m)
+        !This subroutine used to provide initial allocation for the atom and element data arrays
+        integer, intent(in) :: n,m !n-size of element arrays, m-size of atom arrays
+        integer :: allostat
+
+        !Allocate element arrays
+        if (n > 0) then 
+            if (allocated(force_node)) then 
+                deallocate(force_node, virial_node, energy_node)
+            end if 
+            allocate(force_node(3,max_basisnum, max_ng_node, n), &
+                     virial_node(6,max_basisnum, max_ng_node, n), &
+                     energy_node(max_basisnum,max_ng_node,n),    &
+                     stat=allostat)
+            if(allostat > 0) then 
+                print *, "Error allocating element data arrays in mode_metric because of:", allostat
+                stop
+            end if
+
+        end if
+
+        if (m > 0) then
+            if (allocated(force_atom)) then 
+                deallocate(force_atom, virial_atom, energy_atom)
+            end if 
+            allocate(force_atom(3, m), &
+                     virial_atom(6, m), &
+                     energy_atom(m),    &
+                     stat=allostat)
+            if(allostat > 0) then 
+                print *, "Error allocating atom data arrays in mode_metric because of:", allostat
+                stop
+            end if
+        end if
+
+    end subroutine 
+
+    subroutine add_atom_data(ia, eng, force, virial)
+        !Function which sets the atom data arrays 
+        integer, intent(in) :: ia
+        real(kind=dp), intent(in) :: eng, force(3), virial(6)
+
+        energy_atom(ia) = eng
+        force_atom(:,ia) = force(:) 
+        virial_atom(:,ia) = virial(:) 
+        return
+    end subroutine add_atom_data
+
+    subroutine add_element_data(ie, eng, force, virial)
+        !Function which sets the element data arrays
+        integer, intent(in) :: ie
+        real(kind=dp), intent(in) :: eng(max_basisnum, max_ng_node), &
+                                     force(3,max_basisnum, max_ng_node), &
+                                     virial(6,max_basisnum,max_ng_node)
+        energy_node(:,:,ie)   = eng
+        force_node(:,:,:,ie)  = force
+        virial_node(:,:,:,ie) = virial
+        return
+    end subroutine add_element_data
+
+    subroutine reset_data
+        !This function resets all of the data arrays for the elements and atoms
+        atom_num = 0
+        ele_num = 0
+        node_num = 0
+    end subroutine reset_data
+
 end module elements

src/functions.f90

@@ -271,4 +271,124 @@ END FUNCTION StrDnCase
         norm_dis(1:3) = (rk - rl)
         norm_dis(4) = norm2(rk-rl)
     end function
+
+    pure function matinv3(A) result(B)
+        !! Performs a direct calculation of the inverse of a 3×3 matrix.
+        real(kind=dp), intent(in) :: A(3,3)   !! Matrix
+        real(kind=dp)             :: B(3,3)   !! Inverse matrix
+        real(kind=dp)             :: detinv
+
+        if(abs(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)&
+                  - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)&
+                  + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) < lim_zero) then 
+            B(:,:) = 0
+            return
+        else
+            ! Calculate the inverse determinant of the matrix
+            
+            detinv = 1/(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)&
+                      - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)&
+                      + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1))
+
+            ! Calculate the inverse of the matrix
+            B(1,1) = +detinv * (A(2,2)*A(3,3) - A(2,3)*A(3,2))
+            B(2,1) = -detinv * (A(2,1)*A(3,3) - A(2,3)*A(3,1))
+            B(3,1) = +detinv * (A(2,1)*A(3,2) - A(2,2)*A(3,1))
+            B(1,2) = -detinv * (A(1,2)*A(3,3) - A(1,3)*A(3,2))
+            B(2,2) = +detinv * (A(1,1)*A(3,3) - A(1,3)*A(3,1))
+            B(3,2) = -detinv * (A(1,1)*A(3,2) - A(1,2)*A(3,1))
+            B(1,3) = +detinv * (A(1,2)*A(2,3) - A(1,3)*A(2,2))
+            B(2,3) = -detinv * (A(1,1)*A(2,3) - A(1,3)*A(2,1))
+            B(3,3) = +detinv * (A(1,1)*A(2,2) - A(1,2)*A(2,1))
+        end if
+    end function
+
+     pure function transpose3(A) result(B)
+         !!Transposes matrix A
+         real(kind=dp), intent(in) :: A(3,3)
+         real(kind=dp)             :: B(3,3)  
+
+         integer :: i, j
+         forall(i =1:3,j=1:3) B(i,j) = A(j,i)
+
+    end function transpose3
+
+    pure function sqrt3(A) result(B)
+        !This calculates the square of matrix A fulfilling the equation B*B = A according to 
+        !the algorithm by Franca1989
+
+        real(kind=dp), intent(in) :: A(3,3)
+        real(kind=dp)             :: B(3,3)
+
+        real(kind=dp) :: Ione, Itwo, Ithree, l, k, phi, Asq(3,3), lambda, Bone, Btwo, Bthree, p
+    
+        !Step 1 is calculating the invariants of C
+        Ione = A(1,1) + A(2,2)  + A(3,3)
+        Asq = matmul(A,A)
+        Itwo = 0.5_dp *(Ione*Ione - (Asq(1,1) + Asq(2,2) + Asq(3,3)))
+        Ithree = (A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)& 
+                - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)& 
+                + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) 
+
+        if (Ithree < 0) then 
+            B(:,:)=0.0_dp
+            return
+        end if
+        !Check for an isotropic matrix
+        k = Ione*Ione - 3*Itwo
+        if (k < lim_zero) then 
+            lambda = sqrt(Ione/3.0_dp)
+            B = lambda*identity_mat(3)
+        else
+            l = Ione*(Ione*Ione - 9.0_dp/2.0_dp * Itwo) + 27.0_dp/2.0_dp * Ithree
+            p = l/(k**(1.5_dp))
+
+            if (p > 1.0 ) then 
+                B(:,:) = 0.0_dp
+                return
+            end if
+
+            if ((p< -1).or.(p>1)) then 
+                B(:,:) = 0.0_dp
+                return
+            end if
+            phi = acos(p)
+            lambda = sqrt(1.0_dp/3.0_dp * (Ione + 2*sqrt(k)*cos(phi/3)))
+            
+            !Now calculate invariantes of B
+            Bthree = sqrt(Ithree)
+            if((-lambda*lambda + Ione + 2*Ithree/lambda) < 0) then 
+                B(:,:) = 0.0_dp
+                return
+            end if
+            Bone = lambda + sqrt(- lambda*lambda + Ione + 2*Ithree/lambda)
+            Btwo = (Bone*Bone - Ione)/2.0_dp
+            
+            !Now calculate B
+            if(abs(Bone*Btwo -Bthree) < lim_zero) then 
+                B(:,:) = 0.0_dp
+                return
+            end if
+            B = 1/(Bone*Btwo - Bthree) *(Bone*Bthree*identity_mat(3) + (Bone*Bone - Btwo)*A - matmul(A,A))
+        end if
+    end function sqrt3
+
+    pure function permutation(i,j,k) result(e)
+        !Calculates the permutation tensor
+        integer, intent(in) :: i,j,k
+        integer             :: e
+        
+        if ( ((i==1).and.(j==2).and.(k==3)).or. &
+             ((i==2).and.(j==3).and.(k==1)).or. &
+             ((i==3).and.(j==1).and.(k==2))) then
+            e=1
+        else if( ((i==2).and.(j==1).and.(k==3)).or. &
+                 ((i==1).and.(j==3).and.(k==2)).or. &
+                 ((i==3).and.(j==2).and.(k==1))) then
+            e=-1
+        else
+            e=0
+        end if
+    end function permutation
+
 end module functions

src/io.f90

@@ -106,7 +106,7 @@ module io
                 call write_lmpcac(outfiles(i))
             case default
                 print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
-                         " is not accepted for writing. Please select from: xyz and try again"
+                         " is not accepted for writing. Please select from: xyz, lmp, vtk, mb, restart, cac and try again"
                 stop 
 
             end select
@@ -583,34 +583,41 @@ module io
             temp_infile = filename
         end if
 
-        !Infinite loop which only exists if user provides valid filetype
-        do while(.true.)
-
         !Check to see if file exists, if it does then ask user if they would like to overwrite the file
         inquire(file=trim(temp_infile), exist=file_exists)
         if (.not.file_exists) then 
-                print *, "The file ", trim(adjustl(filename)), " does not exist. Please input a filename that exists"
+            print *, "The file ", trim(adjustl(filename)), " does not exist. Please input an existing file to read."
-                read(*,*) temp_infile
+            stop 3
-                cycle
         end if
 
         select case(temp_infile(scan(temp_infile,'.',.true.)+1:))
         case('restart', 'mb', 'cac')
             infilenum=infilenum+1
             infiles(infilenum) =  temp_infile
-                exit
+        case('out')
+            if(atom_types == 0) then 
+                print *, "Please run -set_types command prior to running code requiring reading in pycac_*.out files"
+                stop 3
+            end if
+            select case(trim(adjustl(mode)))
+            case('--convert','--metric')
+                infilenum = infilenum+1
+                infiles(infilenum) = temp_infile
+            case default
+                print *, "Files of type .out cannot be used with mode ", trim(adjustl(mode))
+                stop 3
+            end select
+
         case default    
             print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), "not currently accepted. ", &
-                         "please input a filename with extension from following list: mb, restart."
+                     "please input a filename with extension from following list: mb, restart, cac, or out."
-                read(*,*) temp_infile
+            stop 3
-
         end select
-        end do
 
     end subroutine get_in_file 
 
     subroutine read_in(i, displace, temp_box_bd)
-        !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine
+        !This subroutine reads in file i
 
         integer, intent(in) :: i 
         real(kind=dp), dimension(3), intent(in) :: displace
@@ -624,9 +631,11 @@ module io
             call read_pycac(infiles(i), displace, temp_box_bd)
         case('cac')
             call read_lmpcac(infiles(i), displace, temp_box_bd)
+        case('out')
+            call read_pycac_out(infiles(i), displace, temp_box_bd)
         case default
             print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
-                     " is not accepted for writing. Please select from: mb and try again"
+                     " is not accepted for reading. Please select from: mb,restart,cac,out and try again"
             stop 
 
         end select
@@ -960,6 +969,105 @@ module io
         end if
     end subroutine read_pycac
 
+    subroutine read_pycac_out(file, displace, temp_box_bd)
+        !This subroutine reads in the pyCAC dump file
+    
+
+        !Arguments
+        character(len=100), intent(in) :: file
+        real(kind=dp), dimension(3), intent(in) :: displace
+        real(kind=dp), dimension(6), intent(out) :: temp_box_bd
+
+        !Internal Variables
+        integer :: i, in_eles, in_atoms, inbtypes(10), lat_type, ia, ie, inod, &
+            id, type_node, ilat, esize, tag, type
+        real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), &
+                         ee(1,8), fe(3,1,8), ve(3,1,8), re(3,1,8)
+        character(len=100) :: textholder, fcc
+        character(len=1000) :: line
+
+
+        open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') 
+
+        !Now initialize some important variables if they aren't defined
+        if (max_basisnum==0) max_basisnum = 1
+        if (max_ng_node==0) max_ng_node=8
+        fcc="fcc"
+
+        !Skip header comment lines
+        read(11, *) textholder
+        read(11, *) textholder
+
+        !Read atom number and element number and grow element arrays by needed amount
+        read(11,*) textholder, in_atoms, textholder, in_eles
+        call grow_ele_arrays(in_eles, in_atoms)
+        call alloc_dat_arrays(in_eles, in_atoms)
+
+        !Read boundary information
+        read(11,*) textholder, box_bc(1:1), box_bc(2:2), box_bc(3:3), temp_box_bd(:)
+        do i = 1, 3
+            if (abs(displace(i)) > lim_zero) then 
+                newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
+            else
+                newdisplace(i)=displace(i)
+            end if 
+            temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i)
+            temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i)
+        end do
+
+        call grow_box(temp_box_bd)
+
+        !Allocate sub_box boundaries
+        if (sub_box_num == 0) then 
+            call alloc_sub_box(1)
+        else
+            call grow_sub_box(1)
+        end if
+
+        !Because orientations and other needed sub_box information isn't really 
+        !present within the .cac file we just default a lot of it.
+        sub_box_ori(:,:,sub_box_num+1) = identity_mat(3)
+        sub_box_bd(:, sub_box_num+1) = temp_box_bd
+        sub_box_num = sub_box_num + 1
+
+        if(in_atoms > 0 ) then
+            !Read atom header
+            read(11,*) textholder
+            do ia = 1, in_atoms
+                read(11,'(a)')  line(:)
+                read(line,*) tag, type, ra(:), ea, fa(:), va(:)         
+                call add_atom(tag, type, sub_box_num, ra)
+                call add_atom_data(atom_num, ea, fa, va)
+            end do
+
+        end if
+
+        if(in_eles > 0) then 
+            !Add the lattice_types based on the atom types
+            inbtypes=0
+            do i = 1, maxval(type_atom)
+                inbtypes(1) = i 
+                call lattice_map(1, inbtypes,  8 , 1.0_dp, ilat) !Please check documentation on pycac.out formats
+            end do 
+            !Read element and node headers
+            read(11,*) textholder
+            read(11,*) textholder
+
+            !read element information, currently only 8 node elements with 1 basis
+            do ie =1, in_eles
+                    read(11,*)  tag, lat_type, textholder, textholder, esize
+                do inod =1, 8
+                    read(11,*)  textholder, textholder, textholder, re(:,1,inod), ee(1,inod), fe(:,1,inod), ve(:,1,inod)
+                end do
+                call add_element(tag, fcc, esize+1, lat_type, sub_box_num, re)
+                call add_element_data(ele_num, ee, fe, ve)
+                node_num = node_num + 8
+            end do
+        end if
+        call set_max_esize
+        return
+  end subroutine
+
     subroutine read_lmpcac(file, displace, temp_box_bd)
         !This subroutine is used to read .cac files which are used with the lammpsCAC format
         character(len=100), intent(in) :: file
@@ -980,7 +1088,7 @@ module io
         !Open the file
         open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
         
-        !Now initialiaze some important variables if they aren't defined
+        !Now initialize some important variables if they aren't defined
         if (max_basisnum==0) max_basisnum = 10
         if (max_ng_node==0) max_ng_node=8
 
@@ -1109,4 +1217,25 @@ module io
         print *, in_lattice_type
 
     end subroutine set_cac
+
+    subroutine set_types(apos)
+        !This code 
+        integer, intent(in) :: apos
+        integer :: i, j,arglen, arg_pos, ntypes
+
+        character(len=100) :: textholder
+
+        arg_pos = apos + 1
+        call get_command_argument(arg_pos, textholder, arglen)
+        if (arglen==0) stop "Missing numtypes in io"
+        read(textholder,*) ntypes
+
+        do i=1,ntypes
+            arg_pos = arg_pos + 1
+            call get_command_argument(arg_pos, textholder, arglen)
+            call add_atom_type(textholder, j)
+        end do
+
+        return
+    end subroutine set_types
 end module io

src/main.f90

@@ -63,6 +63,8 @@ program main
         
         case('-set_cac')
             call set_cac(i)
+        case('-set_types')
+            call set_types(i)
         end select 
     end do 
     !Determine if a mode is being used and what it is. The first argument has to be the mode
@@ -71,7 +73,8 @@ program main
 
     argument = trim(adjustl(argument))
     if (argument(1:2) == '--') then 
-        call call_mode(end_mode_arg, argument)
+        mode = argument
+        call call_mode(end_mode_arg)
     end if
 
     !Check to make sure a mode was called
@@ -109,10 +112,12 @@ program main
 
     !Check to make sure a file was passed to be written out and then write out 
     ! Before building do a check on the file
-    if (outfilenum == 0) then
+    if (trim(adjustl(mode)) /= "--metric") then    
+        if ((outfilenum == 0)) then
             argument = 'none'
             call get_out_file(argument)
         end if
         call write_out
+    end if
 
 end program main

src/mode_metric.f90

@@ -0,0 +1,247 @@
+module mode_metric
+    !This mode is used to calculate continuum metrics for the j
+
+    use parameters
+    use io
+    use elements
+    use neighbors
+
+    implicit none
+    
+    integer :: nfiles
+    character(len=100) :: metric_type
+    real(kind=dp) :: rc_off
+
+    !Save reference positions
+    integer :: np, npreal, nmet
+    real(kind=dp), allocatable :: r_zero(:,:), r_curr(:,:), met(:,:)
+
+    public 
+    contains
+    subroutine metric(arg_pos)
+       !This is the main calling subroutine for the metric code 
+        integer, intent(out) :: arg_pos
+        character(len=100) :: infile, outfile
+
+        integer :: i, ibasis, inod, np_temp, ppos
+        real(kind=dp), dimension(6) :: temp_box_bd
+        
+        !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(:,:,:,:)
+
+        !Parse the command arguments
+        call parse_command(arg_pos)
+    
+        !Now read the first file 
+        call read_in(1, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd)
+        np = atom_num + max_basisnum*max_ng_node*ele_num
+        allocate(r_zero(3,atom_num+max_basisnum*max_ng_node*ele_num), &
+                 r_curr(3,atom_num+max_basisnum*max_ng_node*ele_num)) 
+        r_zero(:,:) = -huge(1.0_dp)
+
+        !Set up the met variable for the user desired metric
+        select case(trim(adjustl(metric_type)))
+        case('def_grad')
+            allocate(met(9, np))
+        case('microrotation')
+            allocate(met(3,np))
+        end select
+
+        !Now set the reference positions
+        call convert_positions(r_zero, npreal)
+
+        !Now calculate the neighbor list for the reference configuration
+        call calc_neighbor(5.0_dp, r_zero, np)
+
+        !Reset element and box
+        call reset_data
+        call reset_box 
+
+        !Now loop over new files
+        do i = 2, nfiles
+            call read_in(i, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd)
+            call convert_positions(r_curr, np_temp)
+            if (npreal /= np_temp) then 
+                print *, "Error in mode_metric where number of points in ", i, "th file is ", np_temp, " and number of points in" &
+                // "reference file is", npreal
+            end if
+            call calc_metric
+            !Now create the output file num and write out to xyz format
+            ppos = scan(trim(infiles(i)),".", BACK= .true.)
+            if ( ppos > 0 ) then 
+                outfile = infiles(i)(1:ppos)//'xyz'
+            else 
+                outfile = infiles(i)//'.xyz'
+            end if
+            call write_metric_xyz(outfile)
+            call reset_data
+            call reset_box
+        end do
+    end subroutine metric
+
+    subroutine parse_command(arg_pos)
+        !This subroutine parses the arguments for mode metric
+        integer, intent(out) :: arg_pos
+
+        integer :: i, arglen
+        character(len=100) :: textholder
+        logical :: file_exists
+
+        !First read the metric to be used
+        call get_command_argument(2,metric_type,arglen)
+        if (arglen == 0) stop "Incomplete mode metric command, check documentation"
+        select case(trim(adjustl(metric_type))) 
+        case("microrotation")
+            continue
+        case default
+            print *, "Mode metric does not accept metric ", metric_type, ". Please select from: microrotation"
+            stop 3
+        end select
+
+        !Now read the cutoff radius
+        call get_command_argument(3,textholder,arglen)
+        if (arglen == 0) stop "Incomplete mode metric command, check documentation"
+        read(textholder, *) rc_off
+
+        !Now read the number of files to read and allocate the variables
+        call get_command_argument(4, textholder, arglen)
+        if (arglen == 0) stop "Incomplete mode metric command, check documentation"
+        read(textholder, *) nfiles
+
+        !Now read the files to be read 
+        do i = 1, nfiles
+            call get_command_argument(4+i, textholder, arglen)
+            call get_in_file(textholder)
+        end do 
+
+        arg_pos = 5+nfiles
+        return
+    end subroutine parse_command
+
+    subroutine calc_metric
+        !This subroutine calculates the continuum metric that we require
+        integer :: i, j, k, nei, ip, jp
+        real(kind=dp) :: def_grad(3,3), omega(3,3), eta(3,3), rij(3), eta_inv(3,3), ftf(3,3), &
+                         U(3,3), R(3,3), Rskew(3,3), oldrij(3)
+
+        !Loop over all points
+        do ip = 1, np
+            eta(:,:) = 0.0_dp
+            omega(:,:) = 0.0_dp
+            def_grad(:,:) = 0.0_dp
+            do jp = 1, nei_num(ip)
+                !Calculate the neighbor vec in current configuration
+                nei = nei_list(jp, ip)
+                rij = r_curr(:,nei) - r_curr(:,ip)
+                oldrij = r_zero(:,nei) - r_zero(:,ip)
+
+                !Calculate  eta and omega
+                do i = 1,3
+                    do j = 1,3
+                        omega(i,j) = omega(i,j) + rij(i) * oldrij(j)
+                        eta(i,j) = eta(i,j) + oldrij(i) * oldrij(j)
+                    end do
+                end do 
+            end do
+            
+            eta_inv=matinv3(eta)
+            def_grad=matmul(omega,eta_inv)
+
+            select case(trim(adjustl(metric_type)))
+                case('def_grad')
+                    k = 1
+                    do i = 1,3
+                        do j = 1, 3
+                            met(k, ip) = def_grad(i,j)   
+                        end do
+                    end do
+                case('microrotation')
+                    met(:,ip) = 0.0_dp
+                    if(.not.all(def_grad == 0)) then 
+                        !Now calculate microrotation
+                        ftf = matmul(transpose3(def_grad), def_grad)
+                        U = sqrt3(ftf)
+                        if(.not.all(abs(U) < lim_zero)) then 
+                            R = matmul(def_grad, matinv3(U))
+                            Rskew = 0.5_dp * ( R - transpose3(R))
+                            do k =1,3
+                                do j = 1,3
+                                    do i = 1,3
+                                        met(k,ip) = met(k,ip) -0.5*permutation(i,j,k)*Rskew(i,j)
+                                    end do
+                                end do
+                            end do
+                        end if
+                    end if
+            end select
+        end do
+        return 
+    end subroutine
+
+    subroutine convert_positions(rout, npoints)
+        !This subroutine just converts current atom and element arrays to a single point based form
+        real(kind=dp), dimension(3,atom_num+max_ng_node*max_basisnum*ele_num), intent(inout)  :: rout
+        integer, intent(out) :: npoints
+
+        integer :: i, inod, ibasis
+
+        npoints=0
+
+        if(atom_num > 0) then 
+            do i = 1, atom_num
+                rout(:,tag_atom(i))  = r_atom(:,i)
+                npoints = npoints + 1
+            end do
+        end if
+
+        if (ele_num > 0) then 
+            do i = 1, ele_num
+                do inod = 1, ng_node(lat_ele(i))
+                    do ibasis = 1, basisnum(lat_ele(i))
+                        rout(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis) &
+                            = r_node(:,ibasis,inod,i)
+                        npoints = npoints + 1
+                    end do
+                end do
+            end do
+        end if
+
+    end subroutine convert_positions
+
+    subroutine write_metric_xyz(outfile)
+        character(len=100), intent(in) :: outfile
+
+        integer :: i, inod, ibasis
+        real(kind = dp) :: r(3), eng
+        open (unit=11, file=trim(adjustl(outfile)), action='write', position='rewind', status = 'replace')
+        !Write the header
+        write(11,*) npreal
+        
+        select case(metric_type)
+        case('def_grad')
+            write(11,*) "type x y z F11 F12 F13 F21 F22 F23 F31 F32 F33"
+        case('microrotation')
+            write(11,*) "type x y z micro1 micro2 micro3"
+        end select
+
+        if(atom_num > 0) then 
+            do i = 1, atom_num
+                write(11,*) type_atom(i), r_atom(:,i), met(:,tag_atom(i))
+            end do
+        end if
+
+        if (ele_num > 0) then 
+            do i = 1, ele_num
+                do inod = 1, ng_node(lat_ele(i))
+                    do ibasis = 1, basisnum(lat_ele(i))
+                        write(11,*) basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i), &
+                                    met(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis) 
+                    end do
+                end do
+            end do
+        end if
+    end subroutine write_metric_xyz
+
+end module mode_metric

src/neighbors.f90

@@ -0,0 +1,142 @@
+module neighbors
+
+    use parameters
+    use elements
+    use subroutines
+    use functions
+
+    integer, allocatable :: nei_list(:,:), nei_num(:)
+    real(kind=dp), allocatable :: init_vec(:,:,:), output(:), microrotation(:,:)
+    public 
+    contains
+ 
+    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
+            !Check to see if the current point is a filler point and if so just skip it
+            if(r_list(1,i) < -huge(1.0_dp)+1) cycle
+
+            !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
+
+    subroutine calc_neighbor(rc_off, r_list, n)
+        !This function populates the neighbor list in this module
+        
+        real(kind=dp), intent(in) :: rc_off
+        integer, intent(in) :: n
+        real(kind=dp), dimension(3,n) :: r_list
+
+        integer :: i, c(3), ci, cj, ck, num_nei, nei
+        !Variables for cell list code
+        integer, dimension(3) ::cell_num
+        integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:)
+        integer :: which_cell(3,n)
+        
+        !First reallocate the neighbor list codes
+        if (allocated(nei_list)) then
+            deallocate(nei_list,nei_num)
+        end if
+
+        allocate(nei_list(100,n),nei_num(n))
+
+        !Now first pass the position list and and point num to the cell algorithm
+        call build_cell_list(n, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell)
+
+        !Now loop over every point and find it's neighbors
+        pointloop: do i = 1, n
+            
+            !First check to see if the point is a filler point, if so then skip it 
+            if(r_list(1,i) < -Huge(-1.0_dp)+1) cycle
+
+            !c is the position of the cell that the point
+            c = which_cell(:,i)
+
+            !Loop over all neighboring cells
+            do ci = -1, 1, 1
+                do cj = -1, 1, 1
+                    do ck = -1, 1, 1
+                        !First check to make sure that the neighboring cell exists
+                        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)) then 
+
+                                !Now check to see if it is in the cutoff radius, if it is add it to the neighbor list for that
+                                !atom and calculate the initial neighbor vector
+                                if (norm2(r_list(:,nei)-r_list(:,i)) < rc_off) then 
+                                    
+                                    nei_num(i) = nei_num(i) + 1
+                                    nei_list(nei_num(i), i) = nei
+
+                                end if
+                            end if
+                        end do 
+                    end do
+                end do
+            end do
+
+        end do pointloop
+
+        return
+    end subroutine calc_neighbor
+
+end module neighbors

src/opt_delete.f90

@@ -3,6 +3,7 @@ module opt_delete
     use parameters
     use subroutines
     use elements
+    use neighbors
 
     implicit none 
 

src/opt_disl.f90

@@ -561,7 +561,6 @@ module opt_disl
         !Now reset the list for the scanning algorithm
 
         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

src/opt_group.f90

@@ -87,7 +87,7 @@ module opt_group
             continue
         case default
             print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", &
-                     "Please select from atoms, nodes, or both."
+                     "Please select from atoms, elements, or both."
         end select
 
         arg_pos = arg_pos + 1

src/parameters.f90

@@ -14,4 +14,6 @@ module parameters
     !Numeric constants
     real(kind=dp), parameter :: pi = 3.14159265358979323846_dp
     
+    !Mode type that is being called 
+    character(len=100) :: mode
 end module parameters

src/subroutines.f90

@@ -198,66 +198,6 @@ 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
-
-
     subroutine check_right_ortho(ori, isortho, isrighthanded)
         !This subroutine checks whether provided orientations in the form:
         ! | x1  x2  x3 |
@@ -301,4 +241,5 @@ module subroutines
         return
     end subroutine check_right_ortho
 
+
 end module subroutines