地域密着型エリア広告配信リクルートの地域密着型広告ドコイク?アドネットワーク Ads by ドコイク?[無料でホームページを作成] [通報・削除依頼]
[無料でホームページを作成] [通報・削除依頼]

GPAW

(準備中)ここではGPAWを用いた基本的な計算からXANESの計算までを解説する。
-------------------------------------------------------------------------------
下記のHPを参照すると良い。
http://www.csc.fi/english/csc/courses/archive/gpaw-2008-01
-------------------------------------------------------------------------------
■ インストール(Ubuntu 9.10以降)
  下記のコマンドを順番に入力することでインストールされる。
□ ASE [1,2] & GPAW [1,3]( GPAW 0.8.0 )
  sudo add-apt-repository ppa:campos-dev/campos
  sudo apt-get update
  sudo apt-get install gpaw
■ インストール(Ubuntu 12.04 64bit)(2012/9/16時点では原子位置の構造最適化はagでエラーが出る)
□ ASE [1,2] & GPAW [1,3]( GPAW 0.9.0 )
  sudo bash -c 'echo "deb http://widehat.opensuse.org/repositories/home:/dtufys/xUbuntu_12.04 /" > /etc/apt/sources.list.d/home_dtufys.sources.list'
  wget http://widehat.opensuse.org/repositories/home:/dtufys/xUbuntu_12.04/Release.key && sudo apt-key add Release.key && rm Release.key
  sudo apt-get update
  sudo apt-get install python-ase
  sudo apt-get install python-matplotlib
  sudo apt-get install gpaw
※ もし、動作しない場合は、下記に記載されているパッケージもインストールしてください。
[1] https://wiki.fysik.dtu.dk/gpaw/install/Linux/Ubuntu_ppa.html#ubuntupackage
[2] https://wiki.fysik.dtu.dk/ase/download.html
[3] https://wiki.fysik.dtu.dk/gpaw/install/installationguide.html#installationguide
-------------------------------------------------------------------------------
■ Ubuntuでのインストール[1](上記で動作しなかった場合)
・必要なパッケージ
  sudo apt-get install python-dev
  sudo apt-get install python-numpy
  sudo apt-get install liblapack-dev
・ 推薦するパッケージ
  sudo apt-get install ptyhon-scientific
  sudo apt-get install python-matplotlib
・ 並列化用パッケージ
  sudo apt-get install openmpi-bin
  sudo apt-get install libopenmpi-dev
・ 推薦する並列化用のライブラリ*
  sudo apt-get install libscalapack-mpi-dev
・ ドキュメント
  sudo apt-get install python-sphinx
  sudo apt-get install texlive-latex-base
  sudo apt-get install texlive-latex-extra
  sudo apt-get install dvipng
  sudo apt-get install povray
[1] https://wiki.fysik.dtu.dk/gpaw/install/Linux/Ubuntu.html
* 上記のHPにあるように設定が必要。
-------------------------------------------------------------------------------
■ インストール(Ubuntu以外のOS)
□ ASE[1]
・RHEL/CentOS:
・Fedora:
・openSUSE:
・Debian:
□ GPAW[2]
・RHEL/CentOS/Fedora:
   yum install gpaw
・openSUSE:
   yast -i gpaw
・Debian:
   sudo apt-get update
   sudo apt-get install gpaw
※ 各OSでのインストールは文献[1]に書かれている。
[1] https://wiki.fysik.dtu.dk/ase/download.html?highlight=ubuntu
[2] https://wiki.fysik.dtu.dk/gpaw/install/installationguide.html#installationguide
-------------------------------------------------------------------------------
■ コンパイルする場合(プログラム開発者向け)
pythonによって書かれているので、NumPy及びpython-devが必要になる。
 
  NumPy : http://www.scipy.org/Download
    tar zxvf numpy-1.6.0.b1.tar.gz
    cd numpy-1.6.0b1
    python setup.py build
    python setup.py install
  以上でNumPyをセットアップしたら下記を行います。

  ASEのインストール https://wiki.fysik.dtu.dk/ase/
    tar zxvf python-ase-3.5.0.2101.tar.gz
    cd python-ase-3.5.0.2101
    python setup.py install

  Atomic PAW Setupsの設定
  https://wiki.fysik.dtu.dk/gpaw/setups/setups.html#setups
    tar zxvf gpaw-setups-0.8.7929.tar.gz
    cd gpaw-setups-0.8.7929
    gunzip *.gz
    .bashrcに下記を書き込みます。
    export GPAW_SETUP_PATH=${HOME}/gpaw-setups-0.8.7929
    ${HOME}はgpaw-setups-0.8.7929.tar.gz を解凍したディレクトリ

  GPAWのインストール
    tar zxvf gpaw-0.8.0.7895.tar.gz
    cd gpaw-0.8.0.7895
    python setup.py install
