How to calculate molecular dipole moment from a known wavefunction?
The necessary formal derivation has already been nicely done by AngusTheMan. I'll start from the last equation:

$$
\langle \mu_{z} \rangle = \langle \Psi | \hat{\mu}_{z} | \Psi \rangle
$$

where $\Psi$ is the variational wavefunction; it can be any molecular state. It's important that it's variational, otherwise the expectation value approach is not exact. So, this works for SCF, CI, and MCSCF wavefunctions, but extra derivatives need to be taken for Moller-Plesset and coupled cluster wavefunctions. More work needs to be done for multideterminental wavefunctions like CI and MCSCF, but the complexity is no different for a single state in each wavefunction. There may be some MO space partitioning I'm neglecting that's required for MCSCF, so I'll restrict my work to a single-determinental wavefunction.

Expand the wavefunction as a linear combination of molecular orbitals (MOs)

$$
\Psi = \sum_{i} \psi_{i},
$$

where each molecular orbital is a linear combination of atomic orbitals (AOs)

$$
\psi_{i} = \sum_{\mu} C_{\mu i} \phi_{\mu},
$$

where $C_{\mu i}$ is the MO coefficient matrix, so our expectation value now looks like this:

$$
\langle \mu_{z} \rangle = \sum_{i}^{\textrm{occ MOs}} \sum_{\mu\nu}^{\textrm{AOs}} C_{\mu i} C_{\nu i} \langle \phi_{\mu} | \hat{\mu}_{z} | \phi_{\nu} \rangle.
$$

The indices $\mu,\nu$ run over all AOs, and the index $i$ runs over the occupied MOs. There's only one index because this is a one-electron operator. I'm also neglecting any complex values here, since we almost always work with real-valued AOs and MO coefficients.

We do one last rearrangement. Replace the MO coefficients with the density matrix

$$
P_{\mu\nu} = \sum_{i}^{\textrm{occ MOs}} C_{\mu i} C_{\nu i}
$$

to give the first explicit "working equation":

$$
\langle \mu_{z} \rangle = \sum_{\mu\nu}^{\textrm{AOs}} P_{\mu\nu} \langle \phi_{\mu} | \hat{\mu}_{z} | \phi_{\nu} \rangle
$$

I say first for two reasons: we usually try and avoid explicit loops like this, and the expression can be broken down further, depending on what molecular properties are of interest; I'll be more clear about this later. Inside the sum there are two terms:


The density matrix $P_{\mu\nu}$, which comes from a converged SCF calculation.
The integral of the dipole operator over two basis functions, $\langle \phi_{\mu} | \hat{\mu}_{z} | \phi_{\nu} \rangle$. Atomic orbitals are represented as atom-centered basis functions. These can be calculated once and at any time, since the quantities here don't change over the course of a calculation.


Each of these terms is represented as a matrix. Since the index $\mu$ runs along the rows and $\nu$ runs along the columns for each matrix, "contraction" involves either a matrix product followed by the trace, or an elementwise product followed by an accumulation sum over all matrix elements. There are other details one needs to be careful about, such as what units the result should be in, which changes the prefactor (programs work internally in atomic units), and what the origin for the dipole operator is, but that's really it.

Well, sort of. I'm actually treating some of the program internals as a black box. If you're familiar with Hartree-Fock, it should be clear where $P_{\mu\nu}$ comes from, but what about the integral? For a general expectation value $\langle A \rangle$ with its corresponding operator, where does $\langle \phi_{\mu} | \hat{A} | \phi_{\nu} \rangle$ come from? If it's already available in the code, then you call a wrapper function that then calls the integral engine to do all the nasty work, and you get back a tidy matrix without having to worry about the details. If $\langle \phi_{\mu} | \hat{A} | \phi_{\nu} \rangle$ isn't present, depending on the complexity of $\hat{A}$, there can be a non-trivial amount of derivation required for the working integral equation, followed by the implementation.

Ignoring any possible contraction of primitive basis functions, expand $\langle \phi_{\mu} | \hat{A} | \phi_{\nu} \rangle$ using the definition of $\phi$ in Cartesian coordinates:

