¿Cómo calcular el momento del dipolo molecular a partir de una función de onda conocida?

Question

Digamos que tengo una función de onda molecular como un conjunto de orbitales moleculares y quiero calcular el momento del dipolo de la molécula, ¡pero no sé cómo! He buscado mucho pero no he encontrado ningún ejemplo práctico.

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

Answer 1

La derivación formal necesaria ya ha sido bien hecha por AngusTheMan. Empezaré por la última ecuación:

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

donde $\Psi$ es la función de onda variacional; puede ser cualquier estado molecular. Es importante que sea variacional, de lo contrario el enfoque del valor de expectativa no es exacto. Por lo tanto, esto funciona para las funciones de onda SCF, CI y MCSCF, pero es necesario tomar derivadas adicionales para las funciones de onda Moller-Plesset y de clústeres acoplados. Hay que trabajar más para las funciones de onda multideterminadas como CI y MCSCF, pero la complejidad no es diferente para un solo estado en cada función de onda. Es posible que haya alguna partición del espacio de MO que esté descuidando y que sea necesaria para MCSCF, por lo que limitaré mi trabajo a una función de onda monodeterminada.

Expandir la función de onda como una combinación lineal de orbitales moleculares (MOs)

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

donde cada orbital molecular es una combinación lineal de orbitales atómicos (AOs)

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

donde $C_{\mu i}$ es la matriz de coeficientes de MO, por lo que nuestro valor de expectativa ahora se ve así:

$$ \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. $$

Los índices $\mu,\nu$ se ejecuta sobre todos los AOs, y el índice $i$ recorre las MO ocupadas. Sólo hay un índice porque este es un operador de un solo electrón. También estoy descuidando cualquier valor complejo aquí, ya que casi siempre trabajamos con AOs y coeficientes de MO de valor real.

Hacemos un último reordenamiento. Sustituir los coeficientes de MO por la matriz de densidad

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

para dar la primera "ecuación de trabajo" explícita:

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

Digo lo primero por dos razones: normalmente intentamos evitar bucles explícitos como éste, y la expresión puede descomponerse más, dependiendo de qué propiedades moleculares sean de interés; seré más claro sobre esto más adelante. Dentro de la suma hay dos términos:

• La matriz de densidad $P_{\mu\nu}$ que proviene de un cálculo convergente del SCF.
• La integral del operador dipolar sobre dos funciones base, $\langle \phi_{\mu} | \hat{\mu}_{z} | \phi_{\nu} \rangle$ . Los orbitales atómicos se representan como funciones base centradas en el átomo. Se pueden calcular una vez y en cualquier momento, ya que las cantidades aquí no cambian en el transcurso de un cálculo.

Cada uno de estos términos se representa como una matriz. Dado que el índice $\mu$ corre a lo largo de las filas y $\nu$ recorre las columnas de cada matriz, la "contracción" implica un producto matricial seguido de la traza, o un producto por elementos seguido de una suma de acumulación sobre todos los elementos de la matriz. Hay otros detalles con los que hay que tener cuidado, como en qué unidades debe estar el resultado, que cambia el prefactor (los programas trabajan internamente en unidades atómicas), y cuál es el origen para el operador dipolar, pero eso es todo.

Bueno, más o menos. En realidad, estoy tratando algunas de las partes internas del programa como una caja negra. Si estás familiarizado con Hartree-Fock, debería estar claro dónde $P_{\mu\nu}$ pero ¿qué pasa con la integral? Para un valor de expectativa general $\langle A \rangle$ con su correspondiente operador, donde hace $\langle \phi_{\mu} | \hat{A} | \phi_{\nu} \rangle$ ¿viene? Si ya está disponible en el código, entonces se llama a una función envolvente que luego llama al motor integral para hacer todo el trabajo sucio, y se obtiene una matriz ordenada sin tener que preocuparse por los detalles. Si $\langle \phi_{\mu} | \hat{A} | \phi_{\nu} \rangle$ no está presente, dependiendo de la complejidad de $\hat{A}$ En el caso de la ecuación integral de trabajo, puede ser necesaria una cantidad no trivial de derivación, seguida de la implementación.

