====== Dacapoマニュアル ====== ===== はじめに ===== このドキュメントではASE/Dacapoの使い方を簡単に説明する。ただしASEは最新の第三版ではなく第二版に準拠しているので注意が必要である。pythonスクリプトの中ではASE2.xのクラス名およびクラスメソッドは''ListOfAtoms''のように英単語頭大文字の文字列で、クラス内の予約語およびpython予約語は小文字英単語で、ASEクラスオブジェクトのインスタンスと物理的・数学的に意味のある変数は''calculator2''のように英単語+数字で、forループのカウンタのような制御用変数やその場限りの作業用変数は意味のない英文字列で表すことにする。 ===== 原子配置の設定 ===== ==== ASE ListOfAtoms オブジェクト ==== 平行六面体の中にいくつかの原子を配置したものが''ListOfAtoms''オブジェクトだと思えば良い。Dacapoでの計算（周期境界条件を用いる計算）では、この平行六面体を上下左右前後に並べて計算をする。従って結晶のようにもともと周期性を持っているものを計算する際にはこの平行六面体として単位格子をとり、単位格子内の原子だけをその中に配置すれば良いし、分子の計算をする際にはこの平行六面体を周期的に並べても分子と分子の間に働く相互作用ができるだけ小さくなるように平行六面体の大きくとっておけば良い。 例1. 一酸化炭素分子を大きめの箱の中心に配置する。 from ASE import Atom, ListOfAtoms atom1 = Atom('C', (0.50, 0.50, 0.50)) atom2 = Atom('O', (0.62, 0.50, 0.50)) model1 = ListOfAtoms([atom1, atom2]) model1.SetUnitCell((10.0, 10.0, 10.0)) 10Å×10Å×10Åの立方体の箱の中心に炭素を 例2. アルミニウムの結晶 from ASE import Atom, ListOfAtoms model1 = ListOfAtoms([Atom('Al')]) vector1 = (0.0, 2,5, 2,5) vector2 = (2.5, 0.0, 2,5) vector3 = (2.5, 2,5, 0.0) model1.SetUnitCell((vector1, vector2, vector3)) ListOfAtomsクラスの詳細はASE マニュアルを参照のこと。 ASEでは全エネルギー等を計算するためのオブジェクトをcalculatorオブジェクトと称している。CAMdではここで取り上げるDacapoに加えてVASPやAbinitなど有償無償を含めて様々なcalculatorが使用できるようだが、ここで利用するcalculatorはDacapoのみであるので、以下Dacapoに絞って記述する。定義済みのListOfAtomsオブジェクトであるインスタンス''model1''とDacapoオブジェクトのインスタンス''calculator1''とを結合させるためには''SetCalculator''メソッドを用いて、 model1.SetCalculator(calculator1) とする。 ====== Dacapo ASE calculator interface ====== Dacapoを使って計算する際に、最低限必要Dacapoオブジェクトに与えなければならない情報は計算するコーン・シャム方程式の解の個数である。これは計算する系の価電子数の二分の一以上でなければならない。 上記のように定義されたListOfAtomsクラスのインスタンス''atoms1''とDacapoクラスのインスタンス''calculator1''を結合させ、計算すべき解の個数6を与えるには以下のように記述する。 from Dacapo import Dacapo calculator1 = Dacapo() calculator1.SetNumberOfBands(6) atoms1.SetCalculator(calculator1) Dacapoオブジェクトに与える他の計算条件は適当なデフォルト値が設定されているが、これを変更するためには二つの方法がある。 一つは calculator = Dacapo(nbands = 6) もしくは calculator = Dacapo( ... nbands = 6, ... planewavecutoff = 340) のようにインスタンス生成時にキーワードとともに指定する方法である。一般的な記述は インスタンス名 = Dacapo(キーワード = 値, キーワード = 値, ...) のようになる。もう一つは calculator = Dacapo() calculator.SetNumberOfBands(6) calculator.SetPlaneWaveCutoff(340) のようにSet...系のクラスメソッドを使用する方法である。 一般的な記述は インスタンス名 = Dacapo() インスタンス名.Setメソッド名(値) の様になる。以下に計算する解の個数以外の設定できる計算条件を、設定に使用するクラスメソッド名、インスタンス生成時に指定するキーワード、デフォルト設定値を クラスメソッド名(キーワード = デフォルト値) の形式で列挙する。#以下はその値の単位である。 SetPlaneWaveCutoff(planewavecutoff = 300) #[eV] SetDensityCutoff(densitycutoff = 300) # [eV] SetXCFunctional(xc = ”PW91”) SetBZKPoints(kpts = ((0,0,0))) SetSpinPolarized(spin = False) SetNetCDFFile(out = ”out”) SetTxtFile(txtout = ”out.text”) DipoleCorrection(dipole = False) SetOccupationStatistics(method = FermiDirac) SetElectronicTemperature(temp = 0.1) # [eV] ===== 計算する状態数 ===== ''SetNumberOfBands(nbands)''によって設定する。設定値は自然数である。コーン・シャム方程式を自己無撞着場中の一つの電子に対する波動方程式と見なすと、これは計算すべき一電子状態の数ということになる。したがって計算対象のListOfAtomsオブジェクトに含まれる価電子の総数をスピンの多重度、2で割ったものである。ただし、大きな分子や結晶ではフェルミ準位付近に多くの状態が存在することが多く、状態の占有数を単純な０（非占有状態）、１（スピン分極を考慮する計算の際の占有状態）、２（スピン分極を考慮しない計算の際の占有状態）といった整数値だけでなく後述する''SetOccupationStatistics''メソッドで指定される方法によって適切な実数値で設定する。そのためには価電子の総数から期待される値よりも若干大きく設定する。これは自己無撞着計算の繰り返し中に特定の対称性の状態が計算されたり、計算されなかったりといったことを避けるためにも必要である。適切な値を設定できたかどうかは、計算結果の非占有状態の占有数がゼロになっているかどうかで確認する。また、密度汎関数理論では励起状態のエネルギーの計算は保証されないこと、正確にはコーン・シャム方程式の固有値を一電子エネルギーと見なすことはできないことから非占有状態の固有値に物理的意味を求めることはできない。例えば半導体や絶縁体のエネルギーギャップが実験結果と合わないことは有名である。しかしながら、非占有状態の固有値であってもその変化、例えば波数に対する変化（分散関係）や不純物位置に対する不純物準位の変化などが実験結果とよく一致するのも事実であり、非占有状態の固有値を計算する場合も多い。そのような場合は必要な数だけの状態数を''SetNumberOfBands''で指定する。 ===== 電子系の温度 ===== SetOccupationStatistics(method)を用いると1電子状態の占有数の決め方を変更することができる。Dacapoでは次の二種類の方法から選択できる。 ==== FermiDirac ==== 
 Fermi-Dirac分布を用いる。電子系の温度はSetElectronicTemperature(temp)メソッドで与える。デフォルトでは 0.1 eV　に設定されている. ==== MethfesselPaxton ==== 
 1次のMethfessel-Paxton法((Phys. Rev. **B40** (1989) 3616.))に従って占有数を決定する。 ===== カットオフエネルギー ===== 平面波展開で用いる波数の大きさの上限をエネルギー（eV）の単位で与える。波動関数(planewavecutoff)と電荷密度(densitycutoff)、それぞれに異なるカットオフエネルギーを与えることができる。それぞれ、SetPlaneWaveCutoff(planewavecutoff)メソッドとSetDensityCutoff(densitycutoff)メソッドで定義する。 通常、planewavecutoffとdensitycutoffは同じ値にとられるが、Cuのようなd電子などの局在性の強い原子を含む場合はdensitycutoffをより大きな値にとらなければならない。これは、Cuなどに使われているウルトラソフト擬ポテンシャルが非常にソフトなので、その結果得られる擬波動関数も非常にソフトなのに対して、電荷密度の方は本物の波動関数の局在性のために擬波動関数ほどにはソフトではないからである。 ===== 出力ファイルの指定 ===== Dacapoの出力ファイルには電荷密度など全ての情報が含まれるバイナリ形式のNetCDFファイルと サマリだけがテキスト形式で出力されるASCIIファイルの二つが指定できる。デフォルトではそれぞれ 「out.nc」と「out.txt」である。これを変更するにはSetNetCDFFileメソッドとSetTxtFileメソッドを 用いる。 SetNetCDFFile(filename) SetTxtFile(filename) ===== 交換相関汎関数 ===== SetXCFunctionalメソッドにより計算に使用する交換相関ポテンシャルを選択することができる。 SetXCFunctional(exc) 引数excは以下のものから選択する。 ; LDA : Vosko Wilk Nusair LDA-parametrization. ; PW91 : Perdew Wang 91 GGA-parametrization. ; PBE : Perdew Burke Ernzerhof GGA-parametrization. ; RPBE : Revised PBE GGA-parametrization, Hammer, Hansen and Nørskov. 現在の交換相関ポテンシャルが何かを知るためにはGetXCFunctionalメソッドを用いる。また、GetXCEnergiesメソッドを用いることで、計算に用いた交換相関ポテンシャルとは異なる相関交換汎関数を使って相関交換エネルギーを見積もることもできる。 GetXCEnergies(xcname) 引数 xcname は上の (PBE,RPBE,PW91,VWN)から選択できる。ただし、計算に用いた交換相関ポテンシャル以外で求めた相関交換エネルギーは自己無撞着な計算にはなっていないことに留意せよ。 ===== スピン分極計算 ===== ''SetSpinPolarized(True)''を指定することによりスピン分極をさせた計算を行うことができる。この場合、 from ASE import Atom, ListOfAtoms atom1 = Atom('O', (0, 0, 0), magmom = 1.0) atom2 = Atom('O', (0, 0, 1.2), magmom = 1.0) molecule1 = ListOfAtoms([atom1, atom1], cell = (4, 4, 4)) のように、ListOfAtomsオブジェクトを構成するAtomオブジェクトに磁気モーメントの初期値を設定する必要がある。 これは from ASE import Atom, ListOfAtoms atom1 = Atom('O', (0, 0, 0)) atom2 = Atom('O', (0, 0, 1.2)) molecule1 = ListOfAtoms([atom1, atom1], cell = (4, 4, 4)) molecule1.SetMagneticMoments((1.0, 1.0)) のように、後からListOfAtomsオブジェクトに対して''SetMagneticMoments''メソッドを用いて設定したり、 from ASE import Atom, ListOfAtoms atom1 = Atom('O', (0, 0, 0)) atom2 = Atom('O', (0, 0, 1.2)) atom1.SetMagneticMoment(1.0) atom2.SetMagneticMoment(1.0) molecule1 = ListOfAtoms([atom1, atom1], cell = (4, 4, 4)) のように、個々のAtomオブジェクトに''SetMagneticMoment''メソッドを用いて設定してもよい。 この例では、どれも一つの酸素原子に対して1ボーア磁子、すなわち酸素分子に対して2ボーア磁子の磁気モーメントを設定している。実際にはそれで正しい結果なのだが、上の例で atom1.SetMagneticMoment( 1.0) atom2.SetMagneticMoment(-1.0) などとすると全磁気モーメントゼロの反強磁性の状態から計算を始めることもできる。この場合も最終的には全磁気モーメントが2ボーア磁子になる正しい計算結果が得られるが、自己無撞着な解に到達するまでの繰り返し回数は大きくなる。また、以下の例のようにキーワード''fixmagmom''を用いれば磁気モーメントをある値(ここでは2ボーア磁子)に固定して計算をすることもできる。 calc1.SetSpinPolarized(spinpol = True, fixmagmom = 2.0) ===== K点 ===== ''SetBZKPoints(kpts)''メソッドを使って計算に用いるブリリュアンゾーン内のK点を定義できる。 省略時のデフォルト値は(0,0,0)のみの点、すなわちΓ点一点だけの計算である。 これは''SetBZKpoint((1,1,1))''を指定するのと同等である。 再構成表面を扱う場合など実空間のユニットセルの体積が大きい、すなわち逆空間のウィグナー・サイツセルであるブリリュアンゾーンの体積が小さいときにはΓ点一点だけでもそこそこの計算になる。一般にはkptsは逆格子の単位格子を単位として表わした座標をpythonのlistの形で表現したものであるが、''SetBZKpoint((5,4,3))''等とすることによって後述する5×4×3等のMonkhorst-Packのスペシャルポイント((H.J.Monkhorst and J.D.Pack, Physical Review B, vol. 13, page 5188, 1976.))を用いることもできる。また、表面系を計算する場合にはChadi-Cohenのスペシャルポイント((Chadi and Cohen, Phys. Review B., vol 8,page 5747))を使うことも多い。 ==== MonkhorstPack ==== calc1.SetBZKPoints((n1, n2, n3)) と設定することによって''n1''×''n2''×''n3''のMonkhorst-Packの定義した三次元の特別なk点のセットを用いて計算する。このとき各逆格子ベクトルに沿っての大きさはそれぞれn1, n2, n3で与えられる。 ==== ChadiCohen ==== ASEのUtilitiesモジュールにはChadiとCohenによって定義された主に表面・薄膜系で用いられる二次元のK点が名前付きで定義されている。これは from ASE.Utilities.ChadiCohen import CC6_1x1 calc.SetBZKPoints(CC6_1x1) のように使用することができる。定義済みのChadi-CohenのK点は ''CC6_1x1'', ''CC12_2x3'', ''CC18_sq3xsq3'', ''CC18_1x1'', ''CC54_sq3xsq3'', ''CC54_1x1'', ''CC162_sq3xsq3'', ''CC162_1x1'' などである。詳細は、''ASE.Utilities.ChadiCohen''のソースコードを参照せよ。 実際にどのようなK点の座標が設定されたかは''GetBZKPoints''メソッドでその座標のリストを得ることができる。計算対象の対称性によって計算するK点を減少させた場合には、''GetIBZKPoints''メソッドで対称性を考慮して独立なK点のリストを得ることができる。 ===== 擬ポテンシャル ===== Dacapoで使用される擬ポテンシャルの一覧は[[http://www/|Pseudopotential Library]]を参照のこと。 標準で用意されている擬ポテンシャルのファイルの場所(パス名)はGetPseudoPotentialメソッドで参照できる GetPseudoPotential(elementnumber) 引数は元素の原子番号である。 このファイル以外のファイルから擬ポテンシャルを読み込んで計算を実行するにはSetPseudoPotentialメソッドを用いる。 ``SetPseudoPotential(elementnumber, path)``: 引数elementnumberは元素の原子番号、引数pathは擬ポテンシャルファイルのファイル名（フルパス名）である。 ====== 電荷密度の混成 ====== 自己無撞着計算の説明の際、「新しく得られた波動関数から電荷密度を計算し、これを再度コーン・シャム方程式の交換相関ポテンシャルに代入して、・・・」のように述べることが多いが、実際の実装に際しては毎回電荷密度を完全に置き換えてしまうのではなく古い電荷密度と新しく計算された電荷密度を適当な比率で平均化して次の繰り返し計算を行うことがよく行われる。この混合を全く行わないこともできる。この場合には自己無撞着な計算にはならない訳だが、少ないK点に対して自己無撞着な電荷密度を求めておいて、分散関係を描くために細かくK点を採り直した計算をするような用途に用いられる。 また、平均化の方法も指定できる。 The method ``SetChargeMixing(value)`` can be used to enable or disable charge density mixing. Using ``SetChargeMixing(True)`` the program will use Pulay mixing for the density (default). You can change the Pulay density mixing coeffficient (default 0.1) and the Pulay mixing history (default 10) using:: calc.SetNetCDFEntry("ChargeMixing",att="Pulay_DensityMixingCoeff",value=0.1) calc.SetNetCDFEntry("ChargeMixing",att="Pulay_MixingHistory",value=10) Using ``SetChargeMixing(False)`` correspond to running a calculation using the Harris-Foulkes functional using the input density. See also the `Harris calculation example`_. The method ``SetKerkerPreconditioning(value)`` can be used to set Kerker preconditioning for the density mixing. For value=True Kerker preconditiong is used, i.e. q is different from zero, see eq. 82 in Kresse/Furthmller: Comp. Mat. Sci. 6 (1996). The value of q is fix to give a damping of 20 of the lowest q vector. For value=False q is zero and mixing is linear (default). ====== 固有値の計算法 ====== ''SetEigenvalueSolver(method)''により固有値の計算アルゴリズムを指定できる。現在''method''に指定できるのは次の二つである。 ==== eigsolve ==== Block Davidson algorithm (Claus Bendtsen et al). これが既定値である。 ==== rmm-diis ==== Residual minimization method (RMM), using DIIS (direct inversion in the iterate subspace). The implementation follows closely the algorithm outlined in Kresse and Furthmuller, Comp. Mat. Sci, III.G/III.H ====== 双極子補正 ====== 表面や薄膜を計算する目的でいわゆるスーパースラブ模型を用いて、半導体や絶縁体などの誘電体や、金属であっても電位陰性度やイオン化傾向の大きな原子の吸着を扱う場合、もしくはそもそも有極分子の吸着を扱おうとすると必然的にスラブ自身が分極し付加的な静電ポテンシャルが生じることがある。 DipoleCorrection(MixingParameter　=　0.2, ... InitialValue = 0.0, ... AdditiveDipoleField = 0) と指定することにより、このスラブ間の静電的な相互作用を打ち消すような双極子補正を行うことができる。現在この機能は第三の単位格子ベクトルに沿った方向でのみ有効である。従って、用いるスーパースラブのスーパーセルを適切に選ばなければならない。また、双極子より高次の静電多重極子に対する補正はまだ実装されていない。既定値では真空中の電場が不連続になる位置はスラブの両側のどの原子からも最も遠くなる点が選ばれる。 キーワードAdditiveDipoleFieldにより付加的な電気二重層に相当する一定の静電場を印可する。この電場が双極子補正に加えられる。電場の不連続性は動的双極子補正の場所に追随する。（ここのところ未確認） 詳細は“Dipole correction example”, “workfunction pdf document”, “dacapo netcdf manual”を参照。 双極子補正で付加される電場ではない一定の外場を印可するにはNetCDF entryの ExternalDipolePotentialに値を設定する。例えば、 calc.SetNetCDFEntry("ExternalDipolePotential",　value　=　0.5) により、 0.5[eV/A]の外場を与えることができる。 ====== 対称性 ====== ; SetSymmetryOff : 対称性を考慮すると同等なK点に対しても計算を行う。 ; SetSymmetryOn : 対称性を考慮して計算するK点の数を減らす。 ====== Electrostatic decoupling ====== SetElectrostaticDecoupling(numberofgaussians=3,ecutoff=100)'' Add electrostatic decoupling (Jan Rossmeisl). Decoupling activates the three dimensional electrostatic decoupling. Based on paper by Peter E. Bloechl: JCP **103** page7422 (1995). ====== 応力計算 ====== 標準では応力計算は実行されない。 ''CalculateStress()'' メソッドを使用することで応力計算が実行される。 ====== 局所状態密度 ====== Dacapoでは原子軌道に射影した状態密度を計算することができる。典型的なpythonコードは次のようなものである。 atoms = ListOfAtoms(..) calc = Dacapo(..) calc.CalculateAtomicDOS()   atoms.SetCalculator(calc)   energy = atoms.GetPotentialEnergy() dos = calc.GetLDOS(atoms = [1], angularchannels = ["s"]) dos.GetPlot() この例では最初の原子のs軌道に射影した状態密度を与える。（原子の番号は1から順番に振られるが、必ずしもappendした順番ではないのでどの原子が何番の原子なのかは出力ファイルで事前に確認する必要がある The method:: * ''CalculateAtomicDOS(energywindow = None)''\\ tells the fortran program to calculate the local density of states projected onto atomic orbitals, using ''energywindow'' for calculating the energy resolved LDOS in a specific range (must be within the eigenvalue spectrum). * ''GetLDOS(**keywords )''\\ returns an instance of ''AtomProjectedDOSTool'' (see below). Keyword for the ''GetDOS'' method: * atoms\\ List of atoms numbers (default is all atoms) * angularchannels\\ List of angular channel names.\\ The full list of names are: ''s'', ''p_x'', ''p_y'', ''p_z'', ''d_zz'', ''dxx-yy'', ''d_xy'', ''d_xz'', ''d_yz''. ''p'' and ''d'' can be used as shorthand for all p and all d channels respectively. (default is all d channels) * spin\\ List of spins (default all spins) * cutoffradius\\ ''short'' or ''infinite''. * For ''cutoffradius = 'short''' the integrals are truncated at 1 Angstrom. * For ''cutoffradius = 'infinite''' the integrals are not truncated. (default 'infinite') The resulting DOS are added over the members of the three list (atoms,angularchannels and spin). GetDOS returns a instance of the class AtomProjectedDOSTool (Thanks to John Kitchen for providing grace plot and band moments methods). Methods for AtomProjectedDOSTool: * ''GetPlot''\\ returns a GnuPlotAvartar for the energy resolved DOS A parent can be given too combine plot. * ''GetIntegratedDOS''\\ returns the integral up to the fermi energy * ''GetData''\\ returns the energy resolved DOS * ''GetBandMoment(moments)''\\ returns the moments of the projected DOS. * ''GetBandMoment(0)''\\ return the integrated DOS. * ''GetBandMoment(1, 2)''\\ returns the first and second moment of the projected DOS (center and width) * ''SaveData(filename)''\\ saves the data to the filename * ''GetGracePlot''\\ makes a GracePlot (if GracePlot is installed) The eigen values are not relative to the Fermi level, you can get the fermi level for the calculation using calc.GetFermiLevel() ===== Multicenter orbitals ===== You add the multicenter orbitals using:: mcenters = [(0,0,0,1.0),(1,0,0,-1.0)] calc.SetMultiCenters(mcenters) this will define an atom1_s-atom2_s orbital. The general form is [{(atomnr,l,m,weight)}] You can retrieve the centers using centers = calc.GetMultiCenters() This will return a list of instances of the class MultiCenterProjectedDOSTool. This class has the same method as described above, ie. ''GetBandMoment'', ''GetCutoff'', ''GetData'', ''GetDescription'', ''GetEFermi'', ''GetGracePlot'', ''GetIntegratedDOS'', ''GetPlot'', ''GetSpin'' and ''SaveData''. You can give the following keywords to ''SetMultiCenters'' calc.SetMultiCenters(mcenters, energywindow = (-15, 5), energywidth = 0.2, numberenergypoints = 250, cutoffradius = 1.0) The direction cosines for which the spherical harmonics are set up are using the next different atom in the list (cyclic) as direction pointer, so the z-direction is chosen along the direction to this next atom. At the moment the rotation matrices is only given in the text file, you can use grep 'MUL: Rmatrix' out_o2.txt to get this information. Below is a small example, which defines a s-orbital on atom 0:: #!/usr/bin/env python from Dacapo import Dacapo from ASE import Atom, ListOfAtoms atoms = ListOfAtoms([Atom('O', (2, 2, 2), magmom=1.0), Atom('O', (3.1, 2, 2), magmom=1.0)], cell=(4, 4, 4)) calc = Dacapo(planewavecutoff=350, # in eV nbands=8, # 5 extra empty bands out='o2.nc', spinpol=True ) atoms.SetCalculator(calc) mcenters = [] mcenters.append([(0,0,0,1.0), (1,0,0,0.0)]) # [{(atomnr,l,m,weight)}] calc.SetMultiCenters(mcenters, energywidth=0.1) energy = atoms.GetPotentialEnergy() centers = calc.GetMultiCenters() center1 = centers[0] print 'description:' center1.GetDescription() print 'energy resolved dos:' for d in center1.GetData(): print d[0],d[1] ====== 収束判定条件の設定 ====== 繰り返し計算で自己無撞着な解を得る場合、入力と出力、もしくは繰り返しの一つ前の結果と現在の結果を比較し、変化がなくなったら計算を終了する。しかしながら、実数値が数値計算上一致することはありえない（もし一致していたらプログラムのバグを疑った方がよい）ので ある範囲内で一致しているかどうかで解が収束したかどうかを判定する。DacapoではSetConvergenceParametersメソッドを用いて収束判定パラメータを設定できる。 SetConvergenceParameters(keywords=..., ..., ) 引数keywordsには次のようなパラメータを設定する。 ; energy : 全エネルギー (デフォルト値: 1e-5) ; density : 電荷密度 (デフォルト値: 0.0001) ; occupation : 状態占有数 (デフォルト値: 0.001) SetConvergenceParameters( ... energy = 0.00001, ... density = 0.0001, ... occupation = 0.001) と設定するとデフォルトの設定と同じパラメータを設定したことになる。 GetConvergenceParametersメソッドを用いると、現在の設定を確認することができる。 GetConvergenceParameters() ====== NetCDF specific methods ====== Two methods are defined for general reading and writing of ''NetCDF'' entries. This provides a method for reading and writing any NetCDF variable defined in the dacapo netcdf manual, so that any settings for the dacapo fortran program, not defined by the standard sets of method from the DFT Calculator or from the extra Dacapo method described above, can be defined. * ''SetNetCDFEntry(variable, att = None, value = None)''\\ Define a general netcdf entry for this calculation. If the netcdf entry is allready defined, the value are modified. * ''GetNetCDFEntry(variable, att = None)''\\ Returns the value of NetCDF entry '''' or NetCDF attribute ''att'', if specified. calc = Dacapo() calc.SetNumberOfBands(10) print calc.GetNumberOfBands() calc = Dacapo() calc.SetNetCDFEntry('ElectronicBands', att = 'NumberOfBands', value = 10) print calc.GetNetCDFEntry('ElectronicBands', att = 'NumberOfBands') ====== 計算結果の再利用 ====== 一旦Dacapoの計算が終われば、NetCDFの出力ファイルには自己無撞着な電荷密度をはじめとして計算できるすべての情報が書き込まれている。この情報を再利用するにはDacapoモジュールのReadAtomsメソッドを用いてListOfAtomsクラスのインスタンスを生成しGetCalculatorメソッドを用いてDacapoクラスのインスタンスを得る。このDacapoクラスのインスタンスに対してGet....系のクラスメソッドを用いることにより、計算結果と同時に計算時に設定した計算条件も知ることができる。 from Dacapo import Dacapo atoms = Dacapo.ReadAtoms(filename = 'old_output_file.nc') calc = atoms.GetCalculator() print calc.GetEigenvalues() ただし、このようにして生成したListOfAtomsクラスやDacapoクラスのインスタンスそれ自身を再利用してはならない。 ===== 計算条件の問合せ ===== ^ 問合せメソッド ^ 返値 ^ | ''GetNumberOBands()'' | （整数）計算した状態の数 | | ''GetXCFunctional()'' | （文字列）用いた汎関数の種類、('LDA','PBE' などを返す) | | ''GetBZKPoints()'' | （座標のリスト）計算したブリリュアンゾーン内のK点、単位は逆格子の基底 | | ''GetSpinPolarized()'' | （真偽値）スピン分極の計算であればTrue、非磁性の計算であればFalse | | ''GetMagneticMoment()'' | （実数）全磁気モーメント、単位はボーア磁子 | | ''GetIBZKPoints()'' | （座標のリスト）計算したブリリュアンゾーン内の__独立な__K点、単位は逆格子の基底 | | ''GetIBZKPointWeights()'' | （実数のリスト）ブリリュアンゾーン内のK点についての積分に用いる重み | ===== 計算結果の取得 ===== ==== GetEigenvalues(kpt = 0, spin = 0) ==== K点kpt、スピンspin（省略時kpt = 0, spin = 0）に対して計算されたエネルギー固有値を実数のリストとして返す。エネルギー固有値についてはフェルミエネルギーからの相対値がeV（エレクトロンボルト）の単位で出力される。すなわち0eVのエネルギー固有値の状態がフェルミ準位。 ==== GetOccupationNumbers(kpt = 0, spin = 0) ==== K点kpt、スピンspin（省略時kpt = 0, spin = 0）に対して計算された状態の占有数を実数のリストとして返す。スピン分極を有効にした計算では０から１までの、スピン分極を考慮しない計算では０から２までの値である。この占有数と前述エネルギー固有値の積の総和が全エネルギーである。 ==== GetWavefunctionArray(band, kpt=0, spin=0) ==== Return array of wavefunction values. The wavefunction are returned for the band band, K-point kpt and spin spin. The wavefunction is returned as a three-dimensional Numerical array, corresponding to the grid use in the calculation. ==== GetWavefunctionArrays(kpt=0, spin=0) ==== Return array of wavefunction values for all bands for the given K-point kpt and the given spin spin. The wavefunctions are returned as a three-dimensional Numerical array. ==== GetDensityArray(spin=0) ==== Return array of density values for the given spin spin. ====== 波動関数と電荷密度 ====== 密度汎関数理論では基底状態のエネルギーと[[電子密度|電荷密度]]の関係しか保証しないので、コーン・シャム方程式の固有値と固有関数を一電子エネルギーや電子の波動関数とを安易に対応させてはならない。しかしながら、コーン・シャム方程式は自己無撞着な交換相関ポテンシャルを含んだ一電子のシュレディンガー方程式と同じ形をしているので、その固有値と固有関数は個々の電子のエネルギーや波動関数に対する物理的な洞察を与えてくれるのも確かである。 [[波動関数]]の情報はDacapoの計算結果が出力されたNetCDFファイル（拡張子.nc）に含まれているがこのファイルはテキストエディタで読むことができないファイル（バイナリファイル）なので、直接波動関数の情報を得ることはできない。そもそも波動関数は、３次元のしかも複素関数であるから、もし情報を得られたとしてもそれを可視化するためには何かしらのアプリケーションプログラムが必要になる。 ASEでは波動関数もしくは電荷密度の情報を扱うために二つの方法が用意されている。一つはVTK (Visual Tool Kit) を用いて端末上で直接可視化する方法である。これは、X-windowシステムが動作している端末で作業している際に便利な方法である。もう一つは著名な分子軌道計算パッケージであるGaussianの出力ファイ形式として知られる[[cube]]ファイル形式で三次元の関数の情報を出力しておき、[[VISTA]]や[[VMD]]などの[[cube]]ファイルが扱えるアプリケーションソフトウェアで処理をする方法である。これは、出力結果を一旦MacintoshやWindows PCにダウンロードしてから処理をする際に便利である。 VTKを用いて波動関数を描画するには、 from Dacapo import Dacapo from ASE.Visualization.VTK import VTKPlotWaveFunction atoms = Dacapo.ReadAtoms(filename = 'out.nc') VTKPlotWaveFunction(atoms) NetCDFファイルからcubeファイルを作成するには、別稿に述べるpythonスクリプトを用いる。 ====== 構造最適化と分子動力学 ====== Dacapoプログラムを実行することにより、全エネルギーや電荷密度、エネルギー固有値や波動関数だけでなく、各原子に働く力（ヘルマン・ファイマン力）を計算することができる。出力ファイルに出力された力や、GetCartecianForcesなどのメソッドを使って読み出した力を用いて手動もしくは自作のプログラムで原子を動かすことも可能だが、ASEに含まれるクラスオブジェクトを利用すると、簡便に原子にかかっている力に応じて原子を動かすことができる。一般に原子にかかっている力に応じて、主にその方向に原子を動かすことによって力がゼロになる原子配列を求める計算を「構造最適化」と呼ぶ。一方原子に働く力からその原子の加速度をニュートンの運動法則から見積もり運動方程式を逐次積分することによって原子の位置の時間発展を追う計算を「分子動力学(Molecular dynamics)」と呼ぶ。 ASEでは構造最適化には * 最急降下法(MDMin) * 共役勾配法(Conjugate Gradient) * 擬ニュートン法(Quasi Newton) を、分子動力学法には * 速度ベルレ法(Velocity Verlet) * ランジュバン法(Langevin) を用いることができる。構造最適化および各手法の詳細は別稿もしくは専門書を参照せよ。 このような構造最適化計算を実行する際にはメソッド''StayAliveOn()''を指定する。これは''GetPotentialEnergy''メソッドが呼ばれるたびに、Dacapoプログラムが完全なNetCDFの出力ファイルを書き出して終了しないようにし、そのかわりにpythonプログラムが新しい原子位置を設定するのまで待機するように命令する。これにより、入出力に費やされる時間を短縮し、ネットワークファイル共有のトラフィックを軽減することができる。ただし、後述のNudged Elastic Band計算時にこの設定をしてはならない。メモリ上に数十個のDacapoプログラムがロードされることになってしまう。 ====== Wannier orbitals ====== See the [[Wannier example]]. ====== Nudged Elastic Band ====== See the [[Nudged Elastic Band example]].