[1] http://www.cybernet.co.jp/quantumwise/example/tips/tips-10.html
-------------------------------------------------------------------------------
■ ASE-GUI
□ 起動方法
  ・GPAWやASE用のGUIは下記の ag のコマンドで起動する。
  ag
□ 基本的な使い方
  ・File → Open → cif ファイルを指定すれば結晶構造が読み込まれる。
  ・左クリックで原子を選択 → Edit → Modify → Atom に元素名を入力 → OK で元素置換可能。
  ・右クリックしながらマウスを動かすと座標を回転できる。
  ・Calculate → Set Calculator → Density Functional Theory (GPAW) → Setup → k点を入力する(4 4 4)など → スピンが必要な場合は Spin polarized にチェックを入れる → OK → Apply → OK
  ・File → save → py やその他のファイルを選択 → save 
□ SCF計算
  ・Calculate → Energy and Forces → Run
□ 構造最適化
  ・Calculate → Energy Minimization → Run
□ 構造最適化
  ・Calculate → Scale system → Maximal scale factor を 0.150 など大きめにする → Number of steps も大きな値にする → Load optimal configuration にチェック → Run 
  ・原子位置も最適化するには、Atomic relaxations で On にする。
□ NEB計算[2](調査中)
  ・ ag --interpolate 3 initial.xyz final.xyz -o interpolated_path.traj
[1] https://wiki.fysik.dtu.dk/ase/ase/gui/gui.html#module-gui
[2] https://wiki.fysik.dtu.dk/ase/ase/gui/basics.html#neb-calculations
[3] http://sstxp.ee.ous.ac.jp/dokuwiki/doku.php
-------------------------------------------------------------------------------
※ このHPでの入力ファイルをコピー&ペーストしても動作しないことがあります。GPAWのHPにあるTutorialをコピー&ペーストして、それを元に下記の例のように書き換えていくのが良いでしょう。
■ 入力ファイル作成の基本
  from stuff import f, C
  print f(1, 2)
  print C(1).m(2)
  ・一行目は、stuff.py に記されている関数fとクラスCを用いるための宣言である。
  ・二行目は、関数fでの結果を表示する(print は表示させるためのコマンド)。
  ・三行目は、クラスCで記されている関数mでの結果を表示する。

■ 共通の入力ファイルの作成(原子座標の入力)
  ・terminalで ag を入力 → File → Open → cif ファイル →Open → File → save → py やその他のファイルを選択 → save
  ・terminalで emacs 先ほど作成したpyファイル を入力
  ・image を atom に書き換える(多くのtutorialの入力ファイルを僅かに修正するだけでよくなるため)。
  ・from ase import Atoms の次の行に from gpaw import GPAW を追加する。 
  ・Atoms や GPAW は ワイルドカードである * に書き換える(計算に慣れたらtutorialに書かれているのを真似して、理解を深める為にクラスを指定しよう)。 

□ スーパーセルの作成(pyを作成する場合に成功率が高いようだ)
  ag → setup → Bulk crystal → 結晶構造を入力 → Apply → python で得たデータをコピーする。
下記のように入力ファイルを作成する。x, y, z に格子定数を入力する。XAS計算の場合は from gpaw import GPAW を除いて、XASにある記述に従う。
from ase import Atoms
from gpaw import GPAW
import numpy as np
(python で得たデータをここにペーストする)
atoms.set_cell((x,y,z), scale_atoms=True)
atoms *= (3,3,3)

