การแสดงภาพ Frontier Orbitals ด้วย Psi4 และ Python

Jan 15 2023
Frontier Orbitals หมายถึงวงโคจรระดับโมเลกุลที่ถูกยึดครองสูงสุด (HOMO) และวงโคจรที่ถูกยึดครองต่ำสุด (LUMO) นักเคมีอินทรีย์ส่วนใหญ่สนใจเฉพาะออร์บิทัลเหล่านี้เท่านั้น เนื่องจากเป็นจุดที่มีปฏิสัมพันธ์ซึ่งอิเล็กตรอนที่มีพลังงานสูงสามารถกระโดดเข้าไปในออร์บิทัลถัดไปได้

Frontier Orbitals หมายถึงวงโคจรระดับโมเลกุลที่ถูกยึดครองสูงสุด (HOMO) และวงโคจรที่ถูกยึดครองต่ำสุด (LUMO) นักเคมีอินทรีย์ส่วนใหญ่สนใจเฉพาะออร์บิทัลเหล่านี้เท่านั้น เนื่องจากเป็นจุดที่มีปฏิสัมพันธ์ซึ่งอิเล็กตรอนที่มีพลังงานสูงสามารถกระโดดเข้าไปในออร์บิทัลถัดไปได้ พวกมันเป็นปฏิสัมพันธ์ที่โดดเด่นของโมเลกุลที่น่าจะเป็นไปได้มากที่สุด ตัวอย่างของไดอะแกรม HOMO/LUMO แสดงในรูปที่ 1

รูปที่ 1: ตัวอย่างการโคจรของโมเลกุล HOMO/LUMO

อย่างไรก็ตาม เราทราบดีว่าอิเล็กตรอนเป็นคลื่นและมีเมฆอิเล็กตรอนเกิดขึ้นที่ออร์บิทัลเหล่านี้ การมองเห็นเมฆเหล่านี้ทำให้เรามีความคิดบางอย่างเกี่ยวกับพฤติกรรมของอิเล็กตรอนและรูปร่างที่น่าจดจำซึ่งเราอาจอนุมานบางสิ่งได้ ในการทำเช่นนี้ ให้กลับไปที่ psi4 เมทานอล z-matrix ของเรา และตามที่พบในโพสต์บล็อก ก่อน หน้า มารับรหัสกัน:

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
)

เราต้องการสองสิ่ง จำนวนอิเล็กตรอนแอลฟาและจำนวนออร์บิทัลของโมเลกุล เพราะเราต้องกำหนดว่าออร์บิทัลใดเป็นคู่โฮโม/ลูโม เราใช้ 2 เพราะเราปฏิบัติตามหลักการกีดกันของ 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()