diff --git a/src/opt_group.f90 b/src/opt_group.f90 index 2176a88..3b40f14 100644 --- a/src/opt_group.f90 +++ b/src/opt_group.f90 @@ -261,7 +261,7 @@ module opt_group subroutine remesh_group !This command is used to remesh the group to a desired element size - integer :: i, j, ix, iy, iz, inod, ibasis, ie, type_interp(max_basisnum*max_esize**3), add_atom_num, orig_atom_num, & + integer :: i, j, k, 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), r(3), & r_new_node(3,max_basisnum, max_ng_node), orient(3,3), group_in_lat(3,8) @@ -313,9 +313,6 @@ module opt_group 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 - print *, i, atom_index(i), r_atom(:,atom_index(i)) - end if if(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then stop "Multiple atoms share same position in lat point array, this shouldn't happen" else @@ -328,7 +325,11 @@ module opt_group ie = element_index(i) call interpolate_atoms(type_ele(ie), size_ele(ie), lat_ele(ie), r_node(:,:,:,ie), type_interp, r_interp) do j = 1, size_ele(ie)**3 * basisnum(lat_ele(ie)) - r_lat = matmul(fcc_inv,matmul(ori_inv, r_interp(:,atom_index(i))/remesh_lat_pam)) + 1 + r = r_interp(:,j)/remesh_lat_pam + r = matmul(fcc_inv,matmul(ori_inv,r)) + do k = 1, 3 + r_lat(k) = nint(r(k)) + 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(lat_points(r_lat(1)-bd_in_lat(1)+5,r_lat(2)-bd_in_lat(3)+5,r_lat(3)-bd_in_lat(5)+5)) then @@ -395,6 +396,7 @@ module opt_group do i = 1, 3 vlat(i) = vlat(i) + bd_in_lat(2*i-1)-5 end do + lat_points(ix,iy,iz) = .false. r = matmul(orient, matmul(fcc_mat, vlat))*remesh_lat_pam call add_atom(remesh_type, sbox_atom(atom_index(1)), r) end if