Compare commits

...

10 Commits

@ -33,7 +33,7 @@ $(OBJDIR)/%.o: %.f90
deps: deps:
@fortdepend -o Makefile.dep -i mpi -b obj -w @makedepf90 -b obj *.f90 > Makefile.dep
#----------------- CLEAN UP -----------------# #----------------- CLEAN UP -----------------#

@ -1,26 +1,190 @@
obj/atoms.o : atoms.f90 obj/functions.o obj/parameters.o # This file is generated automatically. DO NOT EDIT!
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/main : \
obj/elements.o : elements.f90 obj/box.o obj/sorts.o obj/subroutines.o obj/functions.o obj/parameters.o obj/atoms.o \
obj/functions.o : functions.f90 obj/parameters.o obj/box.o \
obj/io.o : io.f90 obj/str.o obj/box.o obj/atoms.o obj/parameters.o obj/elements.o obj/caller.o \
obj/main.o : main.f90 obj/caller.o obj/io.o obj/elements.o obj/parameters.o obj/elements.o \
obj/mode_calc.o : mode_calc.f90 obj/box.o obj/elements.o obj/subroutines.o obj/io.o obj/parameters.o obj/functions.o \
obj/mode_convert.o : mode_convert.f90 obj/io.o obj/elements.o obj/box.o obj/parameters.o obj/io.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/main.o \
obj/mode_da.o : mode_da.f90 obj/neighbors.o obj/elements.o obj/io.o obj/parameters.o obj/mode_calc.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_convert.o \
obj/mode_metric.o : mode_metric.f90 obj/neighbors.o obj/elements.o obj/io.o obj/parameters.o obj/mode_create.o \
obj/neighbors.o : neighbors.f90 obj/functions.o obj/subroutines.o obj/elements.o obj/parameters.o obj/mode_da.o \
obj/opt_bubble.o : opt_bubble.f90 obj/opt_group.o obj/box.o obj/elements.o obj/parameters.o obj/mode_merge.o \
obj/opt_deform.o : opt_deform.f90 obj/box.o obj/elements.o obj/subroutines.o obj/parameters.o obj/mode_metric.o \
obj/opt_delete.o : opt_delete.f90 obj/neighbors.o obj/elements.o obj/subroutines.o obj/parameters.o obj/neighbors.o \
obj/opt_disl.o : opt_disl.f90 obj/box.o obj/subroutines.o obj/elements.o obj/parameters.o obj/opt_bubble.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_deform.o \
obj/opt_orient.o : opt_orient.f90 obj/box.o obj/elements.o obj/subroutines.o obj/parameters.o obj/opt_delete.o \
obj/opt_redef_box.o : opt_redef_box.f90 obj/subroutines.o obj/elements.o obj/box.o obj/opt_disl.o \
obj/opt_slip_plane.o : opt_slip_plane.f90 obj/subroutines.o obj/functions.o obj/elements.o obj/parameters.o obj/opt_group.o \
obj/parameters.o : parameters.f90 obj/opt_orient.o \
obj/sorts.o : sorts.f90 obj/parameters.o obj/opt_redef_box.o \
obj/str.o : str.f90 obj/opt_slip_plane.o \
obj/subroutines.o : subroutines.f90 obj/str.o obj/box.o obj/functions.o obj/parameters.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/atoms.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/box.o \
obj/elements.o \
obj/functions.o \
obj/parameters.o \
obj/subroutines.o
obj/opt_bubble.o : \
obj/atoms.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/atoms.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/str.o

@ -60,6 +60,8 @@ module caller
call group(arg_pos) call group(arg_pos)
case('-ow') case('-ow')
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
case('-novel')
arg_pos = arg_pos+1
case('-wrap') case('-wrap')
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
case('-norefine') case('-norefine')

