It is an interesting combinitorial question since nchoosek(120,56) gives 7.4149e+34 possible ways to select 56 pairs out of the 120. So you shouldn't solve this by a loop over all selections.
To handle the problem I introduce a representation of selected pairs given by a matrix A = (aij), where aij == 1 if and only if (i,j) oder (j,i) is a pair among the 56 selected ones. This matrix represents a selection of 56 pairs as described above if and only if row sums are 7, i.e. A*ones(16,1) = 7*ones(16,1) (symmetric!). It is easy to solve this systematically to print all possible solutions, but this results in about 1e20 candidates...
But according to this high number you can easily select one example randomly by some stupid code like the following:
%% find some random example
A = zeros(16,16);
b1 = ones(16,1);
b7 = 7*b1;
while norm(A*b1 - b7) ~= 0
    A = zeros(16,16);
    change = true;
    while norm(A*b1 - b7) ~= 0 && change
        change = false;
        todo = find(A*b1<7);
        idx = todo(randperm(length(todo)));
        for k = 1:2:length(idx)-1
            if A(idx(k),idx(k+1)) == 0
                change = true;
            end
            A(idx(k),idx(k+1)) = 1;
            A(idx(k+1),idx(k)) = 1;
        end        
    end    
end
%% print it
pIdx = 0;
for lIdx = 1:16
    for cIdx = 1:lIdx-1
        if A(lIdx,cIdx) == 1
            fprintf(['(' num2str(cIdx) ',' num2str(lIdx) ')  ']);
            pIdx = pIdx + 1;
            if mod(pIdx,7) == 0
                fprintf('\n');
            end            
        end
    end
end
Example result:
(1,3)  (2,4)  (1,5)  (2,5)  (2,6)  (3,6)  (4,6)  
(5,6)  (1,7)  (3,7)  (6,7)  (4,8)  (5,8)  (2,9)  
(5,9)  (6,9)  (7,9)  (8,9)  (2,10)  (3,10)  (4,10)  
(5,10)  (7,10)  (1,11)  (3,11)  (6,11)  (7,11)  (9,11)  
(1,12)  (4,12)  (7,12)  (11,12)  (3,13)  (8,13)  (9,13)  
(10,13)  (1,14)  (2,14)  (4,14)  (8,14)  (12,14)  (13,14)  
(1,15)  (3,15)  (5,15)  (8,15)  (11,15)  (12,15)  (13,15)  
(2,16)  (4,16)  (8,16)  (10,16)  (12,16)  (13,16)  (14,16)
Update: Added code to scan all possible solutions...
function main()
    %% provide n-choose-k sequences (compute only once!)
    global ncks maxResult;
    disp('Prepare n choose k sequences...');
    for n = 1:15
        fprintf([num2str(n) '...  ']);
        for k = 1:min(7,n)
            ncks{n}{k} = nchoosek_sequence(n,k);
        end
    end
    fprintf('\n');
    %% start finding solutions
    disp('Scan space of solutions...')
    A = zeros(16,16);
    maxResult = -Inf;
    tic;
    findMats(A,1)
    toc;
end
function findMats(A,curRow)
    global ncks maxResult;
    remain = 7-sum(A+A');
    if curRow == 16
        if remain(16) == 0
            A = A + A';
            if norm(sum(A)-7*ones(1,16)) ~= 0
                error('AHHH');
            end
            %A
            %pause();
            maxResult = max(YOUR_FUNCTION(A),maxResult); % <- add your function here
        end
        return;
    end
    if remain(curRow) == 0
        findMats(A,curRow+1);
    elseif remain(curRow) > 0
        remainingRows = curRow + find(remain(curRow+1:end) > 0);
        if ~isempty(remainingRows) && length(remainingRows) >= remain(curRow)
            for r = 1:nchoosek(length(remainingRows),remain(curRow))
                row = remainingRows(ncks{length(remainingRows)}{remain(curRow)}(r,:));
                Anew = A;
                Anew(curRow,row) = 1;
                findMats(Anew,curRow+1);
            end
        end
    end
end
% very simple and not optimized way to get a sequence of subset-selections
% in increasing order
%
% Note: This function was easy to implement and sufficient for our
% purposes. For more efficient implementations follow e.g.
%  http://stackoverflow.com/questions/127704/algorithm-to-return-all-combin
%  ations-of-k-elements-from-n
function s = nchoosek_sequence(N,k)
    s = zeros(nchoosek(N,k),k);
    idx = 1;
    for n=(2^k-1):(2^N-1)
        if sum(dec2bin(n)=='1') == k
            temp = dec2bin(n);
            s(idx,:) = find([zeros(1,N-length(temp)), temp=='1']);
            idx = idx + 1;
        end
    end    
    % just for convinience ;-)
    s = s(end:-1:1,:);
end
But this can take a while... ;-)