Merge pull request #18 from aselimov/development

alpha v1.2
master
aselimov 5 years ago committed by GitHub
commit 21be35041c
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23

@ -0,0 +1,16 @@
# Install CACmb
In order to build CACmb switch to the src repository and run:
```
make
```
Once this command is complete you can run:
```
sudo make install
```
in order to copy the installation to `/usr/local/bin`.

@ -113,6 +113,26 @@ This mode merges multiple data files and creates one big simulation cell. The pa
`dim` - the dimension they are to be stacked along, can be either `x`, `y`, or `z`. If the argument `none` is passed then the cells are just overlaid. Future options will include a delete overlap command.
**Additional options:**
**Shift**
```
shift x y z
```
If the shift command is passed to mode merge then each file after the first file in the merge command is displaced by the vector `[x, y, z]`. This is additive so if you are merging three files and this command is passed then the second file is shifted by `[x,y,z]` and the third file is shifted by `2*[x,y,z]`.
Example: `cacmb --merge z 2 Cu.mb Cu2.mb Cu3.mb Cumerged.mb shift 2 0 0` will shift the atomic and element positions in the `Cu2.mb` file by `[2,0,0]` and the positions in `Cu3.mb` by `[4,0,0]`.
**Wrap**
```
wrap
```
This will wrap atomic positions back inside the box. Effectively as if periodic boundary conditions are applied so that atoms which exit from one side of the simulation cell enter back in through the other.
## Options
Options are additional portions of code which have additional functionality. Options are performed in the order that they appear in the argument list and can be added to any mode. If wanting to use strictly options use `--convert` to specify input and output files.
@ -120,7 +140,7 @@ Options are additional portions of code which have additional functionality. Opt
### Option dislgen
```
-dislgen [ijk] [hkl] x y z char_angle poisson
-dislgen [ijk] [hkl] x y z burgers char_angle poisson
```
@ -131,6 +151,8 @@ This options adds an arbitrarily oriented dislocation into your model based on u
`[hkl]` - The vector for the slip plane
`burgers` - The magnitude of the burgers vector for the dislocation to be inserted
`x y z` - The position of the dislocation centroid
`char_angle` - Character angle of the dislocation (0 is screw and 90 is edge)
@ -153,4 +175,54 @@ This option deletes vacancies on a plane which when minimized should result in a
`bx by bz` - The burgers vector for the dislocation
`poisson` - Poisson ratio for continuum solution
`poisson` - Poisson ratio for continuum solution
### Option Group
`-group select_type group_shape shape_arguments additional keywords`
This option selects a group of either elements, nodes, atoms and applies some transformation to them.
`select_type` - Either `nodes`, `atoms`, `elements`, `nodes/atoms`, `all`. When using the option `nodes` only nodes which are within the group are selected, `elements` selects elements based on whether the element center is within the group, `nodes/atoms` selects both nodes and atoms for the group. `all` selects elements based on the element center and atoms based on their position.
`group_shape` - Specifies what shape the group takes and dictates which options must be passed. Each shape requires different arguments and these arguments are represented by the placeholder `shape_arguments`. The accepted group shapes and arguments are below:
*Block:*
`-group nodes block xlo xhi ylo yhi zlo zhi`
This selects a group residing in a block with edges perpendicular to the simulation cell. The block boundaries are given by `xlo xhi ylo yhi zlo zhi`.
`additional keywords`- Represents the various transformations which can be performed on a group. These additional keywords are given below.
**Displace**
```
displace x y z
```
This is similar to the mode merge `-shift` argument. This simply shift atoms and elements within the group by a vector `[x,y,z]`.
**Wrap**
```
wrap
```
This command wraps atoms back into the simulation cell as though periodic boundary conditions are being used.
**Remesh**
```
remesh esize
```
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.
### Option overwrite
```
-ow
```
If this option is passed then all files are automatically overwritten without asking the user.

