# How to calculate molecular dipole moment from a known wavefunction?

Say I have a molecular wavefunction as a set of molecular orbitals and want to calculate the molecule's dipole moment, but don't know how! I searched a lot but couldn't find any practical example.

$$\psi _i=\sum ^N_{i=1}C_i\mathrm e^{-\alpha _ir^2}$$

• You could use the variational theorem to determine the molecular orbital coefficients of the molecular orbitals (one electron wavefunctions). If you let you molecular wavefunctions be a linear superposition of basis atomic wavefunctions $\Psi =\sum c_i\psi_i$, with the orbital coefficients you can understand key properties of your molecule. You can rationalise trends in bond polarity too, that can't be explained with other theories! – AngusTheMan May 9 '15 at 1:19

The dipole moment $\mu$ of a molecule is a measure of charge distribution in the molecule and the polarity formed by the nuclei and electron cloud.

We can perturb our system with an external electric field $\vec E$ and gauge the response of the electron cloud and nuclei by the polarisability, i.e how much the dipole moment changes. In practice the nuclei might be so heavy that their motion is not perturbed, while electrons being light are very mobile. If we imagine that the external field is caused by some other species, and that it itself is not changing, we can call this external constant electric field $\vec E$ at least over the volume of the molecule we are considering. Imagine for arbitrary book keeping that we point it down the $z$ axis. We could also investigate how the dipole moment changes with bond vibrations to discuss IR spectroscopy or if the polarisability changes during a vibration to give Raman spectroscopy.

We can use perturbation theory to expand the wavefunction and the molecular energy in terms of small perturbations of the field. We start by Taylor expanding the energy and molecular wavefunction in terms of the electric field which acts as the perturbation parameter. \begin{equation} E(\vec E)=E^0+\bigg(\frac{\partial E}{\partial \vec E}\bigg)_0\vec E+\bigg(\frac{\partial ^2E}{\partial \vec E^2}\bigg)_0\frac{\vec E^2}{2!}+\bigg(\frac{\partial ^3E}{\partial \vec E^3}\bigg)_0\frac{\vec E^3}{3!}+\dots \end{equation} \begin{equation} \psi(\vec E)=\psi^0+\bigg(\frac{\partial \psi}{\partial \vec E}\bigg)_0\vec E+\bigg(\frac{\partial ^2\psi}{\partial \vec E^2}\bigg)_0\frac{\vec E^2}{2!}+\bigg(\frac{\partial ^3\psi}{\partial \vec E^3}\bigg)_0\frac{\vec E^3}{3!}+\dots \end{equation} If we use the notation that the wavefunction derivatives are given by; \begin{equation} \bigg(\frac {1}{i!}\bigg)\frac{\partial ^i\psi}{\partial \vec E^i}=\psi ^{(i)} \end{equation} The Hamiltonian for such a system under the influence of an electric field in the $z$ direction is $\hat H(\vec E)$. \begin{equation} \hat H(\vec E)=\hat H^0-\vec E\hat \mu _z \end{equation} Where $\hat \mu _z$ is the dipole moment operator that is a summation of the charges of the nuclei and electrons in the molecule. This is the result of Hellmann-Feynman theorem of the energy and the electric field, e.g. $\hat H(\vec E)=\hat H^0 +\hat H^1(\vec E)$, \begin{equation} \frac{dE}{d\vec E}=\bigg\langle \frac{d\hat H}{d\vec E}\bigg\rangle=\bigg\langle \frac{d(-\mu _z\vec E)}{d\vec E}\bigg\rangle \end{equation}

The time-independent Schrödinger equation is now, \begin{equation} \hat H(\vec E)\psi(\vec E)=E(\vec E)\psi (\vec E) \end{equation} With an energy; \begin{equation} E(\vec E)=\big \langle \psi (\vec E)\big|\hat H(\vec E)\big|\psi (\vec E)\big \rangle=\big\langle \psi ^{(0)}+\psi ^{(1)}\vec E+\psi ^2 \vec E^2+ \dots \big|\hat H-\vec E\hat \mu _z \big|\psi ^{(0)}+\psi ^{(1)}\vec E+\psi^2\vec E^2+\dots \big \rangle \end{equation}

