Bader電荷の描画 – Matlantis株式会社

Glycine ジペプチドを例にBadar電荷の描画方法を紹介します。

Glycine ジペプチドの初期構造は smiles_to_atoms で生成し、構造最適化しました。主鎖が直線状に伸びた構造です。

from pfcc_extras.structure.ase_rdkit_converter import smiles_to_atoms

diglycine = smiles_to_atoms('NCC(=O)NCC(=O)O')

準備する

Calculator をセットして電荷を取得します。

from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator

calc_mode, model_version = "PBE_PLUS_D3", "v8.0.0"
estimator = Estimator(calc_mode=calc_mode, model_version=model_version)
calculator = ASECalculator(estimator)

diglycine.calc = calculator

charge = diglycine.get_charges()

描画する

(a) 三次元構造の原子の色を値でマッピング

正負のうち電荷の偏りが最大の値をカラーマップの最大値、最小値とするため、絶対値の最大値を abs_charge_max として取得します。

abs_max_charge = max(abs(charge))

電荷に基づくカラーコードを c_hex として生成します。色付けは負を青、中央の0を白、正を赤で、値の大きさに対応した濃さのグラデーションで表示するために Matplot に組み込みのカラーマップから bwr（blue-white-red）を使います。

Matplotlib に組み込みのカラーマップはこちらをご参照ください。https://matplotlib.org/stable/users/explain/colors/colormaps.html#diverging

from matplotlib import cm, colors
import numpy as np

cmap=cm.bwr_r
rgba = cmap((diglycine.get_charges()/abs_max_charge+1)/2)

c_hex = [colors.rgb2hex(i) for i in rgba]

ビューアー（v）に表示する準備をします。

from pfcc_extras import view_ngl
v = view_ngl(diglycine)
v

この時点では、以下のように描画されます。

別のセルで以下のように実行すると上のビューアーに重ねて、カラーマップによる色付けで原子が色付けられます。

from ase.data import vdw_radii
r_scale = 0.3 # vdW半径を好みで調整
radi = [float(vdw_radii[a.number]) * r_scale for a in diglycine] # vdW半径に基づいた大きさで表示
#radi = np.full(num_atoms, 0.5) # 全ての原子を同じ大きさで表示

v.view.clear_representations()
for i,j,k in zip(diglycine, c_hex, radi):
    v.view.add_representation("spacefill", selection=[i.index], color=j, radius=k)
v.view.add_unitcell()
v.view.add_representation('licorice', radius=0.20, color='parlgray')

先ほどのビューアーが次のように更新されます。

最後に、カラーバーを作成します。

import matplotlib.pyplot as plt
import numpy as np

# カラーバーの大きさ：横に長くしたいなら 4 を大きく、太くしたいなら 0.15 を大きくします
fig, ax = plt.subplots(figsize=(4, 0.15))

# 範囲と色の設定
norm = plt.Normalize(vmin=-abs_max_charge, vmax=abs_max_charge)
sm = plt.cm.ScalarMappable(cmap='bwr', norm=norm)

# 目盛りを5点に設定
ticks = np.linspace(-abs_max_charge, abs_max_charge, 5)
cb = fig.colorbar(sm, cax=ax, orientation='horizontal', ticks=ticks)

# 目盛りとラベルの設定
cb.set_ticklabels([f'{t:+.2f}' for t in ticks]) # 数値はスッキリと
cb.set_label('Partial charge (e)', labelpad=5)   # 必要な単位情報を下部に配置

plt.show()
#plt.savefig('colorbar.png', dpi=300, bbox_inches='tight', transparent=True)　#保存する場合

plt.savefig で transparent=True とすると以下のように背景が透明のカラーバーが作成できます。（保存する場合はコメントアウトを解除して実行してください。）

完成

補足：原子のインデックスと原子種も表示したい場合、別のセルで以下のように実行します。

v.view.add_label(color="black", 
                 labelType="text",
                 labelText=[diglycine[i].symbol + str(i) for i in range(diglycine.get_global_number_of_atoms())],
                 zOffset=1.5, 
                 attachment='middle_center',
                 radius=1.0)

(b) 二次元構造上に値を表示

RDkit を使用して二次元構造上に電荷を表示します。

電荷は準備で定義した charge を使います。

from IPython.display import Image, display
from rdkit.Chem.Draw import rdMolDraw2D
from rdkit import Chem
from rdkit.Chem import AllChem

# SMILES から 2D 構造生成
smiles = "NCC(=O)NCC(=O)O" 
mol = Chem.MolFromSmiles(smiles)
mol = Chem.AddHs(mol)
AllChem.Compute2DCoords(mol)

# 各原子に電荷の値を注釈（Note）としてセット
for i, i_charge in enumerate(charge):
    lbl = f"{i_charge:+.2f}"
    mol.GetAtomWithIdx(i).SetProp("atomNote", lbl)

# キャンバス作成
d2d = rdMolDraw2D.MolDraw2DCairo(600, 400)
opts = d2d.drawOptions()

# index を表示に設定
# opts.addAtomIndices = False
# index の文字サイズを少し小さくする（電荷の数値と区別しやすくするため）
# opts.atomIndexFontSize = 10 

# 描画
d2d.DrawMolecule(mol)
# d2d.DrawMolecules([mol], legends=["Diglycine Bader Charge Distribution"]) #Title をつけたいとき
d2d.FinishDrawing()

# 保存
png_data = d2d.GetDrawingText() # 描画データを変数に格納
with open('diglycine_2d_w_bader-charge.png', 'wb') as f:
    f.write(png_data)

# Jupyterの画面で確認
display(Image(png_data))

(c) 棒グラフとしてプロット

全原子の電荷を棒グラフとしてプロットし、数値を定量的に比較します。