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

Merge pull request #28 from aselimov/ft--rand-element-pos

Browse Source 
Ft  rand element pos
This commit is contained in:
aselimov 2020-03-04 21:07:11 -05:00 committed by GitHub
commit 5a1676cd66
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
16 changed files with 1037 additions and 199 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
Numbered_element.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 164 KiB

88
README.md
View File

@ -214,10 +214,26 @@ This command wraps atoms back into the simulation cell as though periodic bounda
**Remesh**
```
remesh esize
remesh esize lattice_parameter lattice_type
```
This command remeshes the atoms/elements within the group to the new element size `esize`. Currently only accepts an `esize` of 2 which refines it to full atomistics.
This command remeshes the atoms/elements within the group to the new element size `esize`. Currently only accepts an `esize` of 2 which refines it to full atomistics. When remeshing to atomistics the group can contain any orientations of elements but when remeshing to different finite elements, the group must contain all atoms/elements with the same orientation. `lattice_parameter` is the lattice parameter for the elements within the group and `lattice_type` is the lattice type (integer) that these new elements will be assigned to.
**Max**
```
max
```
This command attempts to reduce the degrees of freedom in the model by replacing them with graded elements. This code works by starting at elements with size `esize` and then checks all degrees of freedom to see which ones can be replaced by inserting the element. It then iterates over elements of `esize-2` to `esize=2` which is full atomic resolution.
**Delete**
```
delete
```
This command deletes all selected atoms and elements within the group.
### Option overwrite
@ -225,4 +241,70 @@ This command remeshes the atoms/elements within the group to the new element siz
-ow
```
If this option is passed then all files are automatically overwritten without asking the user.
If this option is passed then all files are automatically overwritten without asking the user.
### Option boundary
```
-boundary box_bc
```
This allows the user to specify the boundary conditions for the model being outputted. The format is a 3 character string with `p` indicating periodic and `s` indicating shrink-wrapped.
**Example:** `-boundary psp`
### Option Delete
```
-delete keywords
```
Delete requires the usage of additional keywords to specify which delete action will be taken. These additional keywords are below:
**overlap**
```
-delete overlap rc_off
```
This command will delete all overlapping atoms within a specific cutoff radius `rc_off`. This currently does not affect elements.
****
## Position Specification
Specifying positions in cacmb can be done through a variety of ways. Examples of each format is shown below.
`val` - Where `val` is a number, then that value in Angstroms is used as the position. As an example, `11.1` would be read in as a position of 11.1 $\AA$.
`-inf` - This specifies the lower box boundary in the specific dimension. The vector `-inf -inf -inf` specifies the bottom corner of the simulation cell which also acts as the simulation cell origin. The vector `-inf 10 3` instead puts only the x position at the simulation cell origin.
`inf` - Similar to `-inf` but references the upper boundary of the box in that dimension
`inf-val` - Using a minus sign reduces the position from the **upper boundary** by `val`. `inf-10` would be at a distance of $10 \AA$ from the upper boundary in that dimension.
`inf+val` - This increases the position from the **lower boundary**. `inf+10` would be a position $10\AA$ from the lower boundary within the cell.
`inf*val` - This gives you a fractional position in the simulation cell. As an example `inf*0.5` gives you the center point of the simulation cell.
`rand` - Returns a random position that lies within the simulation cell.
`rand[val1:val2]` - returns a random position that lies within the range
`rande[facenum]` - Returns a random position in an interelement boundary which is offset of the element face `facenum`. Face numbers are based on the which vertices comprise the face. Vertex numbers are shown in the figure below for the primitive fcc unit cell which is what the fcc rhombohedral element is based from. The face numbers are:
Face 1: [1,2,3,4]
Face 2: [1,2,6,5]
Face 3: [2,3,7,6]
Face 4: [3,4,8,7]
Face 5: [1,4,8,5]
Face 6: [5,6,7,8]
Image for vertex numbers is:
![](/home/alexselimov/Documents/CACmb/Numbered_element.png)

8
src/Makefile
View File

@ -1,8 +1,8 @@
FC=ifort
FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin
FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays
#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays
MODES=mode_create.o mode_merge.o mode_convert.o
OPTIONS=opt_disl.o opt_group.o
OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o
OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES) $(OPTIONS) call_option.o
.SUFFIXES:
@ -40,4 +40,4 @@ call_mode.o : $(MODES)
call_option.o : $(OPTIONS)
$(MODES) $(OPTIONS) subroutines.o io.o : atoms.o box.o
$(MODES) main.o : io.o
testfuncs.o elements.o mode_create.o $(OPTIONS): subroutines.o
testfuncs.o elements.o mode_create.o $(OPTIONS) $(MODES): subroutines.o

21
src/call_option.f90
View File

@ -2,6 +2,9 @@ subroutine call_option(option, arg_pos)
use parameters
use opt_disl
use opt_group
use opt_orient
use opt_delete
use box
implicit none
integer, intent(inout) :: arg_pos
@ -14,9 +17,21 @@ subroutine call_option(option, arg_pos)
call group(arg_pos)
case('-ow')
arg_pos = arg_pos + 1
continue
case default
print *, 'Option ', trim(adjustl(option)), ' is not currently accepted. Skipping to next argument'
case('-wrap')
arg_pos = arg_pos + 1
case('-orient')
call orient(arg_pos)
case('-unorient')
call unorient
arg_pos = arg_pos + 1
case('-boundary')
arg_pos=arg_pos+1
call get_command_argument(arg_pos, box_bc)
arg_pos=arg_pos+1
case('-delete')
call run_delete(arg_pos)
case default
print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.'
stop 3
end select
end subroutine call_option

214
src/elements.f90
View File

@ -5,6 +5,7 @@ module elements
use parameters
use functions
use subroutines
use box
implicit none
!Data structures used to represent the CAC elements. Each index represents an element
@ -17,6 +18,7 @@ module elements
!Data structure used to represent atoms
integer, allocatable :: type_atom(:)!atom type
integer, allocatable :: sbox_atom(:)
real(kind =dp),allocatable :: r_atom(:,:) !atom position
integer :: atom_num=0 !Number of atoms
@ -41,6 +43,8 @@ module elements
integer :: basis_type(10,10)
real(kind=dp) :: lapa(10)
!Additional module level variables we need
logical :: wrap_flag
public
contains
@ -60,11 +64,11 @@ module elements
shape(fcc_cell))
!Now we create a list containing the list of vertices needed to describe the 6 cube faces
cubic_faces(:,1) = (/ 1, 4, 8, 5 /)
cubic_faces(:,2) = (/ 2, 3, 7, 6 /)
cubic_faces(:,3) = (/ 1, 2, 6, 5 /)
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, 2, 3, 4 /)
cubic_faces(:,5) = (/ 1, 4, 8, 5 /)
cubic_faces(:,6) = (/ 5, 6, 7, 8 /)
!!Now initialize the fcc primitive cell
@ -151,7 +155,7 @@ module elements
if(m > 0) then
!Allocate atom arrays
allocate(type_atom(m), r_atom(3,m), stat=allostat)
allocate(type_atom(m), sbox_atom(m), r_atom(3,m), stat=allostat)
if(allostat > 0) then
print *, "Error allocating atom arrays in elements.f90 because of: ", allostat
stop
@ -209,6 +213,11 @@ module elements
temp_int(atom_size+1:) = 0
call move_alloc(temp_int, type_atom)
allocate(temp_int(m+atom_num+buffer_size))
temp_int(1:atom_size) = sbox_atom
temp_int(atom_size+1:) = 0
call move_alloc(temp_int, sbox_atom)
allocate(temp_real(3,m+atom_num+buffer_size))
temp_real(:,1:atom_size) = r_atom
temp_real(:, atom_size+1:) = 0.0_dp
@ -235,9 +244,9 @@ module elements
end subroutine add_element
subroutine add_atom(type, r)
subroutine add_atom(type, sbox, r)
!Subroutine which adds an atom to the atom arrays
integer, intent(in) :: type
integer, intent(in) :: type, sbox
real(kind=dp), intent(in), dimension(3) :: r
atom_num = atom_num+1
@ -245,6 +254,7 @@ module elements
call grow_ele_arrays(0,1)
type_atom(atom_num) = type
r_atom(:,atom_num) = r(:)
sbox_atom(atom_num) = sbox
end subroutine add_atom
@ -313,7 +323,7 @@ module elements
real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: r_interp !Interpolated atomic positions
!Internal variables
integer :: i, it, is, ir, ibasis, inod, ia, bnum, lat_type_temp
integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp
real(kind=dp), allocatable :: a_shape(:)
real(kind=dp) :: t, s, r
@ -390,7 +400,7 @@ module elements
integer, intent(in) :: num
integer, intent(inout), dimension(num) :: index
integer :: i, j
integer :: i
call heapsort(index)
@ -413,7 +423,7 @@ module elements
integer, intent(in) :: num
integer, intent(inout), dimension(num) :: index
integer :: i, j
integer :: i
call heapsort(index)
@ -436,5 +446,189 @@ module elements
end if
ele_num = ele_num - 1
end do
end subroutine delete_elements
subroutine wrap_atoms
!This subroutine wraps atoms back into the simulation cell if they have exited for any reason
integer :: i
do i = 1, atom_num
call apply_periodic(r_atom(:,i))
end do
end subroutine wrap_atoms
subroutine def_new_box
!This subroutine calculates new box boundaries based on minimum and maximum nodal/atomic positions
integer :: i, j, inod, ibasis
real(kind=dp) :: max_bd(3), min_bd(3)
max_bd(:) = -huge(1.0_dp)
min_bd(:) = huge(1.0_dp)
do i = 1, atom_num
do j = 1, 3
if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + lim_zero
if (r_atom(j,i) < min_bd(j)) min_bd(j) = r_atom(j,i) - lim_zero
end do
end do
do i = 1, ele_num
do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
do j = 1, 3
if (r_node(j,ibasis,inod,i) > max_bd(j)) max_bd(j) = r_node(j,ibasis,inod,i) + lim_zero
if (r_node(j,ibasis,inod,i) < min_bd(j)) min_bd(j) = r_node(j,ibasis,inod,i) -lim_zero
end do
end do
end do
end do
do j = 1, 3
box_bd(2*j) = max_bd(j)
box_bd(2*j-1) = min_bd(j)
end do
end subroutine
recursive subroutine parse_pos(i, pos_string, pos)
!This subroutine parses the pos command allowing for command which include inf
integer, intent(in) :: i !The dimension of the position
character(len=100), intent(in) :: pos_string !The position string
real(kind=dp), intent(out) :: pos !The output parsed position value
integer :: iospara, face, ele, randsize
real(kind=dp) :: rand, rone, rtwo, rand_ele_pos(3)
character(len=100) :: cone, ctwo
iospara = 0
if(trim(adjustl(pos_string)) == 'inf') then
pos=box_bd(2*i)
else if(trim(adjustl(pos_string)) == '-inf') then
pos=box_bd(2*i-1)
else if (trim(adjustl(pos_string)) == 'rand') then
call random_number(rand)
pos = (box_bd(2*i)-box_bd(2*i-1))*rand + box_bd(2*i-1)
else if (index(pos_string,'rande')>0) then
!First select a random element
call random_number(rand)
ele= 1 + floor(ele_num*rand)
!Now read the rest of the command which specifies the face we need and check
!to make sure it's an accepted face number
cone = pos_string(index(pos_string, '[')+1:index(pos_string, '[')+1)
!Check to see if we also pass an element size, if it was passed then make sure our
!random element is the right size
if(index(pos_string,':') > 0) then
ctwo = pos_string(index(pos_string, ':')+1:index(pos_string,']')-1)
read(ctwo, *) randsize
do while(randsize /= size_ele(ele))
call random_number(rand)
ele= 1 + floor(ele_num*rand)
end do
end if
read(cone, *) face
if ((face < 1).or.(face > 6)) stop "Current face number must be 1,2,3,4,5,6. Please check documentation"
!Now get the position
call offset_pos(ele, face, rand_ele_pos)
pos = rand_ele_pos(i)
else if (index(pos_string,'rand')>0) then
call random_number(rand)
cone = pos_string(index(pos_string, '[')+1:index(pos_string,':')-1)
call parse_pos(i, cone, rone)
ctwo = pos_string(index(pos_string, ':')+1:index(pos_string,']')-1)
call parse_pos(i, ctwo, rtwo)
pos = (rtwo - rone)*rand + rone
else if ((index(pos_string,'-') > 0).and.(index(pos_string,'inf')>0)) then
!Now extract the number we are reducing from infinity
if(index(pos_string,'inf') < index(pos_string,'-')) then
read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
else
read(pos_string(1:index(pos_string,'-')-1), *, iostat=iospara) pos
end if
pos = box_bd(2*i) - pos
else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then
!Now extract the number we are reducing from infinity
if(index(pos_string,'inf') < index(pos_string,'+')) then
read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
else
read(pos_string(1:index(pos_string,'+')-1), *, iostat=iospara) pos
end if
pos = box_bd(2*i-1) + pos
else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then
!Now extract the number we are reducing from infinity
if(index(pos_string,'inf') < index(pos_string,'*')) then
read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos
else
read(pos_string(1:index(pos_string,'*')-1), *, iostat=iospara) pos
end if
pos = (box_bd(2*i)-box_bd(2*i-1))*pos + box_bd(2*i-1)
else
read(pos_string, *, iostat=iospara) pos
end if
if (iospara > 0) then
print *, "Error reading position argument ", trim(adjustl(pos_string)), ". Please reformat and try again."
end if
end subroutine parse_pos
subroutine offset_pos(ie, iface, pos)
!This returns a position slightly offset from the center of an element face
!This is used to return a random position within an element discontinuity
integer, intent(in) :: ie !Element index
integer, intent(in) :: iface !Face number between 1 and 6
real(kind=dp), dimension(3), intent(out) :: pos !Position vector
!Other variables we need
integer :: esize
real(kind=dp) :: orient(3,3), ori_inv(3,3), lp, r_cubic_node(3,8)
!First find the offset vector for the face in the untransformed cubic cell
esize = size_ele(ie)
select case(iface)
case(1)
pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp /)
case(2)
pos = (/ real(esize-1,dp)/2.0_dp, -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /)
case(3)
pos = (/ (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
case(4)
pos = (/ real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp /)
case(5)
pos = (/ -10.0_dp**-2.0_dp, real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp /)
case(6)
pos = (/ real(esize-1,dp)/2.0_dp, real(esize-1,dp)/2.0_dp, (esize-1)+10.0_dp**-2.0_dp /)
end select
!Now transform it to real space and adjust it to the position of the element in the first node.
select case (trim(adjustl(type_ele(ie))))
case('fcc')
!First we have to extract the element lattice parameter
call matrix_inverse(sub_box_ori(:,:,sbox_ele(ie)),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 = pos + r_node(:,1, 1, ie)
case default
print *, trim(adjustl(type_ele(ie))), " is not a currently accepted element type for random element positions"
stop 3
end select
end subroutine
end module elements

86
src/io.f90
View File

@ -8,7 +8,7 @@ module io
implicit none
integer :: outfilenum = 0, infilenum = 0
character(len=100) :: outfiles(10), infiles(10)
character(len=100) :: outfiles(100), infiles(100)
logical :: force_overwrite
public
@ -87,6 +87,9 @@ module io
call set_max_esize
do i = 1, outfilenum
print *, "Writing data out to ", trim(adjustl(outfiles(i)))
!Pull out the extension of the file and call the correct write subroutine
select case(trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))))
case('xyz')
@ -130,14 +133,14 @@ module io
do i = 1, ele_num
do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
write(11, '(a, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
write(11, '(i16, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
end do
end do
end do
!Write atom positions
do i = 1, atom_num
write(11, '(a, 3f23.15)') type_atom(i), r_atom(:,i)
write(11, '(i16, 3f23.15)') type_atom(i), r_atom(:,i)
end do
!Finish writing
@ -307,6 +310,9 @@ module io
20 format('SCALARS lattice_type float', /&
'LOOKUP_TABLE default')
21 format('SCALARS esize float', /&
'LOOKUP_TABLE default')
!First we write the vtk file containing the atoms
open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
@ -355,6 +361,10 @@ module io
do i = 1, ele_num
write(11, '(i16)') lat_ele(i)
end do
write(11,21)
do i = 1, ele_num
write(11, '(i16)') size_ele(i)
end do
close(11)
end subroutine
@ -363,7 +373,8 @@ module io
!NOTE: This code doesn't work for arbitrary number of basis atoms per node. It assumes that the
!each element has only 1 atom type at the node.
character(len=100), intent(in) :: file
integer :: interp_max, i, j, lat_size, inod, ibasis, ip, unique_index(10), unique_num
integer :: interp_max, i, j, inod, ibasis, ip, unique_index(50), unique_size(50), unique_num, &
etype
real(kind=dp) :: box_vec(3)
1 format('time' / i16, f23.15)
@ -395,6 +406,22 @@ module io
!Below writes the header information for the restart file
!First figure out all of the unique element types
unique_num = 0
unique_index(:) = 0
eleloop:do i = 1, ele_num
do j =1 , unique_num
if ( ( size_ele(i) == size_ele( unique_index(j) ) ).and. &
( lat_ele(i) == lat_ele(unique_index(j)) ) ) then
cycle eleloop
end if
end do
unique_num = unique_num + 1
unique_index(unique_num) = i
unique_size(unique_num) = size_ele(i)
end do eleloop
!Calculate the max number of atoms per element
select case(max_ng_node)
case(8)
@ -405,7 +432,7 @@ module io
write(11,20) interp_max
write(11,3) node_num
write(11,19) max_ng_node
write(11,4) lattice_types
write(11,4) unique_num
write(11,5) atom_num
write(11,6) 1 !Grain_num is ignored
write(11,16) lattice_types, atom_types
@ -444,26 +471,21 @@ module io
!write the element information
if(ele_num > 0) then
write(11,12)
!First figure out all of the unique element types
unique_num = 0
unique_index(:) = 0
eleloop:do i = 1, ele_num
do j =1 , unique_num
if ( ( size_ele(i) == size_ele( unique_index(j) ) ).and. &
( lat_ele(i) == lat_ele(unique_index(j)) ) ) then
cycle eleloop
end if
end do
unique_num = unique_num + 1
unique_index(unique_num) = i
end do eleloop
do i = 1, unique_num
write(11,'(3i16)') i, size_ele(i)-1, basis_type(1,i)
write(11,'(3i16)') i, size_ele(unique_index(i))-1, basis_type(1,lat_ele(unique_index(i)))
end do
ip = 0
write(11,13)
do i = 1, ele_num
write(11, '(4i16)') i, lat_ele(i), 1, basis_type(1,lat_ele(i))
!Figure out the ele type
do j = 1, unique_num
if ( unique_size(j) == size_ele(i)) then
etype = j
exit
endif
end do
write(11, '(4i16)') i, etype, 1, basis_type(1,lat_ele(i))
do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
ip = ip + 1
@ -520,7 +542,7 @@ module io
!Write out atoms first
do i = 1, atom_num
write(11,*) i, type_atom(i), r_atom(:,i)
write(11,*) i, type_atom(i), sbox_atom(i), r_atom(:,i)
end do
!Write out the elements, this is written in two stages, one line for the element and then 1 line for
@ -621,12 +643,18 @@ module io
!Read in the box boundary and grow the current active box bd
read(11, *) temp_box_bd(:)
print *, "displace", displace
do i = 1, 3
newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
if (abs(displace(i)) > lim_zero) then
newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
else
newdisplace(i)=displace(i)
end if
temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i)
temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i)
end do
print *, "newdisplace", newdisplace
call grow_box(temp_box_bd)
!Read in the number of sub_boxes and allocate the variables
read(11, *) n
@ -713,8 +741,9 @@ module io
!Read the atoms
do i = 1, in_atoms
read(11,*) j, type, r(:)
call add_atom(new_type_to_type(type), r+newdisplace)
read(11,*) j, type, sbox, r(:)
r = r+newdisplace
call add_atom(new_type_to_type(type), sbox+sub_box_num, r)
end do
!Read the elements
@ -726,7 +755,7 @@ module io
r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace
end do
end do
call add_element(etype, size, new_lattice_map(type), sbox+n, r_innode)
call add_element(etype, size, new_lattice_map(type), sbox+sub_box_num, r_innode)
end do
!Close the file being read
@ -734,10 +763,11 @@ module io
!Only increment the lattice types if there are elements, if there are no elements then we
!just overwrite the arrays
if(in_eles > 0) lattice_types = maxval(new_lattice_map)
lattice_types = maxval(new_lattice_map)
sub_box_num = sub_box_num + n
call set_max_esize
end subroutine read_mb
end module io

19
src/main.f90
View File

@ -37,21 +37,29 @@ program main
!Call initialization functions
call lattice_init
call box_init
call random_seed
force_overwrite=.false.
wrap_flag = .false.
end_mode_arg = 0
! Command line parsing
arg_num = command_argument_count()
!Check to see if overwrite flag is passed
!First check to see if certain commands are passed, these commands must be known before code
!is executed.
do i = 1, arg_num
call get_command_argument(i,argument)
select case(trim(adjustl(argument)))
!This lets us know if we are overwriting all files
case('-ow')
force_overwrite = .true.
print *, "Overwrite flag passed, output files will be overwritten"
!This lets us know if we need to wrap atomic positions back into the cell
case('-wrap')
wrap_flag=.true.
end select
end do
!Determine if a mode is being used and what it is. The first argument has to be the mode
@ -85,11 +93,14 @@ program main
call call_option(argument, arg_pos)
!Otherwise print that the argument is not accepted and move on
else
print *, trim(adjustl(argument)), " is not accepted. Skipping to next argument"
arg_pos = arg_pos + 1
print *, trim(adjustl(argument)), " is not an accepted command."
stop 3
end if
end do
!If wrap flag was passed then call the wrap atoms command
if(wrap_flag) call wrap_atoms
!Check to make sure a file was passed to be written out and then write out
! Before building do a check on the file
if (outfilenum == 0) then
@ -98,4 +109,4 @@ program main
end if
call write_out
end program main
end program main

2
src/mode_convert.f90
View File

@ -11,7 +11,7 @@ module mode_convert
subroutine convert(arg_pos)
!This subroutine converts a single input file from one format to another
integer, intent(out) :: arg_pos
character(len=100) :: infile, outfile
character(len=100) :: infile
real(kind = dp) :: temp_box_bd(6)
!First read in the file
call get_command_argument(2, infile)

53
src/mode_create.f90
View File

@ -24,7 +24,6 @@ module mode_create
subroutine create(arg_pos)
! Main subroutine which controls execution
character(len=100) :: textholder
integer, intent(out) :: arg_pos
integer :: i, ibasis, inod
@ -106,8 +105,8 @@ module mode_create
end do
end do
do i = 1,3
box_bd(2*i) = maxval(r_node_temp(i,:,:))
box_bd(2*i-1) = origin(i)
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(element_type, esize, 1, 1, r_node_temp)
end if
@ -128,14 +127,15 @@ module mode_create
!Now that it is built multiply by the lattice parameter
box_bd = box_bd*lattice_parameter
print *, "Using mode create, ", lat_ele_num, " elements are created and ", lat_atom_num*basisnum(1), " atoms are created."
print *, "Using mode create, ", lat_ele_num, " elements are created and ", lat_atom_num*basisnum(1), &
" atoms are created."
!Allocate variables
call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1))
if(lat_atom_num > 0) then
do i = 1, lat_atom_num
do ibasis = 1, basisnum(1)
call add_atom(basis_type(ibasis, 1), (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
call add_atom(basis_type(ibasis, 1), 1, (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
end do
end do
deallocate(r_atom_lat)
@ -169,6 +169,7 @@ module mode_create
integer :: ori_pos, i, j, arglen, stat
character(len=100) :: textholder
character(len=8) :: orient_string
logical :: isortho, isrighthanded
!Pull out all required positional arguments
@ -205,25 +206,25 @@ module mode_create
do i = 1, 3
call get_command_argument(arg_pos, orient_string, arglen)
if (arglen==0) STOP "Missing orientation in orient command of mode create"
call parse_ori_vec(orient_string, orient(i,:))
arg_pos = arg_pos+1
ori_pos=2
do j = 1,3
if (orient_string(ori_pos:ori_pos) == '-') then
ori_pos = ori_pos + 1
read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
if (stat>0) STOP "Error reading orient value"
orient(i,j) = -orient(i,j)
ori_pos = ori_pos + 1
else
read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
if(stat>0) STOP "Error reading orient value"
ori_pos=ori_pos + 1
end if
end do
! ori_pos=2
! do j = 1,3
! if (orient_string(ori_pos:ori_pos) == '-') then
! ori_pos = ori_pos + 1
! read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
! if (stat>0) STOP "Error reading orient value"
! orient(i,j) = -orient(i,j)
! ori_pos = ori_pos + 1
! else
! read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
! if(stat>0) STOP "Error reading orient value"
! ori_pos=ori_pos + 1
! end if
! end do
end do
!If the duplicate command is passed then we extract the information on the new bounds.
case('duplicate')
if(dim_flag) STOP "Both duplicate and dim options cannot be used in mode_create"
@ -284,6 +285,13 @@ module mode_create
end select
!Now normalize the orientation matrix
orient = matrix_normal(orient,3)
!Now check these to make sure they are right handed and orthogonal
call check_right_ortho(orient, isortho, isrighthanded)
if (.not.isortho) then
stop "Directions in orient are not orthogonal"
else if (.not.isrighthanded) then
stop "Directions in orient are not righthanded"
end if
!Set lattice_num to 1 and add the lattice_parameter to the elements module lattice paramter variable
lattice_types = 1
@ -309,7 +317,6 @@ module mode_create
integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,8), rzero(3), &
vlat(3), temp_lat(3,8), m, n, o
real(kind=dp) :: v(3), temp_nodes(3,1,8)
real(kind=dp), allocatable :: resize_lat_array(:,:)
logical, allocatable :: lat_points(:,:,:)
logical :: node_in_bd(8)
@ -386,7 +393,7 @@ module mode_create
outerloop: do iz = 1, bd_in_array(3)
do iy = 1, bd_in_array(2)
do ix = 1, bd_in_array(1)
node_in_bd(8) = .false.
node_in_bd(:) = .false.
do inod = 1, 8
vlat = ele(:,inod) + (/ ix, iy, iz /)

32
src/mode_merge.f90
View File

@ -10,7 +10,7 @@ module mode_merge
character(len=4) :: dim
integer :: in_num, new_starts(2)
real(kind=dp) :: shift_vec(3)
logical :: wrap, shift_flag
logical :: shift_flag
public
contains
@ -22,7 +22,6 @@ module mode_merge
print *, '-----------------------Mode Merge---------------------------'
wrap = .false.
shift_flag = .false.
shift_vec(:) = 0.0_dp
@ -115,8 +114,6 @@ module mode_merge
if (arglen==0) stop "Missing vector component for shift command"
read(textholder, *) shift_vec(i)
end do
case('wrap')
wrap = .true.
case default
!If it isn't an available option to mode merge then we just exit
exit
@ -126,16 +123,16 @@ module mode_merge
end subroutine parse_command
subroutine shift(array_start, filenum)
!This subroutine applies a shift to newly added atoms and elements. It also wraps the atoms
!if the user provides the wrap flag
!This subroutine applies a shift to newly added atoms and elements.
integer, dimension(2), intent(in) :: array_start
integer, intent(in) :: filenum
integer :: i, j, ibasis, inod
integer :: i, ibasis, inod
real(kind=dp), dimension(3) :: current_shift
character(len=3) :: alldims
alldims = 'xyz'
!Calculate the current shift which is the filenum-1 multiplied by the user specified shift
current_shift = (filenum-1)*shift_vec
@ -155,19 +152,16 @@ module mode_merge
end do
end do
!Now we check if we have to wrap the atoms, nodes are not wrapped. For elements the periodic
!boundary conditions are applied in the actual CAC codes
if(wrap) then
do i = array_start(1), atom_num
call apply_periodic(r_atom(:,i))
end do
!If we don't include the wrap command then we have to increase the size of the box
else
if(.not.(wrap_flag)) then
do i = 1,3
if (current_shift(i) < -lim_zero) then
box_bd(2*i-1) = box_bd(2*i-1) - current_shift(i)
else if (current_shift(i) > lim_zero) then
if (alldims(i:i) /= trim(adjustl(dim))) then
if (current_shift(i) < -lim_zero) then
box_bd(2*i-1) = box_bd(2*i-1) - current_shift(i)
else if (current_shift(i) > lim_zero) then
box_bd(2*i) = box_bd(2*i) + current_shift(i)
end if
else if (alldims(i:i) == trim(adjustl(dim))) then
box_bd(2*i) = box_bd(2*i) + current_shift(i)
end if
end do

125
src/opt_delete.f90 Normal file
View File

@ -0,0 +1,125 @@
module opt_delete
use parameters
use subroutines
use elements
implicit none
real(kind=dp) :: rc_off
public
contains
subroutine run_delete(arg_pos)
integer, intent(inout) :: arg_pos
rc_off = -lim_zero
!Main calling function for delete option.
print *, '-----------------------Option Delete------------------------'
call parse_delete(arg_pos)
if (rc_off > 0.0_dp) call delete_overlap
end subroutine run_delete
subroutine parse_delete(arg_pos)
!Parse the delete command
integer, intent(inout) :: arg_pos
integer :: arg_len
character(len=100) :: textholder
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arg_len)
if(arg_len==0) stop "Missing argument to delete command"
select case(textholder)
case('overlap')
arg_pos=arg_pos + 1
call get_command_argument(arg_pos, textholder, arg_len)
if(arg_len==0) stop "Missing argument to delete overlap command"
read(textholder, *) rc_off
case default
print *, "Command ", trim(adjustl(textholder)), " is not accepted for option delete"
stop 3
end select
arg_pos = arg_pos + 1
end subroutine parse_delete
subroutine delete_overlap
!This subroutine deletes all overlapping atoms, which is defined as atoms which are separated by a distance of
!less then rc_off
integer :: i, c(3), ci, cj, ck, num_nei, nei, delete_num
integer, dimension(atom_num) :: for_delete
!These are the variables containing the cell list information
integer, dimension(3) :: cell_num
integer, allocatable :: num_in_cell(:,:,:), which_cell(:,:)
integer, allocatable :: cell_list(:,:,:,:)
allocate(which_cell(3,atom_num))
!First pass the atom list and atom num to the algorithm which builds the cell list
call build_cell_list(atom_num, r_atom, rc_off, cell_num, num_in_cell, cell_list, which_cell)
!Now loop over every atom and figure out if it has neighbors within the rc_off
delete_num = 0
atom_loop: do i = 1, atom_num
!c is the position of the cell that the atom belongs to
c = which_cell(:,i)
!Check to make sure it hasn't already been deleted
if(all(c /= 0)) then
!Now loop over all neighboring cells
do ci = -1, 1, 1
do cj = -1, 1, 1
do ck = -1, 1, 1
if (any((c + (/ ck, cj, ci /)) == 0)) cycle
if( (c(1) + ck > cell_num(1)).or.(c(2) + cj > cell_num(2)).or. &
(c(3) + ci > cell_num(3))) cycle
do num_nei = 1, num_in_cell(c(1) + ck, c(2) + cj, c(3) + ci)
nei = cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci)
!Check to make sure the atom isn't the same index as the atom we are checking
!and that the neighbor hasn't already been deleted
if((nei /= i).and.(nei/= 0)) then
!Now check to see if it is in the cutoff radius, if it is add it to the delete code
if (norm2(r_atom(:,nei)-r_atom(:,i)) < rc_off) then
delete_num = delete_num + 1
for_delete(delete_num) = max(i,nei)
!Now zero out the larger index
if(i > nei) then
which_cell(:,i) = 0
cycle atom_loop
else
which_cell(:,nei) = 0
cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci) = 0
end if
end if
end if
end do
end do
end do
end do
end if
end do atom_loop
print *, "Overlap command deletes ", delete_num, " atoms"
!Now delete all the atoms
call delete_atoms(delete_num, for_delete(1:delete_num))
end subroutine delete_overlap
end module opt_delete

35
src/opt_disl.f90
View File

@ -47,8 +47,8 @@ module opt_disl
integer, intent(inout) :: arg_pos
integer :: i,arglen
character(len=8) :: ori_string
character(len=100) :: textholder
character(len=8) :: ori_string
!Parse all of the commands
arg_pos = arg_pos + 1
@ -180,6 +180,9 @@ module opt_disl
end do
end if
!Now make sure all atoms are wrapped back into the simulation cell
call wrap_atoms
end subroutine dislgen
subroutine parse_disloop(arg_pos)
@ -187,8 +190,7 @@ module opt_disl
integer, intent(inout) :: arg_pos
integer :: i,arglen, sbox
character(len=8) :: ori_string
integer :: i,arglen
character(len=100) :: textholder
!Parse all of the commands
@ -283,18 +285,18 @@ module opt_disl
if(loop_radius < 0.0_dp) then
ALLOCATE(xLoop(4,3))
xLoop(:,:) = 0.d0
xLoop(1,a1) = centroid(1) - loop_radius
xLoop(1,a2) = centroid(2) - loop_radius
xLoop(1,a3) = centroid(3)
xLoop(2,a1) = centroid(1) + loop_radius
xLoop(2,a2) = centroid(2) - loop_radius
xLoop(2,a3) = centroid(3)
xLoop(3,a1) = centroid(1) + loop_radius
xLoop(3,a2) = centroid(2) + loop_radius
xLoop(3,a3) = centroid(3)
xLoop(4,a1) = centroid(1) - loop_radius
xLoop(4,a2) = centroid(2) + loop_radius
xLoop(4,a3) = centroid(3)
xLoop(1,a1) = centroid(a1) + loop_radius
xLoop(1,a2) = centroid(a2) + loop_radius
xLoop(1,a3) = centroid(a3)
xLoop(2,a1) = centroid(a1) - loop_radius
xLoop(2,a2) = centroid(a2) + loop_radius
xLoop(2,a3) = centroid(a3)
xLoop(3,a1) = centroid(a1) - loop_radius
xLoop(3,a2) = centroid(a2) - loop_radius
xLoop(3,a3) = centroid(a3)
xLoop(4,a1) = centroid(a1) + loop_radius
xLoop(4,a2) = centroid(a2) - loop_radius
xLoop(4,a3) = centroid(a3)
else
!Calculate loop perimeter
perimeter = 2.0_dp*pi*loop_radius
@ -392,6 +394,9 @@ module opt_disl
end do
end do
return
!Now make sure all atoms are wrapped back into the simulation cell
call wrap_atoms
end subroutine
!********************************************************

268
src/opt_group.f90
View File

@ -8,10 +8,10 @@ module opt_group
use box
implicit none
integer :: group_ele_num, group_atom_num, remesh_size
integer :: group_ele_num, group_atom_num, remesh_size, remesh_type
character(len=15) :: type, shape !Type indicates what element type is selected and shape is the group shape
real(kind=dp) :: block_bd(6), disp_vec(3)
logical :: displace, wrap
real(kind=dp) :: block_bd(6), disp_vec(3), remesh_lat_pam
logical :: displace, delete, max_remesh, refine
integer, allocatable :: element_index(:), atom_index(:)
@ -27,6 +27,11 @@ module opt_group
group_ele_num = 0
group_atom_num = 0
remesh_size=0
displace=.false.
delete=.false.
max_remesh=.false.
refine = .false.
if(allocated(element_index)) deallocate(element_index)
if(allocated(atom_index)) deallocate(atom_index)
@ -38,6 +43,11 @@ module opt_group
if(displace) call displace_group
if(remesh_size > 0) call remesh_group
if(delete) call delete_group
if(refine) call refine_group
end subroutine group
subroutine parse_group(arg_pos)
@ -96,13 +106,25 @@ module opt_group
if (arglen==0) stop "Missing vector component for shift command"
read(textholder, *) disp_vec(i)
end do
case('wrap')
wrap = .true.
case('refine')
refine=.true.
case('remesh')
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) stop "Missing remesh element size in group command"
read(textholder, *) remesh_size
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) stop "Missing remesh lattice parameter in group command"
read(textholder, *) remesh_lat_pam
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) stop "Missing remesh type in group command"
read(textholder, *) remesh_type
case('max')
max_remesh =.true.
case('delete')
delete=.true.
case default
!If it isn't an available option to opt_disl then we just exit
exit
@ -192,16 +214,8 @@ module opt_group
end do
end do
!Now either apply periodic boundaries if wrap command was passed or adjust box dimensions
!Now we check if we have to wrap the atoms, nodes are not wrapped. For elements the periodic
!boundary conditions are applied in the actual CAC codes
if(wrap) then
do i = 1, atom_num
call apply_periodic(r_atom(:,i))
end do
!If we don't include the wrap command then we have to increase the size of the box
else
if (.not.(wrap_flag)) then
do i = 1,3
if (disp_vec(i) < -lim_zero) then
box_bd(2*i-1) = box_bd(2*i-1) - disp_vec(i)
@ -213,47 +227,217 @@ module opt_group
end subroutine displace_group
subroutine remesh_group
subroutine refine_group
!This command is used to remesh the group to a desired element size
integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3)
!Refining to atoms and remeshing to elements are different processes so check which code we need to run
select case(remesh_size)
!Refining to atoms
case(2)
if(group_ele_num > 0) then
orig_atom_num = atom_num
!Estimate number of atoms we are adding, this doesn't have to be exact
add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3
call grow_ele_arrays(0,add_atom_num)
do i = 1, group_ele_num
ie = element_index(i)
!Get the interpolated atom positions
call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
if(group_ele_num > 0) then
orig_atom_num = atom_num
!Estimate number of atoms we are adding, this doesn't have to be exact
add_atom_num = group_ele_num*basisnum(lat_ele(element_index(1)))*size_ele(element_index(1))**3
call grow_ele_arrays(0,add_atom_num)
do i = 1, group_ele_num
ie = element_index(i)
!Get the interpolated atom positions
call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
!Loop over all interpolated atoms and add them to the system, we apply periodic boundaries
!here as well to make sure they are in the box
do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3
call apply_periodic(r_interp(:,j))
call add_atom(type_interp(j), sbox_ele(ie), r_interp(:,j))
end do
end do
!Once all atoms are added we delete all of the elements
call delete_elements(group_ele_num, element_index)
print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms."
end if
!Loop over all interpolated atoms and add them to the system, we apply periodic boundaries here as well to make sure
!they are in the box
do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3
call apply_periodic(r_interp(:,j))
call add_atom(type_interp(j), r_interp(:,j))
end do
end subroutine
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), old_ele, old_atom, &
max_loops, working_esize
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)
logical, allocatable :: lat_points(:,:,:)
character(len=100) :: remesh_ele_type
!Right now we just hardcode only remeshing to elements
remesh_ele_type = 'fcc'
!Get the orientations, this assumes that the orientation of the subbox for the first atom is the
!same as the rest of the atoms
!If this assumption is false then the code will break and exit
orient = sub_box_ori(:, :, sbox_atom(atom_index(1)))
call matrix_inverse(orient,3,ori_inv)
!First calculate max position in lattice space to be able to allocate lat_points array, also sum the total numbers of
!degrees of freedom which are added
dof = 0
select case(trim(adjustl(shape)))
case('block')
group_in_lat = reshape((/ block_bd(1),block_bd(3),block_bd(5), &
block_bd(2),block_bd(3),block_bd(5), &
block_bd(2),block_bd(4),block_bd(5), &
block_bd(1),block_bd(4),block_bd(5), &
block_bd(1),block_bd(3),block_bd(6), &
block_bd(2),block_bd(3),block_bd(6), &
block_bd(2),block_bd(4),block_bd(6), &
block_bd(1),block_bd(4),block_bd(6) /), [3,8])
group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/remesh_lat_pam))
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
!Once all atoms are added we delete all of the elements
call delete_elements(group_ele_num, element_index)
end select
print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms."
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
end if
!Remeshing to elements, currently not available
case default
print *, "Remeshing to elements is currently not available. Please refine to atoms by passing a remsh size of 2"
end select
!Now place all group atoms and group interpolated atoms into lat_points
do i = 1, group_atom_num
r = r_atom(:,atom_index(i))/remesh_lat_pam
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
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)/remesh_lat_pam
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
print *, "Group has ", dof, " degrees of freedom to remesh"
!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)
old_atom = atom_num
old_ele = ele_num
!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-2)/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(remesh_type)
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))*remesh_lat_pam
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.
call add_element(remesh_ele_type, working_esize, remesh_type, 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))*remesh_lat_pam
call add_atom(remesh_type, sbox_atom(atom_index(1)), r)
end if
end do
end do
end do
print *, "Remeshing has created ", ele_num-old_ele, " elements and ", atom_num-old_atom, " atoms."
end subroutine remesh_group
subroutine delete_group
!This subroutine deletes all atoms/elements within a group
print *, "Deleting group containing ", group_atom_num, " atoms and ", group_ele_num, " elements."
!Delete atoms
call delete_atoms(group_atom_num, atom_index)
!Delete elements
call delete_elements(group_ele_num, element_index)
end subroutine delete_group
end module opt_group

132
src/opt_orient.f90 Normal file
View File

@ -0,0 +1,132 @@
module opt_orient
!This module contains the orient option which allows for the reorientation
!of simulation cells. This can be used to create arbitrarily oriented dislocation or loops.
use parameters
use subroutines
use elements
use box
implicit none
real(kind=dp), save :: new_orient(3,3)
real(kind=dp), dimension(6) :: orig_box_bd
real(kind=dp), allocatable :: orig_sub_box_ori(:,:,:)
public
contains
subroutine orient(arg_pos)
integer, intent(inout) :: arg_pos
integer :: i, ibasis, inod
logical :: isortho, isrighthanded
real(kind=dp) :: inv_sub_box_ori(3,3,sub_box_num)
!First parse the orient command
call parse_orient(arg_pos)
!Now rotate the basis. To do this we transform the basis to [100] [010] [001] and then
!transform to user specified basis.
!Find all inverse orientation matrices for all sub_boxes
do i = 1, sub_box_num
call matrix_inverse(sub_box_ori, 3, inv_sub_box_ori)
end do
!Now transform all atoms
do i = 1, atom_num
r_atom(:,i) = matmul(new_orient,matmul(inv_sub_box_ori(:,:,sbox_atom(i)),r_atom(:,i)))
end do
!Now transform all elements
do i = 1, ele_num
do inod =1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
r_node(:,ibasis,inod,i) = matmul(new_orient,matmul(inv_sub_box_ori(:,:,sbox_ele(i)),r_node(:,ibasis,inod,i)))
end do
end do
end do
!Now save the original sub_box_ori and overwrite them
if(allocated(orig_sub_box_ori)) deallocate(orig_sub_box_ori)
allocate(orig_sub_box_ori(3,3,sub_box_num))
orig_sub_box_ori = sub_box_ori
!Now overwrite the orientations
do i = 1, sub_box_num
sub_box_ori(:,:,i) = new_orient
end do
!Save original box boundaries
orig_box_bd = box_bd
!Now find new box boundaries
call def_new_box
end subroutine orient
subroutine parse_orient(arg_pos)
!This parses the orient option
integer, intent(inout) :: arg_pos
integer :: i, arg_len
logical :: isortho, isrighthanded
character(len=8) :: ori_string
!Pull out the new user orientation
do i = 1, 3
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, ori_string, arg_len)
if (arg_len == 0) print *, "Missing orientation vector in -orient option"
call parse_ori_vec(ori_string, new_orient(i,:))
end do
!Normalize the orientation matrix
new_orient = matrix_normal(new_orient,3)
!Check right hand rule and orthogonality
call check_right_ortho(new_orient, isortho, isrighthanded)
if (.not.isortho) then
stop "Directions in orient are not orthogonal"
else if (.not.isrighthanded) then
stop "Directions in orient are not righthanded"
end if
arg_pos = arg_pos + 1
end subroutine parse_orient
subroutine unorient
integer :: i, ibasis, inod
real(kind=dp) :: inv_ori(3,3)
!Now rotate the basis. To do this we transform the basis to [100] [010] [001] and then
!transform to the original sbox_ele
!Find the inverse for the new orientation matrix
call matrix_inverse(new_orient, 3, inv_ori)
!Recover original sub_box_ori
sub_box_ori = orig_sub_box_ori
!Now transform all atoms
do i = 1, atom_num
r_atom(:,i) = matmul(sub_box_ori(:,:,sbox_atom(i)),matmul(inv_ori(:,:),r_atom(:,i)))
end do
!Now transform all elements
do i = 1, ele_num
do inod =1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
r_node(:,ibasis,inod,i) = matmul(sub_box_ori(:,:,sbox_ele(i)),matmul(inv_ori,r_node(:,ibasis,inod,i)))
end do
end do
end do
!Restore original box boundaries
box_bd = orig_box_bd
end subroutine unorient
end module opt_orient

2
src/parameters.f90
View File

@ -7,7 +7,7 @@ module parameters
!Parameters for floating point tolerance
real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), &
lim_large = huge(1.0_dp)
logical, save :: lmpcac
logical, save :: lmpcac
!Numeric constants
real(kind=dp), parameter :: pi = 3.14159265358979323846_dp

151
src/subroutines.f90
View File

@ -5,6 +5,7 @@ module subroutines
implicit none
integer :: allostat, deallostat
public
contains
@ -170,52 +171,7 @@ module subroutines
return
end subroutine parse_ori_vec
subroutine parse_pos(i, pos_string, pos)
!This subroutine parses the pos command allowing for command which include inf
integer, intent(in) :: i !The dimension of the position
character(len=100), intent(in) :: pos_string !The position string
real(kind=dp), intent(out) :: pos !The output parsed position value
integer :: iospara
iospara = 0
if(trim(adjustl(pos_string)) == 'inf') then
pos=box_bd(2*i)
else if(trim(adjustl(pos_string)) == '-inf') then
pos=box_bd(2*i-1)
else if ((index(pos_string,'-') > 0).and.(index(pos_string,'inf')>0)) then
!Now extract the number we are reducing from infinity
if(index(pos_string,'inf') < index(pos_string,'-')) then
read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
else
read(pos_string(1:index(pos_string,'-')-1), *, iostat=iospara) pos
end if
pos = box_bd(2*i) - pos
else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then
!Now extract the number we are reducing from infinity
if(index(pos_string,'inf') < index(pos_string,'+')) then
read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
else
read(pos_string(1:index(pos_string,'+')-1), *, iostat=iospara) pos
end if
pos = box_bd(2*i-1) + pos
else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then
!Now extract the number we are reducing from infinity
if(index(pos_string,'inf') < index(pos_string,'*')) then
read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos
else
read(pos_string(1:index(pos_string,'*')-1), *, iostat=iospara) pos
end if
pos = (box_bd(2*i)-box_bd(2*i-1))*pos + box_bd(2*i-1)
else
read(pos_string, *, iostat=iospara) pos
end if
if (iospara > 0) then
print *, "Error reading position argument ", trim(adjustl(pos_string)), ". Please reformat and try again."
end if
end subroutine parse_pos
subroutine heapsort(a)
@ -283,4 +239,107 @@ module subroutines
end do
end subroutine
subroutine build_cell_list(numinlist, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell)
!This subroutine builds a cell list based on rc_off
!----------------------------------------Input/output variables-------------------------------------------
integer, intent(in) :: numinlist !The number of points within r_list
real(kind=dp), dimension(3,numinlist), intent(in) :: r_list !List of points to be used for the construction of
!the cell list.
real(kind=dp), intent(in) :: rc_off ! Cutoff radius which dictates the size of the cells
integer, dimension(3), intent(inout) :: cell_num !Number of cells in each dimension.
integer, allocatable, intent(inout) :: num_in_cell(:,:,:) !Number of points within each cell
integer, allocatable, intent(inout) :: cell_list(:,:,:,:) !Index of points from r_list within each cell.
integer, dimension(3,numinlist), intent(out) :: which_cell !The cell index for each point in r_list
!----------------------------------------Begin Subroutine -------------------------------------------
integer :: i, j, cell_lim, c(3)
real(kind=dp) :: box_len(3)
integer, allocatable :: resize_cell_list(:,:,:,:)
!First calculate the number of cells that we need in each dimension
do i = 1,3
box_len(i) = box_bd(2*i) - box_bd(2*i-1)
cell_num(i) = int(box_len(i)/(rc_off/2))+1
end do
!Initialize/allocate variables
cell_lim = 10
allocate(num_in_cell(cell_num(1),cell_num(2),cell_num(3)), cell_list(cell_lim, cell_num(1), cell_num(2), cell_num(3)))
!Now place points within cell
do i = 1, numinlist
!c is the position of the cell that the point belongs to
do j = 1, 3
c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off/2)) + 1
end do
!Place the index in the correct position, growing if necessary
num_in_cell(c(1),c(2),c(3)) = num_in_cell(c(1),c(2),c(3)) + 1
if (num_in_cell(c(1),c(2),c(3)) > cell_lim) then
allocate(resize_cell_list(cell_lim+10,cell_num(1),cell_num(2),cell_num(3)))
resize_cell_list(1:cell_lim, :, :, :) = cell_list
resize_cell_list(cell_lim+1:, :, :, :) = 0
call move_alloc(resize_cell_list, cell_list)
end if
cell_list(num_in_cell(c(1),c(2),c(3)),c(1),c(2),c(3)) = i
which_cell(:,i) = c
end do
return
end subroutine build_cell_list
subroutine check_right_ortho(ori, isortho, isrighthanded)
!This subroutine checks whether provided orientations in the form:
! | x1 x2 x3 |
! | y1 y2 y3 |
! | z1 z2 z3 |
!are right handed
real(kind=dp), dimension(3,3), intent(in) :: ori
logical, intent(out) :: isortho, isrighthanded
integer :: i, j
real(kind=dp) :: v(3), v_k(3)
!Initialize variables
isortho = .true.
isrighthanded=.true.
do i = 1, 3
do j = i+1, 3
if(abs(dot_product(ori(i,:), ori(j,:))) > lim_zero) then
isortho = .false.
end if
!Check if they are righthanded
if (j == i+1) then
v(:) = cross_product(ori(i,:), ori(j,:))
v_k(:) = v(:) - ori(mod(j, 3)+1,:)
else if ((i==1).and.(j==3)) then
v(:) = cross_product(ori(j,:),ori(i,:))
v_k(:) = v(:) - ori(i+1, :)
end if
if(norm2(v_k) > 10.0_dp**(-8.0_dp)) then
isrighthanded=.false.
end if
end do
end do
return
end subroutine check_right_ortho
end module subroutines