Compare commits

...

10 Commits

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

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

@ -2,6 +2,7 @@ module elements
!This module contains the elements data structures, structures needed for building regions
!and operations that are done on elements
use atoms
use parameters
use functions
use subroutines
@ -31,20 +32,21 @@ module elements
!Mapping atom type to provided name
character(len=2), dimension(10) :: type_to_name
integer :: atom_types = 0
real(kind=dp) :: masses(10)
!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), &
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.
!We currently have a limit of 10 lattice types for simplicities sake but this can be easily increased.
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_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.
!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
@ -104,6 +106,14 @@ module elements
cubic_faces(:,5) = (/ 4, 3, 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
fcc_cell = reshape((/ 0.0_dp, 0.0_dp, 0.0_dp, &
0.5_dp, 0.5_dp, 0.0_dp, &
@ -345,9 +355,9 @@ module elements
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
character(len=2), intent(in) :: type
real(kind=dp), intent(in) :: mass
integer, intent(out) :: inttype
logical, optional, intent(in) :: force_new
@ -364,7 +374,7 @@ module elements
exists = .false.
if(.not.force) then
do i=1, 10
if(type == type_to_name(i)) then
if(is_equal(mass, masses(i))) then
exists = .true.
inttype = i
exit
@ -378,7 +388,8 @@ module elements
print *, "Defined atom types are greater than 10 which is currently not supported."
stop 3
end if
type_to_name(atom_types) = type
call atommassspecies(mass, type_to_name(atom_types))
masses(atom_types) = mass
inttype = atom_types
end if
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), optional, intent(in) :: eng(:,:), f(:,:,:), &
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
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) :: t, s, r
@ -451,16 +462,15 @@ module elements
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)-(esize-1)/2)/(1.0_dp*(esize-1)/2)
t=-1.0_dp +(it-1)*(2.0_dp/(esize-1))
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
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)
do ibasis = 1, bnum
@ -518,6 +528,68 @@ module elements
return
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)
!Shape function for rhombohedral elements
real(kind=8), intent(in) :: r, s, t
@ -651,6 +723,7 @@ module elements
end do
end subroutine wrap_atoms
subroutine wrap_elements
integer :: i, inod, ibasis, j, node_in_bd(3,max_ng_node)
real(kind =dp) :: box_len(3)
@ -1008,11 +1081,36 @@ do i = 1, atom_num
subroutine add_atom_data(ia, eng, force, virial)
!Function which sets the atom data arrays
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
force_atom(:,ia) = force(:)
virial_atom(:,ia) = virial(:)
vflag = .true.
return
end subroutine add_atom_data

@ -168,11 +168,11 @@ END FUNCTION StrDnCase
do i =1 ,3
!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.
exit
!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.
exit
end if

