It sounds like you're trying to find the shifts in the energy levels caused by $H_V$. To find the energy shift that goes as the first order in your small parameter you compute
$$E_n^{(1)} = \langle \psi_n^{(0)}|H_V|\psi_n^{(0)} \rangle \quad (1)$$
as you noted. In this expression $|\psi_n^{(0)}\rangle$ means "the $n^{\textrm{th}}$ eigenstate of the unperturbed Hamiltonian." That superscript (0) means "to zeroth order in the perturbation" ie. without any perturbation.
So to compute the energy level shifts you have to know the eigenstates of the unperturbed Hamiltonian. You say you already found those to be $|1,\pm 1\rangle$, $|1,0\rangle$ and $|0,0\rangle$. To make contact with the notation in equation (1) let's go ahead and number these
$$|\psi_1^{(0)}\rangle = |0,0 \rangle \quad |\psi_2^{(0)}\rangle = |1,0 \rangle \quad |\psi_3^{(0)}\rangle = |1,-1 \rangle \quad |\psi_4^{(0)}\rangle = |1,+1 \rangle$$
Now you want to compute the energy shifts. Let's just do one of them as an example. Let's compute $E_1^{(1)}$.
$$\begin{eqnarray}
E_1^{(1)} &=& \langle \psi_1^{(0)}|H_V|\psi_1^{(0)}\rangle \\
&=& -\mu B \cdot \langle0,0|\sigma_1 + \sigma_2 | 0,0 \rangle \\
&=& -\mu B \cdot \left( \langle0,0|\sigma_1|0,0\rangle + \langle0,0|\sigma_2|0,0\rangle \right) \end{eqnarray}$$
Now, those two operators $\sigma_1$ and $\sigma_2$ are the spin operators for each individual particle. The state $|0,0\rangle$ is expressed in a basis that is not the single particle basis. $|0,0\rangle$ is a state designed to be simple in the basis of total spin. In order to work out the matrix elements you need for the computation you have to either convert the operators to the total spin basis, or re-express the states in the individual spin basis. I think the latter is easier in this case. Just use
$$|0,0\rangle \equiv \left( \uparrow \downarrow - \downarrow \uparrow \right)/\sqrt{2}.$$
Now it's easy to compute what you need. Here's an example
$$B \cdot \langle 0,0|\sigma_1| 0,0\rangle = \frac{1}{2} \langle\uparrow \downarrow-\downarrow \uparrow|B^x\sigma^x_1 + B^y\sigma^y_1 + B^z\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle.$$
You asked about what the matrices for $\sigma$ should be. Now that we've got everything expressed in the single particle spin basis, it's easy. For example,
$$\sigma^x_1 = \left(\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array} \right).$$
The thing to note is that this only acts on the first particle, so you get things like this
$$\begin{eqnarray}
\langle \uparrow \downarrow | \sigma^x_1 | \downarrow \uparrow \rangle &=& \langle \uparrow | \sigma^x_1 | \downarrow \rangle\cdot\langle \downarrow|\uparrow \rangle \\
&=& \langle \uparrow|\uparrow \rangle \cdot \langle \downarrow|\uparrow \rangle \\
&=& 1 \cdot 0 \\
&=& 0.
\end{eqnarray}$$
At this point I'm sure you can figure out what's going on. I think your problem was one of notation and hopefully seeing this written out helps.
Extension:
I've shown how to compute simple terms like
$$\langle \uparrow \downarrow | \sigma^x_1 | \downarrow \uparrow \rangle$$
but the original expression to compute was more complex. In particular we had
$$\langle\uparrow \downarrow-\downarrow \uparrow|B^x\sigma^x_1 + B^y\sigma^y_1 + B^z\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle.$$
The trick here is that this is all linear so we can break it into little parts like this
$$\langle\uparrow \downarrow-\downarrow \uparrow|B^x\sigma^x_1 + B^y\sigma^y_1 + B^z\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle$$
$$=B_x\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^x_1|\uparrow\downarrow - \downarrow\uparrow\rangle + B_y\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^y_1|\uparrow\downarrow - \downarrow\uparrow\rangle + B_z\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle.$$
Each one of these terms can be broken down further:
$$\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^x_1|\uparrow\downarrow - \downarrow\uparrow\rangle$$
$$\langle \uparrow \downarrow|\sigma^x_1|\uparrow\downarrow\rangle - \langle \uparrow\downarrow | \sigma^x_1 | \downarrow \uparrow\rangle - \langle \downarrow \uparrow | \sigma^x_1 | \uparrow \downarrow \rangle + \langle \downarrow \uparrow | \sigma_1^x | \downarrow \uparrow \rangle$$
From here you can surely see how to to calculate the desired quantity. Note that if you think carefully you can tell that many terms are zero without having to calculate them.
Best Answer
I'm not so sure, if this is really, what you're looking for, but you can of course solve this easy problem analytically.
To do this, it is clever to first analyze the easier Hamiltonian $H_0 = 2g (\vec L \cdot \vec S)$, where the $L_i$ and $S_j$ fulfill independent $SU(2)$-algebrae $$ [L_i, L_j] = i \epsilon_{ijk} L_k\\ [S_i, S_j] = i \epsilon_{ijk} S_k. $$ This Hamiltonian can be written as $$H_0 = g(J^2 - L^2 - S^2),$$ where we have defined the following operators: $$L^2 = \sum_{i=1}^3 L_i^2 \otimes \mathbb{1},\\ S^2 = \sum_{i=1}^3 \mathbb{1} \otimes S_i^2, \\ J_i = L_i\otimes \mathbb{1} + \mathbb{1}\otimes S_i,\\J^2 = \sum_{i=1}^3 J_i^2.$$ Now as $L^2$ and $S^2$ are equal to $\frac{1}{2} (1+\frac{1}{2})$ on the subspace we are interested in (namely the one of a spin-1/2-particle), we can write the Hamiltonian as $$H_0 = g(J^2 - 3/2)$$ and by simple addition of angular momenta, we find the eigenstates $|j, m\rangle$: $$|1, 1\rangle = \left|\uparrow\uparrow\right\rangle\\ |1, 0\rangle = \frac{1}{\sqrt{2}}(\left|\uparrow\downarrow\right\rangle+\left|\downarrow\uparrow\right\rangle) \\ |1, -1\rangle = \left|\downarrow\downarrow\right\rangle\\ |0, 0\rangle = \frac{1}{\sqrt{2}}(\left|\uparrow\downarrow\right\rangle-\left|\downarrow\uparrow\right\rangle)$$ with energies $E_0 = -3g/2$ and $E_1 = g/2$ respectively (and obvious notation for the product base of the two particle hilbert space).
Now let's proceed to the real problem and add the second term. As you said, $$[J^2, L_z] = 2i (\vec L\times \vec S)_z$$ and so the eigenvalue of $J^2$ is no good quantum number anymore. But $$[H, J_z] = 0,$$ so we can label the states of the system by their energy and their $J_z$ component. Indeed, if we calculate the action of $H$ on the former eigenbase, we see, that it only mixes $|1, 0\rangle$ and $|0, 0\rangle$. By diagonalizing the full Hamiltonian we find a new base of eigenstates:
$$|m=\pm 1, E_{\pm1} \rangle = |1, \pm 1\rangle\\ |m=0, E_{0,\pm}\rangle = \frac{1}{C_\pm}\left(\left(g\pm\sqrt{g^2+d^2}\right) |1, 0\rangle + d|0, 0\rangle\right),$$ where $C_\pm^2 = 2\left(g^2 \pm g \sqrt{g^2+d^2} + d^2\right)$ and corresponding energies $$E_{\pm1} = \frac{g}{2}\pm d\\ E_{0, \pm} = -\frac{g}{2} \pm \sqrt{g^2 + d^2},$$ which of course reduce to the eigenvalues of $H_0$ if you turn $d$ to zero.