Background
Let us recollect the Goldbach's conjecture. For any $n \geq 2$ there exists primes $p_i$ and $p_j$ such that:
$$ 2n = p_i + p_j $$
where $p_k$ is the $k$'th prime. Now we define a function $G(2n)$ as the number of "unique" solutions for $2n$. For example:
$$ G(4) = 2 + 2 = p_1 + p_1$$
Since $p_i = p_j$ the number of "unique" solutions is $1$. However for $G(8)$:
$$ G(8) = 3 + 5 = p_2 + p_3$$
Since $p_i \neq p_j$ the number of solutions is $2$. What about sometime like $G(14)$?
$$G(14) = 7+ 7 = 3+ 11 $$
Hence the number of unique solutions is $G(14)=1 + 2 = 3$
Essentially for each prime pairs $p_i$ and $p_j$ whose sum is $2n$ if we add all such pairs (twice if $i \neq j$ and once if $i=j$). If there are no solutions for the number $k$ then $G(k) = 0$ (for example $G(2) = 0$).
Question
Is the below conjecture true or false?
$$ \frac{P(s)^2}{\zeta(s)} \sim 2 \sum_{r=3}^\infty \frac{G(2r)}{(2r)^{2s}} $$
As $s \nearrow 1$. Or:
$$ \lim_{s \to 1} \frac{\frac{P(s)^2}{\zeta(s)} }{ \sum_{r=3}^\infty \frac{G(2r)}{(2r)^{2s}}} = 2$$
Heuristic Reasoning
Let's think of the prime zeta function as Dirchlet series:
$$ P(s) = 2^{-s} + 3^{-s} + \dots $$
Choosing a function $f(x) = e^{-x}$ and this formula:
$$ (e^{-2\epsilon} + e^{-3 \epsilon} + e^{-5 \epsilon} + \dots) \epsilon \sim \frac{P(s)}{\zeta(s)} \times \int_0^\infty f(x) dx $$
as $s \nearrow 1$ and $\epsilon$ tends to $0$. Removing $e^{-2 \epsilon}$ since it is small
$$ ( e^{-3 \epsilon} + e^{-5 \epsilon} + \dots) \epsilon \sim \frac{P(s)}{\zeta(s)} $$
since $\int_0^\infty f(x) dx = 1$. Squaring both sides and resuming the L.H.S:
$$ ( G(6)e^{-6 \epsilon} + G(8)e^{-8 \epsilon} + \dots) \epsilon^2 \sim (\frac{P(s)}{\zeta(s)})^2 $$
However recognizing the L.H.S as an integral of a function $g(x) = e^{- \sqrt x}$ in it's own right with $\epsilon^2 \to \epsilon' $:
$$ ( G(6)e^{-6 \sqrt \epsilon '} + G(8)e^{-8 \sqrt \epsilon'} + \dots) \epsilon' \sim \sum_{r=3}^\infty \frac{G(2r)}{(2r)^{2s}} \times \frac{1}{\zeta(s)} \times \int_0^\infty g(x) dx $$
Comparing the above $2$ equations:
$$ \frac{P(s)^2}{\zeta(s)} \sim 2 \sum_{r=3}^\infty \frac{G(2r)}{(2r)^{2s}} $$
since $\int_0^\infty g(x) dx = 2$
It is clearly false as $(\frac{P(s)}{\zeta(s)})^2 \to 0$ as $s\to 1^+$ in contrary to your Dirichlet series with non-negative coefficients.
It follows from the PNT that as $x\to 0^+$
$$f(x)=\sum_{p\ge 3} e^{-px}=(1-e^{-x}) \sum_{n\ge 2} (\pi(n)-1)e^{-nx}$$ $$\sim (1-e^{-x}) \sum_{n\ge 2} \sum_{2\le m\le n} \frac1{\log m} e^{-nx}= \sum_{n\ge 2} \frac{e^{-nx}}{\log n} $$ For $s > 2$ $$\sum_{r \ge 1} G(2r)(2r)^{-2s}=\frac1{\Gamma(s)}\int_0^\infty \sum_{r\ge 1} G(2r) e^{-2rx}x^{s-1}dx=\frac1{\Gamma(s)}\int_0^\infty f(x)^2 x^{s-1}dx$$ From this, since $f(x)=O(\frac{e^{-x}}{1-e^{-x}})$, you get that for some $0<A<B$ $$\sum_{r \ge 1} G(2r)(2r)^{-2s}\in (A,B)\int_0^{1/2} f(x)^2 x^{s-1}dx$$
As $x\to 0$, $$f(x)= O(\frac{1}{x\log x}),\qquad \frac{1}{x\log x}=O(f(x))$$ whence you get that for some $0<a<b$ $$\sum_{r \ge 1} G(2r)(2r)^{-2s} \in (a,b)\int_0^{1/2} \frac1{x^2 \log x^2}x^{s-1}dx$$ which converges for $s\ge 2$. The RHS diverges for $s<2$ and so does the LHS.