使用プログラム・ライブラリ（2026年11月12日現在）（ChatGPTに頼んだら、数回のやり取りでさくっと作成・改良してくれました。驚きです。Version 0.1.6。）:
- CREST 3.0.1（配座探索）
- xtb 6.7.1（半経験的量子化学計算）
- ORCA 6.1.1（量子化学計算）
- Python 3.12.11（Webアプリのバックエンド、計算入力の生成、結果解析。標準ライブラリを使用）
- SQLite（ユーザー情報・ジョブ履歴の管理。Python標準の
sqlite3モジュールから利用） - JavaScript/HTML Canvas API（分子構造・スペクトルの描画。外部描画ライブラリを使わず独自実装）
- Slurm 20.11.9（計算ジョブの投入・状態管理）
- Apache HTTP Server（Webアプリの公開・リバースプロキシ）
学部ネットワークからならアクセスできます（URLは内部専用ページに）。学生さんは簡単な分子で遊んでみてください。
使い方: 酢酸エチルを例に
Step 1. 分子構造の読み込み
分子構造ファイル（XYZ, SDF, MOL, PDB, MOL2）を読み込みます。とりあえず試すだけなら、左下の「エタノールの入力例」ボタンを押すとエタノール分子が表示されます。この例では酢酸エチル分子のXYZファイルを読み込みました。3次元構造は、Avogadroなどで作成するか、化合物がPubChemに登録されていれば、3D conformer→DOWNLOAD COORDINATESからSDF形式のファイルをダウンロードして利用できます。
Step 2. 配座探索の計算条件の設定と入力ファイルの作成、計算の投入、 配座探索結果の確認
計算レベル、溶媒モデル、溶媒、残すコンホマーのエネルギー閾値などを設定した後、「入力ファイルを作成」ボタンをクリックしてください。
次に「計算条件・資源設定・入力ファイルを確認しました。」にチェックを入れてから、「Slurmへ投入」ボタンをクリックしてください。
ジョブ・ログ管理セクションで、「Slurm待機中→Slurm実行中→実行終了・結果を確認」となれば配座探索が終了です。
「解析へ追加」ボタンをクリックすると、下部の「CREST配座とエネルギー」セクションで各配座のエネルギーと構造を確認できます。今回は6種類のコンホマーが得られました。
Step 3. ORCAでの構造最適化、振動解析
ジョブ・ログ管理セクションに戻って、「全配座を次へ」ボタンをクリックすると、6個の配座全てを次のDFT最適化と振動解析へ持っていけます。左側の座標に複数の構造が読み込まれていることを確認してください。
計算レベルは以下の条件を選びました:
- 構造最適化: B3LYP/6-31+G(d,p)、気相
- 振動解析: B3LYP/6-311+G(2d,p)、クロロホルム溶媒
（注: ORCAは計算を高速化するため、ハイブリッド汎関数に対してデフォルトでRIJCOSX近似を使用します）。とりあえず遊んでみるなら、B3LYP/6-31G(d) やB3LYP/def2-SVPレベル（気相）で試してみてください。B3LYP/6-31G(d) レベルだと、分子量400程度の低分子のエネルギー一点計算に30〜50秒ほど（18 CPUの並列計算）かかります。2014年のPierensの論文では、NMR予測の前の構造最適化をB3LYP/6-31+G(d,p), 気相で、さらにコンホマー比を計算するためにB3LYP/6-311++G(2d,p) + 溶媒モデルで再構造最適化と振動解析を行っています。
「最適化後に振動解析・自由エネルギー計算を行う」にチェックを入れてから、「入力ファイルを作成」ボタンをクリックして、「計算条件・資源設定・入力ファイルを確認しました。」にチェックを入れてから、「Slurmへ投入」ボタンをクリックしてください。
ジョブの終了を確認したら「全配座を次へ」ボタンをクリックします。
Step 4. ORCAでのECDスペクトル予測
・TDA: B3LYP/6-311+G(2d,p)、クロロホルム溶媒
に設定して、「最適化後の座標を入力し、収束・虚振動・重複を確認しました。」にチェック→「入力ファイルを作成」をクリック→「計算条件・資源設定・入力ファイルを確認しました。」にチェック→「Slurmへ投入」をクリック、で計算が開始されます。
計算が終了したら、「解析へ追加」をクリックして、「スペクトルを表示・平均を計算」でスペクトルが表示されます。黒い太線がボルツマン平均です。酢酸エチルはアキラルな分子なので、エナンチオメリックなコンホマー対が計算に含まれていれば、相殺されてゼロになるはずですが、なっていませんでした。鏡像配座の取りこぼしや重み付けなどを含め、原因は現在確認中です。
Step 5. ORCAでのNMRスペクトル予測
Pierens (2014) で使われていた計算レベルFBS07（mPW1PW91/6-311+G(2d,p)、chloroform溶媒）でNMRスペクトルを予測します。ただし、Pierensの係数はGaussianのIEF-PCM条件で得られたものです。本例のORCA計算にはSMD（The Solvation Model based on Density）を使用して、Pierensの係数を移用しました。
「最適化後の座標を入力し、収束・虚振動・重複を確認しました。」にチェック→「入力ファイルを作成」をクリック→入力確認のチェック →「Slurmへ投入」をクリック→計算終了後に「解析へ追加」をクリック。
次に、スペクトルを表示するオプションを指定します。NMR補正プリセットは「FBS07」、補正係数の溶媒は「chloroform」を選択します。下部の等価な核グループでは、メチレン基を指定する「CH2@C3: 7,8」という行を手動で追加します。補正係数はGaussian／IEF-PCMで求められた文献値を使用しています。今回は計算プログラム・溶媒モデル・構造条件が異なるため、論文と同等の精度を保証するものではありません。
最後に、「スペクトルを表示・平均を計算」ボタンを押すと、NMRスペクトルが下部に表示されます。
NMR予測結果（in CHCl3; 文献値はCDCl3中のresidual solvent peakの値）
酢酸エチルのNMR予測値を、Fulmerらが報告したCDCl₃中の化学シフトと比較しました。計算ではクロロホルムの連続溶媒モデルを使用しています。
Gregory R. Fulmer, Alexander J. M. Miller, Nathaniel H. Sherden, Hugo E. Gottlieb, Abraham Nudelman, Brian M. Stoltz, John E. Bercaw, Karen I. Goldberg; NMR Chemical Shifts of Trace Impurities: Common Laboratory Solvents, Organics, and Gases in Deuterated Solvents Relevant to the Organometallic Chemist. Organometallics 10 May 2010; 29 (9): 2176–2179. https://doi.org/10.1021/om100106e
| 1H | 予測値 δH / ppm (J in Hz) | 文献値 δH / ppm (J in Hz) | ChemDraw 20.0予測値 / ppm (J in Hz) |
|---|---|---|---|
| CH3CO | 2.13 | 2.05 | 2.04 |
| CH2CH3 | 3.99 (6.4) | 4.12 (7) | 4.12 (8) |
| CH2CH3 | 1.34 (6.4) | 1.26 (7) | 1.26 (8) |
| 13C | 予測値 / ppm | 文献値 / ppm | ChemDraw 20.0予測値 / ppm |
|---|---|---|---|
| CH3CO | 20.76 | 21.04 | 20.7 |
| CO | 173.64 | 171.36 | 170.2 |
| CH2 | 61.64 | 60.49 | 61.0 |
| CH3 | 14.17 | 14.19 | 14.1 |
天然物の1Hおよび13C NMRスペクトル予測: Mexicanin I
Mexicanin Iは、キク科の多年草Helenium mexicanumから単離されるPseudoguaiaine骨格を有するセスキテルペンラクトンです。分子式はC15H18O4、分子量262.30 g/mol。Dračínský et al. (2012) とPierens (2014) で、検証に使われていたので、計算してみました。
計算条件:
- 初期構造: PubChemから3D conformerをSDF形式でダウンロード（https://pubchem.ncbi.nlm.nih.gov/compound/Mexicanin-I）
- 配座探索: GFN2-xTB/ALPB(CHCl3)、-ewin 2.5 (kcal/mol)
- 構造最適化、RIJCOSX-B3LYP/6-31+G(d,p)（気相）; 再構造最適化+振動解析、SMD(chloroform)/RIJCOSX-B3LYP/6-311+G(2d,p); NMR予測、気相での最適化構造を使ってSMD(chloroform)/RIJCOSX-mPW1PW91//6-311+G(2d,p)。
CRESTでの配座探索で最安定から2.5 kcal/mol以内に得られたのは6位のOH基の水素原子の向きが異なる3種類のコンホマーだけでした。自由エネルギー差は0.0000, 0.7211, 1.9205 kcal/mol、存在比は74.894, 22.177, 2.929%（298.15 K, 1 atm）。
計算時間（DFT計算には18CPU使用）は、
- 配座探索、6分
- 構造最適化（気相）、各6〜9分
- 再構造最適化 + 振動解析（CHCl3）、各40分程度
- NMR予測、各21分
で合計3時間半程度でした。
| 1H | 予測値δH (SMD/CHCl3) / ppm | Exp. δH (in CDCl3) / ppm | Δδ (予測 − Exp.) | ChemDraw 20.0 δH / ppm | 13C | 予測値δC (SMD/CHCl3) / ppm | Exp. δC (in CDCl3) / ppm | Δδ (予測 − Exp.) | ChemDraw 20.0 δC / ppm |
|---|---|---|---|---|---|---|---|---|---|
| H-1 | 2.77 | 2.68 | 0.09 | 2.00 | C-1 | 54.12 | 52.87 | 1.25 | 53.3 |
| H-2 | 7.92 | 7.67 | 0.25 | 6.57 | C-2 | 169.55 | 162.49 | 7.06 | 162.0 |
| H-3 | 6.00 | 6.15 | -0.15 | 6.07 | C-3 | 128.05 | 130.24 | -2.19 | 129.6 |
| C-4 | 214.78 | 214.35 | 0.43 | 209.0 | |||||
| C-5 | 58.87 | 56.20 | 2.67 | 57.9 | |||||
| H-6 | 4.43 | 4.54 | -0.11 | 3.55 | C-6 | 67.11 | 65.76 | 1.35 | 72.4 |
| H-7 | 3.12 | 3.11 | 0.01 | 2.57 | C-7 | 53.21 | 52.73 | 0.48 | 49.2 |
| H-8 | 4.75 | 4.81 | -0.06 | 4.50 | C-8 | 76.48 | 75.41 | 1.07 | 77.0 |
| H-9α | 1.31 | 1.39 | -0.08 | 1.38 | C-9 | 43.16 | 44.40 | -1.24 | 42.5 |
| H-9β | 2.34 | 2.56 | -0.22 | 1.48 | |||||
| H-10 | 2.15 | 2.22 | -0.07 | 1.43 | C-10 | 29.36 | 27.18 | 2.18 | 26.5 |
| C-11 | 137.34 | 134.63 | 2.71 | 136.8 | |||||
| C-12 | 168.96 | 169.43 | -0.47 | 170.0 | |||||
| H-13 | 6.42 | 6.41 | 0.01 | 6.33 | C-13 | 123.76 | 122.84 | 0.92 | 122.5 |
| H-13′ | 5.72 | 5.68 | 0.04 | 5.81 | |||||
| H-14 | 1.18 | 1.25 | -0.07 | 0.88 | C-14 | 17.69 | 19.75 | -2.06 | 19.2 |
| H-15 | 1.16 | 1.24 | -0.08 | 1.24 | C-15 | 19.36 | 20.28 | -0.92 | 20.6 |
| MAE | 0.095 | 0.474 | 1.800 | 1.793 |
1Hの化学シフトについての決定係数R2 = 0.9978、MAE = 0.095、RMSD = 0.12、13Cについての決定係数R2 = 0.9987、MAE = 1.80、RMSD = 2.40。
Pierens (2014) での、一番近い計算条件はFBS07 (B) です（Gaussian 09で計算）: 3つのコンホマー、構造最適化はB3LYP/6-31+G(d,p)（気相）、NMR予測はPCM(chloroform)/B3LYP/6-311+G(2d,p)、自由エネルギー差を求めるための再構造最適化+振動解析はPCM(chloroform)/B3LYP/6-311++G(2d,p)。Pierens (2014) のTable 3に掲載されているMAE = 0.091 for 1H、1.935 for 13Cでした。今回はRIJCOSX近似を使っている点と再構造最適化/振動解析のための基底関数が違いますが（H原子にdiffuse関数を追加している）、化学シフト予測はかなり良く再現できていそうです。
Mexicanin Iの1Hでは、MAEは本計算0.095 ppm、ChemDraw 20.0で0.474 ppmでした。13CのMAEはそれぞれ1.800、1.793 ppmで、RMSDは2.40、2.60 ppmでした。計算時間が3時間半と数秒ということを考えると、ChemDrawでの13C化学シフト予測のお手軽さと優秀さが際立ちますね（価格は高い）。
天然物の構造決定の論文ではDP4+で、もっともらしさが評価されているのをよく目にします。DP4+は、候補構造の計算NMR化学シフトと実測値との差を、計算レベルに対応する誤差分布で評価し、候補集合内での相対的な確率を推定する統計手法です。DP4+App（https://pypi.org/project/dp4plus-app/ ）はその実装の一つです。
参考文献
- Martin Dračínský, Miloš Buděšínský, Beata Warżajtis, Urszula Rychlewska; Solution and Solid-State Effects on NMR Chemical Shifts in Sesquiterpene Lactones: NMR, X-ray, and Theoretical Methods. J. Phys. Chem. A 2012; 116 (1): 680–688. https://doi.org/10.1021/jp209408b
- Pierens GK. 1H and 13C NMR scaling factors for the calculation of chemical shifts in commonly used solvents using density functional theory. J Comput Chem. 2014;35(18):1388-94. doi: 10.1002/jcc.23638.
（了）