@ -2,6 +2,7 @@ module elements
!This module contains the elements data structures, structures needed for building regions !This module contains the elements data structures, structures needed for building regions
!and operations that are done on elements !and operations that are done on elements
use atoms
use parameters use parameters
use functions use functions
use subroutines use subroutines
@ -31,20 +32,21 @@ module elements
!Mapping atom type to provided name !Mapping atom type to provided name
character(len=2), dimension(10) :: type_to_name character(len=2), dimension(10) :: type_to_name
integer :: atom_types = 0 integer :: atom_types = 0
real(kind=dp) :: masses(10)
!Variables for creating elements based on primitive cells !Variables for creating elements based on primitive cells
real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3), & real(kind=dp) :: cubic_cell(3,8), fcc_cell(3,8), fcc_mat(3,3), fcc_inv(3,3), bcc_cell(3,8), bcc_mat(3,3), bcc_inv(3,3), &
cube20(3,20) cube20(3,20)
integer :: cubic_faces(4,6) integer :: cubic_faces(4,6), oneonetwopairs(2,6)
!Below are lattice type arrays which provide information on the general form of the elements. !Below are lattice type arrays which provide information on the general form of the elements.
!We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased. !We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased.
integer :: lattice_types = 0 integer :: lattice_types = 0
integer :: max_ng_node, ng_node(10) !Max number of nodes per element and number of nodes per element for each lattice type integer :: max_ng_node, ng_node(10) !Max number of nodes per element and number of nodes per element for each lattice type
integer :: max_esize=0 !Maximum number of atoms per side of element integer :: max_esize=0 !Maximum number of atoms per side of element
real(kind=dp) :: lapa(10), masses(10) real(kind=dp) :: lapa(10)
!These variables contain information on the basis, for simplicities sake we limit !These variables contain informatione on the basis, for simplicities sake we limit
!the user to the definition of 10 lattice types with 10 basis atoms at each lattice point. !the user to the definition of 10 lattice types with 10 basis atoms at each lattice point.
!This can be easily increased with no change to efficiency !This can be easily increased with no change to efficiency
integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type
@ -104,6 +106,14 @@ module elements
cubic_faces(:,5) = (/ 4, 3, 7, 8 /) cubic_faces(:,5) = (/ 4, 3, 7, 8 /)
cubic_faces(:,6) = (/ 5, 6, 7, 8 /) cubic_faces(:,6) = (/ 5, 6, 7, 8 /)
!Now mark which node pairs represent the 112 directions. This is used for the dislocation analysis.
oneonetwopairs(:,1) = (/ 1, 3 /)
oneonetwopairs(:,2) = (/ 1, 6 /)
oneonetwopairs(:,3) = (/ 2, 7 /)
oneonetwopairs(:,4) = (/ 1, 8 /)
oneonetwopairs(:,5) = (/ 4, 7 /)
oneonetwopairs(:,6) = (/ 5, 7 /)
!!Now initialize the fcc primitive cell !!Now initialize the fcc primitive cell
fcc_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, & fcc_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
0.5_dp, 0.5_dp, 0.0_dp, & 0.5_dp, 0.5_dp, 0.0_dp, &
@ -345,9 +355,9 @@ module elements
end subroutine add_atom end subroutine add_atom
subroutine add_atom_type(type, inttype, force_new) subroutine add_atom_type(mass, inttype, force_new)
!This subroutine adds a new atom type to the module list of atom types !This subroutine adds a new atom type to the module list of atom types
character(len=2), intent(in) :: type real(kind=dp), intent(in) :: mass
integer, intent(out) :: inttype integer, intent(out) :: inttype
logical, optional, intent(in) :: force_new logical, optional, intent(in) :: force_new
@ -364,7 +374,7 @@ module elements
exists = .false. exists = .false.
if(.not.force) then if(.not.force) then
do i=1, 10 do i=1, 10
if(type == type_to_name(i)) then if(is_equal(mass, masses(i))) then
exists = .true. exists = .true.
inttype = i inttype = i
exit exit
@ -378,7 +388,8 @@ module elements
print *, "Defined atom types are greater than 10 which is currently not supported." print *, "Defined atom types are greater than 10 which is currently not supported."
stop 3 stop 3
end if end if
type_to_name(atom_types) = type call atommassspecies(mass, type_to_name(atom_types))
masses(atom_types) = mass
inttype = atom_types inttype = atom_types
end if end if
return return
@ -426,10 +437,10 @@ module elements
real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: r_interp !Interpolated atomic positions real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: r_interp !Interpolated atomic positions
real(kind = dp), optional, intent(in) :: eng(:,:), f(:,:,:), & real(kind = dp), optional, intent(in) :: eng(:,:), f(:,:,:), &
v(:,:,:) v(:,:,:)
real(kind=dp), dimension(10, max_basisnum*max_esize**3), optional,intent(out) :: data_interp !Interpolated atomic positions real(kind=dp), dimension(:,:), optional,intent(out) :: data_interp !Interpolated atomic positions
!Internal variables !Internal variables
integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp, adjust
real(kind=dp) :: a_shape(max_ng_node) real(kind=dp) :: a_shape(max_ng_node)
real(kind=dp) :: t, s, r real(kind=dp) :: t, s, r
@ -451,16 +462,15 @@ module elements
lat_type_temp = lat_type lat_type_temp = lat_type
end select end select
select case(type) select case(type)
case('fcc','bcc') case('fcc','bcc')
!Now loop over all the possible sites !Now loop over all the possible sites
do it = 1, esize do it = 1, esize
t = (1.0_dp*(it-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2) t=-1.0_dp +(it-1)*(2.0_dp/(esize-1))
do is =1, esize do is =1, esize
s = (1.0_dp*(is-1)-(esize-1)/2)/(1.0_dp*(esize-1)/2) s=-1.0_dp +(is-1)*(2.0_dp/(esize-1))
do ir = 1, esize do ir = 1, esize
r = (1.0_dp*(ir-1) - (esize-1)/2)/(1.0_dp*(esize-1)/2) r=-1.0_dp +(ir-1)*(2.0_dp/(esize-1))
call rhombshape(r,s,t,a_shape) call rhombshape(r,s,t,a_shape)
do ibasis = 1, bnum do ibasis = 1, bnum
@ -518,6 +528,68 @@ module elements
return return
end subroutine interpolate_atoms end subroutine interpolate_atoms
subroutine interpolate_vel(type, esize, lat_type, vel_in, vel_interp)
!This subroutine returns the interpolated atoms from the elements.
!Arguments
character(len=*), intent(in) :: type !The type of element that it is
integer, intent(in) :: esize !The number of atoms per side
integer, intent(in) :: lat_type !The integer lattice type of the element
real(kind=dp), dimension(:,:,:), intent(in) :: vel_in !Nodal positions
real(kind=dp), dimension(3, max_basisnum*max_esize**3), intent(out) :: vel_interp !Interpolated atomic positions
!Internal variables
integer :: it, is, ir, ibasis, inod, ia, bnum, lat_type_temp, adjust
real(kind=dp) :: a_shape(max_ng_node)
real(kind=dp) :: t, s, r
!Initialize some variables
vel_interp(:,:) = 0.0_dp
ia = 0
!Define bnum based on the input lattice type. If lat_type=0 then we are interpolating lattice points which means
!the basis is 0,0,0, and the type doesn't matter
select case(lat_type)
case(0)
bnum=1
lat_type_temp = 1
case default
bnum = basisnum(lat_type)
lat_type_temp = lat_type
end select
select case(type)
case('fcc','bcc')
!Now loop over all the possible sites
do it = 1, esize
t=-1.0_dp +(it-1)*(2.0_dp/(esize-1))
do is =1, esize
s=-1.0_dp +(is-1)*(2.0_dp/(esize-1))
do ir = 1, esize
r=-1.0_dp +(ir-1)*(2.0_dp/(esize-1))
call rhombshape(r,s,t,a_shape)
do ibasis = 1, bnum
ia = ia + 1
do inod = 1, 8
vel_interp(:,ia) = vel_interp(:,ia) + a_shape(inod) * vel_in(:,ibasis,inod)
end do
end do
end do
end do
end do
end select
if (ia /= bnum*esize**3) then
print *, "Incorrect interpolation"
stop 3
end if
return
end subroutine interpolate_vel
subroutine rhombshape(r,s,t, shape_fun) subroutine rhombshape(r,s,t, shape_fun)
!Shape function for rhombohedral elements !Shape function for rhombohedral elements
real(kind=8), intent(in) :: r, s, t real(kind=8), intent(in) :: r, s, t
@ -651,6 +723,7 @@ module elements
end do end do
end subroutine wrap_atoms end subroutine wrap_atoms
subroutine wrap_elements subroutine wrap_elements
integer :: i, inod, ibasis, j, node_in_bd(3,max_ng_node) integer :: i, inod, ibasis, j, node_in_bd(3,max_ng_node)
real(kind =dp) :: box_len(3) real(kind =dp) :: box_len(3)
@ -1008,11 +1081,36 @@ do i = 1, atom_num
subroutine add_atom_data(ia, eng, force, virial) subroutine add_atom_data(ia, eng, force, virial)
!Function which sets the atom data arrays !Function which sets the atom data arrays
integer, intent(in) :: ia integer, intent(in) :: ia
real(kind=dp), intent(in) :: eng, force(3), virial(6) real(kind=dp), intent(in) :: eng, force(:), virial(:)
integer :: size_atom_array
real(kind=dp), allocatable :: tmp_eng(:), tmp_force(:,:), tmp_virial(:,:)
size_atom_array=size(tag_atom)
if(ia > atom_num) then
stop "ia in add_atom_dat shouldn't be larger than atom_num"
else if(ia > size(energy_atom)) then
allocate(tmp_eng(size(energy_atom)+1024))
allocate(tmp_force(3, size(energy_atom)+1024))
allocate(tmp_virial(6, size(energy_atom)+1024))
tmp_eng=0.0_dp
tmp_eng(1:size(energy_atom)) = energy_atom
call move_alloc(tmp_eng, energy_atom)
tmp_force=0.0_dp
tmp_force(:, 1:size(force_atom, 2)) = force_atom
call move_alloc(tmp_force, force_atom)
tmp_virial=0.0_dp
tmp_virial(:, 1:size(virial_atom,2)) = virial_atom
call move_alloc(tmp_virial, virial_atom)
end if
energy_atom(ia) = eng energy_atom(ia) = eng
force_atom(:,ia) = force(:) force_atom(:,ia) = force(:)
virial_atom(:,ia) = virial(:) virial_atom(:,ia) = virial(:)
vflag = .true. vflag = .true.
return return
end subroutine add_atom_data end subroutine add_atom_data

@ -168,11 +168,11 @@ END FUNCTION StrDnCase
do i =1 ,3 do i =1 ,3
!Check upper bound !Check upper bound
if(v(i) > (box_bd(2*i)+10.0_dp**(-10)) ) then if(v(i) > (box_bd(2*i))) then
in_block_bd =.false. in_block_bd =.false.
exit exit
!Check lower bound !Check lower bound
else if (v(i) < (box_bd(2*i-1)-10.0_dp**(-10))) then else if (v(i) < (box_bd(2*i-1))+1d-6) then
in_block_bd = .false. in_block_bd = .false.
exit exit
end if end if

@ -139,7 +139,7 @@ module io
do i = 1, ele_num do i = 1, ele_num
do inod = 1, ng_node(lat_ele(i)) do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i)) do ibasis = 1, basisnum(lat_ele(i))
write(11, '(2i16, 3f23.15)') basis_type(ibasis,lat_ele(i)), 1, r_node(:,ibasis,inod,i) write(11, *) basis_type(ibasis,lat_ele(i)), 1, r_node(:,ibasis,inod,i)
outn = outn + 1 outn = outn + 1
end do end do
end do end do
@ -151,7 +151,7 @@ module io
!Write atom positions !Write atom positions
do i = 1, atom_num do i = 1, atom_num
write(11, '(2i16, 3f23.15)') type_atom(i), 0, r_atom(:,i) write(11, *) type_atom(i), 0, r_atom(:,i)
outn = outn + 1 outn = outn + 1
end do end do
@ -207,7 +207,7 @@ module io
write(11, '(a)') 'Atoms' write(11, '(a)') 'Atoms'
write(11, '(a)') ' ' write(11, '(a)') ' '
do i = 1, atom_num do i = 1, atom_num
write(11, '(2i16, 3f23.15)') tag_atom(i), type_atom(i), r_atom(:,i) write(11, *) tag_atom(i), type_atom(i), r_atom(:,i)
end do end do
!Write refined element atomic positions !Write refined element atomic positions
@ -219,7 +219,7 @@ module io
do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
interp_num = interp_num+1 interp_num = interp_num+1
call apply_periodic(r_interp(:,iatom)) call apply_periodic(r_interp(:,iatom))
write(11, '(2i16, 3f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom) write(11, *) atom_num+interp_num, type_interp(iatom), r_interp(:,iatom)
end do end do
end select end select
end do end do
@ -277,12 +277,12 @@ module io
if (write_dat) then if (write_dat) then
do i = 1, atom_num do i = 1, atom_num
write(11, '(2i16, 13f23.15)') tag_atom(i), type_atom(i), r_atom(:,i), energy_atom(i), & write(11, *) tag_atom(i), type_atom(i), r_atom(:,i), energy_atom(i), &
force_atom(:,i), virial_atom(:,i) force_atom(:,i), virial_atom(:,i)
end do end do
else else
do i = 1, atom_num do i = 1, atom_num
write(11, '(2i16, 3f23.15)') tag_atom(i), type_atom(i), r_atom(:,i) write(11, *) tag_atom(i), type_atom(i), r_atom(:,i)
end do end do
end if end if
@ -294,12 +294,12 @@ module io
do inod =1, ng_node(lat_ele(i)) do inod =1, ng_node(lat_ele(i))
do ibasis=1, basisnum(lat_ele(i)) do ibasis=1, basisnum(lat_ele(i))
if(write_dat) then if(write_dat) then
write(11, '(2i16, 13f23.15)') atom_num+interp_num, basis_type(ibasis,lat_ele(i)), & write(11, *) atom_num+interp_num, basis_type(ibasis,lat_ele(i)), &
r_node(:, ibasis, inod, i), energy_node(ibasis,inod,i), force_node(:, ibasis, inod, i), & r_node(:, ibasis, inod, i), energy_node(ibasis,inod,i), force_node(:, ibasis, inod, i), &
virial_node(:, ibasis, inod, i) virial_node(:, ibasis, inod, i)
else else
write(11, '(2i16, 3f23.15)') atom_num+interp_num, basis_type(ibasis,lat_ele(i)), & write(11, *) atom_num+interp_num, basis_type(ibasis,lat_ele(i)), &
r_node(:, ibasis, inod, i) r_node(:, ibasis, inod, i)
end if end if
interp_num = interp_num +1 interp_num = interp_num +1
@ -321,14 +321,14 @@ module io
do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
interp_num = interp_num+1 interp_num = interp_num+1
call apply_periodic(r_interp(:,iatom)) call apply_periodic(r_interp(:,iatom))
write(11, '(2i16, 13f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom), & write(11, *) atom_num+interp_num, type_interp(iatom), r_interp(:,iatom), &
data_interp(:,iatom) data_interp(:,iatom)
end do end do
else else
do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3 do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
interp_num = interp_num+1 interp_num = interp_num+1
call apply_periodic(r_interp(:,iatom)) call apply_periodic(r_interp(:,iatom))
write(11, '(2i16, 3f23.15)') atom_num+interp_num, type_interp(iatom), r_interp(:,iatom) write(11, *) atom_num+interp_num, type_interp(iatom), r_interp(:,iatom)
end do end do
end if end if
end select end select
@ -436,6 +436,7 @@ module io
16 format(/'POINT_DATA', i16) 16 format(/'POINT_DATA', i16)
17 format('SCALARS weight float', / & 17 format('SCALARS weight float', / &
'LOOKUP_TABLE default') 'LOOKUP_TABLE default')
18 format('SCALARS atom_type float', / & 18 format('SCALARS atom_type float', / &
'LOOKUP_TABLE default') 'LOOKUP_TABLE default')
@ -445,6 +446,9 @@ module io
21 format('SCALARS esize float', /& 21 format('SCALARS esize float', /&
'LOOKUP_TABLE default') 'LOOKUP_TABLE default')
22 format('SCALARS energy float', /&
'LOOKUP_TABLE default')
!First we write the vtk file containing the atoms !First we write the vtk file containing the atoms
open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind') open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
@ -467,6 +471,14 @@ module io
do i = 1, atom_num do i = 1, atom_num
write(11, '(i16)') type_atom(i) write(11, '(i16)') type_atom(i)
end do end do
!Write the atom data
if(allocated(force_atom)) then
write(11,22)
do i =1, atom_num
write(11, "(f23.15)") energy_atom(i)
end do
end if
close(11) close(11)
!Now we write the vtk file for the elements !Now we write the vtk file for the elements
@ -499,6 +511,18 @@ module io
do i = 1, ele_num do i = 1, ele_num
write(11, '(i16)') size_ele(i) write(11, '(i16)') size_ele(i)
end do end do
if(allocated(force_node)) then
write(11,16) node_num
write(11,22)
do i = 1, ele_num
do inod=1, ng_node(lat_ele(i))
do ibasis=1, basisnum(lat_ele(i))
write(11,'(f23.15)') energy_node(ibasis, inod, i)
end do
end do
end do
end if
close(11) close(11)
end subroutine end subroutine
@ -509,7 +533,7 @@ module io
character(len=100), intent(in) :: file character(len=100), intent(in) :: file
integer :: interp_max, i, j, inod, ibasis, ip, unique_index(50), unique_size(50), unique_type(50), unique_num, & integer :: interp_max, i, j, inod, ibasis, ip, unique_index(50), unique_size(50), unique_type(50), unique_num, &
etype etype
real(kind=dp) :: box_vec(3), masses(10) real(kind=dp) :: box_vec(3)
1 format('time' / i16, f23.15) 1 format('time' / i16, f23.15)
2 format('number of elements' / i16) 2 format('number of elements' / i16)
@ -548,9 +572,6 @@ module io
write(11,3) node_num write(11,3) node_num
write(11,19) max_ng_node, max_basisnum write(11,19) max_ng_node, max_basisnum
write(11,5) atom_num write(11,5) atom_num
do i = 1, atom_types
call atommass(type_to_name(i),masses(i))
end do
write(11,*) "masses: " write(11,*) "masses: "
write(11, *) (masses(i), i = 1, atom_types) write(11, *) (masses(i), i = 1, atom_types)
write(11,7) box_bc(1:1), box_bc(2:2), box_bc(3:3) write(11,7) box_bc(1:1), box_bc(2:2), box_bc(3:3)
@ -603,12 +624,12 @@ module io
if(allocated(vel_atom)) then if(allocated(vel_atom)) then
write(11, 16) write(11, 16)
do i = 1, atom_num do i = 1, atom_num
write(11, '(2i16, 6f23.15)') i, type_atom(i), r_atom(:,i), vel_atom(:, i) write(11, '(2i16, 6f23.15)') tag_atom(i), type_atom(i), r_atom(:,i), vel_atom(:, i)
end do end do
else else
write(11,14) write(11,14)
do i = 1, atom_num do i = 1, atom_num
write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i) write(11, '(2i16, 3f23.15)') tag_atom(i), type_atom(i), r_atom(:,i)
end do end do
end if end if
end if end if
@ -631,7 +652,7 @@ module io
write(11,*) box_bd(:) write(11,*) box_bd(:)
!Write the number of atom types in the current model and all of their names !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) write(11,*) atom_types, (masses(i), i=1, atom_types)
!Write the number of lattice_types, basisnum and number of nodes for each lattice type !Write the number of lattice_types, basisnum and number of nodes for each lattice type
write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types) write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types)
!Now for every lattice type write the basis atom types !Now for every lattice type write the basis atom types
@ -754,9 +775,8 @@ module io
new_type_to_type(10), new_lattice_types, 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) node_map(10)
character(len=100) :: etype character(len=100) :: etype
real(kind=dp) :: r(3), newdisplace(3), lapa_map(10) real(kind=dp) :: r(3), newdisplace(3), lapa_map(10), new_masses(10)
real(kind=dp), allocatable :: r_innode(:,:,:) real(kind=dp), allocatable :: r_innode(:,:,:)
character(len = 2) :: new_type_to_name(10)
!First open the file !First open the file
open(unit=11, file=trim(adjustl(file)), action='read',position='rewind') open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
@ -777,11 +797,11 @@ module io
call grow_box(temp_box_bd) call grow_box(temp_box_bd)
!Read in the number of atom types and all their names !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) read(11, *) new_atom_types, (new_masses(i), i = 1, new_atom_types)
!Now fit these into the global list of atom types, after this new_type_to_type is the actual global !Now fit these into the global list of atom types, after this new_type_to_type is the actual global
!type of the atoms within this file !type of the atoms within this file
do i = 1, new_atom_types do i = 1, new_atom_types
call add_atom_type(new_type_to_name(i), new_type_to_type(i), all_new) call add_atom_type(new_masses(i), new_type_to_type(i), all_new)
end do end do
!Read the number of lattice types, basisnum, and number of nodes for each lattice type !Read the number of lattice types, basisnum, and number of nodes for each lattice type
read(11,*) new_lattice_types, (bnum_map(i), i = 1, new_lattice_types),(node_map(i), i =1, new_lattice_types) read(11,*) new_lattice_types, (bnum_map(i), i = 1, new_lattice_types),(node_map(i), i =1, new_lattice_types)
@ -893,8 +913,7 @@ module io
!Read define atom_types by mass !Read define atom_types by mass
print *, j print *, j
do i = 1, j do i = 1, j
call realatomspecies(atomic_masses(i), atomic_element) call add_atom_type(atomic_masses(i), atom_type_map(i), all_new)
call add_atom_type(atomic_element, atom_type_map(i), all_new)
print *, i, atom_type_map(i) print *, i, atom_type_map(i)
end do end do
@ -979,7 +998,7 @@ module io
!Internal Variables !Internal Variables
integer :: i, j, in_eles, in_atoms, inbtypes(10), lat_type, ia, ie, inod, & integer :: i, j, in_eles, in_atoms, inbtypes(10), lat_type, ia, ie, inod, &
id, type_node, ilat, esize, tag, type, bnum, n, ibasis, ip, atom_type_map(100) id, type_node, ilat, esize, tag, atype, bnum, n, ibasis, ip, atom_type_map(100)
real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), vela(3),& real(kind=dp) :: newdisplace(3), ra(3), in_lapa, ea, fa(3), va(6), vela(3),&
ee(10,20), fe(3,10,20), ve(6,10,20), re(3,10,20), vele(3,10,20), atomic_masses(10) ee(10,20), fe(3,10,20), ve(6,10,20), re(3,10,20), vele(3,10,20), atomic_masses(10)
character(len=100) :: textholder, fcc character(len=100) :: textholder, fcc
@ -1029,8 +1048,7 @@ module io
!Read define atom_types by mass !Read define atom_types by mass
do i = 1, j do i = 1, j
call realatomspecies(atomic_masses(i), atomic_element) call add_atom_type(atomic_masses(i), atom_type_map(i), all_new)
call add_atom_type(atomic_element, atom_type_map(i), all_new)
end do end do
if(in_atoms > 0 ) then if(in_atoms > 0 ) then
@ -1046,11 +1064,11 @@ module io
do ia = 1, in_atoms do ia = 1, in_atoms
read(11,'(a)') line(:) read(11,'(a)') line(:)
if(read_vel) then if(read_vel) then
read(line,*) tag, type, ra(:), ea, fa(:), va(:), vela(:) read(line,*) tag, atype, ra(:), ea, fa(:), va(:), vela(:)
else else
read(line,*) tag, type, ra(:), ea, fa(:), va(:) read(line,*) tag, atype, ra(:), ea, fa(:), va(:)
end if end if
call add_atom(tag, atom_type_map(type), ra) call add_atom(tag, atom_type_map(atype), ra)
call add_atom_data(atom_num, ea, fa, va) call add_atom_data(atom_num, ea, fa, va)
if(read_vel) vel_atom(:,atom_num) = vela if(read_vel) vel_atom(:,atom_num) = vela
end do end do
@ -1152,8 +1170,7 @@ module io
!Read atomic masses !Read atomic masses
do i = 1, type_in do i = 1, type_in
read(11,*) j, mass read(11,*) j, mass
call realatomspecies(mass, atom_species) call add_atom_type(mass, type_map(i), all_new)
call add_atom_type(atom_species, type_map(i), all_new)
end do end do
!Read useless info !Read useless info
@ -1226,8 +1243,7 @@ module io
!Read atomic masses !Read atomic masses
do i = 1, type_in do i = 1, type_in
read(11,*) j, mass, textholder read(11,*) j, mass, textholder
call realatomspecies(mass, atom_species) call add_atom_type(mass, type_map(i), all_new)
call add_atom_type(atom_species, type_map(i), all_new)
end do end do
!Read useless info !Read useless info
@ -1311,6 +1327,7 @@ module io
integer :: i, j,arglen, arg_pos, ntypes integer :: i, j,arglen, arg_pos, ntypes
character(len=100) :: textholder character(len=100) :: textholder
real(kind=dp) :: mass
arg_pos = apos + 1 arg_pos = apos + 1
call get_command_argument(arg_pos, textholder, arglen) call get_command_argument(arg_pos, textholder, arglen)
@ -1320,7 +1337,8 @@ module io
do i=1,ntypes do i=1,ntypes
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen) call get_command_argument(arg_pos, textholder, arglen)
call add_atom_type(textholder, j) call atommass(textholder, mass)
call add_atom_type(mass, j)
end do end do
return return

