diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 9a01f1f..dbc984a 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -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