diff --git a/INSTALL.md b/INSTALL.md new file mode 100644 index 0000000..e760b93 --- /dev/null +++ b/INSTALL.md @@ -0,0 +1,16 @@ +# Install CACmb + +In order to build CACmb switch to the src repository and run: + +``` +make +``` + +Once this command is complete you can run: + +``` +sudo make install +``` + +in order to copy the installation to `/usr/local/bin`. + diff --git a/README.md b/README.md index cf7fcfb..5d45ef3 100644 --- a/README.md +++ b/README.md @@ -113,6 +113,26 @@ This mode merges multiple data files and creates one big simulation cell. The pa `dim` - the dimension they are to be stacked along, can be either `x`, `y`, or `z`. If the argument `none` is passed then the cells are just overlaid. Future options will include a delete overlap command. +**Additional options:** + +**Shift** + +``` +shift x y z +``` + +If the shift command is passed to mode merge then each file after the first file in the merge command is displaced by the vector `[x, y, z]`. This is additive so if you are merging three files and this command is passed then the second file is shifted by `[x,y,z]` and the third file is shifted by `2*[x,y,z]`. + +Example: `cacmb --merge z 2 Cu.mb Cu2.mb Cu3.mb Cumerged.mb shift 2 0 0` will shift the atomic and element positions in the `Cu2.mb` file by `[2,0,0]` and the positions in `Cu3.mb` by `[4,0,0]`. + +**Wrap** + +``` +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. + ## 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. @@ -120,7 +140,7 @@ Options are additional portions of code which have additional functionality. Opt ### Option dislgen ``` --dislgen [ijk] [hkl] x y z char_angle poisson +-dislgen [ijk] [hkl] x y z burgers char_angle poisson ``` @@ -131,6 +151,8 @@ This options adds an arbitrarily oriented dislocation into your model based on u `[hkl]` - The vector for the slip plane +`burgers` - The magnitude of the burgers vector for the dislocation to be inserted + `x y z` - The position of the dislocation centroid `char_angle` - Character angle of the dislocation (0 is screw and 90 is edge) @@ -153,4 +175,54 @@ This option deletes vacancies on a plane which when minimized should result in a `bx by bz` - The burgers vector for the dislocation -`poisson` - Poisson ratio for continuum solution \ No newline at end of file +`poisson` - Poisson ratio for continuum solution + +### Option Group + +`-group select_type group_shape shape_arguments additional keywords` + +This option selects a group of either elements, nodes, atoms and applies some transformation to them. + +`select_type` - Either `nodes`, `atoms`, `elements`, `nodes/atoms`, `all`. When using the option `nodes` only nodes which are within the group are selected, `elements` selects elements based on whether the element center is within the group, `nodes/atoms` selects both nodes and atoms for the group. `all` selects elements based on the element center and atoms based on their position. + +`group_shape` - Specifies what shape the group takes and dictates which options must be passed. Each shape requires different arguments and these arguments are represented by the placeholder `shape_arguments`. The accepted group shapes and arguments are below: + +*Block:* + +`-group nodes block xlo xhi ylo yhi zlo zhi` + +This selects a group residing in a block with edges perpendicular to the simulation cell. The block boundaries are given by `xlo xhi ylo yhi zlo zhi`. + +`additional keywords`- Represents the various transformations which can be performed on a group. These additional keywords are given below. + +**Displace** + +``` +displace x y z +``` + +This is similar to the mode merge `-shift` argument. This simply shift atoms and elements within the group by a vector `[x,y,z]`. + +**Wrap** + +``` +wrap +``` + +This command wraps atoms back into the simulation cell as though periodic boundary conditions are being used. + +**Remesh** + +``` +remesh esize +``` + +This command remeshes the atoms/elements within the group to the new element size `esize`. Currently only accepts an `esize` of 2 which refines it to full atomistics. + +### Option overwrite + +``` +-ow +``` + +If this option is passed then all files are automatically overwritten without asking the user. \ No newline at end of file diff --git a/src/Makefile b/src/Makefile index 0cdd008..a3f926f 100644 --- a/src/Makefile +++ b/src/Makefile @@ -2,8 +2,8 @@ FC=ifort FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin #FFLAGS=-mcmodel=large -Ofast -no-wrap-margin MODES=mode_create.o mode_merge.o mode_convert.o -OPTIONS=opt_disl.o -OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o call_option.o $(MODES) $(OPTIONS) +OPTIONS=opt_disl.o opt_group.o +OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES) $(OPTIONS) call_option.o .SUFFIXES: .SUFFIXES: .c .f .f90 .F90 .o @@ -38,6 +38,6 @@ atoms.o subroutines.o testfuncs.o box.o : functions.o main.o io.o $(MODES) $(OPTIONS) : elements.o call_mode.o : $(MODES) call_option.o : $(OPTIONS) -$(MODES) io.o: atoms.o box.o +$(MODES) $(OPTIONS) subroutines.o io.o : atoms.o box.o $(MODES) main.o : io.o -testfuncs.o elements.o mode_create.o opt_disl.o: subroutines.o +testfuncs.o elements.o mode_create.o $(OPTIONS): subroutines.o diff --git a/src/call_option.f90 b/src/call_option.f90 new file mode 100644 index 0000000..60551af --- /dev/null +++ b/src/call_option.f90 @@ -0,0 +1,22 @@ +subroutine call_option(option, arg_pos) + use parameters + use opt_disl + use opt_group + implicit none + + integer, intent(inout) :: arg_pos + character(len=100), intent(in) :: option + + select case(trim(adjustl(option))) + case('-dislgen', '-disloop') + call dislocation(option, arg_pos) + case('-group') + call group(arg_pos) + case('-ow') + arg_pos = arg_pos + 1 + continue + case default + print *, 'Option ', trim(adjustl(option)), ' is not currently accepted. Skipping to next argument' + arg_pos = arg_pos + 1 + end select +end subroutine call_option \ No newline at end of file diff --git a/src/elements.f90 b/src/elements.f90 index b4e59be..9aa5f54 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -223,6 +223,8 @@ module elements real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node) ele_num = ele_num + 1 + !Check to see if we need to grow the arrays + call grow_ele_arrays(1,0) type_ele(ele_num) = type size_ele(ele_num) = size lat_ele(ele_num) = lat @@ -239,6 +241,8 @@ module elements real(kind=dp), intent(in), dimension(3) :: r atom_num = atom_num+1 + !Check to see if we need to grow the arrays + call grow_ele_arrays(0,1) type_atom(atom_num) = type r_atom(:,atom_num) = r(:) @@ -356,7 +360,7 @@ module elements end do end select - if (ia /= esize**3) then + if (ia /= bnum*esize**3) then print *, "Incorrect interpolation" stop 3 end if @@ -403,4 +407,34 @@ module elements atom_num = atom_num - 1 end do end subroutine + + subroutine delete_elements(num, index) + !This subroutine deletes elements from the element array + integer, intent(in) :: num + integer, intent(inout), dimension(num) :: index + + integer :: i, j + + call heapsort(index) + + !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 + if(index(i) == ele_num) then + r_node(:,:,:,index(i)) = 0.0_dp + type_ele(index(i)) ='' + size_ele(index(i)) = 0 + lat_ele(index(i)) = 0 + sbox_ele(index(i)) = 0 + else + node_num = node_num - ng_node(lat_ele(index(i))) + r_node(:,:,:,index(i)) = r_node(:,:,:,ele_num) + type_ele(index(i)) = type_ele(ele_num) + size_ele(index(i)) = size_ele(ele_num) + lat_ele(index(i)) = lat_ele(ele_num) + sbox_ele(index(i)) = sbox_ele(ele_num) + end if + ele_num = ele_num - 1 + end do + end subroutine delete_elements end module elements \ No newline at end of file diff --git a/src/io.f90 b/src/io.f90 index e5b098d..26d4250 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -9,7 +9,7 @@ module io integer :: outfilenum = 0, infilenum = 0 character(len=100) :: outfiles(10), infiles(10) - + logical :: force_overwrite public contains @@ -38,7 +38,7 @@ module io !Check to see if file exists, if it does then ask user if they would like to overwrite the file inquire(file=trim(temp_outfile), exist=file_exists) - if (file_exists) then + if (file_exists.and.(.not.(force_overwrite))) then if (overwrite == 'r') print *, "File ", trim(temp_outfile), " already exists. Would you like to overwrite? (Y/N)" read(*,*) overwrite if((scan(overwrite, "n") > 0).or.(scan(overwrite, "N") > 0)) then @@ -194,6 +194,7 @@ module io select case(trim(adjustl(type_ele(i)))) case('fcc') do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 + call apply_periodic(r_interp(:,iatom)) write(11, '(2i16, 3f23.15)') atom_num+iatom, type_interp(iatom), r_interp(:,iatom) end do end select @@ -362,7 +363,7 @@ module io !NOTE: This code doesn't work for arbitrary number of basis atoms per node. It assumes that the !each element has only 1 atom type at the node. character(len=100), intent(in) :: file - integer :: interp_max, i, j, lat_size, inod, ibasis, ip + integer :: interp_max, i, j, lat_size, inod, ibasis, ip, unique_index(10), unique_num real(kind=dp) :: box_vec(3) 1 format('time' / i16, f23.15) @@ -398,6 +399,8 @@ module io select case(max_ng_node) case(8) interp_max = (max_esize)**3 + case default + interp_max = 0 end select write(11,20) interp_max write(11,3) node_num @@ -441,14 +444,21 @@ module io !write the element information if(ele_num > 0) then write(11,12) - do i = 1, lattice_types - do j = 1, ele_num - if (lat_ele(j) == i) then - lat_size = size_ele(j)-1 - exit + !First figure out all of the unique element types + unique_num = 0 + unique_index(:) = 0 + eleloop:do i = 1, ele_num + do j =1 , unique_num + if ( ( size_ele(i) == size_ele( unique_index(j) ) ).and. & + ( lat_ele(i) == lat_ele(unique_index(j)) ) ) then + cycle eleloop end if - end do - write(11,'(3i16)') i, lat_size, basis_type(1,i) + end do + unique_num = unique_num + 1 + unique_index(unique_num) = i + end do eleloop + do i = 1, unique_num + write(11,'(3i16)') i, size_ele(i)-1, basis_type(1,i) end do ip = 0 write(11,13) @@ -599,8 +609,8 @@ module io real(kind=dp), dimension(3), intent(in) :: displace real(kind = dp), dimension(6), intent(out) :: temp_box_bd - integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, & - new_type_to_type(10), new_lattice_types, sbox + integer :: i, j, k, l, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, & + new_type_to_type(10), new_lattice_types, sbox, new_lattice_map(10) character(len=100) :: etype real(kind=dp) :: r(3), newdisplace(3) real(kind=dp), allocatable :: r_innode(:,:,:) @@ -617,6 +627,7 @@ module io temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i) end do + call grow_box(temp_box_bd) !Read in the number of sub_boxes and allocate the variables read(11, *) n @@ -632,7 +643,11 @@ module io read(11,*) ((sub_box_ori(j, k, sub_box_num+i), j = 1, 3), k = 1, 3) !Read in subbox boundaries read(11,*) sub_box_bd(:,sub_box_num+i) - sub_box_bd(:,sub_box_num+i) = sub_box_bd(:, sub_box_num+i) + displace(:) + + do j = 1, 3 + sub_box_bd(2*j-1,sub_box_num+i) = sub_box_bd(2*j-1, sub_box_num+i) + displace(j) + sub_box_bd(2*j,sub_box_num+i) = sub_box_bd(2*j, sub_box_num+i) + displace(j) + end do !Read in sub_box_array_bd read(11,*) ((sub_box_array_bd(j, k, sub_box_num+i), j = 1, 2), k = 1, 2) end do @@ -661,14 +676,41 @@ module io do i = 1, basisnum(j) basis_type(i,j) = new_type_to_type(basis_type(i,j)) end do - end do + end do !Read the lattice parameters for every lattice type read(11,*) (lapa(i), i = lattice_types+1, lattice_types+new_lattice_types) + + !Now we loop over all new lattice types and check to see if they are exactly the same as any old lattice types + !If they are then we map the new type to the old type. + k = lattice_types + 1 + new_lattice_map(:) = 0 + new_loop:do i = 1, new_lattice_types + old_loop:do j = 1, lattice_types + !First check all the lattice level variables + if ((basisnum(i) == basisnum(j)).and. & + (ng_node(i) == ng_node(j))) then + !Now check the basis level variables + do ibasis =1, basisnum(i) + if(basis_type(ibasis,lattice_types+i) /= basis_type(ibasis,j)) then + cycle old_loop + end if + end do + new_lattice_map(i) = j + cycle new_loop + end if + end do old_loop + new_lattice_map(i) = k + k = k+1 + end do new_loop + !Read number of elements and atoms and allocate arrays read(11, *) in_atoms, in_eles call grow_ele_arrays(in_eles, in_atoms) allocate(r_innode(3,max_basisnum, max_ng_node)) + print *, "Read in ", in_eles, " elements and ", in_atoms, " atoms from ", trim(adjustl(file)) + print *, "New box dimensions are: ", box_bd + !Read the atoms do i = 1, in_atoms read(11,*) j, type, r(:) @@ -677,15 +719,14 @@ module io !Read the elements do i = 1, in_eles - read(11, *) n, type, size, sbox, etype + read(11, *) l, type, size, sbox, etype do inod = 1, ng_node(type) do ibasis =1, basisnum(type) read(11,*) j, k, r_innode(:, ibasis, inod) r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace end do end do - type = type + lattice_types - call add_element(etype, size, type, sbox+n, r_innode) + call add_element(etype, size, new_lattice_map(type), sbox+n, r_innode) end do !Close the file being read @@ -693,7 +734,7 @@ module io !Only increment the lattice types if there are elements, if there are no elements then we !just overwrite the arrays - if(in_eles > 0) lattice_types = lattice_types + new_lattice_types + if(in_eles > 0) lattice_types = maxval(new_lattice_map) sub_box_num = sub_box_num + n diff --git a/src/main.f90 b/src/main.f90 index ed34c53..31f5f58 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,33 +1,59 @@ program main - !**************************** CACmb ******************************* - !* CAC model building toolkit * - ! ____________ * - ! / / * - ! / / * - ! /___________/ * - ! _|_ _|_ _|____________ * - ! / / * - ! / / * - ! /___________/ * - ! * - !******************************************************************* + !**************************** CACmb ******************** + !* CAC model building toolkit * + !* ____________ * + !* / / * + !* / / * + !* /___________/ * + !* _|_ _|_ _|____________ * + !* / / * + !* / / * + !* /___________/ * + !* * + !******************************************************** use parameters use elements use io - + integer :: i, end_mode_arg, arg_num, arg_pos character(len=100) :: argument + !Print introduction text + print *, '*********************** CACmb *********************' + print *, '* CAC model building toolkit *' + print *, '* _______ *' + print *, '* / / *' + print *, '* / / *' + print *, '* /______ / *' + print *, '* _|_ _|_ _|_______ *' + print *, '* / / *' + print *, '* / / *' + print *, '* /______ / *' + print *, '* *' + print *, '****************************************************' + !Call initialization functions call lattice_init call box_init + force_overwrite=.false. + end_mode_arg = 0 ! Command line parsing arg_num = command_argument_count() + !Check to see if overwrite flag is passed + do i = 1, arg_num + call get_command_argument(i,argument) + select case(trim(adjustl(argument))) + + case('-ow') + force_overwrite = .true. + print *, "Overwrite flag passed, output files will be overwritten" + end select + end do !Determine if a mode is being used and what it is. The first argument has to be the mode !if a mode is being used call get_command_argument(1, argument) diff --git a/src/mode_convert.f90 b/src/mode_convert.f90 index 1f8a6bd..32ded60 100644 --- a/src/mode_convert.f90 +++ b/src/mode_convert.f90 @@ -13,8 +13,6 @@ module mode_convert integer, intent(out) :: arg_pos character(len=100) :: infile, outfile real(kind = dp) :: temp_box_bd(6) - !We have to allocate the element and atom arrays with a size of 1 for the read in code to work - call alloc_ele_arrays(1,1) !First read in the file call get_command_argument(2, infile) call get_in_file(infile) diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 9463e92..c3e6bcc 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -30,6 +30,8 @@ module mode_create integer :: i, ibasis, inod real(kind=dp), allocatable :: r_node_temp(:,:,:) + print *, '-----------------------Mode Create---------------------------' + !Initialize default parameters orient = reshape((/ 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp /), shape(orient)) cell_mat(:,:)=0.0_dp @@ -91,6 +93,8 @@ module mode_create end do else + print *, "Creating 1 element" + call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) !If the user doesn't pass any build instructions than we just put the cell mat into the element_array call alloc_ele_arrays(1,0) @@ -110,6 +114,7 @@ module mode_create !If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays if(dup_flag.or.dim_flag) then + !Call the build function with the correct transformation matrix select case(trim(adjustl(element_type))) case('fcc') @@ -123,6 +128,8 @@ module mode_create !Now that it is built multiply by the lattice parameter box_bd = box_bd*lattice_parameter + print *, "Using mode create, ", lat_ele_num, " elements are created and ", lat_atom_num*basisnum(1), " atoms are created." + !Allocate variables call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1)) if(lat_atom_num > 0) then diff --git a/src/mode_merge.f90 b/src/mode_merge.f90 index cfb4273..14a015d 100644 --- a/src/mode_merge.f90 +++ b/src/mode_merge.f90 @@ -8,7 +8,9 @@ module mode_merge use elements character(len=4) :: dim - integer :: in_num + integer :: in_num, new_starts(2) + real(kind=dp) :: shift_vec(3) + logical :: wrap, shift_flag public contains @@ -18,12 +20,26 @@ module mode_merge integer :: i real(kind=dp) :: displace(3), temp_box_bd(6) + print *, '-----------------------Mode Merge---------------------------' + + wrap = .false. + shift_flag = .false. + + shift_vec(:) = 0.0_dp !First we parse the merge command call parse_command(arg_pos) !Now loop over all files and stack them do i = 1, in_num displace(:) = 0.0_dp + + !The new starts variable dictate where in the atom and element array each new + !file starts. This is used for additional options that can be applied to solely + !these new atoms/elements that are read in. + + new_starts(1) = atom_num + 1 + new_starts(2) = ele_num + 1 + if ((i==1).or.(trim(adjustl(dim)) == 'none')) then call read_in(i, displace, temp_box_bd) call grow_box(temp_box_bd) @@ -40,6 +56,8 @@ module mode_merge call read_in(i, displace, temp_box_bd) call grow_box(temp_box_bd) end if + + if(shift_flag) call shift(new_starts, i) end do return @@ -89,7 +107,16 @@ module mode_merge !Choose what to based on what the option string is select case(trim(textholder)) - + case('shift') + shift_flag = .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, *) shift_vec(i) + end do + case('wrap') + wrap = .true. case default !If it isn't an available option to mode merge then we just exit exit @@ -97,4 +124,54 @@ module mode_merge end do end subroutine parse_command + + subroutine shift(array_start, filenum) + !This subroutine applies a shift to newly added atoms and elements. It also wraps the atoms + !if the user provides the wrap flag + + integer, dimension(2), intent(in) :: array_start + integer, intent(in) :: filenum + + integer :: i, j, ibasis, inod + real(kind=dp), dimension(3) :: current_shift + + + !Calculate the current shift which is the filenum-1 multiplied by the user specified shift + current_shift = (filenum-1)*shift_vec + + print *, "Atoms/elements from file ", trim(adjustl(infiles(filenum))), " are shifted by ", current_shift + + !First shift all the atoms + do i = array_start(1), atom_num + r_atom(:,i) = r_atom(:,i) + current_shift + end do + + !Now shift all the elements + 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) + current_shift + end do + end do + end do + + !Now we check if we have to wrap the atoms, nodes are not wrapped. For elements the periodic + !boundary conditions are applied in the actual CAC codes + if(wrap) then + do i = array_start(1), atom_num + call apply_periodic(r_atom(:,i)) + end do + + !If we don't include the wrap command then we have to increase the size of the box + else + do i = 1,3 + if (current_shift(i) < -lim_zero) then + box_bd(2*i-1) = box_bd(2*i-1) - current_shift(i) + else if (current_shift(i) > lim_zero) then + box_bd(2*i) = box_bd(2*i) + current_shift(i) + end if + end do + end if + + end subroutine shift end module mode_merge \ No newline at end of file diff --git a/src/opt_disl.f90 b/src/opt_disl.f90 index 4510ede..025d680 100644 --- a/src/opt_disl.f90 +++ b/src/opt_disl.f90 @@ -23,12 +23,12 @@ module opt_disl contains subroutine dislocation(option, arg_pos) - !Main calling function for all codes related to dislocations character(len=100), intent(in) :: option integer, intent(inout) :: arg_pos + print *, '--------------------Option Dislocation-----------------------' select case(trim(adjustl(option))) case('-dislgen') @@ -70,7 +70,11 @@ module opt_disl call parse_pos(i, textholder, centroid(i)) end do - print *, centroid + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing parameter in dislgen command" + read(textholder, *) b + arg_pos = arg_pos + 1 call get_command_argument(arg_pos, textholder, arglen) @@ -84,26 +88,6 @@ module opt_disl read(textholder, *) poisson arg_pos = arg_pos + 1 - call get_command_argument(arg_pos, lattice, arglen) - if (arglen==0) STOP "Missing lattice in dislgen command" - - arg_pos = arg_pos + 1 - call get_command_argument(arg_pos, textholder, arglen) - if (arglen==0) STOP "Missing lattice parameter in dislgen command" - read(textholder, *) lattice_parameter - - arg_pos = arg_pos + 1 - !Now set the vurgers vector based on the lattice paarameter and lattice type - select case(lattice) - case('fcc') - b = lattice_parameter / sqrt(2.0_dp) - ! case('bcc') - ! b = lattice_parameter * sqrt(3.0_dp) / 2.0_dp - case default - print *, 'Error: Lattice structure', lattice, ' is not accepted for dislgen option' - STOP - end select - end subroutine parse_dislgen subroutine dislgen @@ -113,6 +97,8 @@ module opt_disl real(kind=dp) :: ss_ori(3,3), ss_inv(3,3), be, bs, slipx(3), disp_transform(3,3), inv_transform(3,3), & actan, r(3), disp(3) + print *, "Dislocation with centroid ", centroid, " is inserted" + !Calculate screw and edge burgers vectors be = sin(char_angle*pi/180.0_dp)*b bs = cos(char_angle*pi/180.0_dp)*b @@ -168,7 +154,7 @@ module opt_disl do i = 1, ele_num do inod=1, ng_node(lat_ele(i)) do ibasis = 1, basisnum(lat_ele(i)) - r = r_node(:,ibasis,inod,i) + r = r_node(:,ibasis,inod,i)-centroid r = matmul(inv_transform, r) if (r(1) == 0) then actan = pi/2 @@ -270,6 +256,9 @@ module opt_disl real(kind = dp) :: perimeter, angle, theta, omega, xA(3), xB(3), xC(3), u(3) real(kind=dp), dimension(:,:), allocatable :: xloop !coordinate of points forming loop + + print *, "Dislocation loop with centroid ", centroid, " is inserted" + if(allocated(xLoop)) deallocate(xLoop) !Define new directions @@ -291,36 +280,52 @@ module opt_disl a3 = 3 end select - !Calculate loop perimeter - perimeter = 2.0_dp*pi*loop_radius - - !Define the number of points forming the loop - ! The following criteria are used as a trade-off between - ! good accuracy and computational efficiency: - ! - each dislocation segment should have a length of 5 angströms; - ! - the loop should contain at least 3 points - ! (for very small loops, this will result in segments shorter than 5 A); - ! - there should not be more than 100 points - ! (for very large loops, this will result in segments longer than 5 A). - Npoints = MAX( 3 , MIN( NINT(perimeter/5.d0) , 100 ) ) - - !angle between two consecutive points - theta = 2.0_dp*pi / dble(Npoints) - - !allocate xLoop - allocate(xLoop(Npoints,3)) - xLoop(:,:) = 0.0_dp - - !Calculate the position of each point in the loop - angle = 0.0_dp - do i = 1, size(xLoop,1) - xLoop(i,a1) = centroid(a1) + loop_radius*dcos(angle) - xLoop(i,a2) = centroid(a2) + loop_radius*dsin(angle) - xLoop(i,a3) = centroid(a3) - ! Increment angle for next point - angle = angle + theta - end do - + if(loop_radius < 0.0_dp) then + ALLOCATE(xLoop(4,3)) + xLoop(:,:) = 0.d0 + xLoop(1,a1) = centroid(1) - loop_radius + xLoop(1,a2) = centroid(2) - loop_radius + xLoop(1,a3) = centroid(3) + xLoop(2,a1) = centroid(1) + loop_radius + xLoop(2,a2) = centroid(2) - loop_radius + xLoop(2,a3) = centroid(3) + xLoop(3,a1) = centroid(1) + loop_radius + xLoop(3,a2) = centroid(2) + loop_radius + xLoop(3,a3) = centroid(3) + xLoop(4,a1) = centroid(1) - loop_radius + xLoop(4,a2) = centroid(2) + loop_radius + xLoop(4,a3) = centroid(3) + else + !Calculate loop perimeter + perimeter = 2.0_dp*pi*loop_radius + + !Define the number of points forming the loop + ! The following criteria are used as a trade-off between + ! good accuracy and computational efficiency: + ! - each dislocation segment should have a length of 5 angströms; + ! - the loop should contain at least 3 points + ! (for very small loops, this will result in segments shorter than 5 A); + ! - there should not be more than 100 points + ! (for very large loops, this will result in segments longer than 5 A). + Npoints = MAX( 3 , MIN( NINT(perimeter/5.d0) , 100 ) ) + + !angle between two consecutive points + theta = 2.0_dp*pi / dble(Npoints) + + !allocate xLoop + allocate(xLoop(Npoints,3)) + xLoop(:,:) = 0.0_dp + + !Calculate the position of each point in the loop + angle = 0.0_dp + do i = 1, size(xLoop,1) + xLoop(i,a1) = centroid(a1) + loop_radius*dcos(angle) + xLoop(i,a2) = centroid(a2) + loop_radius*dsin(angle) + xLoop(i,a3) = centroid(a3) + ! Increment angle for next point + angle = angle + theta + end do + end if !Now actually calculate the displacement created by a loop for every atom do i = 1, atom_num u = 0.0_dp diff --git a/src/opt_group.f90 b/src/opt_group.f90 new file mode 100644 index 0000000..6d129f1 --- /dev/null +++ b/src/opt_group.f90 @@ -0,0 +1,259 @@ +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 + character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape + real(kind=dp) :: block_bd(6), disp_vec(3) + logical :: displace, wrap + + 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 *, '-----------------------Option Group-------------------------' + + group_ele_num = 0 + group_atom_num = 0 + remesh_size=0 + if(allocated(element_index)) deallocate(element_index) + if(allocated(atom_index)) deallocate(atom_index) + + call parse_group(arg_pos) + + call get_group + + !Now call the transformation functions for the group + if(displace) call displace_group + + if(remesh_size > 0) call remesh_group + end subroutine group + + subroutine parse_group(arg_pos) + !Parse the group command + integer, intent(inout) :: arg_pos + + integer :: i,arglen + character(len=100) :: textholder + + !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, nodes, 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 dislgen command" + call parse_pos(int((i+1)/2), textholder, block_bd(i)) + end do + + 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('wrap') + wrap = .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 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, inod, ibasis + integer, allocatable :: resize_array(:) + real(kind=dp) :: r_center(3) + + select case(trim(adjustl(shape))) + case('block') + + print *, "Group has block shape with boundaries: ", block_bd + + !Allocate variables to arbitrary size + 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 + r_center(:) = 0.0_dp + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + r_center = r_center + r_node(:,ibasis,inod,i)/(basisnum(lat_ele(i))*ng_node(lat_ele(i))) + end do + end do + + if (in_block_bd(r_center, block_bd)) 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 + 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_block_bd(r_atom(:,i),block_bd)) 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 + end select + end select + + 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 + + !Now either apply periodic boundaries if wrap command was passed or adjust box dimensions + !Now we check if we have to wrap the atoms, nodes are not wrapped. For elements the periodic + !boundary conditions are applied in the actual CAC codes + if(wrap) then + do i = 1, atom_num + call apply_periodic(r_atom(:,i)) + end do + + !If we don't include the wrap command then we have to increase the size of the box + else + 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 remesh_group + !This command is used to remesh the group to a desired element size + + 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 and remeshing to elements are different processes so check which code we need to run + select case(remesh_size) + + !Refining to atoms + case(2) + 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(type_interp(j), 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 + !Remeshing to elements, currently not available + case default + print *, "Remeshing to elements is currently not available. Please refine to atoms by passing a remsh size of 2" + end select + + + end subroutine remesh_group +end module opt_group \ No newline at end of file diff --git a/src/subroutines.f90 b/src/subroutines.f90 index 330fbca..9fd159c 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -178,22 +178,36 @@ module subroutines real(kind=dp), intent(out) :: pos !The output parsed position value integer :: iospara + + iospara = 0 if(trim(adjustl(pos_string)) == 'inf') then pos=box_bd(2*i) else if(trim(adjustl(pos_string)) == '-inf') then pos=box_bd(2*i-1) else if ((index(pos_string,'-') > 0).and.(index(pos_string,'inf')>0)) then !Now extract the number we are reducing from infinity - read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos + if(index(pos_string,'inf') < index(pos_string,'-')) then + read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos + else + read(pos_string(1:index(pos_string,'-')-1), *, iostat=iospara) pos + end if pos = box_bd(2*i) - pos else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then !Now extract the number we are reducing from infinity - read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos + if(index(pos_string,'inf') < index(pos_string,'+')) then + read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos + else + read(pos_string(1:index(pos_string,'+')-1), *, iostat=iospara) pos + end if pos = box_bd(2*i-1) + pos else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then !Now extract the number we are reducing from infinity - read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos - pos = (box_bd(2*i)-box_bd(2*i-1))*pos + if(index(pos_string,'inf') < index(pos_string,'*')) then + read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos + else + read(pos_string(1:index(pos_string,'*')-1), *, iostat=iospara) pos + end if + pos = (box_bd(2*i)-box_bd(2*i-1))*pos + box_bd(2*i-1) else read(pos_string, *, iostat=iospara) pos end if @@ -249,4 +263,24 @@ module subroutines end do end subroutine siftdown + + + subroutine apply_periodic(r) + !This function checks if an atom is outside the box and wraps it back in. This is generally used when some + !kind of displacement is applied but the simulation cell is desired to be maintained as the same size. + + real(kind=dp), dimension(3), intent(inout) :: r + + integer :: j + real(kind=dp) ::box_len + do j = 1, 3 + box_len = box_bd(2*j) - box_bd(2*j-1) + if (r(j) > box_bd(2*j)) then + r(j) = r(j) - box_len + else if (r(j) < box_bd(2*j-1)) then + r(j) = r(j) + box_len + end if + end do + end subroutine + end module subroutines \ No newline at end of file