Visualisierung von Grenzorbitalen mit Psi4 und Python

Jan 15 2023
Grenzorbitale bedeutet das höchste besetzte Molekülorbital (HOMO) und das niedrigste besetzte Orbital (LUMO). Die meisten Organischen Chemiker interessieren sich meist nur für diese Orbitale, weil dies der Punkt von Wechselwirkungen ist, an dem ein Elektron mit hoher Energie in das nächste Orbital springen könnte.

Grenzorbitale bedeutet das höchste besetzte Molekülorbital (HOMO) und das niedrigste besetzte Orbital (LUMO). Die meisten Organischen Chemiker interessieren sich meist nur für diese Orbitale, weil dies der Punkt von Wechselwirkungen ist, an dem ein Elektron mit hoher Energie in das nächste Orbital springen könnte. Sie sind die dominierenden Wechselwirkungen von Molekülen, die höchstwahrscheinlich miteinander interagieren würden. Ein Beispiel für ein HOMO/LUMO-Diagramm ist in Abbildung 1 dargestellt .

Abbildung 1: Beispiel eines HOMO/LUMO-Molekülorbitals

Wir wissen jedoch, dass das Elektron eine Welle ist und dass in diesen Orbitalen eine Elektronenwolke entsteht. Die Visualisierung dieser Wolken gibt uns eine Vorstellung davon, wie sich das Elektron verhält, und eine einprägsame Form, aus der wir vielleicht einige Dinge ableiten könnten. Kehren wir dazu zu unserer psi4-Methanol-Z-Matrix zurück, die wir im vorherigen Blog-Beitrag gefunden haben . Kommen wir gleich zum 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
)

Wir wollen also zwei Dinge, die Anzahl der Alpha-Elektronen und die Anzahl der Molekülorbitale, weil wir bestimmen müssen, welche Orbitale das HOMO/LUMO-Paar sind. Wir verwenden 2, weil wir uns an das Pauli-Ausschlussprinzip halten.

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()