高次元アルゴリズムを用いた分子構造最適化のPCクラスターによる並列処理

寺前 裕之, 大田原 一成


Return

1 はじめに

1970年代末頃からab initio分子軌道法によってエネルギー勾配法が利用できるようになり、電子状態計算によって分子の構造を求めることが広く行われるようになった。このような電子状態計算による分子構造の最適化においては、一般に最適化手法として、Newton-Raphson法、最大傾斜法、共役勾配法などが用いられている[1, 2]。しかし、これらの手法では、初期構造として、エネルギー極小値の近傍、つまりは、求めたい分子構造の近傍から最適化を始める必要がある。また、一つの極小値から他の極小値を求める事は困難であり、局所解を脱することが難しいという課題がある。例えば、初期構造の近くに極小値が存在した場合、しばしばその極小値に収束して最適化が終了するため、大域的にエネルギー最小値を持つ最安定構造を見つけられないことになる。そのため、想定していない構造を見つけることが困難であり、n-C4Me10のような簡単なありふれた分子であってもその構造異性体が長く見つからないで放置されていた[3]。
局所安定構造を脱する方法として、モンテカルロ法やシミュレーテッドアニーリング法などの方法があるが、モンテカルロ法ではエネルギーの高い不安定構造に計算時間の多くを費やし、シミュレーテッドアニーリング法では、最適化の制御パラメータである温度の設定に任意性が残るため、臨界領域における温度のコントロールに困難さがある[4, 5]。これらの大域的な最適化手法を、分子力場法やTight-Binding近似のバンド計算などに適用し、分子構造の最適化に利用した例はあるが、ab initio分子軌道法と組み合わせて用いられた例は無いと思われる[6]。
我々は、極小値近傍の初期構造を必要としない分子構造最適化法として、新上らによる、古典力学系のダイナミクスを利用した高次元アルゴリズム[7]を用い、ab initio分子軌道計算に応用した。高次元アルゴリズムは、力学系の運動を利用した最適化手法であり、大域的な最適化に適した一つの方法である。高次元アルゴリズムは、光ファイバーの多軸合わせ、通信ネットワークのトラフィックの最適化、画像圧縮における量子化テーブルの最適化など種々の工学的問題に応用が行われている[8]。
我々は以前の論文で高次元アルゴリズムによってポテンシャル面を効率良く探る事ができることを示し、一つの極小値から出発しても他の極小値を求めることさえ可能であることを示した[9]。しかしながら一方で、このように大域的に極小値探索を行うためには、最低でも1000回程度の繰り返し計算を必要とするが、毎回のab initio分子軌道計算において、2電子積分、SCF、エネルギー勾配の計算に計算時間の多くが費やされる。従って、化学的に興味がある薬分子やアミノ酸多量体などの大きな分子では、単一CPUでの処理では事実上計算が不可能となってしまう。
一方で、近年のパーソナルコンピュータ(PC)の高速化と低価格化は従来では高価なスーパーコンピューターなどの科学技術計算専用機を使用しないと不可能であった計算も簡便に行うことを可能とした。我々は分子軌道計算部分について複数台のPCを使用した並列処理を行うことで、大きな基底関数を採用する系や薬剤分子のように、分子が対称性を持たずまた構成原子数が多く自由度が大きい系であって、通常の分子構造最適化が困難な場合であっても、高次元アルゴリズムを適用して大域的に分子構造を最適化する事が可能となることをめざした。
並列計算機としては上で述べたように入手が容易でコスト的に有利なPCを複数台使用する。本論文ではPC8台から構成したPCクラスターを使用した。本論文では高次元アルゴリズムのような動的計算の対象としては比較的大きな分子である、ベンゾジアゼピン系の抗不安薬に対して並列計算を試みた。分子軌道およびエネルギー勾配一点あたりの計算時間を計算機の台数を変化させて調べることにより、実際に高次元アルゴリズムを用いた分子構造の最適化の並列処理を行った場合の、台数を増加させた場合の効果や、実際の計算にかかる実時間の短縮効果について検討を行った。

2 計算方法

高次元アルゴリズムによる最適化の過程を説明する。最初に位置xの粒子、ここでは分子内の原子の仮想的な運動を考える。各粒子のポテンシャルとして、最適化する評価関数V(x)、ここでは、ab initio計算で求めるエネルギーに仮想的な運動による運動エネルギーを加えた式(1)を系のエネルギーとする。

