Lots of updates to cacmb

src/elements.f90

@@ -393,7 +393,7 @@ module elements
         !This subroutine returns the interpolated atoms from the elements. 
 
         !Arguments
-        character(len=100), intent(in) :: type !The type of element that it is
+        character(len=*), intent(in) :: type !The type of element that it is
         integer, intent(in) :: esize !The number of atoms per side
         integer, intent(in) :: lat_type !The integer lattice type of the element
         real(kind=dp), dimension(3,max_basisnum, max_ng_node), intent(in) :: r_in !Nodal positions
@@ -863,4 +863,27 @@ do i = 1, atom_num
         node_num = 0
     end subroutine reset_data
 
+    function ele_in_bounds(bd, etype, esize, lat_type, r_in )
+        real(kind=dp), intent(in) :: bd(6), r_in(3, max_basisnum, max_ng_node)
+        integer, intent(in) ::  lat_type, esize
+        character(len=*), intent(in) :: etype
+
+        integer :: i
+        real(kind=dp) :: rinterp(3, max_basisnum*max_esize**3) 
+        integer :: tinterp(max_basisnum*max_esize**3) 
+        logical :: ele_in_bounds
+
+        ele_in_bounds=.false.
+        call interpolate_atoms(etype, esize, lat_type, r_in, tinterp, rinterp) 
+
+        do i=1, esize**3
+            if(in_block_bd(rinterp(:,i), bd))then 
+                ele_in_bounds = .true.
+                exit
+            end if
+        end do
+
+        return
+    end function ele_in_bounds
+
 end module elements

src/io.f90

@@ -865,7 +865,7 @@ module io
                    atom_type_map(100), etype_map(100),  lat_type, new_lattice_map(100), &
                    atom_type, stat, bnum, nnum, esize, btypes(10), tmp, ip, jp
         real(kind=dp) :: newdisplace(3), r_in(3,10,8), r_in_atom(3), atomic_masses(10)
-        character(len=100) :: textholder, etype
+        character(len=1000) :: textholder, etype
         character(len=2) :: atomic_element
         !First open the file
         open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') 
@@ -897,9 +897,11 @@ module io
         read(textholder, *) (atomic_masses(i), i=1, j)
 
         !Read define atom_types by mass
+        print *, j
         do i = 1, j
             call atommassspecies(atomic_masses(i), atomic_element)
             call add_atom_type(atomic_element, atom_type_map(i))
+            print *, i, atom_type_map(i)
         end do 
 
         !Read in the boundary
@@ -1105,7 +1107,7 @@ module io
         real(kind=dp), dimension(3), intent(in) :: displace
         real(kind = dp), dimension(6), intent(out) :: temp_box_bd
         
-        character(len=100) :: textholder, element_type
+        character(len=1000) :: textholder, element_type
         character(len=2) :: atom_species
         integer :: i, j, k, atom_in, type_in, type_map(10), in_basis, node_types(10,8), &
                    lat_type, id 
@@ -1175,7 +1177,8 @@ module io
         !Start the reading loop
         do i = 1, atom_in
             read(11,*) id, type_in, r_in(:)
-            call add_atom(id, type_in, sub_box_num, r_in) 
+            r_in(:) = r_in + newdisplace
+            call add_atom(id, type_map(type_in), sub_box_num, r_in) 
         end do
 
         close(11)

src/mode_create.f90

@@ -515,8 +515,8 @@ module mode_create
                                 end if
                             end do
 