Ignorando cualquier posible contracción de las funciones base primitivas, expande $\langle \phi_{\mu} | \hat{A} | \phi_{\nu} \rangle$ utilizando la definición de $\phi$ en coordenadas cartesianas:

$$ \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} $$

donde $\mathbf{r} = (x, y, z)$ es la posición del electrón, $\mathbf{A} = (A_x, A_y, A_z)$ es la posición de la función base (casi siempre centrada en el átomo), y $\mathbf{a} = (a_x, a_y, a_z)$ son los momentos angulares para cada coordenada, con $l_{\textrm{max}} = a_x + a_y + a_z$ momento angular total de la función base. $(0,0,0)$ es una función s, $(1,1,0)$ y $(0,0,2)$ son funciones d, y así sucesivamente.

Formando la integral de forma más explícita se obtiene

$$ \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] $$

Antes de ir más lejos, $\hat{A}$ debe definirse. Si $\hat{A} = 1$ se convierte en una integral de solapamiento. El operador dipolar en el z -está dada por $\hat{A} = \hat{\mu}_{z} = -ez = -e(z_3 - C_z)$ , donde $z_3$ es la coordenada de integración y $C_z$ es el origen del dipolo en el z -dirección, que normalmente se toma como cero. Todo se mantiene en unidades atómicas hasta después de la integral/contracción de la densidad, por lo que se deja de lado el prefactor $-e$ . Ahora podemos generalizar esto a un operador de momento multipolar cartesiano arbitrario,

$$ \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} $$

donde $(c_x, c_y, c_z)$ determinan la coordenada de cada multipolar, y su suma es el orden multipolar total; por ejemplo, $(1,0,0), (0,1,0), (0,0,1)$ son los x -, y -, y z -direcciones del operador dipolar. El operador de momento cartesiano se parece a una función base gaussiana donde $\zeta = 0$ .

Una vez que se ha derivado la forma de un operador, es necesario implementarlo como parte de un paquete de integrales, cada uno de los cuales implementa uno o más algoritmos para calcular integrales. Cada algoritmo lleva el nombre de los autores del artículo en el que se introdujeron, y suelen estar abreviados. Por ejemplo, el primero que conozco es el trabajo de Taketa, Huzinaga, O-Ohata (THO, DOI: 10.1143/JPSJ.21.2313 ), donde se dan ecuaciones de trabajo explícitas para las integrales de solapamiento de 2 centros, de energía cinética de 2 centros, de atracción electrónica-nuclear de 2 centros y de repulsión electrónica de 4 centros. Se puede encontrar una implementación de trabajo en el PyQuante paquete. Hice un Traducción del cuaderno IPython de los fragmentos de código en la primera página de la documentación oficial.

Otros algoritmos más complejos son los de 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) y Head-Gordon-Pople (HGP). Estoy seguro de que me estoy olvidando de algunos, incluido el artículo seminal de Boys que introdujo el uso de funciones gaussianas como sustituto de las funciones de tipo Slater en los conjuntos de bases. Una buena revisión de estos algoritmos se encuentra en un artículo de Peter Gill (DOI: 10.1016/S0065-3276(08)60019-2 ); es el autor original del código integral en ambos Gaussiano y Q-Chem .

Para cerrar el círculo, hace unos meses escribí un código para calcular el momento dipolar usando pyquante2 y un wrapper que llama a una implementación del algoritmo integral recursivo de Obara-Saika. Puedes encontrarlo en aquí con algunas comparaciones con programas cuánticos "industriales".

Answer 2

El momento dipolar $\mu$ de una molécula es una medida de la distribución de la carga en la molécula y la polaridad formada por los núcleos y la nube de electrones.