ここでは仮想的な運動として、運動の非対角項を許しておくことにする。そのため、各々の粒子の運動は、互いに混合される(ミキシング)ことになる。ミキシングの大きさをパラメータbで指定する。運動方程式は式(2)で与えられ、ここでfはエネルギー勾配法から計算される力である。

式(2)で与えられる運動方程式に従い、ヴェルレ法[10]により、時刻tにおける速度v(t)を求め(式(3))、各粒子の1つ前の位置座標と速度から次の粒子の位置座標を求める(式(4))ことを繰り返すことで、連続的な運動、即ち、分子構造の変化を生成する。連続的な運動を続ける中で、評価関数であるポテンシャルエネルギーの最小値を記録更新し、充分な時間、運動を繰り返した後、最終的に更新された最小値を最適解の近似値とする。この時、運動が混合性を持ち、充分な運動を追跡したなら大域的な最適解が探索できる。ここで、運動の混合性は、空間中に多数の粒子がそれぞれの初期値で始まる別々の運動を運動方程式に従って行っていくと、運動が時間とともに空間中に一様に分布する状態を指す。運動の軌跡を描くことで混合性を簡易的に確認できるが、混合性がない場合は、運動の軌跡は規則性を残し、混合性がある場合は、時間の経過にともない、軌跡は空間中を隅々まで一様に巡ることになる。
高次元アルゴリズムを用いた分子構造最適化のプログラムは、ab initio分子軌道計算プログラムであるプログHONDO5[11]を元に作成した。プログラムは、2電子積分の計算とFockの行列要素の計算、及び、2電子積分の核座標に関する微分の計算について並列処理化を行うとともに、1電子積分の計算についても並列処理化を行った。並列処理プログラムの概要をFigure 1に示す[12]。


Figure 1. Parallelized program of two-electron integrals.

PCクラスターは、CPUにクロック周波数3.0GHzのIntel Pentium4(2次キャッシュ1024KB)を用い、8台で構成した。CPUを搭載した8台のマザーボードには、メインメモリとして各々1.0GBを備えた。また外部記憶装置として各々160GBのハードディスク(Ultra ATA 133, 7200rpm)を備え、各PCは1000BaseTのイーサーネットによりネットワーク接続した。オペレーティングシステムにはLinux kernel 2.6.9 (CentOS 4.4) [13]を使用し、並列処理を行うためのメッセージパッシングライブラリーはMPICH ver.1.2.7p1 [14]を用い、FORTRANコンパイラはIntel FORTRAN ver. 9.1 [15]を使用した。
計算対象とした分子は、種々の置換基により薬効に違いがあり、ドラッグデザインの対象として化学的にも興味が持たれる抗不安薬とした。具体的には、実際に市販されている薬剤分子でFigure 2に示したベンゾチアゾピン(BZP)系のフルトプラゼパム(1, C19H16ClFN2O)、トリアゾラム(2, C17H12Cl2N4)、フルタゾラム(3, C19H18ClFN2O3)、ロラゼパム(4, C15H10Cl2N2O2)、およびチエノジアゼピン(TZP)系のクロチアゼパム(5, C16H15ClN2OS)、エチゾラム(6, C17H16ClN4S)の各分子とした。チエノジアゼンピン系はベンゾジアゼピン系の骨格中のベンゼン環をチオフェン環へ置換したもので、一般的にはベンゾジアゼンピン系と同種の薬剤と見なされている。基底関数には3-21G基底[16]を使用した。計算時間の計測はSCF計算およびエネルギー勾配の一点計算を各10回行い、最も速い値を採用した。


Figure 2. Calculated minor tranquilizer molecules with benzodiazepin and thienodiazepin framework.

3 結果と考察

PCの台数を増加させた場合の並列処理による計算時間に対する効果を調べるため、フルトプラゼパムの計算について、並列処理に用いたCPU数に対する各々のCPU時間および実時間を計測した結果をTable 1に示す。CPU時間は並列化のプログラミング効率の目安となるが、本研究の場合のように計算機を独占して使用している場合には、実際に計算が終了するまでの時間を示す実時間が特に重要である。CPU数が増えるに従って、CPU時間、実時間ともに減少し、並列処理によって計算時間が短縮されていることが確認できる。ここで、CPU数が1の場合の処理時間を基準とし、各CPU数における処理時間との比を加速率として定義して比較する。主に演算に使われるCPU時間として加算されるUser CPU時間は、CPU数4まではほぼ線形に加速率は増加した。加速率はCPU数の90%程度であり、さらにCPU数が増えるに従って加速率の上昇は緩やかになり、CPU数8では加速率がCPU数の73%となった。このようにCPU数が4台程度までは、ほぼ線形の速度向上が認められ、並列化の効率は比較的高いことが示された。

