diff --git a/README.md b/README.md index b4b690d..cf7fcfb 100644 --- a/README.md +++ b/README.md @@ -140,13 +140,17 @@ This options adds an arbitrarily oriented dislocation into your model based on u ### Option disloop ```` --disloop loop_normal loop_size x y z +-disloop loop_normal radius x y z bx by bz poisson ```` This option deletes vacancies on a plane which when minimized should result in a dislocation loop structure. The arguments are below: -`dim` - The box dimension which defines the normal to the loop plane. As of now this dimension must be a closed back direction, meaning that for fcc a box dimension has to be of the (111) family of planes +`dim` - The box dimension which defines the normal to the loop plane. As of now this dimension must be a closed back direction, meaning that for fcc a box dimension has to be of the (111) family of planes. Either `x`, `y`, or `z`. `n` - The number of atoms to delete on the loop plane -`x y z` - The centroid of the loop. As a note, \ No newline at end of file +`x y z` - The centroid of the loop. + +`bx by bz` - The burgers vector for the dislocation + +`poisson` - Poisson ratio for continuum solution \ No newline at end of file diff --git a/src/io.f90 b/src/io.f90 index 46ece37..ba26062 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -315,7 +315,7 @@ module io write(11,3) node_num write(11,19) max_ng_node write(11,4) lattice_types - write(11,2) atom_num + write(11,5) atom_num write(11,6) 1 !Grain_num is ignored write(11,16) lattice_types, atom_types write(11,21) diff --git a/src/opt_disl.f90 b/src/opt_disl.f90 new file mode 100644 index 0000000..4510ede --- /dev/null +++ b/src/opt_disl.f90 @@ -0,0 +1,520 @@ +module opt_disl + + !This module contains all code associated with dislocations + + use parameters + use elements + use subroutines + use box + implicit none + + + 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 + character(len=10) :: lattice + character(len=1) :: loop_normal + real(kind=dp) :: loop_radius !Dislocation loop radius + + real(kind=dp) :: b !burgers vector + + + public + 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 + + + select case(trim(adjustl(option))) + case('-dislgen') + call parse_dislgen(arg_pos) + call dislgen + case('-disloop') + call parse_disloop(arg_pos) + call disloop + end select + end subroutine dislocation + + subroutine parse_dislgen(arg_pos) + + !Parse dislgen command + + integer, intent(inout) :: arg_pos + + integer :: i,arglen + character(len=8) :: ori_string + character(len=100) :: textholder + + !Parse all of the commands + arg_pos = arg_pos + 1 + line(:) = 0.0_dp + call get_command_argument(arg_pos, ori_string, arglen) + if (arglen==0) STOP "Missing line vector in dislgen command" + call parse_ori_vec(ori_string, line) + + arg_pos=arg_pos + 1 + slip_plane(:) = 0.0_dp + call get_command_argument(arg_pos, ori_string, arglen) + if (arglen==0) STOP "Missing plane vector in dislgen command" + call parse_ori_vec(ori_string, slip_plane) + + do i = 1, 3 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing centroid in dislgen command" + 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 character angle in dislgen 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 dislgen command" + 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 + !This subroutine creates the actual dislocation + + integer :: i, sub_box, inod, ibasis + 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) + + !Calculate screw and edge burgers vectors + be = sin(char_angle*pi/180.0_dp)*b + bs = cos(char_angle*pi/180.0_dp)*b + + !Figure out which sub box you are in so you can access it's orientation, this code will not work + !with overlapping sub_boxes + do i = 1, sub_box_num + if(in_block_bd(centroid, sub_box_bd(:,i))) then + sub_box = i + exit + end if + end do + + !Construct the slip system orientation matrix in an unrotated system + slipx = cross_product(slip_plane, line) + ss_ori(1,:) = unitvec(3,slipx) + ss_ori(2,:) = unitvec(3,slip_plane) + ss_ori(3,:) = unitvec(3,line) + call matrix_inverse(ss_ori, 3, ss_inv) + + !Apply the rotation + disp_transform = matmul(sub_box_ori(:,:,i), ss_inv) + call matrix_inverse(disp_transform,3,inv_transform) + + if(atom_num > 0) then + do i = 1, atom_num + r=r_atom(:,i) - centroid + r=matmul(inv_transform, r) + if (r(1) == 0) then + actan=pi/2 + else + actan = datan2(r(2),r(1)) + end if + + if ((r(1)**2 + r(2)**2) == 0) cycle + + !This is the elastic displacement field for dislocation according to Hirth and Lowe + disp(1) = be/(2.0_dp*pi) * (actan + (r(1)*r(2))/(2.0_dp*(1.0_dp-poisson)*(r(1)**2.0_dp + r(2)**2.0_dp))) + disp(2) = -be/(2.0_dp*pi)*((1.0_dp-2.0_dp*poisson)/(4.0_dp-4.0_dp*poisson) * & + log(r(1)**2.0_dp + r(2)**2.0_dp) & + + (r(1)**2.0_dp - r(2)**2.0_dp)/(4.0_dp*(1.0_dp-poisson)& + *(r(1)**2.0_dp+r(2)**2.0_dp))) + disp(3) = bs/(2.0_dp*pi) * actan + + !This transforms the displacement to the correct orientation + disp = matmul(disp_transform, disp) + + 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) + r = matmul(inv_transform, r) + if (r(1) == 0) then + actan = pi/2 + else + actan = datan2(r(2),r(1)) + end if + + if ((r(1)**2 + r(2)**2) == 0) cycle + + !This is the elastic displacement field for dislocation according to Hirth and Lowe + disp(1) = be/(2.0_dp*pi) * (actan + (r(1)*r(2))/(2.0_dp*(1.0_dp-poisson)*(r(1)**2.0_dp + r(2)**2.0_dp))) + disp(2) = -be/(2.0_dp*pi)*((1.0_dp-2.0_dp*poisson)/(4.0_dp-4.0_dp*poisson) * & + log(r(1)**2.0_dp + r(2)**2.0_dp) & + + (r(1)**2.0_dp - r(2)**2.0_dp)/(4.0_dp*(1.0_dp-poisson)& + *(r(1)**2.0_dp+r(2)**2.0_dp))) + disp(3) = bs/(2.0_dp*pi) * actan + + + disp = matmul(disp_transform, disp) + r_node(:,ibasis,inod,i) = r_node(:,ibasis,inod,i) + disp + end do + end do + end do + end if + + end subroutine dislgen + + subroutine parse_disloop(arg_pos) + !This subroutine parses the disloop command + + integer, intent(inout) :: arg_pos + + integer :: i,arglen, sbox + character(len=8) :: ori_string + character(len=100) :: textholder + + !Parse all of the commands + arg_pos = arg_pos + 1 + loop_normal = ' ' + call get_command_argument(arg_pos, loop_normal, arglen) + if (arglen==0) STOP "Missing loop_normal in disloop command" + !Convert the loop_normal to the dimension + select case(loop_normal) + case('x','X', 'y', 'Y', 'z', 'Z') + continue + case default + print *, "Dimension argument must either be x, y, or z not", loop_normal + stop 3 + end select + + arg_pos = arg_pos + 1 + loop_radius = 0 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing loop_size in disloop command" + read(textholder, *) loop_radius + + do i = 1, 3 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing centroid in disloop command" + call parse_pos(i, textholder, centroid(i)) + end do + + burgers(:) = 0.0_dp + do i = 1, 3 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing burgers vector in disloop command" + read(textholder, *) burgers(i) + end do + + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, textholder, arglen) + if (arglen==0) STOP "Missing poisson ratio in disloop command" + read(textholder, *) poisson + + arg_pos = arg_pos + 1 + + !Now check to make sure that the dimension selected is actually a 1 1 1 direction. + ! call in_sub_box(centroid, sbox) + ! if(.not.((abs(sub_box_ori(loop_normal,1,sbox)) == abs(sub_box_ori(loop_normal,2,sbox))).and. & + ! (abs(sub_box_ori(loop_normal,2,sbox)) == abs(sub_box_ori(loop_normal,3,sbox))).and. & + ! (abs(sub_box_ori(loop_normal,3,sbox)) == abs(sub_box_ori(loop_normal,1,sbox))))) then + + ! print *, "The selected dimension ", loop_normal, " for sub_box ", sbox, " is ", & + ! sub_box_ori(loop_normal,:,sbox), " which is not in the (111) family of planes" + ! STOP 3 + ! end if + + + end subroutine + + !Code for the creation of dislocation loops is based on functions from atomsk + !which is available from https://atomsk.univ-lille.fr/. They have been adapted to fit cacmb. + subroutine disloop + !This subroutine applies the displacement field for a dislocation loop to all atoms and elements. + integer :: a1, a2, a3 !New directions + integer :: i, j, inod, ibasis, Npoints + 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 + + if(allocated(xLoop)) deallocate(xLoop) + + !Define new directions + select case(loop_normal) + case('x','X') + !Loop in (y,z) plane + a1 = 2 + a2 = 3 + a3 = 1 + case('y','Y') + !Loop in (x,z) plane + a1 = 3 + a2 = 1 + a3 = 2 + case default + !Loop in (x,y) plane + a1 = 1 + a2 = 2 + 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 + + !Now actually calculate the displacement created by a loop for every atom + do i = 1, atom_num + u = 0.0_dp + omega = 0.0_dp + xC(:) = centroid - r_atom(:,i) + + !Loop over all dislocation segments + do j = 1, size(xLoop,1) + !Coordinates of point A + if (j ==1) then + xA(:) = xLoop(size(xLoop,1),:) - r_atom(:,i) + else + xA(:) = xLoop(j-1, :) - r_atom(:,i) + end if + + !Coordinates of point B + xB(:) = xLoop(j,:) - r_atom(:,i) + + !Displacement due to solid angle + omega = omega + SolidAngle(xA, xB, xC) + !Displacement due to elasticity + u(:) = u(:) + DisloSeg_displacement_iso(xA, xB, burgers(:), poisson) + end do + + !Total displacement + u=u(:) + burgers(:) * omega + + !Apply displacement to atom + r_atom(:,i) = r_atom(:,i) + u(:) + end do + + !Repeat for element nodes + do i = 1, ele_num + do inod =1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + u = 0.0_dp + omega = 0.0_dp + xC(:) = centroid - r_node(:,ibasis,inod,i) + + !Loop over all dislocation segments + do j = 1, size(xLoop,1) + !Coordinates of point A + if (j ==1) then + xA(:) = xLoop(size(xLoop,1),:) - r_node(:,ibasis,inod,i) + else + xA(:) = xLoop(j-1, :) - r_node(:,ibasis,inod,i) + end if + + !Coordinates of point B + xB(:) = xLoop(j,:) - r_node(:,ibasis,inod,i) + + !Displacement due to solid angle + omega = omega + SolidAngle(xA, xB, xC) + !Displacement due to elasticity + u(:) = u(:) + DisloSeg_displacement_iso(xA, xB, burgers(:), poisson) + end do + + !Total displacement + u=u(:) + burgers(:) * omega + + !Apply displacement to atom + r_node(:,ibasis,inod,i) = r_node(:,ibasis,inod,i) + u(:) + end do + end do + end do + return + end subroutine + + !******************************************************** + ! SOLIDANGLE + ! Calculate solid angle (normalized by 4*pi) used to obtain the + ! displacement created by a dislocation triangle loop ABC at the origin + ! (field point = origin) + ! Ref.: Van Oosterom, A. and Strackee, J., + ! The Solid Angle of a Plane Triangle, + ! IEEE Transactions on Biomedical Engineering BME-30, 125 (1983). + !******************************************************** + FUNCTION SolidAngle(xA, xB, xC) RESULT(Omega) + ! + IMPLICIT NONE + ! + ! Extremities of the triangle loop + REAL(dp),DIMENSION(3),INTENT(IN) :: xA, xB, xC + ! + ! Solid angle (normalized by 4*pi) + REAL(dp):: omega + REAL(dp),PARAMETER:: factor=1.d0/(2.d0*pi) + ! + REAL(dp) :: rA, rB, rC, numerator, denominator + ! + rA = norm2(xA) + rB = norm2(xB) + rC = norm2(xC) + ! + numerator = TRIPLE_PRODUCT( xA, xB, xC ) + denominator = rA*rB*rC + DOT_PRODUCT( xA, xB )*rC & + & + DOT_PRODUCT( xB, xC )*rA + DOT_PRODUCT( xC, xA )*rB + ! + omega = factor*ATAN2( numerator, denominator ) + ! + END FUNCTION SolidAngle + + !******************************************************** + ! DISLOSEG_DISPLACEMENT_ISO + ! Calculate displacement created by dislocation segment AB + ! once the solid angle part has been removed + ! Isotropic elastic calculation with nu Poisson coef. + ! Ref.: Eq. (1) in Barnett, D. M. + ! The Displacement Field of a Triangular Dislocation Loop + ! Philos. Mag. A, 1985, 51, 383-387 + !******************************************************** + FUNCTION DisloSeg_displacement_iso(xA, xB, b, nu) RESULT(u) + ! + IMPLICIT NONE + REAL(dp),DIMENSION(3),INTENT(IN):: xA, xB ! Extremities of the segment + REAL(dp),DIMENSION(3),INTENT(IN):: b ! Burgers vector + REAL(dp),INTENT(IN):: nu ! Poisson coefficient + REAL(dp),DIMENSION(3):: u ! Displacement + REAL(dp):: rA, rB + REAL(dp),DIMENSION(1:3):: tAB, nAB + ! + rA = norm2(xA) + rB = norm2(xB) + ! + ! Tangent vector + tAB(:) = xB(:) - xA(:) + tAB(:) = tAB(:)/norm2(tAB) + ! + ! Normal vector + nAB(:) = CROSS_PRODUCT(xA,xB) + nAB(:) = nAB(:)/norm2(nAB) + ! + u(:) = ( -(1.d0-2.d0*nu)*CROSS_PRODUCT(b, tAB)* & + & LOG( (rB + DOT_PRODUCT(xB,tAB))/(rA + DOT_PRODUCT(xA,tAB)) ) & + & + DOT_PRODUCT(b,nAB)*CROSS_PRODUCT(xB/rB-xA/rA,nAB) ) & + & /(8.d0*pi*(1.d0-nu)) + ! + END FUNCTION DisloSeg_displacement_iso +! + !This code simply creates a planar vacancy cluster and does not apply the dislocation loop displacement field. + ! subroutine disloop + ! !This subroutine actually creates the dislocation loop. + + ! real(kind=dp) :: neighbor_dis(loop_size), temp_box(6), dis + ! integer :: i, j, index(loop_size) + + ! neighbor_dis(:) = HUGE(1.0_dp) + ! index(:) = 0 + + ! !First find the nearest atom to the centroid + ! do i = 1, atom_num + ! if(norm2(r_atom(:,i) - centroid) < neighbor_dis(1)) then + ! neighbor_dis(1) = norm2(r_atom(:,i) - centroid) + ! index(1) = i + ! end if + ! end do + + ! !Define a new box, this box tries to isolate all atoms on the plane of the atom + ! !closest to the user defined centroid. + ! temp_box(:) = box_bd(:) + ! temp_box(2*loop_normal) = r_atom(loop_normal,index(1)) + 10.0_dp**(-2.0_dp) + ! temp_box(2*loop_normal-1) = r_atom(loop_normal,index(1)) - 10.0_dp**(-2.0_dp) + + ! !Now reset the list for the scanning algorithm + ! index(1) = 0 + ! neighbor_dis(1) = HUGE(1.0_dp) + + ! !Now scan over all atoms again and find the closest loop_size number of atoms to the initial atom + ! !that reside on the same plane. + + ! do i = 1, atom_num + ! !Check to see if it is on the same plane + ! if (in_block_bd(r_atom(:,i), temp_box)) then + ! dis = norm2(r_atom(:,i) - centroid) + ! do j = 1, loop_size + ! !Check to see if it is closer than other atoms + ! if (dis < neighbor_dis(j)) then + ! !Move values in the neighbor array and add the new neighbor + ! if(j < loop_size) then + ! neighbor_dis(j+1:loop_size) = neighbor_dis(j:loop_size-1) + ! index(j+1:loop_size) = index(j:loop_size-1) + ! end if + ! neighbor_dis(j) = dis + ! index(j) = i + ! exit + ! end if + ! end do + ! end if + ! end do + + ! !Now delete the atoms + ! call delete_atoms(loop_size, index) + + ! return + ! end subroutine disloop + +end module opt_disl \ No newline at end of file