Podemos perturbar nuestro sistema con un campo eléctrico externo $\vec E$ y medir la respuesta de la nube de electrones y los núcleos mediante la polarizabilidad, es decir, cuánto cambia el momento dipolar. En la práctica, los núcleos pueden ser tan pesados que su movimiento no se ve perturbado, mientras que los electrones, al ser ligeros, son muy móviles. Si imaginamos que el campo externo es causado por alguna otra especie, y que él mismo no cambia, podemos llamar a este campo eléctrico externo constante $\vec E$ al menos sobre el volumen de la molécula que estamos considerando. Imaginemos, para llevar una contabilidad arbitraria, que apuntamos hacia abajo el $z$ eje. También podríamos investigar cómo cambia el momento dipolar con las vibraciones de los enlaces para hablar de la espectroscopia IR o si la polarizabilidad cambia durante una vibración para dar la espectroscopia Raman.

Podemos utilizar teoría de la perturbación para expandir la función de onda y la energía molecular en términos de pequeñas perturbaciones del campo. Empezamos por Taylor expandiéndose la energía y la función de onda molecular en términos del campo eléctrico que actúa como parámetro de perturbación. \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} Si utilizamos la notación de que las derivadas de la función de onda están dadas por; \begin{equation} \bigg(\frac {1}{i!}\bigg)\frac{\partial ^i\psi}{\partial \vec E^i}=\psi ^{(i)} \end{equation} El Hamiltoniano para tal sistema bajo la influencia de un campo eléctrico en el $z$ dirección es $\hat H(\vec E)$ . \begin{equation} \hat H(\vec E)=\hat H^0-\vec E\hat \mu _z \end{equation} Dónde $\hat \mu _z$ es el operador de momento dipolar que es una suma de las cargas de los núcleos y los electrones de la molécula. Este es el resultado de Teorema de Hellmann-Feynman de la energía y el campo eléctrico, por ejemplo $\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}

La ecuación de Schrödinger independiente del tiempo es ahora, \begin{equation} \hat H(\vec E)\psi(\vec E)=E(\vec E)\psi (\vec E) \end{equation} Con una energía; \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}

Con un poco de álgebra y el uso de $E^{(0)}=\langle \psi ^{(0)}|\hat H^0 |\psi ^{(0)}\rangle$ así como las propiedades hermitianas del hamiltoniano. \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} Utilizando la ecuación de Schrödinger $H\psi =E\psi$ y sacando la energía escalar de la 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} Desde $\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} Por lo tanto, el valor esperado del momento dipolar a lo largo del $z$ eje para un estado molecular $\psi ^{(0)}$ es $\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} Para comprender la fuerza de la interacción que provoca la transición entre los estados $\psi ^{(0)}$ y $\psi ^{(1)}$ utilizamos el momento dipolar de transición, que es básicamente el mismo, pero hay una función de onda de los dos estados implicados, ¡inicial y final!

Si se repitiera este proceso pero se mantuviera un orden más alto (¡desordenado!) se obtendría (2º orden) la polarizabilidad de la molécula que, en esencia, es la susceptibilidad de la nube de electrones a cambiar con respecto a un campo eléctrico externo (por tanto, cómo cambia el momento dipolar. El tercer orden daría la hiperpolarizabilidad, etc., etc...

Como he dicho, también se podría enfocar esto desde un ángulo realmente diferente, interpretando el diagrama de orbitales moleculares y utilizando la química computacional (así principio variacional etc.) para encontrar los coeficientes orbitales moleculares. Eso te daría una buena idea de lo que está pasando.

Answer 3

Un enfoque que se utiliza a menudo, especialmente para los métodos semiempíricos, es utilizar un conjunto de cargas centradas en el átomo (la mayoría de las veces las cargas de Mulliken) para calcular el momento dipolar. En ese caso, el momento dipolar molecular viene dado por:

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

Tenga en cuenta que si el sistema tiene una carga finita, esta ecuación depende de la posición del sistema molecular con respecto al origen del espacio de coordenadas. La mayoría de los programas utilizan el centro de masa o el centro de carga nuclear como origen en estos casos.