Say we have $n$ initial pairs of people where the people within a pair are identical (i.e., identical twins) but people in different pairs are distinguishable. What is the number of ways to form $k\leq n$ pairs such that none of these new pairs consists of twin pairs?
Motivation/Context
Problem is related to the statistical physics of the dimerization of biomolecules, in particular a collection of homodimer proteins which are identical within the dimer but distinguishable from proteins in other dimers.
Related Problems
I have been able to solve two related problem using the principle of inclusion-exclusion. Here are the statements and answers for reference:
Problem I: Say we have $n$ initial man-woman couples where all the men and all the women are distinguishable. What is the number of ways to form $k\leq n$ man-woman couples such that none of these new couples coincides with any initial pairings?
Answer: The number of ways to form $k\leq n$ man-woman couples such that none is from the initial set of couples is \begin{equation} b_{n, k} = \sum_{m=0}^{k} (-1)^{m} \binom{n}{m} \binom{n-m}{k-m}^2 (k-m)!. \end{equation}
Problem II: Say we have $n$ initial couples where the all $2n$ people are distinguishable. Assuming each person can couple with any other person, what is the number of ways to form $k$ couples such that none of these new couples coincides with any initial pairings?
Answer: The number of ways to form $k\leq n$ the man-woman couples such that none is from the initial set of couples is \begin{equation} a_{n, k} = \sum_{j=0}^{k}(-1)^{j} \binom{n}{j} \binom{2n-2j}{2k-2j} (2k-2j-1)!!. \label{eq:anl_fin} \end{equation}
Let's call the answer to your original problem $c_{n,k}$. We can compute $c_{n,k}$ by dividing $a_{n,k}$ out by the symmetries given by swapping within each initial pair. More precisely, in the context of Problem II, let's say the initial pairs of people were $\{x_1,y_1\},\dots,\{x_n,y_n\}$. Problem II then counts the set $S_{n,k}$ of all sets of $k$ disjoint two-element subsets of $\{x_1,y_1,\dots,x_n,y_n\}$, such that none of the subsets are among the original pairs. Your original problem then counts the quotient of $S_{n,k}$ by the action of the group $G=(\mathbb{Z}/2\mathbb{Z})^n$ where the $j$th coordinate says whether to swap $x_j$ and $y_j$.
So, we just have to figure out how large the stabilizers of this action are. If $P\in S_{n,k}$, let its width be the number of different $j$ such that either $x_j$ or $y_j$ is in some element of $P$. Let $S_{n,k,w}$ be the set of elements of $S_{n,k}$ of width $w$ and let $a_{n,k,w}=|S_{n,k,w}|$. Note that the stabilizer any $P\in S_{n,k,w}$ has at least $2^{n-w}$ elements: for each of the $n-w$ pairs $\{x_j,y_j\}$ are disjoint from all elements of $P$, swapping $x_j$ and $y_j$ will not change $P$. However, the stabilizer could be larger: there are distinct $i,j$ such that both $\{x_i,x_j\}$ and $\{y_i,y_j\}$ (or both $\{x_i,y_j\}$ and $\{y_i,x_j\}$) are in $P$, then simultaneously swapping $x_i$ with $y_i$ and $x_j$ with $y_j$ will not change $P$, increasing the size of the stabilizer by a factor of $2$.
Let us call such a pair $\{i,j\}$ a bad pair, and let $S_{n,k,w,\ell}$ denote the set of those elements of $S_{n,k,w}$ that have exactly $\ell$ bad pairs. Then the stabilizer of the action of $G$ for each $P\in S_{n,k,w,\ell}$ has $2^{n-w+\ell}$ elements, so each orbit of $G$ on $S_{n,k,w,\ell}$ has $2^{w-\ell}$ elements. Thus we have $$c_{n,k}=\sum_{\ell=0}^{\lfloor k/2\rfloor}\sum_{w=k}^{\min(n,2k-2\ell)}\frac{a_{n,k,w,\ell}}{2^{w-\ell}}$$ where $a_{n,k,w,\ell}=|S_{n,k,w,\ell}|$.
Of course, for this formula to be of any use, we still need to compute $a_{n,k,w,\ell}$. First, note that the number of ways to choose and rearrange $\ell$ bad pairs is $$\frac{n!}{(n-2\ell)!\ell!}.$$ This is because there are $\frac{n!}{(n-2\ell)!}$ ways to choose $2\ell$ ordered elements of $\{1,\dots,n\}$ to form the pairs, then you divide by $2^\ell$ since the order within each pair does not matter, and divide by $\ell!$ since the order of the $\ell$ pairs does not matter. However, each bad pair can be rearranged in two different configurations ($\{x_i,x_j\}$ and $\{y_i,y_j\}$ or $\{x_i,y_j\}$ and $\{y_i,x_j\}$), so we multiply by $2^\ell$.
Now, given a specific rearrangement of $\ell$ pairs to be bad as above, the number of elements of $S_{n,k,w}$ for which exactly those pairs are bad is $a_{n-2\ell,k-2\ell,w-2\ell,0}$. We can then count $a_{n-2\ell,k-2\ell,w-2\ell,0}$ by starting from $a_{n-2\ell,k-2\ell,w-2\ell}$ and then using inclusion-exclusion to exclude all cases with any bad pairs, similar to how you solved your Problems I and II. This gives $$\begin{align*} a_{n,k,\ell,w}&=\frac{n!}{(n-2\ell)!\ell!}\sum_{m=0}^{\lfloor k/2\rfloor-\ell}(-1)^m\frac{(n-2\ell)!}{(n-2\ell-2m)!m!}a_{n-2\ell-2m,k-2\ell-2m,w-2\ell-2m} \\ &=\frac{n!}{\ell!}\sum_{m=0}^{\lfloor k/2\rfloor-\ell}(-1)^m\frac{a_{n-2\ell-2m,k-2\ell-2m,w-2\ell-2m}}{(n-2\ell-2m)!m!} \end{align*}$$
Now we just have to compute $a_{n,k,w}$. We can do this by an inclusion-exclusion process similar to how you computed $a_{n,k}$. First, note that the number of ways to pick $2k$ elements of $\{x_1,y_1,\dots,x_n,y_n\}$ which intersect $w$ of the original pairs is $$2^{2w-2k}\binom{n}{2w-2k}\binom{n-2w+2k}{2k-w}.$$ Here we are first picking $2w-2k$ elements of distinct original pairs, and then $2k-w$ more full pairs, to give a total of $(2w-2k)+2(2k-w)=2k$ elements spanning $(2w-2k)+(2k-w)=w$ different pairs. Having picked these elements, there are then $(2k-1)!!$ ways to arrange them into new pairs, without regard for whether all the new pairs are mismatched. Using inclusion-exclusion to remove the cases where some new pairs are not mismatched, we get $$a_{n,k,w}=\sum_{j=0}^{2k-w}(-1)^j\binom{n}{j}2^{2w-2k}\binom{n-j}{2w-2k}\binom{n-2w+2k-j}{2k-w-j}(2k-2j-1)!!.$$
Combining all these formulas gives the following rather horrendous formula for $c_{n,k}$ that does not fit on one line:
$$c_{n,k}=\sum_{\ell=0}^{\lfloor k/2\rfloor}\sum_{w=k}^{\min(n,2k-2\ell)}\sum_{m=0}^{\lfloor k/2\rfloor-\ell}\sum_{j=0}^{2k-w-2\ell-2m} \\ (-1)^{m+j} 2^{w-2k+\ell} \frac{n!(2k-4\ell-4m-2j-1)!!}{\ell!m!j!(2w-2k)!(2k-w-2\ell-2m-j)!(n-w)!}$$