module opt_group !This module contains all code associated with dislocations use parameters use elements use subroutines use box implicit none integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness logical :: displace, delete, max_remesh, refine, group_nodes, flip, efill, refinefill integer, allocatable :: element_index(:), atom_index(:) public contains subroutine group(arg_pos) !Main calling function for the group option integer, intent(inout) :: arg_pos print *, '------------------------------------------------------------' print *, 'Option Group' print *, '------------------------------------------------------------' group_ele_num = 0 group_atom_num = 0 remesh_size=0 random_num=0 group_type=0 notsize=0 displace=.false. delete=.false. max_remesh=.false. refine = .false. flip = .false. if(allocated(element_index)) deallocate(element_index) if(allocated(atom_index)) deallocate(atom_index) call parse_group(arg_pos) !Now call the transformation functions for the group if(refine) then call get_group call refine_group end if if(refinefill) then call get_group call refinefill_group end if if(displace)then call get_group call displace_group end if if(delete)then call get_group call delete_group end if if(remesh_size > 0)then call get_group call remesh_group end if if(group_type > 0) then call get_group call change_group_type end if end subroutine group subroutine parse_group(arg_pos) !Parse the group command integer, intent(inout) :: arg_pos integer :: i, j, arglen, in_num character(len=100) :: textholder, type_spec real(kind=dp) H !Parse type and shape command arg_pos = arg_pos + 1 call get_command_argument(arg_pos, type, arglen) if (arglen==0) STOP "Missing select_type in group command" select case(trim(adjustl(type))) case('atoms', 'elements', 'both') continue case default print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", & "Please select from atoms, elements, or both." end select arg_pos = arg_pos + 1 call get_command_argument(arg_pos, shape, arglen) if (arglen==0) STOP "Missing group_shape in group command" !Now parse the arguments required by the user selected shape select case(trim(adjustl(shape))) case('block') do i= 1, 6 arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing block boundary in group command" call parse_pos(int((i+1)/2), textholder, block_bd(i)) end do case('wedge') arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing normal dim in group wedge command" read(textholder,*) dim1 arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing normal dim in group wedge command" read(textholder,*) dim2 do i = 1, 3 arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing centroid in group wedge command" call parse_pos(i, textholder, centroid(i)) end do arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing base width in group wedge command" read(textholder,*) bwidth !Calculate the vertex positions vertices(:,1) = centroid vertices(dim2,1) = 0.0_dp do i = 1, 3 if (i == dim1) then if (bwidth > 0) then vertices(i,2) = box_bd(2*i) vertices(i,3) = box_bd(2*i) else if (bwidth < 0) then vertices(i,2) = box_bd(2*i-1) vertices(i,3) = box_bd(2*i-1) else print *, "bwidth cannot be 0 in wedge shaped group" stop 3 end if else if (i == dim2) then vertices(i,2) = 0.0_dp vertices(i,3) = 0.0_dp else vertices(i,2) = centroid(i) + bwidth vertices(i,3) = centroid(i) - bwidth end if end do case('notch') arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing normal dim in group notch command" read(textholder,*) dim1 arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing normal dim in group notch command" read(textholder,*) dim2 do i = 1, 3 arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing centroid in group notch command" call parse_pos(i, textholder, centroid(i)) end do arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing base width in group notch command" read(textholder,*) bwidth arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing tip radius in group notch command" read(textholder,*) radius !Calculate the vertex positions vertices(:,1) = centroid vertices(dim2,1) = 0.0_dp do i = 1, 3 if (i == dim1) then if (bwidth > 0) then vertices(i,2) = box_bd(2*i) vertices(i,3) = box_bd(2*i) H= (box_bd(2*i) - centroid(i)) !Calculate the height of the wedge else if (bwidth < 0) then vertices(i,2) = box_bd(2*i-1) vertices(i,3) = box_bd(2*i-1) H = (centroid(i) - box_bd(2*i-1)) else print *, "bwidth cannot be 0 in wedge shaped group" stop 3 end if else if (i == dim2) then vertices(i,2) = 0.0_dp vertices(i,3) = 0.0_dp else vertices(i,2) = centroid(i) + 0.5_dp*bwidth vertices(i,3) = centroid(i) - 0.5_dp*bwidth end if end do !Now update the centroid so that the desire tip diameter matches the wedge with if (bwidth > 0) then centroid(dim1) = centroid(dim1) + 2*radius*(H)/bwidth end if !Read the ID type shape for group case('id') arg_pos = arg_pos + 1 !For this type we have to call different options if type is atoms, elements, or both. Both is the most complex. select case(trim(adjustl(type))) case('atoms') !Read number of ids call get_command_argument(arg_pos, textholder, arglen) if(arglen == 0) then print *, "Missing number of input atoms for group id" end if read(textholder,*) in_num !allocate arrays allocate(atom_index(in_num)) !Read ids do i = 1, in_num arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if(arglen == 0) then print *, "Missing atom id in group atom id" stop 3 end if read(textholder,*) atom_index(i) end do group_atom_num = in_num case('elements') !Read number of ids call get_command_argument(arg_pos, textholder, arglen) if(arglen == 0) then print *, "Missing number of input elements for group id" end if read(textholder, *) in_num !allocate arrays allocate(element_index(in_num)) !Read ids do i = 1, in_num arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if(arglen == 0) then print *, "Missing element id in group element id" stop 3 end if read(textholder,*) element_index(i) end do group_ele_num = in_num case('both') !We repeat this code twice, once for the atoms and once for the elements allocate(element_index(1024),atom_index(1024)) do i = 1, 2 !Read the first id type (either atoms or elements) call get_command_argument(arg_pos, type_spec, arglen) if (arglen == 0) then print *, "Missing type specifier in group id command" stop 3 end if !Now read the first in_num arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if(arglen == 0) then print *, "Missing input number in group_id" stop 3 end if read(textholder, *) in_num select case(trim(adjustl(type_spec))) case('atoms','atom') if (group_atom_num > 0) then print *, "Atoms specifier used more than once in group id command with type both, either use type ", & "atoms or include elements identifier" stop 3 do j = 1, in_num arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen == 0) then print *, "Missing atom id in group atom id" stop 3 end if read(textholder, *) atom_index(j) end do group_atom_num = in_num end if case('elements','element') if (group_ele_num > 0) then print *, "Elements specifier used more than once in group id command with type both, either use type ",& "elements or include atoms identifier" stop 3 do j = 1, in_num arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen == 0) then print *, "Missing element id in group element id" stop 3 end if read(textholder, *) element_index(j) end do group_ele_num = in_num end if end select if(i ==1) arg_pos = arg_pos + 1 end do end select case('sphere') !First extract the sphere centroid do i = 1, 3 arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing sphere centroid in group command" call parse_pos(i, textholder, centroid(i)) end do !Now get the tip radius arg_pos=arg_pos+1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing sphere radius in group command" read(textholder, *) radius case('shell') !First extract the shell centroid do i = 1, 3 arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing shell centroid in group command" call parse_pos(i, textholder, centroid(i)) end do !Now get the radius arg_pos=arg_pos+1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing shell radius in group command" read(textholder, *) radius !Now get the shell thickness arg_pos=arg_pos+1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) STOP "Missing shell thickness in group command" read(textholder, *) shell_thickness case('all') !Do nothing if the shape is all continue case default print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", & "for accepted group shapes." end select !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('displace') displace = .true. do i = 1,3 arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) stop "Missing vector component for shift command" read(textholder, *) disp_vec(i) end do case('refine') refine=.true. case('refinefill') refinefill=.true. case('remesh') arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) stop "Missing remesh element size in group command" read(textholder, *) remesh_size case('max') max_remesh =.true. case('delete') delete=.true. case('nodes') group_nodes=.true. case('random') arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) stop "Missing number of random atoms in group command" read(textholder, *) random_num case('flip') flip=.true. case('efill') efill=.true. case('type') arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if (arglen==0) stop "Missing atom type for group" call add_atom_type(textholder, group_type) case('notsize') arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) if(arglen ==0) stop "Missing notsize size" read(textholder, *) notsize print *, "Ignoring elements with size ", notsize case default !If it isn't an available option to opt_disl then we just exit exit end select end do end subroutine parse_group subroutine get_group !This subroutine finds all elements and/or atoms within the group boundaries !specified by the user. integer :: i, j, inod, ibasis, temp, node_in_out(max_ng_node) integer, allocatable :: resize_array(:) real(kind=dp) :: r_center(3), rand select case(trim(adjustl(shape))) case('block') print *, "Group has block shape with boundaries: ", block_bd case ('wedge') print *, "Group has wedge shape with dim1", dim1, "and dim2", dim2, "and vertices ", vertices case ('notch') print *, "Group has notch shape with dim1", dim1, "and dim2", dim2, " tip radius ", radius, "and vertices ", & vertices case('id') print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms." return case('sphere') print *, "Group has sphere shape with centroid ", centroid(:), " and radius ", radius case('all') print *, "Group has all of type ", type end select !Reset group if needed if(allocated(element_index)) deallocate(element_index,atom_index) group_ele_num = 0 group_atom_num = 0 allocate(element_index(1024), atom_index(1024)) !Check the type to see whether we need to find the elements within the group select case(trim(adjustl(type))) case('elements', 'both') do i = 1, ele_num if(in_group_ele(size_ele(i), lat_ele(i), r_node(:,:,:,i))) then group_ele_num = group_ele_num + 1 if(group_ele_num > size(element_index)) then allocate(resize_array(size(element_index) + 1024)) resize_array(1:group_ele_num-1) = element_index resize_array(group_ele_num:) = 0 call move_alloc(resize_array, element_index) end if element_index(group_ele_num) = i end if end do if(random_num > 0) then !If we have the random option enabled then we select random_num number of elements from the group and overwrite !the group with those elements do i = 1, random_num call random_number(rand) j = i + floor((group_ele_num+1-i)*rand) temp = element_index(j) element_index(j) = element_index(i) element_index(i) = temp end do group_ele_num = random_num end if end select !Check the type to see if we need to find the atoms within the group select case(trim(adjustl(type))) case('atoms','both') do i = 1, atom_num if(in_group(r_atom(:,i)).neqv.flip) then group_atom_num = group_atom_num + 1 if (group_atom_num > size(atom_index)) then allocate(resize_array(size(atom_index) + 1024)) resize_array(1:group_atom_num -1) = atom_index resize_array(group_atom_num:) = 0 call move_alloc(resize_array, atom_index) end if atom_index(group_atom_num) = i end if end do if(random_num > 0) then !If we have the random option enabled then we select random_num number of atom from the group and overwrite !the group with those atom do i = 1, random_num call random_number(rand) j = i + floor((group_atom_num+1-i)*rand) temp = atom_index(j) atom_index(j) = atom_index(i) atom_index(i) = temp end do group_atom_num = random_num end if end select j = 0 do i = 1, group_atom_num if (atom_index(i) == 23318348) then j = j + 1 end if end do if (j > 1) stop "Code broken" print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms." end subroutine get_group subroutine displace_group !This subroutine applies a displacement to elements/atoms in the groups integer :: i, inod, ibasis print *, "Elements/atoms in group displaced by ", disp_vec !Displace atoms do i = 1, group_atom_num r_atom(:,atom_index(i)) = r_atom(:,atom_index(i)) + disp_vec end do !Displace elements do i = 1, group_ele_num do inod = 1, ng_node(lat_ele(element_index(i))) do ibasis = 1, basisnum(lat_ele(element_index(i))) r_node(:,ibasis,inod,element_index(i)) = r_node(:,ibasis,inod,element_index(i)) + disp_vec end do end do end do !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 (disp_vec(i) < -lim_zero) then box_bd(2*i-1) = box_bd(2*i-1) - disp_vec(i) else if (disp_vec(i) > lim_zero) then box_bd(2*i) = box_bd(2*i) + disp_vec(i) end if end do end if end subroutine displace_group subroutine refine_group !This command is used to refine the group to full atomistics integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3) !Refining to atoms if(group_ele_num > 0) then orig_atom_num = atom_num !Estimate number of atoms we are adding, this doesn't have to be exact add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3 call grow_ele_arrays(0,add_atom_num) do i = 1, group_ele_num ie = element_index(i) !Get the interpolated atom positions call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) !Loop over all interpolated atoms and add them to the system, we apply periodic boundaries !here as well to make sure they are in the box do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3 call apply_periodic(r_interp(:,j)) call add_atom(0,type_interp(j), sbox_ele(ie), r_interp(:,j)) end do end do !Once all atoms are added we delete all of the elements call delete_elements(group_ele_num, element_index) print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms." end if end subroutine refine_group subroutine refinefill_group !This command is used to refine the group to full atomistics integer :: i, 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, ibasis, inod, vlat(3), nump_ele, 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 :: lat_points(max_esize, max_esize, max_esize) !Refining to atoms if(group_ele_num > 0) then orig_atom_num = atom_num new_ele_num = 0 !Estimate number of atoms we are adding, this doesn't have to be exact add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3 call grow_ele_arrays(0,add_atom_num) do i = 1, group_ele_num ie = element_index(i) !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 calculate the number of elemets which are available for remeshing nump_ele = size_ele(ie)**3 do o =1, size_ele(ie) do n = 1, size_ele(ie) do m =1, size_ele(ie) call get_interp_pos(m,n,o,i,rfill(:,:,1)) ravg(:) = 0 do ibasis = 1, basisnum(lat_ele(ie)) ravg = ravg + rfill(:,ibasis, 1)/basisnum(lat_ele(ie)) end do if( in_group(ravg)) then nump_ele = nump_ele - 1 end if end do end do end do !Now start the remeshing loop for the element esize = size_ele(ie) - 2 added_points=0 do while(esize > min_efillsize) if(nump_ele < min_efillsize**3) then exit else if (nump_ele < esize**3) then esize = esize - 2 else ele = cubic_cell*(esize-1) do o = 1, size_ele(ie) - esize do n = 1, size_ele(ie) - esize latloop:do m = 1, size_ele(ie) - esize do inod = 1, ng_node(lat_ele(ie)) vlat = ele(:,inod) + (/ m, n, o /) if (.not.lat_points(vlat(1), vlat(2),vlat(3))) cycle latloop call get_interp_pos(vlat(1), vlat(2), vlat(3), ie, rfill(:,:,inod)) end do !Check to make sure all lattice points exist for the current element if(any(.not.lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1))) cycle latloop if (.not.in_group_ele(esize, lat_ele(ie), rfill)) then nump_ele=nump_ele - esize**3 lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) = .false. call add_element(0,type_ele(ie), esize, lat_ele(ie), sbox_ele(ie), rfill) new_ele_num = new_ele_num + 1 added_points = added_points + esize**3 end if end do latloop end do end do esize=esize-2 end if end do !Now add the leftover lattice points as atoms 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)) call add_atom(0, basis_type(ibasis,lat_ele(ie)), sbox_ele(ie), ratom(:,ibasis)) added_points=added_points + 1 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 do !Once all atoms are added we delete all of the elements call delete_elements(group_ele_num, element_index) print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms and ", new_ele_num, & " elements." end if end subroutine refinefill_group subroutine remesh_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), new_ele, new_atom, & 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(:,:,:), new_sbox, new_lat character(len=100) :: remesh_ele_type !Right now we just hardcode only remeshing to elements remesh_ele_type = 'fcc' ! Determine which sub_boxes and lattices types are within in the 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))) 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 if do j = 1, group_lat_num 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 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))) 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_ele(element_index(i)) end if do j = 1, group_lat_num 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 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)) call matrix_inverse(orient,3,ori_inv) do ilat = 1, group_lat_num !First calculate max position in lattice space to be able to allocate lat_points array, also sum the total number !of degrees of freedom which are added dof = 0 do j = 1, 3 group_bd(2*j) = -huge(1.0_dp) group_bd(2*j-1) = huge(1.0_dp) end do do i = 1, group_atom_num if ((type_atom(atom_index(i)) == basis_type(1,ilat)).and.(sbox_atom(atom_index(i)) == is)) then do j =1 ,3 if (r_atom(j,atom_index(i)) > group_bd(2*j)) group_bd(2*j) = r_atom(j,atom_index(i)) if (r_atom(j,atom_index(i)) < group_bd(2*j-1)) group_bd(2*j-1) = r_atom(j,atom_index(i)) end do dof = dof + 1 end if end do do i = 1, group_ele_num if ((lat_ele(element_index(i)) == ilat).and.(sbox_ele(element_index(i)) == is)) then do inod =1, ng_node(ilat) do ibasis = 1, basisnum(ilat) do j = 1, 3 r =r_node(j,ibasis,inod,element_index(i)) if (r(j) > group_bd(2*j)) group_bd(2*j) = r(j) if (r(j) < group_bd(2*j-1)) group_bd(2*j-1) = r(j) end do end do end do dof = dof + size_ele(element_index(i))**3 end if end do !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), & group_bd(2),group_bd(4),group_bd(5), & group_bd(1),group_bd(4),group_bd(5), & group_bd(1),group_bd(3),group_bd(6), & group_bd(2),group_bd(3),group_bd(6), & group_bd(2),group_bd(4),group_bd(6), & group_bd(1),group_bd(4),group_bd(6) /), [3,8]) group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/lapa(ilat))) do i = 1, 3 bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:))) bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:))) end do allocate(lat_points(bd_in_lat(2)-bd_in_lat(1)+10, bd_in_lat(4)-bd_in_lat(3)+10, bd_in_lat(6)-bd_in_lat(5)+10)) lat_points(:,:,:) = .false. dof = 0 !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 r_lat(j) = nint(r(j)) end do !Do a check to make sure the code is working and that lattice points aren't being written on top of each other. !This is primarily a debugging statement if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then stop "Multiple atoms share same position in lat point array, this shouldn't happen" else lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true. dof = dof + 1 end if end do !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)) r = r_interp(:,j)/lapa(ilat) r = matmul(fcc_inv,matmul(ori_inv,r)) do k = 1, 3 r_lat(k) = nint(r(k)) end do !Do a check to make sure the code is working and that lattice points aren't being written on top of each !other. This is primarily a debugging statement if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then stop "Multiple atoms/interpolated atoms share same position in lat point array, this shouldn't happen" else lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true. dof = dof + 1 end if end do end do !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 bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10 bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10 if (max_remesh) then max_loops = (remesh_size-3)/2 else max_loops = 1 end if do j = 1, max_loops working_esize = remesh_size - 2*(j-1) ele = (working_esize-1)*cubic_cell zloop: do iz = 1, bd_in_array(3) yloop: do iy = 1, bd_in_array(2) xloop: do ix = 1, bd_in_array(1) if (lat_points(ix, iy,iz)) then r_new_node(:,:,:) = 0.0_dp !Check to see if the element overshoots the bound if (iz+working_esize-1 > bd_in_array(3)) then exit zloop else if (iy+working_esize-1 > bd_in_array(2)) then cycle zloop else if (ix+working_esize-1 > bd_in_array(1)) then cycle yloop end if if (all(lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1))) then do inod = 1, ng_node(ilat) vlat = ele(:,inod) + (/ix, iy, iz /) do i = 1, 3 vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 end do r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat) end do 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(0,remesh_ele_type, working_esize, ilat, & sbox_atom(atom_index(1)),r_new_node) end if end if end do xloop end do yloop end do zloop end do !Now we have to add any leftover lattice points as atoms do iz = 1, bd_in_array(3) do iy=1, bd_in_array(2) do ix = 1, bd_in_array(1) if(lat_points(ix,iy,iz)) then vlat = (/ ix, iy, iz /) do i = 1, 3 vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 end do lat_points(ix,iy,iz) = .false. r = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat) new_atom=new_atom+1 call add_atom(0,basis_type(1,ilat), is, r) end if end do end do end do deallocate(lat_points) end do end do !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 !This subroutine deletes all atoms/elements within a group print *, "Deleting group containing ", group_atom_num, " atoms and ", group_ele_num, " elements." !Delete atoms call delete_atoms(group_atom_num, atom_index) !Delete elements call delete_elements(group_ele_num, element_index) end subroutine delete_group subroutine change_group_type !This subroutine changes all atoms and nodes at atoms within a group to a specific type integer :: i, j, ltype,ibasis, inod, basis_type(10) print *, "Changing ", group_atom_num, " atoms and ", group_ele_num, " elements to atom type ", group_type !Change all atom group types do i = 1, group_atom_num j = atom_index(i) type_atom(j) = group_type end do !Map to a new lattice type for all element do i =1, group_ele_num j = element_index(i) ltype = lat_ele(j) do ibasis = 1, basisnum(ltype) basis_type(ibasis) = group_type end do call lattice_map(basisnum(ltype), basis_type, ng_node(ltype), lapa(ltype), lat_ele(j)) end do end subroutine change_group_type subroutine split_group_elements ! end subroutine split_group_elements function in_group(r) !This subroutine determines if a point is within the group boundaries real(kind=dp), intent(in) :: r(3) real(kind=dp) :: rnorm integer :: dim3, i logical :: in_group select case(trim(adjustl(shape))) case('block') in_group=in_block_bd(r,block_bd) case('wedge') in_group = in_wedge_bd(r,vertices) case('notch') do i = 1, 3 if (.not.((dim1==i).or.(dim2==i))) dim3 = i end do in_group = in_wedge_bd(r,vertices) !Do a check to make sure the wedge isn't used if it should be the tip radius if (bwidth>0) then if (r(dim1) < centroid(dim1)) in_group = .false. end if rnorm = sqrt((r(dim1) - centroid(dim1))**2 + (r(dim3)-centroid(dim3))**2) in_group = in_group.or.(rnorm < radius) case('sphere') rnorm = norm2(r(:) - centroid(:)) if (rnorm <= radius) then in_group = .true. else in_group = .false. end if case('shell') rnorm = norm2(r(:) - centroid(:)) if ((rnorm >= radius).and.(rnorm<=(radius + shell_thickness))) then in_group = .true. else in_group = .false. end if case('all') in_group = .true. end select end function in_group function in_group_ele(esize, elat, rn) !This figures out if the elements are in the group boundaries real(kind=dp), intent(in) :: rn(3,max_basisnum, max_ng_node) integer, intent(in) :: esize, elat logical :: in_group_ele integer :: i, inod, ibasis, node_in_out(max_ng_node) real(kind=dp) :: r_center(3) in_group_ele=.false. if(trim(adjustl(shape)) == 'shell') then node_in_out(:) = -1 !First calculate whether each element node is within the shell region, inside the shell sphere, or outside the !shell region nodeloop:do inod = 1, ng_node(elat) r_center(:)=0.0_dp do ibasis = 1, basisnum(elat) r_center(:)= r_center(:) + rn(:,ibasis,inod)/basisnum(elat) end do if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(size_ele(i)/=notsize)) then node_in_out(inod) = 2 exit nodeloop end if shape ='sphere' if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then node_in_out(inod) = 1 else node_in_out(inod) = 0 end if shape='shell' end do nodeloop !If any nodes are within the shell region, or if the shell region interescts an element then add it to the group if(any(node_in_out == 2).or.(any(node_in_out==1).and.(any(node_in_out==0)))) then in_group_ele=.true. return end if else if(.not.(group_nodes)) then r_center(:) = 0.0_dp do inod = 1, ng_node(elat) do ibasis = 1, basisnum(elat) r_center = r_center + rn(:,ibasis,inod)/(basisnum(elat)*ng_node(elat)) end do end do if ((in_group(r_center).neqv.flip).and.(esize/= notsize)) then in_group_ele=.true. return end if else if(group_nodes) then r_center(:) = 0.0_dp do inod = 1, ng_node(elat) do ibasis = 1, basisnum(elat) if ((in_group(rn(:,ibasis,inod)).neqv.flip).and.(esize/=notsize)) then in_group_ele=.true. return end if end do end do end if end function in_group_ele end module opt_group