diff --git a/src/call_option.f90 b/src/call_option.f90 index a195a0e..92dbe12 100644 --- a/src/call_option.f90 +++ b/src/call_option.f90 @@ -14,7 +14,7 @@ subroutine call_option(option, arg_pos) character(len=100), intent(in) :: option select case(trim(adjustl(option))) - case('-dislgen', '-disloop','-vacancydisloop') + case('-disl','-dislgen', '-disloop','-vacancydisloop') call dislocation(option, arg_pos) case('-group') call group(arg_pos) diff --git a/src/opt_disl.f90 b/src/opt_disl.f90 index dc0e72c..351b850 100644 --- a/src/opt_disl.f90 +++ b/src/opt_disl.f90 @@ -9,6 +9,7 @@ module opt_disl implicit none integer :: vloop_size ! Number of atoms to remove in planar vacancy cluster + integer :: iline, islip_plane real(kind=dp), dimension(3) :: line, slip_plane, centroid!dislocation line, slip plane vectors, centroid, real(kind=dp) :: burgers(3) !burgers vector of loop real(kind=dp) :: poisson, char_angle, lattice_parameter!Poisson ratio and character angle, lattice_parameter for burgers vector @@ -31,6 +32,9 @@ module opt_disl print *, '--------------------Option Dislocation-----------------------' select case(trim(adjustl(option))) + case('-disl') + call parse_disl(arg_pos) + call disl case('-dislgen') call parse_dislgen(arg_pos) call dislgen @@ -43,6 +47,148 @@ module opt_disl end select end subroutine dislocation + subroutine parse_disl(arg_pos) + + !Parse disl command + + integer, intent(inout) :: arg_pos + + integer :: i,arglen + character(len=100) :: textholder + character(len=1) :: parse_dim + + !Parse all of the commands + arg_pos = arg_pos + 1 + line(:) = 0.0_dp + + call get_command_argument(arg_pos, parse_dim, arglen) + if (arglen==0) STOP "Missing line dim in disl command" + select case(parse_dim) + case('x','X') + iline=1 + case('y','Y') + iline=2 + case('z','Z') + iline=3 + end select + + arg_pos=arg_pos + 1 + call get_command_argument(arg_pos, parse_dim, arglen) + if (arglen==0) STOP "Missing plane dim in disl command" + select case(parse_dim) + case('x','X') + islip_plane=1 + case('y','Y') + islip_plane=2 + case('z','Z') + islip_plane=3 + end select + do i = 1, 3 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing centroid in disl command" + call parse_pos(i, textholder, centroid(i)) + end do + + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing parameter in disl command" + read(textholder, *) b + + + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing character angle in disl command" + read(textholder, *) char_angle + + + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing poisson in disl command" + read(textholder, *) poisson + + arg_pos = arg_pos + 1 + end subroutine parse_disl + + subroutine disl + !This subroutine creates the actual dislocation + + integer :: i, sub_box, inod, ibasis, ix, iy, iz + 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 + + !Map box dimensions to dislocation dimensions, iz is along the dislocation line and iy is along the slip plane + iz = iline + iy = islip_plane + do i = 1, 3 + if ((i /= iy).and.(i /=iz)) then + ix = i + exit + end if + end do + + if(atom_num > 0) then + do i = 1, atom_num + r=r_atom(:,i) - centroid + if (r(ix) == 0) then + actan=pi/2 + else + actan = datan2(r(iy),r(ix)) + end if + + if ((r(ix)**2 + r(iy)**2) == 0) cycle + + !This is the elastic displacement field for dislocation according to Hirth and Lowe + disp(ix) = be/(2.0_dp*pi) * (actan + (r(ix)*r(iy))/(2.0_dp*(1.0_dp-poisson)*(r(ix)**2.0_dp + r(iy)**2.0_dp))) + disp(iy) = -be/(2.0_dp*pi)*((1.0_dp-2.0_dp*poisson)/(4.0_dp-4.0_dp*poisson) * & + log(r(ix)**2.0_dp + r(iy)**2.0_dp) & + + (r(ix)**2.0_dp - r(iy)**2.0_dp)/(4.0_dp*(1.0_dp-poisson)& + *(r(ix)**2.0_dp+r(iy)**2.0_dp))) + disp(3) = bs/(2.0_dp*pi) * actan + + r_atom(:,i) = r_atom(:,i) + disp + end do + end if + + if(ele_num > 0) then + 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)-centroid + if (r(ix) == 0) then + actan = pi/2 + else + actan = datan2(r(iy),r(ix)) + end if + + if ((r(ix)**2 + r(iy)**2) == 0) cycle + + !This is the elastic displacement field for dislocation according to Hirth and Lowe + disp(ix) = be/(2.0_dp*pi)*(actan + (r(ix)*r(iy))/(2.0_dp*(1.0_dp-poisson)*(r(ix)**2.0_dp + r(iy)**2.0_dp))) + disp(iy) = -be/(2.0_dp*pi)*((1.0_dp-2.0_dp*poisson)/(4.0_dp-4.0_dp*poisson) * & + log(r(ix)**2.0_dp + r(iy)**2.0_dp) & + + (r(ix)**2.0_dp - r(iy)**2.0_dp)/(4.0_dp*(1.0_dp-poisson)& + *(r(ix)**2.0_dp+r(iy)**2.0_dp))) + disp(3) = bs/(2.0_dp*pi) * actan + + + r_node(:,ibasis,inod,i) = r_node(:,ibasis,inod,i) + disp + end do + end do + end do + end if + + !Now make sure all atoms are wrapped back into the simulation cell + call wrap_atoms + + end subroutine disl + subroutine parse_dislgen(arg_pos) !Parse dislgen command diff --git a/src/opt_redef_box.f90 b/src/opt_redef_box.f90 index da86342..f6330e0 100644 --- a/src/opt_redef_box.f90 +++ b/src/opt_redef_box.f90 @@ -78,7 +78,7 @@ module opt_redef_box end subroutine redef_box - subroutine parse_redef_box(arg_pos) + subroutine parse_redef_box(arg_pos) !Parse the command integer, intent(inout) :: arg_pos