@ -2,8 +2,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
MODES=mode_create.o mode_merge.o mode_convert.o
OPTIONS=opt_disl.o
OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o call_option.o $(MODES) $(OPTIONS)
OPTIONS=opt_disl.o opt_group.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:
.SUFFIXES: .c .f .f90 .F90 .o
@ -38,6 +38,6 @@ atoms.o subroutines.o testfuncs.o box.o : functions.o
main.o io.o $(MODES) $(OPTIONS) : elements.o
call_mode.o : $(MODES)
call_option.o : $(OPTIONS)
$(MODES) io.o: atoms.o box.o
$(MODES) $(OPTIONS) subroutines.o io.o : atoms.o box.o
$(MODES) main.o : io.o
testfuncs.o elements.o mode_create.o opt_disl.o: subroutines.o
testfuncs.o elements.o mode_create.o $(OPTIONS): subroutines.o

@ -0,0 +1,22 @@
subroutine call_option(option, arg_pos)
use parameters
use opt_disl
use opt_group
implicit none
integer, intent(inout) :: arg_pos
character(len=100), intent(in) :: option
select case(trim(adjustl(option)))
case('-dislgen', '-disloop')
call dislocation(option, arg_pos)
case('-group')
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'
arg_pos = arg_pos + 1
end select
end subroutine call_option

@ -223,6 +223,8 @@ module elements
real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node)
ele_num = ele_num + 1
!Check to see if we need to grow the arrays
call grow_ele_arrays(1,0)
type_ele(ele_num) = type
size_ele(ele_num) = size
lat_ele(ele_num) = lat
@ -239,6 +241,8 @@ module elements
real(kind=dp), intent(in), dimension(3) :: r
atom_num = atom_num+1
!Check to see if we need to grow the arrays
call grow_ele_arrays(0,1)
type_atom(atom_num) = type
r_atom(:,atom_num) = r(:)
@ -356,7 +360,7 @@ module elements
end do
end select
if (ia /= esize**3) then
if (ia /= bnum*esize**3) then
print *, "Incorrect interpolation"
stop 3
end if
@ -403,4 +407,34 @@ module elements
atom_num = atom_num - 1
end do
end subroutine
subroutine delete_elements(num, index)
!This subroutine deletes elements from the element array
integer, intent(in) :: num
integer, intent(inout), dimension(num) :: index
integer :: i, j
call heapsort(index)
!We go from largest index to smallest index just to make sure that we don't miss
!accidentally overwrite values which need to be deleted
do i = num, 1, -1
if(index(i) == ele_num) then
r_node(:,:,:,index(i)) = 0.0_dp
type_ele(index(i)) =''
size_ele(index(i)) = 0
lat_ele(index(i)) = 0
sbox_ele(index(i)) = 0
else
node_num = node_num - ng_node(lat_ele(index(i)))
r_node(:,:,:,index(i)) = r_node(:,:,:,ele_num)
type_ele(index(i)) = type_ele(ele_num)
size_ele(index(i)) = size_ele(ele_num)
lat_ele(index(i)) = lat_ele(ele_num)
sbox_ele(index(i)) = sbox_ele(ele_num)
end if
ele_num = ele_num - 1
end do
end subroutine delete_elements
end module elements