@ -139,7 +139,7 @@ module io
do i = 1, ele_num
do inod = 1, ng_node(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
end do
end do
@ -151,7 +151,7 @@ module io
!Write atom positions
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
end do
@ -207,7 +207,7 @@ module io
write(11, '(a)') 'Atoms'
write(11, '(a)') ' '
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
!Write refined element atomic positions
@ -219,7 +219,7 @@ module io
do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
interp_num = interp_num+1
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 select
end do
@ -277,12 +277,12 @@ module io
if (write_dat) then
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)
end do
else
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 if
@ -294,12 +294,12 @@ module io
do inod =1, ng_node(lat_ele(i))
do ibasis=1, basisnum(lat_ele(i))
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), &
virial_node(:, ibasis, inod, i)
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)
end if
interp_num = interp_num +1
@ -321,14 +321,14 @@ module io
do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
interp_num = interp_num+1
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)
end do
else
do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
interp_num = interp_num+1
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 if
end select
@ -436,6 +436,7 @@ module io
16 format(/'POINT_DATA', i16)
17 format('SCALARS weight float', / &
'LOOKUP_TABLE default')
18 format('SCALARS atom_type float', / &
'LOOKUP_TABLE default')
@ -445,6 +446,9 @@ module io
21 format('SCALARS esize float', /&
'LOOKUP_TABLE default')
22 format('SCALARS energy float', /&
'LOOKUP_TABLE default')
!First we write the vtk file containing the atoms
open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
@ -467,6 +471,14 @@ module io
do i = 1, atom_num
write(11, '(i16)') type_atom(i)
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)
!Now we write the vtk file for the elements
@ -499,6 +511,18 @@ module io
do i = 1, ele_num
write(11, '(i16)') size_ele(i)
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)
end subroutine
@ -509,7 +533,7 @@ module io
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, &
etype
real(kind=dp) :: box_vec(3), masses(10)
real(kind=dp) :: box_vec(3)
1 format('time' / i16, f23.15)
2 format('number of elements' / i16)
@ -548,9 +572,6 @@ module io
write(11,3) node_num
write(11,19) max_ng_node, max_basisnum
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(i), i = 1, atom_types)
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
write(11, 16)
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
else
write(11,14)
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 if
end if
@ -631,7 +652,7 @@ module io
write(11,*) box_bd(:)
!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(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
@ -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), &
node_map(10)
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(:,:,:)
character(len = 2) :: new_type_to_name(10)
!First open the file
open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
@ -777,11 +797,11 @@ module io
call grow_box(temp_box_bd)
!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
!type of the atoms within this file
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
!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)
@ -893,8 +913,7 @@ module io
!Read define atom_types by mass
print *, j
do i = 1, j
call realatomspecies(atomic_masses(i), atomic_element)
call add_atom_type(atomic_element, atom_type_map(i), all_new)
call add_atom_type(atomic_masses(i), atom_type_map(i), all_new)
print *, i, atom_type_map(i)
end do
@ -979,7 +998,7 @@ module io
!Internal Variables
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),&
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
@ -1029,8 +1048,7 @@ module io
!Read define atom_types by mass
do i = 1, j
call realatomspecies(atomic_masses(i), atomic_element)
call add_atom_type(atomic_element, atom_type_map(i), all_new)
call add_atom_type(atomic_masses(i), atom_type_map(i), all_new)
end do
if(in_atoms > 0 ) then
@ -1046,11 +1064,11 @@ module io
do ia = 1, in_atoms
read(11,'(a)') line(:)
if(read_vel) then
read(line,*) tag, type, ra(:), ea, fa(:), va(:), vela(:)
read(line,*) tag, atype, ra(:), ea, fa(:), va(:), vela(:)
else
read(line,*) tag, type, ra(:), ea, fa(:), va(:)
read(line,*) tag, atype, ra(:), ea, fa(:), va(:)
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)
if(read_vel) vel_atom(:,atom_num) = vela
end do
@ -1152,8 +1170,7 @@ module io
!Read atomic masses
do i = 1, type_in
read(11,*) j, mass
call realatomspecies(mass, atom_species)
call add_atom_type(atom_species, type_map(i), all_new)
call add_atom_type(mass, type_map(i), all_new)
end do
!Read useless info
@ -1226,8 +1243,7 @@ module io
!Read atomic masses
do i = 1, type_in
read(11,*) j, mass, textholder
call realatomspecies(mass, atom_species)
call add_atom_type(atom_species, type_map(i), all_new)
call add_atom_type(mass, type_map(i), all_new)
end do
!Read useless info
@ -1311,6 +1327,7 @@ module io
integer :: i, j,arglen, arg_pos, ntypes
character(len=100) :: textholder
real(kind=dp) :: mass
arg_pos = apos + 1
call get_command_argument(arg_pos, textholder, arglen)
@ -1320,7 +1337,8 @@ module io
do i=1,ntypes
arg_pos = arg_pos + 1
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
return

@ -25,6 +25,9 @@ module mode_calc
case('tot_virial')
allocate(calculated(6))
call calc_tot_virial
case('tot_energy')
allocate(calculated(1))
call calc_tot_energy
case default
print *, trim(adjustl(calc_opt)), " is not accepted as a calc option in mode_calc"
stop 3
@ -87,5 +90,26 @@ module mode_calc
print *, "Total virial is calculated as : (v11, v22, v33, v32, v31, v21)"
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

