pythonがインタープリタであることを活かして、一つ一つ説明しながら全エネルギーの計算を実行してみよう。

計算は、

1. Atomオブジェクトを用いてListOfAtomsオブジェクトに原子配列の情報を設定する。
2. Dacapoオブジェクトに計算条件を設定する。
3. 実際の計算を実行する。

の順で行う。

まず、python を起動する。

$ python

単位格子 (COを入れて周期的に並べる立方体) の一辺の長さ (4Å) と、C-O間の距離 (1.1Å) を変数に代入しておく。 長さの単位はオングストローム (1Å=0.1[nm]) 。

>>> a = 4.0; d = 1.1

Atomic Simulation Environment (ASE)のライブラリからListOfAtomsモジュールを読み込む。

>>> from ASE import ListOfAtoms

これによりListOfAtomsオブジェクトが定義できるようになる。 定義の仕方は、

適当な変数名 = ListOfAtoms ([Atomオブジェクト, Atomオブジェクト, ...]) 

このListOfAtomsは、後に述べるAtomオブジェクトを[]で囲んだ「リスト」を引数にとる。 したがって最も簡単な定義は「空」の「リスト」をつかって、

>>> s1 = ListOfAtoms([])

とすることである。 このListOfAtomsオブジェクトに対していろいろな操作を行うことによって計算が進行する。 このListOfAtomsという箱に原子を配置し、その箱を周期的に並べたものを計算対象にしていると考えておけばよい。 上の例ならば、s1という原子も入っていなければ大きさも決まっていない（実はデフォルトで1Å×1Å×1Åの箱が用意されている）箱を用意したことになる。

Atomic Simulation Environment (ASE)のライブラリからAtomモジュールを読み込む。これでAtomオブジェクトが定義できるようになる。

>>> from ASE import Atom

C(炭素)原子とO(酸素)原子とに対応するAtomオブジェクトを定義する。

適当な変数名 = Atom ('元素記号', 位置ベクトル) 

後で箱（単位格子）の形と大きさを設定すると自動的に拡大縮小されるので、位置ベクトルは箱に対する相対座標だと考えておけばよい。

>>> p1 = Atom('C', (0.5, 0.5, 0.5))
>>> s1.append(p1)

とすることで、s1の中心に炭素原子を配置したことになる。Atomオブジェクトを一々変数に代入しなくても直接

>>> s1.append(Atom('O', (0.5 + d/a, 0.5, 0.5)))

とすれば、中心からd/aだけずれた位置に酸素原子を配置することができる。 このとき閉じ括弧の数を間違えやすいのでよく注意すること。 後から箱のx方向の大きさをaに設定すれば、酸素原子と炭素原子の原子間距離はdになる。

>>> from ASE import ListOfAtoms
>>> s1 = ListOfAtoms([p1, p2])
>>> s1.SetUnitCell([(a,0,0),(0,a,0),(0,0,a)])

>>> from Dacapo import Dacapo
>>> c1.SetNumberOfBands(10)
>>> s1.SetCalculator(c1) 

>>> c1.WriteAsNetCDFFile('input_data.nc')

まず、

>> a = 4.0; d = 1.1

三つの基本単位格子ベクトルを変数に代入しましょう。ベクトルは()で囲んだ「タプル」を使います。

>> a1 = (a, 0, 0) 
>> a2 = (0, a, 0) 
>> a3 = (0, 0, a)

炭素原子と酸素原子とを配置する座標（位置ベクトル）を変数に代入しましょう。座標は単位格子に対する相対座標で設定します。

炭素原子の座標を単位格子の中心に

>> v1 = (0, 0, 0)

酸素原子の座標を中心からx方向に 1.1[Å]/4[Å] だけずれた位置に設定します。

>> v2 = (d/a, 0, 0)

>> from ASE import Atom

です。

>> p1 = Atom('C', v1) 
>> p2 = Atom('O', v2)

C(炭素)原子ひとつとO(酸素)原子一つとからなるListOfAtomsを定義します。定義の仕方は、

>> s1 = ListOfAtoms([p1, p2])

基本単位格子ベクトルを設定します。

>> s1.SetUnitCell([a1, a2, a3])

Dacapoモジュールを読み込みます。これでDacapoオブジェクトが定義できるようになります。Dacapoオブジェクトの定義は、

適当な変数名 = Dacapo ()

のように何も指定しなくても定義できます。必要なパラメータは後から設定します。

>> from Dacapo import Dacapo 
>> c1 = Dacapo()

展開に用いる平面波基底の最大エネルギーを300eVに設定

>> c1.SetPlaneWaveCutoff(300)

計算する電子状態数、 (4(Cの価電子数)+6(Oの価電子数))/2(スピン重複度)+5(非占有状態の余裕)

>> c1.SetNumberOfBands(10)

「計算モデル」に「計算方法」の適用

>> s1.SetCalculator(c1)

「計算」の詳細をファイルinput_data.ncに保存

>> c1.WriteAsNetCDFFile('input_data.nc')

pythonを終了後

% dacapo.run input_data.nc

または

% dacapo.run input_data.nc -out output_data.text

のようにしてDacapoを実行します。