@ -25,6 +25,9 @@ module mode_calc
case('tot_virial') case('tot_virial')
allocate(calculated(6)) allocate(calculated(6))
call calc_tot_virial call calc_tot_virial
case('tot_energy')
allocate(calculated(1))
call calc_tot_energy
case default case default
print *, trim(adjustl(calc_opt)), " is not accepted as a calc option in mode_calc" print *, trim(adjustl(calc_opt)), " is not accepted as a calc option in mode_calc"
stop 3 stop 3
@ -87,5 +90,26 @@ module mode_calc
print *, "Total virial is calculated as : (v11, v22, v33, v32, v31, v21)" print *, "Total virial is calculated as : (v11, v22, v33, v32, v31, v21)"
print *, calculated print *, calculated
end subroutine end subroutine calc_tot_virial
subroutine calc_tot_energy
integer :: i, inod, ibasis, j
calculated=0
!Sum atom energies
do i = 1, atom_num
calculated(1) = calculated(1) + energy_atom(i)
end do
!Sum element energies
do i =1, ele_num
do inod=1, ng_node(lat_ele(i))
do ibasis=1, basisnum(lat_ele(i))
calculated(1) = calculated(1) + energy_node(ibasis, inod, i)*(size_ele(i)**3)/ng_node(lat_ele(i))
end do
end do
end do
print *, 'Total energy in eV is:'
print *, calculated(1)
end subroutine calc_tot_energy
end module mode_calc end module mode_calc