@ -12,9 +12,9 @@ module mode_create
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), &
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), &
basis_pos(3,10), esize_nums, esize_index(10)
esize_nums, esize_index(10)
logical :: dup_flag, dim_flag, efill, crossb(3)
real(kind=dp), allocatable :: r_lat(:,:,:), r_atom_lat(:,:)
@ -186,6 +186,7 @@ module mode_create
character(len=100) :: orient_string
character(len=2) :: btype
logical :: isortho, isrighthanded, bool
real(kind=dp) :: mass
!Pull out all required positional arguments
@ -265,7 +266,8 @@ module mode_create
do i = 1, max_basisnum
call get_command_argument(arg_pos, btype)
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
call get_command_argument(arg_pos, textholder)
read(textholder, *) basis_pos(j,i)
@ -334,7 +336,8 @@ module mode_create
if (basisnum(1) == 0) then
max_basisnum = 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
end if
@ -371,6 +374,8 @@ module mode_create
numlatpoints = numlatpoints*(bd_in_lat(2*i)-bd_in_lat(2*i-1))
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
select case(esize)
!Atomistics

@ -8,7 +8,9 @@ module mode_da
implicit none
integer :: num_d_points
integer, allocatable :: is_slipped(:,:,:)
real(kind=dp), allocatable :: disl_line(:,:), disl_data(:)
private :: write_xyz
public
@ -31,16 +33,28 @@ module mode_da
if(allocated(is_slipped)) deallocate(is_slipped)
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
!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.)
if ( ppos > 0 ) then
outfile = 'da_'//infiles(1)(1:ppos)//'xyz'
else
outfile = 'da_'//infiles(1)//'.xyz'
end if
call write_da_xyz(outfile)
return
@ -65,16 +79,39 @@ module mode_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)
real(kind = dp) :: face_centroids(3, ele_num*6), r(3), rc, vnode(3, max_basisnum, 4), vnorm(max_basisnum, 4), &
max_node_dist, ndiff, rmax(3), rmin(3), rnorm
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_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
l = 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
do i = 1, ele_num
do j = 1, 6
@ -101,8 +138,6 @@ module mode_da
end do
!Calculate the largest difference between nodes for each face
end do
!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
do j = 1, 3
vnode(j,1:basisnum(lat_ele(face_ele(i))),:) = vnode(j,1:basisnum(lat_ele(face_ele(i))),:) - &
minval(vnode(j,1:basisnum(lat_ele(face_ele(i))),:))
end do
!Now calculate the maximum difference between nodes
vnorm = 0
do j = 1, 4
do ibasis = 1, basisnum(lat_ele(face_ele(i)))
vnorm(ibasis,j) = norm2(vnode(:, ibasis, j))
do k=j,4
v1=vnode(:,1,j) - vnode(:,1,k)
v1norm=norm2(v1)
if (vnorm < v1norm) then
vnorm = v1norm
diffvec= v1
end if
end do
end do
!Now calculate the difference between the largest norm and the smallest norm, if it's larger than 0.5 then mark it
!as slipped. This value probably can be converted to a variable value that depends on the current lattice parameter
!I think 0.5 works ok though.
if (any(vnorm > 0.5_dp)) then
if (vnorm>0.5_dp) then
print *, "Element number ", face_ele(i), " is dislocated along face ", face_types(i), &
" with neighbor ", face_ele(nei), " with max displacement of ", maxval(vnorm)
" with neighbor ", face_ele(nei), " with max displacement of ", vnorm
l=0
do j = 1, 4
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
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 do
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)
!This subroutine write the element positions to a .xyz file and marks whether they are slipped or not
character(len=*), intent(in) :: outfile
@ -184,4 +309,4 @@ module mode_da
close(11)
end subroutine write_da_xyz
end module
end module mode_da

