teaching:progappchim:matplotlib_gallery:potentiel_energy_surface

Ceci est une ancienne révision du document !


Surface d'énergie potentielle

Eyring et Polanyi ont publié en 1931 l'article On Simple Gas Reactions dans lequel ils décrivent les trajets des atomes dans la réaction H2 + H –> H + H2 (échange d'atomes). Ces travaux aboutiront au développement des notions de complexe activé (activated complex) ou état de transition (transition state).

L'article “On a New Method of Drawing the Potential Energy Surface” (Shin Sato, J. Chem. Phys. 23, 592, 1955) présente une simplification relativement facile à mettre en oeuvre dans le cas où les 3 atomes d'hydrogène sont alignés.

Des expression analytiques sont proposées pour un état d'énergie liant et un état d'énergie non-liant :

  • $E_{bond}= D_e [\exp(-2\beta(r-r_e))-2\exp(-\beta(r-r_e))]$
  • $E_{ant}= \frac{D_e}{2} [\exp(-2\beta(r-r_e))+2\exp(-\beta(r-r_e))]$

$r_e$ est la distance interatomique d'équilibre de H2, $D_e$ la profondeur du puits de potentiel et $\beta$ un paramètre pour ajuster sa largeur (voir le Potentiel de Morseplugin-autotooltip__default plugin-autotooltip_bigPotentiel de Morse

Potentiel de Morse et approximation harmonique, avec représentation des niveaux d'énergie des modèles quantiques correspondants.

Code source :

#! /usr/bin/env python # -*- coding: utf-8 -*- """ Représentation du potentiel de Morse pour H2 http://en.wikipedia.org/wiki/Morse_potential http://en.wikipedia.org/wiki/Quantum_harmonic_oscillator approximation harmonique D_e = 7.6E-19 J a = 19.3E-15 m r_e= 74.1E-12 m dérivée de seconde d2V/dr2 = 2 * D_e * a**2. """ import matplot…
, et l'approximation harmonique).

Pour 2 atomes d'hydrogène A et B, une approximation est :

  • $E_{bond}= \frac{Q_{AB}+\alpha_{AB}}{1+S^2_{AB}} = \frac{Q_{AB}+\alpha_{AB}}{1+k}$
  • $E_{ant}= \frac{Q_{AB}-\alpha_{AB}}{1-S^2_{AB}} = \frac{Q_{AB}-\alpha_{AB}}{1-k}$

Où $k=S^2_{AB}$ et $Q_{AB}$, $\alpha_{AB}$ et $S_{AB}$ sont respectivement les intégrales de coulomb, d'échange et de recouvrement, toutes fonctions de la distance $r_{AB}$ entre les atomes A et B.

La solution proposée par Sato pour 3 atomes A, B, C, avec l'hypothèse $S^2_{AB}=S^2_{BC}=S^2_{CA}=k$ est :

  • $E = \frac{1}{1+k} \{ Q_{AB} + Q_{BC} + Q_{CA} - \sqrt{\frac{2}{1}[(\alpha_{AB} - \alpha_{BC})^2 + (\alpha_{BC} - \alpha_{CA})^2 + (\alpha_{CA} - \alpha_{AB})^2 ]} \}$

On obtient facilement $Q_{AB}$ et $\alpha_{AB}$ :

  • $Q_{AB} = ((1+k)E_{bond} + (1-k)E_{ant}) / 2$
  • $\alpha_{AB} = ((1+k)E_{bond} - (1-k)E_{ant}) / 2$

Sato présente des PES avec l'hypothèse k = 0.18 pour des distances jusque 0.5 nm.

<sxh python; title : PES-contour-01.py> #!/usr/bin/env python # -*- coding: utf-8 -*- “”“ Tracés de lignes de niveau ou isolignes Application : Potentiel Energy Surface de la réaction H + H2 –> H2 + H

”“” # ref : http://bulldog2.redlands.edu/facultyfolder/deweerd/tutorials/Tutorial-ContourPlot.pdf

import matplotlib.pyplot as plt # directive d'importation standard de Matplotlib import numpy as np # directive d'importation standard de numpy from mpl_toolkits.mplot3d import Axes3D # Axes3D

def Ebond(rAB):

  return D_e * (np.exp(-2.*beta*(rAB-r_e)) - 2.*np.exp(-beta*(rAB-r_e)))

def Eant(rAB):

  return 0.5 * D_e * (np.exp(-2.*beta*(rAB-r_e)) + 2.*np.exp(-beta*(rAB-r_e)))

def Q(rAB):

  return 0.5 * ((1.+k)*Ebond(rAB) + (1.-k)*Eant(rAB))

def a(rAB):

  return 0.5 * ((1.+k)*Ebond(rAB) - (1.-k)*Eant(rAB))

beta=19.3E-3 # pm-1 r_e=74.1 # pm D_e = .76 # E-18 J k=0.18 rmin=10. rmax=400. num=100 x_1d = np.linspace(rmin,rmax, num) print x_1d.shape, x_1d.dtype, x_1d.ndim y_1d = np.linspace(rmin,rmax, num) print y_1d.shape, y_1d.dtype, y_1d.ndim X, Y = np.meshgrid(x_1d, y_1d) print X.shape, X.dtype, X.ndim, Y.shape, Y.dtype, Y.ndim E=(Q(X)+Q(Y)+Q(X+Y)-np.sqrt(2.*1)2.+(a(Y)-a(X+Y))2.+(a(X+Y)-a(X))**2.)))/(1.+k) print np.min(E) #valeur minimale de E

fig = plt.figure(figsize=(12, 12), dpi=80) ax = fig.add_subplot(111) # cf. http://stackoverflow.com/questions/7965743/how-can-i-set-the-aspect-ratio-in-matplotlib ax.set_aspect(“equal”) levels = np.linspace(-1.7, 1.0, 53) CS1 = plt.contour(X, Y, E, levels, colors='k') plt.clabel(CS1, colors = 'k', fmt = '%2.2f', fontsize=14) CS2 = plt.contourf(X, Y, E, levels) #plt.colorbar(CS2) # visualisation éventuelle de l'échelle de couleur

plt.title('Isolignes') plt.xlabel('x (pm)') plt.ylabel('y (pm)')

fig = plt.figure(2,figsize=(15, 15)) ax = Axes3D(fig) ax.plot_surface(X,Y,E, rstride=1,cstride=1 ,cmap=plt.cm.jet) ax.set_xlabel('X') ax.set_ylabel('Y') ax.set_zlabel('E') plt.show() </sxh>

Avec les paramètres essayés, la valeur minimale de E est environ -1.603

Lignes de contour

Surface 3D

Voir aussi :


1)
a(X)-a(Y
Ce site web utilise des cookies. En utilisant le site Web, vous acceptez le stockage de cookies sur votre ordinateur. Vous reconnaissez également que vous avez lu et compris notre politique de confidentialité. Si vous n'êtes pas d'accord, quittez le site.En savoir plus
  • teaching/progappchim/matplotlib_gallery/potentiel_energy_surface.1462885327.txt.gz
  • Dernière modification: 2016/05/10 15:02
  • de villersd