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

Lots of accumulated fixes to various functionality with cacmb

Browse Source
This commit is contained in:
Alex Selimov 2022-01-13 12:16:50 -05:00
parent 1114a46c60
commit d9920bbef0
8 changed files with 251 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
src/elements.f90
View File

@ -143,6 +143,7 @@ module elements
max_basisnum = 0 max_basisnum = 0
basisnum(:) = 0 basisnum(:) = 0
basis_type(:,:) = 0
ng_node(:) = 0 ng_node(:) = 0
node_num = 0 node_num = 0
node_atoms = 0 node_atoms = 0
@ -308,7 +309,7 @@ module elements
!Subroutine which adds an element to the element arrays !Subroutine which adds an element to the element arrays
integer, intent(in) :: size, lat, sbox, tag integer, intent(in) :: size, lat, sbox, tag
character(len=100), intent(in) :: type character(len=100), intent(in) :: type
real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node) real(kind=dp), intent(in) :: r(:,:,:)
integer :: newtag integer :: newtag
@ -435,8 +436,8 @@ module elements
real(kind=dp), dimension(:,:,:), intent(in) :: r_in !Nodal positions real(kind=dp), dimension(:,:,:), intent(in) :: r_in !Nodal positions
integer, dimension(:), intent(out) :: type_interp !Interpolated atomic positions integer, dimension(:), intent(out) :: type_interp !Interpolated atomic positions
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(max_basisnum, max_ng_node), f(3, max_basisnum, max_ng_node), & real(kind = dp), optional, intent(in) :: eng(:,:), f(:,:,:), &
v(6, max_basisnum, max_ng_node) v(:,:,:)
real(kind=dp), dimension(10, max_basisnum*max_esize**3), optional,intent(out) :: data_interp !Interpolated atomic positions real(kind=dp), dimension(10, max_basisnum*max_esize**3), optional,intent(out) :: data_interp !Interpolated atomic positions
!Internal variables !Internal variables
@ -463,7 +464,7 @@ module elements
end select end select
select case(trim(adjustl(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
@ -657,11 +658,45 @@ module elements
subroutine wrap_atoms subroutine wrap_atoms
!This subroutine wraps atoms back into the simulation cell if they have exited for any reason !This subroutine wraps atoms back into the simulation cell if they have exited for any reason
integer :: i integer :: i
do i = 1, atom_num do i = 1, atom_num
call apply_periodic(r_atom(:,i)) call apply_periodic(r_atom(:,i))
end do end do
end subroutine wrap_atoms 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)
do j = 1, 3
box_len(j) = box_bd(2*j) - box_bd(2*j-1)
end do
print *, box_bd
do i = 1, ele_num
node_in_bd = 0
do inod = 1, ng_node(lat_ele(i))
do ibasis = 1, basisnum(lat_ele(i))
do j = 1, 3
if(r_node(j,ibasis,inod,i) < box_bd(2*j-1)) then
node_in_bd(j,inod) = 1
else if(r_node(j,ibasis,inod,i) > box_bd(2*j)) then
node_in_bd(j,inod) = -1
end if
end do
end do
end do
do j = 1, 3
if(all(node_in_bd(j,:) == 1)) then
r_node(j,:,:,i) = r_node(j,:,:,i) + box_len(j)
else if(all(node_in_bd(j,:) == -1)) then
r_node(j,:,:,i) = r_node(j,:,:,i) - box_len(j)
end if
end do
end do
end subroutine wrap_elements
subroutine def_new_box subroutine def_new_box
!This subroutine calculates new box boundaries based on minimum and maximum nodal/atomic positions !This subroutine calculates new box boundaries based on minimum and maximum nodal/atomic positions
@ -1002,10 +1037,12 @@ do i = 1, atom_num
real(kind=dp), intent(in) :: eng(max_basisnum, max_ng_node), & real(kind=dp), intent(in) :: eng(max_basisnum, max_ng_node), &
force(3,max_basisnum, max_ng_node), & force(3,max_basisnum, max_ng_node), &
virial(6,max_basisnum,max_ng_node) virial(6,max_basisnum,max_ng_node)
integer :: inod
vflag = .true. vflag = .true.
energy_node(:,:,ie) = eng energy_node(:,:,ie) = eng
force_node(:,:,:,ie) = force force_node(:,:,:,ie) = force
virial_node(:,:,:,ie) = virial virial_node(:,:,:,ie) = virial
return return
end subroutine add_element_data end subroutine add_element_data

6
src/io.f90
View File

@ -275,11 +275,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)') i, type_atom(i), r_atom(:,i), energy_atom(i), force_atom(:,i), virial_atom(:,i) write(11, '(2i16, 13f23.15)') tag_atom(i), type_atom(i), r_atom(:,i), energy_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)') 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
@ -1146,6 +1147,7 @@ module io
!Now set the nodes and basisnum !Now set the nodes and basisnum
max_ng_node = maxval(ng_node(1:lattice_types)) max_ng_node = maxval(ng_node(1:lattice_types))
max_basisnum = maxval(basisnum(1:lattice_types)) max_basisnum = maxval(basisnum(1:lattice_types))
return return
end subroutine read_pycac_out end subroutine read_pycac_out

