From 44efb4be4a4ce38c09fe4d563661c9987b9ad88a Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 3 Feb 2020 17:30:41 -0500 Subject: [PATCH 1/2] Added orient and unorient options --- src/Makefile | 4 ++-- src/call_option.f90 | 9 ++++++--- src/elements.f90 | 48 ++++++++++++++++++++++++++++++++++++++++++--- src/io.f90 | 12 ++++++------ src/mode_create.f90 | 32 ++++++++++++++++-------------- src/opt_group.f90 | 2 +- 6 files changed, 77 insertions(+), 30 deletions(-) diff --git a/src/Makefile b/src/Makefile index a3f926f..4a9f60c 100644 --- a/src/Makefile +++ b/src/Makefile @@ -2,7 +2,7 @@ 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 MODES=mode_create.o mode_merge.o mode_convert.o -OPTIONS=opt_disl.o opt_group.o +OPTIONS=opt_disl.o opt_group.o opt_orient.o OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES) $(OPTIONS) call_option.o .SUFFIXES: @@ -40,4 +40,4 @@ call_mode.o : $(MODES) call_option.o : $(OPTIONS) $(MODES) $(OPTIONS) subroutines.o io.o : atoms.o box.o $(MODES) main.o : io.o -testfuncs.o elements.o mode_create.o $(OPTIONS): subroutines.o +testfuncs.o elements.o mode_create.o $(OPTIONS) $(MODES): subroutines.o diff --git a/src/call_option.f90 b/src/call_option.f90 index 0d7f9c7..e67454b 100644 --- a/src/call_option.f90 +++ b/src/call_option.f90 @@ -2,6 +2,7 @@ subroutine call_option(option, arg_pos) use parameters use opt_disl use opt_group + use opt_orient use box implicit none @@ -15,14 +16,16 @@ subroutine call_option(option, arg_pos) call group(arg_pos) case('-ow') arg_pos = arg_pos + 1 - continue case('-wrap') arg_pos = arg_pos + 1 - continue + case('-orient') + call orient(arg_pos) + case('-unorient') + call unorient + arg_pos = arg_pos + 1 case('-boundary') arg_pos=arg_pos+1 call get_command_argument(arg_pos, box_bc) - print *, box_bc arg_pos=arg_pos+1 case default print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.' diff --git a/src/elements.f90 b/src/elements.f90 index 068d188..b4b8814 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -5,6 +5,7 @@ module elements use parameters use functions use subroutines + use box implicit none !Data structures used to represent the CAC elements. Each index represents an element @@ -17,6 +18,7 @@ module elements !Data structure used to represent atoms integer, allocatable :: type_atom(:)!atom type + integer, allocatable :: sbox_atom(:) real(kind =dp),allocatable :: r_atom(:,:) !atom position integer :: atom_num=0 !Number of atoms @@ -153,7 +155,7 @@ module elements if(m > 0) then !Allocate atom arrays - allocate(type_atom(m), r_atom(3,m), stat=allostat) + allocate(type_atom(m), sbox_atom(m), r_atom(3,m), stat=allostat) if(allostat > 0) then print *, "Error allocating atom arrays in elements.f90 because of: ", allostat stop @@ -211,6 +213,11 @@ module elements temp_int(atom_size+1:) = 0 call move_alloc(temp_int, type_atom) + allocate(temp_int(m+atom_num+buffer_size)) + temp_int(1:atom_size) = sbox_atom + temp_int(atom_size+1:) = 0 + call move_alloc(temp_int, sbox_atom) + allocate(temp_real(3,m+atom_num+buffer_size)) temp_real(:,1:atom_size) = r_atom temp_real(:, atom_size+1:) = 0.0_dp @@ -237,9 +244,9 @@ module elements end subroutine add_element - subroutine add_atom(type, r) + subroutine add_atom(type, sbox, r) !Subroutine which adds an atom to the atom arrays - integer, intent(in) :: type + integer, intent(in) :: type, sbox real(kind=dp), intent(in), dimension(3) :: r atom_num = atom_num+1 @@ -247,6 +254,7 @@ module elements call grow_ele_arrays(0,1) type_atom(atom_num) = type r_atom(:,atom_num) = r(:) + sbox_atom(atom_num) = sbox end subroutine add_atom @@ -449,4 +457,38 @@ module elements end do end subroutine wrap_atoms + subroutine def_new_box + !This subroutine calculates new box boundaries based on minimum and maximum nodal/atomic positions + integer :: i, j, inod, ibasis + + real(kind=dp) :: max_bd(3), min_bd(3) + + max_bd(:) = -huge(1.0_dp) + min_bd(:) = huge(1.0_dp) + + do i = 1, atom_num + do j = 1, 3 + if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + lim_zero + if (r_atom(j,i) < min_bd(j)) min_bd(j) = r_atom(j,i) - lim_zero + end do + end do + + do i = 1, ele_num + do inod = 1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + do j = 1, 3 + if (r_node(j,ibasis,inod,i) > max_bd(j)) max_bd(j) = r_node(j,ibasis,inod,i) + lim_zero + if (r_node(j,ibasis,inod,i) < min_bd(j)) min_bd(j) = r_node(j,ibasis,inod,i) -lim_zero + + end do + end do + end do + end do + + do j = 1, 3 + box_bd(2*j) = max_bd(j) + box_bd(2*j-1) = min_bd(j) + end do + end subroutine + end module elements \ No newline at end of file diff --git a/src/io.f90 b/src/io.f90 index c4ced7f..088459a 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -8,7 +8,7 @@ module io implicit none integer :: outfilenum = 0, infilenum = 0 - character(len=100) :: outfiles(10), infiles(10) + character(len=100) :: outfiles(100), infiles(100) logical :: force_overwrite public @@ -130,14 +130,14 @@ module io do i = 1, ele_num do inod = 1, ng_node(lat_ele(i)) do ibasis = 1, basisnum(lat_ele(i)) - write(11, '(a, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i) + write(11, '(i16, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i) end do end do end do !Write atom positions do i = 1, atom_num - write(11, '(a, 3f23.15)') type_atom(i), r_atom(:,i) + write(11, '(i16, 3f23.15)') type_atom(i), r_atom(:,i) end do !Finish writing @@ -520,7 +520,7 @@ module io !Write out atoms first do i = 1, atom_num - write(11,*) i, type_atom(i), r_atom(:,i) + write(11,*) i, type_atom(i), sbox_atom(i), r_atom(:,i) end do !Write out the elements, this is written in two stages, one line for the element and then 1 line for @@ -713,8 +713,8 @@ module io !Read the atoms do i = 1, in_atoms - read(11,*) j, type, r(:) - call add_atom(new_type_to_type(type), r+newdisplace) + read(11,*) j, type, sbox, r(:) + call add_atom(new_type_to_type(type), sbox, r+newdisplace ) end do !Read the elements diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 607f6d0..934895a 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -135,7 +135,7 @@ module mode_create if(lat_atom_num > 0) then do i = 1, lat_atom_num do ibasis = 1, basisnum(1) - call add_atom(basis_type(ibasis, 1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis)) + call add_atom(basis_type(ibasis, 1), 1, (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis)) end do end do deallocate(r_atom_lat) @@ -205,21 +205,23 @@ module mode_create do i = 1, 3 call get_command_argument(arg_pos, orient_string, arglen) if (arglen==0) STOP "Missing orientation in orient command of mode create" + call parse_ori_vec(orient_string, orient(i,:)) 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 + + ! 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 end do diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 08167a6..8e0fc08 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -229,7 +229,7 @@ module opt_group !here as well to make sure they are in the box do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3 call apply_periodic(r_interp(:,j)) - call add_atom(type_interp(j), r_interp(:,j)) + call add_atom(type_interp(j), sbox_ele(ie), r_interp(:,j)) end do end do From f37f19cbcf427c87b155308fa381c33a138649d0 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Mon, 3 Feb 2020 17:35:22 -0500 Subject: [PATCH 2/2] Forgot to add the option file --- src/opt_orient.f90 | 121 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 121 insertions(+) create mode 100644 src/opt_orient.f90 diff --git a/src/opt_orient.f90 b/src/opt_orient.f90 new file mode 100644 index 0000000..21d9d3f --- /dev/null +++ b/src/opt_orient.f90 @@ -0,0 +1,121 @@ +module opt_orient + !This module contains the orient option which allows for the reorientation + !of simulation cells. This can be used to create arbitrarily oriented dislocation or loops. + use parameters + use subroutines + use elements + use box + implicit none + + + real(kind=dp), save :: new_orient(3,3) + real(kind=dp), dimension(6) :: orig_box_bd + real(kind=dp), allocatable :: orig_sub_box_ori(:,:,:) + + public + contains + + subroutine orient(arg_pos) + + integer, intent(inout) :: arg_pos + + integer :: i, ibasis, inod + real(kind=dp) :: inv_sub_box_ori(3,3,sub_box_num) + + !First parse the orient command + call parse_orient(arg_pos) + + !Now rotate the basis. To do this we transform the basis to [100] [010] [001] and then + !transform to user specified basis. + + !Find all inverse orientation matrices for all sub_boxes + do i = 1, sub_box_num + call matrix_inverse(sub_box_ori, 3, inv_sub_box_ori) + end do + + !Now transform all atoms + do i = 1, atom_num + r_atom(:,i) = matmul(new_orient,matmul(inv_sub_box_ori(:,:,sbox_atom(i)),r_atom(:,i))) + end do + + !Now transform all elements + do i = 1, ele_num + do inod =1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + r_node(:,ibasis,inod,i) = matmul(new_orient,matmul(inv_sub_box_ori(:,:,sbox_ele(i)),r_node(:,ibasis,inod,i))) + end do + end do + end do + + !Now save the original sub_box_ori and overwrite them + if(allocated(orig_sub_box_ori)) deallocate(orig_sub_box_ori) + + allocate(orig_sub_box_ori(3,3,sub_box_num)) + orig_sub_box_ori = sub_box_ori + + !Now overwrite the orientations + do i = 1, sub_box_num + sub_box_ori(:,:,i) = new_orient + end do + + !Save original box boundaries + orig_box_bd = box_bd + + !Now find new box boundaries + call def_new_box + end subroutine orient + + subroutine parse_orient(arg_pos) + !This parses the orient option + integer, intent(inout) :: arg_pos + + integer :: i, arg_len + character(len=8) :: ori_string + + !Pull out the new user orientation + do i = 1, 3 + arg_pos = arg_pos + 1 + call get_command_argument(arg_pos, ori_string, arg_len) + if (arg_len == 0) print *, "Missing orientation vector in -orient option" + call parse_ori_vec(ori_string, new_orient(i,:)) + end do + + !Normalize the orientation matrix + new_orient = matrix_normal(new_orient,3) + arg_pos = arg_pos + 1 + + end subroutine parse_orient + + subroutine unorient + + integer :: i, ibasis, inod + real(kind=dp) :: inv_ori(3,3) + + !Now rotate the basis. To do this we transform the basis to [100] [010] [001] and then + !transform to the original sbox_ele + + !Find the inverse for the new orientation matrix + call matrix_inverse(new_orient, 3, inv_ori) + + !Recover original sub_box_ori + sub_box_ori = orig_sub_box_ori + + !Now transform all atoms + do i = 1, atom_num + r_atom(:,i) = matmul(sub_box_ori(:,:,sbox_atom(i)),matmul(inv_ori(:,:),r_atom(:,i))) + end do + + !Now transform all elements + do i = 1, ele_num + do inod =1, ng_node(lat_ele(i)) + do ibasis = 1, basisnum(lat_ele(i)) + r_node(:,ibasis,inod,i) = matmul(sub_box_ori(:,:,sbox_ele(i)),matmul(inv_ori,r_node(:,ibasis,inod,i))) + end do + end do + end do + + !Restore original box boundaries + box_bd = orig_box_bd + end subroutine unorient + +end module opt_orient \ No newline at end of file