Arbitrary Phase FIR design using FFT
You are trying to design a FIR allpass filter. That's impossible (other than a simple delay) since true allpass filters consist of zero/pole pairs that are inverse of each others and the FIR poles are always at the origin.

However, with enough resources, you can always get "close enough" for your specific application. Of course, you would have to define what "close enough" means. 

The actual design is simple enough


Sample your filter target on a very dense FFT grid
Do an inverse FFT
Window/truncate the resulting impulse response until you meet your "close enough" requirements. Typically the impulse response will decay exponentially, so truncation often gives better results that any type windowing (for the same filter length)  
In most cases the impulse response will be non-causal, so you will need to add bulk delay.


Below is a piece of Matlab code that does that

%% start with a Butter filter
fs = 44100; % sample rate
bwCut = 1000; % 1000 Hz cut off frequency
bwOrder = 3;  % thrid order
%% design the filter
[b,a] = butter(bwOrder,bwCut/(fs/2));

%% sample on a large FFT grid
nFFT = 16384;
delta = zeros(nFFT,1); delta(1) = 1;
bwIR = filter(b,a,delta); % impulse response
bwTF = fft(bwIR);         % transer fundction

%% design the target
targetTF = exp(-1i.*angle(bwTF)); % that works every where except for Nyquist, where
                                  % the phase is undefined, we'll manually smootht this out
targetTF(nFFT/2+1) = 1;           % Nyquist phase repair                                  
h0 = real(ifft(targetTF));        % target impulse response

%% now we have to cut it to a practical size
% let's see if we can capture all samples that are above a certain
% thershold. That's in essence a rectangular window
thresh = max(abs(h0)).*0.01;      % -40dB re max
tmp = circshift(h0,nFFT/2);       % shift to the middle of FFT window
i1 = find(abs(tmp) &gt; thresh);     % find the area above threshold
h1 = tmp(min(i1):max(i1));        % grab that area

%% analyse the result. This really depends on the application.
% In this example we will look at the difference of the resulting
% transfer function to the "ideal" answer. The difference looks at both the
% real and imaginary parts, so it's a decent metric that combines magnitude
% and phase. It coould be interpreted as the amount of "noise" induced by
% the truncation error.
idealTF = abs(bwTF);

% model the results, rectangular window
h2 = conv(h1,bwIR);                   % convolution
h2 = h2(1:nFFT);                      % cut off the excess length
h2 = circshift(h2,min(i1)-nFFT/2-1);  % time align
diffTFRect = fft(h2)-idealTF;         % calculate difference
% same length hanning window
h1w = h1.*hanning(length(h1));
h2 = conv(h1w,bwIR);                  % convolution
h2 = h2(1:nFFT);                      % cut off the excess length
h2 = circshift(h2,min(i1)-nFFT/2-1);  % time align
diffTFHanning = fft(h2)-idealTF;


%% plot it
clf;
i = (1:nFFT/2);
fx = [diffTFRect diffTFHanning];
plot(i*fs/nFFT, 20*log10(abs(fx(1+i,:))));
xlabel('Frequency in Hz -&gt;');
set(gca,'xlim',fs*[1/nFFT 1/2]);
set(gca,'xscale','log');
set(gca,'ylim',[-120, 10]);
ylabel('Difference in dB');
grid('on');
legend('Rectangle','Hanning');
title('Error');