Table 1. The measured CPU time and wall clock time in seconds by parallel processing of one SCF and gradient calculation of flutoprazepam.
# CPUAcceleration ratiowallAcceleration ratio
totalusrsystotalusrsysclockof wall clock
1445.9257.3188.61.001.001.002509.11.00
2231.8137.294.61.921.881.991430.31.75
3158.195.962.22.822.683.031001.02.51
4109.973.936.04.063.485.23588.14.27
580.361.618.75.554.1810.0985.829.26
670.254.216.06.354.7511.7673.234.26
762.348.414.07.155.3213.5266.137.95
856.644.112.67.875.8415.0259.642.09

一方、実時間では、CPU数4まではCPU時間と同様にほぼ線形に加速率が増加するが、CPU数が5を超えると、CPU数5で29.3倍、CPU数8では42.1倍と、CPU数以上の大きな加速率、いわゆるスーパーリニアリティが得られた。また、実時間とCPU時間の差を比較すると、CPU数4までは、実時間がCPU時間の5.3倍程度と、CPU時間に対して実時間が大きくなるのに対し、CPU数5以上では、1.1倍程度でありCPU時間と実時間の差が少ないことがわかる。
このCPU数4と5の間を境にした大幅な実時間の短縮の理由は、CPU数4までは、計算した2電子積分のファイルをハードディスク上に記録してSCF計算の際に繰り返し読み出すために必要な時間が、CPU数5以上では事実上不要となることが考えられる。つまりPC台数の増加により、CPU数が増えるとともにメモリー容量が1台につき1GBずつ増加し、計算に用いる2電子積分のファイルが、全てLinuxオペレーティングシステムにより自動的にメモリー上にバッファーされるので、最初の1回以外はハードディスクの読み出しをする必要がなくなるためと考えられる。
また、CPU時間は、主に演算に使われるuser CPU (usr)時間と、通信パケット処理やファイル入出力処理などに使われるsystem CPU (sys)時間の和として算出されるが、そのうちsys時間がCPU数5以上では大幅に短縮されており、ファイル入出力処理が減少したことによると理解できる。
Table 2に計算した各分子についての並列度に対するCPU時間および実時間の計測結果を示す。フルトプラゼパムと同様に、トリアゾラムは、CPU数3と4の間で実時間の大きな加速率の違いがあり、フルタゾラムは、CPU数6と7の間、ロラゼパムは、CPU数2と3の間、クロチアゼパムは、CPU数2と3の間、エチゾラムは、CPU数3と4の間に実時間の大きな加速率の違いが見られることがわかる。Table 3に各分子の3-21G基底における2電子積分の保存に必要なメモリー容量を示す。トリアゾラムでは、2電子積分の必要容量は3.3GBであり、それを超えるCPU数4以上で実時間の大きな加速率の違いがある。またフルタゾラムでは、2電子積分の必要容量は6.0GBであり、充分なメモリー容量となるCPU数7以上で実時間が大きく加速されている。

Table 2. Comparison of the CPU time and the wall clock time by parallel processing.
<
MoleculeNumber of CPUCPU timeAcceleration ratio of CPU timeWall timeAcceleration ratio of wall time
Flutoprazepam (1)1445.91.002509.11.00
2231.81.921430.31.75
3158.12.821001.02.51
4109.94.06588.14.27
580.35.5585.829.26
670.26.3573.234.26
762.37.1566.137.95
856.67.8759.642.09
Triazolam (2)1614.11.00 1817.31.00
2312.71.961007.31.80
3217.52.82712.62.55
4147.04.18150.412.08
5118.55.18122.414.85
6105.35.83107.716.87
787.77.0090.520.08
878.67.8281.222.38
Flutazolam (3)11220.11.003637.31.00
2620.81.972066.81.76
3437.12.791416.92.57
4320.33.811093.13.33
5254.74.97875.34.16
6217.25.62589.56.17
7171.87.10177.720.47
8154.17.92158.023.02
Lorazepam (4)1454.91.001471.91.00
2231.21.97816.31.80
3146.93.10149.49.85
4108.24.20110.313.34
587.75.1990.216.31
677.65.8680.118.38
764.37.0767.421.84
858.27.8160.524.33
Clotiazepam (5)1567.51.001715.41.00
2287.31.981015.61.69
3179.13.17191.08.98
4134.24.23137.212.51
5107.85.26110.715.49
694.56.0196.917.71
782.66.8785.320.11
872.27.8674.622.98
Etizolam (6)1663.61.001785.81.00
2337.51.97987.21.81
3238.22.79705.42.53
4160.24.14164.110.88
5129.95.11133.213.41
6115.55.75117.715.17
795.16.9898.518.13
885.07.8187.520.40

