From 8fd2a7f65b84da14584cf82d4a1d4f29f2275f43 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Tue, 10 May 2022 09:46:46 -0400 Subject: [PATCH] Updated dislocation analysis and fix to alloy opt_group option --- src/mode_da.f90 | 108 ++++++++++++++++++++++++++++++++++++++-------- src/opt_group.f90 | 14 +++++- 2 files changed, 103 insertions(+), 19 deletions(-) diff --git a/src/mode_da.f90 b/src/mode_da.f90 index d290066..4798cab 100644 --- a/src/mode_da.f90 +++ b/src/mode_da.f90 @@ -8,7 +8,9 @@ module mode_da implicit none - integer, allocatable :: is_slipped(:,:,:) + integer :: num_d_points + integer, allocatable :: is_slipped(:,:,:) + real(kind=dp), allocatable :: disl_line(:,:), disl_data(:) private :: write_xyz public @@ -31,16 +33,28 @@ module mode_da if(allocated(is_slipped)) deallocate(is_slipped) allocate(is_slipped(max_basisnum, max_ng_node, ele_num)) + if(allocated(disl_line)) deallocate(disl_line) + num_d_points=0 + allocate(disl_line(3,1024)) + allocate(disl_data(1024)) + 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)//'vtk' + else + outfile = 'da_'//infiles(1)//'.vtk' + end if + call write_da_vtk(outfile) + 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 @@ -68,9 +82,10 @@ module mode_da integer :: i, j, k, l, ibasis, nei, close_neighbors(2,4,6), far_neighbor(4,6), minindex, minindices(2), linevec(3) integer :: face_types(ele_num*6), face_ele(6*ele_num), diff_pair(2,6) - real(kind = dp) :: face_centroids(3, ele_num*6), r(3), rc, vnode(3, max_basisnum, 4), vnorm(max_basisnum, 4), & + real(kind = dp) :: face_centroids(3, ele_num*6), r(3), rc, vnode(3, max_basisnum, 4), vnorm, & max_node_dist, ndiff, rmax(3), rmin(3), rnorm, v1(3), v2(3), v1norm, v2norm, theta1, theta2, & - diff_mag(6) + diff_mag(6), diffvec(3) + real(kind=dp), allocatable :: disl_line_temp(:,:), disl_data_temp(:) logical :: is_edge !Initialize variables @@ -123,8 +138,6 @@ module mode_da end do - !Calculate the largest difference between nodes for each face - end do !Now set the cut off distance as half the distance from opposite nodes in the largest element @@ -153,22 +166,24 @@ module mode_da 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 + !Now calculate the maximum difference between nodes + vnorm = 0 do j = 1, 4 - do ibasis = 1, basisnum(lat_ele(face_ele(i))) - vnorm(ibasis,j) = norm2(vnode(:, ibasis, j)) + do k=j,4 + v1=vnode(:,1,j) - vnode(:,1,k) + v1norm=norm2(v1) + if (vnorm < v1norm) then + vnorm = v1norm + diffvec= v1 + end if 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 > 1_dp)) then + if (vnorm>0.5_dp) then print *, "Element number ", face_ele(i), " is dislocated along face ", face_types(i), & - " with neighbor ", face_ele(nei), " with max displacement of ", maxval(vnorm) - + " with neighbor ", face_ele(nei), " with max displacement of ", vnorm l=0 do j = 1, 4 is_slipped(:, cubic_faces(j,face_types(i)), face_ele(i)) = 1 @@ -185,14 +200,15 @@ module mode_da end do end do minindex = minloc(diff_mag,1) - !Now figure out if the min differnce between nodes is associated with a 112 direction + !Now figure out if the min difference between nodes is associated with a 112 direction is_edge = .false. - do j = 1,6 if((diff_pair(1,minindex) == oneonetwopairs(1,j)).and.(diff_pair(2,minindex)==oneonetwopairs(2,j))) then is_edge=.true. + exit else if((diff_pair(2,minindex)==oneonetwopairs(1,j)).and.(diff_pair(1,minindex)==oneonetwopairs(2,j)))then is_edge=.true. + exit end if end do @@ -201,12 +217,70 @@ module mode_da else print *, 'Dislocation has primarily screw character' end if + + if( i < nei) then + num_d_points=num_d_points+2 + if(num_d_points > size(disl_line,2)) then + allocate(disl_line_temp(3,size(disl_line,2)+1024)) + disl_line_temp=0.0_dp + disl_line_temp(:, 1:size(disl_line,2))=disl_line + call move_alloc(disl_line_temp, disl_line) + end if + if(num_d_points/2 > size(disl_data)) then + allocate(disl_data_temp(size(disl_data)+1024)) + disl_data_temp=0 + disl_data_temp(1:size(disl_data)) = disl_data + call move_alloc(disl_data_temp, disl_data) + end if + print *, num_d_points/2, vnorm + disl_data(num_d_points/2)=vnorm + disl_line(:,num_d_points-1) =r_node(:, 1, diff_pair(1,minindex), face_ele(i)) + disl_line(:,num_d_points) =r_node(:, 1, diff_pair(2,minindex), face_ele(i)) + end if end if end if end do end subroutine cga_da + subroutine write_da_vtk(outfile) + character (len=*), intent(in) :: outfile + + integer :: i + open(unit=11, file=outfile, action='write', status='replace', position='rewind') + + write(11, '(a)') '# vtk DataFile Version 4.0.1' + write(11, '(a)') 'Dislocation line analyzed via cacmb dislocation analysis' + write(11, '(a)') 'ASCII' + write(11, '(a)') 'DATASET UNSTRUCTURED_GRID' + write(11, *) 'POINTS', num_d_points, 'float' + + !Write the points to the file + do i = 1, num_d_points + write(11, '(3f23.15)') disl_line(:, i) + end do + + write(11,*) 'CELLS', num_d_points/2, num_d_points/2*3 + do i=1, num_d_points, 2 + write(11,'(3i16)') 2, i-1, i + end do + + write(11,*) 'CELL_TYPES', num_d_points/2 + do i=1, num_d_points, 2 + write(11,'(i16)') 3 + end do + + write(11,*) 'CELL_DATA', num_d_points/2 + write(11, '(a)') 'SCALARS burgers_vec, float' + write(11, '(a)') 'LOOKUP_TABLE default' + do i=1, num_d_points, 2 + write(11, '(f23.15)') disl_data((i+1)/2) + end do + + close(11) + + end subroutine write_da_vtk + 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 diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 58b81bb..b8d6fe9 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -1187,17 +1187,27 @@ module opt_group subroutine alloy_group !This subroutine randomizes the atom types to reach desired concentrations, this only operates on atoms integer :: i, j, ia, type_map(10), added_types(num_species) - real(kind=dp) :: rand, mass + real(kind=dp) :: rand, mass, old_fractions(num_species) + character(len=2) :: sorted_species_type(10) print *, "Alloying group with desired fractions", s_fractions(1:num_species) + old_fractions=s_fractions(1:num_species) !First we sort the fractions call quicksort(s_fractions(1:num_species)) + + do i = 1,num_species + do j=1, num_species + if (is_equal(s_fractions(i), old_fractions(j))) then + sorted_species_type(i)=species_type(j) + end if + end do + end do !Now get the atom type maps for all the atoms and make the fractions a running sum do i = 1, num_species - call atommass(species_type(i), mass) + call atommass(sorted_species_type(i), mass) call add_atom_type(mass, type_map(i)) if(i > 1) s_fractions(i) = s_fractions(i) + s_fractions(i-1) end do