5
src/main.f90
View File

@ -113,7 +113,10 @@ program main
end do end do
!If wrap flag was passed then call the wrap atoms command !If wrap flag was passed then call the wrap atoms command
if(wrap_flag) call wrap_atoms if(wrap_flag) then
call wrap_atoms
call wrap_elements
end if
!If we called the boundary command then we adjust the box bounds !If we called the boundary command then we adjust the box bounds
if(bound_called) call def_new_box if(bound_called) call def_new_box

117
src/mode_create.f90
View File

@ -15,7 +15,7 @@ module mode_create
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)
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) basis_pos(3,10), esize_nums, esize_index(10)
logical :: dup_flag, dim_flag, efill 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(:,:)
integer, allocatable :: elat(:) integer, allocatable :: elat(:)
@ -79,10 +79,10 @@ module mode_create
do i = 1, 3 do i = 1, 3
box_bd(2*i) = maxval(box_vert(i,:)) - 0.25_dp*lattice_space(i) box_bd(2*i) = maxval(box_vert(i,:)) - 0.25_dp*lattice_space(i)
box_bd(2*i-1) = origin(i)-0.25_dp*lattice_space(i) box_bd(2*i-1) = origin(i)-0.25_dp*lattice_space(i)
box_len(i) = box_bd(2*i) - box_bd(2*i-1)
end do end do
else if(dim_flag) then else if(dim_flag) then
!As a note everything is defined so that the lattice parameter is multiplied in at the end !As a note everything is defined so that the lattice parameter is multiplied in at the end
!so we have to divide all the real Angstroms units by the lattice parameter !so we have to divide all the real Angstroms units by the lattice parameter
@ -131,11 +131,11 @@ module mode_create
!Call the build function with the correct transformation matrix !Call the build function with the correct transformation matrix
select case(trim(adjustl(element_type))) select case(trim(adjustl(element_type)))
case('fcc') case('fcc')
call build_with_rhomb(box_lat_vert, fcc_mat, 8) call build_with_rhomb(box_lat_vert, fcc_mat, 8, fcc_inv)
case('bcc') case('bcc')
call build_with_rhomb(box_lat_vert, bcc_mat, 8) call build_with_rhomb(box_lat_vert, bcc_mat, 8, bcc_inv)
case('20fcc') case('20fcc')
call build_with_rhomb(box_lat_vert, fcc_mat, 20) call build_with_rhomb(box_lat_vert, fcc_mat, 20, fcc_inv)
case default case default
print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ",& print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ",&
"element type" "element type"
@ -176,6 +176,11 @@ module mode_create
sub_box_num = 1 sub_box_num = 1
sub_box_ori(:,:,1) = orient sub_box_ori(:,:,1) = orient
sub_box_bd(:,1) = box_bd sub_box_bd(:,1) = box_bd
!If any elements are fully outside the box then wrap them back in
if (any(crossb)) then
call wrap_elements
end if
end subroutine create end subroutine create
!This subroutine parses the command and pulls out information needed for mode_create !This subroutine parses the command and pulls out information needed for mode_create
subroutine parse_command(arg_pos) subroutine parse_command(arg_pos)
@ -183,9 +188,9 @@ module mode_create
integer, intent(out) :: arg_pos integer, intent(out) :: arg_pos
integer :: ori_pos, i, j, arglen, stat integer :: ori_pos, i, j, arglen, stat
character(len=100) :: textholder character(len=100) :: textholder
character(len=20) :: orient_string character(len=100) :: orient_string
character(len=2) :: btype character(len=2) :: btype
logical :: isortho, isrighthanded logical :: isortho, isrighthanded, bool
!Pull out all required positional arguments !Pull out all required positional arguments
@ -203,6 +208,7 @@ module mode_create
call get_command_argument(5, textholder, arglen) call get_command_argument(5, textholder, arglen)
if(arglen==0) STOP "Esize missing in mode create" if(arglen==0) STOP "Esize missing in mode create"
read(textholder, *, iostat=stat) esize read(textholder, *, iostat=stat) esize
max_esize = esize
if(stat > 0) STOP "Error reading esize" if(stat > 0) STOP "Error reading esize"
arg_pos = 6 arg_pos = 6
@ -250,6 +256,12 @@ module mode_create
read(textholder, *) origin(i) read(textholder, *) origin(i)
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
end do end do
case('crossb')
do i = 1, 3
call get_command_argument(arg_pos, textholder)
read(textholder, *) crossb(i)
arg_pos = arg_pos + 1
end do
case('basis') case('basis')
call get_command_argument(arg_pos, textholder) call get_command_argument(arg_pos, textholder)
read(textholder, *) basisnum(1) read(textholder, *) basisnum(1)
@ -275,11 +287,13 @@ module mode_create
end select end select
end do end do
!Calculate the lattice periodicity length in lattice units !Calculate the lattice periodicity length in lattice units
do i = 1, 3 do i = 1, 3
lattice_space(i) = norm2(orient(i,:)) lattice_space(i) = norm2(orient(i,:))
end do end do
!Now normalize the orientation matrix
orient = matrix_normal(orient,3)
!Check special periodicity relations !Check special periodicity relations
select case(trim(adjustl(element_type))) select case(trim(adjustl(element_type)))
case('fcc', '20fcc') case('fcc', '20fcc')
@ -309,8 +323,6 @@ module mode_create
end if end if
end do end do
end select end select
!Now normalize the orientation matrix
orient = matrix_normal(orient,3)
!Now check these to make sure they are right handed and orthogonal !Now check these to make sure they are right handed and orthogonal
call check_right_ortho(orient, isortho, isrighthanded) call check_right_ortho(orient, isortho, isrighthanded)
if (.not.isortho) then if (.not.isortho) then
@ -330,21 +342,22 @@ module mode_create
call add_atom_type(name, basis_type(1,1)) !If basis command not defined then we use name as the atom_name call add_atom_type(name, 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
!
end subroutine end subroutine
subroutine build_with_rhomb(box_in_lat, transform_matrix, nn) subroutine build_with_rhomb(box_in_lat, transform_matrix, nn, transform_inverse)
!This subroutine returns all the lattice points in the box in r_lat !This subroutine returns all the lattice points in the box in r_lat
!Inputs !Inputs
integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space integer, dimension(3,8), intent(in) :: box_in_lat !The box vertices transformed to lattice space
real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space real(kind=dp), dimension(3,3), intent(in) :: transform_matrix !The transformation matrix from lattice_space to real space
real(kind=dp), dimension(3,3), intent(in) :: transform_inverse !The inverse transform
integer, intent(in) :: nn integer, intent(in) :: nn
!Internal variables !Internal variables
integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,nn), rzero(3), efill_size, & integer :: i, inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz, numlatpoints, ele(3,nn), rzero(3), efill_size, &
vlat(3), temp_lat(3,nn), m, n, o, j, k, nump_ele, efill_temp_lat(3,nn), filzero(3), & vlat(3), temp_lat(3,nn), m, n, o, j, k, nump_ele, efill_temp_lat(3,nn), filzero(3), &
bd_ele_lat(6), efill_ele(3,nn), ebd(6) bd_ele_lat(6), efill_ele(3,nn), ebd(6), shift_vec(3), type_interp(max_basisnum*max_esize**3)
real(kind=dp) :: v(3), temp_nodes(3,1,nn), r(3), centroid_bd(6) real(kind=dp) :: v(3), temp_nodes(3,1,nn), r(3), centroid_bd(6), vreal(3), r_interp(3, max_basisnum*max_esize**3)
logical, allocatable :: lat_points(:,:,:) logical, allocatable :: lat_points(:,:,:)
logical :: node_in_bd(nn), add, lat_points_ele(esize,esize,esize), intersect_bd(3) logical :: node_in_bd(nn), add, lat_points_ele(esize,esize,esize), intersect_bd(3)
@ -422,6 +435,7 @@ module mode_create
end do end do
end do end do
!Now we redefine bd_in_lat The first 3 indices contains limits for the lat_points array !Now we redefine bd_in_lat The first 3 indices contains limits for the lat_points array
bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10 bd_in_array(1) = bd_in_lat(2) - bd_in_lat(1) + 10
bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10 bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10
@ -493,6 +507,31 @@ module mode_create
!If within array boundaries check to see if it is a lattice point !If within array boundaries check to see if it is a lattice point
else if(lat_points(vlat(1),vlat(2),vlat(3))) then else if(lat_points(vlat(1),vlat(2),vlat(3))) then
node_in_bd(inod) = .true. node_in_bd(inod) = .true.
else if(any(crossb)) then
vreal=0
do i = 1, 3
if(crossb(i)) then
if(temp_nodes(i,1,inod) < box_bd(2*i-1)) then
vreal(i) = temp_nodes(i,1,inod)+box_len(i)
else if(temp_nodes(i,1,inod) > box_bd(2*i)) then
vreal(i) = temp_nodes(i,1,inod)-box_len(i)
else
vreal(i) = temp_nodes(i,1,inod)
end if
else
vreal(i) = temp_nodes(i,1,inod)
end if
end do
v = matmul(transform_inverse, matmul(orient_inv, vreal))
do i = 1, 3
vlat(i) = nint(v(i) - bd_in_lat(2*i-1)+5)
end do
if(any(vlat > shape(lat_points)).or.any(vlat < 1)) then
continue
!If within array boundaries check to see if it is a lattice point
else if(lat_points(vlat(1),vlat(2),vlat(3))) then
node_in_bd(inod) = .true.
end if
end if end if
end do end do
@ -502,14 +541,50 @@ module mode_create
lat_ele_num = lat_ele_num+1 lat_ele_num = lat_ele_num+1
r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
elat(lat_ele_num) = esize elat(lat_ele_num) = esize
!Now set all the lattice points contained within an element to false if(any(crossb)) then
do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:)) call interpolate_atoms('fcc', esize, 0, temp_nodes, type_interp, r_interp)
do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:)) j= 0
do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:)) do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
lat_points(m,n,o) = .false. do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))
do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:))
j=j+1
do i = 1, 3
if(crossb(i)) then
if(r_interp(i,j) < box_bd(2*i-1)) then
vreal(i) = r_interp(i,j)+box_len(i)
else if(r_interp(i,j) > box_bd(2*i)) then
vreal(i) = r_interp(i,j)-box_len(i)
else
vreal(i) = r_interp(i,j)
end if
else
vreal(i) = r_interp(i,j)
end if
end do
v = matmul(transform_inverse, matmul(orient_inv, vreal))
do i = 1, 3
vlat(i) = nint(v(i) - bd_in_lat(2*i-1)+5)
end do
if(lat_points(vlat(1), vlat(2), vlat(3))) then
lat_points(vlat(1), vlat(2), vlat(3)) = .false.
else
print *, "Lat points should be true not false"
stop 2
end if
end do
end do end do
end do end do
end do else
!Now set all the lattice points contained within an element to false
do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))
do m = minval(temp_lat(1,:)), maxval(temp_lat(1,:))
lat_points(m,n,o) = .false.
end do
end do
end do
end if
!If any nodes are in the boundary and we want to efill then run the efill code !If any nodes are in the boundary and we want to efill then run the efill code
else if(any(node_in_bd).and.efill) then else if(any(node_in_bd).and.efill) then

