Visualiser les orbitales frontières avec Psi4 et Python

Jan 15 2023
Frontier Orbitals désigne l'orbite moléculaire la plus occupée (HOMO) et l'orbite la plus basse occupée (LUMO). La plupart des chimistes organiques ne s'intéressent principalement qu'à ces orbitales car ce sont le point d'interaction où un électron à haute énergie pourrait sauter dans l'orbitale suivante.

Frontier Orbitals désigne l'orbite moléculaire la plus occupée (HOMO) et l'orbite la plus basse occupée (LUMO). La plupart des chimistes organiques ne s'intéressent principalement qu'à ces orbitales car ce sont le point d'interaction où un électron à haute énergie pourrait sauter dans l'orbitale suivante. Ce sont les interactions dominantes des molécules qui seraient très probablement les unes avec les autres. Un exemple de diagramme HOMO/LUMO est illustré à la figure 1 .

Figure 1 : Exemple d'une orbitale moléculaire HOMO/LUMO

Cependant, nous savons que l'électron est une onde et qu'un nuage d'électrons est produit dans ces orbitales. La visualisation de ces nuages ​​nous donne une idée du comportement de l'électron et une forme mémorable que nous pourrions peut-être déduire. Pour ce faire, revenons à notre matrice z de méthanol psi4, et comme trouvé dans le billet de blog précédent . Entrons directement dans le code :

import psi4

psi4.set_memory('1000mb')
psi4.core.set_num_threads(1)

zmatrix = '''\
    H11
    O11  H11  0.9600
    C11  O11  1.4000  H11  108.0000
    H12  C11  1.1000  O11  112.0000  H11  -60.0000
    H13  C11  1.1000  O11  112.0000  H11   60.0000
    H14  C11  1.1000  O11  111.8699  H11 -180.0000
'''

universe = psi4.geometry(zmatrix)
universe.update_geometry()
universe.print_out()

energy, wave_function = psi4.optimize(
  'hf/6-31g*',
   return_wfn =True,
  molecule=universe
)

Nous voulons donc deux choses, le nombre d'électrons alpha et le nombre d'orbitales moléculaires car nous devons déterminer quelles orbitales sont la paire HOMO/LUMO. Nous utilisons 2 parce que nous adhérons au principe d'exclusion de Pauli.

number_of_occupied = 2
number_of_virutal = 2

alpha_electrons = wave_function.nalpha()
molecule = wave_function.nmo()

min_orb = max(1, alpha_electrons + 1 - number_of_occupied)
max_orb = min(molecule, alpha_electrons + number_of_virutal)
orbitals = [k for k in range(min_orb, max_orb + 1)]

cubeprop_tasks = []
  cubeprop_tasks.append('DENSITY')
  cubeprop_tasks.append('ORBITALS')

  psi4.set_options(
    {
      'scf_type': 'df',
      'g_convergence': 'gau_tight',
      'freeze_core': 'true',
      'CUBEPROP_TASKS': cubeprop_tasks,
      'CUBEPROP_FILEPATH': './',
      'CUBEPROP_ORBITALS': orbitals,
    }
  )

  psi4.cubeprop(wave_function)

Da.cube
Psi_a_9_2-A.cube

pip install moly
pip install opt_einsum

import moly

fig = moly.Figure()
fig.add_cube(
    'Psi_a_9_2-A.cube',
    iso=0.03,
    colorscale="rdbu",
    opacity=0.1
)

fig.show()