diff --git a/Numbered_element.png b/Numbered_element.png new file mode 100644 index 0000000..9930917 Binary files /dev/null and b/Numbered_element.png differ diff --git a/README.md b/README.md index 5d45ef3..d387a8d 100644 --- a/README.md +++ b/README.md @@ -214,10 +214,26 @@ This command wraps atoms back into the simulation cell as though periodic bounda **Remesh** ``` -remesh esize +remesh esize lattice_parameter lattice_type ``` -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. +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. When remeshing to atomistics the group can contain any orientations of elements but when remeshing to different finite elements, the group must contain all atoms/elements with the same orientation. `lattice_parameter` is the lattice parameter for the elements within the group and `lattice_type` is the lattice type (integer) that these new elements will be assigned to. + +**Max** + +``` +max +``` + +This command attempts to reduce the degrees of freedom in the model by replacing them with graded elements. This code works by starting at elements with size `esize` and then checks all degrees of freedom to see which ones can be replaced by inserting the element. It then iterates over elements of `esize-2` to `esize=2` which is full atomic resolution. + +**Delete** + +``` +delete +``` + +This command deletes all selected atoms and elements within the group. ### Option overwrite @@ -225,4 +241,70 @@ This command remeshes the atoms/elements within the group to the new element siz -ow ``` -If this option is passed then all files are automatically overwritten without asking the user. \ No newline at end of file +If this option is passed then all files are automatically overwritten without asking the user. + +### Option boundary + +``` +-boundary box_bc +``` + +This allows the user to specify the boundary conditions for the model being outputted. The format is a 3 character string with `p` indicating periodic and `s` indicating shrink-wrapped. + +**Example:** `-boundary psp` + +### Option Delete + +``` +-delete keywords +``` + +Delete requires the usage of additional keywords to specify which delete action will be taken. These additional keywords are below: + +**overlap** + +``` +-delete overlap rc_off +``` + +This command will delete all overlapping atoms within a specific cutoff radius `rc_off`. This currently does not affect elements. + +**** + +## Position Specification + +Specifying positions in cacmb can be done through a variety of ways. Examples of each format is shown below. + +`val` - Where `val` is a number, then that value in Angstroms is used as the position. As an example, `11.1` would be read in as a position of 11.1 $\AA$. + +`-inf` - This specifies the lower box boundary in the specific dimension. The vector `-inf -inf -inf` specifies the bottom corner of the simulation cell which also acts as the simulation cell origin. The vector `-inf 10 3` instead puts only the x position at the simulation cell origin. + +`inf` - Similar to `-inf` but references the upper boundary of the box in that dimension + +`inf-val` - Using a minus sign reduces the position from the **upper boundary** by `val`. `inf-10` would be at a distance of $10 \AA$ from the upper boundary in that dimension. + +`inf+val` - This increases the position from the **lower boundary**. `inf+10` would be a position $10\AA$ from the lower boundary within the cell. + +`inf*val` - This gives you a fractional position in the simulation cell. As an example `inf*0.5` gives you the center point of the simulation cell. + +`rand` - Returns a random position that lies within the simulation cell. + +`rand[val1:val2]` - returns a random position that lies within the range + +`rande[facenum]` - Returns a random position in an interelement boundary which is offset of the element face `facenum`. Face numbers are based on the which vertices comprise the face. Vertex numbers are shown in the figure below for the primitive fcc unit cell which is what the fcc rhombohedral element is based from. The face numbers are: + +Face 1: [1,2,3,4] + +Face 2: [1,2,6,5] + +Face 3: [2,3,7,6] + +Face 4: [3,4,8,7] + +Face 5: [1,4,8,5] + +Face 6: [5,6,7,8] + +Image for vertex numbers is: + +![](/home/alexselimov/Documents/CACmb/Numbered_element.png) \ No newline at end of file diff --git a/src/Makefile b/src/Makefile index a3f926f..9a3fced 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,8 +1,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 +FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays +#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays MODES=mode_create.o mode_merge.o mode_convert.o -OPTIONS=opt_disl.o opt_group.o +OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.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: @@ -40,4 +40,4 @@ call_mode.o : $(MODES) call_option.o : $(OPTIONS) $(MODES) $(OPTIONS) subroutines.o io.o : atoms.o box.o $(MODES) main.o : io.o -testfuncs.o elements.o mode_create.o $(OPTIONS): subroutines.o +testfuncs.o elements.o mode_create.o $(OPTIONS) $(MODES): subroutines.o diff --git a/src/call_option.f90 b/src/call_option.f90 index 60551af..1755d89 100644 --- a/src/call_option.f90 +++ b/src/call_option.f90 @@ -2,6 +2,9 @@ subroutine call_option(option, arg_pos) use parameters use opt_disl use opt_group + use opt_orient + use opt_delete + use box implicit none integer, intent(inout) :: arg_pos @@ -14,9 +17,21 @@ subroutine call_option(option, arg_pos) 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' + case('-wrap') + arg_pos = arg_pos + 1 + case('-orient') + call orient(arg_pos) + case('-unorient') + call unorient arg_pos = arg_pos + 1 + case('-boundary') + arg_pos=arg_pos+1 + call get_command_argument(arg_pos, box_bc) + arg_pos=arg_pos+1 + case('-delete') + call run_delete(arg_pos) + case default + print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.' + stop 3 end select end subroutine call_option \ No newline at end of file diff --git a/src/elements.f90 b/src/elements.f90 index 9aa5f54..f63b6b2 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -5,6 +5,7 @@ module elements use parameters use functions use subroutines + use box implicit none !Data structures used to represent the CAC elements. Each index represents an element @@ -17,6 +18,7 @@ module elements !Data structure used to represent atoms integer, allocatable :: type_atom(:)!atom type + integer, allocatable :: sbox_atom(:) real(kind =dp),allocatable :: r_atom(:,:) !atom position integer :: atom_num=0 !Number of atoms @@ -41,6 +43,8 @@ module elements integer :: basis_type(10,10) real(kind=dp) :: lapa(10) + !Additional module level variables we need + logical :: wrap_flag public contains @@ -60,11 +64,11 @@ module elements shape(fcc_cell)) !Now we create a list containing the list of vertices needed to describe the 6 cube faces - cubic_faces(:,1) = (/ 1, 4, 8, 5 /) - cubic_faces(:,2) = (/ 2, 3, 7, 6 /) - cubic_faces(:,3) = (/ 1, 2, 6, 5 /) + cubic_faces(:,1) = (/ 1, 2, 3, 4 /) + cubic_faces(:,2) = (/ 1, 2, 6, 5 /) + cubic_faces(:,3) = (/ 2, 3, 7, 6 /) cubic_faces(:,4) = (/ 3, 4, 8, 7 /) - cubic_faces(:,5) = (/ 1, 2, 3, 4 /) + cubic_faces(:,5) = (/ 1, 4, 8, 5 /) cubic_faces(:,6) = (/ 5, 6, 7, 8 /) !!Now initialize the fcc primitive cell @@ -151,7 +155,7 @@ module elements if(m > 0) then !Allocate atom arrays - allocate(type_atom(m), r_atom(3,m), stat=allostat) + allocate(type_atom(m), sbox_atom(m), r_atom(3,m), stat=allostat) if(allostat > 0) then print *, "Error allocating atom arrays in elements.f90 because of: ", allostat stop @@ -209,6 +213,11 @@ module elements temp_int(atom_size+1:) = 0 call move_alloc(temp_int, type_atom) + allocate(temp_int(m+atom_num+buffer_size)) + temp_int(1:atom_size) = sbox_atom + temp_int(atom_size+1:) = 0 + call move_alloc(temp_int, sbox_atom) + allocate(temp_real(3,m+atom_num+buffer_size)) temp_real(:,1:atom_size) = r_atom temp_real(:, atom_size+1:) = 0.0_dp @@ -235,9 +244,9 @@ module elements end subroutine add_element - subroutine add_atom(type, r) + subroutine add_atom(type, sbox, r) !Subroutine which adds an atom to the atom arrays - integer, intent(in) :: type + integer, intent(in) :: type, sbox real(kind=dp), intent(in), dimension(3) :: r atom_num = atom_num+1 @@ -245,6 +254,7 @@ module elements call grow_ele_arrays(0,1) type_atom(atom_num) = type r_atom(:,atom_num) = r(:) + sbox_atom(atom_num) = sbox end subroutine add_atom @@ -313,7 +323,7 @@ module elements real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: r_interp !Interpolated atomic positions !Internal variables - integer :: i, it, is, ir, ibasis, inod, ia, bnum, lat_type_temp + integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp real(kind=dp), allocatable :: a_shape(:) real(kind=dp) :: t, s, r @@ -390,7 +400,7 @@ module elements integer, intent(in) :: num integer, intent(inout), dimension(num) :: index - integer :: i, j + integer :: i call heapsort(index) @@ -413,7 +423,7 @@ module elements integer, intent(in) :: num integer, intent(inout), dimension(num) :: index - integer :: i, j + integer :: i call heapsort(index) @@ -436,5 +446,189 @@ module elements end if ele_num = ele_num - 1 end do + end subroutine delete_elements + + subroutine wrap_atoms + !This subroutine wraps atoms back into the simulation cell if they have exited for any reason + integer :: i + do i = 1, atom_num + call apply_periodic(r_atom(:,i)) + end do + end subroutine wrap_atoms + + subroutine def_new_box + !This subroutine calculates new box boundaries based on minimum and maximum nodal/atomic positions + integer :: i, j, inod, ibasis + + real(kind=dp) :: max_bd(3), min_bd(3) + + max_bd(:) = -huge(1.0_dp) + min_bd(:) = huge(1.0_dp) + + do i = 1, atom_num + do j = 1, 3 + if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + lim_zero + if (r_atom(j,i) < min_bd(j)) min_bd(j) = r_atom(j,i) - lim_zero + end do + end do + + do i = 1, ele_num + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + do j = 1, 3 + if (r_node(j,ibasis,inod,i) > max_bd(j)) max_bd(j) = r_node(j,ibasis,inod,i) + lim_zero + if (r_node(j,ibasis,inod,i) < min_bd(j)) min_bd(j) = r_node(j,ibasis,inod,i) -lim_zero + + end do + end do + end do + end do + + do j = 1, 3 + box_bd(2*j) = max_bd(j) + box_bd(2*j-1) = min_bd(j) + end do + end subroutine + + recursive subroutine parse_pos(i, pos_string, pos) + !This subroutine parses the pos command allowing for command which include inf + integer, intent(in) :: i !The dimension of the position + character(len=100), intent(in) :: pos_string !The position string + real(kind=dp), intent(out) :: pos !The output parsed position value + + integer :: iospara, face, ele, randsize + real(kind=dp) :: rand, rone, rtwo, rand_ele_pos(3) + character(len=100) :: cone, ctwo + + 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 (trim(adjustl(pos_string)) == 'rand') then + call random_number(rand) + pos = (box_bd(2*i)-box_bd(2*i-1))*rand + box_bd(2*i-1) + + else if (index(pos_string,'rande')>0) then + !First select a random element + call random_number(rand) + ele= 1 + floor(ele_num*rand) + !Now read the rest of the command which specifies the face we need and check + !to make sure it's an accepted face number + cone = pos_string(index(pos_string, '[')+1:index(pos_string, '[')+1) + + !Check to see if we also pass an element size, if it was passed then make sure our + !random element is the right size + if(index(pos_string,':') > 0) then + ctwo = pos_string(index(pos_string, ':')+1:index(pos_string,']')-1) + read(ctwo, *) randsize + do while(randsize /= size_ele(ele)) + call random_number(rand) + ele= 1 + floor(ele_num*rand) + end do + end if + + read(cone, *) face + if ((face < 1).or.(face > 6)) stop "Current face number must be 1,2,3,4,5,6. Please check documentation" + !Now get the position + call offset_pos(ele, face, rand_ele_pos) + pos = rand_ele_pos(i) + + else if (index(pos_string,'rand')>0) then + call random_number(rand) + cone = pos_string(index(pos_string, '[')+1:index(pos_string,':')-1) + call parse_pos(i, cone, rone) + ctwo = pos_string(index(pos_string, ':')+1:index(pos_string,']')-1) + call parse_pos(i, ctwo, rtwo) + pos = (rtwo - rone)*rand + rone + + else if ((index(pos_string,'-') > 0).and.(index(pos_string,'inf')>0)) then + !Now extract the number we are reducing from infinity + 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 + 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 + 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 + + if (iospara > 0) then + print *, "Error reading position argument ", trim(adjustl(pos_string)), ". Please reformat and try again." + end if + end subroutine parse_pos + + subroutine offset_pos(ie, iface, pos) + !This returns a position slightly offset from the center of an element face + !This is used to return a random position within an element discontinuity + + integer, intent(in) :: ie !Element index + integer, intent(in) :: iface !Face number between 1 and 6 + real(kind=dp), dimension(3), intent(out) :: pos !Position vector + + !Other variables we need + integer :: esize + real(kind=dp) :: orient(3,3), ori_inv(3,3), lp, r_cubic_node(3,8) + + !First find the offset vector for the face in the untransformed cubic cell + esize = size_ele(ie) + select case(iface) + case(1) + pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp /) + case(2) + pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /) + case(3) + pos = (/ (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) + case(4) + pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /) + case(5) + pos = (/ -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /) + case(6) + pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp /) + end select + + !Now transform it to real space and adjust it to the position of the element in the first node. + select case (trim(adjustl(type_ele(ie)))) + + case('fcc') + !First we have to extract the element lattice parameter + call matrix_inverse(sub_box_ori(:,:,sbox_ele(ie)),3,ori_inv) + r_cubic_node = r_node(:,1,:,ie) + r_cubic_node = matmul(fcc_inv,matmul(ori_inv,r_cubic_node)) + lp = (maxval(r_cubic_node(1,:))-minval(r_cubic_node(1,:)))/(esize-1) + pos = matmul(sub_box_ori(:,:,sbox_ele(ie)),matmul(fcc_mat,pos))*lp + pos = pos + r_node(:,1, 1, ie) + + case default + print *, trim(adjustl(type_ele(ie))), " is not a currently accepted element type for random element positions" + stop 3 + + end select + end subroutine + end module elements \ No newline at end of file diff --git a/src/io.f90 b/src/io.f90 index 26d4250..4aab723 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -8,7 +8,7 @@ module io implicit none integer :: outfilenum = 0, infilenum = 0 - character(len=100) :: outfiles(10), infiles(10) + character(len=100) :: outfiles(100), infiles(100) logical :: force_overwrite public @@ -87,6 +87,9 @@ module io call set_max_esize do i = 1, outfilenum + + print *, "Writing data out to ", trim(adjustl(outfiles(i))) + !Pull out the extension of the file and call the correct write subroutine select case(trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:)))) case('xyz') @@ -130,14 +133,14 @@ module io do i = 1, ele_num do inod = 1, ng_node(lat_ele(i)) do ibasis = 1, basisnum(lat_ele(i)) - write(11, '(a, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i) + write(11, '(i16, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i) end do end do end do !Write atom positions do i = 1, atom_num - write(11, '(a, 3f23.15)') type_atom(i), r_atom(:,i) + write(11, '(i16, 3f23.15)') type_atom(i), r_atom(:,i) end do !Finish writing @@ -307,6 +310,9 @@ module io 20 format('SCALARS lattice_type float', /& 'LOOKUP_TABLE default') +21 format('SCALARS esize float', /& + 'LOOKUP_TABLE default') + !First we write the vtk file containing the atoms open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind') @@ -355,6 +361,10 @@ module io do i = 1, ele_num write(11, '(i16)') lat_ele(i) end do + write(11,21) + do i = 1, ele_num + write(11, '(i16)') size_ele(i) + end do close(11) end subroutine @@ -363,7 +373,8 @@ 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, unique_index(10), unique_num + integer :: interp_max, i, j, inod, ibasis, ip, unique_index(50), unique_size(50), unique_num, & + etype real(kind=dp) :: box_vec(3) 1 format('time' / i16, f23.15) @@ -395,6 +406,22 @@ module io !Below writes the header information for the restart file + + !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 + unique_num = unique_num + 1 + unique_index(unique_num) = i + unique_size(unique_num) = size_ele(i) + end do eleloop + !Calculate the max number of atoms per element select case(max_ng_node) case(8) @@ -405,7 +432,7 @@ module io write(11,20) interp_max write(11,3) node_num write(11,19) max_ng_node - write(11,4) lattice_types + write(11,4) unique_num write(11,5) atom_num write(11,6) 1 !Grain_num is ignored write(11,16) lattice_types, atom_types @@ -444,26 +471,21 @@ module io !write the element information if(ele_num > 0) then write(11,12) - !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 - 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) + write(11,'(3i16)') i, size_ele(unique_index(i))-1, basis_type(1,lat_ele(unique_index(i))) end do ip = 0 write(11,13) do i = 1, ele_num - write(11, '(4i16)') i, lat_ele(i), 1, basis_type(1,lat_ele(i)) + !Figure out the ele type + do j = 1, unique_num + if ( unique_size(j) == size_ele(i)) then + etype = j + exit + endif + end do + write(11, '(4i16)') i, etype, 1, basis_type(1,lat_ele(i)) do inod = 1, ng_node(lat_ele(i)) do ibasis = 1, basisnum(lat_ele(i)) ip = ip + 1 @@ -520,7 +542,7 @@ module io !Write out atoms first do i = 1, atom_num - write(11,*) i, type_atom(i), r_atom(:,i) + write(11,*) i, type_atom(i), sbox_atom(i), r_atom(:,i) end do !Write out the elements, this is written in two stages, one line for the element and then 1 line for @@ -621,12 +643,18 @@ module io !Read in the box boundary and grow the current active box bd read(11, *) temp_box_bd(:) + print *, "displace", displace do i = 1, 3 - newdisplace(i) = displace(i) - temp_box_bd(2*i-1) + if (abs(displace(i)) > lim_zero) then + newdisplace(i) = displace(i) - temp_box_bd(2*i-1) + else + newdisplace(i)=displace(i) + end if temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i) temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i) end do - + + print *, "newdisplace", newdisplace call grow_box(temp_box_bd) !Read in the number of sub_boxes and allocate the variables read(11, *) n @@ -713,8 +741,9 @@ module io !Read the atoms do i = 1, in_atoms - read(11,*) j, type, r(:) - call add_atom(new_type_to_type(type), r+newdisplace) + read(11,*) j, type, sbox, r(:) + r = r+newdisplace + call add_atom(new_type_to_type(type), sbox+sub_box_num, r) end do !Read the elements @@ -726,7 +755,7 @@ module io r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace end do end do - call add_element(etype, size, new_lattice_map(type), sbox+n, r_innode) + call add_element(etype, size, new_lattice_map(type), sbox+sub_box_num, r_innode) end do !Close the file being read @@ -734,10 +763,11 @@ 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 = maxval(new_lattice_map) + lattice_types = maxval(new_lattice_map) sub_box_num = sub_box_num + n - + + call set_max_esize end subroutine read_mb end module io diff --git a/src/main.f90 b/src/main.f90 index 31f5f58..85b3b49 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -37,21 +37,29 @@ program main !Call initialization functions call lattice_init call box_init + call random_seed force_overwrite=.false. + wrap_flag = .false. end_mode_arg = 0 ! Command line parsing arg_num = command_argument_count() - !Check to see if overwrite flag is passed + !First check to see if certain commands are passed, these commands must be known before code + !is executed. do i = 1, arg_num call get_command_argument(i,argument) select case(trim(adjustl(argument))) + !This lets us know if we are overwriting all files case('-ow') force_overwrite = .true. print *, "Overwrite flag passed, output files will be overwritten" + + !This lets us know if we need to wrap atomic positions back into the cell + case('-wrap') + wrap_flag=.true. end select end do !Determine if a mode is being used and what it is. The first argument has to be the mode @@ -85,11 +93,14 @@ program main call call_option(argument, arg_pos) !Otherwise print that the argument is not accepted and move on else - print *, trim(adjustl(argument)), " is not accepted. Skipping to next argument" - arg_pos = arg_pos + 1 + print *, trim(adjustl(argument)), " is not an accepted command." + stop 3 end if end do + !If wrap flag was passed then call the wrap atoms command + if(wrap_flag) call wrap_atoms + !Check to make sure a file was passed to be written out and then write out ! Before building do a check on the file if (outfilenum == 0) then @@ -98,4 +109,4 @@ program main end if call write_out -end program main \ No newline at end of file +end program main diff --git a/src/mode_convert.f90 b/src/mode_convert.f90 index 32ded60..3bb8a3b 100644 --- a/src/mode_convert.f90 +++ b/src/mode_convert.f90 @@ -11,7 +11,7 @@ module mode_convert subroutine convert(arg_pos) !This subroutine converts a single input file from one format to another integer, intent(out) :: arg_pos - character(len=100) :: infile, outfile + character(len=100) :: infile real(kind = dp) :: temp_box_bd(6) !First read in the file call get_command_argument(2, infile) diff --git a/src/mode_create.f90 b/src/mode_create.f90 index c3e6bcc..3cbafa2 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -24,7 +24,6 @@ module mode_create subroutine create(arg_pos) ! Main subroutine which controls execution - character(len=100) :: textholder integer, intent(out) :: arg_pos integer :: i, ibasis, inod @@ -106,8 +105,8 @@ module mode_create end do end do do i = 1,3 - box_bd(2*i) = maxval(r_node_temp(i,:,:)) - box_bd(2*i-1) = origin(i) + box_bd(2*i) = maxval(r_node_temp(i,:,:))+10.0_dp**-6.0_dp + box_bd(2*i-1) = minval(r_node_temp(i,:,:)) - 10.0_dp**-6.0_dp end do call add_element(element_type, esize, 1, 1, r_node_temp) end if @@ -128,14 +127,15 @@ 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." + 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 do i = 1, lat_atom_num do ibasis = 1, basisnum(1) - call add_atom(basis_type(ibasis, 1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis)) + call add_atom(basis_type(ibasis, 1), 1, (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis)) end do end do deallocate(r_atom_lat) @@ -169,6 +169,7 @@ module mode_create integer :: ori_pos, i, j, arglen, stat character(len=100) :: textholder character(len=8) :: orient_string + logical :: isortho, isrighthanded !Pull out all required positional arguments @@ -205,24 +206,24 @@ module mode_create do i = 1, 3 call get_command_argument(arg_pos, orient_string, arglen) if (arglen==0) STOP "Missing orientation in orient command of mode create" + call parse_ori_vec(orient_string, orient(i,:)) arg_pos = arg_pos+1 - ori_pos=2 - do j = 1,3 - if (orient_string(ori_pos:ori_pos) == '-') then - ori_pos = ori_pos + 1 - read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j) - if (stat>0) STOP "Error reading orient value" - orient(i,j) = -orient(i,j) - ori_pos = ori_pos + 1 - else - read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j) - if(stat>0) STOP "Error reading orient value" - ori_pos=ori_pos + 1 - end if - end do - end do - + ! ori_pos=2 + ! do j = 1,3 + ! if (orient_string(ori_pos:ori_pos) == '-') then + ! ori_pos = ori_pos + 1 + ! read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j) + ! if (stat>0) STOP "Error reading orient value" + ! orient(i,j) = -orient(i,j) + ! ori_pos = ori_pos + 1 + ! else + ! read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j) + ! if(stat>0) STOP "Error reading orient value" + ! ori_pos=ori_pos + 1 + ! end if + ! end do + end do !If the duplicate command is passed then we extract the information on the new bounds. case('duplicate') @@ -284,6 +285,13 @@ module mode_create end select !Now normalize the orientation matrix orient = matrix_normal(orient,3) + !Now check these to make sure they are right handed and orthogonal + call check_right_ortho(orient, isortho, isrighthanded) + if (.not.isortho) then + stop "Directions in orient are not orthogonal" + else if (.not.isrighthanded) then + stop "Directions in orient are not righthanded" + end if !Set lattice_num to 1 and add the lattice_parameter to the elements module lattice paramter variable lattice_types = 1 @@ -309,7 +317,6 @@ module mode_create integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,8), rzero(3), & vlat(3), temp_lat(3,8), m, n, o real(kind=dp) :: v(3), temp_nodes(3,1,8) - real(kind=dp), allocatable :: resize_lat_array(:,:) logical, allocatable :: lat_points(:,:,:) logical :: node_in_bd(8) @@ -386,7 +393,7 @@ module mode_create outerloop: do iz = 1, bd_in_array(3) do iy = 1, bd_in_array(2) do ix = 1, bd_in_array(1) - node_in_bd(8) = .false. + node_in_bd(:) = .false. do inod = 1, 8 vlat = ele(:,inod) + (/ ix, iy, iz /) diff --git a/src/mode_merge.f90 b/src/mode_merge.f90 index 14a015d..afc9bba 100644 --- a/src/mode_merge.f90 +++ b/src/mode_merge.f90 @@ -10,7 +10,7 @@ module mode_merge character(len=4) :: dim integer :: in_num, new_starts(2) real(kind=dp) :: shift_vec(3) - logical :: wrap, shift_flag + logical :: shift_flag public contains @@ -22,7 +22,6 @@ module mode_merge print *, '-----------------------Mode Merge---------------------------' - wrap = .false. shift_flag = .false. shift_vec(:) = 0.0_dp @@ -115,8 +114,6 @@ module mode_merge 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 @@ -126,16 +123,16 @@ module mode_merge 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 + !This subroutine applies a shift to newly added atoms and elements. integer, dimension(2), intent(in) :: array_start integer, intent(in) :: filenum - integer :: i, j, ibasis, inod + integer :: i, ibasis, inod real(kind=dp), dimension(3) :: current_shift + character(len=3) :: alldims - + alldims = 'xyz' !Calculate the current shift which is the filenum-1 multiplied by the user specified shift current_shift = (filenum-1)*shift_vec @@ -155,19 +152,16 @@ module mode_merge 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 + if(.not.(wrap_flag)) then 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 + if (alldims(i:i) /= trim(adjustl(dim))) then + 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 + else if (alldims(i:i) == trim(adjustl(dim))) then box_bd(2*i) = box_bd(2*i) + current_shift(i) end if end do diff --git a/src/opt_delete.f90 b/src/opt_delete.f90 new file mode 100644 index 0000000..09cefb1 --- /dev/null +++ b/src/opt_delete.f90 @@ -0,0 +1,125 @@ +module opt_delete + + use parameters + use subroutines + use elements + + implicit none + + real(kind=dp) :: rc_off + + public + contains + + subroutine run_delete(arg_pos) + + integer, intent(inout) :: arg_pos + + rc_off = -lim_zero + !Main calling function for delete option. + print *, '-----------------------Option Delete------------------------' + + call parse_delete(arg_pos) + + if (rc_off > 0.0_dp) call delete_overlap + end subroutine run_delete + + subroutine parse_delete(arg_pos) + !Parse the delete command + + integer, intent(inout) :: arg_pos + + integer :: arg_len + character(len=100) :: textholder + arg_pos = arg_pos + 1 + + call get_command_argument(arg_pos, textholder, arg_len) + if(arg_len==0) stop "Missing argument to delete command" + + select case(textholder) + case('overlap') + arg_pos=arg_pos + 1 + call get_command_argument(arg_pos, textholder, arg_len) + if(arg_len==0) stop "Missing argument to delete overlap command" + read(textholder, *) rc_off + case default + print *, "Command ", trim(adjustl(textholder)), " is not accepted for option delete" + stop 3 + end select + + arg_pos = arg_pos + 1 + end subroutine parse_delete + + subroutine delete_overlap + !This subroutine deletes all overlapping atoms, which is defined as atoms which are separated by a distance of + !less then rc_off + + integer :: i, c(3), ci, cj, ck, num_nei, nei, delete_num + integer, dimension(atom_num) :: for_delete + + !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(:,:,:,:) + + allocate(which_cell(3,atom_num)) + + !First pass the atom list and atom num to the algorithm which builds the cell list + call build_cell_list(atom_num, r_atom, 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 + delete_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 + + delete_num = delete_num + 1 + for_delete(delete_num) = max(i,nei) + + !Now zero out the larger index + if(i > nei) then + which_cell(:,i) = 0 + cycle atom_loop + else + 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 *, "Overlap command deletes ", delete_num, " atoms" + !Now delete all the atoms + call delete_atoms(delete_num, for_delete(1:delete_num)) + end subroutine delete_overlap +end module opt_delete diff --git a/src/opt_disl.f90 b/src/opt_disl.f90 index 025d680..4281d12 100644 --- a/src/opt_disl.f90 +++ b/src/opt_disl.f90 @@ -47,8 +47,8 @@ module opt_disl integer, intent(inout) :: arg_pos integer :: i,arglen - character(len=8) :: ori_string character(len=100) :: textholder + character(len=8) :: ori_string !Parse all of the commands arg_pos = arg_pos + 1 @@ -180,6 +180,9 @@ module opt_disl end do end if + !Now make sure all atoms are wrapped back into the simulation cell + call wrap_atoms + end subroutine dislgen subroutine parse_disloop(arg_pos) @@ -187,8 +190,7 @@ module opt_disl integer, intent(inout) :: arg_pos - integer :: i,arglen, sbox - character(len=8) :: ori_string + integer :: i,arglen character(len=100) :: textholder !Parse all of the commands @@ -283,18 +285,18 @@ module opt_disl 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) + xLoop(1,a1) = centroid(a1) + loop_radius + xLoop(1,a2) = centroid(a2) + loop_radius + xLoop(1,a3) = centroid(a3) + xLoop(2,a1) = centroid(a1) - loop_radius + xLoop(2,a2) = centroid(a2) + loop_radius + xLoop(2,a3) = centroid(a3) + xLoop(3,a1) = centroid(a1) - loop_radius + xLoop(3,a2) = centroid(a2) - loop_radius + xLoop(3,a3) = centroid(a3) + xLoop(4,a1) = centroid(a1) + loop_radius + xLoop(4,a2) = centroid(a2) - loop_radius + xLoop(4,a3) = centroid(a3) else !Calculate loop perimeter perimeter = 2.0_dp*pi*loop_radius @@ -392,6 +394,9 @@ module opt_disl end do end do return + + !Now make sure all atoms are wrapped back into the simulation cell + call wrap_atoms end subroutine !******************************************************** diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 6d129f1..03b3063 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -8,10 +8,10 @@ module opt_group use box implicit none - integer :: group_ele_num, group_atom_num, remesh_size + integer :: group_ele_num, group_atom_num, remesh_size, remesh_type 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 + real(kind=dp) :: block_bd(6), disp_vec(3), remesh_lat_pam + logical :: displace, delete, max_remesh, refine integer, allocatable :: element_index(:), atom_index(:) @@ -27,6 +27,11 @@ module opt_group group_ele_num = 0 group_atom_num = 0 remesh_size=0 + displace=.false. + delete=.false. + max_remesh=.false. + refine = .false. + if(allocated(element_index)) deallocate(element_index) if(allocated(atom_index)) deallocate(atom_index) @@ -38,6 +43,11 @@ module opt_group if(displace) call displace_group if(remesh_size > 0) call remesh_group + + if(delete) call delete_group + + if(refine) call refine_group + end subroutine group subroutine parse_group(arg_pos) @@ -96,13 +106,25 @@ module opt_group if (arglen==0) stop "Missing vector component for shift command" read(textholder, *) disp_vec(i) end do - case('wrap') - wrap = .true. + case('refine') + refine=.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 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) stop "Missing remesh lattice parameter in group command" + read(textholder, *) remesh_lat_pam + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) stop "Missing remesh type in group command" + read(textholder, *) remesh_type + case('max') + max_remesh =.true. + case('delete') + delete=.true. case default !If it isn't an available option to opt_disl then we just exit exit @@ -192,16 +214,8 @@ module opt_group 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 + 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) @@ -213,47 +227,217 @@ module opt_group end subroutine displace_group - subroutine remesh_group + subroutine refine_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 + + 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), 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 + + 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), old_ele, old_atom, & + max_loops, working_esize + + 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) + logical, allocatable :: lat_points(:,:,:) + character(len=100) :: remesh_ele_type + + + !Right now we just hardcode only remeshing to elements + remesh_ele_type = 'fcc' + + !Get the orientations, this assumes that the orientation of the subbox for the first atom is the + !same as the rest of the atoms + !If this assumption is false then the code will break and exit + orient = sub_box_ori(:, :, sbox_atom(atom_index(1))) + call matrix_inverse(orient,3,ori_inv) + + !First calculate max position in lattice space to be able to allocate lat_points array, also sum the total numbers of + !degrees of freedom which are added + dof = 0 + + select case(trim(adjustl(shape))) + case('block') + group_in_lat = reshape((/ block_bd(1),block_bd(3),block_bd(5), & + block_bd(2),block_bd(3),block_bd(5), & + block_bd(2),block_bd(4),block_bd(5), & + block_bd(1),block_bd(4),block_bd(5), & + block_bd(1),block_bd(3),block_bd(6), & + block_bd(2),block_bd(3),block_bd(6), & + block_bd(2),block_bd(4),block_bd(6), & + block_bd(1),block_bd(4),block_bd(6) /), [3,8]) + + group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/remesh_lat_pam)) + 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 - !Once all atoms are added we delete all of the elements - call delete_elements(group_ele_num, element_index) + end select - print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms." + 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 - 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 + !Now place all group atoms and group interpolated atoms into lat_points + do i = 1, group_atom_num + r = r_atom(:,atom_index(i))/remesh_lat_pam + 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 + 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)/remesh_lat_pam + 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 + + print *, "Group has ", dof, " degrees of freedom to remesh" + + !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) + + old_atom = atom_num + old_ele = ele_num + + + + !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-2)/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(remesh_type) + 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))*remesh_lat_pam + 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. + call add_element(remesh_ele_type, working_esize, remesh_type, 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))*remesh_lat_pam + call add_atom(remesh_type, sbox_atom(atom_index(1)), r) + end if + end do + end do + end do + + print *, "Remeshing has created ", ele_num-old_ele, " elements and ", atom_num-old_atom, " atoms." 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 end module opt_group \ No newline at end of file diff --git a/src/opt_orient.f90 b/src/opt_orient.f90 new file mode 100644 index 0000000..c1bef1e --- /dev/null +++ b/src/opt_orient.f90 @@ -0,0 +1,132 @@ +module opt_orient + !This module contains the orient option which allows for the reorientation + !of simulation cells. This can be used to create arbitrarily oriented dislocation or loops. + use parameters + use subroutines + use elements + use box + implicit none + + + real(kind=dp), save :: new_orient(3,3) + real(kind=dp), dimension(6) :: orig_box_bd + real(kind=dp), allocatable :: orig_sub_box_ori(:,:,:) + + public + contains + + subroutine orient(arg_pos) + + integer, intent(inout) :: arg_pos + + integer :: i, ibasis, inod + logical :: isortho, isrighthanded + real(kind=dp) :: inv_sub_box_ori(3,3,sub_box_num) + + !First parse the orient command + call parse_orient(arg_pos) + + !Now rotate the basis. To do this we transform the basis to [100] [010] [001] and then + !transform to user specified basis. + + !Find all inverse orientation matrices for all sub_boxes + do i = 1, sub_box_num + call matrix_inverse(sub_box_ori, 3, inv_sub_box_ori) + end do + + !Now transform all atoms + do i = 1, atom_num + r_atom(:,i) = matmul(new_orient,matmul(inv_sub_box_ori(:,:,sbox_atom(i)),r_atom(:,i))) + end do + + !Now transform all elements + do i = 1, ele_num + do inod =1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + r_node(:,ibasis,inod,i) = matmul(new_orient,matmul(inv_sub_box_ori(:,:,sbox_ele(i)),r_node(:,ibasis,inod,i))) + end do + end do + end do + + !Now save the original sub_box_ori and overwrite them + if(allocated(orig_sub_box_ori)) deallocate(orig_sub_box_ori) + + allocate(orig_sub_box_ori(3,3,sub_box_num)) + orig_sub_box_ori = sub_box_ori + + !Now overwrite the orientations + do i = 1, sub_box_num + sub_box_ori(:,:,i) = new_orient + end do + + !Save original box boundaries + orig_box_bd = box_bd + + !Now find new box boundaries + call def_new_box + end subroutine orient + + subroutine parse_orient(arg_pos) + !This parses the orient option + integer, intent(inout) :: arg_pos + + integer :: i, arg_len + logical :: isortho, isrighthanded + character(len=8) :: ori_string + + !Pull out the new user orientation + do i = 1, 3 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, ori_string, arg_len) + if (arg_len == 0) print *, "Missing orientation vector in -orient option" + call parse_ori_vec(ori_string, new_orient(i,:)) + end do + + !Normalize the orientation matrix + new_orient = matrix_normal(new_orient,3) + + !Check right hand rule and orthogonality + call check_right_ortho(new_orient, isortho, isrighthanded) + if (.not.isortho) then + stop "Directions in orient are not orthogonal" + else if (.not.isrighthanded) then + stop "Directions in orient are not righthanded" + end if + + arg_pos = arg_pos + 1 + + end subroutine parse_orient + + subroutine unorient + + integer :: i, ibasis, inod + real(kind=dp) :: inv_ori(3,3) + + !Now rotate the basis. To do this we transform the basis to [100] [010] [001] and then + !transform to the original sbox_ele + + !Find the inverse for the new orientation matrix + call matrix_inverse(new_orient, 3, inv_ori) + + !Recover original sub_box_ori + sub_box_ori = orig_sub_box_ori + + !Now transform all atoms + do i = 1, atom_num + r_atom(:,i) = matmul(sub_box_ori(:,:,sbox_atom(i)),matmul(inv_ori(:,:),r_atom(:,i))) + end do + + !Now transform all elements + do i = 1, ele_num + do inod =1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + r_node(:,ibasis,inod,i) = matmul(sub_box_ori(:,:,sbox_ele(i)),matmul(inv_ori,r_node(:,ibasis,inod,i))) + end do + end do + end do + + !Restore original box boundaries + box_bd = orig_box_bd + end subroutine unorient + +end module opt_orient \ No newline at end of file diff --git a/src/parameters.f90 b/src/parameters.f90 index 6054140..2fa3e37 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -7,7 +7,7 @@ module parameters !Parameters for floating point tolerance real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), & lim_large = huge(1.0_dp) - logical, save :: lmpcac + logical, save :: lmpcac !Numeric constants real(kind=dp), parameter :: pi = 3.14159265358979323846_dp diff --git a/src/subroutines.f90 b/src/subroutines.f90 index 9fd159c..0be8d0e 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -5,6 +5,7 @@ module subroutines implicit none integer :: allostat, deallostat + public contains @@ -170,52 +171,7 @@ module subroutines return end subroutine parse_ori_vec - - subroutine parse_pos(i, pos_string, pos) - !This subroutine parses the pos command allowing for command which include inf - integer, intent(in) :: i !The dimension of the position - character(len=100), intent(in) :: pos_string !The position string - 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 - 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 - 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 - 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 - - if (iospara > 0) then - print *, "Error reading position argument ", trim(adjustl(pos_string)), ". Please reformat and try again." - end if - end subroutine parse_pos + subroutine heapsort(a) @@ -283,4 +239,107 @@ module subroutines end do end subroutine + subroutine build_cell_list(numinlist, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell) + !This subroutine builds a cell list based on rc_off + + !----------------------------------------Input/output variables------------------------------------------- + + integer, intent(in) :: numinlist !The number of points within r_list + + real(kind=dp), dimension(3,numinlist), intent(in) :: r_list !List of points to be used for the construction of + !the cell list. + + real(kind=dp), intent(in) :: rc_off ! Cutoff radius which dictates the size of the cells + + integer, dimension(3), intent(inout) :: cell_num !Number of cells in each dimension. + + integer, allocatable, intent(inout) :: num_in_cell(:,:,:) !Number of points within each cell + + integer, allocatable, intent(inout) :: cell_list(:,:,:,:) !Index of points from r_list within each cell. + + integer, dimension(3,numinlist), intent(out) :: which_cell !The cell index for each point in r_list + + !----------------------------------------Begin Subroutine ------------------------------------------- + + integer :: i, j, cell_lim, c(3) + real(kind=dp) :: box_len(3) + integer, allocatable :: resize_cell_list(:,:,:,:) + + !First calculate the number of cells that we need in each dimension + do i = 1,3 + box_len(i) = box_bd(2*i) - box_bd(2*i-1) + cell_num(i) = int(box_len(i)/(rc_off/2))+1 + end do + + !Initialize/allocate variables + cell_lim = 10 + allocate(num_in_cell(cell_num(1),cell_num(2),cell_num(3)), cell_list(cell_lim, cell_num(1), cell_num(2), cell_num(3))) + + !Now place points within cell + do i = 1, numinlist + !c is the position of the cell that the point belongs to + do j = 1, 3 + c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off/2)) + 1 + end do + + !Place the index in the correct position, growing if necessary + num_in_cell(c(1),c(2),c(3)) = num_in_cell(c(1),c(2),c(3)) + 1 + if (num_in_cell(c(1),c(2),c(3)) > cell_lim) then + allocate(resize_cell_list(cell_lim+10,cell_num(1),cell_num(2),cell_num(3))) + resize_cell_list(1:cell_lim, :, :, :) = cell_list + resize_cell_list(cell_lim+1:, :, :, :) = 0 + call move_alloc(resize_cell_list, cell_list) + end if + + cell_list(num_in_cell(c(1),c(2),c(3)),c(1),c(2),c(3)) = i + which_cell(:,i) = c + end do + + return + end subroutine build_cell_list + + + subroutine check_right_ortho(ori, isortho, isrighthanded) + !This subroutine checks whether provided orientations in the form: + ! | x1 x2 x3 | + ! | y1 y2 y3 | + ! | z1 z2 z3 | + !are right handed + + real(kind=dp), dimension(3,3), intent(in) :: ori + logical, intent(out) :: isortho, isrighthanded + + integer :: i, j + real(kind=dp) :: v(3), v_k(3) + !Initialize variables + isortho = .true. + isrighthanded=.true. + + do i = 1, 3 + do j = i+1, 3 + + if(abs(dot_product(ori(i,:), ori(j,:))) > lim_zero) then + isortho = .false. + end if + + !Check if they are righthanded + if (j == i+1) then + v(:) = cross_product(ori(i,:), ori(j,:)) + v_k(:) = v(:) - ori(mod(j, 3)+1,:) + else if ((i==1).and.(j==3)) then + v(:) = cross_product(ori(j,:),ori(i,:)) + v_k(:) = v(:) - ori(i+1, :) + end if + + if(norm2(v_k) > 10.0_dp**(-8.0_dp)) then + isrighthanded=.false. + end if + + end do + + end do + + return + end subroutine check_right_ortho + end module subroutines \ No newline at end of file