With a little algebra and use of $E^{(0)}=\langle \psi ^{(0)}|\hat H^0 |\psi ^{(0)}\rangle$, as well as the Hermitian properties of the Hamiltonian. \begin{equation} E(\vec E)=E^{(0)}+2\vec E\big \langle \psi ^{(1)}\big|\hat H^{(0)}\big|\psi ^{(0)}\rangle -\vec E\big\langle \psi ^{(0)}\big|\hat \mu _z\big|\psi ^{(0)}\big \rangle +\mathcal O(\vec E^2) \end{equation} Using the Schrödinger equation $H\psi =E\psi$ and pulling the scalar energy out of the integral; \begin{equation} \big\langle \psi ^{(1)}\big|\hat H^0\big|\psi ^{(0)}\rangle =E^{(0)}\big\langle \psi ^{(1)}\big|\psi ^{(0)}\big \rangle \end{equation} Since $\langle \phi |\psi \rangle=0$, \begin{equation} E(\vec E)=E^{(0)}-\vec E\big \langle \psi ^{(0)}\big|\hat \mu _z\big|\psi ^{(0)}\big\rangle \end{equation} Therefore the expectation value of the dipole moment along the $z$ axis for a molecular state $\psi ^{(0)}$ is $\langle \mu _z \rangle$; \begin{equation} \langle \mu _z\rangle =\big\langle \psi ^{(0)}\big|\hat \mu _z\big|\psi ^{(0)}\big\rangle \end{equation} To understand the strength of the interaction that causes the transition between the states $\psi ^{(0)}$ and $\psi ^{(1)}$ we use the transition dipole moment which is basically exactly the same but there is a wavefunction from both of the states involved, initial and final!

If you were to repeat this process but retain higher orders (Messy!) you would get (2nd order) the polarisability of the molecule which in essence is the susceptibility of the electron cloud to change with respect to an external electric field (so how the dipole moment changes. Third order would give the hyperpolarisability etc., etc...

As I said, you could also approach this from a really different angle by interpreting the molecular orbital diagram and using computational chemistry (so variational principle etc) to find the molecular orbital coefficients! That would give you a good idea of what is going on!

• Hrm. This is a quite thorough formal elaboration of the theory involved, but it doesn't actually describe how one would implement practically the dipole moment calculation. (At least, I am no closer to understanding how I would code a dipole moment calculation given a set of MO coefficients and basis functions.) – hBy2Py Sep 3 '15 at 11:57
• This whole derivation is unnecessary. Per definition the dipole moment is the expectation value of the dipole moment operator with the given wave-function. – Greg Sep 3 '15 at 14:38
• Thank you both for your comments, I will take this into account and update my answer shortly. I agree with @Brian, this answer does not give a non-specialist sufficient information on how to perform calculations to obtain values for molecules etc. However, I do not feel it is obvious that the dipole moment is the expectation value to someone who is new to this material and not as knowledgable as others. It is my experience that undergrads can struggle in this area if they can not see where something comes from. That's why I always like to give people a little background material. :) – AngusTheMan Sep 3 '15 at 15:39
• @Greg For those like myself (background in chemical engineering, and not applied mathematics / quantum physical chemistry; but with an interest in a general understanding of the inner workings of quantum computation) it is perhaps not quite as sad...? I very much appreciate the efforts of both AngusTheMan and pentavalentcarbon to lay out the details. Also, frankly: isn't this sort of exposition the entire purpose of StackExchange? – hBy2Py Sep 22 '15 at 13:56
• @Greg Practical quantum computation was the topic of the original question -- so, no, it's not at all off-topic. Also, if you're looking for mathematically rigorous developments from axiom to theorem or whatever, you're on the wrong SE site. There are Math.SE and MathOverflow for that. – hBy2Py Sep 22 '15 at 14:11

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.

• Thanks all. I really recommend reading the Obara-Saika paper if you're at all interested in integral evaluation. The DALTON package (not sure what the primary integral algorithm is) has an impressive list of one-electron integrals that can be calculated, many of which correspond directly (expectation value times some constants) to molecular properties. – pentavalentcarbon Sep 20 '15 at 1:29

An often used approach, especially for semi-empirical methods, is to use a set of atom-centered charges (most often the Mulliken charges) to calculate the dipole moment. In that case, the molecular dipole moment is given by:

$\vec{\mu} = \sum_a \vec{r} \times q_a$

Note that if the system has a finite charge, this equation has a dependence on the position of the molecular system relative to the origin of the coordinate space. Most program use the center of mass or center of nuclear charge as origin in such cases.