-                            lat_points_ele(1:(bd_ele_lat(2)-bd_ele_lat(1)),1:(bd_ele_lat(4)-bd_ele_lat(3)),&
+                            lat_points_ele(1:(bd_ele_lat(2)-bd_ele_lat(1)+1),1:(bd_ele_lat(4)-bd_ele_lat(3)+1),&
-                            1:(bd_ele_lat(6)-bd_ele_lat(5)))= lat_points(bd_ele_lat(1):bd_ele_lat(2), &
+                            1:(bd_ele_lat(6)-bd_ele_lat(5))+1)= lat_points(bd_ele_lat(1):bd_ele_lat(2), &
                                                                          bd_ele_lat(3):bd_ele_lat(4), &
                                                                          bd_ele_lat(5):bd_ele_lat(6)) 
                             !Now start looping through elements and try to fit as many as you can

src/mode_merge.f90

@@ -6,18 +6,23 @@ module mode_merge
     use io
     use subroutines
     use elements
+    use neighbors
+
+    implicit none
 
     character(len=4) :: dim
     integer :: in_num, new_starts(2)
-    real(kind=dp) :: shift_vec(3)
+    real(kind=dp) :: shift_vec(3), replace_vec(3) 
-    logical :: shift_flag
+    character(len=100) :: replace_str(3)
+    logical :: shift_flag, replace_flag
+    real(kind=dp), private, save :: rc_off 
 
     public 
     contains
     subroutine merge(arg_pos)
         
         integer, intent(out) :: arg_pos 
-        integer :: i
+        integer :: i, j 
         real(kind=dp) :: displace(3), temp_box_bd(6)
 
         print *, '-----------------------Mode Merge---------------------------'
@@ -55,9 +60,20 @@ module mode_merge
                 call read_in(i, displace, temp_box_bd)
             end if 
 
+
             if(shift_flag) call shift(new_starts, i)
+
+            if(replace_flag.and.(i>1)) then 
+                !Parse the replace vector
+                do j = 1, 3
+                    call parse_pos(j, replace_str(j), replace_vec(j))
+                end do
+                call replace(new_starts, temp_box_bd)
+            end if
+
         end do
 
+
         !Now reset tags
         do i = 1, atom_num
             tag_atom(i) = i
@@ -122,6 +138,17 @@ module mode_merge
                     if (arglen==0) stop "Missing vector component for shift command"
                     read(textholder, *) shift_vec(i)
                 end do
+            case('replace')
+                replace_flag = .true.
+                do i = 1,3
+                    arg_pos = arg_pos + 1
+                    call get_command_argument(arg_pos, replace_str(i), arglen)
+                    if (arglen==0) stop "Missing vector component for shift command"
+                end do
+                arg_pos = arg_pos+1
+                call get_command_argument(arg_pos, textholder, arglen)
+                read(textholder,*) rc_off
+
             case default
                 !If it isn't an available option to mode merge then we just exit
                 exit 
@@ -176,4 +203,147 @@ module mode_merge
         end if
 
     end subroutine shift
+
+    subroutine replace(array_start, rbox_bd)
+        integer, intent(in) :: array_start(2)
+        real(kind = dp), intent(in) :: rbox_bd(6)
+
+        integer ::  ibasis, inod, del_num, del_index(atom_num), nump_ele, interp_start
+        integer ::  j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num, m, n, o, esize, &
+                   ele(3,8), new_ele_num, vlat(3), added_points
+        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum,max_ng_node), ravg(3), ratom(3,max_basisnum)
+        logical :: in_bd, lat_points(max_esize, max_esize, max_esize)
+        real(kind=dp) :: del_bd(6)
+        integer :: i, c(3), ci, cj, ck, num_nei, nei, delete_num
+
+        !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(:,:,:,:)
+
+        !First apply the replace vec to all new nodes and elements
+        do i = array_start(1), atom_num
+            r_atom(:,i) = r_atom(:, i) + replace_vec
+        end do  
+
+        do i = array_start(2), ele_num
+            do inod = 1, ng_node(lat_ele(i))
+                do ibasis=1, basisnum(lat_ele(i))
+                    r_node(:, ibasis,inod, i) = r_node(:, ibasis,inod, i) + replace_vec 
+                end do
+            end do
+
+        end do
+
+        !Calculate new boundary
+        do i = 1, 6
+            del_bd(i) = rbox_bd(i) + replace_vec((i-1)/2 + 1)
+        end do
+        
+        del_num = 0
+        del_index=0
+        interp_start = atom_num +1
+        !Now loop over all old elements, 
+        do ie = 1, array_start(2)-1
+            !If any element points are within the boundary then we run the refine code
+            if(ele_in_bounds(del_bd, type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie))) then 
+                added_points=0
+                del_num = del_num + 1
+                del_index(del_num) = ie
+
+                !Find all possible elements that we can make while making sure they aren't in the group
+                lat_points(1:size_ele(ie),1:size_ele(ie),1:size_ele(ie)) = .true.
+
+                !Now add the leftover lattice points as atoms, only if they aren't within the new boundaries
+                do o = 1, size_ele(ie)
+                    do n = 1, size_ele(ie)
+                        do m = 1, size_ele(ie)
+                            if(lat_points(m,n,o)) then 
+                                call get_interp_pos(m,n,o, ie, ratom(:,:))
+                                do ibasis = 1, basisnum(lat_ele(ie))
+                                    call apply_periodic(ratom(:,ibasis))
+                                    added_points=added_points + 1
+                                    call add_atom(0, basis_type(ibasis,lat_ele(ie)), sbox_ele(ie), ratom(:,ibasis))
+                                end do
+                            end if
+                        end do
+                    end do
+                end do
+
+                if (added_points /= (size_ele(ie)**3)) then 
+                    print *, "Element ", ie, " is refined incorrectly in refinefill"
+                end if
+            end if
+        end do
+
+        !Once all atoms are added we delete all of the elements
+        call delete_elements(del_num, del_index)
+
+        !Now delete overlapping atoms
+        allocate(which_cell(3,atom_num))
+
+        !First pass the atom list and atom num to the algorithm which builds the cell list
+        print *, rc_off
+        call build_cell_list(atom_num, r_atom, 4*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
+        del_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 
+                                        
+                                        del_num = del_num + 1
+
+                                        !Make sure to delete the older value
+                                        if( (i < array_start(1)).or.(i > interp_start)) then 
+                                            del_index(del_num) = i
+                                            which_cell(:,i) = 0
+                                            cycle atom_loop
+                                        else
+                                            del_index(del_num) = nei
+                                            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 *, "Replace command deletes ", del_num, " atoms"
+        !Now delete all the atoms
+        call delete_atoms(del_num, del_index(1:del_num))
+
+
+        return
+    end subroutine replace
+
 end module mode_merge

src/neighbors.f90

@@ -41,12 +41,14 @@ module neighbors
             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 
+        print *, box_len, cell_num
 
         !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
+        num_in_cell = 0
         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
@@ -60,9 +62,10 @@ module neighbors
             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(1:cell_lim, :, :, :) = cell_list(1:cell_lim, :, :, :)
                 resize_cell_list(cell_lim+1:, :, :, :) = 0
                 call move_alloc(resize_cell_list, cell_list)
+                cell_lim = cell_lim + 10
             end if
 
             cell_list(num_in_cell(c(1),c(2),c(3)),c(1),c(2),c(3)) = i

src/opt_deform.f90

@@ -10,6 +10,7 @@ module opt_deform
 
     real(kind=dp), save :: applied_strain
     integer, save :: sdim
+    logical :: percent 
 
     public 
     contains
@@ -21,9 +22,11 @@ module opt_deform
 
         character(len=1) :: dims(3)
         integer :: i, j, k
-        real(kind=dp) :: frac_atom(atom_num), frac_node(max_basisnum, max_ng_node, ele_num)
+        real(kind=dp) :: frac_atom(atom_num), frac_node(max_basisnum, max_ng_node, ele_num), blen
+
 
         !initialize some variables
+        percent=.false.
         dims(1) = 'x'
         dims(2) = 'y'
         dims(3) = 'z'
@@ -44,7 +47,13 @@ module opt_deform
         end do
 
         print *, "Original box bounds in ", dims(sdim), " are ", box_bd(2*sdim-1:2*sdim)
-        box_bd(2*sdim) = box_bd(2*sdim) + applied_strain
+        if(percent) then 
+            blen = box_bd(2*sdim) - box_bd(2*sdim-1) 
+            box_bd(2*sdim) = box_bd(2*sdim-1) + blen*applied_strain
+        else
+            box_bd(2*sdim) = box_bd(2*sdim) + applied_strain
+        end if
+
         print *, "New box bounds are ", box_bd(2*sdim-1:2*sdim)
 
         !Now reassign the positions
@@ -90,6 +99,25 @@ module opt_deform
         if (arg_len == 0) stop "Missing strain in deform command"
         read(textholder, *) applied_strain
 
+        !Now parse the additional options which may be present
+        do while(.true.)
+            if(arg_pos > command_argument_count()) exit
+            !Pull out the next argument which should either be a keyword or an option
+            arg_pos=arg_pos+1
+            call get_command_argument(arg_pos, textholder)
+            textholder=adjustl(textholder)
+
+            !Choose what to based on what the option string is 
+            select case(trim(textholder))
+            case('percent')
+                percent=.true. 
+            case default
+                arg_pos=arg_pos-1
+                !if it isn't an available option to opt_group then we just exit
+                exit 
+            end select
+        end do
+
         arg_pos = arg_pos + 1
 
     end subroutine parse_deform

src/opt_group.f90

@@ -6,15 +6,18 @@ module opt_group
     use elements
     use subroutines
     use box 
+    use sorts
+
     implicit none
 
     integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize, insert_type, &
-               insert_site
+               insert_site, num_species
     character(len=15) :: type, gshape!Type indicates what element type is selected and shape is the group shape
+    character(len=2) :: species_type(10)
     real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness, insert_conc, &
-    insert_lattice
+    insert_lattice, s_fractions(10)
 
-    logical :: displace, delete, max_remesh, refine, group_nodes, flip, efill, refinefill
+    logical :: displace, delete, max_remesh, refine, group_nodes, flip, efill, refinefill, alloy
 
     integer, allocatable :: element_index(:), atom_index(:) 
 
@@ -41,6 +44,8 @@ module opt_group
         max_remesh=.false.
         refine = .false.
         flip = .false.
+        group_nodes=.false.
+        num_species = 0
 
         if(allocated(element_index)) deallocate(element_index)
         if(allocated(atom_index)) deallocate(atom_index)
@@ -83,6 +88,12 @@ module opt_group
             call insert_group
         end if
 
+        if(num_species > 0) then 
+            call get_group
+            call alloy_group
+        end if
+
+
     end subroutine group
 
     subroutine parse_group(arg_pos)
@@ -434,6 +445,23 @@ module opt_group
                 if(arglen ==0) stop "Missing notsize size"
                 read(textholder, *) notsize
                 print *, "Ignoring elements with size ", notsize
+            case('alloy')
+                alloy=.true.
+                arg_pos=arg_pos+1
+                call get_command_argument(arg_pos, textholder, arglen)
+                print *, textholder
+                if(arglen == 0) stop "Missing alloy num"
+                read(textholder, *) num_species
+                do i = 1, num_species
+                    arg_pos=arg_pos + 1
+                    call get_command_argument(arg_pos, textholder, arglen)
+                    print *, textholder
+                    read(textholder, *) species_type(i)
+                    arg_pos=arg_pos + 1
+                    call get_command_argument(arg_pos, textholder, arglen)
+                    print *, textholder
+                    read(textholder, *) s_fractions(i)
+                end do
             case('insert')
                 arg_pos=arg_pos+1
                 call get_command_argument(arg_pos, textholder, arglen)
@@ -454,14 +482,14 @@ module opt_group
                 end select
                 arg_pos=arg_pos+1
                 call get_command_argument(arg_pos, textholder, arglen)
-                if(arglen ==0) stop "Missing lattice_type in insert command"
+                if(arglen ==0) stop "missing lattice_type in insert command"
                 read(textholder, *) insert_lattice
                 arg_pos=arg_pos+1
                 call get_command_argument(arg_pos, textholder, arglen)
-                if(arglen ==0) stop "Missing concentration in insert command"
+                if(arglen ==0) stop "missing concentration in insert command"
                 read(textholder, *) insert_conc
             case default
-                !If it isn't an available option to opt_group then we just exit
+                !if it isn't an available option to opt_group then we just exit
                 exit 
             end select
         end do
@@ -1107,6 +1135,43 @@ module opt_group
         end do addloop
     end subroutine insert_group
 
+    subroutine alloy_group
+        !This subroutine randomizes the atom types to reach desired concentrations, this only operates on atoms
+        integer :: i, j, ia, type_map(10), added_types(num_species)
+        real(kind=dp) ::  rand
+
+
+        print *, "Alloying group with desired fractions", s_fractions(1:num_species)
+        !First we sort the fractions
+        call quicksort(s_fractions(1:num_species))
+
+
+        !Now get the atom type maps for all the atoms and make the fractions a running sum
+        do i = 1, num_species
+            call add_atom_type(species_type(i), type_map(i))
+            if(i > 1) s_fractions(i) = s_fractions(i) + s_fractions(i-1)
+        end do
+        print *, s_fractions
+        !Now randomize the atom types
+        added_types = 0
+        print *, type_map(j)
+        do i = 1, group_atom_num
+            ia = atom_index(i)
+            call random_number(rand)
+            do j = 1, num_species
+                if(rand < s_fractions(j)) then 
+                    type_atom(ia) = type_map(j)
+                    added_types(j) = added_types(j) +1 
+                    exit
+                end if
+            end do
+        end do
+
+        print *, "Converted ", added_types(1:num_species), " atoms"
+
+        return
+    end subroutine alloy_group
+
     function in_group(r)
         !This subroutine determines if a point is within the group boundaries
         real(kind=dp), intent(in) :: r(3)
@@ -1218,7 +1283,6 @@ module opt_group
             end do
         end if
 
-
     end function in_group_ele
 
 end module opt_group