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

Add group shape id

Browse Source
master
Alex Selimov 4 years ago
parent f6c0cf93f2
commit 284fce3c17
1 changed files with 152 additions and 39 deletions
Show Stats Download Patch File Download Diff File

191
src/opt_group.f90
Unescape View File

@ -54,8 +54,8 @@ module opt_group
!Parse the group command
integer, intent(inout) :: arg_pos
integer :: i,arglen
character(len=100) :: textholder
integer :: i, j, arglen, in_num
character(len=100) :: textholder, type_spec
real(kind=dp) bwidth, wheight
!Parse type and shape command
@ -85,51 +85,162 @@ module opt_group
end do
case('wedge')
arg_pos = arg_pos + 1
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing normal dim in group wedge 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 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 normal dim in group wedge command"
read(textholder,*) dim1
arg_pos = arg_pos + 1
if (arglen==0) STOP "Missing centroid in group wedge 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 command"
read(textholder,*) bwidth
!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 shaped 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
!Read the ID type shape for group
case('id')
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.
select case(trim(adjustl(type)))
case('atoms')
!Read number of ids
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing normal dim in group wedge command"
read(textholder,*) dim2
if(arglen == 0) then
print *, "Missing number of input atoms for group id"
end if
read(textholder,*) in_num
!allocate arrays
allocate(atom_index(in_num))
!Read ids
do i = 1, in_num
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if(arglen == 0) then
print *, "Missing atom id in group atom id"
stop 3
end if
read(textholder,*) atom_index(i)
end do
group_atom_num = in_num
do i = 1, 3
case('elements')
!Read number of ids
call get_command_argument(arg_pos, textholder, arglen)
if(arglen == 0) then
print *, "Missing number of input elements for group id"
end if
read(textholder, *) in_num
!allocate arrays
allocate(element_index(in_num))
!Read ids
do i = 1, in_num
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing centroid in group wedge command"
call parse_pos(i, textholder, centroid(i))
end do
if(arglen == 0) then
print *, "Missing element id in group element id"
stop 3
end if
read(textholder,*) element_index(i)
end do
group_ele_num = in_num
case('both')
!We repeat this code twice, once for the atoms and once for the elements
allocate(element_index(1024),atom_index(1024))
do i = 1, 2
!Read the first id type (either atoms or elements)
call get_command_argument(arg_pos, type_spec, arglen)
if (arglen == 0) then
print *, "Missing type specifier in group id command"
stop 3
end if
!Now read the first in_num
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if(arglen == 0) then
print *, "Missing input number in group_id"
stop 3
end if
read(textholder, *) in_num
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen==0) STOP "Missing base width in group wedge command"
read(textholder,*) bwidth
select case(trim(adjustl(type_spec)))
case('atoms','atom')
if (group_atom_num > 0) then
print *, "Atoms specifier used more than once in group id command with type both, either use type ", &
"atoms or include elements identifier"
stop 3
!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 shaped group"
do j = 1, in_num
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen == 0) then
print *, "Missing atom id in group atom id"
stop 3
end if
read(textholder, *) atom_index(j)
end do
group_atom_num = in_num
end if
case('elements','element')
if (group_ele_num > 0) then
print *, "Elements specifier used more than once in group id command with type both, either use type ", &
"elements or include atoms identifier"
stop 3
do j = 1, in_num
arg_pos = arg_pos + 1
call get_command_argument(arg_pos, textholder, arglen)
if (arglen == 0) then
print *, "Missing element id in group element id"
stop 3
end if
read(textholder, *) element_index(j)
end do
group_ele_num = in_num
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
end select
if(i ==1) arg_pos = arg_pos + 1
end do
end select
case default
print *, "Group shape ", trim(adjustl(shape)), " is not currently accepted. Please check documentation ", &
"for accepted group shapes."
@ -183,6 +294,9 @@ module opt_group
print *, "Group has block shape with boundaries: ", block_bd
case ('crack')
print *, "Group has crack shape with dim1", dim1, "and dim2", dim2, "and vertices ", vertices
case('id')
print *, 'Group contains ', group_ele_num, " elements and ", group_atom_num, " atoms."
return
end select
!Allocate variables to arbitrary size
@ -274,7 +388,6 @@ module opt_group
real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3)
!Refining to atoms
if(group_ele_num > 0) then
orig_atom_num = atom_num
!Estimate number of atoms we are adding, this doesn't have to be exact

Write Preview
Loading…
Cancel
Save