@ -9,7 +9,7 @@ module io
integer :: outfilenum = 0, infilenum = 0
character(len=100) :: outfiles(10), infiles(10)
logical :: force_overwrite
public
contains
@ -38,7 +38,7 @@ module io
!Check to see if file exists, if it does then ask user if they would like to overwrite the file
inquire(file=trim(temp_outfile), exist=file_exists)
if (file_exists) then
if (file_exists.and.(.not.(force_overwrite))) then
if (overwrite == 'r') print *, "File ", trim(temp_outfile), " already exists. Would you like to overwrite? (Y/N)"
read(*,*) overwrite
if((scan(overwrite, "n") > 0).or.(scan(overwrite, "N") > 0)) then
@ -194,6 +194,7 @@ module io
select case(trim(adjustl(type_ele(i))))
case('fcc')
do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
call apply_periodic(r_interp(:,iatom))
write(11, '(2i16, 3f23.15)') atom_num+iatom, type_interp(iatom), r_interp(:,iatom)
end do
end select
@ -362,7 +363,7 @@ 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
integer :: interp_max, i, j, lat_size, inod, ibasis, ip, unique_index(10), unique_num
real(kind=dp) :: box_vec(3)
1 format('time' / i16, f23.15)
@ -398,6 +399,8 @@ module io
select case(max_ng_node)
case(8)
interp_max = (max_esize)**3
case default
interp_max = 0
end select
write(11,20) interp_max
write(11,3) node_num
@ -441,14 +444,21 @@ module io
!write the element information
if(ele_num > 0) then
write(11,12)
do i = 1, lattice_types
do j = 1, ele_num
if (lat_ele(j) == i) then
lat_size = size_ele(j)-1
exit
!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
write(11,'(3i16)') i, lat_size, basis_type(1,i)
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)
end do
ip = 0
write(11,13)
@ -599,8 +609,8 @@ module io
real(kind=dp), dimension(3), intent(in) :: displace
real(kind = dp), dimension(6), intent(out) :: temp_box_bd
integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, &
new_type_to_type(10), new_lattice_types, sbox
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)
character(len=100) :: etype
real(kind=dp) :: r(3), newdisplace(3)
real(kind=dp), allocatable :: r_innode(:,:,:)
@ -617,6 +627,7 @@ module io
temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i)
end do
call grow_box(temp_box_bd)
!Read in the number of sub_boxes and allocate the variables
read(11, *) n
@ -632,7 +643,11 @@ module io
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)
sub_box_bd(:,sub_box_num+i) = sub_box_bd(:, sub_box_num+i) + displace(:)
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
!Read in sub_box_array_bd
read(11,*) ((sub_box_array_bd(j, k, sub_box_num+i), j = 1, 2), k = 1, 2)
end do
@ -661,14 +676,41 @@ module io
do i = 1, basisnum(j)
basis_type(i,j) = new_type_to_type(basis_type(i,j))
end do
end do
end do
!Read the lattice parameters for every lattice type
read(11,*) (lapa(i), i = lattice_types+1, lattice_types+new_lattice_types)
!Now we loop over all new lattice types and check to see if they are exactly the same as any old lattice types
!If they are then we map the new type to the old type.
k = lattice_types + 1
new_lattice_map(:) = 0
new_loop:do i = 1, new_lattice_types
old_loop:do j = 1, lattice_types
!First check all the lattice level variables
if ((basisnum(i) == basisnum(j)).and. &
(ng_node(i) == ng_node(j))) then
!Now check the basis level variables
do ibasis =1, basisnum(i)
if(basis_type(ibasis,lattice_types+i) /= basis_type(ibasis,j)) then
cycle old_loop
end if
end do
new_lattice_map(i) = j
cycle new_loop
end if
end do old_loop
new_lattice_map(i) = k
k = k+1
end do new_loop
!Read number of elements and atoms and allocate arrays
read(11, *) in_atoms, in_eles
call grow_ele_arrays(in_eles, in_atoms)
allocate(r_innode(3,max_basisnum, max_ng_node))
print *, "Read in ", in_eles, " elements and ", in_atoms, " atoms from ", trim(adjustl(file))
print *, "New box dimensions are: ", box_bd
!Read the atoms
do i = 1, in_atoms
read(11,*) j, type, r(:)
@ -677,15 +719,14 @@ module io
!Read the elements
do i = 1, in_eles
read(11, *) n, type, size, sbox, etype
read(11, *) l, type, size, sbox, etype
do inod = 1, ng_node(type)
do ibasis =1, basisnum(type)
read(11,*) j, k, r_innode(:, ibasis, inod)
r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace
end do
end do
type = type + lattice_types
call add_element(etype, size, type, sbox+n, r_innode)
call add_element(etype, size, new_lattice_map(type), sbox+n, r_innode)
end do
!Close the file being read
@ -693,7 +734,7 @@ 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 = lattice_types + new_lattice_types
if(in_eles > 0) lattice_types = maxval(new_lattice_map)
sub_box_num = sub_box_num + n

