薬の設計のための量子化学シミュレーション

IBMが開催しているQuantum Challengeという量子コンピューターのプログラミングコンテストがあります。2019年からすでに５、６回開催されていますが、2021年の秋に、**[IBM Quantum Challenge Africa 2021](https://github.com/qiskit-community/ibm-quantum-challenge-africa-2021)** というアフリカ大陸の社会問題を取り扱った大会がありました。１問目は最適化、２問目は金融の問題、３問目は量子化学となっており、量子コンピューターの応用事例としておもしろい例題が並んでいました。
今回、３問目を翻訳したので、概要を紹介したいと思います。（3問目の全編翻訳は[こちら](https://github.com/quantum-tokyo/introduction/blob/main/materials/QuantumChallengeAfrica2021/lab3/lab3_ja.ipynb)です。１問目は、[こちら](https://github.com/purepureclub/IFCO2021DEC/blob/main/IFCO2021Dec_qiskit_handson.ipynb)にリョウコさんによる日本語解説があります。）
化学についてすっかり忘れてしまった方にも、原子や分子について思い出しながら読んでいただけるような内容になるように心がけました！

# 目標：量子コンピューターを使ってHIVの薬の効果を確認する
皆さんもご存知の通り、HIVは公衆衛生上、大きな課題となっているウイルスです。このHIVの薬の効果をVQE（変分量子固有値ソルバー）を使って調べるのが、このチャレンジの課題です。(ただし実際の計算にはトイモデルを使います。)
量子コンピューターによって、これまでシミュレーションの難しかった医薬品の設計が飛躍的に向上する可能性があります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151117/6e5d83f5-9b25-9746-4742-f83010ace2b5.png)

HIV（ヒト免疫不全ウイルス、Human Immunodeficiency Virus）は、人の免疫細胞に感染して破壊し、最終的に後天性免疫不全症候群 (AIDS)を発症させるウイルスです。増えていく過程で、ウイルスの組み立てに必要なタンパク質を、感染した細胞に作らせますが、大きなタンパク質のままではウイルスが組み立てられません。
HIVに特有の「プロテアーゼ」という酵素によって、そのタンパク質が適当な大きさに切断され、増殖されていきます。この「プロテアーゼ」が今回のシミュレーションの対象となる分子その１で、いわゆる悪者の分子です。
これに対し、2つ目のシミュレーションの対象となる分子が「抗レトロウイルス薬」で、プロテアーゼに結合してブロックする薬剤です。まとめると、

シミュレーションの対象となる分子：
１）「プロテアーゼ」：ウイルスタンパクをハサミでカットしてウイルスを増やす悪者
２）「抗レトロウイルス薬」：プロテアーゼにくっついてハサミの能力をブロックする薬

# VQE概要
今回使う、VQE (Variational Quantum Eigensolver、変分量子固有値ソルバー）アルゴリズムは、量子化学の分子エネルギー固有値問題において使うことができる、量子・古典のハイブリッドアルゴリズムです。特にNISQ(ノイズあり中規模量子コンピューター)向けの量子アルゴリズムとして活用が期待されていて、以下のような手順で計算を行います。
1. 対象とする分子のエネルギーに相当するハミルトニアン$H$(演算子)を準備する。
1. パラメーター$\theta$の入った試行状態(ansatz)$|\psi(\theta)\rangle$を作る。
1. $|\psi(\theta)\rangle$に$H$演算をして測定し、エネルギー期待値$\langle H(\theta)\rangle$ を求める。
1. 古典コンピューターで$\theta$を調整して 2.、3.を繰り返し、エネルギーが最小になるように最適化する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151117/1619b6c5-7aeb-ecef-6515-5ffa425ddd70.png)

図の出典：Alberto Peruzzo,他 ["A variational eigenvalue solver on a quantum processor"](https://arxiv.org/pdf/1304.3061.pdf)
具体的な計算方法は、この後に水素分子の例が載っています。さらに詳しい解説は、[QiskitテキストブックのVQEの章](https://qiskit.org/textbook/ja/ch-applications/vqe-molecules.html)をご覧ください。

# 化学の用語
原子の一般的な構造の模式図は、[wikipedia](https://ja.wikipedia.org/wiki/%E5%8E%9F%E5%AD%90)によると以下のようになっています。化学反応の決め手となるのは、原子の外側にある電子の反応であるため、電子軌道について調べることが重要となります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151117/6f24e891-024b-0aec-ba27-f2a31a309641.png)

上の図は模式図であり、もう少し詳しく電子軌道を見てみると、[原子軌道のWikipedia](https://ja.wikipedia.org/wiki/%E5%8E%9F%E5%AD%90%E8%BB%8C%E9%81%93)に掲載されている以下の図からも見て取れるように、原子核の周りに存在する電子の波動関数が電子軌道であり、大雑把にいうと電子の存在確率が高い場所が電子軌道になります。（絶対値の二乗が存在確率に比例。）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151117/1c89ad59-c47f-ad7a-4599-bdad208d3ea7.png)
（左から順に1s、2s、2px、2py、2pz軌道。）

ここで、量子コンピューターにおける量子化学分野での頻出用語をまとめておくと、
- 電子軌道：電子の波動関数（電子の存在確率が高い場所）。一つの軌道には↑スピンと↓スピンの2個の電子が入る。
- 原子軌道：原子核のまわりの電子軌道。
- 分子軌道：分子（原子が複数集まったもの）の電子軌道。
- スピン軌道：電子の個数分ある(電子軌道の２倍)スピンの軌道。(量子化学計算で使われるJordan Wigner変換ではこの数だけ量子ビットが必要になる。)
- アルファ電子：↑(+)スピンの電子。
- ベータ電子：↓(-)スピンの電子。

原子軌道や分子軌道といわれると原子や分子が回っている軌道のことかと勘違いするかもしれませんが、原子軌道、分子軌道とも電子の軌道のことを指しています。

[元素の周期表](https://ja.wikipedia.org/wiki/%E5%91%A8%E6%9C%9F%E8%A1%A8)では、1から順番に番号が振られていますが、それぞれの元素は、その原子番号の数だけ電子を持っています。例えば、原子番号1の水素原子の電子の数は1個、原子番号6の炭素原子の電子の数は6個です。
また、基底状態の電子は入る軌道が決まっていて、エネルギーの低い順に順番に席が埋まっていきます。下の図のように、↑または↓で表された電子が、1sから順番に2個ずつペアで配置されていきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151117/043dcb7c-0f68-cee5-3afd-5c2b75c114d5.png)


# 例題：水素分子の基底エネルギーを求める
目標のHIVの薬の効果を調べる前に、[本文(和訳版)](https://github.com/quantum-tokyo/introduction/blob/main/materials/QuantumChallengeAfrica2021/lab3/lab3_ja.ipynb)では、水素分子の基底エネルギーを例に、VQEの手順を解説しています。以下のようなステップで行います。(ステップの番号は本文の番号と同じです。)

1  ：ドライバーで、原子の位置を指定して、古典計算から必要な情報を取得。
2a：取得した古典情報からフェルミ演算子のハミルトニアンを作成。(Hartree-Fock近似(電子1個ずつの運動に近似)。電子(フェルミオン)の存在する場所を量子化(第二量子化)。)
2b：ハミルトニアンを量子ビット(パウリ演算子)に変換。(Jordan Wigner変換では、量子ビットをスピン軌道として表現：$|0\rangle$、 $|1\rangle$ がスピン軌道の占有,非占有に相当)
3a：初期状態をセット。
3a：試行状態（パラメーター付き量子状態、Ansatz）をセット。
3b：量子計算の実行環境を設定。
3c：(古典オプティマイザーを設定。今回は設定していませんでした。)
4  ：VQEを実行し、最適なパラメーターを探し、エネルギーの近似解を求める。

# HIVチャレンジ：薬の効果を調べる
さて、とうとうHIVの「プロテアーゼ(悪者)」と「抗レトロウイルス薬」の反応を調べます。基本の考え方は、二つの分子が結合すると全体のエネルギーが下がるという自然の現象を使います。具体的には、「プロテアーゼ(悪者)」分子と「抗レトロウイルス薬」分子の位置を近づけていって、その位置ごとにVQEでエネルギーを計算し、[エネルギー曲面](https://ja.wikipedia.org/wiki/%E3%83%9D%E3%83%86%E3%83%B3%E3%82%B7%E3%83%A3%E3%83%AB%E3%82%A8%E3%83%8D%E3%83%AB%E3%82%AE%E3%83%BC%E6%9B%B2%E9%9D%A2)を作り、局所的にエネルギーが低いところがあるかどうか確認します。ある一定の距離内でエネルギーの最小値が見つかれば、二つの分子はその位置で結合した状態が安定であり、薬の効果が現れることになります。
以下は、水素分子(水素原子が2個くっついて水素分子になる)の2つの原子間距離を変えたときのエネルギー曲面の例です。エネルギーが最小となる原子核の位置が安定な水素分子の構造です。(水素の基底状態の構造はよく知られているため今回は既知の原子核の位置を0に設定したので0が最小値になっています。)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151117/d2b583ec-f5bb-3488-c7f5-86194014e812.png)

ここで、「プロテアーゼ(悪者)」と「抗レトロウイルス薬」について、本物の分子は大きすぎてそのまま計算にかけることはできないため、トイモデルを使います。

# トイモデルの設定
本物の「プロテアーゼ(悪者)」は、約100個のアミノ酸からなる2本のポリペプチド鎖で構成されています。ペプチド鎖とは、ペプチド結合したアミノ酸がたくさんつながったものです。アミノ酸とペプチド結合について振り返っておくと、下の図のように
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151117/26b2c8b9-a2c3-d5a7-2473-55a6d08191cb.png)

アミノ酸とはアミノ基($-NH_2$)とカルボキシル基($-COOH$)を持つ有機化合物のことです。2つのアミノ酸がペプチド結合すると、ジペプチドと水($H_2O$)ができます。「プロテアーゼ(悪者)」は、このペプチド結合したアミノ酸が100個くらいつながった鎖が2本で作られているとのことです。
ペプチド結合は、一般的なタンパク質の折りたたみやHIVの「プロテアーゼ(悪者)」のハサミの切断能力などタンパク質の化学的な性質を決定する最も重要な要素の一つです。このペプチド結合に注目して、トイモデルを設定します。「プロテアーゼ(悪者)」のトイモデルを、ペプチド結合の$O=C-N$に水素原子($H$)を足した$HCONH_2$ (ホルムアミド)にします。ホルムアミドは、実際にイオン性の溶媒でもあるので、イオン結合を「切る」ことができる物質でもあるとのことです。以下の図がそのモデルです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151117/59d2601d-8624-8f72-ff39-fe77b8bfded8.png)


また、「抗レトロウイルス薬」として、今回は、1個の炭素原子($C$)を使います。
この２つのトイモデルを使って、窒素($N$)と炭素($C$)の間の距離を変えてエネルギーが局所的に最小となるところがあるか調べます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151117/9e39b039-86ef-a717-b153-5e50f25b546e.png)

# 具体的な計算ステップ
以下は、実際のHIVチャレンジの手順と設問をまとめたものです。上から順に解いていくことで、トイモデルを使って薬の効果があるかどうか調べ、また、最後は独自の試行状態(Ansatz)を作り、最も精度の良いシミュレーションモデルを作成するオープンチャレンジになっています。

1：分子の形を変化させる設定
- $N$と$C$の原子間距離を変化させます。
- **Ex 3a：原子の座標と動かす原子を設定する問題。**

2：量子ビット数を削減
- そのままでは60量子ビット必要になり、ビット数が多すぎて実行できないため、Active Space Transformerで内殻電子をただのポテンシャルとみなし、反応に必要な外殻電子のみ計算で扱う設定をします。
- **Ex 3b：分子構造についての化学的な質問。**
    - Q1: 内側の電子より外側の電子を先に量子処理する方が高分子の計算において、経験的に意味を持つのはなぜですか？ 
    - Q2: 量子で扱う電子の数が一定の場合、それらの電子がアクセスできる軌道の数を増やすと、計算される基底状態のエネルギーは、上からも下からも漸近的なエネルギーに近づくのでしょうか？

3：VQEを実行し、原子間距離を変化させた際のエネルギー曲面を表示
- 結合するかどうかを全体のエネルギーが最小となる状態があるかどうかで判定。
- **Ex 3c：薬の効果があるか調べる問題。**
    - 有限の距離の間に明確な最小値はありますか？結合は起こりますか？ 

4：試行状態（パラメーター付き量子状態）の調整。
- **Ex 3d：パラメーターの繰り返し数と計算精度の関係について調べる問題。**
- **Ex 3e：$HCONH_2-C$の分子について、8量子ビットまでのVQEで、誤差が最も小さくなるハードウェア由来の試行状態を見つける問題。**
   - スコアは以下で計算されます：
$$\text{score} = -10 \times \log_{10}{\left(\left\lvert{\frac{E_{true} - E_{computed}}{E_{true}}}\right\rvert\right)}$$


# まとめ
IBM Quantum Challenge Africa 2021の3問目は、HIVのトイモデルを使って薬の効果をVQEで調べる問題でした。まだモデルレベルでの計算ですが、量子コンピューターの具体的な応用例の一つを学べる良い例題になっていると思います。
ぜひ皆さんも、薬の効果があるのかどうか、実際にVQEを実行して調べてみてください！