Chemistry - How to calculate molecular dipole moment from a known wavefunction?

Solution 1:

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.

Solution 2:

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!

---

Solution 3:

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.