@ -1,33 +1,59 @@
program main
!**************************** CACmb *******************************
!* CAC model building toolkit *
! ____________ *
! / / *
! / / *
! /___________/ *
! _|_ _|_ _|____________ *
! / / *
! / / *
! /___________/ *
! *
!*******************************************************************
!**************************** CACmb ********************
!* CAC model building toolkit *
!* ____________ *
!* / / *
!* / / *
!* /___________/ *
!* _|_ _|_ _|____________ *
!* / / *
!* / / *
!* /___________/ *
!* *
!********************************************************
use parameters
use elements
use io
integer :: i, end_mode_arg, arg_num, arg_pos
character(len=100) :: argument
!Print introduction text
print *, '*********************** CACmb *********************'
print *, '* CAC model building toolkit *'
print *, '* _______ *'
print *, '* / / *'
print *, '* / / *'
print *, '* /______ / *'
print *, '* _|_ _|_ _|_______ *'
print *, '* / / *'
print *, '* / / *'
print *, '* /______ / *'
print *, '* *'
print *, '****************************************************'
!Call initialization functions
call lattice_init
call box_init
force_overwrite=.false.
end_mode_arg = 0
! Command line parsing
arg_num = command_argument_count()
!Check to see if overwrite flag is passed
do i = 1, arg_num
call get_command_argument(i,argument)
select case(trim(adjustl(argument)))
case('-ow')
force_overwrite = .true.
print *, "Overwrite flag passed, output files will be overwritten"
end select
end do
!Determine if a mode is being used and what it is. The first argument has to be the mode
!if a mode is being used
call get_command_argument(1, argument)

@ -13,8 +13,6 @@ module mode_convert
integer, intent(out) :: arg_pos
character(len=100) :: infile, outfile
real(kind = dp) :: temp_box_bd(6)
!We have to allocate the element and atom arrays with a size of 1 for the read in code to work
call alloc_ele_arrays(1,1)
!First read in the file
call get_command_argument(2, infile)
call get_in_file(infile)

@ -30,6 +30,8 @@ module mode_create
integer :: i, ibasis, inod
real(kind=dp), allocatable :: r_node_temp(:,:,:)
print *, '-----------------------Mode Create---------------------------'
!Initialize default parameters
orient = reshape((/ 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp /), shape(orient))
cell_mat(:,:)=0.0_dp
@ -91,6 +93,8 @@ module mode_create
end do
else
print *, "Creating 1 element"
call cell_init(lattice_parameter, esize, element_type, orient, cell_mat)
!If the user doesn't pass any build instructions than we just put the cell mat into the element_array
call alloc_ele_arrays(1,0)
@ -110,6 +114,7 @@ module mode_create
!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
if(dup_flag.or.dim_flag) then
!Call the build function with the correct transformation matrix
select case(trim(adjustl(element_type)))
case('fcc')
@ -123,6 +128,8 @@ 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."
!Allocate variables
call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1))
if(lat_atom_num > 0) then