@ -12,9 +12,9 @@ module mode_create
character(len=100) :: name, element_type character(len=100) :: name, element_type
real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), & real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), &
orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3) orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3), basis_pos(3,10)
integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), & integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
basis_pos(3,10), esize_nums, esize_index(10) esize_nums, esize_index(10)
logical :: dup_flag, dim_flag, efill, crossb(3) logical :: dup_flag, dim_flag, efill, crossb(3)
real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:) real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
@ -186,6 +186,7 @@ module mode_create
character(len=100) :: orient_string character(len=100) :: orient_string
character(len=2) :: btype character(len=2) :: btype
logical :: isortho, isrighthanded, bool logical :: isortho, isrighthanded, bool
real(kind=dp) :: mass
!Pull out all required positional arguments !Pull out all required positional arguments
@ -265,7 +266,8 @@ module mode_create
do i = 1, max_basisnum do i = 1, max_basisnum
call get_command_argument(arg_pos, btype) call get_command_argument(arg_pos, btype)
arg_pos = arg_pos+1 arg_pos = arg_pos+1
call add_atom_type(btype, basis_type(i,1)) call atommass(btype, mass)
call add_atom_type(mass, basis_type(i,1))
do j = 1, 3 do j = 1, 3
call get_command_argument(arg_pos, textholder) call get_command_argument(arg_pos, textholder)
read(textholder, *) basis_pos(j,i) read(textholder, *) basis_pos(j,i)
@ -334,7 +336,8 @@ module mode_create
if (basisnum(1) == 0) then if (basisnum(1) == 0) then
max_basisnum = 1 max_basisnum = 1
basisnum(1) = 1 basisnum(1) = 1
call add_atom_type(name, basis_type(1,1)) !If basis command not defined then we use name as the atom_name call atommass(name, mass)
call add_atom_type(mass, basis_type(1,1)) !If basis command not defined then we use name as the atom_name
basis_pos(:,1) = 0.0_dp basis_pos(:,1) = 0.0_dp
end if end if
@ -371,6 +374,8 @@ module mode_create
numlatpoints = numlatpoints*(bd_in_lat(2*i)-bd_in_lat(2*i-1)) numlatpoints = numlatpoints*(bd_in_lat(2*i)-bd_in_lat(2*i-1))
end do end do
if(numlatpoints < 0) stop "number of lattice points is negative, this can occur if the model is too big"
!Allocate the correct lat variable !Allocate the correct lat variable
select case(esize) select case(esize)
!Atomistics !Atomistics

