オブジェクト指向設計による分子動力学計算プログラム

オブジェク ト指向設計 による分子動力学計算 プログラム

石 井

晃 (平成 6年 6月21日受理)

1.序

最近, ソフ トウェア開発の現場 に終 いてオブジェク ト指 向 とい う考えに基づいたプログラム開発 が行 われ る ようになって きている。 これ は従 来 の構造化 プログラ ミングと違 って ソフ トウェアで シュ ミレー トす る対象 をまず「 オブジェク ト」 とい う構成要素 に分解 し

,こ

のオブジェ ク トについ てデー タとそのデー タを処理す る関数 を付加 して一つの まとまったプログラム単位 にす る方法であ る。 しか し

,昔

か ら計算機 を多用 して来 た科学技術分野ではこのオブジェク ト指 向 に基づ いたプロ グラム開発 の理解及 び実行 は十分 とはいえない。それはオブジェク ト指向 を標榜 したプログラム開 発環境のほ とん どが

,主

に大量 の文書 や顧客 リス トな どの処理 を念頭 に置いて整備 されて きている ことか らもわかる。 しか しなが ら

,元

々オブジェク ト指 向 プログラ ミングが構造化 プログラ ミングを批判す る形でで て きたのは

,大

規模 プログラムの開発 を容易 にす る とい う理 由であ り

,そ

のィ点では特 に対象分野 を 限 る もの で は ない。批 判 の対 象 に され た の は

FORTRANに

附 い て言 え ば

COMMON文

やsub_ rouineの _{引数 とい う形でのデー タのや り取 りは甚 だデバ ックが しに くい とい う点 であ る。 た とえ} ば対象 とす る問題 自体 をい くつかの部分 に分割 し

,そ

れ をそれぞれサ ブルーチ ンとい う形 で分割管 理す るのは可能であるが,メ イ ンプログラムの中で,それ らのサ ブルーチ ンの引数 はむ き出 しになっ て しまい

,従

って どの部分で トラブルが起 こって もメイ ンプログラムを経由 してその影響 は全 プロ グラム に及ぶ。 それ を防 ぐため に

COMMON文

の中で変数 を受 け渡 しすれば

,一

応 その変数 を操 作で きるプログラムは局所 的 に限定 され るが

,

この場合 は どの プログラムが どの

COMMON文

を 使用 しているのか を掴 むのが甚 だ難 しくなる。つ ま り変数の受 け渡 しは地下 に潜 って しまい,エデ ィ タの検索機能で探 し回ることになる。 さらに

,構

造化 プログラミングは現場 での大規模 なソフ トウェア開発 に於 いて致命的 な欠陥があ る。それは通常 開発現場では一つの ソフ トウェアが複数の人間に よつて開発 される とい うことであ る。従 って

,一

人の人間が書 いたプログラムに別 の人間が コー ドをつけ加 えて機 能 を拡 張す ること が 日常茶飯事 であることを念頭 に置 く必要があ る。 これは特 にソフ トウェア会社 だけの問題 ではな

石 井 晃

く

,大

学の研究室においても毎年新 しい学生が事業論文制作

.修

士論文制作

,あ

るいは博士論文の

準備として入れ替わ り立ち替わ り開発プロジェク トに携わるのであるから

,大

学研究室におけるソ

フトウェアの開発の状況はソフトウェア開発会社よりさらに意思の疎通という意味では厳

.し

いとい

える。そうした常に入れ替わる複数の人間が一つのプログラムをいじっていくとすると

,た

とえば

COMMON文

を用いてどのプログラムからどこヘデータが流れていったのか

,掴

むのは非常に難 し

くなる。

例を挙げよう。たとえば次のようなプログラムがあったとする