@ -8,7 +8,9 @@ module mode_merge
use elements
character(len=4) :: dim
integer :: in_num
integer :: in_num, new_starts(2)
real(kind=dp) :: shift_vec(3)
logical :: wrap, shift_flag
public
contains
@ -18,12 +20,26 @@ module mode_merge
integer :: i
real(kind=dp) :: displace(3), temp_box_bd(6)
print *, '-----------------------Mode Merge---------------------------'
wrap = .false.
shift_flag = .false.
shift_vec(:) = 0.0_dp
!First we parse the merge command
call parse_command(arg_pos)
!Now loop over all files and stack them
do i = 1, in_num
displace(:) = 0.0_dp
!The new starts variable dictate where in the atom and element array each new
!file starts. This is used for additional options that can be applied to solely
!these new atoms/elements that are read in.
new_starts(1) = atom_num + 1
new_starts(2) = ele_num + 1
if ((i==1).or.(trim(adjustl(dim)) == 'none')) then
call read_in(i, displace, temp_box_bd)
call grow_box(temp_box_bd)
@ -40,6 +56,8 @@ module mode_merge
call read_in(i, displace, temp_box_bd)
call grow_box(temp_box_bd)
end if
if(shift_flag) call shift(new_starts, i)
end do
return
@ -89,7 +107,16 @@ module mode_merge
!Choose what to based on what the option string is
select case(trim(textholder))
case('shift')
shift_flag = .true.
do i = 1,3
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
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
@ -97,4 +124,54 @@ module mode_merge
end do
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
integer, dimension(2), intent(in) :: array_start
integer, intent(in) :: filenum
integer :: i, j, ibasis, inod
real(kind=dp), dimension(3) :: current_shift
!Calculate the current shift which is the filenum-1 multiplied by the user specified shift
current_shift = (filenum-1)*shift_vec
print *, "Atoms/elements from file ", trim(adjustl(infiles(filenum))), " are shifted by ", current_shift
!First shift all the atoms
do i = array_start(1), atom_num
r_atom(:,i) = r_atom(:,i) + current_shift
end do
!Now shift all the elements
do i = array_start(2), ele_num
do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
r_node(:,ibasis, inod, i) = r_node(:,ibasis, inod, i) + current_shift
end do
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
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
box_bd(2*i) = box_bd(2*i) + current_shift(i)
end if
end do
end if
end subroutine shift
end module mode_merge

@ -23,12 +23,12 @@ module opt_disl
contains
subroutine dislocation(option, arg_pos)
!Main calling function for all codes related to dislocations
character(len=100), intent(in) :: option
integer, intent(inout) :: arg_pos
print *, '--------------------Option Dislocation-----------------------'
select case(trim(adjustl(option)))
case('-dislgen')
@ -70,7 +70,11 @@ module opt_disl
call parse_pos(i, textholder, centroid(i))
end do
print *, centroid
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing parameter in dislgen command"
read(textholder, *) b
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
@ -84,26 +88,6 @@ module opt_disl
read(textholder, *) poisson
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, lattice, arglen)
if (arglen==0) STOP "Missing lattice in dislgen command"
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing lattice parameter in dislgen command"
read(textholder, *) lattice_parameter
arg_pos = arg_pos + 1
!Now set the vurgers vector based on the lattice paarameter and lattice type
select case(lattice)
case('fcc')
b = lattice_parameter / sqrt(2.0_dp)
! case('bcc')
! b = lattice_parameter * sqrt(3.0_dp) / 2.0_dp
case default
print *, 'Error: Lattice structure', lattice, ' is not accepted for dislgen option'
STOP
end select
end subroutine parse_dislgen
subroutine dislgen
@ -113,6 +97,8 @@ module opt_disl
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)
print *, "Dislocation with centroid ", centroid, " is inserted"
!Calculate screw and edge burgers vectors
be = sin(char_angle*pi/180.0_dp)*b
bs = cos(char_angle*pi/180.0_dp)*b
@ -168,7 +154,7 @@ module opt_disl
do i = 1, ele_num
do inod=1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
r = r_node(:,ibasis,inod,i)
r = r_node(:,ibasis,inod,i)-centroid
r = matmul(inv_transform, r)
if (r(1) == 0) then
actan = pi/2
@ -270,6 +256,9 @@ module opt_disl
real(kind = dp) :: perimeter, angle, theta, omega, xA(3), xB(3), xC(3), u(3)
real(kind=dp), dimension(:,:), allocatable :: xloop !coordinate of points forming loop
print *, "Dislocation loop with centroid ", centroid, " is inserted"
if(allocated(xLoop)) deallocate(xLoop)
!Define new directions
@ -291,36 +280,52 @@ module opt_disl
a3 = 3
end select
!Calculate loop perimeter
perimeter = 2.0_dp*pi*loop_radius
!Define the number of points forming the loop
! The following criteria are used as a trade-off between
! good accuracy and computational efficiency:
! - each dislocation segment should have a length of 5 angströms;
! - the loop should contain at least 3 points
! (for very small loops, this will result in segments shorter than 5 A);
! - there should not be more than 100 points
! (for very large loops, this will result in segments longer than 5 A).
Npoints = MAX( 3 , MIN( NINT(perimeter/5.d0) , 100 ) )
!angle between two consecutive points
theta = 2.0_dp*pi / dble(Npoints)
!allocate xLoop
allocate(xLoop(Npoints,3))
xLoop(:,:) = 0.0_dp
!Calculate the position of each point in the loop
angle = 0.0_dp
do i = 1, size(xLoop,1)
xLoop(i,a1) = centroid(a1) + loop_radius*dcos(angle)
xLoop(i,a2) = centroid(a2) + loop_radius*dsin(angle)
xLoop(i,a3) = centroid(a3)
! Increment angle for next point
angle = angle + theta
end do
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)
else
!Calculate loop perimeter
perimeter = 2.0_dp*pi*loop_radius
!Define the number of points forming the loop
! The following criteria are used as a trade-off between
! good accuracy and computational efficiency:
! - each dislocation segment should have a length of 5 angströms;
! - the loop should contain at least 3 points
! (for very small loops, this will result in segments shorter than 5 A);
! - there should not be more than 100 points
! (for very large loops, this will result in segments longer than 5 A).
Npoints = MAX( 3 , MIN( NINT(perimeter/5.d0) , 100 ) )
!angle between two consecutive points
theta = 2.0_dp*pi / dble(Npoints)
!allocate xLoop
allocate(xLoop(Npoints,3))
xLoop(:,:) = 0.0_dp
!Calculate the position of each point in the loop
angle = 0.0_dp
do i = 1, size(xLoop,1)
xLoop(i,a1) = centroid(a1) + loop_radius*dcos(angle)
xLoop(i,a2) = centroid(a2) + loop_radius*dsin(angle)
xLoop(i,a3) = centroid(a3)
! Increment angle for next point
angle = angle + theta
end do
end if
!Now actually calculate the displacement created by a loop for every atom
do i = 1, atom_num
u = 0.0_dp

