diff --git a/src/Makefile.dep b/src/Makefile.dep index b1ecea9..861d4a4 100644 --- a/src/Makefile.dep +++ b/src/Makefile.dep @@ -1,185 +1,26 @@ -# This file is generated automatically. DO NOT EDIT! - -obj/main : \ - obj/atoms.o \ - obj/box.o \ - obj/caller.o \ - obj/elements.o \ - obj/functions.o \ - obj/io.o \ - obj/main.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/neighbors.o \ - obj/opt_bubble.o \ - obj/opt_deform.o \ - obj/opt_delete.o \ - obj/opt_disl.o \ - obj/opt_group.o \ - obj/opt_orient.o \ - obj/opt_redef_box.o \ - obj/opt_slip_plane.o \ - obj/parameters.o \ - obj/sorts.o \ - obj/str.o \ - obj/subroutines.o - -obj/atoms.o : \ - obj/functions.o \ - obj/parameters.o - -obj/box.o : \ - obj/functions.o \ - obj/parameters.o - -obj/caller.o : \ - obj/box.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_bubble.o \ - obj/opt_deform.o \ - obj/opt_delete.o \ - obj/opt_disl.o \ - obj/opt_group.o \ - obj/opt_orient.o \ - obj/opt_redef_box.o \ - obj/opt_slip_plane.o \ - obj/parameters.o - -obj/elements.o : \ - obj/box.o \ - obj/functions.o \ - obj/parameters.o \ - obj/sorts.o \ - obj/subroutines.o - -obj/functions.o : \ - obj/parameters.o - -obj/io.o : \ - obj/atoms.o \ - obj/box.o \ - obj/elements.o \ - obj/parameters.o \ - obj/str.o - -obj/main.o : \ - obj/caller.o \ - obj/elements.o \ - obj/io.o \ - obj/parameters.o - -obj/mode_calc.o : \ - obj/box.o \ - obj/elements.o \ - obj/io.o \ - obj/parameters.o \ - obj/subroutines.o - -obj/mode_convert.o : \ - obj/box.o \ - obj/elements.o \ - obj/io.o \ - obj/parameters.o - -obj/mode_create.o : \ - obj/atoms.o \ - obj/box.o \ - obj/elements.o \ - obj/io.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 \ - obj/io.o \ - obj/neighbors.o \ - obj/parameters.o \ - obj/subroutines.o - -obj/mode_metric.o : \ - obj/elements.o \ - obj/io.o \ - obj/neighbors.o \ - obj/parameters.o - -obj/neighbors.o : \ - obj/elements.o \ - obj/functions.o \ - obj/parameters.o \ - obj/subroutines.o - -obj/opt_bubble.o : \ - obj/box.o \ - obj/elements.o \ - obj/opt_group.o \ - obj/parameters.o - -obj/opt_deform.o : \ - obj/box.o \ - obj/elements.o \ - obj/parameters.o \ - obj/subroutines.o - -obj/opt_delete.o : \ - obj/elements.o \ - obj/neighbors.o \ - obj/parameters.o \ - obj/subroutines.o - -obj/opt_disl.o : \ - obj/box.o \ - obj/elements.o \ - obj/parameters.o \ - obj/subroutines.o - -obj/opt_group.o : \ - obj/box.o \ - obj/elements.o \ - obj/parameters.o \ - obj/sorts.o \ - obj/subroutines.o - -obj/opt_orient.o : \ - obj/box.o \ - obj/elements.o \ - obj/parameters.o \ - obj/subroutines.o - -obj/opt_redef_box.o : \ - obj/box.o \ - obj/elements.o \ - obj/subroutines.o - -obj/opt_slip_plane.o : \ - obj/elements.o \ - obj/functions.o \ - obj/parameters.o \ - obj/subroutines.o - -obj/parameters.o : - -obj/sorts.o : \ - obj/parameters.o - -obj/str.o : - -obj/subroutines.o : \ - obj/box.o \ - obj/functions.o \ - obj/parameters.o +obj/atoms.o : atoms.f90 obj/functions.o obj/parameters.o +obj/box.o : box.f90 obj/functions.o obj/parameters.o +obj/caller.o : caller.f90 obj/box.o obj/opt_bubble.o obj/opt_slip_plane.o obj/opt_redef_box.o obj/opt_delete.o obj/opt_deform.o obj/opt_orient.o obj/opt_group.o obj/opt_disl.o obj/parameters.o obj/mode_calc.o obj/mode_metric.o obj/mode_merge.o obj/mode_convert.o obj/mode_create.o obj/mode_da.o +obj/elements.o : elements.f90 obj/box.o obj/sorts.o obj/subroutines.o obj/functions.o obj/parameters.o +obj/functions.o : functions.f90 obj/parameters.o +obj/io.o : io.f90 obj/str.o obj/box.o obj/atoms.o obj/parameters.o obj/elements.o +obj/main.o : main.f90 obj/caller.o obj/io.o obj/elements.o obj/parameters.o +obj/mode_calc.o : mode_calc.f90 obj/box.o obj/elements.o obj/subroutines.o obj/io.o obj/parameters.o +obj/mode_convert.o : mode_convert.f90 obj/io.o obj/elements.o obj/box.o obj/parameters.o +obj/mode_create.o : mode_create.f90 obj/box.o obj/elements.o obj/subroutines.o obj/io.o obj/atoms.o obj/parameters.o +obj/mode_da.o : mode_da.f90 obj/neighbors.o obj/elements.o obj/io.o obj/parameters.o +obj/mode_merge.o : mode_merge.f90 obj/neighbors.o obj/elements.o obj/subroutines.o obj/io.o obj/atoms.o obj/parameters.o +obj/mode_metric.o : mode_metric.f90 obj/neighbors.o obj/elements.o obj/io.o obj/parameters.o +obj/neighbors.o : neighbors.f90 obj/functions.o obj/subroutines.o obj/elements.o obj/parameters.o +obj/opt_bubble.o : opt_bubble.f90 obj/opt_group.o obj/box.o obj/elements.o obj/parameters.o +obj/opt_deform.o : opt_deform.f90 obj/box.o obj/elements.o obj/subroutines.o obj/parameters.o +obj/opt_delete.o : opt_delete.f90 obj/neighbors.o obj/elements.o obj/subroutines.o obj/parameters.o +obj/opt_disl.o : opt_disl.f90 obj/box.o obj/subroutines.o obj/elements.o obj/parameters.o +obj/opt_group.o : opt_group.f90 obj/sorts.o obj/box.o obj/subroutines.o obj/elements.o obj/parameters.o +obj/opt_orient.o : opt_orient.f90 obj/box.o obj/elements.o obj/subroutines.o obj/parameters.o +obj/opt_redef_box.o : opt_redef_box.f90 obj/subroutines.o obj/elements.o obj/box.o +obj/opt_slip_plane.o : opt_slip_plane.f90 obj/subroutines.o obj/functions.o obj/elements.o obj/parameters.o +obj/parameters.o : parameters.f90 +obj/sorts.o : sorts.f90 obj/parameters.o +obj/str.o : str.f90 +obj/subroutines.o : subroutines.f90 obj/str.o obj/box.o obj/functions.o obj/parameters.o diff --git a/src/box.f90 b/src/box.f90 index 9f2c7b2..8b985b3 100644 --- a/src/box.f90 +++ b/src/box.f90 @@ -5,14 +5,11 @@ module box implicit none real(kind=dp) :: box_bd(6) !Global box boundaries + real(kind=dp) :: box_ori(3,3) !Box orientation character(len=3) :: box_bc !Box boundary conditions (periodic or shrinkwrapped) logical :: bound_called !The subbox variables contain values for each subbox, being the boxes read in through some - !command. Currently only mode_merge will require sub_boxes, for mode_create it will always - !allocate to only 1 sub_box - integer :: sub_box_num = 0 - real(kind=dp), allocatable :: sub_box_ori(:,:,:)!Orientations for each of the subboxes - real(kind=dp), allocatable :: sub_box_bd(:,:)!Boundaries for each of the sub_boxes + !command. !Below are some simulation parameters which may be adjusted if reading in restart files integer :: timestep=0 @@ -29,42 +26,6 @@ module box bound_called=.false. end subroutine box_init - subroutine alloc_sub_box(n) - !Allocate the sub_box variables - - integer, intent(in) :: n - - integer :: i - - allocate(sub_box_ori(3,3,n), sub_box_bd(6,n)) - do i = 1, n - sub_box_ori(:,:,i) = identity_mat(3) - sub_box_bd(:,i) = 0.0_dp - end do - end subroutine alloc_sub_box - - subroutine grow_sub_box(n) - !Grows sub box arrays, this is only called when a new file is read in - integer, intent(in) :: n - - integer, allocatable :: temp_array_bd(:,:,:), temp_file(:) - real(kind=dp), allocatable :: temp_ori(:,:,:), temp_bd(:,:) - !Allocate temporary arrays - allocate(temp_ori(3,3,sub_box_num+n),temp_bd(6,sub_box_num+n), & - temp_array_bd(2,2,sub_box_num+n), temp_file(sub_box_num+n)) - - !Move allocation for all sub_box_arrays - temp_ori(:,:,1:sub_box_num) = sub_box_ori - temp_ori(:,:,sub_box_num+1:) = 0.0_dp - call move_alloc(temp_ori, sub_box_ori) - - temp_bd(:, 1:sub_box_num) = sub_box_bd - temp_bd(:, sub_box_num+1:) = 0.0_dp - call move_alloc(temp_bd, sub_box_bd) - - return - end subroutine grow_sub_box - subroutine grow_box(temp_box_bd) !This function takes in a temporary box boundary and adjusts the overall box boundaries !to include it @@ -85,24 +46,6 @@ 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 - subroutine reset_box !This subroutine just resets the box boundary and position box_bc = "ppp" @@ -115,4 +58,5 @@ module box box_volume = (box_bd(2) - box_bd(1)) * (box_bd(4) - box_bd(3))*(box_bd(6) - box_bd(5)) return end function + end module box diff --git a/src/caller.f90 b/src/caller.f90 index 2eeca41..a95fbb3 100644 --- a/src/caller.f90 +++ b/src/caller.f90 @@ -76,8 +76,8 @@ module caller call get_command_argument(arg_pos, box_bc) arg_pos=arg_pos+1 bound_called = .true. - case('-sbox_ori') - call sbox_ori(arg_pos) + case('-box_ori') + call set_box_ori(arg_pos) case('-deform') call deform(arg_pos) case('-delete') diff --git a/src/elements.f90 b/src/elements.f90 index 1dc58bf..655f4e5 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -11,7 +11,7 @@ module elements !Data structures used to represent the CAC elements. Each index represents an element character(len=100), allocatable :: type_ele(:) !Element type - integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:), tag_ele(:) !Element size + integer, allocatable :: size_ele(:), lat_ele(:), tag_ele(:) !Element size real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array !Element result data structures real(kind=dp), allocatable :: force_node(:,:,:,:), virial_node(:,:,:,:), energy_node(:,:,:), vel_node(:,:,:,:) @@ -22,7 +22,7 @@ module elements !Data structure used to represent atoms integer, allocatable :: type_atom(:)!atom type - integer, allocatable :: sbox_atom(:), tag_atom(:) + integer, allocatable :: tag_atom(:) real(kind =dp),allocatable :: r_atom(:,:) !atom position integer :: atom_num=0 !Number of atoms !Atom result data structures information @@ -100,8 +100,8 @@ module elements cubic_faces(:,1) = (/ 1, 2, 3, 4 /) cubic_faces(:,2) = (/ 1, 2, 6, 5 /) cubic_faces(:,3) = (/ 2, 3, 7, 6 /) - cubic_faces(:,4) = (/ 3, 4, 8, 7 /) - cubic_faces(:,5) = (/ 1, 4, 8, 5 /) + cubic_faces(:,4) = (/ 1, 4, 8, 5 /) + cubic_faces(:,5) = (/ 4, 3, 7, 8 /) cubic_faces(:,6) = (/ 5, 6, 7, 8 /) !!Now initialize the fcc primitive cell @@ -204,7 +204,7 @@ module elements !Allocate element arrays if(n > 0) then - allocate(type_ele(n), tag_ele(n), size_ele(n), lat_ele(n), sbox_ele(n), r_node(3,max_basisnum, max_ng_node,n), & + allocate(type_ele(n), tag_ele(n), size_ele(n), lat_ele(n), r_node(3,max_basisnum, max_ng_node,n), & stat=allostat) if(allostat > 0) then print *, "Error allocating element arrays in elements.f90 because of: ", allostat @@ -214,7 +214,7 @@ module elements if(m > 0) then !Allocate atom arrays - allocate(type_atom(m), sbox_atom(m), tag_atom(m), r_atom(3,m), stat=allostat) + allocate(type_atom(m), tag_atom(m), r_atom(3,m), stat=allostat) if(allostat > 0) then print *, "Error allocating atom arrays in elements.f90 because of: ", allostat stop @@ -258,11 +258,6 @@ module elements temp_int(ele_size+1:) = 0 call move_alloc(temp_int, size_ele) - allocate(temp_int(n+ele_size+buffer_size)) - temp_int(1:ele_size) = sbox_ele(1:ele_size) - temp_int(ele_size+1:) = 0 - call move_alloc(temp_int, sbox_ele) - allocate(char_temp(n+ele_size+buffer_size)) char_temp(1:ele_size) = type_ele(1:ele_size) call move_alloc(char_temp, type_ele) @@ -290,11 +285,6 @@ module elements temp_int(atom_size+1:) = 0 call move_alloc(temp_int, tag_atom) - allocate(temp_int(m+atom_size+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_size+buffer_size)) temp_real(:,1:atom_size) = r_atom temp_real(:, atom_size+1:) = 0.0_dp @@ -305,9 +295,9 @@ module elements end if end subroutine - subroutine add_element(tag, type, size, lat, sbox, r) + subroutine add_element(tag, type, size, lat, r) !Subroutine which adds an element to the element arrays - integer, intent(in) :: size, lat, sbox, tag + integer, intent(in) :: size, lat, tag character(len=100), intent(in) :: type real(kind=dp), intent(in) :: r(:,:,:) @@ -329,15 +319,14 @@ module elements type_ele(ele_num) = type size_ele(ele_num) = size lat_ele(ele_num) = lat - sbox_ele(ele_num) = sbox r_node(:,:,:,ele_num) = r(:,:,:) end subroutine add_element - subroutine add_atom(tag, type, sbox, r) + subroutine add_atom(tag, type, r) !Subroutine which adds an atom to the atom arrays - integer, intent(in) :: type, sbox, tag + integer, intent(in) :: type, tag real(kind=dp), intent(in), dimension(3) :: r integer :: newtag @@ -353,7 +342,6 @@ module elements tag_atom(atom_num) = newtag type_atom(atom_num) = type r_atom(:,atom_num) = r(:) - sbox_atom(atom_num) = sbox end subroutine add_atom @@ -641,14 +629,12 @@ module elements type_ele(sorted_index(i)) ='' size_ele(sorted_index(i)) = 0 lat_ele(sorted_index(i)) = 0 - sbox_ele(sorted_index(i)) = 0 tag_ele(sorted_index(i)) = 0 else r_node(:,:,:,sorted_index(i)) = r_node(:,:,:,ele_num) type_ele(sorted_index(i)) = type_ele(ele_num) size_ele(sorted_index(i)) = size_ele(ele_num) lat_ele(sorted_index(i)) = lat_ele(ele_num) - sbox_ele(sorted_index(i)) = sbox_ele(ele_num) tag_ele(sorted_index(i)) = tag_ele(ele_num) end if ele_num = ele_num - 1 @@ -878,11 +864,11 @@ do i = 1, atom_num case('fcc') !First we have to extract the element lattice parameter - call matrix_inverse(sub_box_ori(:,:,sbox_ele(ie)),3,ori_inv) + call matrix_inverse(box_ori(:,:),3,ori_inv) r_cubic_node = r_node(:,1,:,ie) r_cubic_node = matmul(fcc_inv,matmul(ori_inv,r_cubic_node)) lp = (maxval(r_cubic_node(1,:))-minval(r_cubic_node(1,:)))/(esize-1) - pos = matmul(sub_box_ori(:,:,sbox_ele(ie)),matmul(fcc_mat,pos))*lp + pos = matmul(box_ori,matmul(fcc_mat,pos))*lp pos = pos + r_node(:,1, 1, ie) case default diff --git a/src/io.f90 b/src/io.f90 index de8a7bc..a81f67c 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -104,7 +104,9 @@ module io case('restart') call write_pycac(outfiles(i)) case('cac') - call write_lmpcac(outfiles(i)) + print *, "Compatibility with lammpscac has been broken" + stop 1 +! call write_lmpcac(outfiles(i)) case('dump') call write_ldump(outfiles(i)) case default @@ -334,86 +336,86 @@ module io end if end subroutine write_ldump - subroutine write_lmpcac(file) - !This subroutine writes out a .lmp style dump file - character(len=100), intent(in) :: file - integer :: write_num, i, inod, ibasis - real(kind=dp) :: mass, fcc_adjust(3,8), local_adjust(3,8), rout(3) - -1 format(i16, ' Eight_Node', 4i16) -2 format(i16, ' Atom', 4i16) -3 format(3i16,3f23.15) - - open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') - - !Comment line - write(11, '(a)') '# CAC input file made with cacmb' - write(11, '(a)') - !Calculate total atom number - write_num = atom_num + ele_num - - !Write total number of atoms + elements - write(11, '(i16, a)') write_num, ' cac elements' - !Write number of atom types - write(11, '(i16, a)') atom_types, ' atom types' - - write(11,'(a)') ' ' - !Write box bd - write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi' - write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi' - write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi' - - !Masses - write(11, '(a)') 'Masses' - - write(11, '(a)') ' ' - do i =1, atom_types - call atommass(type_to_name(i),mass) - write(11, '(i16, f23.15, 2a)') i, mass, ' # ', type_to_name(i) - end do - write(11, '(a)') ' ' - - write(11, '(a)') 'CAC Elements' - write(11, '(a)') ' ' - - !Set up the nodal adjustment variables for all the different element types. This adjusts the node centers - !from the center of the unit cell (as formulated in this code) to the corners of the unit cells - do inod = 1, 8 - do i = 1,3 - if(is_equal(cubic_cell(i, inod),0.0_dp)) then - fcc_adjust(i,inod) = -0.5_dp - else - fcc_adjust(i, inod) = 0.5_dp - end if - end do - end do - fcc_adjust = matmul(fcc_mat, fcc_adjust) - - !Write element nodal positions - do i = 1, ele_num - select case(trim(adjustl(type_ele(i)))) - case('fcc') - !Now orient the current adjustment vector to the correct orientation - local_adjust = matmul(sub_box_ori(:,:,sbox_ele(i)), fcc_adjust) * lapa(lat_ele(i)) - !The first entry is the element specifier - write(11,1) i, basisnum(lat_ele(i)), size_ele(i), size_ele(i), size_ele(i) - do ibasis = 1, basisnum(lat_ele(i)) - do inod = 1, 8 - !Nodal information for every node - rout = r_node(:,ibasis,inod,i) + local_adjust(:,inod) - write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), rout - end do - end do - end select - end do - - do i = 1, atom_num - !Element specifier dictating that it is an atom - write(11,2) ele_num+i, 1, 1, 1, 1 - !Write the atomic information - write(11,3) 1, 1, type_atom(i), r_atom(:,i) - end do - end subroutine write_lmpcac +! subroutine write_lmpcac(file) +! !This subroutine writes out a .lmp style dump file +! character(len=100), intent(in) :: file +! integer :: write_num, i, inod, ibasis +! real(kind=dp) :: mass, fcc_adjust(3,8), local_adjust(3,8), rout(3) +! +!1 format(i16, ' Eight_Node', 4i16) +!2 format(i16, ' Atom', 4i16) +!3 format(3i16,3f23.15) +! +! open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') +! +! !Comment line +! write(11, '(a)') '# CAC input file made with cacmb' +! write(11, '(a)') +! !Calculate total atom number +! write_num = atom_num + ele_num +! +! !Write total number of atoms + elements +! write(11, '(i16, a)') write_num, ' cac elements' +! !Write number of atom types +! write(11, '(i16, a)') atom_types, ' atom types' +! +! write(11,'(a)') ' ' +! !Write box bd +! write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi' +! write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi' +! write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi' +! +! !Masses +! write(11, '(a)') 'Masses' +! +! write(11, '(a)') ' ' +! do i =1, atom_types +! call atommass(type_to_name(i),mass) +! write(11, '(i16, f23.15, 2a)') i, mass, ' # ', type_to_name(i) +! end do +! write(11, '(a)') ' ' +! +! write(11, '(a)') 'CAC Elements' +! write(11, '(a)') ' ' +! +! !Set up the nodal adjustment variables for all the different element types. This adjusts the node centers +! !from the center of the unit cell (as formulated in this code) to the corners of the unit cells +! do inod = 1, 8 +! do i = 1,3 +! if(is_equal(cubic_cell(i, inod),0.0_dp)) then +! fcc_adjust(i,inod) = -0.5_dp +! else +! fcc_adjust(i, inod) = 0.5_dp +! end if +! end do +! end do +! fcc_adjust = matmul(fcc_mat, fcc_adjust) +! +! !Write element nodal positions +! do i = 1, ele_num +! select case(trim(adjustl(type_ele(i)))) +! case('fcc') +! !Now orient the current adjustment vector to the correct orientation +! local_adjust = matmul(sub_box_ori(:,:,sbox_ele(i)), fcc_adjust) * lapa(lat_ele(i)) +! !The first entry is the element specifier +! write(11,1) i, basisnum(lat_ele(i)), size_ele(i), size_ele(i), size_ele(i) +! do ibasis = 1, basisnum(lat_ele(i)) +! do inod = 1, 8 +! !Nodal information for every node +! rout = r_node(:,ibasis,inod,i) + local_adjust(:,inod) +! write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), rout +! end do +! end do +! end select +! end do +! +! do i = 1, atom_num +! !Element specifier dictating that it is an atom +! write(11,2) ele_num+i, 1, 1, 1, 1 +! !Write the atomic information +! write(11,3) 1, 1, type_atom(i), r_atom(:,i) +! end do +! end subroutine write_lmpcac subroutine write_vtk(file) !This subroutine writes out a vtk style dump file @@ -624,16 +626,9 @@ module io !Open the .mb file for writing open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind') - !First write the box boundary information - !Write the global box boundaries + !First write the box boundary and orientation + write(11,*) box_ori(:,:) write(11,*) box_bd(:) - !Write the number of sub_boxes in the system - write(11,*) sub_box_num - !For every subbox write the orientation, sub box boundar - do i = 1, sub_box_num - write(11,*) sub_box_ori(:,:,i) - write(11,*) sub_box_bd(:,i) - 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) @@ -649,13 +644,13 @@ module io !Write out atoms first do i = 1, atom_num - write(11,*) tag_atom(i), type_atom(i), sbox_atom(i), r_atom(:,i) + write(11,*) tag_atom(i), type_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 !every basis at every node do i = 1, ele_num - write(11, *) tag_ele(i), lat_ele(i), size_ele(i), sbox_ele(i), trim(adjustl(type_ele(i))) + write(11, *) tag_ele(i), lat_ele(i), size_ele(i), trim(adjustl(type_ele(i))) do inod = 1, ng_node(lat_ele(i)) do ibasis =1, basisnum(lat_ele(i)) write(11,*) inod, ibasis, r_node(:, ibasis, inod, i) @@ -734,7 +729,9 @@ module io case('lmp') call read_lmp(infiles(i), displace, temp_box_bd) case('cac') - call read_lmpcac(infiles(i), displace, temp_box_bd) + print *, "Compatibility for lammpsCAC has been broken" + stop 1 +! call read_lmpcac(infiles(i), displace, temp_box_bd) case('out') call read_pycac_out(infiles(i), displace, temp_box_bd) case default @@ -754,7 +751,7 @@ module io real(kind = dp), dimension(6), intent(out) :: temp_box_bd integer :: i, j, k, l, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, & - new_type_to_type(10), new_lattice_types, sbox, new_lattice_map(10), btype_map(10,10), bnum_map(10), & + new_type_to_type(10), new_lattice_types, new_lattice_map(10), btype_map(10,10), bnum_map(10), & node_map(10) character(len=100) :: etype real(kind=dp) :: r(3), newdisplace(3), lapa_map(10) @@ -764,6 +761,7 @@ module io open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') !Read in the box boundary and grow the current active box bd + read(11, *) box_ori(:,:) read(11, *) temp_box_bd(:) do i = 1, 3 @@ -777,27 +775,6 @@ module io end do call grow_box(temp_box_bd) - !Read in the number of sub_boxes and allocate the variables - read(11, *) n - - if (sub_box_num == 0) then - call alloc_sub_box(n) - else - call grow_sub_box(n) - end if - - !Read in subbox orientations and boundaries - do i = 1, n - !Read in orientation with column major ordering - read(11,*) ((sub_box_ori(j, k, sub_box_num+i), j = 1, 3), k = 1, 3) - !Read in subbox boundaries - read(11,*) sub_box_bd(:,sub_box_num+i) - - do j = 1, 3 - sub_box_bd(2*j-1,sub_box_num+i) = sub_box_bd(2*j-1, sub_box_num+i) + displace(j) - sub_box_bd(2*j,sub_box_num+i) = sub_box_bd(2*j, sub_box_num+i) + displace(j) - end do - end do !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) @@ -840,21 +817,21 @@ module io !Read the atoms do i = 1, in_atoms - read(11,*) j, type, sbox, r(:) + read(11,*) j, type, r(:) r = r+newdisplace - call add_atom(j, new_type_to_type(type), sbox+sub_box_num, r) + call add_atom(j, new_type_to_type(type), r) end do !Read the elements do i = 1, in_eles - read(11, *) j, type, size, sbox, etype + read(11, *) j, type, size, etype do inod = 1, ng_node(type) do ibasis =1, basisnum(type) read(11,*) k, l, r_innode(:, ibasis, inod) r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace end do end do - call add_element(j, etype, size, new_lattice_map(type), sbox+sub_box_num, r_innode) + call add_element(j, etype, size, new_lattice_map(type), r_innode) end do !Close the file being read @@ -865,8 +842,6 @@ module io lattice_types = maxval(new_lattice_map) - sub_box_num = sub_box_num + n - call set_max_esize end subroutine read_mb @@ -942,18 +917,6 @@ module io call grow_box(temp_box_bd) - !Allocate sub_box - if (sub_box_num == 0) then - call alloc_sub_box(1) - else - call grow_sub_box(1) - end if - - !Because orientations and other needed sub_box information isn't really - !present within the restart file we just default a lot of it. - sub_box_ori(:,:,sub_box_num+1) = identity_mat(3) - sub_box_bd(:, sub_box_num+1) = temp_box_bd - !Read in more useless info do i = 1, 4 read(11,*) textholder @@ -976,7 +939,7 @@ module io end do ip = ip + 8 call lattice_map(bnum, btypes, 8, 1.0_dp, lat_type) - call add_element(j, etype, esize+1, lat_type, sub_box_num + 1, r_in(:,1:max_basisnum, 1:max_ng_node)) + call add_element(j, etype, esize+1, lat_type, r_in(:,1:max_basisnum, 1:max_ng_node)) end do end if @@ -993,14 +956,13 @@ module io stop end if r_in_atom = r_in_atom + newdisplace - call add_atom(j,atom_type_map(atom_type), sub_box_num + 1, r_in_atom) + call add_atom(j,atom_type_map(atom_type), r_in_atom) end do !Close file close(11) lattice_types = maxval(new_lattice_map) - sub_box_num = sub_box_num + 1 call set_max_esize end if @@ -1060,13 +1022,6 @@ module io call grow_box(temp_box_bd) - !Allocate sub_box boundaries - if (sub_box_num == 0) then - call alloc_sub_box(1) - else - call grow_sub_box(1) - end if - !Now read in masses for atoms read(11, '(a)') line j = tok_count(line) @@ -1078,14 +1033,6 @@ module io call add_atom_type(atomic_element, atom_type_map(i), all_new) end do - - - !Because orientations and other needed sub_box information isn't really - !present within the .cac file we just default a lot of it. - sub_box_ori(:,:,sub_box_num+1) = identity_mat(3) - sub_box_bd(:, sub_box_num+1) = temp_box_bd - sub_box_num = sub_box_num + 1 - if(in_atoms > 0 ) then !Read atom header read(11,'(a)') line @@ -1103,7 +1050,7 @@ module io else read(line,*) tag, type, ra(:), ea, fa(:), va(:) end if - call add_atom(tag, atom_type_map(type), sub_box_num, ra) + call add_atom(tag, atom_type_map(type), ra) call add_atom_data(atom_num, ea, fa, va) if(read_vel) vel_atom(:,atom_num) = vela end do @@ -1138,7 +1085,7 @@ module io else if (n == 20) then fcc = "20fcc" end if - call add_element(tag, fcc, esize+1, lat_type, sub_box_num, re) + call add_element(tag, fcc, esize+1, lat_type, re) call add_element_data(ele_num, ee, fe, ve) if(read_vel) vel_node(:, :, :, ele_num) = vele(:,1:max_basisnum, 1:max_ng_node) end do @@ -1199,20 +1146,6 @@ module io !Grow box boundaries call grow_box(temp_box_bd) - !Allocate sub_box - if (sub_box_num == 0) then - call alloc_sub_box(1) - else - call grow_sub_box(1) - end if - - !Because orientations and other needed sub_box information isn't really - !present within the .cac file we just default a lot of it. - sub_box_ori(:,:,sub_box_num+1) = identity_mat(3) - sub_box_bd(:, sub_box_num+1) = temp_box_bd - sub_box_num = sub_box_num + 1 - - !Read useless information read(11,*) textholder @@ -1230,7 +1163,7 @@ module io do i = 1, atom_in read(11,*) id, type_in, r_in(:) r_in(:) = r_in + newdisplace - call add_atom(id, type_map(type_in), sub_box_num, r_in) + call add_atom(id, type_map(type_in), r_in) end do close(11) @@ -1287,20 +1220,6 @@ module io !Grow box boundaries call grow_box(temp_box_bd) - !Allocate sub_box - if (sub_box_num == 0) then - call alloc_sub_box(1) - else - call grow_sub_box(1) - end if - - !Because orientations and other needed sub_box information isn't really - !present within the .cac file we just default a lot of it. - sub_box_ori(:,:,sub_box_num+1) = identity_mat(3) - sub_box_bd(:, sub_box_num+1) = temp_box_bd - sub_box_num = sub_box_num + 1 - - !Read useless information read(11,*) textholder @@ -1352,11 +1271,11 @@ module io call lattice_map(in_basis, in_basis_types, 8, in_lapa, lat_type) !Now add the element - call add_element(0,in_lattice_type, esize, lat_type, sub_box_num, r_in(:,1:max_basisnum,1:max_ng_node)) + call add_element(0,in_lattice_type, esize, lat_type, r_in(:,1:max_basisnum,1:max_ng_node)) case('Atom') read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,1,1) - call add_atom(0,in_basis_types(ibasis), sub_box_num, r_in(:,1,1)) + call add_atom(0,in_basis_types(ibasis), r_in(:,1,1)) end select end do diff --git a/src/mode_create.f90 b/src/mode_create.f90 index bbebaa2..34e1bf0 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -122,7 +122,7 @@ module mode_create box_bd(2*i) = maxval(r_node_temp(i,:,:))+10.0_dp**(-6.0_dp) box_bd(2*i-1) = minval(r_node_temp(i,:,:)) - 10.0_dp**(-6.0_dp) end do - call add_element(0,element_type, esize, 1, 1, r_node_temp) + call add_element(0,element_type, esize, 1, r_node_temp) end if !If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays @@ -152,7 +152,7 @@ module mode_create if(lat_atom_num > 0) then do i = 1, lat_atom_num do ibasis = 1, basisnum(1) - call add_atom(0,basis_type(ibasis, 1), 1, (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis)) + call add_atom(0,basis_type(ibasis, 1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis)) end do end do deallocate(r_atom_lat) @@ -166,21 +166,16 @@ module mode_create end do end do - call add_element(0,element_type, elat(i), 1, 1, r_node_temp) + call add_element(0,element_type, elat(i), 1, r_node_temp) end do end if end if - !The last thing we do is setup the sub_box_boundaries - call alloc_sub_box(1) - sub_box_num = 1 - sub_box_ori(:,:,1) = orient - sub_box_bd(:,1) = box_bd - !If any elements are fully outside the box then wrap them back in if (any(crossb)) then call wrap_elements end if + box_ori = orient end subroutine create !This subroutine parses the command and pulls out information needed for mode_create subroutine parse_command(arg_pos) diff --git a/src/mode_da.f90 b/src/mode_da.f90 index 82a9bc9..8074d17 100644 --- a/src/mode_da.f90 +++ b/src/mode_da.f90 @@ -68,21 +68,29 @@ module mode_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) + real(kind = dp) :: face_centroids(3, ele_num*6), r(3), rc, vnode(3, max_basisnum, 4), vnorm(max_basisnum, 4), & + max_node_dist, ndiff, rmax(3), rmin(3), rnorm !Initialize variables l = 0 + max_node_dist = 0 !First calculate all of the face centroids do i = 1, ele_num do j = 1, 6 r(:) = 0.0_dp + rmax= -Huge(1.0) + rmin= Huge(1.0) do k = 1, 4 do ibasis = 1, basisnum(lat_ele(i)) r = r + r_node(:, ibasis, cubic_faces(k,j), i) + rmax= max(rmax, r_node(:, ibasis, cubic_faces(k,j), i)) + rmin= min(rmin, r_node(:, ibasis, cubic_faces(k,j), i)) end do end do - + + rnorm=norm2(rmax-rmin) + max_node_dist = max(max_node_dist, rnorm) r = r/(basisnum(lat_ele(i))*4) !add the face centroids, the type, and map the elements faces to the face arrays @@ -90,11 +98,15 @@ module mode_da face_centroids(:, l) = r face_types(l) = j face_ele(l) = i + + end do + !Calculate the largest difference between nodes for each face + end do - !Now calculate the nearest faces - rc = max_esize*maxval(lapa) + !Now set the cut off distance as half the distance from opposite nodes in the largest element + rc = 0.5*max_node_dist 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 @@ -102,8 +114,9 @@ module mode_da do i = 1, l nei = nn(i) - !Skip if it's 0 - if (nei == 0) cycle + !Skip if it's 0 or the closest face belongs to the same element + if ((nei == 0).or.(face_ele(i) == face_ele(nei))) 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 @@ -124,15 +137,18 @@ module mode_da end do do j = 1, 4 do ibasis = 1, basisnum(lat_ele(face_ele(i))) - vnorm = norm2(vnode(:, ibasis, j)) + vnorm(ibasis,j) = 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 + print *, "Element number ", face_ele(i), " is dislocated along face ", face_types(i), & + " with neighbor ", face_ele(nei), " with max displacement of ", maxval(vnorm) do j = 1, 4 is_slipped(:, cubic_faces(j,face_types(i)), face_ele(i)) = 1 + print *, j, vnode(:,1,j) end do end if end if diff --git a/src/mode_merge.f90 b/src/mode_merge.f90 index 1e8d291..da8eaa5 100644 --- a/src/mode_merge.f90 +++ b/src/mode_merge.f90 @@ -263,7 +263,7 @@ module mode_merge do ibasis = 1, basisnum(lat_ele(ie)) call apply_periodic(ratom(:,ibasis)) added_points=added_points + 1 - call add_atom(0, basis_type(ibasis,lat_ele(ie)), sbox_ele(ie), ratom(:,ibasis)) + call add_atom(0, basis_type(ibasis,lat_ele(ie)), ratom(:,ibasis)) end do end if end do diff --git a/src/neighbors.f90 b/src/neighbors.f90 index 98799ee..e66392f 100644 --- a/src/neighbors.f90 +++ b/src/neighbors.f90 @@ -174,7 +174,7 @@ module neighbors c = which_cell(:,i) !Initialize the min vec - rmin=Huge(1.0_dp) + rmin=rc_off !loop over all neighboring cells do ci = -1, 1, 1 @@ -201,7 +201,6 @@ module neighbors end do end do pointloop - print *, nn(1) return end subroutine calc_NN diff --git a/src/opt_bubble.f90 b/src/opt_bubble.f90 index 5eaf2f5..0dce055 100644 --- a/src/opt_bubble.f90 +++ b/src/opt_bubble.f90 @@ -68,7 +68,7 @@ module opt_bubble if (norm2(p-c) < br) exit end do - call add_atom(0, new_type, 1, p) + call add_atom(0, new_type, p) end do end subroutine bubble diff --git a/src/opt_disl.f90 b/src/opt_disl.f90 index e3f7aa9..7e5d13e 100644 --- a/src/opt_disl.f90 +++ b/src/opt_disl.f90 @@ -30,7 +30,7 @@ module opt_disl integer, intent(inout) :: arg_pos print *, '--------------------Option Dislocation-----------------------' - + select case(trim(adjustl(option))) case('-disl') call parse_disl(arg_pos) @@ -113,7 +113,7 @@ module opt_disl subroutine disl !This subroutine creates the actual dislocation - integer :: i, sub_box, inod, ibasis, ix, iy, iz + integer :: i, 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) @@ -240,7 +240,7 @@ module opt_disl subroutine dislgen !This subroutine creates the actual dislocation - integer :: i, sub_box, inod, ibasis + integer :: i, 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) @@ -250,14 +250,6 @@ module opt_disl 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) @@ -267,7 +259,7 @@ module opt_disl call matrix_inverse(ss_ori, 3, ss_inv) !Apply the rotation - disp_transform = matmul(sub_box_ori(:,:,i), ss_inv) + disp_transform = matmul(box_ori(:,:), ss_inv) call matrix_inverse(disp_transform,3,inv_transform) if(atom_num > 0) then @@ -380,18 +372,6 @@ module opt_disl 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 @@ -651,18 +631,6 @@ module opt_disl 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 parse_vacancydisloop diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 5e30420..776d161 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -70,8 +70,10 @@ module opt_group end if if(remesh_size > 0)then + print *, "Remesh command has been dropped" + stop 1 call get_group - call remesh_group +! call remesh_group end if if(group_type > 0) then @@ -84,8 +86,10 @@ module opt_group call displace_group end if if(insert_type > 0) then + print *, "Insert command has been dropped" + stop 1 call get_group - call insert_group +! call insert_group end if if(num_species > 0) then @@ -705,7 +709,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(0,type_interp(j), sbox_ele(ie), r_interp(:,j)) + call add_atom(0,type_interp(j), r_interp(:,j)) end do end do !Once all atoms are added we delete all of the elements @@ -778,7 +782,7 @@ module opt_group if (.not.in_group_ele(esize, lat_ele(ie), rfill)) then nump_ele=nump_ele - esize**3 lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) = .false. - call add_element(0,type_ele(ie), esize, lat_ele(ie), sbox_ele(ie), rfill) + call add_element(0,type_ele(ie), esize, lat_ele(ie), rfill) new_ele_num = new_ele_num + 1 added_points = added_points + esize**3 end if @@ -797,7 +801,7 @@ module opt_group call get_interp_pos(m,n,o, ie, ratom(:,:)) do ibasis = 1, basisnum(lat_ele(ie)) call apply_periodic(ratom(:,ibasis)) - call add_atom(0, basis_type(ibasis,lat_ele(ie)), sbox_ele(ie), ratom(:,ibasis)) + call add_atom(0, basis_type(ibasis,lat_ele(ie)), ratom(:,ibasis)) added_points=added_points + 1 end do end if @@ -818,270 +822,270 @@ module opt_group end subroutine refinefill_group - subroutine remesh_group - !This command is used to remesh the group to a desired element size - - integer :: i, j, k, ix, iy, iz, inod, ibasis, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num, & - current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3), new_ele, new_atom, & - max_loops, working_esize, group_lat_num, lat_list(10), group_sbox_num, sbox_list(100), is, ilat, tot_dof - - real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), ori_inv(3,3), r(3), & - r_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8), group_bd(6) - logical, allocatable :: lat_points(:,:,:), new_sbox, new_lat - character(len=100) :: remesh_ele_type - - - !Right now we just hardcode only remeshing to elements - remesh_ele_type = 'fcc' - - ! Determine which sub_boxes and lattices types are within in the group - group_sbox_num = 0 - group_lat_num = 0 - do i = 1, group_atom_num - new_sbox=.true. - new_lat=.true. - do j = 1, group_sbox_num - if (sbox_list(j) == sbox_atom(atom_index(i))) then - new_sbox=.false. - exit - end if - end do - - if(new_sbox) then - group_sbox_num = group_sbox_num + 1 - sbox_list(group_sbox_num) = sbox_atom(atom_index(i)) - end if - - do j = 1, group_lat_num - if (basis_type(1,lat_list(j)) == type_atom(atom_index(i))) then - new_lat = .false. - exit - end if - end do - - if (new_lat) then - group_lat_num = group_lat_num + 1 - do k = 1, lattice_types - if (basis_type(1,k) == type_atom(atom_index(i))) lat_list(group_lat_num) = k - end do - end if - - end do - - do i = 1, group_ele_num - new_sbox=.true. - new_lat = .true. - do j = 1, group_sbox_num - if (sbox_list(j) == sbox_ele(element_index(i))) then - new_sbox = .false. - exit - end if - end do - - if (new_sbox) then - group_sbox_num = group_sbox_num + 1 - sbox_list(group_sbox_num) = sbox_ele(element_index(i)) - end if - - do j = 1, group_lat_num - if (lat_list(group_lat_num) == lat_ele(element_index(i))) then - new_lat=.false. - exit - end if - end do - - if (new_lat) then - group_lat_num = group_lat_num + 1 - lat_list(group_lat_num) = lat_ele(element_index(i)) - end if - end do - - new_atom = 0 - new_ele=0 - tot_dof=0 - do is = 1, group_sbox_num - - orient = sub_box_ori(:, :, sbox_list(is)) - call matrix_inverse(orient,3,ori_inv) - - do ilat = 1, group_lat_num - - !First calculate max position in lattice space to be able to allocate lat_points array, also sum the total number - !of degrees of freedom which are added - dof = 0 - do j = 1, 3 - group_bd(2*j) = -huge(1.0_dp) - group_bd(2*j-1) = huge(1.0_dp) - end do - do i = 1, group_atom_num - if ((type_atom(atom_index(i)) == basis_type(1,ilat)).and.(sbox_atom(atom_index(i)) == is)) then - do j =1 ,3 - if (r_atom(j,atom_index(i)) > group_bd(2*j)) group_bd(2*j) = r_atom(j,atom_index(i)) - if (r_atom(j,atom_index(i)) < group_bd(2*j-1)) group_bd(2*j-1) = r_atom(j,atom_index(i)) - end do - dof = dof + 1 - end if - end do - - do i = 1, group_ele_num - if ((lat_ele(element_index(i)) == ilat).and.(sbox_ele(element_index(i)) == is)) then - do inod =1, ng_node(ilat) - do ibasis = 1, basisnum(ilat) - do j = 1, 3 - r =r_node(j,ibasis,inod,element_index(i)) - if (r(j) > group_bd(2*j)) group_bd(2*j) = r(j) - if (r(j) < group_bd(2*j-1)) group_bd(2*j-1) = r(j) - end do - end do - end do - dof = dof + size_ele(element_index(i))**3 - end if - end do - - !If for some reason there are no dof in this loop then cycle out - if(dof == 0) cycle - tot_dof = tot_dof+dof - - group_in_lat = reshape((/ group_bd(1),group_bd(3),group_bd(5), & - group_bd(2),group_bd(3),group_bd(5), & - group_bd(2),group_bd(4),group_bd(5), & - group_bd(1),group_bd(4),group_bd(5), & - group_bd(1),group_bd(3),group_bd(6), & - group_bd(2),group_bd(3),group_bd(6), & - group_bd(2),group_bd(4),group_bd(6), & - group_bd(1),group_bd(4),group_bd(6) /), [3,8]) - - group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/lapa(ilat))) - do i = 1, 3 - bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:))) - bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:))) - end do - - allocate(lat_points(bd_in_lat(2)-bd_in_lat(1)+10, bd_in_lat(4)-bd_in_lat(3)+10, bd_in_lat(6)-bd_in_lat(5)+10)) - lat_points(:,:,:) = .false. - dof = 0 - - !Now place all group atoms and group interpolated atoms into lat_points - do i = 1, group_atom_num - if (.not.((sbox_atom(atom_index(i)) == is).and.(basis_type(1,ilat) == type_atom(atom_index(i))))) cycle - r = r_atom(:,atom_index(i))/lapa(ilat) - r = matmul(fcc_inv,matmul(ori_inv,r)) - do j = 1, 3 - r_lat(j) = nint(r(j)) - end do - !Do a check to make sure the code is working and that lattice points aren't being written on top of each other. - !This is primarily a debugging statement - if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then - stop "Multiple atoms share same position in lat point array, this shouldn't happen" - else - lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true. - dof = dof + 1 - end if - end do - - !Now place interpolated atoms within lat_points array - do i =1, group_ele_num - if (.not.((sbox_ele(element_index(i)) == is).and.( lat_ele(element_index(i)) == ilat))) cycle - ie = element_index(i) - call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) - do j = 1, size_ele(ie)**3 * basisnum(lat_ele(ie)) - r = r_interp(:,j)/lapa(ilat) - r = matmul(fcc_inv,matmul(ori_inv,r)) - do k = 1, 3 - r_lat(k) = nint(r(k)) - end do - !Do a check to make sure the code is working and that lattice points aren't being written on top of each - !other. This is primarily a debugging statement - if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then - stop "Multiple atoms/interpolated atoms share same position in lat point array, this shouldn't happen" - else - lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true. - dof = dof + 1 - end if - end do - end do - - !Now run remeshing algorithm, not the most optimized or efficient but gets the job done - !Figure out new looping boundaries - bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10 - bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10 - bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10 - - - if (max_remesh) then - max_loops = (remesh_size-3)/2 - else - max_loops = 1 - end if - do j = 1, max_loops - working_esize = remesh_size - 2*(j-1) - ele = (working_esize-1)*cubic_cell - zloop: do iz = 1, bd_in_array(3) - yloop: do iy = 1, bd_in_array(2) - xloop: do ix = 1, bd_in_array(1) - if (lat_points(ix, iy,iz)) then - r_new_node(:,:,:) = 0.0_dp - - !Check to see if the element overshoots the bound - if (iz+working_esize-1 > bd_in_array(3)) then - exit zloop - else if (iy+working_esize-1 > bd_in_array(2)) then - cycle zloop - else if (ix+working_esize-1 > bd_in_array(1)) then - cycle yloop - end if - - if (all(lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1))) then - do inod = 1, ng_node(ilat) - vlat = ele(:,inod) + (/ix, iy, iz /) - do i = 1, 3 - vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 - end do - r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat) - end do - - lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1) = .false. - !Add the element, for the sbox we just set it to the same sbox that we get the orientation - !from. In this case it is from the sbox of the first atom in the group. - new_ele = new_ele+1 - call add_element(0,remesh_ele_type, working_esize, ilat, & - sbox_atom(atom_index(1)),r_new_node) - - end if - end if - end do xloop - end do yloop - end do zloop - end do - - !Now we have to add any leftover lattice points as atoms - do iz = 1, bd_in_array(3) - do iy=1, bd_in_array(2) - do ix = 1, bd_in_array(1) - if(lat_points(ix,iy,iz)) then - vlat = (/ ix, iy, iz /) - do i = 1, 3 - vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 - end do - lat_points(ix,iy,iz) = .false. - r = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat) - new_atom=new_atom+1 - call add_atom(0,basis_type(1,ilat), is, r) - end if - end do - end do - end do - deallocate(lat_points) - end do - end do - - !Delete all elements and atoms to make space for new elements and atoms - call delete_atoms(group_atom_num, atom_index) - call delete_elements(group_ele_num, element_index) - - print *, tot_dof, " degrees of freedom in group" - print *, "remesh_group creates ", new_atom, " atoms and ", new_ele, " elements." - end subroutine remesh_group +! subroutine remesh_group +! !This command is used to remesh the group to a desired element size +! +! integer :: i, j, k, ix, iy, iz, inod, ibasis, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num, & +! current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3), new_ele, new_atom, & +! max_loops, working_esize, group_lat_num, lat_list(10), is, ilat, tot_dof +! +! real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), ori_inv(3,3), r(3), & +! r_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8), group_bd(6) +! logical, allocatable :: lat_points(:,:,:), new_sbox, new_lat +! character(len=100) :: remesh_ele_type +! +! +! !Right now we just hardcode only remeshing to elements +! remesh_ele_type = 'fcc' +! +! ! Determine which sub_boxes and lattices types are within in the group +! group_sbox_num = 0 +! group_lat_num = 0 +! do i = 1, group_atom_num +! new_sbox=.true. +! new_lat=.true. +! do j = 1, group_sbox_num +! if (sbox_list(j) == sbox_atom(atom_index(i))) then +! new_sbox=.false. +! exit +! end if +! end do +! +! if(new_sbox) then +! group_sbox_num = group_sbox_num + 1 +! sbox_list(group_sbox_num) = sbox_atom(atom_index(i)) +! end if +! +! do j = 1, group_lat_num +! if (basis_type(1,lat_list(j)) == type_atom(atom_index(i))) then +! new_lat = .false. +! exit +! end if +! end do +! +! if (new_lat) then +! group_lat_num = group_lat_num + 1 +! do k = 1, lattice_types +! if (basis_type(1,k) == type_atom(atom_index(i))) lat_list(group_lat_num) = k +! end do +! end if +! +! end do +! +! do i = 1, group_ele_num +! new_sbox=.true. +! new_lat = .true. +! do j = 1, group_sbox_num +! if (sbox_list(j) == sbox_ele(element_index(i))) then +! new_sbox = .false. +! exit +! end if +! end do +! +! if (new_sbox) then +! group_sbox_num = group_sbox_num + 1 +! sbox_list(group_sbox_num) = sbox_ele(element_index(i)) +! end if +! +! do j = 1, group_lat_num +! if (lat_list(group_lat_num) == lat_ele(element_index(i))) then +! new_lat=.false. +! exit +! end if +! end do +! +! if (new_lat) then +! group_lat_num = group_lat_num + 1 +! lat_list(group_lat_num) = lat_ele(element_index(i)) +! end if +! end do +! +! new_atom = 0 +! new_ele=0 +! tot_dof=0 +! do is = 1, group_sbox_num +! +! orient = sub_box_ori(:, :, sbox_list(is)) +! call matrix_inverse(orient,3,ori_inv) +! +! do ilat = 1, group_lat_num +! +! !First calculate max position in lattice space to be able to allocate lat_points array, also sum the total number +! !of degrees of freedom which are added +! dof = 0 +! do j = 1, 3 +! group_bd(2*j) = -huge(1.0_dp) +! group_bd(2*j-1) = huge(1.0_dp) +! end do +! do i = 1, group_atom_num +! if ((type_atom(atom_index(i)) == basis_type(1,ilat)).and.(sbox_atom(atom_index(i)) == is)) then +! do j =1 ,3 +! if (r_atom(j,atom_index(i)) > group_bd(2*j)) group_bd(2*j) = r_atom(j,atom_index(i)) +! if (r_atom(j,atom_index(i)) < group_bd(2*j-1)) group_bd(2*j-1) = r_atom(j,atom_index(i)) +! end do +! dof = dof + 1 +! end if +! end do +! +! do i = 1, group_ele_num +! if ((lat_ele(element_index(i)) == ilat).and.(sbox_ele(element_index(i)) == is)) then +! do inod =1, ng_node(ilat) +! do ibasis = 1, basisnum(ilat) +! do j = 1, 3 +! r =r_node(j,ibasis,inod,element_index(i)) +! if (r(j) > group_bd(2*j)) group_bd(2*j) = r(j) +! if (r(j) < group_bd(2*j-1)) group_bd(2*j-1) = r(j) +! end do +! end do +! end do +! dof = dof + size_ele(element_index(i))**3 +! end if +! end do +! +! !If for some reason there are no dof in this loop then cycle out +! if(dof == 0) cycle +! tot_dof = tot_dof+dof +! +! group_in_lat = reshape((/ group_bd(1),group_bd(3),group_bd(5), & +! group_bd(2),group_bd(3),group_bd(5), & +! group_bd(2),group_bd(4),group_bd(5), & +! group_bd(1),group_bd(4),group_bd(5), & +! group_bd(1),group_bd(3),group_bd(6), & +! group_bd(2),group_bd(3),group_bd(6), & +! group_bd(2),group_bd(4),group_bd(6), & +! group_bd(1),group_bd(4),group_bd(6) /), [3,8]) +! +! group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/lapa(ilat))) +! do i = 1, 3 +! bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:))) +! bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:))) +! end do +! +! allocate(lat_points(bd_in_lat(2)-bd_in_lat(1)+10, bd_in_lat(4)-bd_in_lat(3)+10, bd_in_lat(6)-bd_in_lat(5)+10)) +! lat_points(:,:,:) = .false. +! dof = 0 +! +! !Now place all group atoms and group interpolated atoms into lat_points +! do i = 1, group_atom_num +! if (.not.((sbox_atom(atom_index(i)) == is).and.(basis_type(1,ilat) == type_atom(atom_index(i))))) cycle +! r = r_atom(:,atom_index(i))/lapa(ilat) +! r = matmul(fcc_inv,matmul(ori_inv,r)) +! do j = 1, 3 +! r_lat(j) = nint(r(j)) +! end do +! !Do a check to make sure the code is working and that lattice points aren't being written on top of each other. +! !This is primarily a debugging statement +! if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then +! stop "Multiple atoms share same position in lat point array, this shouldn't happen" +! else +! lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true. +! dof = dof + 1 +! end if +! end do +! +! !Now place interpolated atoms within lat_points array +! do i =1, group_ele_num +! if (.not.((sbox_ele(element_index(i)) == is).and.( lat_ele(element_index(i)) == ilat))) cycle +! ie = element_index(i) +! call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) +! do j = 1, size_ele(ie)**3 * basisnum(lat_ele(ie)) +! r = r_interp(:,j)/lapa(ilat) +! r = matmul(fcc_inv,matmul(ori_inv,r)) +! do k = 1, 3 +! r_lat(k) = nint(r(k)) +! end do +! !Do a check to make sure the code is working and that lattice points aren't being written on top of each +! !other. This is primarily a debugging statement +! if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then +! stop "Multiple atoms/interpolated atoms share same position in lat point array, this shouldn't happen" +! else +! lat_points(r_lat(1)-bd_in_lat(1)+5, r_lat(2)-bd_in_lat(3)+5, r_lat(3)-bd_in_lat(5)+5) = .true. +! dof = dof + 1 +! end if +! end do +! end do +! +! !Now run remeshing algorithm, not the most optimized or efficient but gets the job done +! !Figure out new looping boundaries +! bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10 +! bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10 +! bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10 +! +! +! if (max_remesh) then +! max_loops = (remesh_size-3)/2 +! else +! max_loops = 1 +! end if +! do j = 1, max_loops +! working_esize = remesh_size - 2*(j-1) +! ele = (working_esize-1)*cubic_cell +! zloop: do iz = 1, bd_in_array(3) +! yloop: do iy = 1, bd_in_array(2) +! xloop: do ix = 1, bd_in_array(1) +! if (lat_points(ix, iy,iz)) then +! r_new_node(:,:,:) = 0.0_dp +! +! !Check to see if the element overshoots the bound +! if (iz+working_esize-1 > bd_in_array(3)) then +! exit zloop +! else if (iy+working_esize-1 > bd_in_array(2)) then +! cycle zloop +! else if (ix+working_esize-1 > bd_in_array(1)) then +! cycle yloop +! end if +! +! if (all(lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1))) then +! do inod = 1, ng_node(ilat) +! vlat = ele(:,inod) + (/ix, iy, iz /) +! do i = 1, 3 +! vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 +! end do +! r_new_node(:,1,inod) = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat) +! end do +! +! lat_points(ix:ix+working_esize-1,iy:iy+working_esize-1,iz:iz+working_esize-1) = .false. +! !Add the element, for the sbox we just set it to the same sbox that we get the orientation +! !from. In this case it is from the sbox of the first atom in the group. +! new_ele = new_ele+1 +! call add_element(0,remesh_ele_type, working_esize, ilat, & +! sbox_atom(atom_index(1)),r_new_node) +! +! end if +! end if +! end do xloop +! end do yloop +! end do zloop +! end do +! +! !Now we have to add any leftover lattice points as atoms +! do iz = 1, bd_in_array(3) +! do iy=1, bd_in_array(2) +! do ix = 1, bd_in_array(1) +! if(lat_points(ix,iy,iz)) then +! vlat = (/ ix, iy, iz /) +! do i = 1, 3 +! vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 +! end do +! lat_points(ix,iy,iz) = .false. +! r = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat) +! new_atom=new_atom+1 +! call add_atom(0,basis_type(1,ilat), is, r) +! end if +! end do +! end do +! end do +! deallocate(lat_points) +! end do +! end do +! +! !Delete all elements and atoms to make space for new elements and atoms +! call delete_atoms(group_atom_num, atom_index) +! call delete_elements(group_ele_num, element_index) +! +! print *, tot_dof, " degrees of freedom in group" +! print *, "remesh_group creates ", new_atom, " atoms and ", new_ele, " elements." +! end subroutine remesh_group subroutine delete_group !This subroutine deletes all atoms/elements within a group @@ -1120,71 +1124,71 @@ module opt_group end subroutine change_group_type - subroutine insert_group - !This code inserts atoms into interstitial sites. This only works on atoms within the group due to the limitations with the - !Coarse-graining methodology which doesn't allow for concentration fluctuations. - real(kind=dp) interstitial_sites(3,14), rand, rinsert(3) - integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi, sbox - integer, allocatable :: used_sites(:,:) - - !First save all of the displacement vectors from a lattice site to interstitial site - !The first 6 are the octohedral sites and the next 8 are the tetrahedral sites - interstitial_sites= reshape( (/ -0.5_dp, 0.0_dp, 0.0_dp, & - 0.5_dp, 0.0_dp, 0.0_dp, & - 0.0_dp,-0.5_dp, 0.0_dp, & - 0.0_dp, 0.5_dp, 0.0_dp, & - 0.0_dp, 0.0_dp,-0.5_dp, & - 0.0_dp, 0.0_dp, 0.5_dp, & - -0.25_dp,-0.25_dp,-0.25_dp, & - -0.25_dp,-0.25_dp, 0.25_dp, & - -0.25_dp, 0.25_dp,-0.25_dp, & - -0.25_dp, 0.25_dp, 0.25_dp, & - 0.25_dp,-0.25_dp,-0.25_dp, & - 0.25_dp,-0.25_dp, 0.25_dp, & - 0.25_dp, 0.25_dp,-0.25_dp, & - 0.25_dp, 0.25_dp, 0.25_dp /), & - shape(interstitial_sites)) - - !First we calculate the number of atoms needed based on the number of atoms in the group and the concentration - interstitial_sites=interstitial_sites*insert_lattice - - add_num = (insert_conc*group_atom_num)/(1-insert_conc) - allocate(used_sites(2,add_num)) - - print *, "Inserting ", add_num, " atoms as atom type ", insert_type - - !Now set up the random number generator for the desired interstitial type - select case(insert_site) - case(1) - rlo=1 - rhi=6 - case(2) - rlo=7 - rhi = 14 - case(3) - rlo=1 - rhi=14 - end select - - !Now add the atoms - i = 1 - addloop:do while ( i < add_num) - call random_number(rand) - rindex = int(1+rand*(group_atom_num-1)) - ia=atom_index(rindex) - call random_number(rand) - sindex = int(rlo+rand*(rhi-rlo)) - do j = 1, i - if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop - end do - rinsert = r_atom(:,ia) + matmul(sub_box_ori(:,:,sbox_atom(ia)),interstitial_sites(:,sindex)) - sbox = sbox_atom(ia) - call add_atom(0, insert_type, sbox, rinsert) - used_sites(1,i) = ia - used_sites(2,i) = sindex - i = i + 1 - end do addloop - end subroutine insert_group +! subroutine insert_group +! !This code inserts atoms into interstitial sites. This only works on atoms within the group due to the limitations with the +! !Coarse-graining methodology which doesn't allow for concentration fluctuations. +! real(kind=dp) interstitial_sites(3,14), rand, rinsert(3) +! integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi +! integer, allocatable :: used_sites(:,:) +! +! !First save all of the displacement vectors from a lattice site to interstitial site +! !The first 6 are the octohedral sites and the next 8 are the tetrahedral sites +! interstitial_sites= reshape( (/ -0.5_dp, 0.0_dp, 0.0_dp, & +! 0.5_dp, 0.0_dp, 0.0_dp, & +! 0.0_dp,-0.5_dp, 0.0_dp, & +! 0.0_dp, 0.5_dp, 0.0_dp, & +! 0.0_dp, 0.0_dp,-0.5_dp, & +! 0.0_dp, 0.0_dp, 0.5_dp, & +! -0.25_dp,-0.25_dp,-0.25_dp, & +! -0.25_dp,-0.25_dp, 0.25_dp, & +! -0.25_dp, 0.25_dp,-0.25_dp, & +! -0.25_dp, 0.25_dp, 0.25_dp, & +! 0.25_dp,-0.25_dp,-0.25_dp, & +! 0.25_dp,-0.25_dp, 0.25_dp, & +! 0.25_dp, 0.25_dp,-0.25_dp, & +! 0.25_dp, 0.25_dp, 0.25_dp /), & +! shape(interstitial_sites)) +! +! !First we calculate the number of atoms needed based on the number of atoms in the group and the concentration +! interstitial_sites=interstitial_sites*insert_lattice +! +! add_num = (insert_conc*group_atom_num)/(1-insert_conc) +! allocate(used_sites(2,add_num)) +! +! print *, "Inserting ", add_num, " atoms as atom type ", insert_type +! +! !Now set up the random number generator for the desired interstitial type +! select case(insert_site) +! case(1) +! rlo=1 +! rhi=6 +! case(2) +! rlo=7 +! rhi = 14 +! case(3) +! rlo=1 +! rhi=14 +! end select +! +! !Now add the atoms +! i = 1 +! addloop:do while ( i < add_num) +! call random_number(rand) +! rindex = int(1+rand*(group_atom_num-1)) +! ia=atom_index(rindex) +! call random_number(rand) +! sindex = int(rlo+rand*(rhi-rlo)) +! do j = 1, i +! if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop +! end do +! rinsert = r_atom(:,ia) + matmul(sub_box_ori(:,:,sbox_atom(ia)),interstitial_sites(:,sindex)) +! sbox = sbox_atom(ia) +! call add_atom(0, insert_type, sbox, rinsert) +! used_sites(1,i) = ia +! used_sites(2,i) = sindex +! i = i + 1 +! end do addloop +! end subroutine insert_group subroutine alloy_group !This subroutine randomizes the atom types to reach desired concentrations, this only operates on atoms diff --git a/src/opt_orient.f90 b/src/opt_orient.f90 index fd42cb1..c538a55 100644 --- a/src/opt_orient.f90 +++ b/src/opt_orient.f90 @@ -10,7 +10,7 @@ module opt_orient real(kind=dp), save :: new_orient(3,3) real(kind=dp), dimension(6) :: orig_box_bd - real(kind=dp), allocatable :: orig_sub_box_ori(:,:,:) + real(kind=dp), dimension(3,3) :: orig_box_ori character(len=3) :: orig_box_bc public @@ -22,7 +22,7 @@ module opt_orient integer :: i, j, k, ibasis, inod logical :: isortho, isrighthanded, matching - real(kind=dp) :: inv_sub_box_ori(3,3,sub_box_num) + real(kind=dp) :: inv_box_ori(3,3) !First parse the orient command @@ -31,36 +31,25 @@ module opt_orient !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(:,:,i), 3, inv_sub_box_ori(:,:,i)) - end do + call matrix_inverse(box_ori(:,:), 3, inv_box_ori(:,:)) !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))) + r_atom(:,i) = matmul(new_orient,matmul(inv_box_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(new_orient,matmul(inv_sub_box_ori(:,:,sbox_ele(i)),r_node(:,ibasis,inod,i))) + r_node(:,ibasis,inod,i) = matmul(new_orient,matmul(inv_box_ori(:,:),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 - + !Now save the original box_ori + orig_box_ori = box_ori + box_ori = new_orient !Save original box boundaries orig_box_bd = box_bd @@ -70,18 +59,16 @@ module opt_orient orig_box_bc = box_bc do i = 1,3 matching=.true. - sbox_loop:do j = 1, sub_box_num - do k = 1, 3 - if(.not.is_equal(orig_sub_box_ori(i,k,j), new_orient(i,k))) then - matching = .false. - exit sbox_loop - end if - end do - end do sbox_loop + do k = 1, 3 + if(.not.is_equal(orig_box_ori(i,k), new_orient(i,k))) then + matching = .false. + exit + end if + end do if(.not.matching) box_bc(i:i)='s' end do - call def_new_box + call def_new_box end subroutine orient_opt subroutine parse_orient(arg_pos) @@ -121,24 +108,24 @@ module opt_orient 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 + !transform to the original box_ori !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 + !Recover original box_ori and box_bd + box_ori = orig_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))) + r_atom(:,i) = matmul(box_ori(:,:),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))) + r_node(:,ibasis,inod,i) = matmul(box_ori(:,:),matmul(inv_ori,r_node(:,ibasis,inod,i))) end do end do end do @@ -148,29 +135,4 @@ module opt_orient box_bc = orig_box_bc end subroutine unorient - subroutine sbox_ori(arg_pos) - integer, intent(inout) :: arg_pos - - integer :: i, sbox_in, arg_len - real(kind = dp) :: new_orient(3,3) - character(len=100) :: textholder - character(len=20) :: ori_string - - arg_pos = arg_pos + 1 - call get_command_argument(arg_pos, textholder,arg_len) - if (arg_len== 0) stop 'Missing sbox in sbox_ori command' - read(textholder,*) sbox_in - - 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,:)) - new_orient(i,:) = new_orient(i,:) / norm2(new_orient(i,:)) - end do - - sub_box_ori(:,:,sbox_in) = new_orient - arg_pos = arg_pos + 1 - - end subroutine sbox_ori end module opt_orient diff --git a/src/opt_redef_box.f90 b/src/opt_redef_box.f90 index f6330e0..0d55b47 100644 --- a/src/opt_redef_box.f90 +++ b/src/opt_redef_box.f90 @@ -3,6 +3,7 @@ module opt_redef_box use box use elements use subroutines + implicit none character(len=3) :: redef_dim, new_bc @@ -13,7 +14,8 @@ module opt_redef_box subroutine redef_box(arg_pos) !This is the main calling function for opt_redef_box integer, intent(inout) :: arg_pos - integer :: i, inod, ibasis, iatom, delete_list(atom_num), delete_num, type_interp(max_basisnum*max_esize**3) + integer :: i, j, inod, ibasis, iatom, delete_list(atom_num), delete_num, type_interp(max_basisnum*max_esize**3), & + new_snum real(kind=dp):: r_interp(3, max_basisnum*max_esize**3) logical :: node_out(8) @@ -64,7 +66,7 @@ module opt_redef_box !loop over all interpolated atoms and add them to the system do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 if(in_block_bd(r_interp(:,iatom), new_bd)) then - call add_atom(0,type_interp(iatom), sbox_ele(i), r_interp(:,iatom)) + call add_atom(0,type_interp(iatom), r_interp(:,iatom)) end if end do end if @@ -128,4 +130,26 @@ module opt_redef_box arg_pos = arg_pos + 1 end subroutine parse_redef_box + subroutine set_box_ori(arg_pos) + + integer, intent(inout) :: arg_pos + + integer :: arg_len, i + real(kind=dp) :: new_orient(3,3) + character(len=100) :: textholder + character(len=20) :: ori_string + + 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,:)) + new_orient(i,:) = new_orient(i,:) / norm2(new_orient(i,:)) + end do + + box_ori(:,:) = new_orient + arg_pos = arg_pos + 1 + + end subroutine set_box_ori + end module opt_redef_box diff --git a/src/opt_slip_plane.f90 b/src/opt_slip_plane.f90 index 80da930..e000460 100644 --- a/src/opt_slip_plane.f90 +++ b/src/opt_slip_plane.f90 @@ -61,7 +61,7 @@ module opt_slip_plane call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3 call apply_periodic(r_interp(:,ia)) - call add_atom(0, type_interp(ia), sbox_ele(ie), r_interp(:,ia)) + call add_atom(0, type_interp(ia), r_interp(:,ia)) end do !If we are efilling then the code is slightly more complex else @@ -93,7 +93,7 @@ module opt_slip_plane if(.not.(spos < maxp).and.(spos > minp))then nump_ele = nump_ele - esize**3 lat_points(m:m+esize-1, n:n+esize-1, o:o+esize-1) = .false. - call add_element(0, type_ele(ie), esize, lat_ele(ie), sbox_ele(ie), rfill) + call add_element(0, type_ele(ie), esize, lat_ele(ie), rfill) new_ele_num = new_ele_num + 1 end if end do latloop @@ -110,7 +110,7 @@ module opt_slip_plane call get_interp_pos(m,n,o, ie, ratom(:,:)) do ibasis = 1, basisnum(lat_ele(ie)) call apply_periodic(r_atom(:,ibasis)) - call add_atom(0, basis_type(ibasis,lat_ele(ie)), sbox_ele(ie), ratom(:,ibasis)) + call add_atom(0, basis_type(ibasis,lat_ele(ie)), ratom(:,ibasis)) end do end if end do