aselimov/CACmb
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Refactor to remove sub boxes from CACmb code, functionality was not useful and became detrimental for generation of new models

Browse Source
This commit is contained in:
Alex Selimov 2022-01-26 20:48:16 -05:00
parent 971e0a5e8d
commit de8b4e168a
15 changed files with 567 additions and 909 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

211
src/Makefile.dep
View File

@ -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

62
src/box.f90
View File

@ -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

4
src/caller.f90
View File

@ -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')

38
src/elements.f90
View File

@ -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

287
src/io.f90
View File

@ -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

13
src/mode_create.f90
View File

@ -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)

30
src/mode_da.f90
View File

@ -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

2
src/mode_merge.f90
View File

@ -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

3
src/neighbors.f90
View File

@ -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

2
src/opt_bubble.f90
View File

@ -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

40
src/opt_disl.f90
View File

@ -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

672
src/opt_group.f90
View File

@ -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

78
src/opt_orient.f90
View File

@ -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

28
src/opt_redef_box.f90
View File

@ -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

6
src/opt_slip_plane.f90
View File

@ -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