diff --git a/src/main.f90 b/src/main.f90 index d73270f..0f2d0f9 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -37,7 +37,7 @@ program main !Call initialization functions call lattice_init call box_init - call random_seed + call init_random_seed force_overwrite=.false. wrap_flag = .false. diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 603bf6e..c181e53 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -8,9 +8,12 @@ module opt_group use box implicit none - integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize - character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape - real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness + integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize, insert_type, & + insert_site + character(len=15) :: type, gshape!Type indicates what element type is selected and shape is the group shape + real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness, insert_conc, & + insert_lattice + logical :: displace, delete, max_remesh, refine, group_nodes, flip, efill, refinefill integer, allocatable :: element_index(:), atom_index(:) @@ -29,6 +32,7 @@ module opt_group group_ele_num = 0 group_atom_num = 0 remesh_size=0 + insert_type = 0 random_num=0 group_type=0 notsize=0 @@ -75,6 +79,11 @@ module opt_group call change_group_type end if + if(insert_type > 0) then + call get_group + call insert_group + end if + end subroutine group subroutine parse_group(arg_pos) @@ -98,11 +107,11 @@ module opt_group end select arg_pos = arg_pos + 1 - call get_command_argument(arg_pos, shape, arglen) + call get_command_argument(arg_pos, gshape, 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))) + select case(trim(adjustl(gshape))) case('block') do i= 1, 6 arg_pos = arg_pos + 1 @@ -146,7 +155,7 @@ module opt_group vertices(i,2) = box_bd(2*i-1) vertices(i,3) = box_bd(2*i-1) else - print *, "bwidth cannot be 0 in wedge shaped group" + print *, "bwidth cannot be 0 in wedge gshaped group" stop 3 end if else if (i == dim2) then @@ -369,7 +378,7 @@ module opt_group continue case default - print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", & + print *, "Group shape ", trim(adjustl(gshape)), " is not currently accepted. Please check documentation ", & "for accepted group shapes." end select @@ -426,8 +435,34 @@ module opt_group if(arglen ==0) stop "Missing notsize size" read(textholder, *) notsize print *, "Ignoring elements with size ", notsize + case('insert') + arg_pos=arg_pos+1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) stop "Missing element type for insert command" + call add_atom_type(textholder, insert_type) + arg_pos=arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + select case(trim(adjustl(textholder))) + case('tetra') + insert_site=1 + case('octa') + insert_site=2 + case('mixed') + insert_site=3 + case default + print *, "site value must be tetra, octa, or mixed and not ", trim(adjustl(textholder)) + stop 3 + end select + arg_pos=arg_pos+1 + call get_command_argument(arg_pos, textholder, arglen) + if(arglen ==0) stop "Missing lattice_type in insert command" + read(textholder, *) insert_lattice + arg_pos=arg_pos+1 + call get_command_argument(arg_pos, textholder, arglen) + if(arglen ==0) stop "Missing concentration in insert command" + read(textholder, *) insert_conc case default - !If it isn't an available option to opt_disl then we just exit + !If it isn't an available option to opt_group then we just exit exit end select end do @@ -440,7 +475,7 @@ module opt_group integer, allocatable :: resize_array(:) real(kind=dp) :: r_center(3), rand - select case(trim(adjustl(shape))) + select case(trim(adjustl(gshape))) case('block') print *, "Group has block shape with boundaries: ", block_bd case ('wedge') @@ -1007,10 +1042,71 @@ module opt_group end subroutine change_group_type + subroutine insert_group + !This code inserts atoms into interstitial sites. This only works on atoms within the group due to the limitations with the + !Coarse-graining methodology which doesn't allow for concentration fluctuations. + real(kind=dp) interstitial_sites(3,14), rand, rinsert(3) + integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi, sbox + integer, allocatable :: used_sites(:,:) + + !First save all of the displacement vectors from a lattice site to interstitial site + !The first 6 are the octohedral sites and the next 8 are the tetrahedral sites + interstitial_sites= reshape( (/ -0.5_dp, 0.0_dp, 0.0_dp, & + 0.5_dp, 0.0_dp, 0.0_dp, & + 0.0_dp,-0.5_dp, 0.0_dp, & + 0.0_dp, 0.5_dp, 0.0_dp, & + 0.0_dp, 0.0_dp,-0.5_dp, & + 0.0_dp, 0.0_dp, 0.5_dp, & + -0.25_dp,-0.25_dp,-0.25_dp, & + -0.25_dp,-0.25_dp, 0.25_dp, & + -0.25_dp, 0.25_dp,-0.25_dp, & + -0.25_dp, 0.25_dp, 0.25_dp, & + 0.25_dp,-0.25_dp,-0.25_dp, & + 0.25_dp,-0.25_dp, 0.25_dp, & + 0.25_dp, 0.25_dp,-0.25_dp, & + 0.25_dp, 0.25_dp, 0.25_dp /), & + shape(interstitial_sites)) + + !First we calculate the number of atoms needed based on the number of atoms in the group and the concentration + interstitial_sites=interstitial_sites*insert_lattice + + add_num = (insert_conc*group_atom_num)/(1-insert_conc) + allocate(used_sites(2,add_num)) + + print *, "Inserting ", add_num, " atoms as atom type ", insert_type + + !Now set up the random number generator for the desired interstitial type + select case(insert_site) + case(1) + rlo=1 + rhi=6 + case(2) + rlo=7 + rhi = 14 + case(3) + rlo=1 + rhi=14 + end select - subroutine split_group_elements - ! - end subroutine split_group_elements + !Now add the atoms + i = 1 + addloop:do while ( i < add_num) + call random_number(rand) + rindex = int(1+rand*(group_atom_num-1)) + ia=atom_index(rindex) + call random_number(rand) + sindex = int(rlo+rand*(rhi-rlo)) + do j = 1, i + if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop + end do + rinsert = r_atom(:,ia) + matmul(sub_box_ori(:,:,sbox_atom(ia)),interstitial_sites(:,sindex)) + sbox = sbox_atom(ia) + call add_atom(0, insert_type, sbox, rinsert) + used_sites(1,i) = ia + used_sites(2,i) = sindex + i = i + 1 + end do addloop + end subroutine insert_group function in_group(r) !This subroutine determines if a point is within the group boundaries @@ -1019,7 +1115,7 @@ module opt_group integer :: dim3, i logical :: in_group - select case(trim(adjustl(shape))) + select case(trim(adjustl(gshape))) case('block') in_group=in_block_bd(r,block_bd) case('wedge') @@ -1068,7 +1164,7 @@ module opt_group in_group_ele=.false. - if(trim(adjustl(shape)) == 'shell') then + if(trim(adjustl(gshape)) == 'shell') then node_in_out(:) = -1 !First calculate whether each element node is within the shell region, inside the shell sphere, or outside the !shell region @@ -1083,13 +1179,13 @@ module opt_group exit nodeloop end if - shape ='sphere' + gshape ='sphere' if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then node_in_out(inod) = 1 else node_in_out(inod) = 0 end if - shape='shell' + gshape='shell' end do nodeloop !If any nodes are within the shell region, or if the shell region interescts an element then add it to the group diff --git a/src/subroutines.f90 b/src/subroutines.f90 index b272c16..5cf147f 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -241,5 +241,52 @@ module subroutines return end subroutine check_right_ortho + subroutine init_random_seed() + implicit none + integer, allocatable :: seed(:) + integer :: i, n, un, istat, dt(8), pid, t(2), s + integer(8) :: count, tms + + call random_seed(size = n) + allocate(seed(n)) + ! First try if the OS provides a random number generator + open(newunit=un, file="/dev/urandom", access="stream", & + form="unformatted", action="read", status="old", iostat=istat) + if (istat == 0) then + read(un) seed + close(un) + else + ! Fallback to XOR:ing the current time and pid. The PID is + ! useful in case one launches multiple instances of the same + ! program in parallel. + call system_clock(count) + if (count /= 0) then + t = transfer(count, t) + else + call date_and_time(values=dt) + tms = (dt(1) - 1970) * 365_8 * 24 * 60 * 60 * 1000 & + + dt(2) * 31_8 * 24 * 60 * 60 * 1000 & + + dt(3) * 24 * 60 * 60 * 60 * 1000 & + + dt(5) * 60 * 60 * 1000 & + + dt(6) * 60 * 1000 + dt(7) * 1000 & + + dt(8) + t = transfer(tms, t) + end if + s = ieor(t(1), t(2)) + pid = getpid() + 1099279 ! Add a prime + s = ieor(s, pid) + if (n >= 3) then + seed(1) = t(1) + 36269 + seed(2) = t(2) + 72551 + seed(3) = pid + if (n > 3) then + seed(4:) = s + 37 * (/ (i, i = 0, n - 4) /) + end if + else + seed = s + 37 * (/ (i, i = 0, n - 1 ) /) + end if + end if + call random_seed(put=seed) + end subroutine init_random_seed end module subroutines