====== 初めての第一原理計算（インタープリタで一歩一歩）(ase3/Jacapo) ======
pythonがインタープリタであることを活かして、一つ一つ説明しながら全エネルギーの計算を実行してみよう。

計算は、
  - [[Atomオブジェクト]]を用いて[[Atomsオブジェクト]]に原子配列の情報を設定する。
  - [[Jacapoオブジェクト]]に計算条件を設定する。
  - 実際の計算を実行する。
の順で行う。


===== 使用する変数・定数の設定・定義 =====
まず、python を起動する。

  $ python

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

  >>> a = 4.0
  >>> d = 1.1

このような場合

  >>> a = 4.0; d = 1.1

のようにセミコロンを使って一行で入力することもできる。

===== 原子配列の設定 =====
Atomic Simulation Environment (ASE)のライブラリかAtomsモジュールを読み込む。

  >>> from ASE import Atoms

Atomsと複数形になっていることに注意すること。これにより[[Atomsオブジェクト]]が定義できるようになる。<note>定義の仕方は、

  適当な変数名 = Atoms ([Atomオブジェクト, Atomオブジェクト, …], cell = …, pbc = …) 

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

  >>> s1 = ListOfAtoms([])

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

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

  >>> from ASE import Atom

C(炭素)原子とO(酸素)原子とに対応するAtomオブジェクトを定義する。<note>定義の仕方は、

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

箱（単位格子）の形と大きさをよく考えて原子の座標をA単位で設定する。箱が直方体でない場合は斜方座標になるので注意。

  >>> p1 = Atom('C', (a/2, a/2, a/2))
  >>> s1.append(p1)

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

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

とすれば、中心からaだけずれた位置に酸素原子を配置することができる。
このとき閉じ括弧の数を間違えやすいのでよく注意すること。



 

  >>> 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を実行します。


[[ase2dacapo_2|旧バージョン（ASE2/Dacapo編）]]