Table 3. The amount of memory requirement to store two-electron integrals.
MoleculeNumber of basis functionsAmount of two-electron integrals [Gbyte]
1Flutoprazepam2484.2
2Triazolam2393.3
3Flutazolam2746.0
4Lorazepam2172.5
5Clotiazepam2272.9
6Etizolam2453.7

このように2電子積分の計算に必要なメモリー容量に対して、システム全体のメモリー容量が充分に大きくなるCPU数となった時点で、実時間の大きな加速率の違いが生じることがわかる。このような効果はメモリー共有型ではない複数台の計算機を用いた並列化において特徴的なことである。対象とする分子軌道計算の2電子積分のファイルがメモリー上に保存可能となるまで台数を増加させるだけで大幅な実時間の短縮が実現すると考えられる。実時間を短縮するのに、プログラムを書き直す必要が無いのは大きな利点である。
一般の計算対象に対して、どの程度の台数を用意すれば良いのであろうか。効率良く構造最適化計算するのに必要なPCクラスターのシステム規模を知るため、上で述べた2電子積分ファイルのメモリーによるバッファー効果で高速化を実現するのに必要となるPCクラスターのCPU数をモデル計算により見積もることを試みた。モデル計算の対象としてポリアミノ酸であるグリシン多量体を選び、3-21G、6-31G、6-31G*および6-31G**基底での2電子積分数、2電子積分容量を計算した結果をTable 4に示す。グリシン多量体の構造は平面構造とし分子モデリングソフトであるChem3D[17]上で3次元構造を作成し、計算に用いる座標を生成した。

Table 4. Number of two-electron integrals and the amount of memory requirement on calculation of glycine oligomer.
Basis setNumber of glycine unitsNumber of atomsNumber of basis functionsNumber of two-electron integralsAmount of two-electelectron integrals [Gbyte]
3-21G110551,134,8870.017
217978,912,4020.133
32413927,134,6060.404
43118155,662,7780.829
53822393,693,0221.396
10108433421,979,9216.280
15143643970,985,02314.469
201788531,730,850,00425.792
6-31G110551,166,2400.017
217979,677,0550.144
32413930,218,8260.450
43118162,344,1150.929
538223104,989,7481.565
10108433469,559,9756.997
151436431,074,368,60216.009
201788531,907,382,33728.422
6-31G*110856,566,1480.098
21715153,738,9960.801
324217158,833,1372.367
431283316,826,4454.721
538349521,664,9477.773
101086792,226,896,41033.183
6-31G**11010012,338,2250.184
21717589,982,4881.341
324250255,164,6043.802
431325499,632,7397.445
538400816,127,89612.161
101087753,430,410,26351.117

グリシン数が1から10に増えるに従って、急速に2電子積分の量が増加しているのがわかる。
通常SCF計算では、2電子積分を一度計算した後、メモリー上またはハードディスクに保存し、繰り返し計算の中で読み出して再利用する。大きな分子を計算する場合、メモリー・ハードディスク容量の制限やディスク入出力による速度低下を防ぐため、2電子積分を保存せずに毎回計算するdirect SCF法が使われることが一般的である。しかしながら、並列処理においては、前述のとおり台数を確保することでメモリー容量を確保することができ、ほとんどの場合、2電子積分を保存する通常SCF計算を高速に処理することが可能となると考えられる。
グリシン10量体の場合、3-21Gでは433基底、6-31G**では775基底であり、各々の積分容量は、6.3GBおよび51.1GBである。これから2電子積分ファイルの読み書きについて、メモリー上へのバッファー効果が期待できるPCクラスター台数は、PC毎のメモリーを1GBとすると、各々7台、52台となる。なおTable 4より2電子積分数は、基底関数数Nの約3.1〜3.3乗に比例していることがわかるため、これを用いてあらかじめ必要台数を見積もることも可能となる。例えば3.3乗に比例すると仮定すると、433基底では7.5GB、775基底では51.0GBとなるため、各8台、52台となる。
高次元アルゴリズムによる大域的な構造最適化計算では、概ね1000回程度の繰り返し計算が必要であることから、例えばフルトプラゼパムの場合、並列化しないCPU数1では、仮に充分なメモリーを搭載し、CPU時間だけで処理できたとしても計算時間に 257.3秒×1000回 = 257300秒、71時間で、約3日間が必要であると見積もることができる。もしマザーボード上に搭載できるメモリーの制限などにより、充分なメモリーを準備できない場合には、約29日間が必要と見積もられる。今回構成したPCクラスターシステムでは、同じ計算についてCPU数8の場合で16.5時間と見積もることができる。以上のように大幅な計算時間の短縮が期待でき、並列処理の効果が特に大きいことがわかる。