@ -8,7 +8,9 @@ module mode_da
implicit none implicit none
integer :: num_d_points
integer, allocatable :: is_slipped(:,:,:) integer, allocatable :: is_slipped(:,:,:)
real(kind=dp), allocatable :: disl_line(:,:), disl_data(:)
private :: write_xyz private :: write_xyz
public public
@ -31,16 +33,28 @@ module mode_da
if(allocated(is_slipped)) deallocate(is_slipped) if(allocated(is_slipped)) deallocate(is_slipped)
allocate(is_slipped(max_basisnum, max_ng_node, ele_num)) allocate(is_slipped(max_basisnum, max_ng_node, ele_num))
if(allocated(disl_line)) deallocate(disl_line)
num_d_points=0
allocate(disl_line(3,1024))
allocate(disl_data(1024))
call cga_da call cga_da
!Now create the output file num and write out to xyz format !Now create the output file num and write out to xyz format
ppos = scan(trim(infiles(1)),".", BACK= .true.)
if ( ppos > 0 ) then
outfile = 'da_'//infiles(1)(1:ppos)//'vtk'
else
outfile = 'da_'//infiles(1)//'.vtk'
end if
call write_da_vtk(outfile)
ppos = scan(trim(infiles(1)),".", BACK= .true.) ppos = scan(trim(infiles(1)),".", BACK= .true.)
if ( ppos > 0 ) then if ( ppos > 0 ) then
outfile = 'da_'//infiles(1)(1:ppos)//'xyz' outfile = 'da_'//infiles(1)(1:ppos)//'xyz'
else else
outfile = 'da_'//infiles(1)//'.xyz' outfile = 'da_'//infiles(1)//'.xyz'
end if end if
call write_da_xyz(outfile) call write_da_xyz(outfile)
return return
@ -65,16 +79,39 @@ module mode_da
subroutine cga_da subroutine cga_da
integer :: i, j, k, l, ibasis, nei integer :: i, j, k, l, ibasis, nei, close_neighbors(2,4,6), far_neighbor(4,6), minindex, minindices(2), linevec(3)
integer :: face_types(ele_num*6), face_ele(6*ele_num) integer :: face_types(ele_num*6), face_ele(6*ele_num), diff_pair(2,6)
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_node_dist, ndiff, rmax(3), rmin(3), rnorm max_node_dist, ndiff, rmax(3), rmin(3), rnorm, v1(3), v2(3), v1norm, v2norm, theta1, theta2, &
diff_mag(6), diffvec(3)
real(kind=dp), allocatable :: disl_line_temp(:,:), disl_data_temp(:)
logical :: is_edge
!Initialize variables !Initialize variables
l = 0 l = 0
max_node_dist = 0 max_node_dist = 0
!Now save the close and far neighbors of the nodes. This is done to attempt to figure out the character of the dislocation
!by comparing how the burgers vector is distributed over the face.
do j = 1, 6
close_neighbors(1, 1, j) = cubic_faces(4, j)
close_neighbors(2, 1, j) = cubic_faces(2, j)
far_neighbor(1,j) = cubic_faces(3,j)
close_neighbors(1, 2, j) = cubic_faces(1, j)
close_neighbors(2, 2, j) = cubic_faces(3, j)
far_neighbor(2,j) = cubic_faces(4,j)
close_neighbors(1, 3, j) = cubic_faces(2, j)
close_neighbors(2, 3, j) = cubic_faces(4, j)
far_neighbor(3,j) = cubic_faces(1,j)
close_neighbors(1, 4, j) = cubic_faces(3, j)
close_neighbors(2, 4, j) = cubic_faces(1, j)
far_neighbor(4,j) = cubic_faces(2,j)
end do
!First calculate all of the face centroids !First calculate all of the face centroids
do i = 1, ele_num do i = 1, ele_num
do j = 1, 6 do j = 1, 6
@ -101,8 +138,6 @@ module mode_da
end do end do
!Calculate the largest difference between nodes for each face
end do end do
!Now set the cut off distance as half the distance from opposite nodes in the largest element !Now set the cut off distance as half the distance from opposite nodes in the largest element
@ -131,31 +166,121 @@ module mode_da
end do end do
end do end do
do j = 1, 3 !Now calculate the maximum difference between nodes
vnode(j,1:basisnum(lat_ele(face_ele(i))),:) = vnode(j,1:basisnum(lat_ele(face_ele(i))),:) - & vnorm = 0
minval(vnode(j,1:basisnum(lat_ele(face_ele(i))),:))
end do
do j = 1, 4 do j = 1, 4
do ibasis = 1, basisnum(lat_ele(face_ele(i))) do k=j,4
vnorm(ibasis,j) = norm2(vnode(:, ibasis, j)) v1=vnode(:,1,j) - vnode(:,1,k)
v1norm=norm2(v1)
if (vnorm < v1norm) then
vnorm = v1norm
diffvec= v1
end if
end do end do
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 !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 !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. !I think 0.5 works ok though.
if (any(vnorm > 0.5_dp)) then if (vnorm>0.5_dp) then
print *, "Element number ", face_ele(i), " is dislocated along face ", face_types(i), & print *, "Element number ", face_ele(i), " is dislocated along face ", face_types(i), &
" with neighbor ", face_ele(nei), " with max displacement of ", maxval(vnorm) " with neighbor ", face_ele(nei), " with max displacement of ", vnorm
l=0
do j = 1, 4 do j = 1, 4
is_slipped(:, cubic_faces(j,face_types(i)), face_ele(i)) = 1 is_slipped(:, cubic_faces(j,face_types(i)), face_ele(i)) = 1
print *, j, vnode(:,1,j)
!This portion of the code is used to determine what the character of the dislocation most likely is
!effectively how this works is we look for the two nodes with the most similar burgers vector.
!The vector between these two nodes is close to the line direction and as a result we can probably estimate
!if it's close to edge, screw, or if it's pretty fairly mixed.
do k = j+1, 4
l=l+1
diff_mag(l) = norm2(vnode(:, 1, j)-vnode(:,1,k))
diff_pair(1,l)=cubic_faces(j, face_types(i))
diff_pair(2,l)=cubic_faces(k, face_types(i))
end do
end do
minindex = minloc(diff_mag,1)
!Now figure out if the min difference between nodes is associated with a 112 direction
is_edge = .false.
do j = 1,6
if((diff_pair(1,minindex) == oneonetwopairs(1,j)).and.(diff_pair(2,minindex)==oneonetwopairs(2,j))) then
is_edge=.true.
exit
else if((diff_pair(2,minindex)==oneonetwopairs(1,j)).and.(diff_pair(1,minindex)==oneonetwopairs(2,j)))then
is_edge=.true.
exit
end if
end do end do
if(is_edge) then
print *, 'Dislocation has primarily edge character'
else
print *, 'Dislocation has primarily screw character'
end if
if( i < nei) then
num_d_points=num_d_points+2
if(num_d_points > size(disl_line,2)) then
allocate(disl_line_temp(3,size(disl_line,2)+1024))
disl_line_temp=0.0_dp
disl_line_temp(:, 1:size(disl_line,2))=disl_line
call move_alloc(disl_line_temp, disl_line)
end if
if(num_d_points/2 > size(disl_data)) then
allocate(disl_data_temp(size(disl_data)+1024))
disl_data_temp=0
disl_data_temp(1:size(disl_data)) = disl_data
call move_alloc(disl_data_temp, disl_data)
end if
print *, num_d_points/2, vnorm
disl_data(num_d_points/2)=vnorm
disl_line(:,num_d_points-1) =r_node(:, 1, diff_pair(1,minindex), face_ele(i))
disl_line(:,num_d_points) =r_node(:, 1, diff_pair(2,minindex), face_ele(i))
end if
end if end if
end if end if
end do end do
end subroutine cga_da end subroutine cga_da
subroutine write_da_vtk(outfile)
character (len=*), intent(in) :: outfile
integer :: i
open(unit=11, file=outfile, action='write', status='replace', position='rewind')
write(11, '(a)') '# vtk DataFile Version 4.0.1'
write(11, '(a)') 'Dislocation line analyzed via cacmb dislocation analysis'
write(11, '(a)') 'ASCII'
write(11, '(a)') 'DATASET UNSTRUCTURED_GRID'
write(11, *) 'POINTS', num_d_points, 'float'
!Write the points to the file
do i = 1, num_d_points
write(11, '(3f23.15)') disl_line(:, i)
end do
write(11,*) 'CELLS', num_d_points/2, num_d_points/2*3
do i=1, num_d_points, 2
write(11,'(3i16)') 2, i-1, i
end do
write(11,*) 'CELL_TYPES', num_d_points/2
do i=1, num_d_points, 2
write(11,'(i16)') 3
end do
write(11,*) 'CELL_DATA', num_d_points/2
write(11, '(a)') 'SCALARS burgers_vec, float'
write(11, '(a)') 'LOOKUP_TABLE default'
do i=1, num_d_points, 2
write(11, '(f23.15)') disl_data((i+1)/2)
end do
close(11)
end subroutine write_da_vtk
subroutine write_da_xyz(outfile) subroutine write_da_xyz(outfile)
!This subroutine write the element positions to a .xyz file and marks whether they are slipped or not !This subroutine write the element positions to a .xyz file and marks whether they are slipped or not
character(len=*), intent(in) :: outfile character(len=*), intent(in) :: outfile
@ -184,4 +309,4 @@ module mode_da
close(11) close(11)
end subroutine write_da_xyz end subroutine write_da_xyz
end module end module mode_da

@ -10,7 +10,7 @@ module mode_metric
integer :: nfiles integer :: nfiles
character(len=100) :: metric_type character(len=100) :: metric_type
real(kind=dp) :: rc_off real(kind=dp) :: rc_off, old_box_len(3)
!Save reference positions !Save reference positions
integer :: np, npreal, nmet integer :: np, npreal, nmet
@ -57,6 +57,10 @@ module mode_metric
print *, "Getting initial neighbor list" print *, "Getting initial neighbor list"
call calc_neighbor(rc_off, r_zero, np) call calc_neighbor(rc_off, r_zero, np)
!Save box lengths for applying periodic boundaries
do i=1,3
old_box_len(i) = box_bd(2*i) - box_bd(2*i-1)
end do
!Reset element and box !Reset element and box
call reset_data call reset_data
call reset_box call reset_box
@ -128,8 +132,11 @@ module mode_metric
!This subroutine calculates the continuum metric that we require !This subroutine calculates the continuum metric that we require
integer :: i, j, k, nei, ip, jp integer :: i, j, k, nei, ip, jp
real(kind=dp) :: def_grad(3,3), omega(3,3), eta(3,3), rij(3), eta_inv(3,3), ftf(3,3), & real(kind=dp) :: def_grad(3,3), omega(3,3), eta(3,3), rij(3), eta_inv(3,3), ftf(3,3), &
U(3,3), R(3,3), Rskew(3,3), oldrij(3) U(3,3), R(3,3), Rskew(3,3), oldrij(3), box_len(3)
do i = 1, 3
box_len(i) = box_bd(2*i) - box_bd(2*i-1)
end do
!Loop over all points !Loop over all points
do ip = 1, np do ip = 1, np
eta(:,:) = 0.0_dp eta(:,:) = 0.0_dp
@ -138,8 +145,17 @@ module mode_metric
do jp = 1, nei_num(ip) do jp = 1, nei_num(ip)
!Calculate the neighbor vec in current configuration !Calculate the neighbor vec in current configuration
nei = nei_list(jp, ip) nei = nei_list(jp, ip)
rij = r_curr(:,nei) - r_curr(:,ip) rij = r_curr(:,nei) - r_curr(:,ip) +periodvec(:,jp,ip)*box_len
oldrij = r_zero(:,nei) - r_zero(:,ip) if (norm2(rij) > 5*rc_off) then
do j = 1, 3
if (r_curr(j,nei) > r_curr(j,ip)) then
rij(j) = rij(j) - box_len(j)
else if(r_curr(j,nei) < r_curr(j,ip)) then
rij(j) = rij(j) + box_len(j)
end if
end do
end if
oldrij = r_zero(:,nei) - r_zero(:,ip) + periodvec(:,jp,ip)*old_box_len
!Calculate eta and omega !Calculate eta and omega
do i = 1,3 do i = 1,3