@ -10,7 +10,7 @@ module mode_metric
integer :: nfiles
character(len=100) :: metric_type
real(kind=dp) :: rc_off
real(kind=dp) :: rc_off, old_box_len(3)
!Save reference positions
integer :: np, npreal, nmet
@ -57,6 +57,10 @@ module mode_metric
print *, "Getting initial neighbor list"
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
call reset_data
call reset_box
@ -128,8 +132,11 @@ module mode_metric
!This subroutine calculates the continuum metric that we require
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), &
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
do ip = 1, np
eta(:,:) = 0.0_dp
@ -138,8 +145,17 @@ module mode_metric
do jp = 1, nei_num(ip)
!Calculate the neighbor vec in current configuration
nei = nei_list(jp, ip)
rij = r_curr(:,nei) - r_curr(:,ip)
oldrij = r_zero(:,nei) - r_zero(:,ip)
rij = r_curr(:,nei) - r_curr(:,ip) +periodvec(:,jp,ip)*box_len
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
do i = 1,3

@ -4,10 +4,12 @@ module neighbors
use elements
use subroutines
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(:,:)
public
contains
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
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
integer, dimension(3) ::cell_num
integer, allocatable :: num_in_cell(:,:,:), cell_list(:,:,:,:)
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
deallocate(nei_list,nei_num)
deallocate(nei_list,nei_num, periodvec)
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
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
pointloop: do i = 1, n
!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
!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)
!Loop over all neighboring cells
do ci = -1, 1, 1
do cj = -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( (c(1) + ck > cell_num(1)).or.(c(2) + cj > cell_num(2)).or. &
(c(3) + ci > cell_num(3))) cycle
do num_nei = 1, num_in_cell(c(1) + ck, c(2) + cj, c(3) + ci)
nei = cell_list(num_nei,c(1) + ck, c(2) + cj, c(3) + ci)
do num_nei = 1, num_in_cell(cn(1),cn(2),cn(3))
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
!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
!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_list(nei_num(i), i) = nei
periodvec(:,nei_num(i),i) = period_dir
end if
end if

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

@ -87,7 +87,7 @@ module opt_delete
allocate(which_cell(3,atom_num))
!First pass the atom list and atom num to the algorithm which builds the cell list
call build_cell_list(atom_num, r_atom, rc_off, cell_num, num_in_cell, cell_list, which_cell)
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
delete_num = 0

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

