Visualizzare gli orbitali di frontiera con Psi4 e Python

Per orbitali di frontiera si intendono l'orbitale molecolare più alto occupato (HOMO) e l'orbitale più basso occupato (LUMO). La maggior parte dei chimici organici è per lo più interessata solo a questi orbitali perché questi sono il punto di interazione in cui un elettrone ad alta energia potrebbe saltare nell'orbitale successivo. Sono le interazioni dominanti di molecole che molto probabilmente sarebbero tra loro. Un esempio di diagramma HOMO/LUMO è mostrato in Figura 1 .

Tuttavia, sappiamo che l'elettrone è un'onda e c'è una nuvola di elettroni che viene prodotta in questi orbitali. Visualizzare queste nuvole ci dà un'idea di come si comporta l'elettrone e una forma memorabile potremmo forse dedurre alcune cose. Per fare questo torniamo alla nostra matrice z di metanolo psi4, e come trovato nel precedente post sul blog . Entriamo subito nel codice:

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
)

Quindi vogliamo due cose, il numero di elettroni alfa e il numero di orbitali molecolari perché dobbiamo determinare quali orbitali sono la coppia HOMO/LUMO. Usiamo 2 perché aderiamo al principio di esclusione di 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()