□ スーパーセルの作成(ag上で構造最適化する場合に良いだろう)
  ・ ag を起動する。
  ・ File → Open でcif ファイルを読み取るか、Setup → Bulk Crystal で結晶構造を入力(unit cells の数値を大きくするとスーパーセルが作成可能)→ Apply → OK
  ・ View → Repeat → Repeat atoms: の下に各方位で拡張する倍率を入力 → Set unit cell を押すとスーパーセルが形成される。(View → Zoom out で構造を見やすくしよう)
  ・元素置換したい場合は、 置換したい原子を図中から左クリックして選択し、Edit → Modify → atom の下の欄に置換する原子を入力する。画面が見づらい場合は、マウスを右クリックしながら上下左右に動かして図形を回転させるとよいだろう)。「正孔を導入する場合は何も変更は必要ない。しかし、単位胞に原子数が多くて混乱しやすい場合は置換する原子の原子番号が1つ大きなものを指定し、正孔を導入した擬ポテンシャルファイルのファイル名も原子番号を1つ大きくしたものにしよう」
  ・File → save → py で save する。

■ SCF計算
  共通の入力ファイルの作成後、そのファイルの一番最後に、下記を追加する。各結晶構造で書き換えるのはk点数を指定するkptsの数値。case.out のcase は好きな名称にすれば良い。その他は特に書き換えなくて良いだろう。下記のXXXに結晶名を入力。
calc = GPAW(h=0.18,
             kpts=(4,4,4),
             xc='PBE',
             txt='case.out',
             )
atom.set_calculator(calc)
e1 = atom.get_potential_energy()
print "total energy:", e1, "eV"
calc.write("XXX.gpw")

□ DOS(調査中)
  Tutorial[5] にある dos.py と pdos.py をそのまま用いればDOSの結果が得られそうだ。
[5] https://wiki.fysik.dtu.dk/gpaw/exercises/dos/dos.html

□ バンド分散(調査中)
  Tutorial[6]にあるPBE band structure が分かりやすい記述になっている。SCF計算の後に、「The band structure can be calculated like this:」と「and plotted like this:」をコピー&ペーストして実行すれば良さそうに見える。ただし、例はFCCのものなので、BCCやHCPも調べておく必要がある。
[6] https://wiki.fysik.dtu.dk/gpaw/tutorials/pbe0/pbe0.html
[7] https://wiki.fysik.dtu.dk/gpaw/tutorials/bandstructures/bandstructures.html#bandstructures
[8] https://wiki.fysik.dtu.dk/gpaw/exercises/band_structure/bands.html#band-exercise

■ 構造最適化 (原子位置の最適化)
  共通の入力ファイルの作成後、そのファイルの一番上に、下記を追加する。
from ase.optimize import QuasiNewton
  入力ファイルの一番最後に、下記を追加する。各結晶構造で書き換えるのはk点数を指定するkptsの数値。case.out のcase は好きな名称にすれば良い。その他は特に書き換えなくて良いだろう。他の原子位置の最適化の手法を採用したい場合は、入力ファイル中のQuasiNewtonをBFGSに換える。
calc = GPAW(h=0.18,
            kpts=(4,4,4),
            xc='PBE',
            txt='case.out',
            )
atom.set_calculator(calc)
relax = QuasiNewton(atom)
relax.run(fmax=0.05)
※ fmax = 0.05 程度にすると、原子位置の誤差は実験値に比べ、格子定数に対して1%程度になる。これ以内であれば対称性の高い構造に換えてみるという考えもどうだろうか?

■ NEB計算

[9] https://wiki.fysik.dtu.dk/ase/tutorials/neb/diffusion.html
[10] https://wiki.fysik.dtu.dk/gpaw/tutorials/neb/neb.html
[11] https://wiki.fysik.dtu.dk/gpaw/exercises/neb/neb.html

■ XAS計算
□ 正孔を導入した擬ポテンシャルの作成
  下記のXXの部分を元素名に換える。正孔1.0にしたい場合は、下記の0.5を1.0にする(Na, Mg は1.0の正孔を導入しようとすると失敗する)。hch1s はおそらく half core-hole 1s の略だと思われる。一方、fch1s は full core-hole 1s の略。
gpaw-setup -f PBE XX --name hch1s --core-hole=1s, 0.5
gpaw-setup -f PBE XX --name fch1s --core-hole=1s, 1.0
  成功すると元素名.hch1s.PBE という名称のファイルが作成される。(Terminal で ghost の記述が表示されればそのファイルは用いてはならない)

□ 基底状態の計算(調査中)
  共通の入力ファイルの作成後、fromのある最後の行に、下記を追加する。
from gpaw import GPAW, setup_paths
setup_paths.insert(0, '.')
   下記のXXXに結晶名を、Yに正孔を導入する原子に対する番号(上から順番{0から数える}に数えたときに空孔を導入する原子の番号)を記入する(Yは通常原点にある原子を指定することになる 0 の方が便利だろう)。入力ファイルの一番最後に、下記を追加する。(nbandsで指定する負の値が大きいほど非占有準位側の表示範囲が広くなる)