14
src/mode_metric.f90
View File

@ -63,6 +63,7 @@ module mode_metric
!Now loop over new files !Now loop over new files
do i = 2, nfiles do i = 2, nfiles
call read_in(i, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd) call read_in(i, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd)
call convert_positions(r_curr, np_temp) call convert_positions(r_curr, np_temp)
if (npreal /= np_temp) then if (npreal /= np_temp) then
@ -196,6 +197,7 @@ module mode_metric
rout = -huge(1.0_dp) rout = -huge(1.0_dp)
if(atom_num > 0) then if(atom_num > 0) then
do i = 1, atom_num do i = 1, atom_num
call apply_periodic(r_atom(:,i))
rout(:,tag_atom(i)) = r_atom(:,i) rout(:,tag_atom(i)) = r_atom(:,i)
npoints = npoints + 1 npoints = npoints + 1
end do end do
@ -205,8 +207,11 @@ module mode_metric
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))
call apply_periodic(r_node(:,ibasis,inod,i))
rout(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis) & rout(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis) &
= r_node(:,ibasis,inod,i) = r_node(:,ibasis,inod,i)
npoints = npoints + 1 npoints = npoints + 1
end do end do
end do end do
@ -226,14 +231,14 @@ module mode_metric
select case(metric_type) select case(metric_type)
case('def_grad') case('def_grad')
write(11,*) "type element x y z F11 F12 F13 F21 F22 F23 F31 F32 F33" write(11,*) "id type element x y z F11 F12 F13 F21 F22 F23 F31 F32 F33"
case('microrotation') case('microrotation')
write(11,*) "type element x y z micro1 micro2 micro3 norm2(micro)" write(11,*) "id type element x y z micro1 micro2 micro3 norm2(micro)"
end select end select
if(atom_num > 0) then if(atom_num > 0) then
do i = 1, atom_num do i = 1, atom_num
write(11,*) type_atom(i), 0, r_atom(:,i), met(:,tag_atom(i)) write(11,*) tag_atom(i), type_atom(i), 0, r_atom(:,i), met(:,tag_atom(i))
end do end do
end if end if
@ -241,7 +246,8 @@ module mode_metric
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,*) basis_type(ibasis,lat_ele(i)), 1, r_node(:,ibasis,inod,i), & write(11,*) atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis, &
basis_type(ibasis,lat_ele(i)), 1, r_node(:,ibasis,inod,i), &
met(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis) met(:, atom_num+(tag_ele(i)-1)*max_ng_node*max_basisnum + (inod-1)*max_basisnum + ibasis)
end do end do
end do end do