@ -0,0 +1,259 @@
module opt_group
!This module contains all code associated with dislocations
use parameters
use elements
use subroutines
use box
implicit none
integer :: group_ele_num, group_atom_num, remesh_size
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
integer, allocatable :: element_index(:), atom_index(:)
public
contains
subroutine group(arg_pos)
!Main calling function for the group option
integer, intent(inout) :: arg_pos
print *, '-----------------------Option Group-------------------------'
group_ele_num = 0
group_atom_num = 0
remesh_size=0
if(allocated(element_index)) deallocate(element_index)
if(allocated(atom_index)) deallocate(atom_index)
call parse_group(arg_pos)
call get_group
!Now call the transformation functions for the group
if(displace) call displace_group
if(remesh_size > 0) call remesh_group
end subroutine group
subroutine parse_group(arg_pos)
!Parse the group command
integer, intent(inout) :: arg_pos
integer :: i,arglen
character(len=100) :: textholder
!Parse type and shape command
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, type, arglen)
if (arglen==0) STOP "Missing select_type in group command"
select case(trim(adjustl(type)))
case('atoms', 'elements', 'both')
continue
case default
print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", &
"Please select from atoms, nodes, or both."
end select
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, shape, arglen)
if (arglen==0) STOP "Missing group_shape in group command"
!Now parse the arguments required by the user selected shape
select case(trim(adjustl(shape)))
case('block')
do i= 1, 6
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing block boundary in dislgen command"
call parse_pos(int((i+1)/2), textholder, block_bd(i))
end do
case default
print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", &
"for accepted group shapes."
end select
!Now parse the additional options which may be present
do while(.true.)
if(arg_pos > command_argument_count()) exit
!Pull out the next argument which should either be a keyword or an option
arg_pos=arg_pos+1
call get_command_argument(arg_pos, textholder)
textholder=adjustl(textholder)
!Choose what to based on what the option string is
select case(trim(textholder))
case('displace')
displace = .true.
do i = 1,3
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) stop "Missing vector component for shift command"
read(textholder, *) disp_vec(i)
end do
case('wrap')
wrap = .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
case default
!If it isn't an available option to opt_disl then we just exit
exit
end select
end do
end subroutine parse_group
subroutine get_group
!This subroutine finds all elements and/or atoms within the group boundaries
!specified by the user.
integer :: i, inod, ibasis
integer, allocatable :: resize_array(:)
real(kind=dp) :: r_center(3)
select case(trim(adjustl(shape)))
case('block')
print *, "Group has block shape with boundaries: ", block_bd
!Allocate variables to arbitrary size
allocate(element_index(1024), atom_index(1024))
!Check the type to see whether we need to find the elements within the group
select case(trim(adjustl(type)))
case('elements', 'both')
do i = 1, ele_num
r_center(:) = 0.0_dp
do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
r_center = r_center + r_node(:,ibasis,inod,i)/(basisnum(lat_ele(i))*ng_node(lat_ele(i)))
end do
end do
if (in_block_bd(r_center, block_bd)) then
group_ele_num = group_ele_num + 1
if(group_ele_num > size(element_index)) then
allocate(resize_array(size(element_index) + 1024))
resize_array(1:group_ele_num-1) = element_index
resize_array(group_ele_num:) = 0
call move_alloc(resize_array, element_index)
end if
element_index(group_ele_num) = i
end if
end do
end select
!Check the type to see if we need to find the atoms within the group
select case(trim(adjustl(type)))
case('atoms','both')
do i = 1, atom_num
if(in_block_bd(r_atom(:,i),block_bd)) then
group_atom_num = group_atom_num + 1
if (group_atom_num > size(atom_index)) then
allocate(resize_array(size(atom_index) + 1024))
resize_array(1:group_atom_num -1) = atom_index
resize_array(group_atom_num:) = 0
call move_alloc(resize_array, atom_index)
end if
atom_index(group_atom_num) = i
end if
end do
end select
end select
print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms."
end subroutine get_group
subroutine displace_group
!This subroutine applies a displacement to elements/atoms in the groups
integer :: i, inod, ibasis
print *, "Elements/atoms in group displaced by ", disp_vec
!Displace atoms
do i = 1, group_atom_num
r_atom(:,atom_index(i)) = r_atom(:,atom_index(i)) + disp_vec
end do
!Displace elements
do i = 1, group_ele_num
do inod = 1, ng_node(lat_ele(element_index(i)))
do ibasis = 1, basisnum(lat_ele(element_index(i)))
r_node(:,ibasis,inod,element_index(i)) = r_node(:,ibasis,inod,element_index(i)) + disp_vec
end do
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
do i = 1,3
if (disp_vec(i) < -lim_zero) then
box_bd(2*i-1) = box_bd(2*i-1) - disp_vec(i)
else if (disp_vec(i) > lim_zero) then
box_bd(2*i) = box_bd(2*i) + disp_vec(i)
end if
end do
end if
end subroutine displace_group
subroutine remesh_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)
!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 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
!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
end subroutine remesh_group
end module opt_group

@ -178,22 +178,36 @@ module subroutines
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
read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
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
read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
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
read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos
pos = (box_bd(2*i)-box_bd(2*i-1))*pos
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
@ -249,4 +263,24 @@ module subroutines
end do
end subroutine siftdown
subroutine apply_periodic(r)
!This function checks if an atom is outside the box and wraps it back in. This is generally used when some
!kind of displacement is applied but the simulation cell is desired to be maintained as the same size.
real(kind=dp), dimension(3), intent(inout) :: r
integer :: j
real(kind=dp) ::box_len
do j = 1, 3
box_len = box_bd(2*j) - box_bd(2*j-1)
if (r(j) > box_bd(2*j)) then
r(j) = r(j) - box_len
else if (r(j) < box_bd(2*j-1)) then
r(j) = r(j) + box_len
end if
end do
end subroutine
end module subroutines
Loading…
Cancel
Save