calc = GPAW( 
            nbands=-30,
            kpts=(4,4,4),
            h=0.2,
            xc='PBE',
            txt='case_xas.txt',
            setups={0: 'hch1s'})
atoms.set_calculator(calc)
e = atoms.get_potential_energy()
calc.write('case_xas.gpw')

□ dks_energy の計算
  共通の入力ファイルの作成後、fromのある最後の行に、下記を追加する。
from gpaw import GPAW, setup_paths
setup_paths.insert(0, '.')
   下記のXXXに結晶名を、Yに正孔を導入する原子に対する番号(上から順番{0から数える}に数えたときに空孔を導入する原子の番号)を記入する(Yは通常原点にある原子を指定することになる 0 の方が便利だろう)。入力ファイルの一番最後に、下記を追加する。
calc1 = GPAW(
             h=0.2,
             kpts=(4,4,4),
             txt='case_gs.txt',
             xc='PBE')
atoms.set_calculator(calc1)
e1 = atoms.get_potential_energy() + calc1.get_reference_energy()

calc2 = GPAW(
             h=0.2,
             kpts=(4,4,4),
             txt='case_exc.txt',
             xc='PBE',
             charge=-1,
             spinpol=True,
             occupations=FermiDirac(0.0, fixmagmom=True),
             setups={0: 'fch1s'})
atoms[0].magmom = 1
atoms.set_calculator(calc2)
e2 = atoms.get_potential_energy() + calc2.get_reference_energy()

print 'Energy difference' , e2 - e1

□ ブロードニングしたXASスペクトルのプロット
  先の計算で得られた Energy difference を dks_energy の値に記述する。
from gpaw import GPAW, setup_paths
from gpaw.xas import XAS
import pylab as plt
setup_paths.insert(0, '.')

dks_energy = 532.774  # from dks calcualtion

calc = GPAW('case_xas.gpw')
calc.set_positions()

xas = XAS(calc, mode='xas')
x, y = xas.get_spectra(fwhm=0.5, linbroad=[4.5, -1.0, 5.0])
x_s, y_s = xas.get_spectra(stick=True)

shift = dks_energy - x_s[0]  # shift the first transition

y_tot = y[0] + y[1] + y[2]
y_tot_s = y_s[0] + y_s[1] + y_s[2]

plt.plot(x + shift, y_tot)
plt.bar(x_s + shift, y_tot_s, width=0.001)
plt.show()

■ XAS計算(Hydock recursion method)
□ 正孔を導入した擬ポテンシャルの作成
   下記のXXの部分を元素名に換える。正孔1.0にしたい場合は、下記の0.5を1.0にする(Na, Mg は1.0の正孔を導入しようとすると失敗する)。hch1s はおそらく half core-hole 1s の略だと思われる。一方、fch1s は full core-hole 1s の略。
gpaw-setup -f PBE XX --name hch1s --core-hole=1s, 0.5
gpaw-setup -f PBE XX --name fch1s --core-hole=1s, 1.0
   成功すると元素名.hch1s.PBE という名称のファイルが作成される。(Terminal で ghost の記述が表示されればそのファイルは用いてはならない)

□ 基底状態の計算
  共通の入力ファイルの作成後、fromのある最後の行に、下記を追加する。
from gpaw import GPAW, setup_paths
setup_paths.insert(0, '.')
  下記のXXXに結晶名を、Yに正孔を導入する原子に対する番号(上から順番{0から数える}に数えたときに空孔を導入する原子の番号)を記入する(Yは通常原点にある原子を指定することになる 0 の方が便利だろう)。入力ファイルの一番最後に、下記を追加する。
name= 'XXX_hch'
calc = GPAW( h=0.2, kpts=(4,4,4,), txt = name + '.txt', xc='PBE', setups={Y:'hch1s'})
atom.set_calculator(calc)
e = atom.get_potential_energy(calc)
calc.write(name+'.gpw')

□ recursion coefficients の計算(調査中)
from ase import Atoms
import numpy as np
from gpaw.xas import RecursionMethod
from gpaw import GPAW, setup_paths
setup_paths.insert(0, '.')
  下記のXXXに結晶名を記入する。kptsの数値にk点を指定する。入力ファイルの一番最後に、下記を追加する。
