diff --git a/src/elements.f90 b/src/elements.f90 index f371074..1dc58bf 100644 --- a/src/elements.f90 +++ b/src/elements.f90 @@ -143,6 +143,7 @@ module elements max_basisnum = 0 basisnum(:) = 0 + basis_type(:,:) = 0 ng_node(:) = 0 node_num = 0 node_atoms = 0 @@ -308,7 +309,7 @@ module elements !Subroutine which adds an element to the element arrays integer, intent(in) :: size, lat, sbox, tag 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 @@ -435,8 +436,8 @@ module elements real(kind=dp), dimension(:,:,:), intent(in) :: r_in !Nodal 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), optional, intent(in) :: eng(max_basisnum, max_ng_node), f(3, max_basisnum, max_ng_node), & - v(6, max_basisnum, max_ng_node) + 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 !Internal variables @@ -463,7 +464,7 @@ module elements end select - select case(trim(adjustl(type))) + select case(type) case('fcc','bcc') !Now loop over all the possible sites do it = 1, esize @@ -657,11 +658,45 @@ module elements subroutine wrap_atoms !This subroutine wraps atoms back into the simulation cell if they have exited for any reason + integer :: i do i = 1, atom_num call apply_periodic(r_atom(:,i)) 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) + + 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 !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), & force(3,max_basisnum, max_ng_node), & virial(6,max_basisnum,max_ng_node) + integer :: inod vflag = .true. energy_node(:,:,ie) = eng force_node(:,:,:,ie) = force virial_node(:,:,:,ie) = virial + return end subroutine add_element_data diff --git a/src/io.f90 b/src/io.f90 index 6f6fec6..de8a7bc 100644 --- a/src/io.f90 +++ b/src/io.f90 @@ -275,11 +275,12 @@ module io if (write_dat) then 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 else 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 @@ -1146,6 +1147,7 @@ module io !Now set the nodes and basisnum max_ng_node = maxval(ng_node(1:lattice_types)) max_basisnum = maxval(basisnum(1:lattice_types)) + return end subroutine read_pycac_out diff --git a/src/main.f90 b/src/main.f90 index 92afba1..73c838d 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -113,7 +113,10 @@ program main end do !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(bound_called) call def_new_box diff --git a/src/mode_create.f90 b/src/mode_create.f90 index 12bb2b8..bbebaa2 100644 --- a/src/mode_create.f90 +++ b/src/mode_create.f90 @@ -15,7 +15,7 @@ module mode_create 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), & 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(:,:) integer, allocatable :: elat(:) @@ -79,10 +79,10 @@ module mode_create do i = 1, 3 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_len(i) = box_bd(2*i) - box_bd(2*i-1) end do - else if(dim_flag) then !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 @@ -131,11 +131,11 @@ module mode_create !Call the build function with the correct transformation matrix select case(trim(adjustl(element_type))) 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') - call build_with_rhomb(box_lat_vert, bcc_mat, 8) + call build_with_rhomb(box_lat_vert, bcc_mat, 8, bcc_inv) 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 print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ",& "element type" @@ -176,6 +176,11 @@ module mode_create sub_box_num = 1 sub_box_ori(:,:,1) = orient 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 !This subroutine parses the command and pulls out information needed for mode_create subroutine parse_command(arg_pos) @@ -183,9 +188,9 @@ module mode_create integer, intent(out) :: arg_pos integer :: ori_pos, i, j, arglen, stat character(len=100) :: textholder - character(len=20) :: orient_string + character(len=100) :: orient_string character(len=2) :: btype - logical :: isortho, isrighthanded + logical :: isortho, isrighthanded, bool !Pull out all required positional arguments @@ -203,6 +208,7 @@ module mode_create call get_command_argument(5, textholder, arglen) if(arglen==0) STOP "Esize missing in mode create" read(textholder, *, iostat=stat) esize + max_esize = esize if(stat > 0) STOP "Error reading esize" arg_pos = 6 @@ -250,6 +256,12 @@ module mode_create read(textholder, *) origin(i) arg_pos = arg_pos + 1 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') call get_command_argument(arg_pos, textholder) read(textholder, *) basisnum(1) @@ -275,11 +287,13 @@ module mode_create end select end do + !Calculate the lattice periodicity length in lattice units do i = 1, 3 lattice_space(i) = norm2(orient(i,:)) end do - + !Now normalize the orientation matrix + orient = matrix_normal(orient,3) !Check special periodicity relations select case(trim(adjustl(element_type))) case('fcc', '20fcc') @@ -309,8 +323,6 @@ module mode_create end if end do end select - !Now normalize the orientation matrix - orient = matrix_normal(orient,3) !Now check these to make sure they are right handed and orthogonal call check_right_ortho(orient, isortho, isrighthanded) 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 basis_pos(:,1) = 0.0_dp end if - ! + 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 !Inputs 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_inverse !The inverse transform integer, intent(in) :: nn !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), & - bd_ele_lat(6), efill_ele(3,nn), ebd(6) - real(kind=dp) :: v(3), temp_nodes(3,1,nn), r(3), centroid_bd(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), vreal(3), r_interp(3, max_basisnum*max_esize**3) logical, allocatable :: lat_points(:,:,:) 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 + !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(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 else if(lat_points(vlat(1),vlat(2),vlat(3))) then 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 do @@ -502,14 +541,50 @@ module mode_create lat_ele_num = lat_ele_num+1 r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:) elat(lat_ele_num) = esize - !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 + if(any(crossb)) then + call interpolate_atoms('fcc', esize, 0, temp_nodes, type_interp, r_interp) + j= 0 + 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,:)) + 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 + 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 do + end if !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 diff --git a/src/mode_metric.f90 b/src/mode_metric.f90 index 8d0b6a5..dab29f5 100644 --- a/src/mode_metric.f90 +++ b/src/mode_metric.f90 @@ -63,6 +63,7 @@ module mode_metric !Now loop over new files do i = 2, nfiles + call read_in(i, (/ 0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd) call convert_positions(r_curr, np_temp) if (npreal /= np_temp) then @@ -196,6 +197,7 @@ module mode_metric rout = -huge(1.0_dp) if(atom_num > 0) then do i = 1, atom_num + call apply_periodic(r_atom(:,i)) rout(:,tag_atom(i)) = r_atom(:,i) npoints = npoints + 1 end do @@ -205,8 +207,11 @@ module mode_metric do i = 1, ele_num do inod = 1, ng_node(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) & = r_node(:,ibasis,inod,i) + + npoints = npoints + 1 end do end do @@ -226,14 +231,14 @@ module mode_metric select case(metric_type) 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') - 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 if(atom_num > 0) then 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 if @@ -241,7 +246,8 @@ module mode_metric do i = 1, ele_num do inod = 1, ng_node(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) end do end do diff --git a/src/neighbors.f90 b/src/neighbors.f90 index 27cd608..98799ee 100644 --- a/src/neighbors.f90 +++ b/src/neighbors.f90 @@ -54,7 +54,6 @@ module neighbors if(r_list(1,i) < box_bd(1)) cycle !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 c(j) = int((r_list(j,i)-box_bd(2*j-1))/(rc_off)) + 1 end do diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 69cb818..5e30420 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -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=2) :: species_type(10) 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 @@ -109,11 +109,11 @@ module opt_group call get_command_argument(arg_pos, type, arglen) if (arglen==0) STOP "Missing select_type in group command" select case(trim(adjustl(type))) - case('atoms', 'elements', 'both') + case('atoms', 'elements', 'all') continue case default 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 arg_pos = arg_pos + 1 @@ -177,6 +177,57 @@ module opt_group end if 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') arg_pos = arg_pos + 1 @@ -239,7 +290,7 @@ module opt_group !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. + !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))) case('atoms') !Read number of ids @@ -285,7 +336,7 @@ module opt_group end do group_ele_num = in_num - case('both') + case('all') !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 @@ -308,7 +359,7 @@ module opt_group 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 ", & + print *, "Atoms specifier used more than once in group id command with type all, either use type ", & "atoms or include elements identifier" stop 3 @@ -326,7 +377,7 @@ module opt_group 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 ",& + print *, "Elements specifier used more than once in group id command with type all, either use type ",& "elements or include atoms identifier" stop 3 @@ -527,7 +578,7 @@ module opt_group allocate(element_index(1024), atom_index(1024)) !Check the type to see whether we need to find the elements within the group select case(trim(adjustl(type))) - case('elements', 'both') + case('elements', 'all') do i = 1, ele_num if(in_group_ele(size_ele(i), lat_ele(i), r_node(:,:,:,i))) then group_ele_num = group_ele_num + 1 @@ -559,7 +610,7 @@ module opt_group end select !Check the type to see if we need to find the atoms within the group select case(trim(adjustl(type))) - case('atoms','both') + case('atoms','all') do i = 1, atom_num if(in_group(r_atom(:,i)).neqv.flip) then group_atom_num = group_atom_num + 1 @@ -1151,10 +1202,8 @@ module opt_group call add_atom_type(species_type(i), type_map(i)) if(i > 1) s_fractions(i) = s_fractions(i) + s_fractions(i-1) end do - print *, s_fractions !Now randomize the atom types added_types = 0 - print *, type_map(j) do i = 1, group_atom_num ia = atom_index(i) call random_number(rand) @@ -1184,6 +1233,9 @@ module opt_group in_group=in_block_bd(r,block_bd) case('wedge') 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') do i = 1, 3 if (.not.((dim1==i).or.(dim2==i))) dim3 = i diff --git a/src/subroutines.f90 b/src/subroutines.f90 index 5cf147f..3fec73f 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -2,6 +2,8 @@ module subroutines use parameters use functions use box + use str + implicit none integer :: allostat, deallostat @@ -149,29 +151,43 @@ module subroutines subroutine parse_ori_vec(ori_string, ori_vec) !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) :: prefactor + integer :: start, fin integer :: i, ori_pos, stat - ori_pos=2 - 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 + do i = 1, len(ori_string) + if(ori_string(i:i) == "[") then + start=i+1 + else if(ori_string(i:i) == "]") then + fin = i-1 end if 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 end subroutine parse_ori_vec