diff --git a/src/opt_group.f90 b/src/opt_group.f90 index cfd6499..8ff789a 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -263,7 +263,7 @@ module opt_group integer :: i, j, ix, iy, iz, inod, ibasis, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num, & current_esize, dof, max_lat(3), r_lat(3), ele(3,8), vlat(3), bd_in_lat(6), bd_in_array(3) - real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), ori_inv(3,3), & + real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), ori_inv(3,3), r(3), & r_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8) logical, allocatable :: lat_points(:,:,:) character(len=100) :: remesh_ele_type @@ -295,8 +295,8 @@ module opt_group group_in_lat = matmul(fcc_inv, matmul(ori_inv, group_in_lat/remesh_lat_pam)) do i = 1, 3 - bd_in_lat(2*i-1) = minval(group_in_lat(i,:)) - bd_in_lat(2*i) = maxval(group_in_lat(i,:)) + bd_in_lat(2*i-1) = nint(minval(group_in_lat(i,:))) + bd_in_lat(2*i) = nint(maxval(group_in_lat(i,:))) end do end select @@ -306,8 +306,11 @@ module opt_group !Now place all group atoms and group interpolated atoms into lat_points do i = 1, group_atom_num - r_lat = r_atom(:,atom_index(i))/remesh_lat_pam - r_lat = matmul(fcc_inv,matmul(ori_inv,r_lat)) + r = r_atom(:,atom_index(i))/remesh_lat_pam + r = matmul(fcc_inv,matmul(ori_inv,r)) + do j = 1, 3 + r_lat(j) = nint(r(j)) + end do !Do a check to make sure the code is working and that lattice points aren't being written on top of each other. !This is primarily a debugging statement if((r_lat(1)==13).and.(r_lat(2)==13).and.(r_lat(3)==-13)) then @@ -349,10 +352,20 @@ module opt_group bd_in_array(2) = bd_in_lat(4) - bd_in_lat(3) + 10 bd_in_array(3) = bd_in_lat(6) - bd_in_lat(5) + 10 - do iz = 1, bd_in_array(3) - do iy = 1, bd_in_array(2) - do ix = 1, bd_in_array(1) + zloop: do iz = 1, bd_in_array(3) + yloop: do iy = 1, bd_in_array(2) + xloop: do ix = 1, bd_in_array(1) if (lat_points(ix, iy,iz)) then + + !Check to see if the element overshoots the bound + if (iz+remesh_size-1 > bd_in_array(3)) then + exit zloop + else if (iy+remesh_size-1 > bd_in_array(2)) then + cycle zloop + else if (ix+remesh_size-1 > bd_in_array(1)) then + cycle yloop + end if + if (all(lat_points(ix:ix+remesh_size-1,iy:iy+remesh_size-1,iz:iz+remesh_size-1))) then do inod = 1, 8 vlat = ele(:,inod) + (/ix, iy, iz /) @@ -369,11 +382,25 @@ module opt_group end if end if - end do - end do - end do - + end do xloop + end do yloop + end do zloop + !Now we have to add any leftover lattice points as atoms + do iz = 1, bd_in_array(3) + do iy=1, bd_in_array(2) + do ix = 1, bd_in_array(1) + if(lat_points(ix,iy,iz)) then + vlat = (/ ix, iy, iz /) + do i = 1, 3 + vlat = vlat + bd_in_lat(2*i-1)-5 + end do + r = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam + call add_atom(remesh_type, sbox_atom(atom_index(1)), r) + end if + end do + end do + end do end subroutine remesh_group