name='XXX_hch'
calc = GPAW(name +'.gpw', kpts=(6,6,6), txt=name + '_rec.txt')
calc.set_positions()
r = RecursionMethod(calc)
r.run(600)
r.run(1400, inverse_overlap='approximate')
r.write(name + '_600_1400a.rec', mode='all')

□ ブロードニングしたXASスペクトルのプロット(調査中)
下記のXXXに結晶名を記入する。ZZZは上記にある 「dks_energy の計算」で得る。
import sys
import numpy as np
import pylab as plt
from ase import Atoms
from gpaw.xas import RecursionMethod
from gpaw import GPAW, setup_paths
setup_paths.insert(0, '.')

sys.setrecursionlimit(10000)
name="XXX_hch_600_1400a.rec"
x_start=-20
x_end=100
dx=0.01
x_rec = x_start + np.arange(0, x_end - x_start, dx)
r = RecursionMethod(filename=name)
y = r.get_spectra(x_rec, delta=0.4, fwhm=0.4)
y2 = sum(y)
plt.plot(x_rec + ZZZ, y2)
plt.show()

■ EELS計算
□ EELSスペクトル
  Tutorial[4] の記述が非常に良く出来ている。A simple startup: Bulk aluminum の例をコピー&ペーストして、atoms の行を消して、ag で得られる原子構造のデータに書き換えれば良い。この他に書き換える部分はkpts での数値、そして、 w=np.linspace (0, 40, 401)と大きめの範囲にすれば良いだろう。計算が終了したら、下記のように例にあるpyをコピー&ペースとしてplot_EELS.py として保存し、python plot_EELS.py で plot すればよい。
from pylab import *
import numpy as np

d = np.loadtxt('EELS.dat')
plot(d[:,0], d[:,1], '-k')
show()

□ q点の位置を徐々に変えた場合のEELSスペクトル(調査中)
  Tutorial[4] の記述が非常に良く出来ている。A more sophisticated example: graphite の例をコピー&ペーストして、atoms の行を消して、ag で得られる原子構造のデータに書き換えれば良い。この他に書き換える部分はxc="PBE"、kpts での数値、nbandsでの数値、そして、convergence の数値だろう。計算が終了したら、python plot_EELS.py でplot すればよい。
[4] https://wiki.fysik.dtu.dk/gpaw/tutorials/dielectric_response/dielectric_response.html#example-1-optical-absorption-of-bulk-silicon

■ TDDFT計算

■ STM計算

■ GLLB-sc計算
  共通の入力ファイルの作成後、そのファイルの一番最後に、下記を追加する。各結晶構造で書き換えるのはk点数を指定するkptsの数値。case.out のcase は好きな名称にすれば良い。その他は特に書き換えなくて良いだろう。
calc = GPAW(h=0.16,
            kpts=(10,10,10),
            xc='GLLBSC',
            txt='case.out',
            eigensolver='cg',
            occupations=FermiDirac(width=0.05),
            )
atom.set_calculator(calc)
atom.get_potential_energy()
response = calc.hamiltonian.xc.xcs['RESPONSE']
response.calculate_delta_xc()
EKs, Dxc = response.calculate_delta_xc_perturbation()
Gap = EKs+Dxc
print "Calculated band gap:", Gap

■ PBE0計算

■ RPA計算

■ GW計算

[1] http://www.csc.fi/english/csc/courses/archive/gpaw-2008-01 
[2] https://wiki.fysik.dtu.dk/gpaw/tutorials/tutorials.html
[3] http://www.quantumwise.co.jp/example/tips/tips-10.html 
-------------------------------------------------------------------------------
GPAW : https://wiki.fysik.dtu.dk/gpaw/index.html
チュートリアルを参考に入力用のスクリプト case.py を作成し、python case.py で計算する。
並列計算はmpirun -np 4 gpaw-python case.py で計算できる。
[1] http://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&frm=1&source=web&cd=4&cad=rja&ved=0CDwQFjAD&url=http%3A%2F%2Fwww.csc.fi%2Fenglish%2Fcsc%2Fcourses%2Farchive%2Fmaterial%2FGPAW-2008-1%2Fparallel&ei=SplUUNHADsf1mAWspYGoCQ&usg=AFQjCNHKQiO97RBB51q0gMeohPddYO16Og&sig2=D7IZ_TwlI570eKmg51yuYQ
-------------------------------------------------------------------------------
アクセス数
ページビュー数