diff --git a/src/Makefile.dep b/src/Makefile.dep index 54ee834..e561193 100644 --- a/src/Makefile.dep +++ b/src/Makefile.dep @@ -11,6 +11,7 @@ obj/main : \ obj/mode_calc.o \ obj/mode_convert.o \ obj/mode_create.o \ + obj/mode_da.o \ obj/mode_merge.o \ obj/mode_metric.o \ obj/neighbors.o \ @@ -39,6 +40,7 @@ obj/caller.o : \ obj/mode_calc.o \ obj/mode_convert.o \ obj/mode_create.o \ + obj/mode_da.o \ obj/mode_merge.o \ obj/mode_metric.o \ obj/opt_deform.o \ @@ -94,6 +96,12 @@ obj/mode_create.o : \ obj/parameters.o \ obj/subroutines.o +obj/mode_da.o : \ + obj/elements.o \ + obj/io.o \ + obj/neighbors.o \ + obj/parameters.o + obj/mode_merge.o : \ obj/atoms.o \ obj/elements.o \ diff --git a/src/caller.f90 b/src/caller.f90 index 54c520e..18e91cf 100644 --- a/src/caller.f90 +++ b/src/caller.f90 @@ -1,6 +1,7 @@ module caller !this module just calls modes and options + use mode_da use mode_create use mode_convert use mode_merge @@ -37,6 +38,8 @@ module caller call metric(arg_pos) case('--calc') call calc(arg_pos) + case('--da') + call da(arg_pos) case default print *, "Mode ", trim(adjustl(mode)), " currently not accepted. Please check documentation for ", & "accepted modes and rerun." diff --git a/src/io.f90 b/src/io.f90 index 9798c51..d5dbab4 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -1121,7 +1121,7 @@ module io !Read header information read(11, *) textholder - !Read number of elements + !Read number of atoms read(11, *) atom_in, textholder read(11, *) type_in, textholder @@ -1164,7 +1164,7 @@ module io !Read atomic masses do i = 1, type_in - read(11,*) j, mass, textholder + read(11,*) j, mass call ATOMMASSSPECIES(mass, atom_species) call add_atom_type(atom_species, type_map(i)) end do diff --git a/src/main.f90 b/src/main.f90 index d4e7895..a41f396 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -117,7 +117,7 @@ program main !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 ((trim(adjustl(mode)) /= "--metric").and.(trim(adjustl(mode)) /= "--calc"))then + if ((trim(adjustl(mode)) /= "--metric").and.(trim(adjustl(mode)) /= "--calc").and.(trim(adjustl(mode)) /= "--da"))then if ((outfilenum == 0)) then argument = 'none' call get_out_file(argument) @@ -125,4 +125,5 @@ program main call write_out end if + return end program main diff --git a/src/mode_da.f90 b/src/mode_da.f90 new file mode 100644 index 0000000..82a9bc9 --- /dev/null +++ b/src/mode_da.f90 @@ -0,0 +1,171 @@ +module mode_da + !This mode is used to calculate the dislocation analysis algorithm for both CG and atomistic regions + + use parameters + use io + use elements + use neighbors + + implicit none + + integer, allocatable :: is_slipped(:,:,:) + + private :: write_xyz + public + contains + + subroutine da(arg_pos) + !This is the main calling subroutine for the dislocation analysis code + integer, intent(out) :: arg_pos + real(kind=dp), dimension(6) :: temp_box_bd + + integer :: ppos + character(len=100) :: outfile + + !Now call the dislocation analysis code for the coarse-grained elements + call parse_command(arg_pos) + + call read_in(1, (/0.0_dp,0.0_dp,0.0_dp/), temp_box_bd) + + !Now allocate the necessary variables + if(allocated(is_slipped)) deallocate(is_slipped) + allocate(is_slipped(max_basisnum, max_ng_node, ele_num)) + + call cga_da + + !Now create the output file num and write out to xyz format + ppos = scan(trim(infiles(1)),".", BACK= .true.) + if ( ppos > 0 ) then + outfile = 'da_'//infiles(1)(1:ppos)//'xyz' + else + outfile = 'da_'//infiles(1)//'.xyz' + end if + + call write_da_xyz(outfile) + + return + end subroutine da + + + subroutine parse_command(arg_pos) + !This subroutine parses the arguments for mode command + integer, intent(out) :: arg_pos + integer :: i, arglen + character(len = 100) :: textholder + logical :: file_exists + + !Read the input file + call get_command_argument(2,textholder, arglen) + if (arglen == 0) stop "Missing file for dislocation analysis" + call get_in_file(textholder) + + arg_pos = 3 + return + end subroutine parse_command + + subroutine cga_da + + integer :: i, j, k, l, ibasis, nei + + integer :: face_types(ele_num*6), face_ele(6*ele_num) + real(kind = dp) :: face_centroids(3, ele_num*6), r(3), rc, vnode(3, max_basisnum, 4), vnorm(max_basisnum, 4) + + !Initialize variables + l = 0 + + !First calculate all of the face centroids + do i = 1, ele_num + do j = 1, 6 + r(:) = 0.0_dp + do k = 1, 4 + do ibasis = 1, basisnum(lat_ele(i)) + r = r + r_node(:, ibasis, cubic_faces(k,j), i) + end do + end do + + r = r/(basisnum(lat_ele(i))*4) + + !add the face centroids, the type, and map the elements faces to the face arrays + l = l + 1 + face_centroids(:, l) = r + face_types(l) = j + face_ele(l) = i + end do + end do + + !Now calculate the nearest faces + rc = max_esize*maxval(lapa) + call calc_NN(l, face_centroids(:,1:l), rc) + + !Now loop overall the faces and make sure the nearest neighbor is the opposite type. If it isn't than we dscard it + is_slipped = 0 + do i = 1, l + nei = nn(i) + + !Skip if it's 0 + if (nei == 0) cycle + + !Check the face types, the way that the faces are laid out in the cubic_faces array face 1's opposite is 6 and + ! face 2's opposite is 5 and etc + vnode = 0 + if(face_types(i) == (7-face_types(nei))) then + vnorm = 0 + do j = 1, 4 + do ibasis = 1, basisnum(lat_ele(face_ele(i))) + !Compute the vectors between all nodes at the face. + vnode(:,ibasis,j) = r_node(:,ibasis, cubic_faces(j,face_types(i)), face_ele(i)) - & + r_node(:,ibasis, cubic_faces(j,face_types(nei)), face_ele(nei)) + end do + end do + + do j = 1, 3 + vnode(j,1:basisnum(lat_ele(face_ele(i))),:) = vnode(j,1:basisnum(lat_ele(face_ele(i))),:) - & + minval(vnode(j,1:basisnum(lat_ele(face_ele(i))),:)) + end do + do j = 1, 4 + do ibasis = 1, basisnum(lat_ele(face_ele(i))) + vnorm = norm2(vnode(:, ibasis, j)) + end do + end do + !Now calculate the difference between the largest norm and the smallest norm, if it's larger than 0.5 then mark it + !as slipped. This value probably can be converted to a variable value that depends on the current lattice parameter + !I think 0.5 works ok though. + if (any(vnorm > 0.5_dp)) then + do j = 1, 4 + is_slipped(:, cubic_faces(j,face_types(i)), face_ele(i)) = 1 + end do + end if + end if + end do + + end subroutine cga_da + + subroutine write_da_xyz(outfile) + !This subroutine write the element positions to a .xyz file and marks whether they are slipped or not + character(len=*), intent(in) :: outfile + integer :: i, ibasis, inod, outn + + open(unit=11, file=outfile, action='write', status='replace', position='rewind') + + !Write number of node_atoms + write(11, '(i16)') node_atoms + !Write comment + write(11, '(a)') "is_slipped x y z" + !Write nodal positions + outn = 0 + do i = 1, ele_num + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + write(11, '(1i16, 3f23.15)') is_slipped(ibasis,inod,i), r_node(:,ibasis,inod,i) + outn = outn + 1 + end do + end do + end do + + if(outn /= node_atoms) then + print *, "outn", outn, " doesn't equal node_atoms ", node_atoms + end if + + close(11) + end subroutine write_da_xyz +end module diff --git a/src/neighbors.f90 b/src/neighbors.f90 index 0faebca..f9d4d48 100644 --- a/src/neighbors.f90 +++ b/src/neighbors.f90 @@ -5,7 +5,7 @@ module neighbors use subroutines use functions - integer, allocatable :: nei_list(:,:), nei_num(:) + integer, allocatable :: nei_list(:,:), nei_num(:), nn(:) real(kind=dp), allocatable :: init_vec(:,:,:), output(:), microrotation(:,:) public contains @@ -139,4 +139,67 @@ module neighbors return end subroutine calc_neighbor + subroutine calc_NN(n,points, rc_off) + + integer, intent(in) :: n + real(kind=dp), intent(in) :: points(3,n) + real(kind=dp), intent(in) :: rc_off + + integer :: i, c(3), ci, cj, ck, nei + !cell arrays + integer, dimension(3) ::cell_num + integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:) + integer :: which_cell(3,n) + real(kind = dp) :: rmin + + !First reallocate the neighbor list codes + if (allocated(nn)) then + deallocate(nn) + end if + + allocate(nn(n)) + nn=0 + + !First build the cell lists + call build_cell_list(n, points, rc_off, cell_num, num_in_cell, cell_list, which_cell) + + pointloop: do i = 1, n + !First check to see if the point is a filler point, if so then skip it + if(points(1,i) < -Huge(-1.0_dp)+1) cycle + + !Get the positon of the cell + c = which_cell(:,i) + + !Initialize the min vec + rmin=Huge(1.0_dp) + + !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 + if((nei /= i)) then + !If it's the minimum position than we add it to the nearest neighbor list and updat e the min vec + + if (norm2(points(:,nei)-points(:,i)) < rmin) then + rmin = norm2(points(:, nei) - points(:,i)) + nn(i)=(nei) + end if + end if + end do + end do + end do + end do + end do pointloop + + print *, nn(1) + return + end subroutine calc_NN + end module neighbors