Menggambar Vibrasi sebuah Molekul Air

April 19, 2018   kimia , python3

Halo! 

Kali ini saya ingin berbagi tulisan tentang pembuatan grafik lewat paket kode pyplot dalam bahasa Python. Kode yang saya tulis ini berfungsi untuk menyajikan mode vibrasi molekul secara visual dalam dua dimensi. Data vektor tiap atom dari molekul didapatkan dari hasil perhitungan ab initio via Gaussian 09. Batasan metode ini adalah molekul dan mode vibrasi harus terdefinisi pada sumbu x dan y saja. Di sini saya menggunakan paket Python mathplotlib.pyplot dan numpy. Data yang saya gunakan adalah data penghitungan struktur molekul air (H₂O) menggunakan perangkat Gaussian 09. 

"""
 program to plot vibration modes
 planar, triatomic system
 exp/physchem/54E3
 Alwan Darusslam
 2018/4/18
"""

import numpy as np
import matplotlib.pyplot as plt

# molecule: H2O
# stdo: standard orientation #
stdo = np.array([[0.,0.,0.,0.106840],
                [0.,0.,0.785221,-0.427359],
                [0.,0.,-0.785221,-0.427359]])

# vibration modes (A11, A12, B11)
A11 = np.array([[0.,0.,0.,0.07],
              [0.,0.,-0.38,-0.59],
              [0.,0.,0.38,-0.59]])
A12 = np.array([[0.,0.,0.,0.04],
              [0.,0.,0.61,-0.35],
              [0.,0.,-0.61,-0.35]])
B11 = np.array([[0.,0.,0.,0.],
              [0.,0.,-0.58,0.40],
              [0.,0.,-0.58,-0.40]])
   
# molecular model and coordinates
def molecule(m):
    X, Y, U, V = zip(*m)
    
 #C: atom centre 
    C = plt.Circle((U[0],V[0]), 0.3,
                    color='black', alpha=0.75,
                    hatch='//', fill=False)
    M1 = plt.Circle((U[1],V[1]), 0.2,
                    color='r', alpha=0.75,
                    hatch='//', fill=False)   
    M2 = plt.Circle((U[2],V[2]), 0.2,
                    color='r', alpha=0.75,
                    hatch='//', fill=False)

    bond_y = [U[1], U[0], U[2]]
    bond_z = [V[1], V[0], V[2]]

    return C, M1, M2, bond_y, bond_z

def drawing(m, op):
    fig = plt.figure()
 
    ax = fig.add_subplot(111, aspect='equal')
    ax.set_xlim(-2,2)
    ax.set_ylim(-1.5,1.5)
    ax.set_xlabel('y-axis', size='20')
    ax.set_ylabel('z-axis', size='20')

    C, M1, M2, bond_y, bond_z = molecule(op)
    X1, Y1, U1, V1 = zip(*op)
    X2, Y2, U2, V2 = zip(*m)
    
    n = np.array([[U1[0],V1[0],U2[0],V2[0]],
                  [U1[1],V1[1],U2[1],V2[1]],
                  [U1[2],V1[2],U2[2],V2[2]]])

    mol = ax.quiver(X1, Y1, U1, V1,
                    angles='xy', scale_units='xy', scale=1,
                    label='position vector', color='red',alpha=0.5)
       
    X3, Y3, U3, V3 = zip(*n) 
    vib = ax.quiver(X3, Y3, U3, V3,
                    angles='xy', scale_units='xy', scale=1,
                    label='vibration modes', color='blue')
    
 plt.plot(bond_y, bond_z,'black',alpha=1.0,label='covalent bond')
    ax.add_patch(C)
    ax.add_patch(M1)
    ax.add_patch(M2)

    plt.grid(alpha=0.5)
    plt.legend(loc='upper right')
    plt.title('standard coordinate')

M1 = drawing(A11, stdo) # First 
M2 = drawing(A12, stdo) # Second
M3 = drawing(B11, stdo) # Third

plt.show()

Dengan menginput masukan data berupa koordinat standard stdo dan mode vibrasi A11, A12, dan B11; grafik yang dihasilkan sebagai berikut:

Mode A1(1)

Mode A1(2)

Mode B1

Silakan akses kode menggunakan Google Colaboratory di sini (pemutakhiran: 12/06/2020).

Penulisan kode ini sangatlah banyak kekurangannya. Apabila ada kesempatan, saya ingin mengembangkan algoritma ini menjadi lebih umum terkhusus untuk  molekul anorganik. Tentunya kritik dan saran akan sangat membantu. Semoga bisa bermanfaat.