Calculating the sum $\frac{1}{2} \sum x^T \Sigma x$ for all $x \in \{0,1\}^n$ Note: the equation inside the sum is related to Boltzmann Machines / Hopfield Networks, the energy function of these functions are similar. For further info, i.e. on how to derive the maximum likelihood estimator for $\Sigma$ you could take a look at David MacKay's book, section 43 MacKay calls it $W$)
http://www.inference.phy.cam.ac.uk/itprnn/ps/521.526.pdf
Original Question:
I am trying to calculate the sum 
$$ 
\frac{1}{2} \sum_{\forall x} x^T \Sigma x 
 $$
where $\Sigma$ is a symmetric, positive definite matrix, $x$ is a column vector with binary values. I can assume normalization of $\Sigma$ as another step (so that all rows sum to 1), if that will give me convergence. $\Sigma$ is the result of cooccurence calculation $A^TA$ on sample data that contains only binary values, 0 or 1. 
I ran some numerical experiments on randomly generated samples, 
from sklearn.preprocessing import normalize
import numpy as np
np.random.seed(0)
A = np.random.randint(0,2,(100,4))
cooc = A.T.dot(A).astype(float)
cc = normalize(cooc, norm='l1', axis=1) 

When I apply the formula, that is, run through all possible values of $x \in \{0,1\}^n$, and calculate the sum,
import itertools, numpy.linalg as lin
lst = np.array(list(itertools.product([0, 1], repeat=A.shape[1])))
s = 0
for x in lst: 
    s += np.dot(np.dot(x.T,cc), x) / 2
print s

I get 11.15 on data with 4 dimensions. For n=5 I get around 26.6, on any sample I generate with these dimensions. From this I concluded this number is directly tied to the dimension I am working with, and is a constant, so I was wondering if it could be calculated as a limit, somehow. 
Note: 
I need this number as a normalization constant, I plan to use $p(x;\Sigma) = \frac{1}{2C} x^T \Sigma x$ as a probability mass function. Here is how I arrived to this question. 
I was trying to capture frequency of each single variable and dependency between all variables of multivariate binary samples. I used a simple cooccurence calculation on the sample, 
import numpy as np
A = np.array([\
[0.,1.,1.,0],
[1.,1.,0, 0],
[1.,1.,1.,0],
[0, 1.,1.,1.],
[0, 0, 1.,0]
])
c = A.T.dot(A).astype(float)
print c 

Result
[[ 2.  2.  1.  0.]
 [ 2.  4.  3.  1.]
 [ 1.  3.  4.  1.]
 [ 0.  1.  1.  1.]]

Now for any new data point $x$, if I wanted to calculate a "score", ex:
x = np.array([[0,1,1,0]])
print np.dot(np.dot(x.T,c), x) / 2

would give me 7. The formula basically picks numbers 4,3,3,4 in the middle block of the cooc matrix, and sums them, which was what I wanted because new data point $x=[0,1,1,0]$ has binary variables 2 and 3 turned on (is 1) so I am interested in the dependency between $x_2$ and $x_3$, as well as the frequency of the variables by themselves.
Once I had this score, I started wondering if I could turn this function into a PMF, hence the need for normalization constant and the need to integrate the function for all possible values of $x$.
I toyed with the idea of dividing the sum by $x^Tx$, thereby causing the equation to look like the Rayleigh Quotient, 
$$ 
= \frac{1}{2} \sum_{\forall x} \frac{x^T \Sigma x }{x^Tx}
 $$
then, if I assumed "x=all eigenvalues" instead of "x=all possible values" then perhaps summing all eigenvalues would give me something. But the summation must be for all x. Representing all x's using eigenvectors as basis maybe.. ? 
 A: Let $X_{1}, \cdots, X_{n}$ be i.i.d. Bernoulli random variables. Then with $\Sigma = (\sigma_{ij})$,
\begin{align*}
\frac{1}{2} \sum_{\mathrm{x} \in \{0, 1\}^{n}} \mathrm{x}^{T} \Sigma \mathrm{x}
&= 2^{n-1}\Bbb{E} \sum_{i,j} \sigma_{ij}X_{i}X_{j} \\
&= 2^{n-1}\sum_{i,j} \sigma_{ij} \Bbb{E} (X_{i}X_{j}) \\
&= 2^{n-1}\bigg( \sum_{i} \sigma_{ii} \Bbb{E} (X_{i}) + \sum_{i \neq j} \sigma_{ij} \Bbb{E} (X_{i}X_{j}) \bigg) \\
&= 2^{n-3}\bigg( \operatorname{tr} \Sigma + \sum_{i, j} \sigma_{ij} \bigg).
\end{align*}
Assuming the normalization step so that $\Sigma$ is right stochastic matrix, then the sum $\sum_{i, j} \sigma_{ij}$ reduces to $n$ and hence
$$ \frac{1}{2} \sum_{\mathrm{x} \in \{0, 1\}^{n}} \mathrm{x}^{T} \Sigma \mathrm{x} = 2^{n-3}(n + \operatorname{tr}\Sigma). $$
The expectation of this quantity over all the possible $\Sigma$ depends on the distribution of $A$, but we always have $1 \leq \operatorname{tr}\Sigma \leq n$, so the value $(n+1) 2^{n-3}$ approximates this value well up to magnitude.
A: Based on @sos404's answer, the final equation is
$$
p(x;\Sigma,n) = \frac{ \mathrm{x}^{T} \Sigma \mathrm{x}}{2^{n-2}(n + tr \ \Sigma)} 
$$
Some numeric checks
from sklearn.preprocessing import normalize
import itertools
n = 4
A = np.random.randint(0,2,(200,n))
c = A.T.dot(A).astype(float)
cc = normalize(c, norm='l1', axis=1) 
lst = np.array(list(itertools.product([0, 1], repeat=n)))
res = [np.dot(np.dot(x.T,cc), x) / \
         ((2**(n-2)) * (n+np.sum(np.diag(cc)))) for x in lst]
print np.array(res).sum()

No matter what value is chosed for D, the result is always 1.0.
The itertools.product is the function call that allows the iteration over all possible values of $x \in \{0,1\}^n$. 
