From 0e41c7acc9a1d1dc4b3febb335a29282ff8f8907 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 14 Jun 2021 10:23:26 -0400 Subject: [PATCH] Lots of updates to cacmb --- src/elements.f90 | 25 ++++++- src/io.f90 | 9 ++- src/mode_create.f90 | 4 +- src/mode_merge.f90 | 176 +++++++++++++++++++++++++++++++++++++++++++- src/neighbors.f90 | 5 +- src/opt_deform.f90 | 32 +++++++- src/opt_group.f90 | 78 ++++++++++++++++++-- 7 files changed, 310 insertions(+), 19 deletions(-) diff --git a/src/elements.f90 b/src/elements.f90 index 73807ab..3e73d37 100644 --- a/src/elements.f90 +++ b/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 diff --git a/src/io.f90 b/src/io.f90 index d5dbab4..a7947a5 100644 --- a/src/io.f90 +++ b/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) diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 7994f55..29e4d05 100644 --- a/src/mode_create.f90 +++ b/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)),& - 1:(bd_ele_lat(6)-bd_ele_lat(5)))= lat_points(bd_ele_lat(1):bd_ele_lat(2), & + 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))+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 diff --git a/src/mode_merge.f90 b/src/mode_merge.f90 index 2052022..8a812b1 100644 --- a/src/mode_merge.f90 +++ b/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) - logical :: shift_flag + real(kind=dp) :: shift_vec(3), replace_vec(3) + 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 diff --git a/src/neighbors.f90 b/src/neighbors.f90 index f9d4d48..d3683e6 100644 --- a/src/neighbors.f90 +++ b/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 diff --git a/src/opt_deform.f90 b/src/opt_deform.f90 index c9be654..bba35cf 100644 --- a/src/opt_deform.f90 +++ b/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 diff --git a/src/opt_group.f90 b/src/opt_group.f90 index ad63a15..78f8029 100644 --- a/src/opt_group.f90 +++ b/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