.。 program main common/A/文 ,y,z subroutine alpha(tl subroutine beta ttl end subroutine alpha(tl

common/A/x,y,z

Feturn― end subroutine beta(1) co■

lmon/A/x,y,2

Fetuen end

ここにあげたプログラムでは

,Aと

いう

commOn block

で使われている。ところが

,こ

こで別の学生が

,次

のよ

subroutine crash(t) common/A/x,y,多

x=x+1

t==t return end

が

talpha,betaと

_{いう二つのサブルーチン}

うなプログラムを作ったとしよう。

オブジェク ト指向設計 による分子動力学計算 プログラム

このプログラムは基礎 デー タであるcommOn block Aの 中身 を勝手 に書 き換 える ものであ るか ら,

その学生 が この プログラム を全体 の プログラムの どこに挟 むのかで

,先

の alpha,betaと い うサ ブ ルーチ ンの実行結果 は全 く保証 され な くなって しまう。 この学生 に悪気 はなかったとはいえ

,こ

の プログラムを全体 のプログラムの中の随所 にいれて しまった と した ら

,第

三者が この ことに気づ く のはかな り難 しい。 つ ま り

,従

来の構造化 プログラ ミングの欠陥 とは

,こ

の ような不完全 なサ ブルーチ ンを挟み込 む ことによるデー タの破壊 がい とも簡単 に行 われ ることにある。特 に大量のエラーメッセージを抱 え て, ともか くも動 かそ うと妙 に悪 あが きした ときに

,つ

い に犯 しが ちなことで もある。 オブジェク ト指向設計 はこう したこ とを未然 に防 ごうとい うプログラム設計 の思想 である。従 って科学技術系 の研究室では研 究室で大切 に開発 しているプログラムを保護す るため に も

,是

非必要な設計法 とい える。 そ こで

,本

論 の 目的は

,そ

う したオブジェク ト指 向 によるプログラムの設計・ 開発 をいかに科学 技術

,特

に物理科学 の分野 に導入す るのか を分子動力学 のプログラムを例 に とって考察することに ある。

2.設

計方針

2.1

オブジ ェク ト指向 オブジェク ト指 向設計 と構造化 プログラ ミングとは必ず しも相反す る設計法でないことをまず明 言 してお こう。構造化 プログラ ミングが一人の プログラマーで十分 目が届 く範 囲

,た

とえば千行 か ら二千行 のプログラムを書 く場合 には良 く機能す るのに対 し

,オ

ブジェク ト指 向設計 は とうてい一 人では見渡せ ない範 囲,たとえば数万か ら数十万行 のプログラムの コーデ ィングで威力 を発揮す る。 オブジェ ク ト指 向 は プ ロ グラム設計 の思想 で あ るので

,

どの プ ログ ラム言語 で も

,た

とえ ば

FORTRAN77で

も可能であ る。 しか し

,先

ほ どの ような意 図 しないデー タ書 き換 えの ような もの を防 ぐ意味ではそれな りのデー タ保護機能の備 わったプログラム言語であ るこ とが望 ま しい。 オブ ジェク ト指向言語 を標榜 す る ものはふつ うオブジェク トに相 当す る もの と して クラス とい うもの を 定義 し

,こ

の クラスが

,い

わば変数 とその変数の取 り扱 いに関す る関数の集合体 となる。それぞれ の クラス内の局所 変数 を定義 す る こ とによって クラス内で はい ちいち関数 の引数や

COMMON文

などを使 わな くて も変数 を引用 す ることがで きる。 また

,自

分 の クラス以外 の プログラムか ら参照 で きる変数 とそ うで ない もの とに分 けることがで きるため

,先

に述べ た意図せ ぬデー タの破壊 をか ヽミなり防 ぐことがで きる。 なぜ なら, クラスの外 か らは参照 で きないデー タを参照す るため にはその 対象 となるクラス 自身の コーデ イングを変更す る必要があ り

,そ

れ を実行 す るのは明 らかに「意 図

石 井 晃 した破壊」 しかあ り得 ないか らである。 この ようなデー タ参照 の具護 を情報 隠蔽 と呼 び

,多

くの オ ブジェク ト指向言語 に備 わった機能である。つ ま り

,情

報 隠蔽 は研究室の学生 に

,自

分がい じつて いい部分 といけない部分 を明確 に区分 させ

,

もしその学生 のプロジェク トが失敗 して もいつで も元 のプログラムは 自動 的 に (バックア ップ とい うことでな くて も

)守

られるこ とになる。 さらに

,オ

ブジェク トとしてなるべ く対象 とす る現実 に合 わせ てオブジェク ト化 してい くと

,プ

ログラムが読みやす くなる とい うメ リッ トもある。た とえば本論文で対象 とす る分子動力学 なら, 一つの原子 を一つの クラス として表現す る とい うのはか な りわか りやすい クラス化 であろう。その クラスのデータメ ンバ として波動 関数があ り

,ま

た

,そ

の クラスの関数 と して原子 に働 く力 を計算 す る関数がある

,

とい う具合である。 オブジェク ト指 向言語 の重要 な特徴 として

,派

生 クラス とい う概念がある。 これはあるクラス に 一部機能 をつけ加 えた り

,

また

,変

数 を追加 した りして

,そ

のスーパ ーセ ッ トを派生 クラス として 新 たに定義す ることである。 た とえば

,学

生 に研究室で開発 しているプログラムの クラス 自身はい じらせず, もっぱ らその派生 クラスについてコーデ イングさせ ていれば

,元

のプログラムを全 くい じ らず に い ろ い ろ プ ロ グ ラ ム 開 発 に伴 う テ ス トコー デ イ ン グ が で き る。 も し

,こ

れ が

FORTRAN77で

あ れ ば

,テ

ス ト と し て 動 作 させ る た め に は 今 ま で に 作 ら れ て い る サ ブ ル ー チ ン に つ い て もつ1数 を追 力日した り

,COMMON BLOCKを

追 力日した り しな け れ ば な らず

,バ

ッ ク ア ッ プ フ ァ イル以外で加 えた変更点 を元 に戻 すのはかな り大変 な作業 となる。 しか も

,こ

の機能の追加 の作業 が何週 間 も続 くとなる と

,同

じ名前 のサ ブルーチ ンで引数や

COMMON文

の異 なるバー ジ ョンが 研究室内に同時 に存在す ることとにな り

,そ

のバージ ョン管理 は きわめてや っかいなことになる。 派生 クラス を用 いる とい うや り方 を取 る と

,元

になっているクラスはその ままで派生 クラスにのみ 変更 を加 えるので

,そ

の派生 クラスにその学生 の名前 にちなんだ名前 で も仮 につ けさせておけば, バージ ョン管理の心配 な どは全 く無用 で

,か

な り思 い切 った機 能変更の試 みな ども安心 してで きる ことになる。

2.2

使用言語 オブジェク ト指 向のプ ログラ ミング言語 としてはSmalltalk,C十+,Effelな どいるいろな ものが存 在 してい る。最近処 理系 が出回 り始 めた Fortran90な ども

FORTRAN77と

比べ てか な リオブジェ ク ト指向 を意識 した仕様変更 を している。分子動力学 の ように完成 まで に長期 の開発作業が必要 に なるプロジェク トの場合 は開発途 中で変更で きない プログラ ミング言語の選択 はかな り慎重 に行 わ なければな らない。従 来 か らあ る財 産 を生かす とい う立場 では

FORTRAN系

統 の処理系 の選択 も 十分 に考 え られる。 しか し

,こ

こでの分子動力学 プログラム開発 は

1か

らの開発 であるため

,逆

に 過去の遺産 を気 にせず に最 良の選択 をす ることがで きる。 そ こで

,

こうした処理系の中で

,

ここで

オブジェク ト指向設計 による分子動力学計算 プログラム

151

は

,C十

+を

選択す ることにす る。

C十十は

C言

語の上位互換 にあた り

,C言

語 に関す る過去の遺産 を継承 で きる利点がある。さらに,

FORTRAN同

様 に複素数 を実数 同様 にサポー トす るな ど

,C言

語 に比べ る と科学技術計算 に関す る 機能 は大 幅 に強化 されてい る。 また,M A Ellis and B,Straustrapの “The Annotated C十 +Refer― ence Manu』

"[1]と

ANSIに

よる基準 があ るため

,そ

の範 囲内で プログラ ミングを していれば機 種依存性 は考 えな くていい。 現 在 の コ ン ピュー タ使 用 で は何 らか の ウ イ ン ドウ シス テ ム

(X window,MS Windows,OS/2,

Machintosh OS,etc)の_{使 用 を念頭 に置 い た方 がいいが}

_,こ

_{のユ ーザー イ ンター フェイス と しての} ウイ ン ドウ自体 がC十十 で プロ グラ ミングされ てい るため

,分

子動力 学計算 の部分 もC十十で コー デ ィング しておけば

,あ

らゆ る ウイ ン ドウシステムに移植が容易 になる。現在 の ところ,C十十は パ ソコンか らハ イエ ン ドの ワー クステー シ ョンまで動作 しているが, ダウンサ イジングとともに事 実上あ らゆる科学技術計算 を行 う上でc十十は利用可能である と考 え られる。 本論文 に於 いて以下 に掲載す るプログラムについては

,ANS1 2.1以

降の基準 のC十十であれば, 動作 す る。 しか し

,将

来 的 な発 展 を考 えれ ば

,テ

ンプ レー トや例 外 処 理 の利 用 が望 ま し く, ANS13 0の 基準 に達 した コンパ イ ラの利用が望 ま しい。具体 的に開発 に用 いた開発環境 は

Bonand

C十

+4 0 for DOS,Windows and WindowsNT(英

語 版

)で

あ る。特 に

,ユ

ー ザ ー イ ン ター フェ イス の 部 分 は Borland c十

+40に

付 け ら れ て い る ク ラ ス ラ イ ブ ラ リー Object Windows Library 2 0

(OWL20)を

用いて書かれている。 しか し

,分

子動力学計算本体の プログラム コー ドの動作 は

MicrOsOft VisuЛ C+十 Development System for Windows 1 5(英 _語版

_)及

_{び Symantec C十 十}PrOfes sbnЛ 6 1 for Windows(英語版

)の

上 で も確認 してい る。 また

,ユ

ーザー イ ンター フェイスの部

分 をVisual C十十やSymanted C十十に添付 されている MidrOsoft Foundaion Classを 用 いて コーデ イ

ング して もほぼ同様 に実現す るこ とがで きる。

以上, プログラムの開発 に直接用 いた機種 は

DELL OphPlex 466/LVで

あるが

,ANS13 0の

基準 を満 た している

IBM AⅨ

CI十

/600010及

び

IBM C set++/600021上

で も計算 本体 の部分 は動作す るので

,

これ に

AⅨ

_{上 のユ ーザ ー インターフェイス開発 環境 を組 み合 わせ れ ば}

_,IBMの

ワークス テーシ ョン上への移植が可能で ある。 さらにほか の機種へ も ANS1 3 0の 基準 を満 たす

C‖

コ ンパ

石 井 晃

3.分

子動力学の理論 3.l Pseudopotential,去 近年注 目されてい る

,第

一原理 か らの分子動力学計算 にCar Parrinello法[2]と呼 ばれる方法が ある。 この方法 は計算規模 は大 きいがか な り精密 に計算 す る こ とが で きる。 このCar―Parrinello 法 において電子状態 はpSeudopotendЛ _{法 で計算 す る。一般 に結 晶の電子状 態 を計算 す るには結 晶} 内の電子 の波動 関数 を何等 かの基底 関数で展 開 して計算 す る。数学的 に扱 い易 い平面波 で展 開す る のが一番 いいのだけれ ども

,電

子 の感 じる原子 ポテ ンシャルは中心部 で急激 に大 き く(負に大 き く) なっているので

,そ

こを正確 に しようとす る と非常 にた くさんの平面波 を加 える必要が出て きて実 用 的で はな くなる。pseudopotendal法は原子 の 中心部分 のポテ ンシャル を もっ と滑 らかで

,か

つ 波動関数 はなるべ く本物 に近 い もの を再現す る よう偽者 のポテ ンシャルを作 り電子状態 を計算す る 方法である。古 くは固有値 さえ再現 されていればいい として箱型 な どかな り適 当な擬 ポテ ンシャル を用 いていたが

,最

近では第一原理的 な きちん とした擬 ポテ ンシャル を作 つて計算す る ようになっ て きた。 Hartree Fock近似 による波動 関数 は次 の方程式で計算 され る,

身

″φ

ν

三

_[―

_弓

_井

▽

2 φν=Cッφν こ こ で,

フ

盟

Xア

)=暑

_[2免(ア)一

え

(ア)] はコア電子 との クロー ン相互作用 お よび交換 ・相関相互作用である。

免

①

=μ

②考先②

'テ 2

え①Ψ

①

=撫

)∫

家の

■ ②房

2

Philllps―

Kldnman[3]は

このHartree Fock方程式 を次 の ように書 き改め る。

呵

十 Ｚ 一 γ 一

[孝

▽

2¥十

フ

特″

刊

χ

=Ъ

χ

オブジェク ト指向設計による分子動力学計算プログラム

153

こ こ で,

ド

K=Σ

(ε

ν― ε

♪

1先)(究￨

は擬ポテンシャルであり, χは価電子的で

,且

つコア電子の軌道と直交する事を要求されないある

関数である。従ってχは次のうように表す事が出来るであろう。

χ

=φ

ν

ttΣ

″θ

ここで α♂はパ ラメー タであ り

,こ

れ を適 当に選 んで関数

Xが

ノー ドを持 たない ようにす る。 しか しなが ら

,こ

の方法では擬波動 関数 χが一意 にはな らない。 よ り良い擬 ポテ ンシャル を選ぶため には

,擬

ポテ ンシャル 自身がそれ を求め る ときの価電子 の波 動関数 に依存 して しまってはな らない。 なぜ な ら

,

もし依存 している と

,対

象 とす るエ ネルギ ーを 変 え る た び にい ちい ち擬 ポ テ ンシャル を計 算 しなお さな けれ ば な らな くな る。 そ こで改 め て Hartree―

Fock方

_{程式 を見 てみ ると}

_,交

_{換 ・相 関相互作用 の項が価電子 の波動関数 に依存 して しまっ} ている。特 にコア電子 の軌道 のある領域 での振 る舞 いに依存 して しまっている訳である。

Phillips Kldnmanの 計算 したLiの_{波動 関数 の図 を見 る と}

_,原

_{点付近での挙動}

_,即

_{ち コア電子軌} 道の領域 での振 る舞 いは

S,P,D,の

別があ り

,か

つ

_{, 2Sと}

3S, 2Pと

3Pの

振 るまいが よ く似 ている事 に気付 く。つ ま り

,角

運動量別 に交換相互作用 の部分の振 る舞 いを考 えれば

,い

ろい ろな エネルギー領域で同 じ擬 ポテ ンシャルが成 り立つであろ うとい うとい う予想が出来 る。 そ こで

,有

効 なポテ ンシャル として次 の形 に書 く事 に しよう。

レ

4万=f跳71(つIJ〉

〈′

￨ ここで

_,

かは実際には￨,ので

,必

要なだけ球面調和関数の和 を取る。 この￨′め 〈肋￨は射影演算子 になってお り

,こ

ちらで選ぶ波動関数の基底 (普通は平面波に採 るが

)か

ら

,そ

れぞれの角運動量 成分 を抜 き出す形 になっている。ある程度以上大 きな角運動量成分については同 じポテンシャルを 用いて近似する。 T/J=堵

+1

カ″ι>λ 従って, △‰

=71-7八

十二とか くと

,有

効ポテンシャルは次のように書ける事 になる。

154

石 井

晃 x万三7λ十二(つ十 `乳 △ 71(DIサ )(】￨ 実際の計算 では Sと

p,場

合 によつてはdまで擬 ポテ ンシャル を作 り

,そ

れ よ り大 きな角運動量成 分は

dあ

るいはfについてのポテ ンシャルで近似す る事 になる。

3.2 Norm―

Conserving Pseudopotentiol

Harmann,Schliter Chiang[4,5]に よつて提 示 され たNorm― Conserving Pseudopotentialで は

,擬

ポテ ンシャル を決定す る条件 として

,次

の関係式 に着 目 した。 21●

ψ

)2子

チ

lnψ

h=4π

lη 2″2,7 これは適当な球面 γ=父 より内側の波動関数のノルムさえ元の真の波動 関数 と一致 していれば

,対

数導関数の ″=貴 での一致 は自動 的に保証 されるという関係式である。そこで

,

この

Rよ

り外側で は真のポテンシャルに一致 して内側ではずつと滑 らかになっている擬 ポテンシャルを作 って擬波動 関数 も外側では本物 と完全 に一致するように作 り

,そ

の時にノルムが真の波動関数 と擬波動関数で Rよ り外側では一致するように する。そ うすればこの関係式か ら

,バ

ン ド構造の計算 に必要な ″

=

Fに

おける波動 関数の対数導関数 を一致 させ る事が出来る。Hamann,Schliter and Chiangに よろ

Moの

場合の真 の波動関数 と擬波動関数の図に見 るように内側で節が無 く滑 らかになってお り

,

し か もノルムは一致 している。 具体的には真のポテンシャルをまず用意 し

,そ

れ を ″

=Rよ

り内側で滑 らかになるように変 えた 擬ポテンシャルに直す。その時の関数 として

とする。ここで

7(の

は真のポテンシャルであり

, /(死

)=exp(π

35) がカッ トオフ関数 になる。

ここからノルムを合わせる手続きに入る。まず上の擬ポテンシャル 71'(Dか ら計算された擬波

動関数

Vrl(つ

に対し

, γ 一 η √ 十 の ７ ７ 一_η ／ 一 〓 つ

吟

オブジェク ト指向設計 による分子動力学計算 プログラム

155

吻⑦

=灼

_防ι

⑦

+始

ど

(羽

]

なる第二の擬波動関数を計算する。そしてこれが

/J/栃

ι

(の

α

″

=1

となるように規格化する。そうしておいて 晩 ′

(の

について動径方向の

Schrёdinger方

程式を逆に

解く事よって擬ポテンシャルカ猟か を求める訳である。

1

3.3 Separable pSeudopotential

1982年_{,Kleinman a♀}

d Bylander[6]は

Hamann,schluter and Chian3つ 擬ポテ ンシャルを改良 し て記憶容量が著 しく少 な くて済むようにした。Hamann,Schluter and Chね

ngの

ままの形であると, 擬ポテンシャルを平面波で挟んだとき 〈′T〒 J物_)7(つ (Jηlが生す_{〉∝ ∫ヵ}(力つ7(Dヵ(力の ん 7Pl(cosみた

) i

という形になるが

,こ

れでは波数 について二重 に和 を取 らなければな らな くなる。平面波の数 を

4

個

,M酌

R zone岬

の点 を

m点

_{取 るとす るとこれは 堅寧考}三堂 となるが

,酎

nlkltt mは

29ぐらいであるので

,

これは著 しいオーバ

_=ヘ

ッドを与 える。Kleinman―Bylanderは この積分 を二 つに分けてこの問題 を解決することを提案 した。つ まり

, ￨

∫九(力つ曳 (つ 72Jγ _ル (ヵの ■ (つん TPI(Cos先″) と分離す るわけである。 そのため

,擬

ポテンシャルを 7(つ=71°

σ

煎

(つ

十Σ△

71(つ￨,,初)(Й

η

l

と書いたときに

,こ

こであるエネルギー

Eι

に付いて擬波動関数を η

(の

とすると

,

△

『

ざ

と

(つ =Eιl島)(ξ

′

￨

156

石 井 晃 とする。このとき,

家の

=絆

とする。こうするとす ぐわかるように

△

71配(つ

切

(71=△

y(D9(つ

となる。しかも, この定義だと積分は分離され

,波

数の和はそれぞれについて単独に取ればいいこ

とになる。 しか し, この巧妙 な方法 に も落 とし穴がある。す ぐに気付 くのは

,

しか るべ きエ ネルギーのは波 動 関数 に依存す るため

,対

象 とす る電子 のエ ネルギー をいろいろ変 えた ときに きちん と真 のポテ ン シャル をエ ミュ レー シ ョン してい るか どうか

,つ

ま り,transferbilityが甚 だ怪 しくな る と言 うこ とで ある。 この点 はGonze,Kackel,Schefaer[7]に よって指摘 され

,対

数導 関数 を見 る とghost stateがとんで もない ところに出現す る。この対処法 は,その′ミに気 を付 けて擬 ポテ ンシャル を作 れ, とい う事である。 3.4 Uitrason PseudOpotential

い よい よ本題 であるUItrasoft Pseudopotendalに 話 を移そ う。 これは1990年 にD Vanderbilt[8] によって提案 された ものである。 これ まで見 て きたNorm―Conservihg Pseudopotendal(NCP)で は ノルムを保存 させ るため

,前

に掲 げた

Moの

例 で もわか るように決 して擬 ポテ ンシャル 自身それほ どソフ トにはなつていない。従 つて

, d電

子 や

f電

子 な どを扱 お うと して もまった くソフ トになっ てお らず

,従

ってfun pOtentialの計算 とさ して違 わない平面波 の数が必 要 となって事 実上遷移金 属 を交 えた計算が出来 なか った。 しか し

,そ

の制限は外 か ら見 た原子 の電荷密度 を本物 と同 じに見 せ るため にや む を得 ない制 限 で あった。 そ れ はself―consねtentに 擬 波 動 関数 を計算 す る と きや

Car Par ello法で動力学 をす る ときの原子 間の

Helmann Feynmann力

を正確 に計算す るため には やむ を得 ない制限 と思 われて きた。 ところが

,D Vanderbiltは

この制 限 は必ず しも不 可欠 な もの ではないことを示 したのであ る。

ではまずD,Vanderbiltの論文 に添 って

,NCPを

書 き直す ことか ら入 ろ う。運動エ ネルギー演算 子 を

Tと

書 き

,

まずall―electronの場合の波動 関数 (真の波動 関数

)は

オブジェク ト指向設計 による分子動力学計算 プログラム

157

と書 ける。一つのエネルギー固有値 の ところだけで参照す る場合

,先

のKleinman―Bylanderの や り方 に従 えば lχ,〉

≡

(勇

一

T―T/脇a)1陀)

と定義 して

,nonlocalポ

テンシャルを

晩≡

―

_場粥弊

と定義する。γ

='の

タ

ト

側では 駒=″どであることから

lχ

か≡

(ε ,一

T-7ヵσ

,)￨?ι

〉は明らかにRの 内

側だけに局在した関数である。すぐわかるように

yNL劣_)=￨‰) であ り

,従

って (党―T―‰翔―ン待L)lη

)=0

なるSchrёdinger方 程式 を満 たす。

もし

,複

数の固有値を参照 して擬ポテンシャルを導 くなら

(通

常はせいぜい 2つ であるが

),そ

の固有値をε

ぅ号として,ま ず

B方_=(7ιl χサ) とい う行列 を定義 し

,更

にそれ を使 って l

βl)=Σ

(Bl)方lχ

ブ

〉

という関数を定義する。これを使って

nonioc』 pseudopotenЫ

』 を次のように定義する。

ン

^√

と

≡Σβげ

l

β

τ

)(β

ブ

￨

158

石 井

晃

こうすれ

I出

すぐに計算してわかるように

ン

符

ι

l?,〉

=Σ

Bヴ_{lβ ,.)(β ,1望}

″

〉

=尋

島

lβ '〉

号

(Jl)'(χ

″

lη

″

)

=静

畢(Jl)力lγ″ ― )澪(す1)葬B協 =lχ″) となる。従 つて,

巧屹

￨,″

〉

=(ett― T―

覧づ

￨マ.)

であるから

, (鳥,一T― 駐 -7NL)IP″ 〉=0

となる。ここで

Bj=(71(身 ―

T―

晩♪

￨?ブ )

ず

'ア

劣

の

_卜孝移一

二

_霧上

―

晦⑦

_lη

⑦

である

.。

ただし

,雅

_(つ

=7η

_.(つ

っ

NOrm COnserving Pseudopotentialで

あれば

,こ

こで本物の波動関数と擬波動関数とで

QI」=●恥IWつP―(?,1狩 )P

なる量を

1定

義すれば

,Q方

=0で

ある。しかしュ

Vallderbildの

擬ポテンシャルでは 庁 々においての

波動関数と擬波動関数の対数導関数の一致のみを問題にし

,ノ

ルムは保存させない。従って

,当

然

Q・7は

ゼロにならない。まず

,重

なり

1演

算子として

d-1+Σ

_鋳

lβ,)(β

す

￨

オブジェク ト指向設計 による分子動力学計算プログラム

レ

^ば L=Σ

(B,十

考

@づlβ:)(β

ブ

と定義す る。す る と, これ に対 して, (T十陥じ十1ィ能―ε101η,〉

=0

という擬波動方程式が成立する。ここでエネルギー固有値 に相当するところに重なり演算子がかけ られていることに注意 されたい。 従 って

,ノ

ルムにこだわらずに滑 らかな擬ポテンシャルを勝手に作 って

,そ

れに対 して擬波動関 数を計算すればよいのだが

,ノ

ルムが保存 していないと

,電

荷の保存が怪 しくな り

,従

って分子動 力学で もっとも重要な原子 と原子の間に働 く力の計算が問題 になって くる。 こうした場合 にノルム は 〈ψザldlηデ〉々

=(町

IVl〉F と計算 される。従 って

,こ

の重 な り演算子

Sを

入れていれば

,外

か ら見た電荷 は保証 されるわけで ある。 以上 をとりまとめれば

,一

個の原子 に対する真のポテ ンシャルがあった として

,そ

れ とある適当 な半径 γ

=2で

対数導関数が真の波動関数に (関心のあるエネルギー領域で

)等

しい擬波動関数 を 与える適当な滑 らかな擬ポテンシャルを見つけておけば

,後

はそれについてノルムを計算すれば本 質的にはCar Parrinelloの 方法が可能になる[9,10]。

4.分

子動力学プログラム

4.1

クラスライブラ リー C十_{十で は クラス ライブ ラ リー とい う形 で一度作 った プ ログラムコー ドの再利用 をす る こ とが で} きる。先 に述べ た ようにオブジェク ト指向言語 ではデー タとその手続 きをひ とかたま りに して クラ ス とい う

,い

わば

Cに

おける構造体 やFortran90に終 け るモ ジュール を拡張 した もの をプ ログラ ムの基本的な単位 とす る。そのクラスに適当な初期条件 を与 えてメモ リ上 に展 開 した ものをそのク ラスのイ ンス タンス とい う。 クラス ライブラ リー とは

,そ

うした これ まで に作 ったプログラム コー ドを良 く設計 された複数の クラスの まとま りとして保存 し

,新

しいプログラムの製作 の時 には これ

石 井 晃 まで作 つたプログラムの うち必要 な もの を必要 な場所 で イ ンス タンス として展 開 し

,再

利用す る。 従 つて

,分

子動力学 の クラス ライブラリーは開発 とともに次 第 に多 くの クラス を含 む ように発展 さ せ ることがで き

,か

つ古 い クラス をイ ンス タンス として生成す ることで

,以

前 のバ ージ ョンの プロ グラムコー ドも変更 な しに実行で きる。 この特徴 は分子動力学 クラス ライブ ラリーの移植性 に も貢献す る。数値計算 をそれだけで ひ とま とまりの クラス ライブラ リー とし

,そ

れぞれの ウイ ン ドウシステムでは用 いた部分 の クラス をイ ン ス タンス として生成すれば

,そ

の ままそのウイ ン ドウシステム上で分子動力学計算が稼働す るこ と になる。 グラフイカル・ユーザー・イ ンター フェースは市販 の表計算 ソフ トな どとは違 ってそれほ ど複雑 な もの とはな らないか ら

,そ

れぞれの ウイン ドウシス テム ご とに多少 の プログラ ミングに よって

,分

子動力学 クラスライブラリーを用 いた分子動力学計算 ソフ トウェアがで きあが るこ とに なる。 分子動力学の場合

,た

とえば原子 を

Atomと

い うクラスに して しまうのは一つの方法である。 し か し

,そ

れではあ ま り広す ぎるので

,こ

こでは一個 の原子 をい くつかの クラス に どう分 けるか を考 える。 これが この論文 の主題 であ り

,そ

して これか ら発展 してい く分子動力学計算 ライブラ リーの 全体 を決め る骨格部分 となる。 ここで は一つの原子 につ いて

,そ

のnorm― conserving pseudotential

やultrasoft pseudopotend』 を計算 す ることを念頭 に置いて

,

まず真 のポテ ンシャルに対 して動径 方向の波動関数 を計算す る とい う計算 としては簡単 な問題 に絞 って考察 してい くことにす る。

4.2

グ リッ ドクラス Schr6dinger方 程式 で

,動

径方 向の波動関数 はポテ ンシャルが決 まっている とすれば解 くのはそ れほど難 しくない。求める量 は波動 関数であ り

,必

要 なのは動径方向の座標 の値 とポテ ンシャルの 値である。 ここで

,一

つの クラスの中にこの計算 のすべ てを押 し込 むことも可能である。た とえば 次の ようにクラス を作 ることがで きる。 class Wavefunction ￨ public i complex wavefunction[N]; private i

double grid[N];

double potential[N]; ptlblic: SolveEq();

オブジェク ト指向設計による分子動力学計算プログラム 16ユ ￨ このクラスでは波動関数の値 と,それを求める命令である

SolveEq()が

外 か ら参照で きるパブリッ クであ り

,グ

_{リッドとポテ ンシャルはプライベー トメンバになっている。たとえば学生に対する演} 示 として波動関数 を計算 してみせるだけならこれで十分であろう。 しか し

,一

度真の波動関数につ いて波動関数を計算 した後

,そ

れ とある適当な半径

r=Rの

ところで対数導関数が等 しい擬波動関 数を求めてい く作業 を考えると

,動

径方向の波動関数 を求める部分 はなるべ く共通のコー ドで

,つ

まり

,コ

ー ドの再利用ができることが望 ましい。その意味でいえば

,上

の クラスは落第で

,違

うポ テンシャルについてはこのクラスのコー ド全体 をエディタでコピー した上で異なるポテンシャルに 合わせて一部 コー ドを書 き換 える必要がある。それ も

,SЛ

veEq()本

体 の一部 を書 き換 えること になる。 これでは再利用 とはいえない。 従 って

,グ

リットもポテンシャル もこうしたクラスのメンバーとしてではな く

,独

立 したクラス としてこれをインスタンス として生成することにする。さらに

,グ

リッ ドを定数で刻む場合 と指数 関数で刻 む場合

,さ

らには Harmann Skillmannグ リットのようなものな ど

,い

ろいろなグリッ ド が考えられるので, グリッ トとポテンシャルも切 り放すことにする。そうすると

,グ

リッ トの違い を引数 としてグリッ ドのインス タンスを取 るときのsignitureで_{区別 して扱 う polymorphttmが 可} 能になる。 まず

,等

間隔で刻むグリッ ドのクラスは次のようなものである。

const unsigned int Ngridこ 三500;

const unsigned int N2grid=Ngrid*2;

class SCOnstGrid l private:

int T; pubhc:

int Z,L,N; //charge,angular momentum, //quantum number

double maxR,minR,dR,hR;//max,min,and difference of //the grid _ duble gridR [Ngrid]; //grid

double grid2[N2grid]; //half grid fOr Runge― Kutta

pubhc:

1鬱

1

石 井

晃

SConstGrid()‖

; //dos,uctor

void makeGFid(): //procedure to lllake grid

void maxGFid(); //proceduFe tO get maxinun

l;

4.4

ポテレー

シャルクラスー

￨ポ

テンシャルのクラスはポテンシャルを作るに当たってのグリッドをグリツドクラスから取 り

,

その刻みに従ったポテンシャルをデータとして蓄えるクラスである

.。

その際に

,フ

ルポテンシャル

以外に

nOFn「 COnservitt pseudopotentiЛ

〕

uttrasoft,seudopotialな

ども同じように解きたいため

,

まず

,抽

象基底クラスとして

SPotelatiolを 1定

め

,そ

れからの派生クラスとしてフルポテンシャルを

扱う

SF■1lPOtを

作る。こうすると

,別

に

SNOrmoonservc,SUltrasoftな

どが派生クラスとしてつく

れることになる。こうした場合のメウットは

,実

際にこれらのクラスのインスタンスを生成したと

きに

,記

憶場所として基底クラスの

SPotenialの

メン′ヾ

が先に並帆 その後にそれぞれの派生タラ

ス特有のメンバが並ぶ。従って

,別

のクラスで

S'otenЫ al*▼

とインスタンスを生成したときには

Vの 中身が実際には

I SFu1lPotで

_{あろう力潮}

Jの

であろう力｀

同じように扱えることになる。

実際の定義を示そう。

dtts SPotential i 〃 ポ テ ンシャ ル基 底 クラス 0よbliC: int type i doible V EN2grid]; public:

舒

otend』

01;

rtual∼SPotenia1 0 1 delete□

v;￨;

cottt SPotelltial&operator=(const SPottnial&pot),

//insert oporation

li

tlass SFillPOt i pubHc SPOtenhal l〃 フル ポ テ ンシ ャル の ク ラス

publiC i

SConstGFid*r; //グ リンドタラスのインスタンス

pllЫ ic i

オブジェク ト指向設計 による分子動力学計算プログラム

163

SFuHPot(SConstGrid*);//constructor

∼

SFullPot O ldelete r i delete□ V‖ ; //destructor

void makePotential(SConstGrid*r); _{〃 ポテ ンシ ャル を作 る}

￨;

class SNormconserv i public SPotential l //norm_conser ng pseudopotential

public i SConstGrid*ri 〃 グ リ ッ ドク ラ ス の イ ンス タ ンス pubhc i

SNormconserv(SConstGrid*r) ‖

, ∼

SNormconserv() 1 ;

void makePotential(SConstGrid*r); _{〃 ポ テ ン シ ャ ル を作 る}

￨;

dass SUltrasoft:pubhc SPotential l //ultrasO■ pSeudopotential

public i SConstGrid*r;

〃 グ リ ッ ドク ラ ス の イ ンス タ ンス

pubhc:

SUltrasoft(SConstGrid*r) ‖ ; ∼ SUitrasoft()‖ ; void makePotenhal(SConstGrid*r); 〃ポテンシャルを作る

￨;

まず

,基

底 クラスのSPotentialでポテンシャルの値 をメンバ として持つ。こうしておけば

,実

装 で どのポテ ンシャルについての派生 クラスでインスタンスを構成 して も

,利

用するコー ドの方で SPotenti』 を用いると宣言 しておけば

,そ

のコー ドはどのポテンシャルについて も共通に使 える。 また

,実

際 にポテ ンシャルを生成す る関数であ るmakePoteatialは

,す

べ て引数 としてSCOn_ stGrid*を とっている。こうしておけば

,グ

リッ ドとして等間隔ではな く

,た

とえば↓旨数関数 グリッ ドな どを用いるときで も

,別

に

makePotenhal(SExpGrid*r)を

クラスメンバー関数 として追加す れば,signitureか らC十十コンパイラは二つのmakePotentialを_{区別 して くれる。}

4.5

波動関数クラス グリッ ドクラス とポテンシャルクラスのインスタンスを含む形で波動関数クラスを構成する。波 動関数 自体はRunge―Kutta方 で計算するものとする。 class S VaveFunction ￨

pubhc:

164

石 井

晃 double E; 〃エネルギー protected: dOuble P,Q,dR,R; pubhc i

SConstGrid*r,

SFullPot*v; public:

double ,ave[Ngnd];〃

波動関数

double dwave[Ngnd];

〃波動関数の一階導関数

prOtectcd i

_{〃内部処理用}

doubに

wave fOr[NgFid]; 〃波動関数

double dwavettr[Ngr趙

];

〃波動関数の一階導関数

public:

SWaveFunction() ‖

; SWaVeFunction(int,int,int int); SWaveFunction(SConstGrid*,SFtlllPot*);

∼

SWaveFunction() ‖

i

virtttal void makeWave O,

protected:

_{〃波動関数の計算で用いる関数}

doubl definetttargetttE();

vOid getlnitValue_for(), void solveWave();

void runge_kutta_for(1五 ti); double dwavel(int i); double dwave2(int i)i

l ;

実装に当たっては

pubhcメ

ンバと

protectedメ

ンバの使い分けに注意されたい。一般に外部から

利用 されるもののみを

piblicメ

ンバ とし

}そ

れ以外は

prOtectedか privateと

する。ここで波動関

数の計算をさせる

makeWave Oは

lpuЫ

たであり

,外

から波動関数の計算を促すことが出来るが

,

その中身は

lsolveWave()そ

のものであり

,つ

まり実際の計算を させるプログラムコー ドの実体は

オブジェク ト指向設計 による分子動力学計算プログラム

4.6

ユーザ ーインターフ ェイス

ユーザーイ ンター フェイスはその コー ドがC十十で記述 されてい る限 り

,

どの ウイ ン ドウシス テ

ムにつ いて もほぼ似 た ように実装 で きる。上で挙 げた例 で言 えば

,SWaveFuncionの

イ ンス タ ン スを生成 して

,そ

の中のw ave[Ngrid]につ いて グラフを描 くな り

,適

当な処理 を して行 けばいい のであ る。幸 い

,MS WindOws,OS/2,Machintosh OS,X― Window,MS―

WindowsNTな

ど

,現

在 身 の回 りで メジャー に使 われてい るウイ ン ドウシステムの ほ とん どはC十十に よる コーデ ィングが可 能であ る。 ここではその一例 と して

,Ms windows 31及

び

WindOwsNT(Win32s)で

動作 す る,

波動 関数 の グ ラフを描 くユーザー イ ンター フェースの例 を Object Windows Library 2 0に よって コーデ イング した もの を付録 に示す。ほぼ同 じ構造で,Microsoft Foundadon Class 2 0(2.5)で コー デ ィングすることもで きる。

4.7

今後の拡張

ここで問題 としたのは一個 の原子 についての波動関数の計算 プログラム と

,そ

れ を擬 ポテ ンシャ ルに再利用 してい く場合 の クラス構成である。実際 に分子動力学計算 を行 うにはこれ らをまとめて

SAtomと

い う一 つ の クラス と し

,そ

こにその原 子 の位 置 と擬 波動 関数

,及

びそ の原子 に働 く力

(Hcllmaan Feynmann force)のデー タを集 めてお くのが便利 であ る。その うえで

,実

際 に単位格 子 な り

,ス

ラブ計算 で電子状 態 を計算 してい く部分 では

SAtomク

ラスのイ ンス タンスを扱 う原子 の数だけ生成す ることになる。その際 にすべ ての原子 にわたって共通 の操作 を繰 り返す必要性か ら, クラスの イ ンス タンス を扱 う一般 的なコンテナクラスがC十十コンパ イラーにつけ られていること, それに付 随 して

,不

可欠 とい うほ どで はないに して もテ ンプ レー ト機能の装備 が望 ま しい。その意 味で

,現

在 メジャーな コ ンパ イラの一つ,Microsoft Visu』 C十+1,0/1.51ま 望 ま しい開発 環境 で は ないが

,今

年 の うちに発 売が予定 されてい る次のver2,0におい て対応 す る ことが アナ ウ ンス され ているので これか ら開発 を進 めてい く分 には大 きな問題 ではないであ ろう。 メモ リの問題 について も次第 に解決 されてい く方 向であ る。一般 に小 さなプロ トタイプをパ ソコンで動作確認 した上 で ワークステーシ ョン上 で本格的 に稼働 させ ることが望 ま しいので

,そ

の意味で も分子動力学計算 の 本体 をユ ーザーイ ンター フェースか ら完全 に切 り放 してお くこ とが肝 要である。

5.結

論 本研究 においては

,オ

ブジェク ト指 向 に基づいた完全 なクロス プラ ッ トフォームによる分子計算 のプログラム開発 の重要性 を指摘 し

,そ

れ を実現す る分子動力学計算 プログラムの基本設計 を提示 した。

C++を

用 いてコーデ イング し

,原

子 の波動 関数 の計算 と擬波動 関数 の計算 を

pdymOrphism

166

石 井

晃

を利用 してクラス化する方法を示 した。

謝

辞

本研

―

究で示 したプログラムコー ドのうち

,原

子の動径方向の波動関数 を計算する

.Runge―Kutta法

の部分は逢坂豪先生が以前

BAMCを

用いてコーデイングしたものを C料 に移植 したものである。

参考文献

(1)M.A Ellis and BI ttrOustrup.Tlle Anmtated C++RettFelaCe Mantlal(Addttn Wetty,19921(邦 訳 :注解 C十 十

リファレンスマニュアル, トッパン,1992)

(2)R,Car and M Parrine19,Phぃ Rev Le世 55(19弱)2471

(3)」 c.Phillips and L Kleinman,Phys Rev.116(1959)287

(4)D,R Haman■ ,M Schidter and C.Chiang,Phys.Rev.Lett,43(1970)1494 (5) G.B.Bachelet,DR.Himann and M.Slhlater,P,ys Rev,B26(1982)4199 (6)L klこinmall and D.MI B,lander,Phys Rev Lett 43(1982)1425 (7)X,Co■ ze,Peter Kacken and M stheffler,Phys,Rev.B41(1990)12264 (8.)D Vanderbilt,Phytt RevI B41(1990)7892

(9)K.La鶴o,en,良.Car,C.Lee l■d D Vanderbilt,Phys Revi B43(1991)6796

1101 K Lassoneュ I A Pasquare■9,Ri Car,C.Lee and D Vanderbi比 ,Phys.Rev.B47(1993)10142

付録 分子動力学計算のプログラム実装部分

本文で紹介 しているグリッ ドクラス

,ポ

テンシャルクラス

,波

動関数 クラスの各メ ンバ関数の実 装部分 を示す。 また

,ユ

ーザーインターフェースのtf711を Borland Obj℃ct Windows LibFary 2,0ク

ラスライブラリを用いて示す。

//Sakyu ObieCtiVeヽ lolec,lar Dyta― alnics Class Library

// versio■ 1.00 May 9, 1994

// Calculate radial gridす Ootential of half― gride,

// and wave..`

#indude<iOstFeam.h>

#include<math.h>

#include“

atomlib.h"

オブジェタ ト指向設計による分子動力学計算プログラム

167

SConstGrid ￨

Z=z

L=1

N=n

T=t

i SConstGFid(int 2j int l,int n,int t=4)

for(int i=o,1(Ngrid l i十 +)

￨ gridR[司=0.0, ￨ for(i=0;i(N2grid i ll十) ￨ gFid2[￨]=αO, ￨

minR=0,00001,

￨

void lSConstGrid i i mattGrid()

￨ double anl,a9,dd,ai

a4=2,0*T*2.30259259;/れ

og(10)=2,30259259

aO=am i

dd=2.*N;

a=amttdd*log(aO)i

ヤ

hile((a―aO)〉 9,000oo01*aO)

￨ aO三=a

aO=amttdd*log(aO)i

￨

maxR=0.5*a*N/Z;

￨

163

石 井

晃

dR=(maxR―

ninR)/Ngrid,

hR=三dR*0.5; //incFement clf the grid

gridR[0]=grid2[0]=lalinR;

br(ht-1=1;i(Ngrid;i‖

) ￨

gridR[i]=gndR[i-1]十

dR;

￨ fOF(1=1;i(N2grid i二

十十

) ￨

grid2[i]=grid2[i-1]+hR,

￨ ￨

const SPotential&I SPotential i :― operator=const sPotential&p90

1

if(&,ot!=thi0

￨

for(int i=o;i(N2grid i i十

十

)

I

V[i] =pot,V[i];

￨ type=三pOt.type i ￨ return *this i l

SFullPto: :SFullPto(SCOnstGFid*CD:r(cg)

￨

type=FVLLPOT,

r―

)maxGrid();

r―

〉makcGrid O;

￨

オブジ■タ ト指向設計による分子動力学計算プログラム

169

for(int i=o,1(N2grid,i+ll i

V[i]=―

(ィ

ー

)2)/r―

〉

grid2[i]; ￨ ￨

SWaveFunc■9n i i SWavoFtlュcuon(int z,int L,int N,lnt T=4)

￨

SCorlstGrid*r=new SConstCrid(Z,L,N,T)￨

SFullPot*

▼=new SFullPot(1)￨

// v―

〉

makePotential(r),

// E=deine/target/E(r);

// makcwave(4,V)i

￨

SWaveFunctiOII:i SWaveFmction(SConstttid*rr,sF■ 11'Ot*w)ir(Fr),v lhrv)

ν

/You shOuid call this constructor aFter you construct r■ alad v■

￨

vOid SWaVeFunctio4::nakeWave()

￨

E =define_target_E()￨

getlnitValue_fOF O I .solveWave (), for(int i=0;i(Ngrid:i++) ￨ ヤ

ave[i] 100*wave_fbF[i]i

dw―

ave[I] 100ヤ

dwavefor[i]: ￨ // nOrlllalize O; ￨

double SWaveFunctiOn i:derineitaFgettE 0

￨

170

石 井

晃

double Energy=-9,5*(r―

〉

Z)*(r―

_〉

Z)/(r―)N)/(■

―

〉

N)*Er:

// E=-0.5*pow(・

―

)Z,2)/pOw(r―

〉

N,2.)*EFi

remrn EneFgy,

￨

void SWaveruncti04: :getlnitVJueifOr()

￨

〃Iniunl cOndi on of Runge_Kutta,loFWard

douわle Rl:

Rl=sqrt(2.*功

■

(r・

〉

minR),

P =exp←

Rl)*OOw((21*Rl),r→

二

十

1.);

Q =―

P*Sqrt(-2,燕

El}

dR =(lr― ).maxR)―

(■―minR))/Ngrid, return‐ i ￨

void SWaveFunctiOn i:sdveWave 0

￨

R =r―

〉

.gridR[9];

dR=r)dR;

wave/for[0]=P;

dwavefor[0]=Q;

fOF(lllt i=0,1(Ngrid,1+ll ￨

wave_for[i]=P;

dwavefor[i]=Q,

ringe_kutt広

_{_for(1*2),}

R =Rttdk,

￨ : ￨

void SWaveFunction I:ruIIge/kutta/for(int i)

i

double PP =P;

do!Ы

c QQ =Q;

オブジーェクト指向設計による分子動力学計算プログラム

171

double qO=dR*d市

な

V021i),

P =PP+0.5*p01

Q =QQ+0.5*電

9;

doible pl=dR*dwavel(i+1)

doubに

ql=dR ttdwave2(1+1)

P =PP+0.5*pl,

Q =QQ+0お

*■1; dou心

に

p2‐

dR*dwavel t+1)

doible qZ=dR*dwave2 1+1)

P =PP十

o21

Q =QQ+範

2, d6ibに

p3=dR*dttavel(1+2)

d。

■

b―

le q3=dR*dwave2(1+21

P =PP十

(pO+2-*(p二

十

p2)十o31./6.:

Q =QQttlpO+2,■

(01+。2)十

p3)/6i

￨

dollbに

熱

￨・

幾

veFttnctioll i l dwavel(IIltつ

￨

dow―bに f=Qtt lrl〉

七十

11/1r‐

)8rid2岡

)IP,

Feturn fi

l

doiblo SWaveFtln on i:dOubre SwaveFunction i:dwave12(intつ

￨

dOuble f=2.*(ヤ

ー

〉

V[i]―

E)*P―

(r‐

〉

L+1)/(1→ grid2[1])*Qi

Feturln fi

￨

void SWaveFun,■o■

￨:nermarize o

￨

double s=O,9;

釣T(illt i=1,1(Ngrm●

2,1=i+2)

￨

172

石 井

晃 ￨

S=S*dR/8,;

double C=sqrt(S)i for(1=0,1(Ngrid i i十

+)

￨

wave[1]=wave[i]/Ci

￨ ￨

double SWaveFuncion i:GetMaxValue()

￨

doible max=wa中

[9]f

forOnt i=0;i(Ngrid/2;i十 十

) ￨ if(wave[i]〉 max)

lmax=wave[1]#

￨ retuFn max; ￨ OWL.2-.0に よ,るユーザとインターフェース部分のコーデイング ヘ ッダー部分 と実装本体を合わせてあるが

,現

実は別のファイルにしている。

/* atom4app.clpp OVERVIEW

CItts defini10n for atonttn4Ap。 (TApplicatiO■

)

半/

1 include(。wI¥owlpch.h〉

十

prag4a hdrstop

#include(Owl¥翻

itfle.h)

#include(。

WI¥opensaveh〉

オ

inCl■de``atom4app,rhイ //Deiniton or all rttources.

〃 キ

ITAppllcadOn=atomwh4app‖

Class at。■

wh4App:publ■

TApplicadon i

p.rlvate:

オブジェク ト指向設計 による分子動力学計算プログラム

173

TOpenSaveDialog: :TD― ata FileData, //Data to control opcn/さaVeモ嗚standard dialog puわHc i

atomwin4App()i

ViFtual∼

atomvin4App()￨

vold OpenFllc(const char*fileName=0),

// 11atomwin4AppVIRTUAL_BECIN‖

public i

vltutt void ln

MainWindOw();

// ltatomwin4AppVIRTUALiEND‖

// 11atomwin4AppRSP_TBL_BEGINII

pFOteCted i void CmHelpAbout O;

void DrawWave(),

// 11atomwin4AppRSP_TBL_ENDII

DECLARE_RESPONSE_TABLE(atonwin4App);

￨ , // 11atomwin4Appll # inolude(Owl¥owipch.h〉 オ pragma hdlstop #i―nclude(owl¥vbxctl,h〉

# include f`trnwn4abd,h" //Definition of about dialogi

# indude

“tmywindw.h''

〃 ￨latomwil14App lmp的mcntatio対 ￨

//Build a response table for all lnessages/co■ lmands‐ handled

//by the apohcration //

DEFINE_RESPONSE_TABLE-1(atomwin4App,TApolicatiOn).

ん/ 11atomwil14AppRSP_TBL_BEGINII

EV_COMMAND(CMttHELPABOUT,CmlHelpAbout),

EV COMMAND(CM DRAWWAVL DrawWave),

// 11atomwin4AppRSP_TBL_ENDII

END_RESPONSE_TABLEI

174

石 井

晃

//

class SDIDecFFane:0■ blic TDecoratedFrame l

public I

SDIDecFrame(Tれ

た

_{indOW*parent,const.chaF far*titleJ TWindOw*dientWnd,BOOL}

trackMenuSOlection=FALSE,TModile*module=0):

TDecoratedFranle(paFent,tide,91ientWid,trackMenuSlec,on,InOdule) ￨ ￨ ∼

SDIDecFrame()

￨ ￨

￨;

atomwin4App: l atomwin4App()i TApplicati釘 ュ(“S‐akyu Obiect e MOlecular Dynalllics") ￨

//CO血InOn flle■le naaS and fllters for open/Save As dialogsi Flename and diFeCtOrd aFe

//computed in the member ftlnctions CmFneOpen,and CmFileSaveAs

FllcData Fla―

gs=OFN_FILEMUSTEXIST 1 0FN_HIpEREADONLY I

OFN OVERWRITEPROMPT;

FileDataISetFilter(“

ALL Files(*.*)￨ *,* ￨'り

,

//1NSERT〉)YouF COnStructOr cOdc heFe.

￨

atomwin4App l i atOttwin4App O

￨

//1NSERT〉〉Your destFuCtOr code here,

￨ //atOmwin4App //EEE== //Apphcatioll intlaliZation, //

vdd atomwh4App::I

th4anWhdO寸

() ￨

Clie■t=new TEditFlle(0,0,0);

SDIDesFrame*frane=ne4・

SDIDecFralne(Oi GetNalne(),new TMyWindow): //TMyM/indow in this new opeFattOn is a key paraneter to hnk

オブジェク ト指向設計 による分子動力学計算プログラム

175

//this Mainframe with the TMyttrindow cl ass of tmywindwicpp

nCmdShow=nCmdShow!=SW SHOWMINIMIZED?SW SHOWNORMALinCmdShowi

//Assign ICON w/this apphcation,

frame―

)Setlcon(this,IDI SDIAPPLICATTON);

//

// lenu associated with window and acceletor table associated with table, //

//AssOciate、vith the accelerator table.

frame―

_{)Attr,AccelTable=SDI MENUi}

MainWindow=fFame;

￨

//====二==ご==

//Menu Help AbOu atomwin4 exe conamand void atomwin4App: i CnHelpAbout()

￨

//

//ShOw the modal dialog= //

atomwin4AboutDlg(Mainwindow)Execute();

￨

int Owlmain(int,char*

_岡

)

￨

TBIVbxLibrary vbxSupporti //This apphcatiOn haSヽ /BX controls,

atomwin4App App i

int result i

result=App RunO,

retura result i

￨

void atomwin4App i l Drawれた_ave()

￨

//1NSERT)〉 YOur code heFe.

176

石 井

晃

/■TMアヽVi■

dowicpp OVERVIEW

Cl■

ss deAnitiOll fOF TMyWindOw(TWindOw), *ノ

#include(Owl¥。

wlpchih)

# praglEla hdstop

#includc(owl¥window,れ 〉

#include(Owl¥dch〉

//deす

ice context

#include(dasslib¥arraysilly

#itlclude rtlxp■ 1_dih"

#inclJdc ttatom4app.rh″

//Dゞ

in■ion 9f all resouFeeS・

イ

dude“

ム

tOm■b,■

'

〃分子動力学計算プログラムを引用指定

typedef TAFray(TPoin→ TPbints,

tyOedef TArraylteratoF〈TP10in→ _{TPointsheratOr,}

//11TWindow=TMyWindowll

cttss TMyWindow i public TWindow i

public.:

TDC*宮

rapllDC, TP01■偽* Lirle i

TPoints*Wavei

SColls崎

d*ri

〃分子動力1学_{計算 のグ リツ ドクラスの インス タンス}

SF,1lPot* v: //ポ

テ ンシヤル クラスの イ ンス タンス SWavё

Function* wi //波

動 関数 タラスのイ ンス タンス OrOttCted i i,tZ, illt L i int N; public:

TMyWindow(TWindO帝

*parent=0,ooist char far*titleO,TModuに 工

mOdile=o):

rm』

=TMyWindOw O i

Void Paint(TDC&,B00L.TRect&);

void drawXaxis(TDC■ )i

オブジェク ト指向設計 による分子動力学計算 プログラム

177

//

‖TMylVindoヤ

RSP_TBLiBECIN‖

protected i

void DraWave();

void EvLBItto■

Down(UINT modKeys,TPoint&pointl,

void inOutZ();

void EvRBu,oADo,n(UINT mё dKeys,TPOint&point)i

VOid lnpi■ ()￨ void lnputN()￨ voidi Setlnputs();

//

‖ThrlyЛ

WindowRSP/TBL/ENつ

_‖

DECLARE_RESPONSEI_TABLE(TMyWindOw),

￨!

〃 ￨ITMyWindowll

#i■

cludc(Owl¥owlpchih) #inctude(owi¥inputdia.lly #pragala hdFS‐tOp

■

include“ atOmllbih" //

//B■ 1ld a FespOnse tablc foF all lne.ssages/comlllandS handledl //by the application.

//

DEFINEttRttPONE_TABLEl(TMy訥

「indO帯

_,TWindOw)

//11TMyWindO.wRSPttTBL_BEGINII

EVttCoMMAND(C.M_DRAwwAVE,Dra

Wave),

EVttWM_LBUTTONDOWN,

EV_COMMAND(CM_INPUT_を

,Input2).

EV_WM_LBIUTTONDOWN,

Ev_COMMAND(CM_INPUT_￨,InputL),

EV_COMMAND ICM_INPUT_N,InputN),

EVttCOMMAND(CM_SETINPUT,SetlnputS),

// 11TMttindOWRSP TBL ENDII

END_IRIESPONSE_TABLE;

〃 ￨ITMyWindow lmplemelltatねnl l

178

石 井

晃

TMyWindow i r TMyWindow(TWindow*―

paFent,Const ellar far*tl■

c,TMbdule*modulc):

TWindow(parent,itiei modull)

￨

//Override the defa■lt window style fOr TWindOwi

Attr.Style l=WS_BORDER I WS_CAPT10N II〃

S_CHILD I WS_A/1▲

XIMIZEBOX I

WS_MINMIZEBOX I WS_VIttBLE I

〃INSERT〉〉YouF COnSttu.Ct6r codeにre

z=251

L=1,

N=21

1n五

●

arent,0,① ; graphDC.=01 Line=new TPOints(10i9,19);

Wave=new TPoin榛

(10や,lo)￨

￨

TMttWindOw::TMァ

windOw()

￨

DcstFOy();

//1NSERT〉)Your destruc‐tor cOde here.

delete gFaphDC i

delete Linё , delete Wave;

￨

void TMyWindow i r DrawWave 0

￨

//11NSERT))YOur code heFC. if(l graphDC)

SttCapture(),

graphDG inev TClientDC(*血

is)i gFaphDC■

_{)SotMapMbde(MM_Ⅲ}

METRICl i

オブジエークト指向設計による1分子動力学計算プログラムー

179

gFaOhDC‐ 〉籠tWind。,Org(pO)・

dravXaxis(graphDC);

dFaWYa `(grapllDC),

SConstCrid*r=hew SC6■ sOrid tt L,N,4)￨

r「)mattG

d O;

r」

)makecrid O i

SF■1lPot*v=new SFullPot.(F), v⇒ nakePotentid(う : br(int―

I=0;i(N2grid,i‖

) 1

iF(v―

_)V[i]〉

‐

6000&&v―

〉

v[i](-6900)

￨

int px=(100oO*(r―

〉

gユid2[i])),

int,y=1.*(v⇒ V[』

)￨ TPOint p=TPOlnt(px,ay);

switcll(i)

￨

.case O :

graphDC―

_)Moヤ

eT6 1px,py)i

bFeakl; default i

graphDO―_{)LineTo lpx,py)￨} bFeak i

￨

Line―

〉

Add(pl i

￨

SWttveFunCttn*w

〒

new SWaveFunction(F,Vl,

v―

〉

makさ

Wave O;

double max =wiCetMaxVarue();

foF(1=0,i(Ngrid i i+十)

180

石 井

晃

1■

七

―

px=(10000*●

―

)gridR[i]));

int py=6000*(w〉 wave[iン max)i

TPoint p=T'Ointlpx― ,py), switch(1)

￨

case O :

餅

aphDC―

_〉

LineTO(。

4py);

bFeak,

default i

graphDC・_{)LineTO lpx,py),}

bFeak;

￨

Wave―

Add(o)￨

￨

void TMyWindow.: :Paint(TDC&dc,BOOL,TRect&)

￨

B00L first =τ

RUE,

TPointslteratOr i(*Line), TPointtlteFatOr i(*Wave),

dc.SO出

rapMOde(MM_HIMETRICl i

TPottt pO =TPOint(-500,5000)￨〃

oFttin of tte coordi■ ate

dc.SetWindovorg(pO)i

drawXaxis(&dc);

drawYa

s(&dc)i

whne(il TPoint p=母+i if(1 4rNst) dc.LineTo lp); else

オブジエタ ト指向設計による分子動力学計算プログラム

131

dc,MoveTo(p);

first=FAL.SE; ￨ wlh■ e li) ￨ TPoint q=〕+十, if(l irstl dc.LineTo(q), else ￨

dc,MovOTo(q);

first=FALSE;

￨ ￨

void TMyWindow i l draWXaxis(TDC*D⑥

￨ TPOint,o咆

=TPOint(0,0),

TPoint pend=TPOint(6000,0)i DC―

〉

MOveT。 (porgl, DC―

_〉

LineTo(pelld); ￨

void TMyWindow i :drattYatts(TDC*DC)

￨

TPOintブ

bOttom =TPoint(91-6000),

TPoint ytop =TPoint(0,6000),

DC―

)MoveT。

(yわotto五),

DC―

_〉

LinoTo(ytoo),

￨

石 井

TWれ

dow::E▼

LButtOnDowa(“

odKeys,po 0 1 〃INSERT〉_{)Y.our ω}de―heFe,

r lgreaphDC) ￨ ReleaSeCapture 0 1 Lille→

FIush O;

Waver)FIush(),

hvalidate O, dclete griphDC i l ￨

void TMyWindω

tt I14putZ 0

￨

//1NSERT))YouF COde here, thar inputZ[6];

if((ThputDialog(Fhお_{,'Atomic NImberⅢ ,・}Input a rlew z numbeF"

ino■tZ,sizωf(inputz))).Exec

圭

_{atOi Onputzl i}

￨

void TMyWilldow i t EvRButtonDOwn(UINT modKeぃ _,TPoi

―

&poiltt)

￨

TWindOw i :EvRBitto■

Down(modkeブ

s,point, //1NSERT))YOur code here`

Invalidate O;

￨

void TMyWindow::IttutL.()

￨

〃

INttT)〉

YouF COde here,

O=IDOK)

オブジェク ト指向設計 による分子動力学計算プログラム

188

if((TInputDialog(this,“ A■

gdar Momentum",

_{“Illput a new L number!,}

inputi sizeOf OnputLl)).Ex℃

cute()=IDOK)

￨

L =atoi(inputLl,

￨ ￨

void TMyWindOw: :IュputN()

￨

//1NSERT〉

_〉

Your codc here.

Char input躊];

if((TI■,■tDialog(this,“Main Quantum NultheF",`Input a llew N numbeF",

inputN,sizeof(inputN))).Execu俺

_0=IDOK)

￨

N =atoi(inputN),

￨ vdd TMywindOw i i setlnputs 0 1

//1NSERT)〉 Ybur codc heFei

TINPUTALL_DIALOG adia10g=new TINPUTALL_DIALOG lthis,IDD_INPUTALL,0)￨

adi』og SetLvaluc(L),今 adialog.SetLvaluc(N); if(adialog,Execute O=IDOK) ￨ L=adialog,GetLvalue O; N〒adittogi CetNvalue()i l