@ -2,6 +2,7 @@ module opt_group
!This module contains all code associated with dislocations
use atoms
use parameters
use elements
use subroutines
@ -11,7 +12,7 @@ module opt_group
implicit none
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=2) :: species_type(10)
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
random_num=0
group_type=0
group_atom_types=0
notsize=0
displace=.false.
delete=.false.
@ -46,6 +48,7 @@ module opt_group
flip = .false.
group_nodes=.false.
num_species = 0
exclude_num = 0
if(allocated(element_index)) deallocate(element_index)
if(allocated(atom_index)) deallocate(atom_index)
@ -86,10 +89,10 @@ module opt_group
call displace_group
end if
if(insert_type > 0) then
print *, "Insert command has been dropped"
stop 1
!print *, "Insert command has been dropped"
!stop 1
call get_group
! call insert_group
call insert_group
end if
if(num_species > 0) then
@ -106,7 +109,7 @@ module opt_group
integer :: i, j, arglen, in_num
character(len=100) :: textholder, type_spec
real(kind=dp) H
real(kind=dp) H, mass
!Parse type and shape command
arg_pos = arg_pos + 1
@ -469,6 +472,24 @@ module opt_group
refine=.true.
case('refinefill')
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')
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
@ -493,7 +514,8 @@ module opt_group
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
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')
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
@ -521,7 +543,8 @@ module opt_group
arg_pos=arg_pos+1
call get_command_argument(arg_pos, textholder, arglen)
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
call get_command_argument(arg_pos, textholder, arglen)
select case(trim(adjustl(textholder)))
@ -616,7 +639,7 @@ module opt_group
select case(trim(adjustl(type)))
case('atoms','all')
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
if (group_atom_num > size(atom_index)) then
allocate(resize_array(size(atom_index) + 1024))
@ -643,15 +666,6 @@ module opt_group
end if
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."
end subroutine get_group
@ -693,8 +707,19 @@ module opt_group
!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
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
if(group_ele_num > 0) then
orig_atom_num = atom_num
@ -704,17 +729,38 @@ module opt_group
do i = 1, group_ele_num
ie = element_index(i)
!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)
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
!here as well to make sure they are in the box
do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3
call apply_periodic(r_interp(:,j))
call add_atom(0,type_interp(j), 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
!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
end subroutine refine_group
@ -751,7 +797,7 @@ module opt_group
ravg = ravg + rfill(:,ibasis, 1)/basisnum(lat_ele(ie))
end do
if( in_group(ravg)) then
if( in_group(0,ravg)) then
nump_ele = nump_ele - 1
end if
end do
@ -1124,86 +1170,102 @@ module opt_group
end subroutine change_group_type
! subroutine insert_group
! !This code inserts atoms into interstitial sites. This only works on atoms within the group due to the limitations with the
! !Coarse-graining methodology which doesn't allow for concentration fluctuations.
! real(kind=dp) interstitial_sites(3,14), rand, rinsert(3)
! integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi
! integer, allocatable :: used_sites(:,:)
!
! !First save all of the displacement vectors from a lattice site to interstitial site
! !The first 6 are the octohedral sites and the next 8 are the tetrahedral sites
! interstitial_sites= reshape( (/ -0.5_dp, 0.0_dp, 0.0_dp, &
! 0.5_dp, 0.0_dp, 0.0_dp, &
! 0.0_dp,-0.5_dp, 0.0_dp, &
! 0.0_dp, 0.5_dp, 0.0_dp, &
! 0.0_dp, 0.0_dp,-0.5_dp, &
! 0.0_dp, 0.0_dp, 0.5_dp, &
! -0.25_dp,-0.25_dp,-0.25_dp, &
! -0.25_dp,-0.25_dp, 0.25_dp, &
! -0.25_dp, 0.25_dp,-0.25_dp, &
! -0.25_dp, 0.25_dp, 0.25_dp, &
! 0.25_dp,-0.25_dp,-0.25_dp, &
! 0.25_dp,-0.25_dp, 0.25_dp, &
! 0.25_dp, 0.25_dp,-0.25_dp, &
! 0.25_dp, 0.25_dp, 0.25_dp /), &
! shape(interstitial_sites))
!
! !First we calculate the number of atoms needed based on the number of atoms in the group and the concentration
! interstitial_sites=interstitial_sites*insert_lattice
!
! add_num = (insert_conc*group_atom_num)/(1-insert_conc)
! allocate(used_sites(2,add_num))
!
! print *, "Inserting ", add_num, " atoms as atom type ", insert_type
!
! !Now set up the random number generator for the desired interstitial type
! select case(insert_site)
! case(1)
! rlo=1
! rhi=6
! case(2)
! rlo=7
! rhi = 14
! case(3)
! rlo=1
! rhi=14
! end select
!
! !Now add the atoms
! i = 1
! addloop:do while ( i < add_num)
! call random_number(rand)
! rindex = int(1+rand*(group_atom_num-1))
! ia=atom_index(rindex)
! call random_number(rand)
! sindex = int(rlo+rand*(rhi-rlo))
! do j = 1, i
! if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop
! end do
! rinsert = r_atom(:,ia) + matmul(sub_box_ori(:,:,sbox_atom(ia)),interstitial_sites(:,sindex))
! sbox = sbox_atom(ia)
! call add_atom(0, insert_type, sbox, rinsert)
! used_sites(1,i) = ia
! used_sites(2,i) = sindex
! i = i + 1
! end do addloop
! end subroutine insert_group
subroutine insert_group
!This code inserts atoms into interstitial sites. This only works on atoms within the group due to the limitations with the
!Coarse-graining methodology which doesn't allow for concentration fluctuations.
real(kind=dp) interstitial_sites(3,14), rand, rinsert(3)
integer :: add_num, i, j, rindex, sindex, ia, rlo, rhi
integer, allocatable :: used_sites(:,:)
!First save all of the displacement vectors from a lattice site to interstitial site
!The first 6 are the octohedral sites and the next 8 are the tetrahedral sites
interstitial_sites= reshape( (/ -0.5_dp, 0.0_dp, 0.0_dp, &
0.5_dp, 0.0_dp, 0.0_dp, &
0.0_dp,-0.5_dp, 0.0_dp, &
0.0_dp, 0.5_dp, 0.0_dp, &
0.0_dp, 0.0_dp,-0.5_dp, &
0.0_dp, 0.0_dp, 0.5_dp, &
-0.25_dp,-0.25_dp,-0.25_dp, &
-0.25_dp,-0.25_dp, 0.25_dp, &
-0.25_dp, 0.25_dp,-0.25_dp, &
-0.25_dp, 0.25_dp, 0.25_dp, &
0.25_dp,-0.25_dp,-0.25_dp, &
0.25_dp,-0.25_dp, 0.25_dp, &
0.25_dp, 0.25_dp,-0.25_dp, &
0.25_dp, 0.25_dp, 0.25_dp /), &
shape(interstitial_sites))
!First we calculate the number of atoms needed based on the number of atoms in the group and the concentration
interstitial_sites=interstitial_sites*insert_lattice
!Now set up the random number generator for the desired interstitial type
select case(insert_site)
case(1)
rlo=1
rhi=6
case(2)
rlo=7
rhi = 14
case(3)
rlo=1
rhi=14
end select
if ((insert_conc > 1).or.(insert_conc < 0)) stop "Insert concentration should be between 0 and 1"
add_num = (insert_conc*group_atom_num)
allocate(used_sites(2,add_num))
print *, "Inserting ", add_num, " atoms as atom type ", insert_type
!Now add the atoms
i = 1
addloop:do while ( i < add_num)
call random_number(rand)
rindex = int(1+rand*(group_atom_num-1))
ia=atom_index(rindex)
call random_number(rand)
sindex = int(rlo+rand*(rhi-rlo))
do j = 1, i
if((ia == used_sites(1,i)).and.(sindex == used_sites(2,i))) cycle addloop
end do
rinsert = r_atom(:,ia) + matmul(box_ori,interstitial_sites(:,sindex))
call add_atom(0, insert_type, rinsert)
used_sites(1,i) = ia
used_sites(2,i) = sindex
i = i + 1
end do addloop
end subroutine insert_group
subroutine alloy_group
!This subroutine randomizes the atom types to reach desired concentrations, this only operates on atoms
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)
old_fractions=s_fractions(1:num_species)
!First we sort the fractions
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
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)
end do
!Now randomize the atom types
@ -1225,8 +1287,9 @@ module opt_group
return
end subroutine alloy_group
function in_group(r)
function in_group(atype, r)
!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) :: rnorm
@ -1271,6 +1334,19 @@ module opt_group
case('all')
in_group = .true.
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
function in_group_ele(esize, elat, rn)
@ -1294,13 +1370,13 @@ module opt_group
r_center(:)= r_center(:) + rn(:,ibasis,inod)/basisnum(elat)
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
exit nodeloop
end if
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
else
node_in_out(inod) = 0
@ -1322,7 +1398,7 @@ module opt_group
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.
return
end if
@ -1331,7 +1407,7 @@ module opt_group
r_center(:) = 0.0_dp
do inod = 1, ng_node(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.
return
end if

@ -17,11 +17,13 @@ module opt_slip_plane
!Main calling function for the slip_plane option
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(:)
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)
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(.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)
do ia = 1, basisnum(lat_ele(ie)) * size_ele(ie)**3
call apply_periodic(r_interp(:,ia))
call add_atom(0, type_interp(ia), r_interp(:,ia))
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
else
!First populate the lat points array

Loading…
Cancel
Save