Merge branch 'ft--opt-dislloop' into development
This commit is contained in:
Alex Selimov
commit
c89b7d3b59
11 changed files with 788 additions and 35 deletions
23
src/Makefile
23
src/Makefile
|
|
@ -1,9 +1,9 @@
|
|||
|
||||
FC=ifort
|
||||
FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
|
||||
#FFLAGS=-mcmodel=large -Ofast
|
||||
#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin
|
||||
MODES=mode_create.o mode_merge.o mode_convert.o
|
||||
OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES)
|
||||
OPTIONS=opt_disl.o
|
||||
OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o call_option.o $(MODES) $(OPTIONS)
|
||||
|
||||
.SUFFIXES:
|
||||
.SUFFIXES: .c .f .f90 .F90 .o
|
||||
|
|
@ -19,16 +19,25 @@ clean:
|
|||
$(RM) cacmb *.o
|
||||
|
||||
testfuncs: testfuncs.o functions.o subroutines.o
|
||||
$(FC) testfuncs.o functions.o subroutines.o elements.o -o $@
|
||||
$(FC) testfuncs.o functions.o subroutines.o box.o elements.o -o $@
|
||||
|
||||
.PHONY: cleantest
|
||||
cleantest:
|
||||
$(RM) testfuncs testfuncs.o
|
||||
|
||||
.PHONY: test
|
||||
test: testfuncs
|
||||
./testfuncs
|
||||
|
||||
.PHONY: install
|
||||
install: cacmb
|
||||
cp ./cacmb /usr/local/bin
|
||||
|
||||
$(OBJECTS) : parameters.o
|
||||
atoms.o subroutines.o testfuncs.o : functions.o
|
||||
main.o io.o build_subroutines.o: elements.o
|
||||
atoms.o subroutines.o testfuncs.o box.o : functions.o
|
||||
main.o io.o $(MODES) $(OPTIONS) : elements.o
|
||||
call_mode.o : $(MODES)
|
||||
call_option.o : $(OPTIONS)
|
||||
$(MODES) io.o: atoms.o box.o
|
||||
$(MODES) main.o : io.o
|
||||
testfuncs.o elements.o mode_create.o: subroutines.o
|
||||
testfuncs.o elements.o mode_create.o opt_disl.o: subroutines.o
|
||||
|
|
|
|||
29
src/box.f90
29
src/box.f90
|
|
@ -1,7 +1,7 @@
|
|||
module box
|
||||
!This module contains information on the properties of the current box.
|
||||
use parameters
|
||||
|
||||
use functions
|
||||
implicit none
|
||||
|
||||
real(kind=dp) :: box_bd(6) !Global box boundaries
|
||||
|
|
@ -34,8 +34,14 @@ module box
|
|||
|
||||
integer, intent(in) :: n
|
||||
|
||||
allocate(sub_box_ori(3,3,n), sub_box_bd(6,n), sub_box_array_bd(2,2,n))
|
||||
integer :: i
|
||||
|
||||
allocate(sub_box_ori(3,3,n), sub_box_bd(6,n), sub_box_array_bd(2,2,n))
|
||||
do i = 1, n
|
||||
sub_box_ori(:,:,i) = identity_mat(3)
|
||||
sub_box_bd(:,i) = 0.0_dp
|
||||
sub_box_array_bd(:,:,i) = 1
|
||||
end do
|
||||
end subroutine alloc_sub_box
|
||||
|
||||
subroutine grow_sub_box(n)
|
||||
|
|
@ -58,7 +64,7 @@ module box
|
|||
call move_alloc(temp_bd, sub_box_bd)
|
||||
|
||||
temp_array_bd(:,:,1:sub_box_num) = sub_box_array_bd
|
||||
temp_array_bd(:,:,sub_box_num+1:) = 0.0_dp
|
||||
temp_array_bd(:,:,sub_box_num+1:) = 1
|
||||
call move_alloc(temp_array_bd, sub_box_array_bd)
|
||||
|
||||
return
|
||||
|
|
@ -79,4 +85,21 @@ module box
|
|||
return
|
||||
end subroutine grow_box
|
||||
|
||||
|
||||
subroutine in_sub_box(r, which_sub_box)
|
||||
!This returns which sub_box a point is in. It returns the first sub_box with boundaries which
|
||||
!contain the point.
|
||||
real(kind=dp), dimension(3), intent(in) :: r
|
||||
integer, intent(out) :: which_sub_box
|
||||
|
||||
integer :: i
|
||||
|
||||
do i = 1, sub_box_num
|
||||
if( in_block_bd(r, sub_box_bd(:,i))) then
|
||||
which_sub_box = i
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
return
|
||||
end subroutine in_sub_box
|
||||
end module box
|
||||
|
|
@ -3,6 +3,7 @@ module elements
|
|||
!This module contains the elements data structures, structures needed for building regions
|
||||
!and operations that are done on elements
|
||||
use parameters
|
||||
use functions
|
||||
use subroutines
|
||||
implicit none
|
||||
|
||||
|
|
@ -373,4 +374,26 @@ module elements
|
|||
return
|
||||
end subroutine rhombshape
|
||||
|
||||
subroutine delete_atoms(num, index)
|
||||
!This subroutine deletes atoms from the atom arrays
|
||||
integer, intent(in) :: num
|
||||
integer, intent(inout), dimension(num) :: index
|
||||
|
||||
integer :: i, j
|
||||
|
||||
call heapsort(index)
|
||||
|
||||
!We go from largest index to smallest index just to make sure that we don't miss
|
||||
!accidentally overwrite values which need to be deleted
|
||||
do i = num, 1, -1
|
||||
if(index(i) == atom_num) then
|
||||
r_atom(:,index(i)) = 0.0_dp
|
||||
type_atom(index(i)) = 0
|
||||
else
|
||||
r_atom(:,index(i)) = r_atom(:, atom_num)
|
||||
type_atom(index(i)) = type_atom(atom_num)
|
||||
end if
|
||||
atom_num = atom_num - 1
|
||||
end do
|
||||
end subroutine
|
||||
end module elements
|
||||
|
|
@ -232,4 +232,22 @@ END FUNCTION StrDnCase
|
|||
end if
|
||||
return
|
||||
end function is_equal
|
||||
|
||||
pure function unitvec(n,vec)
|
||||
integer, intent(in) :: n
|
||||
real(kind=dp), intent(in) :: vec(n)
|
||||
real(kind=dp) :: unitvec(n)
|
||||
|
||||
unitvec = vec/norm2(vec)
|
||||
return
|
||||
end function unitvec
|
||||
|
||||
pure function norm_dis(rl,rk)
|
||||
!This just returns the magnitude of the vector between two points
|
||||
real(kind=dp), dimension(3), intent(in) :: rl, rk
|
||||
real(kind=dp) :: norm_dis(4)
|
||||
|
||||
norm_dis(1:3) = (rk - rl)
|
||||
norm_dis(4) = norm2(rk-rl)
|
||||
end function
|
||||
end module functions
|
||||
|
|
|
|||
14
src/io.f90
14
src/io.f90
|
|
@ -387,7 +387,7 @@ module io
|
|||
write(11,3) node_num
|
||||
write(11,19) max_ng_node
|
||||
write(11,4) lattice_types
|
||||
write(11,2) atom_num
|
||||
write(11,5) atom_num
|
||||
write(11,6) 1 !Grain_num is ignored
|
||||
write(11,16) lattice_types, atom_types
|
||||
write(11,21)
|
||||
|
|
@ -481,7 +481,7 @@ module io
|
|||
end do
|
||||
|
||||
!Write the number of atom types in the current model and all of their names
|
||||
write(11,*) atom_types, (type_to_name(i), i=1, atom_types)
|
||||
write(11,*) atom_types, (type_to_name(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
|
||||
|
|
@ -616,10 +616,14 @@ module io
|
|||
read(11,*) sub_box_bd(:,sub_box_num+i)
|
||||
sub_box_bd(:,sub_box_num+i) = sub_box_bd(:, sub_box_num+i) + displace(:)
|
||||
!Read in sub_box_array_bd
|
||||
read(11,*) ((sub_box_ori(j, k, sub_box_num+i), j = 1, 2), k = 1, 2)
|
||||
|
||||
read(11,*) ((sub_box_array_bd(j, k, sub_box_num+i), j = 1, 2), k = 1, 2)
|
||||
end do
|
||||
sub_box_num = sub_box_num + n
|
||||
|
||||
!Add the existing element boundaries
|
||||
sub_box_array_bd(:,1,sub_box_num+1:) = sub_box_array_bd(:,1,sub_box_num+1:) + atom_num
|
||||
sub_box_array_bd(:,2,sub_box_num+1:) = sub_box_array_bd(:,2,sub_box_num+1:) + ele_num
|
||||
|
||||
sub_box_num = sub_box_num + n
|
||||
|
||||
!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)
|
||||
|
|
|
|||
24
src/main.f90
24
src/main.f90
|
|
@ -17,12 +17,13 @@ program main
|
|||
use io
|
||||
|
||||
|
||||
integer :: i, end_mode_arg, arg_num
|
||||
integer :: i, end_mode_arg, arg_num, arg_pos
|
||||
character(len=100) :: argument
|
||||
|
||||
!Call initialization functions
|
||||
call lattice_init
|
||||
call box_init
|
||||
end_mode_arg = 0
|
||||
|
||||
! Command line parsing
|
||||
arg_num = command_argument_count()
|
||||
|
|
@ -47,13 +48,30 @@ program main
|
|||
call call_mode(end_mode_arg, argument)
|
||||
end if
|
||||
|
||||
!Check to make sure a mode was called
|
||||
if (end_mode_arg==0) then
|
||||
stop "Nothing to do, please run cacmb using an available mode"
|
||||
end if
|
||||
!Now we loop through all of the arguments and check for passed options or for a filename to be written out
|
||||
do i = end_mode_arg-1, arg_num
|
||||
call get_command_argument(i, argument)
|
||||
arg_pos = end_mode_arg
|
||||
do while(.true.)
|
||||
!Exit the loop if we are done reading
|
||||
if(arg_pos > arg_num) exit
|
||||
|
||||
call get_command_argument(arg_pos, argument)
|
||||
|
||||
!Check to see if a filename was passed
|
||||
if(scan(argument,'.',.true.) > 0) then
|
||||
call get_out_file(argument)
|
||||
arg_pos = arg_pos + 1
|
||||
|
||||
!Check to see if an option has been passed
|
||||
else if(argument(1:1) == '-') then
|
||||
call call_option(argument, arg_pos)
|
||||
!Otherwise print that the argument is not accepted and move on
|
||||
else
|
||||
print *, trim(adjustl(argument)), " is not accepted. Skipping to next argument"
|
||||
arg_pos = arg_pos + 1
|
||||
end if
|
||||
end do
|
||||
|
||||
|
|
|
|||
|
|
@ -199,20 +199,7 @@ module mode_create
|
|||
call get_command_argument(arg_pos, orient_string, arglen)
|
||||
if (arglen==0) STOP "Missing orientation in orient command of mode create"
|
||||
arg_pos = arg_pos+1
|
||||
ori_pos=2
|
||||
do j = 1,3
|
||||
if (orient_string(ori_pos:ori_pos) == '-') then
|
||||
ori_pos = ori_pos + 1
|
||||
read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
|
||||
if (stat>0) STOP "Error reading orient value"
|
||||
orient(i,j) = -orient(i,j)
|
||||
ori_pos = ori_pos + 1
|
||||
else
|
||||
read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
|
||||
if(stat>0) STOP "Error reading orient value"
|
||||
ori_pos=ori_pos + 1
|
||||
end if
|
||||
end do
|
||||
call parse_ori_vec(orient_string, orient(i,:))
|
||||
end do
|
||||
|
||||
|
||||
|
|
@ -243,10 +230,10 @@ module mode_create
|
|||
|
||||
case default
|
||||
!If it isn't an option then you have to exit
|
||||
arg_pos = arg_pos -1
|
||||
exit
|
||||
end select
|
||||
end do
|
||||
|
||||
!Calculate the lattice periodicity length in lattice units
|
||||
do i = 1, 3
|
||||
lattice_space(i) = norm2(orient(i,:))
|
||||
|
|
|
|||
520
src/opt_disl.f90
Normal file
520
src/opt_disl.f90
Normal file
|
|
@ -0,0 +1,520 @@
|
|||
module opt_disl
|
||||
|
||||
!This module contains all code associated with dislocations
|
||||
|
||||
use parameters
|
||||
use elements
|
||||
use subroutines
|
||||
use box
|
||||
implicit none
|
||||
|
||||
|
||||
real(kind=dp), dimension(3) :: line, slip_plane, centroid!dislocation line, slip plane vectors, centroid,
|
||||
real(kind=dp) :: burgers(3) !burgers vector of loop
|
||||
real(kind=dp) :: poisson, char_angle, lattice_parameter!Poisson ratio and character angle, lattice_parameter for burgers vector
|
||||
character(len=10) :: lattice
|
||||
character(len=1) :: loop_normal
|
||||
real(kind=dp) :: loop_radius !Dislocation loop radius
|
||||
|
||||
real(kind=dp) :: b !burgers vector
|
||||
|
||||
|
||||
public
|
||||
contains
|
||||
|
||||
subroutine dislocation(option, arg_pos)
|
||||
|
||||
!Main calling function for all codes related to dislocations
|
||||
|
||||
character(len=100), intent(in) :: option
|
||||
integer, intent(inout) :: arg_pos
|
||||
|
||||
|
||||
select case(trim(adjustl(option)))
|
||||
case('-dislgen')
|
||||
call parse_dislgen(arg_pos)
|
||||
call dislgen
|
||||
case('-disloop')
|
||||
call parse_disloop(arg_pos)
|
||||
call disloop
|
||||
end select
|
||||
end subroutine dislocation
|
||||
|
||||
subroutine parse_dislgen(arg_pos)
|
||||
|
||||
!Parse dislgen command
|
||||
|
||||
integer, intent(inout) :: arg_pos
|
||||
|
||||
integer :: i,arglen
|
||||
character(len=8) :: ori_string
|
||||
character(len=100) :: textholder
|
||||
|
||||
!Parse all of the commands
|
||||
arg_pos = arg_pos + 1
|
||||
line(:) = 0.0_dp
|
||||
call get_command_argument(arg_pos, ori_string, arglen)
|
||||
if (arglen==0) STOP "Missing line vector in dislgen command"
|
||||
call parse_ori_vec(ori_string, line)
|
||||
|
||||
arg_pos=arg_pos + 1
|
||||
slip_plane(:) = 0.0_dp
|
||||
call get_command_argument(arg_pos, ori_string, arglen)
|
||||
if (arglen==0) STOP "Missing plane vector in dislgen command"
|
||||
call parse_ori_vec(ori_string, slip_plane)
|
||||
|
||||
do i = 1, 3
|
||||
arg_pos = arg_pos + 1
|
||||
call get_command_argument(arg_pos, textholder, arglen)
|
||||
if (arglen==0) STOP "Missing centroid in dislgen command"
|
||||
call parse_pos(i, textholder, centroid(i))
|
||||
end do
|
||||
|
||||
print *, centroid
|
||||
|
||||
arg_pos = arg_pos + 1
|
||||
call get_command_argument(arg_pos, textholder, arglen)
|
||||
if (arglen==0) STOP "Missing character angle in dislgen command"
|
||||
read(textholder, *) char_angle
|
||||
|
||||
|
||||
arg_pos = arg_pos + 1
|
||||
call get_command_argument(arg_pos, textholder, arglen)
|
||||
if (arglen==0) STOP "Missing poisson in dislgen command"
|
||||
read(textholder, *) poisson
|
||||
|
||||
arg_pos = arg_pos + 1
|
||||
call get_command_argument(arg_pos, lattice, arglen)
|
||||
if (arglen==0) STOP "Missing lattice in dislgen command"
|
||||
|
||||
arg_pos = arg_pos + 1
|
||||
call get_command_argument(arg_pos, textholder, arglen)
|
||||
if (arglen==0) STOP "Missing lattice parameter in dislgen command"
|
||||
read(textholder, *) lattice_parameter
|
||||
|
||||
arg_pos = arg_pos + 1
|
||||
!Now set the vurgers vector based on the lattice paarameter and lattice type
|
||||
select case(lattice)
|
||||
case('fcc')
|
||||
b = lattice_parameter / sqrt(2.0_dp)
|
||||
! case('bcc')
|
||||
! b = lattice_parameter * sqrt(3.0_dp) / 2.0_dp
|
||||
case default
|
||||
print *, 'Error: Lattice structure', lattice, ' is not accepted for dislgen option'
|
||||
STOP
|
||||
end select
|
||||
|
||||
end subroutine parse_dislgen
|
||||
|
||||
subroutine dislgen
|
||||
!This subroutine creates the actual dislocation
|
||||
|
||||
integer :: i, sub_box, inod, ibasis
|
||||
real(kind=dp) :: ss_ori(3,3), ss_inv(3,3), be, bs, slipx(3), disp_transform(3,3), inv_transform(3,3), &
|
||||
actan, r(3), disp(3)
|
||||
|
||||
!Calculate screw and edge burgers vectors
|
||||
be = sin(char_angle*pi/180.0_dp)*b
|
||||
bs = cos(char_angle*pi/180.0_dp)*b
|
||||
|
||||
!Figure out which sub box you are in so you can access it's orientation, this code will not work
|
||||
!with overlapping sub_boxes
|
||||
do i = 1, sub_box_num
|
||||
if(in_block_bd(centroid, sub_box_bd(:,i))) then
|
||||
sub_box = i
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
!Construct the slip system orientation matrix in an unrotated system
|
||||
slipx = cross_product(slip_plane, line)
|
||||
ss_ori(1,:) = unitvec(3,slipx)
|
||||
ss_ori(2,:) = unitvec(3,slip_plane)
|
||||
ss_ori(3,:) = unitvec(3,line)
|
||||
call matrix_inverse(ss_ori, 3, ss_inv)
|
||||
|
||||
!Apply the rotation
|
||||
disp_transform = matmul(sub_box_ori(:,:,i), ss_inv)
|
||||
call matrix_inverse(disp_transform,3,inv_transform)
|
||||
|
||||
if(atom_num > 0) then
|
||||
do i = 1, atom_num
|
||||
r=r_atom(:,i) - centroid
|
||||
r=matmul(inv_transform, r)
|
||||
if (r(1) == 0) then
|
||||
actan=pi/2
|
||||
else
|
||||
actan = datan2(r(2),r(1))
|
||||
end if
|
||||
|
||||
if ((r(1)**2 + r(2)**2) == 0) cycle
|
||||
|
||||
!This is the elastic displacement field for dislocation according to Hirth and Lowe
|
||||
disp(1) = be/(2.0_dp*pi) * (actan + (r(1)*r(2))/(2.0_dp*(1.0_dp-poisson)*(r(1)**2.0_dp + r(2)**2.0_dp)))
|
||||
disp(2) = -be/(2.0_dp*pi)*((1.0_dp-2.0_dp*poisson)/(4.0_dp-4.0_dp*poisson) * &
|
||||
log(r(1)**2.0_dp + r(2)**2.0_dp) &
|
||||
+ (r(1)**2.0_dp - r(2)**2.0_dp)/(4.0_dp*(1.0_dp-poisson)&
|
||||
*(r(1)**2.0_dp+r(2)**2.0_dp)))
|
||||
disp(3) = bs/(2.0_dp*pi) * actan
|
||||
|
||||
!This transforms the displacement to the correct orientation
|
||||
disp = matmul(disp_transform, disp)
|
||||
|
||||
r_atom(:,i) = r_atom(:,i) + disp
|
||||
end do
|
||||
end if
|
||||
|
||||
if(ele_num > 0) then
|
||||
do i = 1, ele_num
|
||||
do inod=1, ng_node(lat_ele(i))
|
||||
do ibasis = 1, basisnum(lat_ele(i))
|
||||
r = r_node(:,ibasis,inod,i)
|
||||
r = matmul(inv_transform, r)
|
||||
if (r(1) == 0) then
|
||||
actan = pi/2
|
||||
else
|
||||
actan = datan2(r(2),r(1))
|
||||
end if
|
||||
|
||||
if ((r(1)**2 + r(2)**2) == 0) cycle
|
||||
|
||||
!This is the elastic displacement field for dislocation according to Hirth and Lowe
|
||||
disp(1) = be/(2.0_dp*pi) * (actan + (r(1)*r(2))/(2.0_dp*(1.0_dp-poisson)*(r(1)**2.0_dp + r(2)**2.0_dp)))
|
||||
disp(2) = -be/(2.0_dp*pi)*((1.0_dp-2.0_dp*poisson)/(4.0_dp-4.0_dp*poisson) * &
|
||||
log(r(1)**2.0_dp + r(2)**2.0_dp) &
|
||||
+ (r(1)**2.0_dp - r(2)**2.0_dp)/(4.0_dp*(1.0_dp-poisson)&
|
||||
*(r(1)**2.0_dp+r(2)**2.0_dp)))
|
||||
disp(3) = bs/(2.0_dp*pi) * actan
|
||||
|
||||
|
||||
disp = matmul(disp_transform, disp)
|
||||
r_node(:,ibasis,inod,i) = r_node(:,ibasis,inod,i) + disp
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
|
||||
end subroutine dislgen
|
||||
|
||||
subroutine parse_disloop(arg_pos)
|
||||
!This subroutine parses the disloop command
|
||||
|
||||
integer, intent(inout) :: arg_pos
|
||||
|
||||
integer :: i,arglen, sbox
|
||||
character(len=8) :: ori_string
|
||||
character(len=100) :: textholder
|
||||
|
||||
!Parse all of the commands
|
||||
arg_pos = arg_pos + 1
|
||||
loop_normal = ' '
|
||||
call get_command_argument(arg_pos, loop_normal, arglen)
|
||||
if (arglen==0) STOP "Missing loop_normal in disloop command"
|
||||
!Convert the loop_normal to the dimension
|
||||
select case(loop_normal)
|
||||
case('x','X', 'y', 'Y', 'z', 'Z')
|
||||
continue
|
||||
case default
|
||||
print *, "Dimension argument must either be x, y, or z not", loop_normal
|
||||
stop 3
|
||||
end select
|
||||
|
||||
arg_pos = arg_pos + 1
|
||||
loop_radius = 0
|
||||
call get_command_argument(arg_pos, textholder, arglen)
|
||||
if (arglen==0) STOP "Missing loop_size in disloop command"
|
||||
read(textholder, *) loop_radius
|
||||
|
||||
do i = 1, 3
|
||||
arg_pos = arg_pos + 1
|
||||
call get_command_argument(arg_pos, textholder, arglen)
|
||||
if (arglen==0) STOP "Missing centroid in disloop command"
|
||||
call parse_pos(i, textholder, centroid(i))
|
||||
end do
|
||||
|
||||
burgers(:) = 0.0_dp
|
||||
do i = 1, 3
|
||||
arg_pos = arg_pos + 1
|
||||
call get_command_argument(arg_pos, textholder, arglen)
|
||||
if (arglen==0) STOP "Missing burgers vector in disloop command"
|
||||
read(textholder, *) burgers(i)
|
||||
end do
|
||||
|
||||
arg_pos = arg_pos + 1
|
||||
call get_command_argument(arg_pos, textholder, arglen)
|
||||
if (arglen==0) STOP "Missing poisson ratio in disloop command"
|
||||
read(textholder, *) poisson
|
||||
|
||||
arg_pos = arg_pos + 1
|
||||
|
||||
!Now check to make sure that the dimension selected is actually a 1 1 1 direction.
|
||||
! call in_sub_box(centroid, sbox)
|
||||
! if(.not.((abs(sub_box_ori(loop_normal,1,sbox)) == abs(sub_box_ori(loop_normal,2,sbox))).and. &
|
||||
! (abs(sub_box_ori(loop_normal,2,sbox)) == abs(sub_box_ori(loop_normal,3,sbox))).and. &
|
||||
! (abs(sub_box_ori(loop_normal,3,sbox)) == abs(sub_box_ori(loop_normal,1,sbox))))) then
|
||||
|
||||
! print *, "The selected dimension ", loop_normal, " for sub_box ", sbox, " is ", &
|
||||
! sub_box_ori(loop_normal,:,sbox), " which is not in the (111) family of planes"
|
||||
! STOP 3
|
||||
! end if
|
||||
|
||||
|
||||
end subroutine
|
||||
|
||||
!Code for the creation of dislocation loops is based on functions from atomsk
|
||||
!which is available from https://atomsk.univ-lille.fr/. They have been adapted to fit cacmb.
|
||||
subroutine disloop
|
||||
!This subroutine applies the displacement field for a dislocation loop to all atoms and elements.
|
||||
integer :: a1, a2, a3 !New directions
|
||||
integer :: i, j, inod, ibasis, Npoints
|
||||
real(kind = dp) :: perimeter, angle, theta, omega, xA(3), xB(3), xC(3), u(3)
|
||||
real(kind=dp), dimension(:,:), allocatable :: xloop !coordinate of points forming loop
|
||||
|
||||
if(allocated(xLoop)) deallocate(xLoop)
|
||||
|
||||
!Define new directions
|
||||
select case(loop_normal)
|
||||
case('x','X')
|
||||
!Loop in (y,z) plane
|
||||
a1 = 2
|
||||
a2 = 3
|
||||
a3 = 1
|
||||
case('y','Y')
|
||||
!Loop in (x,z) plane
|
||||
a1 = 3
|
||||
a2 = 1
|
||||
a3 = 2
|
||||
case default
|
||||
!Loop in (x,y) plane
|
||||
a1 = 1
|
||||
a2 = 2
|
||||
a3 = 3
|
||||
end select
|
||||
|
||||
!Calculate loop perimeter
|
||||
perimeter = 2.0_dp*pi*loop_radius
|
||||
|
||||
!Define the number of points forming the loop
|
||||
! The following criteria are used as a trade-off between
|
||||
! good accuracy and computational efficiency:
|
||||
! - each dislocation segment should have a length of 5 angströms;
|
||||
! - the loop should contain at least 3 points
|
||||
! (for very small loops, this will result in segments shorter than 5 A);
|
||||
! - there should not be more than 100 points
|
||||
! (for very large loops, this will result in segments longer than 5 A).
|
||||
Npoints = MAX( 3 , MIN( NINT(perimeter/5.d0) , 100 ) )
|
||||
|
||||
!angle between two consecutive points
|
||||
theta = 2.0_dp*pi / dble(Npoints)
|
||||
|
||||
!allocate xLoop
|
||||
allocate(xLoop(Npoints,3))
|
||||
xLoop(:,:) = 0.0_dp
|
||||
|
||||
!Calculate the position of each point in the loop
|
||||
angle = 0.0_dp
|
||||
do i = 1, size(xLoop,1)
|
||||
xLoop(i,a1) = centroid(a1) + loop_radius*dcos(angle)
|
||||
xLoop(i,a2) = centroid(a2) + loop_radius*dsin(angle)
|
||||
xLoop(i,a3) = centroid(a3)
|
||||
! Increment angle for next point
|
||||
angle = angle + theta
|
||||
end do
|
||||
|
||||
!Now actually calculate the displacement created by a loop for every atom
|
||||
do i = 1, atom_num
|
||||
u = 0.0_dp
|
||||
omega = 0.0_dp
|
||||
xC(:) = centroid - r_atom(:,i)
|
||||
|
||||
!Loop over all dislocation segments
|
||||
do j = 1, size(xLoop,1)
|
||||
!Coordinates of point A
|
||||
if (j ==1) then
|
||||
xA(:) = xLoop(size(xLoop,1),:) - r_atom(:,i)
|
||||
else
|
||||
xA(:) = xLoop(j-1, :) - r_atom(:,i)
|
||||
end if
|
||||
|
||||
!Coordinates of point B
|
||||
xB(:) = xLoop(j,:) - r_atom(:,i)
|
||||
|
||||
!Displacement due to solid angle
|
||||
omega = omega + SolidAngle(xA, xB, xC)
|
||||
!Displacement due to elasticity
|
||||
u(:) = u(:) + DisloSeg_displacement_iso(xA, xB, burgers(:), poisson)
|
||||
end do
|
||||
|
||||
!Total displacement
|
||||
u=u(:) + burgers(:) * omega
|
||||
|
||||
!Apply displacement to atom
|
||||
r_atom(:,i) = r_atom(:,i) + u(:)
|
||||
end do
|
||||
|
||||
!Repeat for element nodes
|
||||
do i = 1, ele_num
|
||||
do inod =1, ng_node(lat_ele(i))
|
||||
do ibasis = 1, basisnum(lat_ele(i))
|
||||
u = 0.0_dp
|
||||
omega = 0.0_dp
|
||||
xC(:) = centroid - r_node(:,ibasis,inod,i)
|
||||
|
||||
!Loop over all dislocation segments
|
||||
do j = 1, size(xLoop,1)
|
||||
!Coordinates of point A
|
||||
if (j ==1) then
|
||||
xA(:) = xLoop(size(xLoop,1),:) - r_node(:,ibasis,inod,i)
|
||||
else
|
||||
xA(:) = xLoop(j-1, :) - r_node(:,ibasis,inod,i)
|
||||
end if
|
||||
|
||||
!Coordinates of point B
|
||||
xB(:) = xLoop(j,:) - r_node(:,ibasis,inod,i)
|
||||
|
||||
!Displacement due to solid angle
|
||||
omega = omega + SolidAngle(xA, xB, xC)
|
||||
!Displacement due to elasticity
|
||||
u(:) = u(:) + DisloSeg_displacement_iso(xA, xB, burgers(:), poisson)
|
||||
end do
|
||||
|
||||
!Total displacement
|
||||
u=u(:) + burgers(:) * omega
|
||||
|
||||
!Apply displacement to atom
|
||||
r_node(:,ibasis,inod,i) = r_node(:,ibasis,inod,i) + u(:)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
return
|
||||
end subroutine
|
||||
|
||||
!********************************************************
|
||||
! SOLIDANGLE
|
||||
! Calculate solid angle (normalized by 4*pi) used to obtain the
|
||||
! displacement created by a dislocation triangle loop ABC at the origin
|
||||
! (field point = origin)
|
||||
! Ref.: Van Oosterom, A. and Strackee, J.,
|
||||
! The Solid Angle of a Plane Triangle,
|
||||
! IEEE Transactions on Biomedical Engineering BME-30, 125 (1983).
|
||||
!********************************************************
|
||||
FUNCTION SolidAngle(xA, xB, xC) RESULT(Omega)
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
! Extremities of the triangle loop
|
||||
REAL(dp),DIMENSION(3),INTENT(IN) :: xA, xB, xC
|
||||
!
|
||||
! Solid angle (normalized by 4*pi)
|
||||
REAL(dp):: omega
|
||||
REAL(dp),PARAMETER:: factor=1.d0/(2.d0*pi)
|
||||
!
|
||||
REAL(dp) :: rA, rB, rC, numerator, denominator
|
||||
!
|
||||
rA = norm2(xA)
|
||||
rB = norm2(xB)
|
||||
rC = norm2(xC)
|
||||
!
|
||||
numerator = TRIPLE_PRODUCT( xA, xB, xC )
|
||||
denominator = rA*rB*rC + DOT_PRODUCT( xA, xB )*rC &
|
||||
& + DOT_PRODUCT( xB, xC )*rA + DOT_PRODUCT( xC, xA )*rB
|
||||
!
|
||||
omega = factor*ATAN2( numerator, denominator )
|
||||
!
|
||||
END FUNCTION SolidAngle
|
||||
|
||||
!********************************************************
|
||||
! DISLOSEG_DISPLACEMENT_ISO
|
||||
! Calculate displacement created by dislocation segment AB
|
||||
! once the solid angle part has been removed
|
||||
! Isotropic elastic calculation with nu Poisson coef.
|
||||
! Ref.: Eq. (1) in Barnett, D. M.
|
||||
! The Displacement Field of a Triangular Dislocation Loop
|
||||
! Philos. Mag. A, 1985, 51, 383-387
|
||||
!********************************************************
|
||||
FUNCTION DisloSeg_displacement_iso(xA, xB, b, nu) RESULT(u)
|
||||
!
|
||||
IMPLICIT NONE
|
||||
REAL(dp),DIMENSION(3),INTENT(IN):: xA, xB ! Extremities of the segment
|
||||
REAL(dp),DIMENSION(3),INTENT(IN):: b ! Burgers vector
|
||||
REAL(dp),INTENT(IN):: nu ! Poisson coefficient
|
||||
REAL(dp),DIMENSION(3):: u ! Displacement
|
||||
REAL(dp):: rA, rB
|
||||
REAL(dp),DIMENSION(1:3):: tAB, nAB
|
||||
!
|
||||
rA = norm2(xA)
|
||||
rB = norm2(xB)
|
||||
!
|
||||
! Tangent vector
|
||||
tAB(:) = xB(:) - xA(:)
|
||||
tAB(:) = tAB(:)/norm2(tAB)
|
||||
!
|
||||
! Normal vector
|
||||
nAB(:) = CROSS_PRODUCT(xA,xB)
|
||||
nAB(:) = nAB(:)/norm2(nAB)
|
||||
!
|
||||
u(:) = ( -(1.d0-2.d0*nu)*CROSS_PRODUCT(b, tAB)* &
|
||||
& LOG( (rB + DOT_PRODUCT(xB,tAB))/(rA + DOT_PRODUCT(xA,tAB)) ) &
|
||||
& + DOT_PRODUCT(b,nAB)*CROSS_PRODUCT(xB/rB-xA/rA,nAB) ) &
|
||||
& /(8.d0*pi*(1.d0-nu))
|
||||
!
|
||||
END FUNCTION DisloSeg_displacement_iso
|
||||
!
|
||||
!This code simply creates a planar vacancy cluster and does not apply the dislocation loop displacement field.
|
||||
! subroutine disloop
|
||||
! !This subroutine actually creates the dislocation loop.
|
||||
|
||||
! real(kind=dp) :: neighbor_dis(loop_size), temp_box(6), dis
|
||||
! integer :: i, j, index(loop_size)
|
||||
|
||||
! neighbor_dis(:) = HUGE(1.0_dp)
|
||||
! index(:) = 0
|
||||
|
||||
! !First find the nearest atom to the centroid
|
||||
! do i = 1, atom_num
|
||||
! if(norm2(r_atom(:,i) - centroid) < neighbor_dis(1)) then
|
||||
! neighbor_dis(1) = norm2(r_atom(:,i) - centroid)
|
||||
! index(1) = i
|
||||
! end if
|
||||
! end do
|
||||
|
||||
! !Define a new box, this box tries to isolate all atoms on the plane of the atom
|
||||
! !closest to the user defined centroid.
|
||||
! temp_box(:) = box_bd(:)
|
||||
! temp_box(2*loop_normal) = r_atom(loop_normal,index(1)) + 10.0_dp**(-2.0_dp)
|
||||
! temp_box(2*loop_normal-1) = r_atom(loop_normal,index(1)) - 10.0_dp**(-2.0_dp)
|
||||
|
||||
! !Now reset the list for the scanning algorithm
|
||||
! index(1) = 0
|
||||
! neighbor_dis(1) = HUGE(1.0_dp)
|
||||
|
||||
! !Now scan over all atoms again and find the closest loop_size number of atoms to the initial atom
|
||||
! !that reside on the same plane.
|
||||
|
||||
! do i = 1, atom_num
|
||||
! !Check to see if it is on the same plane
|
||||
! if (in_block_bd(r_atom(:,i), temp_box)) then
|
||||
! dis = norm2(r_atom(:,i) - centroid)
|
||||
! do j = 1, loop_size
|
||||
! !Check to see if it is closer than other atoms
|
||||
! if (dis < neighbor_dis(j)) then
|
||||
! !Move values in the neighbor array and add the new neighbor
|
||||
! if(j < loop_size) then
|
||||
! neighbor_dis(j+1:loop_size) = neighbor_dis(j:loop_size-1)
|
||||
! index(j+1:loop_size) = index(j:loop_size-1)
|
||||
! end if
|
||||
! neighbor_dis(j) = dis
|
||||
! index(j) = i
|
||||
! exit
|
||||
! end if
|
||||
! end do
|
||||
! end if
|
||||
! end do
|
||||
|
||||
! !Now delete the atoms
|
||||
! call delete_atoms(loop_size, index)
|
||||
|
||||
! return
|
||||
! end subroutine disloop
|
||||
|
||||
end module opt_disl
|
||||
|
|
@ -2,9 +2,14 @@ module parameters
|
|||
|
||||
implicit none
|
||||
|
||||
!Default precision
|
||||
integer, parameter :: dp= selected_real_kind(15,307)
|
||||
!Parameters for floating point tolerance
|
||||
real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), &
|
||||
lim_large = huge(1.0_dp)
|
||||
logical, save :: lmpcac
|
||||
|
||||
end module parameters
|
||||
!Numeric constants
|
||||
real(kind=dp), parameter :: pi = 3.14159265358979323846_dp
|
||||
|
||||
end module parameters
|
||||
|
|
|
|||
|
|
@ -1,6 +1,7 @@
|
|||
module subroutines
|
||||
use parameters
|
||||
use functions
|
||||
use box
|
||||
implicit none
|
||||
|
||||
integer :: allostat, deallostat
|
||||
|
|
@ -145,4 +146,107 @@ module subroutines
|
|||
return
|
||||
end subroutine matrix_inverse
|
||||
|
||||
subroutine parse_ori_vec(ori_string, ori_vec)
|
||||
!This subroutine parses a string to vector in the format [ijk]
|
||||
character(len=8), intent(in) :: ori_string
|
||||
real(kind=dp), dimension(3), intent(out) :: ori_vec
|
||||
|
||||
integer :: i, ori_pos, stat
|
||||
|
||||
ori_pos=2
|
||||
do i = 1,3
|
||||
if (ori_string(ori_pos:ori_pos) == '-') then
|
||||
ori_pos = ori_pos + 1
|
||||
read(ori_string(ori_pos:ori_pos), *, iostat=stat) ori_vec(i)
|
||||
if (stat>0) STOP "Error reading orientation value"
|
||||
ori_vec(i) = -ori_vec(i)
|
||||
ori_pos = ori_pos + 1
|
||||
else
|
||||
read(ori_string(ori_pos:ori_pos), *, iostat=stat) ori_vec(i)
|
||||
if(stat>0) STOP "Error reading orientation value"
|
||||
ori_pos=ori_pos + 1
|
||||
end if
|
||||
end do
|
||||
|
||||
return
|
||||
end subroutine parse_ori_vec
|
||||
|
||||
subroutine parse_pos(i, pos_string, pos)
|
||||
!This subroutine parses the pos command allowing for command which include inf
|
||||
integer, intent(in) :: i !The dimension of the position
|
||||
character(len=100), intent(in) :: pos_string !The position string
|
||||
real(kind=dp), intent(out) :: pos !The output parsed position value
|
||||
|
||||
integer :: iospara
|
||||
if(trim(adjustl(pos_string)) == 'inf') then
|
||||
pos=box_bd(2*i)
|
||||
else if(trim(adjustl(pos_string)) == '-inf') then
|
||||
pos=box_bd(2*i-1)
|
||||
else if ((index(pos_string,'-') > 0).and.(index(pos_string,'inf')>0)) then
|
||||
!Now extract the number we are reducing from infinity
|
||||
read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
|
||||
pos = box_bd(2*i) - pos
|
||||
else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then
|
||||
!Now extract the number we are reducing from infinity
|
||||
read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
|
||||
pos = box_bd(2*i-1) + pos
|
||||
else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then
|
||||
!Now extract the number we are reducing from infinity
|
||||
read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos
|
||||
pos = (box_bd(2*i)-box_bd(2*i-1))*pos
|
||||
else
|
||||
read(pos_string, *, iostat=iospara) pos
|
||||
end if
|
||||
|
||||
if (iospara > 0) then
|
||||
print *, "Error reading position argument ", trim(adjustl(pos_string)), ". Please reformat and try again."
|
||||
end if
|
||||
end subroutine parse_pos
|
||||
|
||||
subroutine heapsort(a)
|
||||
|
||||
integer, intent(in out) :: a(0:)
|
||||
integer :: start, n, bottom
|
||||
real :: temp
|
||||
|
||||
n = size(a)
|
||||
do start = (n - 2) / 2, 0, -1
|
||||
call siftdown(a, start, n);
|
||||
end do
|
||||
|
||||
do bottom = n - 1, 1, -1
|
||||
temp = a(0)
|
||||
a(0) = a(bottom)
|
||||
a(bottom) = temp;
|
||||
call siftdown(a, 0, bottom)
|
||||
end do
|
||||
|
||||
end subroutine heapsort
|
||||
|
||||
subroutine siftdown(a, start, bottom)
|
||||
|
||||
integer, intent(in out) :: a(0:)
|
||||
integer, intent(in) :: start, bottom
|
||||
integer :: child, root
|
||||
real :: temp
|
||||
|
||||
root = start
|
||||
do while(root*2 + 1 < bottom)
|
||||
child = root * 2 + 1
|
||||
|
||||
if (child + 1 < bottom) then
|
||||
if (a(child) < a(child+1)) child = child + 1
|
||||
end if
|
||||
|
||||
if (a(root) < a(child)) then
|
||||
temp = a(child)
|
||||
a(child) = a (root)
|
||||
a(root) = temp
|
||||
root = child
|
||||
else
|
||||
return
|
||||
end if
|
||||
end do
|
||||
|
||||
end subroutine siftdown
|
||||
end module subroutines
|
||||
Loading…
Add table
Add a link
Reference in a new issue