Added disl option

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)

src/io.f90

@@ -119,7 +119,7 @@ module io
         !This is the simplest visualization subroutine, it writes out all nodal positions and atom positions to an xyz file
         character(len=100), intent(in) :: file 
 
-        integer :: i, inod, ibasis
+        integer :: i, inod, ibasis, p_node, p_atom
 
         open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
 
@@ -128,21 +128,27 @@ module io
          
         !Write comment line
         write(11, '(a)') "#Node + atom file created using cacmb"
+        p_node=0 
+        p_atom = 0
 
         !Write nodal positions
         do i = 1, ele_num
             do inod = 1, ng_node(lat_ele(i))
                 do ibasis = 1, basisnum(lat_ele(i))
                         write(11, '(2i16, 3f23.15)') basis_type(ibasis,lat_ele(i)), 1, r_node(:,ibasis,inod,i)
+                        p_node = p_node + 1
                 end do 
             end do 
         end do 
 
+        if(p_node /= node_num) print *, "Error with node num"
         !Write atom positions
         do i = 1, atom_num
             write(11, '(2i16, 3f23.15)') type_atom(i), 0, r_atom(:,i) 
+            p_atom=p_atom +1 
         end do
 
+        if(p_atom /= atom_num) print *, "Error with atom num"
         !Finish writing 
         close(11)
     end subroutine write_xyz

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

src/opt_redef_box.f90

@@ -68,6 +68,7 @@ module opt_redef_box
 
         call delete_elements(delete_num, delete_list(1:delete_num))
 
+        print *, new_bd
         box_bd=new_bd
         box_bc = new_bc