$$
\phi(\mathbf{r}; \mathbf{A}, \mathbf{a}, \zeta) = (x-A_x)^{a_x} (y-A_y)^{a_y} (z-A_z)^{a_z} e^{-\zeta |\mathbf{r} - \mathbf{A}|^2}
$$

where $\mathbf{r} = (x, y, z)$ is the electron position, $\mathbf{A} = (A_x, A_y, A_z)$ is the position of the basis function (almost always atom-centered), and $\mathbf{a} = (a_x, a_y, a_z)$ are the angular momenta for each coordinate, with $l_{\textrm{max}} = a_x + a_y + a_z$ total angular momentum of the basis function. $(0,0,0)$ is an s-function, $(1,1,0)$ and $(0,0,2)$ are d-functions, and so on.

Forming the integral more explicitly gives

$$
\langle \phi_{\mu} | \hat{A} | \phi_{\nu} \rangle = \int\int d\mathbf{r}_1 d\mathbf{r}_2 \left[ (x_1-A_x)^{a_x} (y_1-A_y)^{a_y} (z_1-A_z)^{a_z} e^{-\zeta_a |\mathbf{r}_1 - \mathbf{A}|^2} \right] \\ \times\left[ \hat{A} \right]\left[ (x_2-B_x)^{b_x} (y_2-B_y)^{b_y} (z_2-B_z)^{b_z} e^{-\zeta_b |\mathbf{r}_2 - \mathbf{B}|^2} \right]
$$

Before going any further, $\hat{A}$ must be defined. If $\hat{A} = 1$, this becomes an overlap integral. The dipole operator in the z-direction is given by $\hat{A} = \hat{\mu}_{z} = -ez = -e(z_3 - C_z)$, where $z_3$ is the integration coordinate and $C_z$ is the origin of the dipole in the z-direction, usually taken to be zero. Everything is kept in atomic units until after the integral/density contraction, so drop the prefactor $-e$. We can now generalize this to an arbitrary Cartesian multipole moment operator,

$$
\hat{A} = \mathfrak{M}(\mathbf{r}_3) = (x_3 - C_x)^{c_x} (y_3 - C_y)^{c_y} (z_3 - C_z)^{c_z}
$$

where $(c_x, c_y, c_z)$ determine the coordinate of each multipole, and their sum is the total multipole order; for example, $(1,0,0), (0,1,0), (0,0,1)$ are the x-, y-, and z-directions of the dipole operator. The Cartesian moment operator looks just like a Gaussian basis function where $\zeta = 0$.

Once the form of an operator has been derived, it needs to be implemented as part of an integral package, each of which implement one or more algorithms for computing integrals. Each algorithm is named after the authors of the paper in which they were introduced, and are usually abbreviated. For example, the first one I know of is the Taketa, Huzinaga, O-Ohata paper (THO, DOI: 10.1143/JPSJ.21.2313), where explicit working equations are given for 2-center overlap, 2-center kinetic energy, 2-center electron-nuclear attraction, and 4-center electron repulsion integrals. A working implementation can be found in the PyQuante package. I made an IPython notebook translation of the code snippets on the front page of the official documentation.

Other, more complex algorithms are from the Pople-Hehre (PH), McMurchie-Davidson (MD, DOI: 10.1016/0021-9991(78)90092-X), Obara-Saika (OS, DOI: 10.1063/1.450106), Dupuis-Rys-King (DRK), and Head-Gordon-Pople (HGP) papers. I'm sure I'm neglecting some, including the seminal paper by Boys which introduced the use of Gaussian functions as a substitute for Slater-type functions in basis sets. A good review of these algorithms is found in a Peter Gill paper (DOI: 10.1016/S0065-3276(08)60019-2); he is the original author of the integral code in both Gaussian and Q-Chem.

To bring these things full circle, I wrote some code a few months ago to calculate the dipole moment using pyquante2 and a wrapper that calls an implementation of the Obara-Saika recursive integral algorithm. You can find it here with some comparisons to "industrial" quantum programs.