4 結論

Ab initio分子軌道計算の構造最適化手法として、大域的に最適値を探索する高次元アルゴリズムによる計算を現実的な時間内で処理するために、並列計算機として、Pentium4 (3.0GHz)を搭載した8台のパーソナルコンピュータから構成したPCクラスターで、1電子積分、2電子積分、SCF計算、エネルギー勾配の計算に対する並列処理を行った。計算対象には、化学的に興味の持たれる分子として、ベンゾチアゾピン系およびチエノジアゼピン系の抗不安薬を選び、3-21G基底を用いたSCF計算とエネルギー勾配計算について、PCの台数を増加させた場合の並列処理時間を調べ、実時間の短縮効果を検討した。
例えばフルトプラゼパムの計算では、CPU数1に対し、CPU数4で、CPU時間4.1倍、実時間4.3倍、CPU数8で、CPU時間7.9倍、実時間42.1倍を得たが、他の分子においても同様な結果が得られた。PC台数を増加させてCPU数を増やすことで、同時に各マザーボード上に搭載するメモリーが有効に利用できることがわかった。PCクラスターのメモリー上に必要な2電子積分ファイルをバッファーできるため、ハードディスクのI/O時間が事実上無視できるようになり、実時間の大幅な短縮がなされる。
このように、ab initio分子軌道計算のPCクラスターによる並列処理で高次元アルゴリズムを使用した構造最適化計算の実時間短縮が大幅に計れることが示された。今後はこれらのPCクラスターを用いた計算により、本論文で示した、抗不安薬の構造最適化計算を行っていく予定である[18]。

参考文献

[ 1] A. R. Learch, Molecular Modelling: Principles and Applications, Longman, England (1996).
[ 2] E. Polak, Computational Methods in Optimization, Academic Press, New York (1971).
[ 3] F. Neumann, H. Teramae, J. W. Downing, J. Michl, J. Am. Chem. Soc., 120, 573 (1998).
[ 4] H. Cohn, M. J. Fielding, SIAM J. Opt., 9, 779 (1999).
[ 5] S. Kirkpatrick, C.D. Gelett, M. P. Vecchi, Science, 220, 671 (1983).
[ 6] M. R. Lemes, C. R. Zacharias, A. Dal Pino, Phys. Rev., B56, 9279 (1997).
[ 7] K. Shinjo, T. Sasada, Phys. Rev., E54, 4686 (1996).
[ 8] K. Shinjo, J. Jpn. Soc. Fuzzy Theo. and Sys., 11, 382 (1999).
[ 9] K. Ohtawara, H. Teramae, Chem. Phys. Lett., 390, 84 (2004).
[10] L. Verlet, Phys. Rev., 159, 98 (1967).
[11] M. Dupuis, J. Rys, and H. F. King, J. Chem. Phys., 65, 111 (1976).
[12] H. Teramae, J. Chem. Software, 4, 73 (1998).
[13] http://www.centos.org/
[14] http://www-unix.mcs.anl.gov/mpi/mpich/
[15] http://www.intel.com/
[16] J. S. Binkley, J. A. Pople, W. J. Hehre, J. Am. Chem. Soc., 102, 939 (1980).
[17] 10 CambridgeSoft Corporation, 100 CambridgePark Drive, Cambridge, MA 02140 USA.
[18] H. Teramae, K. Ohtawara, T. Ishimoto, U. Nagashima, Bull. Chem. Soc. Jpn., 81, 1094 (2008).


Return