דוגמאות לחישובים נומריים בכימיה

במחברת זו נשתמש ב־NumPy כדי לבצע שני חישובים נומריים קצרים שמופיעים לעיתים קרובות בכימיה: חישוב שטח מתחת לפס בליעה וחישוב אוכלוסיות בולצמן של רמות אנרגיה. המטרה אינה ללמוד כימיה חדשה, אלא לראות איך מערכים של NumPy מאפשרים לבצע חישובים על הרבה נקודות בבת אחת, בלי לולאות מיותרות.

🤖 Ask PyProf.

🗑 −

▶ שלח

⏳ ממתין...

import numpy as np

דוגמה 1: שטח מתחת לפס בליעה

בספקטרוסקופיה, לעיתים קרובות השטח מתחת לפס בליעה חשוב לא פחות מגובה הפס. לדוגמה, השטח יכול להיות קשור לכמות החומר הבולע או לעוצמת המעבר.

נבנה פס בליעה מלאכותי בצורת גאוסיאן:

\[ A(\lambda) = A_0 \exp\left[-\frac{(\lambda-\lambda_0)^2}{2\sigma^2}\right] \]

כאן \(\lambda\) הוא אורך הגל, \(A(\lambda)\) היא הבליעה, \(A_0\) הוא גובה הפס, \(\lambda_0\) הוא מרכז הפס, ו־\(\sigma\) קובע את רוחב הפס.

ל־NumPy יש הרבה פונקציות מתמטיות שימושיות. כבר ראינו למשל את np.sqrt, ויש גם פונקציות טריגונומטריות כמו np.sin ו־np.tan, פונקציות מעריכיות ולוגריתמיות כמו np.exp ו־np.log, ועוד רבות אחרות. הרשימה המלאה מופיעה בתיעוד של NumPy: https://numpy.org/doc/stable/reference/routines.math.html.

נשתמש גם ב־np.linspace(start, stop, num). הפונקציה הזו יוצרת מערך של num ערכים שווי־מרחק בין start ל־stop. למשל, np.linspace(400, 700, 301) יוצר 301 אורכי גל בין 400 ל־700 ננומטר. זה נוח כאשר רוצים לחשב פונקציה על רשת צפופה ומסודרת של ערכי \(x\).

בקוד נשתמש בפונקציה np.exp, שמחשבת את \(e^x\). כמו פונקציות NumPy אחרות, היא עובדת גם על מספר אחד וגם על מערך שלם. כאשר מכניסים לה מערך, כמו מערך אורכי הגל שלנו, היא מחזירה מערך באותו גודל: ערך אחד לכל אורך גל.

wavelength = np.linspace(400, 700, 301)  # nm

A0 = 1.0        # peak height
lambda0 = 520   # nm
sigma = 20      # nm

absorbance = A0 * np.exp(-(wavelength - lambda0)**2 / (2 * sigma**2))

כעת נחשב בקירוב את האינטגרל של הבליעה לפי אורך הגל:

\[\text{area} = \int A(\lambda)\,d\lambda\]

בנתונים ניסיוניים אין לנו פונקציה רציפה אלא נקודות מדידה בדידות. לכן נשתמש בכלל הטרפזים, שממומש ב־NumPy על־ידי np.trapezoid.

area = np.trapezoid(absorbance, wavelength)

print(f"area = {area:.2f} absorbance nm")

area = 50.13 absorbance nm

דוגמה 2: אוכלוסיות בולצמן

במערכת בשיווי־משקל תרמי, לא כל המולקולות חייבות להיות ברמת האנרגיה הנמוכה ביותר. עבור רמות אנרגיה \(E_i\), המשקל התרמי של רמה \(i\) פרופורציוני ל־

\[e^{-E_i/(k_B T)}\]

כדי לקבל הסתברויות, צריך לנרמל את המשקלים כך שסכום האוכלוסיות יהיה 1.

kB = 8.617333262e-5  # eV/K
T = 300              # K

energies = np.array([0.00, 0.03, 0.08, 0.15])  # eV

weights = np.exp(-energies / (kB * T))
populations = weights / np.sum(weights)

print(populations)
print("sum =", np.sum(populations))

[0.73439706 0.23011934 0.0332652  0.00221841]
sum = 0.9999999999999998

החישוב נעשה על כל רמות האנרגיה בבת אחת. כמו בדוגמת פס הבליעה, np.exp מקבלת כאן מערך ומחזירה מערך שלם של משקלים, אחד לכל רמת אנרגיה.

נציג את התוצאה בצורה נוחה יותר.

for i in range(len(energies)):
    energy = energies[i]
    population = populations[i]
    print(f"E = {energy:.2f} eV, population = {population:.3f}")

E = 0.00 eV, population = 0.734
E = 0.03 eV, population = 0.230
E = 0.08 eV, population = 0.033
E = 0.15 eV, population = 0.002

כעת נבדוק מה קורה כאשר הטמפרטורה משתנה. בטמפרטורה נמוכה יותר, רוב האוכלוסייה תהיה ברמה הנמוכה ביותר. בטמפרטורה גבוהה יותר, רמות מעוררות יקבלו אוכלוסייה גדולה יותר.

def boltzmann_populations(energies, temperature):
    """Return normalized Boltzmann populations for energies in eV."""
    weights = np.exp(-energies / (kB * temperature))
    return weights / np.sum(weights)
temperatures = [100, 300, 1000]
for temperature in temperatures:
    populations = boltzmann_populations(energies, temperature)
    print(f"T = {temperature} K")
    print(populations)
    print()

T = 100 K
[9.70065106e-01 2.98447144e-02 9.01526163e-05 2.67382885e-08]

T = 300 K
[0.73439706 0.23011934 0.0332652  0.00221841]

T = 1000 K
[0.43925047 0.31011223 0.17359211 0.07704519]