1
src/neighbors.f90
View File

@ -54,7 +54,6 @@ module neighbors
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 belongs to !c is the position of the cell that the point belongs to
if(i == 11651203) print *, r_list(:,i), (r_list(1,i) < box_bd(1))
do j = 1, 3 do j = 1, 3
c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off)) + 1 c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off)) + 1
end do end do

74
src/opt_group.f90
View File

@ -15,7 +15,7 @@ module opt_group
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, &
insert_lattice, s_fractions(10) insert_lattice, s_fractions(10), max_height
logical :: displace, delete, max_remesh, refine, group_nodes, flip, efill, refinefill, alloy logical :: displace, delete, max_remesh, refine, group_nodes, flip, efill, refinefill, alloy
@ -109,11 +109,11 @@ module opt_group
call get_command_argument(arg_pos, type, arglen) call get_command_argument(arg_pos, type, arglen)
if (arglen==0) STOP "Missing select_type in group command" if (arglen==0) STOP "Missing select_type in group command"
select case(trim(adjustl(type))) select case(trim(adjustl(type)))
case('atoms', 'elements', 'both') case('atoms', 'elements', 'all')
continue continue
case default case default
print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", & print *, "Select_type ", trim(adjustl(type)), " is not an accept group selection criteria. ", &
"Please select from atoms, elements, or both." "Please select from atoms, elements, or all."
end select end select
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
@ -177,6 +177,57 @@ module opt_group
end if end if
end do end do
case('wedge_cut')
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing normal dim in group wedge_cut command"
read(textholder,*) dim1
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing normal dim in group wedge_cut command"
read(textholder,*) dim2
do i = 1, 3
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing centroid in group wedge_cut command"
call parse_pos(i, textholder, centroid(i))
end do
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing base width in group wedge_cut command"
read(textholder,*) bwidth
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing base width in group wedge_cut command"
read(textholder,*) max_height
!Calculate the vertex positions
vertices(:,1) = centroid
vertices(dim2,1) = 0.0_dp
do i = 1, 3
if (i == dim1) then
if (bwidth > 0) then
vertices(i,2) = box_bd(2*i)
vertices(i,3) = box_bd(2*i)
else if (bwidth < 0) then
vertices(i,2) = box_bd(2*i-1)
vertices(i,3) = box_bd(2*i-1)
else
print *, "bwidth cannot be 0 in wedge gshaped group"
stop 3
end if
else if (i == dim2) then
vertices(i,2) = 0.0_dp
vertices(i,3) = 0.0_dp
else
vertices(i,2) = centroid(i) + bwidth
vertices(i,3) = centroid(i) - bwidth
end if
end do
case('notch') case('notch')
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
@ -239,7 +290,7 @@ module opt_group
!Read the ID type shape for group !Read the ID type shape for group
case('id') case('id')
arg_pos = arg_pos + 1 arg_pos = arg_pos + 1
!For this type we have to call different options if type is atoms, elements, or both. Both is the most complex. !For this type we have to call different options if type is atoms, elements, or all. all is the most complex.
select case(trim(adjustl(type))) select case(trim(adjustl(type)))
case('atoms') case('atoms')
!Read number of ids !Read number of ids
@ -285,7 +336,7 @@ module opt_group
end do end do
group_ele_num = in_num group_ele_num = in_num
case('both') case('all')
!We repeat this code twice, once for the atoms and once for the elements !We repeat this code twice, once for the atoms and once for the elements
allocate(element_index(1024),atom_index(1024)) allocate(element_index(1024),atom_index(1024))
do i = 1, 2 do i = 1, 2
@ -308,7 +359,7 @@ module opt_group
select case(trim(adjustl(type_spec))) select case(trim(adjustl(type_spec)))
case('atoms','atom') case('atoms','atom')
if (group_atom_num > 0) then if (group_atom_num > 0) then
print *, "Atoms specifier used more than once in group id command with type both, either use type ", & print *, "Atoms specifier used more than once in group id command with type all, either use type ", &
"atoms or include elements identifier" "atoms or include elements identifier"
stop 3 stop 3
@ -326,7 +377,7 @@ module opt_group
case('elements','element') case('elements','element')
if (group_ele_num > 0) then if (group_ele_num > 0) then
print *, "Elements specifier used more than once in group id command with type both, either use type ",& print *, "Elements specifier used more than once in group id command with type all, either use type ",&
"elements or include atoms identifier" "elements or include atoms identifier"
stop 3 stop 3
@ -527,7 +578,7 @@ module opt_group
allocate(element_index(1024), atom_index(1024)) allocate(element_index(1024), atom_index(1024))
!Check the type to see whether we need to find the elements within the group !Check the type to see whether we need to find the elements within the group
select case(trim(adjustl(type))) select case(trim(adjustl(type)))
case('elements', 'both') case('elements', 'all')
do i = 1, ele_num do i = 1, ele_num
if(in_group_ele(size_ele(i), lat_ele(i), r_node(:,:,:,i))) then if(in_group_ele(size_ele(i), lat_ele(i), r_node(:,:,:,i))) then
group_ele_num = group_ele_num + 1 group_ele_num = group_ele_num + 1
@ -559,7 +610,7 @@ module opt_group
end select end select
!Check the type to see if we need to find the atoms within the group !Check the type to see if we need to find the atoms within the group
select case(trim(adjustl(type))) select case(trim(adjustl(type)))
case('atoms','both') 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(r_atom(:,i)).neqv.flip) then
group_atom_num = group_atom_num + 1 group_atom_num = group_atom_num + 1
@ -1151,10 +1202,8 @@ module opt_group
call add_atom_type(species_type(i), type_map(i)) call add_atom_type(species_type(i), 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
print *, s_fractions
!Now randomize the atom types !Now randomize the atom types
added_types = 0 added_types = 0
print *, type_map(j)
do i = 1, group_atom_num do i = 1, group_atom_num
ia = atom_index(i) ia = atom_index(i)
call random_number(rand) call random_number(rand)
@ -1184,6 +1233,9 @@ module opt_group
in_group=in_block_bd(r,block_bd) in_group=in_block_bd(r,block_bd)
case('wedge') case('wedge')
in_group = in_wedge_bd(r,vertices) in_group = in_wedge_bd(r,vertices)
case('wedge_cut')
in_group = in_wedge_bd(r,vertices)
in_group = in_group.and.((abs(r(dim1) - vertices(dim1,2)) < max_height))
case('notch') case('notch')
do i = 1, 3 do i = 1, 3
if (.not.((dim1==i).or.(dim2==i))) dim3 = i if (.not.((dim1==i).or.(dim2==i))) dim3 = i

48
src/subroutines.f90
View File

@ -2,6 +2,8 @@ module subroutines
use parameters use parameters
use functions use functions
use box use box
use str
implicit none implicit none
integer :: allostat, deallostat integer :: allostat, deallostat
@ -149,29 +151,43 @@ module subroutines
subroutine parse_ori_vec(ori_string, ori_vec) subroutine parse_ori_vec(ori_string, ori_vec)
!This subroutine parses a string to vector in the format [ijk] !This subroutine parses a string to vector in the format [ijk]
character(len=20), intent(in) :: ori_string character(*), intent(in) :: ori_string
real(kind=dp), dimension(3), intent(out) :: ori_vec real(kind=dp), dimension(3), intent(out) :: ori_vec
real(kind=dp) :: prefactor
integer :: start, fin
integer :: i, ori_pos, stat integer :: i, ori_pos, stat
ori_pos=2 do i = 1, len(ori_string)
do i = 1,3 if(ori_string(i:i) == "[") then
do while(ori_string(ori_pos:ori_pos)==' ') start=i+1
ori_pos=ori_pos+1 else if(ori_string(i:i) == "]") then
end do fin = i-1
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 if
end do end do
ori_pos=start
if (tok_count(ori_string(start:fin)) == 3) then
read(ori_string(start:fin),*) ori_vec(:)
else
do i = 1,3
do while(ori_string(ori_pos:ori_pos)==' ')
ori_pos=ori_pos+1
end do
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
end if
return return
end subroutine parse_ori_vec end subroutine parse_ori_vec