Let $A$ be $d\times d$ matrix, $b\in\mathbb R^d$, $X$ a $d\times n$ matrix, $W$ an $n\times n$ positive definite diagonal matrix with the diagonal elements summing to one, $1_n$ a vector of $n$ ones and $a,c\in\mathbb R$. I need to solve the following system of equations for $c$:
$$\begin{align*} a &= \int \exp(-0.5x'Ax + b'x + c)\mathrm dx\\ aXW1_n &=\int x\exp(-0.5x'Ax + b'x + c)\mathrm dx \\ aXWX' &=\int xx'\exp(-0.5x'Ax + b'x + c)\mathrm dx \end{align*}$$ Using the rules for Gaussian integration, I obtain for the right hand sides: $$\int x\exp(-0.5x'Ax + b'x + c)\mathrm dx = \exp(c)\sqrt{\frac{(2\pi)^d}{\det A}}\exp(0.5b'A^{-1}b)A^{-1}b$$ and $$\int xx'\exp(-0.5x'Ax + b'x + c)\mathrm dx = \exp(c)\sqrt{\frac{(2\pi)^d}{\det A}}\exp(0.5b'A^{-1}b)(A^{-1}bb'A^{-1} + A^{-1})$$ where I assumed that $A$ is symmetric and psotive definite. Otherwise I would not be able to perform the change of variables and derive the expressions. Dividing the second and third equation by the first -- note that $$ \int \exp(-0.5x'Ax + b'x + c) = \exp(c)\sqrt{\frac{(2\pi)^d}{\det A}}\exp(0.5b'A^{-1}b)$$ -- I get: $XW1_n = A^{-1}b$ and $XWX' = A^{-1}bb'A^{-1} + A^{-1}$. Thus, $XWX' - XW1_n1_n'WX' = A^{-1}.$ But this matrix is certainly not positive definite since $1_n'W1_n - 1_n'W1_n1_n'W1_n = 0$ (as the diagonal elements sum to one). So how can this happen? What causes this phenomen? Apparently this matrix is $A$ is indefinite as it also posses negative eigenvalues.
I checked the calculations multiple times and could not spot an error.
Your calculations are correct. But the matrix $$M:=XWX' - XW1_n1_n'WX'$$ is indeed positive definite. This is not affected by the fact that $1_n'W1_n=1$.
To show $v'Mv\ge 0$ for every $v$, write $y:=X'v$, and write the matrix $W$ as $\operatorname{diag}(w_1,\ldots,w_n)$ with all $w_i>0$, so that $W1_n=(w_1,\ldots,w_n)'$ and $y'W1_n=\sum w_iy_i$. Then the expression $v'Mv$ becomes $$y'Wy-y'W1_n(y'W1_n)'=\sum w_iy_i^2-\left(\sum w_iy_i\right)^2,$$ which is non-negative by Cauchy-Schwarz: $$ \left(\sum\sqrt{w_i}\sqrt{w_i}y_i\right)^2\le\sum(\sqrt{w_i})^2\sum(\sqrt{w_i}y_i)^2 $$ or by an application of Jensen's inequality (the function $t\mapsto t^2$ is convex).