Visualizando orbitais de fronteira com Psi4 e Python

Jan 15 2023
Frontier Orbitals significa o Orbital Molecular Ocupado Mais Alto (HOMO) e o Orbital Ocupado Mais Baixo (LUMO). A maioria dos químicos orgânicos está interessada apenas nesses orbitais, porque esses são os pontos de interação em que um elétron com alta energia pode pular para o próximo orbital.

Frontier Orbitals significa o Orbital Molecular Ocupado Mais Alto (HOMO) e o Orbital Ocupado Mais Baixo (LUMO). A maioria dos químicos orgânicos está interessada apenas nesses orbitais, porque esses são os pontos de interação em que um elétron com alta energia pode pular para o próximo orbital. Eles são as interações dominantes de moléculas que provavelmente seriam umas com as outras. Um exemplo de um diagrama HOMO/LUMO é mostrado na Figura 1 .

Figura 1: Exemplo de Orbital Molecular HOMO/LUMO

No entanto, sabemos que o elétron é uma onda e há uma nuvem de elétrons que é produzida nesses orbitais. A visualização dessas nuvens nos dá uma ideia de como o elétron se comporta e uma forma memorável que talvez possamos inferir algumas coisas. Para fazer isso, vamos voltar à nossa matriz z de metanol psi4, conforme encontrado na postagem anterior do blog . Vamos direto ao código:

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
)

Então, queremos duas coisas, o número de elétrons alfa e o número de orbitais moleculares, porque precisamos determinar quais orbitais são o par HOMO/LUMO. Usamos 2 porque aderimos ao Princípio de Exclusão 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()