@ -4,10 +4,12 @@ module neighbors
use elements use elements
use subroutines use subroutines
use functions use functions
use box
integer, allocatable :: nei_list(:,:), nei_num(:), nn(:) integer, allocatable :: nei_list(:,:), nei_num(:), nn(:), periodvec(:,:,:)
real(kind=dp), allocatable :: init_vec(:,:,:), output(:), microrotation(:,:) real(kind=dp), allocatable :: init_vec(:,:,:), output(:), microrotation(:,:)
public public
contains contains
subroutine build_cell_list(numinlist, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell) subroutine build_cell_list(numinlist, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell)
@ -80,44 +82,70 @@ module neighbors
real(kind=dp), intent(in) :: rc_off real(kind=dp), intent(in) :: rc_off
integer, intent(in) :: n integer, intent(in) :: n
real(kind=dp), dimension(3,n) :: r_list real(kind=dp), dimension(:, :), intent(in) :: r_list
integer :: i, c(3), ci, cj, ck, num_nei, nei integer :: i, j, c(3),cn(3), ci, cj, ck, num_nei, nei, v(3), period_dir(3)
!Variables for cell list code !Variables for cell list code
integer, dimension(3) ::cell_num integer, dimension(3) ::cell_num
integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:) integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:)
integer :: which_cell(3,n) integer :: which_cell(3,n)
real(kind=dp) :: r(3), box_len(3)
logical :: period_bd(3)
!First reallocate the neighbor list codes !First reallocate the neighbor list variables
if (allocated(nei_list)) then if (allocated(nei_list)) then
deallocate(nei_list,nei_num) deallocate(nei_list,nei_num, periodvec)
end if end if
allocate(nei_list(100,n),nei_num(n)) allocate(nei_list(100,n),nei_num(n), periodvec(3,100,n))
nei_list=0
periodvec=0
nei_num=0
!Now first pass the position list and and point num to the cell algorithm !Now first pass the position list and and point num to the cell algorithm
call build_cell_list(n, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell) call build_cell_list(n, r_list, rc_off, cell_num, num_in_cell, cell_list, which_cell)
do i=1, 3
if (box_bc(i:i) == 'p') then
period_bd(i) = .true.
else
period_bd(i)=.false.
end if
box_len(i) = box_bd(2*i) - box_bd(2*i-1)
end do
!Now loop over every point and find it's neighbors !Now loop over every point and find it's neighbors
pointloop: do i = 1, n pointloop: do i = 1, n
!First check to see if the point is a filler point, if so then skip it !First check to see if the point is a filler point, if so then skip it
if(r_list(1,i) < box_bd(1)) cycle if(r_list(1,i) < box_bd(1)) cycle
!c is the position of the cell that the point !c is the position of the cell that the point belongs to
c = which_cell(:,i) c = which_cell(:,i)
!Loop over all neighboring cells !Loop over all neighboring cells
do ci = -1, 1, 1 do ci = -1, 1, 1
do cj = -1, 1, 1 do cj = -1, 1, 1
do ck = -1, 1, 1 do ck = -1, 1, 1
!First check to make sure that the neighboring cell exists !First try to apply periodic boundaries
v=(/ ck, cj, ci /)
cn=0
period_dir=0
do j=1, 3
if ((c(j) + v(j) == 0).and.period_bd(j)) then
cn(j) = cell_num(j)
period_dir(j) = 1
else if ((c(j) + v(j) > cell_num(j)).and.period_bd(j)) then
cn(j) = 1
period_dir(j) = -1
else if ((c(j)+v(j) >= 1) .and. (c(j)+v(j) <= cell_num(j))) then
cn(j) = c(j) + v(j)
end if
end do
if(any((c + (/ ck, cj, ci /)) == 0)) cycle 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) do num_nei = 1, num_in_cell(cn(1),cn(2),cn(3))
nei = cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci) nei = cell_list(num_nei,cn(1),cn(2),cn(3))
!Check to make sure the atom isn't the same index as the atom we are checking !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 !and that the neighbor hasn't already been deleted
@ -125,10 +153,13 @@ module neighbors
!Now check to see if it is in the cutoff radius, if it is add it to the neighbor list for that !Now check to see if it is in the cutoff radius, if it is add it to the neighbor list for that
!atom and calculate the initial neighbor vector !atom and calculate the initial neighbor vector
if (norm2(r_list(:,nei)-r_list(:,i)) < rc_off) then
r = r_list(:,nei) + period_dir*box_len
if (norm2(r-r_list(:,i)) < rc_off) then
nei_num(i) = nei_num(i) + 1 nei_num(i) = nei_num(i) + 1
nei_list(nei_num(i), i) = nei nei_list(nei_num(i), i) = nei
periodvec(:,nei_num(i),i) = period_dir
end if end if
end if end if

@ -1,6 +1,7 @@
module opt_bubble module opt_bubble
!This module contains the bubble option which can be used to create voids with specific pressures of certain atoms !This module contains the bubble option which can be used to create voids with specific pressures of certain atoms
use atoms
use parameters use parameters
use elements use elements
use box use box
@ -8,7 +9,7 @@ module opt_bubble
implicit none implicit none
real(kind=dp), private :: br, bt, bp, c(3) real(kind=dp), private :: br, brat, c(3)
character(len=2), private :: species character(len=2), private :: species
public public
@ -16,8 +17,8 @@ module opt_bubble
subroutine bubble(arg_pos) subroutine bubble(arg_pos)
integer, intent(inout) :: arg_pos integer, intent(inout) :: arg_pos
integer :: new_type, n, j, i integer :: new_type, n, j, i, atom_num_pre
real(kind = dp) :: p(3), rand, factor, per, vol real(kind = dp) :: p(3), rand, factor, per, vol, mass
print *, '------------------------------------------------------------' print *, '------------------------------------------------------------'
print *, 'Option Bubble' print *, 'Option Bubble'
@ -28,34 +29,36 @@ module opt_bubble
!Now we use the existing group code to delete a sphere which represents the bubble !Now we use the existing group code to delete a sphere which represents the bubble
centroid=c centroid=c
radius = br radius = br
type='both' type='all'
gshape='sphere' gshape='sphere'
group_nodes = .true. group_nodes = .true.
group_atom_types=0
call get_group call get_group
call refine_group call refine_group
call get_group call get_group
atom_num_pre= atom_num
call delete_group call delete_group
!Now we create a new atom type with the desired species !Now we create a new atom type with the desired species
call add_atom_type(species, new_type) call atommass(species, mass)
call add_atom_type(mass, new_type)
!Now we calculate the number of atoms we need to add for the desired pressure !Now we calculate the number of atoms we need to add for the desired pressure
print *, "Creating a bubble with pressure", bp, " at temperature ", bt, " with radius ", br !print *, "Creating a bubble with pressure", bp, " at temperature ", bt, " with radius ", br
!
factor = 1.0e24/6.02214e23 !factor = 1.0e24/6.02214e23
if (bp <= 10.0) then !if (bp <= 10.0) then
per=factor*(3.29674113+4.51777872*bp**(-0.80473167)) ! per=factor*(3.29674113+4.51777872*bp**(-0.80473167))
else if (bp .le. 20.0) then !else if (bp .le. 20.0) then
per=factor*(4.73419689-0.072919483*bp) ! per=factor*(4.73419689-0.072919483*bp)
else !else
per=factor*(4.73419689-0.072919483*bp) ! per=factor*(4.73419689-0.072919483*bp)
print *, 'warning: pressure is too high' ! print *, 'warning: pressure is too high'
print *, 'equation of state is only valid for < 20 GPa' ! print *, 'equation of state is only valid for < 20 GPa'
endif !endif
vol = 4.0*pi/3.0*br**3.0 !vol = 4.0*pi/3.0*br**3.0
n = vol/per n = brat*(atom_num_pre-atom_num)
print *, "Adding ", n, " atoms of species ", species print *, "Adding ", n, " atoms of species ", species
!Now add n atoms randomly within the sphere !Now add n atoms randomly within the sphere
@ -77,39 +80,66 @@ module opt_bubble
integer, intent(inout) :: arg_pos integer, intent(inout) :: arg_pos
integer :: i, arglen integer :: i, arglen
real(kind=dp) :: mass
character(len=100) :: tmptxt character(len=100) :: tmptxt
!First read in the bubble centroid !First read in the bubble centroid
do i = 1, 3 do i = 1, 3
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
call get_command_argument(arg_pos, tmptxt, arglen) call get_command_argument(arg_pos, tmptxt, arglen)
print *, trim(adjustl(tmptxt))
call parse_pos(i, tmptxt, c(i)) call parse_pos(i, tmptxt, c(i))
end do end do
!Now the bubble radius !Now the bubble radius
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
call get_command_argument(arg_pos, tmptxt, arglen) call get_command_argument(arg_pos, tmptxt, arglen)
print *, trim(adjustl(tmptxt))
if(arglen == 0) stop "Missing bubble radius" if(arglen == 0) stop "Missing bubble radius"
read(tmptxt, *) br read(tmptxt, *) br
!Now bubble pressure !Now bubble ratio
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, tmptxt, arglen)
if(arglen == 0) stop "Missing bubble pressure"
read(tmptxt, *) bp
!Now bubble temperature
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
call get_command_argument(arg_pos, tmptxt, arglen) call get_command_argument(arg_pos, tmptxt, arglen)
if(arglen == 0) stop "Missing bubble temperature" print *, trim(adjustl(tmptxt))
read(tmptxt, *) bt if(arglen == 0) stop "Missing bubble ratio"
read(tmptxt, *) brat
!Now the bubble species !Now the bubble species
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
call get_command_argument(arg_pos, species, arglen) call get_command_argument(arg_pos, species, arglen)
print *, species
if(arglen == 0) stop "Missing bubble species" if(arglen == 0) stop "Missing bubble species"
!OPtional arguments
do while(.true.)
if(arg_pos > command_argument_count()) exit
arg_pos=arg_pos+1 arg_pos=arg_pos+1
call get_command_argument(arg_pos, tmptxt)
tmptxt=adjustl(tmptxt)
select case(trim(tmptxt))
case('excludetypes')
arg_pos=arg_pos + 1
call get_command_argument(arg_pos, tmptxt, arglen)
if(arglen==0) stop "Missing number of atoms to exclude"
read(tmptxt, *) exclude_num
do i=1,exclude_num
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, tmptxt, arglen)
if(arglen==0) stop "Missing exclude atom types"
call atommass(tmptxt, mass)
call add_atom_type(mass, exclude_types(i))
end do
case default
exit
end select
end do
return return
end subroutine parse_bubble end subroutine parse_bubble
end module opt_bubble end module opt_bubble

