From 41024b96743df63463630027198d4846054b21c1 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Sun, 15 Dec 2019 18:59:21 -0500 Subject: [PATCH 1/6] Added new parse ori_vec subroutine to parse orientation vectors --- src/mode_create.f90 | 15 +-------------- src/subroutines.f90 | 25 +++++++++++++++++++++++++ 2 files changed, 26 insertions(+), 14 deletions(-) diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 23b9100..bac9b12 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -199,20 +199,7 @@ module mode_create call get_command_argument(arg_pos, orient_string, arglen) if (arglen==0) STOP "Missing orientation in orient command of mode create" arg_pos = arg_pos+1 - ori_pos=2 - do j = 1,3 - if (orient_string(ori_pos:ori_pos) == '-') then - ori_pos = ori_pos + 1 - read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j) - if (stat>0) STOP "Error reading orient value" - orient(i,j) = -orient(i,j) - ori_pos = ori_pos + 1 - else - read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j) - if(stat>0) STOP "Error reading orient value" - ori_pos=ori_pos + 1 - end if - end do + call parse_ori_vec(orient_string, orient(i,:)) end do diff --git a/src/subroutines.f90 b/src/subroutines.f90 index 0682037..f4447d7 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -145,4 +145,29 @@ module subroutines return end subroutine matrix_inverse + subroutine parse_ori_vec(ori_string, ori_vec) + !This subroutine parses a string to vector in the format [ijk] + character(len=8), intent(in) :: ori_string + real(kind=dp), dimension(3), intent(out) :: ori_vec + + integer :: i, ori_pos, stat + + ori_pos=2 + do i = 1,3 + if (ori_string(ori_pos:ori_pos) == '-') then + ori_pos = ori_pos + 1 + read(ori_string(ori_pos:ori_pos), *, iostat=stat) ori_vec(i) + if (stat>0) STOP "Error reading orientation value" + ori_vec(i) = -ori_vec(i) + ori_pos = ori_pos + 1 + else + read(ori_string(ori_pos:ori_pos), *, iostat=stat) ori_vec(i) + if(stat>0) STOP "Error reading orientation value" + ori_pos=ori_pos + 1 + end if + end do + + return + end subroutine parse_ori_vec + end module subroutines \ No newline at end of file From d0e6253d646928397d43fae8752e908d9008dea1 Mon Sep 17 00:00:00 2001 From: Alex Date: Wed, 25 Dec 2019 16:30:18 -0500 Subject: [PATCH 2/6] Working commit for dislocation generation branch --- README.md | 26 +++++++++++++++++++++++++- src/Makefile | 12 +++++++----- src/box.f90 | 12 +++++++++--- src/functions.f90 | 10 ++++++++++ src/io.f90 | 12 ++++++++---- src/main.f90 | 19 ++++++++++++++++--- src/mode_create.f90 | 2 +- src/parameters.f90 | 5 +++++ src/subroutines.f90 | 33 +++++++++++++++++++++++++++++++++ 9 files changed, 114 insertions(+), 17 deletions(-) diff --git a/README.md b/README.md index d390177..c8eb22b 100644 --- a/README.md +++ b/README.md @@ -111,4 +111,28 @@ This mode merges multiple data files and creates one big simulation cell. The pa `N` - The number of files which are being read -`dim` - the dimension they are to be stacked along, can be either `x`, `y`, or `z`. If the argument `none` is passed then the cells are just overlaid. Future options will include a delete overlap command. \ No newline at end of file +`dim` - the dimension they are to be stacked along, can be either `x`, `y`, or `z`. If the argument `none` is passed then the cells are just overlaid. Future options will include a delete overlap command. + +## Options + +Options are additional portions of code which have additional functionality. Options are performed in the order that they appear in the argument list and can be added to any mode. If wanting to use strictly options use `--convert` to specify input and output files. + +### Option dislgen + +``` +-dislgen [ijk] [hkl] x y z char_angle poisson +``` + + + +This options adds an arbitrarily oriented dislocation into your model based on user inputs using the volterra displacement fields. The options are below + +`[ijk]` - The vector for the line direction + +`[hkl]` - The vector for the slip plane + +`x y z` - The position of the dislocation centroid + +`char_angle` - Character angle of the dislocation (0 is screw and 90 is edge) + +`poisson` - Poisson's ratio used for the displacement field. \ No newline at end of file diff --git a/src/Makefile b/src/Makefile index c84186e..a55e152 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,8 +1,9 @@ FC=ifort #FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -FFLAGS=-mcmodel=large -Ofast +FFLAGS=-mcmodel=large -Ofast -no-wrap-margin MODES=mode_create.o mode_merge.o mode_convert.o -OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES) +OPTIONS=opt_disl.o +OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o call_option.o $(MODES) $(OPTIONS) .SUFFIXES: .SUFFIXES: .c .f .f90 .F90 .o @@ -25,9 +26,10 @@ cleantest: $(RM) testfuncs testfuncs.o $(OBJECTS) : parameters.o -atoms.o subroutines.o testfuncs.o : functions.o -main.o io.o build_subroutines.o: elements.o +atoms.o subroutines.o testfuncs.o box.o : functions.o +main.o io.o $(MODES) $(OPTIONS) : elements.o call_mode.o : $(MODES) +call_option.o : $(OPTIONS) $(MODES) io.o: atoms.o box.o $(MODES) main.o : io.o -testfuncs.o elements.o mode_create.o: subroutines.o +testfuncs.o elements.o mode_create.o opt_disl.o: subroutines.o diff --git a/src/box.f90 b/src/box.f90 index 3043851..de34e23 100644 --- a/src/box.f90 +++ b/src/box.f90 @@ -1,7 +1,7 @@ module box !This module contains information on the properties of the current box. use parameters - + use functions implicit none real(kind=dp) :: box_bd(6) !Global box boundaries @@ -34,8 +34,14 @@ module box integer, intent(in) :: n - allocate(sub_box_ori(3,3,n), sub_box_bd(6,n), sub_box_array_bd(2,2,n)) + integer :: i + allocate(sub_box_ori(3,3,n), sub_box_bd(6,n), sub_box_array_bd(2,2,n)) + do i = 1, n + sub_box_ori(:,:,i) = identity_mat(3) + sub_box_bd(:,i) = 0.0_dp + sub_box_array_bd(:,:,i) = 1 + end do end subroutine alloc_sub_box subroutine grow_sub_box(n) @@ -58,7 +64,7 @@ module box call move_alloc(temp_bd, sub_box_bd) temp_array_bd(:,:,1:sub_box_num) = sub_box_array_bd - temp_array_bd(:,:,sub_box_num+1:) = 0.0_dp + temp_array_bd(:,:,sub_box_num+1:) = 1 call move_alloc(temp_array_bd, sub_box_array_bd) return diff --git a/src/functions.f90 b/src/functions.f90 index eb8e813..f340b86 100644 --- a/src/functions.f90 +++ b/src/functions.f90 @@ -232,4 +232,14 @@ END FUNCTION StrDnCase end if return end function is_equal + + pure function unitvec(n,vec) + integer, intent(in) :: n + real(kind=dp), intent(in) :: vec(n) + real(kind=dp) :: unitvec(n) + + unitvec = vec/norm2(vec) + return + end function unitvec + end module functions diff --git a/src/io.f90 b/src/io.f90 index 6e914b2..46ece37 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -409,7 +409,7 @@ module io end do !Write the number of atom types in the current model and all of their names - write(11,*) atom_types, (type_to_name(i), i=1, atom_types) + write(11,*) atom_types, (type_to_name(i)//' ', i=1, atom_types) !Write the number of lattice_types, basisnum and number of nodes for each lattice type write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types) !Now for every lattice type write the basis atom types @@ -544,10 +544,14 @@ module io read(11,*) sub_box_bd(:,sub_box_num+i) sub_box_bd(:,sub_box_num+i) = sub_box_bd(:, sub_box_num+i) + displace(:) !Read in sub_box_array_bd - read(11,*) ((sub_box_ori(j, k, sub_box_num+i), j = 1, 2), k = 1, 2) - + read(11,*) ((sub_box_array_bd(j, k, sub_box_num+i), j = 1, 2), k = 1, 2) end do - sub_box_num = sub_box_num + n + + !Add the existing element boundaries + sub_box_array_bd(:,1,sub_box_num+1:) = sub_box_array_bd(:,1,sub_box_num+1:) + atom_num + sub_box_array_bd(:,2,sub_box_num+1:) = sub_box_array_bd(:,2,sub_box_num+1:) + ele_num + + sub_box_num = sub_box_num + n !Read in the number of atom types and all their names read(11, *) new_atom_types, (new_type_to_name(i), i = 1, new_atom_types) diff --git a/src/main.f90 b/src/main.f90 index a681d33..6bdc0ed 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -17,7 +17,7 @@ program main use io - integer :: i, end_mode_arg, arg_num + integer :: i, end_mode_arg, arg_num, arg_pos character(len=100) :: argument !Call initialization functions @@ -37,12 +37,25 @@ program main end if !Now we loop through all of the arguments and check for passed options or for a filename to be written out - do i = end_mode_arg-1, arg_num - call get_command_argument(i, argument) + arg_pos = end_mode_arg + do while(.true.) + !Exit the loop if we are done reading + if(arg_pos > arg_num) exit + + call get_command_argument(arg_pos, argument) !Check to see if a filename was passed if(scan(argument,'.',.true.) > 0) then call get_out_file(argument) + arg_pos = arg_pos + 1 + + !Check to see if an option has been passed + else if(argument(1:1) == '-') then + call call_option(argument, arg_pos) + !Otherwise print that the argument is not accepted and move on + else + print *, trim(adjustl(argument)), " is not accepted. Skipping to next argument" + arg_pos = arg_pos + 1 end if end do diff --git a/src/mode_create.f90 b/src/mode_create.f90 index bac9b12..ee35488 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -230,10 +230,10 @@ module mode_create case default !If it isn't an option then you have to exit + arg_pos = arg_pos -1 exit end select end do - !Calculate the lattice periodicity length in lattice units do i = 1, 3 lattice_space(i) = norm2(orient(i,:)) diff --git a/src/parameters.f90 b/src/parameters.f90 index 0443622..b54a702 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -2,7 +2,12 @@ module parameters implicit none + !Default precision integer, parameter :: dp= selected_real_kind(15,307) + !Parameters for floating point tolerance real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), & lim_large = huge(1.0_dp) + !Numeric constants + real(kind=dp), parameter :: pi = 3.14159265358979323846_dp + end module parameters \ No newline at end of file diff --git a/src/subroutines.f90 b/src/subroutines.f90 index f4447d7..d0b004b 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -1,6 +1,7 @@ module subroutines use parameters use functions + use box implicit none integer :: allostat, deallostat @@ -170,4 +171,36 @@ module subroutines return end subroutine parse_ori_vec + subroutine parse_pos(i, pos_string, pos) + !This subroutine parses the pos command allowing for command which include inf + integer, intent(in) :: i !The dimension of the position + character(len=100), intent(in) :: pos_string !The position string + real(kind=dp), intent(out) :: pos !The output parsed position value + + integer :: iospara + if(trim(adjustl(pos_string)) == 'inf') then + pos=box_bd(2*i) + else if(trim(adjustl(pos_string)) == '-inf') then + pos=box_bd(2*i-1) + else if ((index(pos_string,'-') > 0).and.(index(pos_string,'inf')>0)) then + !Now extract the number we are reducing from infinity + read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos + pos = box_bd(2*i) - pos + else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then + !Now extract the number we are reducing from infinity + read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos + pos = box_bd(2*i-1) + pos + else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then + !Now extract the number we are reducing from infinity + read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos + pos = (box_bd(2*i)-box_bd(2*i-1))*pos + else + read(pos_string, *, iostat=iospara) pos + end if + + if (iospara > 0) then + print *, "Error reading position argument ", trim(adjustl(pos_string)), ". Please reformat and try again." + end if + end subroutine parse_pos + end module subroutines \ No newline at end of file From 8fc91f4dd208f4993409ec27a69f68b1b66a5728 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 13 Jan 2020 09:42:57 -0500 Subject: [PATCH 3/6] Working planar vacancy cluster algorithm --- README.md | 16 +++++++++++++++- src/Makefile | 10 +++++++--- src/box.f90 | 17 +++++++++++++++++ src/elements.f90 | 23 +++++++++++++++++++++++ src/functions.f90 | 8 ++++++++ src/subroutines.f90 | 46 +++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 116 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index c8eb22b..b4b690d 100644 --- a/README.md +++ b/README.md @@ -135,4 +135,18 @@ This options adds an arbitrarily oriented dislocation into your model based on u `char_angle` - Character angle of the dislocation (0 is screw and 90 is edge) -`poisson` - Poisson's ratio used for the displacement field. \ No newline at end of file +`poisson` - Poisson's ratio used for the displacement field. + +### Option disloop + +```` +-disloop loop_normal loop_size x y z +```` + +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 + +`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 diff --git a/src/Makefile b/src/Makefile index a55e152..acc2066 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,6 +1,6 @@ FC=ifort -#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -FFLAGS=-mcmodel=large -Ofast -no-wrap-margin +FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin +#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin MODES=mode_create.o mode_merge.o mode_convert.o OPTIONS=opt_disl.o OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o call_option.o $(MODES) $(OPTIONS) @@ -19,12 +19,16 @@ clean: $(RM) cacmb *.o testfuncs: testfuncs.o functions.o subroutines.o - $(FC) testfuncs.o functions.o subroutines.o elements.o -o $@ + $(FC) testfuncs.o functions.o subroutines.o box.o elements.o -o $@ .PHONY: cleantest cleantest: $(RM) testfuncs testfuncs.o +.PHONY: test +test: testfuncs + ./testfuncs + $(OBJECTS) : parameters.o atoms.o subroutines.o testfuncs.o box.o : functions.o main.o io.o $(MODES) $(OPTIONS) : elements.o diff --git a/src/box.f90 b/src/box.f90 index de34e23..de8ef59 100644 --- a/src/box.f90 +++ b/src/box.f90 @@ -85,4 +85,21 @@ module box return end subroutine grow_box + + subroutine in_sub_box(r, which_sub_box) + !This returns which sub_box a point is in. It returns the first sub_box with boundaries which + !contain the point. + real(kind=dp), dimension(3), intent(in) :: r + integer, intent(out) :: which_sub_box + + integer :: i + + do i = 1, sub_box_num + if( in_block_bd(r, sub_box_bd(:,i))) then + which_sub_box = i + exit + end if + end do + return + end subroutine in_sub_box end module box \ No newline at end of file diff --git a/src/elements.f90 b/src/elements.f90 index 22e6eff..1c17643 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -3,6 +3,7 @@ module elements !This module contains the elements data structures, structures needed for building regions !and operations that are done on elements use parameters + use functions use subroutines implicit none @@ -355,4 +356,26 @@ module elements return end subroutine rhombshape + subroutine delete_atoms(num, index) + !This subroutine deletes atoms from the atom arrays + integer, intent(in) :: num + integer, intent(inout), dimension(num) :: index + + integer :: i, j + + call heapsort(index) + + !We go from largest index to smallest index just to make sure that we don't miss + !accidentally overwrite values which need to be deleted + do i = num, 1, -1 + if(index(i) == atom_num) then + r_atom(:,index(i)) = 0.0_dp + type_atom(index(i)) = 0 + else + r_atom(:,index(i)) = r_atom(:, atom_num) + type_atom(index(i)) = type_atom(atom_num) + end if + atom_num = atom_num - 1 + end do + end subroutine end module elements \ No newline at end of file diff --git a/src/functions.f90 b/src/functions.f90 index f340b86..82bcbcd 100644 --- a/src/functions.f90 +++ b/src/functions.f90 @@ -242,4 +242,12 @@ END FUNCTION StrDnCase return end function unitvec + pure function norm_dis(rl,rk) + !This just returns the magnitude of the vector between two points + real(kind=dp), dimension(3), intent(in) :: rl, rk + real(kind=dp) :: norm_dis(4) + + norm_dis(1:3) = (rk - rl) + norm_dis(4) = norm2(rk-rl) + end function end module functions diff --git a/src/subroutines.f90 b/src/subroutines.f90 index d0b004b..330fbca 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -203,4 +203,50 @@ module subroutines end if end subroutine parse_pos + subroutine heapsort(a) + + integer, intent(in out) :: a(0:) + integer :: start, n, bottom + real :: temp + + n = size(a) + do start = (n - 2) / 2, 0, -1 + call siftdown(a, start, n); + end do + + do bottom = n - 1, 1, -1 + temp = a(0) + a(0) = a(bottom) + a(bottom) = temp; + call siftdown(a, 0, bottom) + end do + + end subroutine heapsort + + subroutine siftdown(a, start, bottom) + + integer, intent(in out) :: a(0:) + integer, intent(in) :: start, bottom + integer :: child, root + real :: temp + + root = start + do while(root*2 + 1 < bottom) + child = root * 2 + 1 + + if (child + 1 < bottom) then + if (a(child) < a(child+1)) child = child + 1 + end if + + if (a(root) < a(child)) then + temp = a(child) + a(child) = a (root) + a(root) = temp + root = child + else + return + end if + end do + + end subroutine siftdown end module subroutines \ No newline at end of file From 1f1598dec9f4e4cb0178592b97b23f53914ee0e2 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 13 Jan 2020 09:49:42 -0500 Subject: [PATCH 4/6] Added Make install and fix to main --- src/Makefile | 4 ++++ src/main.f90 | 5 +++++ 2 files changed, 9 insertions(+) diff --git a/src/Makefile b/src/Makefile index acc2066..c20e5b2 100644 --- a/src/Makefile +++ b/src/Makefile @@ -29,6 +29,10 @@ cleantest: test: testfuncs ./testfuncs +.PHONY: install +install: cacmb + cp ./cacmb /usr/local/bin + $(OBJECTS) : parameters.o atoms.o subroutines.o testfuncs.o box.o : functions.o main.o io.o $(MODES) $(OPTIONS) : elements.o diff --git a/src/main.f90 b/src/main.f90 index 6bdc0ed..49b1a36 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -23,6 +23,7 @@ program main !Call initialization functions call lattice_init call box_init + end_mode_arg = 0 ! Command line parsing arg_num = command_argument_count() @@ -36,6 +37,10 @@ program main call call_mode(end_mode_arg, argument) end if + !Check to make sure a mode was called + if (end_mode_arg=0) then + stop "Please run cacmb using an available mode" + end if !Now we loop through all of the arguments and check for passed options or for a filename to be written out arg_pos = end_mode_arg do while(.true.) From ec7980ff0b2186adf813172f305ba4e5f3e1aa2d Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 13 Jan 2020 09:50:29 -0500 Subject: [PATCH 5/6] Fix to main when running without a mode --- src/main.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main.f90 b/src/main.f90 index 49b1a36..ed34c53 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -38,8 +38,8 @@ program main end if !Check to make sure a mode was called - if (end_mode_arg=0) then - stop "Please run cacmb using an available mode" + if (end_mode_arg==0) then + stop "Nothing to do, please run cacmb using an available mode" end if !Now we loop through all of the arguments and check for passed options or for a filename to be written out arg_pos = end_mode_arg From 20546dc2672435e6c0b8c73efb4bd60de47affe8 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 13 Jan 2020 16:43:49 -0500 Subject: [PATCH 6/6] Actual dislocation loop generation algorithm alongside dislgen command --- README.md | 10 +- src/io.f90 | 2 +- src/opt_disl.f90 | 520 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 528 insertions(+), 4 deletions(-) create mode 100644 src/opt_disl.f90 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