@ -87,7 +87,7 @@ module opt_delete
allocate(which_cell(3,atom_num)) allocate(which_cell(3,atom_num))
!First pass the atom list and atom num to the algorithm which builds the cell list !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) call build_cell_list(atom_num, r_atom, 5*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 !Now loop over every atom and figure out if it has neighbors within the rc_off
delete_num = 0 delete_num = 0

@ -520,7 +520,7 @@ module opt_disl
end do end do
return return
end subroutine end subroutine disloop
!******************************************************** !********************************************************
! SOLIDANGLE ! SOLIDANGLE

@ -2,6 +2,7 @@ module opt_group
!This module contains all code associated with dislocations !This module contains all code associated with dislocations
use atoms
use parameters use parameters
use elements use elements
use subroutines use subroutines
@ -11,7 +12,7 @@ module opt_group
implicit none implicit none
integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize, insert_type, & integer :: group_ele_num, group_atom_num, remesh_size,normal, dim1, dim2, random_num, group_type, notsize, insert_type, &
insert_site, num_species insert_site, num_species, group_atom_types, exclude_num, exclude_types(10)
character(len=15) :: type, gshape!Type indicates what element type is selected and shape is the group shape character(len=15) :: type, gshape!Type indicates what element type is selected and shape is the group shape
character(len=2) :: species_type(10) character(len=2) :: species_type(10)
real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness, insert_conc, & real(kind=dp) :: block_bd(6), centroid(3), vertices(3,3),disp_vec(3), radius, bwidth, shell_thickness, insert_conc, &
@ -38,6 +39,7 @@ module opt_group
insert_type = 0 insert_type = 0
random_num=0 random_num=0
group_type=0 group_type=0
group_atom_types=0
notsize=0 notsize=0
displace=.false. displace=.false.
delete=.false. delete=.false.
@ -46,6 +48,7 @@ module opt_group
flip = .false. flip = .false.
group_nodes=.false. group_nodes=.false.
num_species = 0 num_species = 0
exclude_num = 0
if(allocated(element_index)) deallocate(element_index) if(allocated(element_index)) deallocate(element_index)
if(allocated(atom_index)) deallocate(atom_index) if(allocated(atom_index)) deallocate(atom_index)
@ -86,10 +89,10 @@ module opt_group
call displace_group call displace_group
end if end if
if(insert_type > 0) then if(insert_type > 0) then
print *, "Insert command has been dropped" !print *, "Insert command has been dropped"
stop 1 !stop 1
call get_group call get_group
! call insert_group call insert_group
end if end if
if(num_species > 0) then if(num_species > 0) then
@ -106,7 +109,7 @@ module opt_group
integer :: i, j, arglen, in_num integer :: i, j, arglen, in_num
character(len=100) :: textholder, type_spec character(len=100) :: textholder, type_spec
real(kind=dp) H real(kind=dp) H, mass
!Parse type and shape command !Parse type and shape command
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
@ -469,6 +472,24 @@ module opt_group
refine=.true. refine=.true.
case('refinefill') case('refinefill')
refinefill=.true. refinefill=.true.
case('selecttype')
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) stop "Missing select type"
call atommass(textholder, mass)
call add_atom_type(mass, group_atom_types)
case('excludetypes')
arg_pos=arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if(arglen==0) stop "Missing number of atoms to exclude"
read(textholder, *) exclude_num
do i=1,exclude_num
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if(arglen==0) stop "Missing exclude atom types"
call atommass(textholder, mass)
call add_atom_type(mass, exclude_types(i))
end do
case('remesh') case('remesh')
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen) call get_command_argument(arg_pos, textholder, arglen)
@ -493,7 +514,8 @@ module opt_group
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen) call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) stop "Missing atom type for group" if (arglen==0) stop "Missing atom type for group"
call add_atom_type(textholder, group_type) call atommass(textholder, mass)
call add_atom_type(mass, group_type)
case('notsize') case('notsize')
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen) call get_command_argument(arg_pos, textholder, arglen)
@ -521,7 +543,8 @@ module opt_group
arg_pos=arg_pos+1 arg_pos=arg_pos+1
call get_command_argument(arg_pos, textholder, arglen) call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) stop "Missing element type for insert command" if (arglen==0) stop "Missing element type for insert command"
call add_atom_type(textholder, insert_type) call atommass(textholder, mass)
call add_atom_type(mass, insert_type)
arg_pos=arg_pos + 1 arg_pos=arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen) call get_command_argument(arg_pos, textholder, arglen)
select case(trim(adjustl(textholder))) select case(trim(adjustl(textholder)))
@ -616,7 +639,7 @@ module opt_group
select case(trim(adjustl(type))) select case(trim(adjustl(type)))
case('atoms','all') case('atoms','all')
do i = 1, atom_num do i = 1, atom_num
if(in_group(r_atom(:,i)).neqv.flip) then if(in_group(type_atom(i), r_atom(:,i)).neqv.flip) then
group_atom_num = group_atom_num + 1 group_atom_num = group_atom_num + 1
if (group_atom_num > size(atom_index)) then if (group_atom_num > size(atom_index)) then
allocate(resize_array(size(atom_index) + 1024)) allocate(resize_array(size(atom_index) + 1024))
@ -643,15 +666,6 @@ module opt_group
end if end if
end select end select
j = 0
do i = 1, group_atom_num
if (atom_index(i) == 23318348) then
j = j + 1
end if
end do
if (j > 1) stop "Code broken"
print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms." print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms."
end subroutine get_group end subroutine get_group
@ -693,8 +707,19 @@ module opt_group
!This command is used to refine the group to full atomistics !This command is used to refine the group to full atomistics
integer :: i, j, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num 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) real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), data_interp(10, max_basisnum*max_esize**3), &
vel_interp(3, max_basisnum*max_esize**3)
real(kind=dp), allocatable :: vel_tmp(:,:)
logical :: refine_dat
if (allocated(force_node)) then
refine_dat=.true.
else
refine_dat=.false.
end if
print *, size(data_interp)
!Refining to atoms !Refining to atoms
if(group_ele_num > 0) then if(group_ele_num > 0) then
orig_atom_num = atom_num orig_atom_num = atom_num
@ -704,17 +729,38 @@ module opt_group
do i = 1, group_ele_num do i = 1, group_ele_num
ie = element_index(i) ie = element_index(i)
!Get the interpolated atom positions !Get the interpolated atom positions
if (refine_dat) then
call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp, &
energy_node(:,:,i), force_node(:,:,:,i), virial_node(:,:,:,i), data_interp)
else
call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp)
end if
if(allocated(vel_atom)) then
call interpolate_vel(type_ele(i), size_ele(i), lat_ele(i), vel_node(:,:,:,i), vel_interp)
end if
!Loop over all interpolated atoms and add them to the system, we apply periodic boundaries !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 !here as well to make sure they are in the box
do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3 do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3
call apply_periodic(r_interp(:,j)) call apply_periodic(r_interp(:,j))
call add_atom(0,type_interp(j), r_interp(:,j)) call add_atom(0,type_interp(j), r_interp(:,j))
if(refine_dat) then
call add_atom_data(atom_num, data_interp(1, j), data_interp(2:4, j), data_interp(5:10, j))
end if
if(allocated(vel_atom)) then
if(size(vel_atom, 2) < size(type_atom)) then
allocate(vel_tmp(3, size(type_atom)))
vel_tmp=0.0_dp
vel_tmp(:, 1:size(vel_atom,2)) = vel_atom
call move_alloc(vel_tmp, vel_atom)
end if
vel_atom(:, atom_num) = vel_interp(:,j)
end if
end do end do
end do end do
!Once all atoms are added we delete all of the elements !Once all atoms are added we delete all of the elements
call delete_elements(group_ele_num, element_index) call delete_elements(group_ele_num, element_index)
print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms." print *, group_ele_num, " elements of group are refined to ", atom_num -orig_atom_num, " atoms."
end if end if
end subroutine refine_group end subroutine refine_group
@ -751,7 +797,7 @@ module opt_group
ravg = ravg + rfill(:,ibasis, 1)/basisnum(lat_ele(ie)) ravg = ravg + rfill(:,ibasis, 1)/basisnum(lat_ele(ie))
end do end do
if( in_group(ravg)) then if( in_group(0,ravg)) then
nump_ele = nump_ele - 1 nump_ele = nump_ele - 1
end if end if
end do end do
@ -1124,86 +1170,102 @@ module opt_group
end subroutine change_group_type end subroutine change_group_type
! 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 !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. !Coarse-graining methodology which doesn't allow for concentration fluctuations.
! real(kind=dp) interstitial_sites(3,14), rand, rinsert(3) real(kind=dp) interstitial_sites(3,14), rand, rinsert(3)
! integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi
! integer, allocatable :: used_sites(:,:) integer, allocatable :: used_sites(:,:)
!
! !First save all of the displacement vectors from a lattice site to interstitial site !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 !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, & interstitial_sites= reshape( (/ -0.5_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.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.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, &
! 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)) shape(interstitial_sites))
!
! !First we calculate the number of atoms needed based on the number of atoms in the group and the concentration !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 interstitial_sites=interstitial_sites*insert_lattice
!
! add_num = (insert_conc*group_atom_num)/(1-insert_conc)
! allocate(used_sites(2,add_num)) !Now set up the random number generator for the desired interstitial type
! select case(insert_site)
! print *, "Inserting ", add_num, " atoms as atom type ", insert_type case(1)
! rlo=1
! !Now set up the random number generator for the desired interstitial type rhi=6
! select case(insert_site) case(2)
! case(1) rlo=7
! rlo=1 rhi = 14
! rhi=6 case(3)
! case(2) rlo=1
! rlo=7 rhi=14
! rhi = 14 end select
! case(3)
! rlo=1 if ((insert_conc > 1).or.(insert_conc < 0)) stop "Insert concentration should be between 0 and 1"
! rhi=14 add_num = (insert_conc*group_atom_num)
! end select allocate(used_sites(2,add_num))
!
! !Now add the atoms print *, "Inserting ", add_num, " atoms as atom type ", insert_type
! i = 1
! addloop:do while ( i < add_num) !Now add the atoms
! call random_number(rand) i = 1
! rindex = int(1+rand*(group_atom_num-1)) addloop:do while ( i < add_num)
! ia=atom_index(rindex) call random_number(rand)
! call random_number(rand) rindex = int(1+rand*(group_atom_num-1))
! sindex = int(rlo+rand*(rhi-rlo)) ia=atom_index(rindex)
! do j = 1, i call random_number(rand)
! if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop sindex = int(rlo+rand*(rhi-rlo))
! end do do j = 1, i
! rinsert = r_atom(:,ia) + matmul(sub_box_ori(:,:,sbox_atom(ia)),interstitial_sites(:,sindex)) if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop
! sbox = sbox_atom(ia) end do
! call add_atom(0, insert_type, sbox, rinsert) rinsert = r_atom(:,ia) + matmul(box_ori,interstitial_sites(:,sindex))
! used_sites(1,i) = ia call add_atom(0, insert_type, rinsert)
! used_sites(2,i) = sindex used_sites(1,i) = ia
! i = i + 1 used_sites(2,i) = sindex
! end do addloop i = i + 1
! end subroutine insert_group end do addloop
end subroutine insert_group
subroutine alloy_group subroutine alloy_group
!This subroutine randomizes the atom types to reach desired concentrations, this only operates on atoms !This subroutine randomizes the atom types to reach desired concentrations, this only operates on atoms
integer :: i, j, ia, type_map(10), added_types(num_species) integer :: i, j, ia, type_map(10), added_types(num_species)
real(kind=dp) :: rand real(kind=dp) :: rand, mass, old_fractions(num_species)
character(len=2) :: sorted_species_type(10)
logical :: species_mapped(10)
print *, "Alloying group with desired fractions", s_fractions(1:num_species) print *, "Alloying group with desired fractions", s_fractions(1:num_species)
old_fractions=s_fractions(1:num_species)
!First we sort the fractions !First we sort the fractions
call quicksort(s_fractions(1:num_species)) call quicksort(s_fractions(1:num_species))
species_mapped = .false.
do i = 1,num_species
do j=1, num_species
if (is_equal(s_fractions(i), old_fractions(j)) .and. .not. species_mapped(j)) then
sorted_species_type(i)=species_type(j)
species_mapped(j) = .true.
exit
end if
end do
end do
!Now get the atom type maps for all the atoms and make the fractions a running sum !Now get the atom type maps for all the atoms and make the fractions a running sum
do i = 1, num_species do i = 1, num_species
call add_atom_type(species_type(i), type_map(i)) call atommass(sorted_species_type(i), mass)
call add_atom_type(mass, type_map(i))
if(i > 1) s_fractions(i) = s_fractions(i) + s_fractions(i-1) if(i > 1) s_fractions(i) = s_fractions(i) + s_fractions(i-1)
end do end do
!Now randomize the atom types !Now randomize the atom types
@ -1225,8 +1287,9 @@ module opt_group
return return
end subroutine alloy_group end subroutine alloy_group
function in_group(r) function in_group(atype, r)
!This subroutine determines if a point is within the group boundaries !This subroutine determines if a point is within the group boundaries
integer, intent(in) :: atype
real(kind=dp), intent(in) :: r(3) real(kind=dp), intent(in) :: r(3)
real(kind=dp) :: rnorm real(kind=dp) :: rnorm
@ -1271,6 +1334,19 @@ module opt_group
case('all') case('all')
in_group = .true. in_group = .true.
end select end select
if(group_atom_types> 0) then
in_group = in_group.and.(atype== group_atom_types)
elseif(exclude_num > 0) then
if(in_group) then
do i=1,exclude_num
if (atype == exclude_types(i)) then
in_group=.false.
exit
end if
end do
end if
end if
end function in_group end function in_group
function in_group_ele(esize, elat, rn) function in_group_ele(esize, elat, rn)
@ -1294,13 +1370,13 @@ module opt_group
r_center(:)= r_center(:) + rn(:,ibasis,inod)/basisnum(elat) r_center(:)= r_center(:) + rn(:,ibasis,inod)/basisnum(elat)
end do end do
if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(size_ele(i)/=notsize)) then if((in_group(0,rn(:, ibasis, inod)).neqv.flip).and.(size_ele(i)/=notsize)) then
node_in_out(inod) = 2 node_in_out(inod) = 2
exit nodeloop exit nodeloop
end if end if
gshape ='sphere' gshape ='sphere'
if((in_group(rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then if((in_group(0,rn(:, ibasis, inod)).neqv.flip).and.(esize/=notsize)) then
node_in_out(inod) = 1 node_in_out(inod) = 1
else else
node_in_out(inod) = 0 node_in_out(inod) = 0
@ -1322,7 +1398,7 @@ module opt_group
end do end do
end do end do
if ((in_group(r_center).neqv.flip).and.(esize/= notsize)) then if ((in_group(0,r_center).neqv.flip).and.(esize/= notsize)) then
in_group_ele=.true. in_group_ele=.true.
return return
end if end if
@ -1331,7 +1407,7 @@ module opt_group
r_center(:) = 0.0_dp r_center(:) = 0.0_dp
do inod = 1, ng_node(elat) do inod = 1, ng_node(elat)
do ibasis = 1, basisnum(elat) do ibasis = 1, basisnum(elat)
if ((in_group(rn(:,ibasis,inod)).neqv.flip).and.(esize/=notsize)) then if ((in_group(0,rn(:,ibasis,inod)).neqv.flip).and.(esize/=notsize)) then
in_group_ele=.true. in_group_ele=.true.
return return
end if end if

@ -17,11 +17,13 @@ module opt_slip_plane
!Main calling function for the slip_plane option !Main calling function for the slip_plane option
integer, intent(inout) :: arg_pos integer, intent(inout) :: arg_pos
integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), ibasis integer :: ie, ia, slip_enum, old_atom_num, esize, new_ele_num, n, m, o, ele(3,8), nump_ele, inod, vlat(3), ibasis,&
starting_anum
integer, allocatable :: slip_eles(:), temp_int(:) integer, allocatable :: slip_eles(:), temp_int(:)
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node), ratom(3,max_basisnum), & real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), rfill(3,max_basisnum, max_ng_node), ratom(3,max_basisnum), &
maxp, minp maxp, minp, vel_interp(3, max_basisnum*max_esize**3)
real(kind=dp), allocatable :: vel_tmp(:,:)
integer :: type_interp(max_basisnum*max_esize**3) integer :: type_interp(max_basisnum*max_esize**3)
logical :: lat_points(max_esize,max_esize, max_esize) logical :: lat_points(max_esize,max_esize, max_esize)
@ -58,11 +60,25 @@ module opt_slip_plane
!If we aren't efilling then just refine the element !If we aren't efilling then just refine the element
if(.not.efill) then if(.not.efill) then
starting_anum=atom_num
call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) 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 do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3
call apply_periodic(r_interp(:,ia)) call apply_periodic(r_interp(:,ia))
call add_atom(0, type_interp(ia), r_interp(:,ia)) call add_atom(0, type_interp(ia), r_interp(:,ia))
end do end do
if(allocated(vel_atom)) then
call interpolate_vel(type_ele(ie), size_ele(ie), lat_ele(ie), vel_node(:,:,:,ie), vel_interp)
if(size(vel_atom,2) < size(type_atom)) then
allocate(vel_tmp(3, size(type_atom)))
vel_tmp=0.0_dp
vel_tmp(:, 1:size(vel_atom,2)) = vel_atom
call move_alloc(vel_tmp, vel_atom)
end if
do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3
vel_atom(:, starting_anum+ia) = vel_interp(:, ia)
end do
end if
!If we are efilling then the code is slightly more complex !If we are efilling then the code is slightly more complex
else else
!First populate the lat points array !First populate the lat points array

Loading…
Cancel
Save