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

MEM解析

 ここではMEMを用いた解析に関連する情報を掲示していきます。
--------------------------------------------------------------------------------
■ 原子による散乱
・ 原子の散乱能は原子番号(電子の数)に比例する。
・ 水素は散乱能が最も低い。
・ 重金属は良く散乱する。

■ 差Fourier合成法
・ モデルによる計算値との差分をフーリエ逆変換する。
  Δρ(x, y, z) = ρobs(x, y, z) - ρcalc(x, y, z)  
・球対称なモデルとの差分をとれば変形成分が現れて結合状態を議論できる。
※ 差フーリエは水素の位置や有無などを調べるのに用いられる

■ MEM解析の特徴
・ フーリエ逆変換の式を使わないため、打ち切り効果の影響が少ないという利点がある。
・ Rietveld解析は自由原子モデルを用いており、結合電子の情報は含まれていない。したがって、結合電子による寄与が現れやすい低角部では残差が大きくなる。
・ MEMは結合電子による寄与もすべて含んだ結晶構造因子(F)に合う電子密度分布を求めている。
・ フーリエ逆変換の式を用いずに、与えた結晶構造因子の値にその誤差の範囲内で一致する結晶構造因子を再現する電子密度分布の中で、最も情報エントロピーが大きなものを解として求めている。
・ それゆえ、与える結晶構造因子のデータの信頼性が非常に重要である。モデルに強く依存する補正はできるだけ避けた方が良い。
・ 整合のとれないたった一つのデータが電子密度分布に影響を与えることがあるので注意する。
・ 誤差の範囲内で一致する解なので、厳密に言えば正しい解ではない。しかし、フーリエ変換を使った方法よりは真の解に近い解が得られると考えられる。
・ 情報エントロピーが最大である解=最も平坦な電子密度分布を求めている。そのため、電子密度が高い・低いという構造をできるだけ出さないようにしており、原子核付近(電子密度分布の鋭いピーク部分)は再現が難しい。
・ 解析に用いた回折線の面間隔の最小値dminよりも細かいスケールの構造は再現できない。面間隔は見たい化学結合の結合距離の1/2のdが一つの目安と成る。
・ タンパク質の結晶構造解析ではよく「X.XÅ分解能で結晶構造を解いた」という言い方をする。
・ MEM電子密度解析のための回折データ測定においては、統計精度の高いデータが得られるように心がけることが重要。
・ MEM電子密度解析とRietveld解析を組み合わせることにより、大きな単位格子や複雑な結晶構造の解析、初期構造モデルの構築が可能である。
・ 通常の構造解析は散乱体を球として計算するために、観測値と計算値の差が大きくなる。MEMの場合は、ユニットセルを非常に細かく分割し、その一つ一つの小さなboxを散乱体としてみるので、細かなところまで表現することができ、観測値と計算値の差が小さくなる。
・ 共有結合電子の部分まで計算ができ、計算値と観測値の差が小さくなると、通常では見えにくい軽元素の電子密度が浮かび上がってくる。一方、差フーリエ合成だと、電子密度の計算において観測値と計算値の差が大きく、打ち切りによる波打の影響でぼやけが生じてしまう。

■ 初期電子密度分布の問題
・ MEM電子密度分布の解は初期電子密度分布の選び方に依存することが知られている。
・ Sakataらの解析法では、均一な電子密度分布を用いている。これは初期電子密度分布を選ぶ任意性が無くなり、誰が解析しても同じ一つの解に行き着く利点がある。
・ 一方、高圧における回折データの場合など、回折データの質が低く、データ数が少ない場合には、均一な電子密度分布からスタートするのは難しい。計算が収束しないこともある。
・ 自由電子モデルの結晶構造因子(原子散乱因子)から計算した電子密度分布を用いる方法もある。

■ 粉末X線結晶構造解析の流れ
1. 指数付け
・ JCPDFなどによるデータベースを用いて検索。
・ 類似物質の結晶データからの推測
・ 指数付けプログラム
2. 構造モデルの構築 ※
・ 類似物質の結晶構造からの推測
・ 未知構造解析プログラム
3. 構造の精密化
・ リートベルト法(格子定数や原子座標の精密化)
4. 電子密度解析 
・ フーリエ変換、マキシマムエントロピー法、多極子展開法などによる化学結合解析
※ 指数付けの後に格子定数の精密化としてLe Bail法やPewley法によるパターン分解を行うこともある。

■ 粉末回折法の注意点
・ MEMはモデルフリーの解析法であるので、モデルに強く依存する補正はできる限り避ける(消衰効果、吸収補正、選択配向など)。
□ 粉末回折法の利点
・ 単結晶を得ることができない、または得ることが困難な物質に対応可能
・ 測定が簡便。新奇物質の研究において素早く対応可能
・ データ全体の統計精度が高い
・ すべての回折線が同時に測定されるので、入射ビームの変動の影響を受けない
・ 測定中の稼働部が少なく、装置による統計誤差が混入しにくい(種々のアタッチメントの装着が比較的容易)
・ 消衰効果の補正が必要ない(吸収など補正が比較的容易)。
・ デバイリングから試料の質の評価が容易
□ 粉末回折法の欠点
・ 3次元の情報が1次元に重なり一部失われる
・ 単結晶に比べてパラメータの精度が悪い。非等方性熱振動の評価が難しい。
・ 電場や磁場に対する応答をみることができない(誘電歪みや磁歪など)。

■ MEM/Rietveld解析[2,3]
  結合電子、ディスオーダーの構造、動的な構造など、通常構造モデルに表現しにくい電子・原子の空間的分布の情報を得ることができる可能性がある。
1. Rietveld解析の結果を用いて回折強度の再分配を行う
2. 個々の反射の積分反射強度を求める
3. Rietveld解析により得られたスケール因子を用いて、積分反射強度から結晶構造因子を求める(位相はRietveld解析に用いた構造モデルの位相を与える)。

■ MEM/Rietveld解析の注意点
・ 予備的な解析のモデルが比較的良い近似モデルであることが重要。
(金属内包フラーレンでは炭素ゲージの形がほぼ分かっていること、多孔性配位高分子ではホスト骨格の構造がほぼ分かっていることが大事)
・ 通常、MEMのR因子は、リートベルトのR因子よりもほぼ同じか小さくなるが、イオン結晶のような閉殻を作るような物質ではR因子が悪くなることもある。
・ MPF解析が上手くいかないときは手作業でEの値を調整する。手作業は手間がかかり大変だが、最小二乗をPCにやらせるか人間にやらせるかの違いがあるだけである。
・ 手作業でEを決める場合には、cycle数が小さく、R因子が小さいものを採用する。さらに、e/Å3を大きくしたり小さくしたりして、ゴーストが生じていないかをチェックする。
・ ゴーストか軽元素か見分けが付かない場合は、マーデルングエネルギーを計算したり、理論計算と比較したり、中性子回折を用いたりして慎重に分析する。
※ cycle数が少ないと良い理由は、男女の関係と同じで、脈ありだとアプローチ回数が少ない、ということと同じだと考えれば良い(そう泉先生から教えて貰いました)。

■ 熱散漫散乱(Thermal Diffuse Scattering: TDS)の影響
・ブラッグ反射には裾野が広がったTDSが重なる
・TDSは格子振動に起因しており、その強度は(singθ/λ)^2に比例する
・高次の反射では、熱振動のためにブラッグ反射の強度が減衰するのに対し、TDSは強くなる。その結果、格子面間隔dが小さくなるにつれてTDS強度/ブラッグ反射強度は増加する
・粉末回折データの処理では、スムーズなバックグラウンドを差し引いてもTDSを取りきれないため、通常、TDSを無視する。そのため、リートベルト解析やMEM解析の結果にはTDSを無視した悪影響が大なり小なり残る
□ リートベルト・MPF解析におけるTDS
・高角(小d)領域では反射同士の重なり、回折強度の低下、TDSの増加が同時に起こり、結晶構造に関する情報が劣化する。そのため、2θ(TOF)の範囲を変えると、Bや密度分布が変わる
・できるだけ低温で強度データを測定する
・短波長を使用した放射光回折データやTOF粉末中性子回折データ(J-PARC)の解析では、格子面間隔dが小さい反射は解析に含めない(MEM解析では、非ゼロの構造因子を推定させる)

■ MEMの計算式[4-7]
□ Information entropy
  S = -Σr ρ'(r) ln {ρ'(r)/τ'(r)}
  ここでρ'(r) = ρ(r)/Σrρ(r), τ'(r)=τ(r)/Σrτ(r) である。電子密度ではなく確率電子密度を用ていることが特徴である。τ(r)は事前電子密度(prior density distribution)を表す。繰り返し計算(iteration)により解を求めるので、τ(r)はひとつ前のサイクルの電子密度となる。
  S = -ρ'(r) ln {ρ'(r)/τ'(r)} = ρ'(r) ln [{1/ρ'(r)}/{1/τ'(r)}] = ρ'(r) ln {1/ρ'(r)} - ρ'(r) ln {1/τ'(r)}
  ここで大事なことは、事前電子密度の情報の分だけ情報量が少なくなっていることがわかる。
□ Constraint function
  C = (1/N) Σk[|Fcalc(k) - Fobs(k)|2 / σF2(k)]
  ここで Fcalc(k) = V Σr ρ(r) exp (-2πik・r) 、σF2(k)はFobs(k)の測定誤差である。実験データFobsと電子密度分布から計算されるFcalcとが誤差の範囲内で一致するという束縛条件をここで付けている。統計学的にはC=1となることが期待される。
□ Information entropy with constraint
  Q(λ) = -Σr ρ'(r) ln {ρ'(r)/τ'(r)} + λ(C-1)/2
  最大値問題をラグランジュの未定乗数法により解く。λはラグランジュの未定乗数。結果はλの値に依存しないことが分かっているが、大きすぎると計算が発散する。
  ∂Q(λ)/∂ρ'(r) = 0 ← エントロピーが最大になる条件。ρ(r) ではなく、ρ'(r)で偏微分することにより、近似なしで解くことができる。
□ Final MEM equation
  ρ(r) = exp[ ln τ(r) + (λ/N) Σk {|Fcalc(k) - Fobs(k)|2F2(k)}exp(-2πik・r)]
□ 信頼度因子
  RFMEM = √[{Σk |Fcalc(k) - Fobs(k)|2Fobs2(k)}/{Σk Fobs2(k)/σFobs2(k)}]

■ 電子密度の変換 (V. G. Tsirelson, Acta. Cryst. B58 (2002) 632-639.)
in ∇ρ(r)=0 case,
Covalent type: ∇2ρ(r) < 0, he(r) < 0
Ionic type: ∇2ρ(r) > 0, he(r) > 0
he(r): electronic energy density

■ MEM電子密度解析のプログラム[8]
・ ENIGMA [1](並列処理MEM構造解析プログラム)
・ PRIMA (F. Izumi et al.) → Dysnomia (パスワードは speed710、泉先生に公開しても良いと言われたので記載。こんなに良いソフトが使えないなんて勿体無い)
■ 電子密度解析支援プログラム
・ GSAS (A. C. Larson & R. B. Von Dreele)
・ Fullprof (J. Rodriguez-Carvajal)
・ Rietica (B. Hunter)
・ JANA 2006 (M. Dusek et al.)
・ RIETAN-FP (F. Izumi)
■ 電子密度の描画
・ VESTA (F. Izumni et al.)
・ Limner (H. Tanaka et al.)
■ Powder Diffraction Software
  http://www.ccp14.ac.uk/

[1] H. Tanaka et al., J. Appl. Cryst. 35 (2002) 282.
[2] M. Tanaka et al., Nature 46 (1995) 377.
[3] M. Tanaka et al., Z. Kristallogr. 216 (2001) 71.
[4] http://www.netsu.org/j+/Jour_J/pdf/31/31-1-29.pdf
[5] http://www2.itc.nagoya-u.ac.jp/pub/pdf/pdf/vol07_01/003_010campus.pdf
[6] http://tdl.libra.titech.ac.jp/z3950/gakuipdf/1706450/170645009.pdf
[7] http://sbsp.jp/sbsp/Sb/sb72/028.pdf
[8] http://crystruct.info/Rietveld/index.html
[9] http://www.spring8.or.jp/ext/ja/iuss/htm/text/11file/glass_ceramic-2/yashima.pdf
[10] Y. Kawamura、johokiko, 2014.

※ 最強回折線のカウントと議論できる物性
・格子定数の決定:最強回折線の強度は500カウント以上
・リートベルト用のXRD:最強回折線の強度は5千から1万カウント程度以上
・MEM用のXRD:リートベルト用のXRDよりもさらに長時間測定したもの。
  放射光だと1回の測定5分程度で30万カウント程度(放射光は一回の測定で済む)。
  例として、放射光と同じ程度の精度で測定しようとすると実験室の装置では2から3日測定が必要
  私の経験だと、実験室の装置では4.5万カウントあたりから理論値と近い結果が得られるようになってくる(私の使っている実験室の装置だと、計数時間が10.0sで測定時間は18時間程度となり最強回折線の強度は4.5万カウント程度)。

※ MEM測定の注意点
・低角度の回折が全て必要(その部分が結合に重要なところになる。それが無いと原子の間が平らになる)
・高角度側は熱などの影響(熱散漫散乱?)を考えて必要であれば除く。私の場合は2θ=20-120で理論値と整合性のある結果が得られた。
・信号/雑音比=S/N比が悪いとぼやけた像しか得られない。S/N比を高くして、ようやく結合の状態が議論できる。
・選択配向が無いこと。それがあると解析できない。
・出来るだけ単体(リートベルト解析したときに成分を指定して、MEMに持っていくことはできますが……)
・議論をしやすくするために可能な限り低温で(といっても、実験室で低い温度にすることは難しい。液体窒素の温度にできれば、放射光と温度に関して同等になる。放射光の場合はHeガス(気体)(液体窒素ではなくHeガスでした)を吹き付けているだけですけど。それで77 K程度まで冷えます)

※ リートベルト解析
  私はリートベルトを相の同定と、置換元素のサイト選択性、占有率を調べるために用いている。
  置換したときのXRDからの(Vegard則からの)原子間距離とXAFSでの局所的な原子間距離が一致しないことがあるのは当然なので、そこを指摘するレフェリーがいたら丁寧に答えてあげよう。H. Sato et al., Jpn. J. Appl. Phys. 31 (1992) 1118.を読んでください。
--------------------------------------------------------------------------------
■ VESTA (unit)
Unit: Please, show "[17.4.2 Volumetric data] and [17.4.3 Structural and volumetric data]" in VESTA_manual.pdf
・ Dysnomia (*.grid): Å-3
・ WIEN2k: bohr-3
・ VASP (CHG, CHGCAR): bohr-3
・ Abinit (*.xsf): bohr-3
・ PWscf (*.xsf): The unit of primitive lattice vectors (PRIMVEC) and Cartesian coordinates is Å in the XSF format.
--------------------------------------------------------------------------------
■ VESTAでの分割数の算出の仕方
1. 粉末法の場合、限界球の直径は 4/λ である。※1
2. Utilities -> Model Electron Densities -> Resoultion (Å) -> 0.042 -> OK
3. 下の部分に Dimentions として記述される。
※ 原子変位パラメーターが小さい場合や核密度の計算時には、分割数を細かくする。ただし、電子密度の計算では、原子散乱因子のパラメーターの信頼領域が sinθ/λ= 6 Å-1 (resolution指定で 0.042 Å相当)しかない。
[1] 中井泉、泉富士夫 編集、「粉末X線解析の実際 第2版」、朝倉書店
※1 例えば、フーリエ変換の場合、短波長のMoターゲットを使って測定した回折データでは 0.2 Å程度の間隔で電子密度を計算すれば十分で、それ以上細かく刻んでも、見栄えは良くなるが新しい情報は増えない。ただし、この話は電子密度図の分解能についてであり、原子間距離がどの精度で求まるかという話とは異なる。
[2] 門馬綱一、「電子・各密度分布の決定」、日本結晶学会講習会 「粉末X線解析の実際」、2013年7月10日
--------------------------------------------------------------------------------
■ Dysnomia [Dys1-3]
□ 解析手順(秀丸を用いた場合)
1. RITAN-FPでリートベルト解析を行う。
  空間群や構造が正しいかなどの確認のための荒い計算をする場合には6割以上解析が終わっていればよいが、RIETAN-FPでリートベルト解析を可能な限り収束させておく方が良い。
2. 原子の占有率や原子変位パラメーターは固定
3. XRDでおかしな部分以外、全ての低角のピークを int に入れる
4. リートベルト解析を実行する。値が収束するまでRIETANを押す
5. ins ファイルで、
   NMEM = 1 を確認し、NPIXAからCを56(VESTAが推奨値を出すのでそれを入力する)へ。
   NPIXA = 56: Number of pixels along the a axis.
   NPIXB = 56: Number of pixels along the b axis.
   NPIXC = 56: Number of pixels along the c axis.
6. lst ファイルから 密度の下に書かれている 全電子数をコピー
7. ins ファイルで、
   TSCAT1 = 471.276: Total number of electrons (X-ray diffraction) or
 とあるTSCAT1に全電子数をペースト
8. RIETANを押す
9. prf を用意する。サンプルファイルにあるものを利用したりすると良い。
10. Dysnomia を押す
11. VESTAで pgrid ファイルを読み込む
  Properties... > Isosurfaces > Render from front to back にチェックを入れる > No. 1の欄を押して、 isosurface level: を変える。
  筆者の例だと、理論計算での出力ファイル case.rho で 0.041 / bohr3  (= 0.276 / Å3 )、Dysnomiaの出力ファイル case.pgrid で 0.45 /Å3 にしている。
※ 他の系での報告でも MEM : 0.3 /Å3、DFT : 0.2 / Å3 で、MEMとDFTで 0.1 〜 0.2 で違いが出ているようなので少し安心している。可能なら、情報やデータをもっと集めたいと思う。

□ case.3ed
 File > Export Data... > ファイルの種類から*.3ed を選択する。
 電子数は当該原子(化学種)に対応するピクセルの電子密度の総和: http://www16.ocn.ne.jp/~structur/qa/qa.html

□ Line Profile
 Utilities > Line Profile... > internal coordinates (格子定数で割って1を最大とした座標で位置を記述)

□ 2D Data Display
1. Utilities > 2D Data Display...
2. Slice > input hkl ( Distance from origin: の値に注意)
 Edit > Lattice Planes... > New > input hkl, and ??? x d (Slice と同じなので、これでチェックする)
3. General > input Max.: and Min.: ( e/Å3)
4. Contours > Draw contour lines ( e/Å3)

※ Dysnomiaでの等価面レベル( Isosurface level)の単位はÅ-3 である。値はÅ3 中の電子の数となる。
※ VASPのCHGCARは[e/bohr3]単位であるので注意(WIEN2kも同じisosurface levelで同様の分布になる)。
※ 1 [e/bohr3] = 6.748333 [e/Å3]
※ ドープした系での等価面レベル = 無ドープの等価面レベル * (無ドープの全電子数 / ドープした系での全電子数) *(無ドープの体積 / ドープした系での体積) (← この式が正しいか確認が必要) であるが、ドープした系ではほとんど違いがなかった(等高線を描くときのラインの線幅程度しか変化がなかった)。
[Dys1] http://fujioizumi.verse.jp/visualization/visualization.html
[Dys2] http://nano.kut.ac.kr/nanolab/images/b/bf/Kim_DH_2009_MRS_Pd_substitution_at_orthorhombic-NiSi_Si_%28010%29_epitaxial_interface_studied_by_density_functional_theory.pdf
[Dys3] http://renqinzhang.weebly.com/uploads/9/6/1/9/9619514/charge_density_difference.pdf

□ パラメーターの設定
A. case.out を見て、高次のCn(n>2)が1より極端に大きい場合、多重F制約条件(MEMの推定する誤差 (|Fo|-|F|)/σ を改良する)であるλn (n>2)の値を大きくする。つまり、Fractions of lambda for generalized constraints での ID を 1 0 0 0 0 0 0 0 から 1 1 0 0 0 0 0 0 0 0 としたりする。
 残差のn次中心積率が表示されている部分で、右側の列はその自然対数となる。Cn≒1, lin(Cn)≒0となることが理想である。
B. case.raw を見て、ガウス分布に近い形になっているかをチェックする。中心から遠くでも値が大きく存在している場合は、多重F制約条件であるλn (n>2)の値を大きくする。
C. case.raw を見て、Number of cycles や、MEMの信頼度に関する出力である CONSTR,やCONSTRw、RFやwRFをチェックする。これらが小さい値を採用する。※1
D. MEM解析における信頼度因子RFがリートベルト解析における信頼度因子RFよりも小さくなることが望ましい。
E. MEM解析がただ収束すればよいというものではない。電子・原子核密度の等値曲面を精査し、必要なら、化学結合論的・結晶化学的に自然な密度分布となるようEを微調整しなければならない。Gaussian, WIEN2k, VASP, ABINITなどの電子状態計算プログラムで決定した電子密度分布とよく似た分布となるようにEを調節するという便法も、ときには有効であろう。[P1]
F. 結合がないところに電子密度があるものは解析をやり直す。ゴーストやリップルに注意する。
G. 高角領域では、反射同士の重なり、回折強度の低下、TDSの増加が同時におこるので、結晶に関する情報が劣化する。そのため、これらの影響が強くない角度までを解析するようにする。
[P1] 中井ら著、粉末X線回折の実際 第二版、朝倉書店
※ 金属でのMEMの論文や報告が少ないのは、X線では結合状態が見えるが、金属で研究している人が少ないこと、中性子では核を見ているのでディスオーダーがなければあまり面白くないという理由がある。金属のMEMでも信頼性は十分にある(例えば、東大院新領域の木村研の論文など)。
※ 固溶体を分析する場合、1%の置換量で占有率や結合状態をMEMで議論しては駄目だとは言えない。慎重に実験し、色々な(他の手法での)データと付き合せて議論する必要がある。
※ ただし、Eを大きくするほどR因子は小さくなるが、ゴースト(物理的におかしな部分に電荷密度があると表示される)が出やすくなる。理論計算と比較するなりして、物理的に正しいもののなかでNumber of cycles 及び 信頼度に関する出力の値が小さいものを探す。Eの最適値は入射ビーム強度、検出器の種類、試料、測定時間にも依存する。Eは通常、正常な範囲のσ|Fo|を与え、かつMEMの反復回数が最小になるように設定する。[P1] 以下の手順も試してみるとよい。
※ 空孔サイトに何が存在しているかを明らかにするのは難しい。特に低角度側のデータが信頼できないものである場合には、空孔サイトに何の原子が存在しているかをMEMで議論するのは避けた方がよい。なぜなら、低角度側のデータを除くと、空孔サイトの電子密度が大きくなったりするからである。この件に関しては私の力不足で学生に多大な苦労をさせてしまうことになった。反省したい。

◇ E値の決定法(Y. Kawamura、johokiko, 2014.)
1) 中性子の場合はσFの値が1桁、放射光の場合はσFが0.1程度となるEの値を探す。
 適当な反射(例えば、最強回折ピークなど)でF0/I0を計算する。この値が1桁になる1/Eを計算する。または、E=1.0としてcycle数を1にしてMEM計算する。実効したときに作成される*.outファイルを開き、その中に記述されているσFを見る。E=100にするとσFは1/10に、E=0.01とするとσFは10倍になる。Eを変えてσFが目的の値になるEを探す。
2) cycle数を10程度に変更して、Eの値を変えてみる。最初は90%程度の値から下がらなかったR値が一気に下がるEの値を見つける。
3) cycle数を大幅に増やしてEの値を変えながら、出力された結果でcycle数とR値が小さくなる最適値を決定する。
※ R値が最小にならず、cycle数の増加でRが小さくなるような場合は、R値とcycle数をプロットして両者が一番小さくなる部分を探す。その後、e/Å3やfm/Å3を変えてゴーストが出にくいE値を採用する。
※ 構造解析でのRと同じになるEの値を見つけるのも一つの手段。

□ 入力ファイル
  下記は、私が研究に用いているHH(No.216)またはFH(No.225)の試料に対して、WIEN2kでの電荷密度分布の”変化の傾向”(しかし絶対値は)とほぼ一致するパラメータを探した場合での結果である。
・ SPring-8 BL02B2 イメージングプレート利用での場合(下記の二つが候補になった。他ではゴーストが出る)
 下記の2つのどれかを選ぶ。恐ろしいほどWIEN2kの結果と合うので驚く。また、電荷分布の微妙な変化も見つけられる(ただし、このパラメータ設定で良いかは分からない。現在、精力的に調査中である)
(※ case.lst において、E、つまり、E(SCIO)が出力される。それを採用した場合、私の用いている試料や測定環境におけるデータでは、MEM計算においけるCycle数が多くなり、加えて、WIEN2kの結果からも大きく離れる結果となった)

※ ちなみに、VESTAでは原子散乱因子を用いて結晶構造からモデル電子密度・核密度を得ることができる。分解能は実空間で0.042Å、逆格子空間で6Å-1までで、それ以上に意味は無い。
※ 本物の水素は丸いし距離も良いところにある。逆にゴーストだと丸くないし距離も変なところに電荷がある。

◆ 候補1
---------------------------------------------------------------
# The initial Lagrangian multiplier, lambda (dummy if LM_type = 1 or 2), and
# A parameter x to impose weighting factors on the basis of
# lattice-plane spacing.
0.001 2
# The coefficient, E, to adjust estimated standard uncertainty of
# structure factors.
0.1 から 1 の値(0.8を推奨か?)
---------------------------------------------------------------

◆ 候補2
---------------------------------------------------------------
# The initial Lagrangian multiplier, lambda (dummy if LM_type = 1 or 2), and
# A parameter x to impose weighting factors on the basis of
# lattice-plane spacing.
0.001 0
# The coefficient, E, to adjust estimated standard uncertainty of
# structure factors.
0.05
---------------------------------------------------------------

◆ 候補2でのprfファイルを検証と参考のために示す。下記はFeやNiなどの遷移金属を含むFm-3mまたはF-43mの物質に対して検証したものである(SPring-8 BL02B2 イメージングプレート利用でのデータを用いている)。
---------------------------------------------------------------
$PREFERENCES

# The name of the MEM data set file, *.fos or *.mem.
case.fos
case.out
case.pgrid
case.fba
case_eps.raw

# The type of the MEM data set file.
# 0. (a) X-ray diffraction data or (b) neutron diffraction data of a compound
#    containing no element with a negative bc value.
# 1. Neutron diffraction data of a compound containing at least one element
#    with a negative bc value.
0

# From which densities will you start?
# 0. Uniform densities.
# 1. Densities recorded in the 3D densities file.
0

# The approximation mode.
# 0. 0th order single-pixel approximation (ZSPA).
# 1. The Cambridge algorithm.
1

# Will you use the G-constraint?.
# (This option is a dummy to keep backward compatibility with PRIMA.)
1

# Which initial Lagrangian multiplier will you use (LM_type)?
# 0. The value written in this file.
# 1. The value calculated by Dysnomia.
# 2. The value written in the MEM data set file.
1

# The initial Lagrangian multiplier, lambda (dummy if LM_type = 1 or 2), and
# A parameter x to impose weighting factors on the basis of
# lattice-plane spacing.
0.001 2

# The coefficient, t, to adjust the Lagrangian multiplier.
0.05

# The coefficient, E, to adjust estimated standard uncertainty of
# structure factors.
0.8

# Will you save a feedback data file?
# 0. Yes (output only individual F data).
# 1. Yes (output all the F data including those for grouped reflections
#    and estimated for unobserved reflections).
# 2. No.
1

# Will you save a e(GAUSS) distribution data file?
# 0. Yes (raw data file, * eps.raw).
# 1. No.
0

# Maximum number of cycles.
100000

# Save output files (dummy)
1

# Fractions of lambda for generalized constraints.
# (8 parameters; l_2, l_4, l_6, l_8, l_10, l_12, l_14, l_16)
1 0 0 0 0 0 0 0

# Which prior densities will you use?
# 0. Uniform densities.
# 1. Densities recorded in the 3D densities file (*_prior.pgrid).
0
---------------------------------------------------------------
■ Dysnomia の入出力ファイル

□ 入力ファイル
・ *.fos: RIETAN-FPから作られる”観測構造因子”。
 (*.itxで緑の線(ycal) でピークが描かれている反射に対する構造因子が作られている)
・ *.prf: MEMでの計算条件のファイル

□ 出力ファイル
・ *.out: 計算結果の概要
・ *.pgrid: 電子・核密度データ
・ *.fba: 結晶構造因子のリスト
・ *.eps.raw: 観測構造因子と計算構造因子の残差

□ ケンブリッジアルゴリズムが計算可能
  理論的に精度が高い方法で、厳密解が得られる方法であると言われている。
 エントロピーが最大になる解がえられる。

□ study (学習方法) J-Stage
[1] RIETAN徹底活用ガイド(3) : https://www.jstage.jst.go.jp/article/jcrsj1959/44/6/44_6_380/_pdf
[2] RIETAN徹底活用ガイド(2) : https://www.jstage.jst.go.jp/article/jcrsj1959/44/5/44_5_311/_pdf
[3] RIETAN徹底活用ガイド(1) : https://www.jstage.jst.go.jp/article/jcrsj1959/44/4/44_4_246/_pdf
---------------------------------------------------------------
■ Charge Density Plots [CDDP1-2]
◇ VESTA & WIEN2k or VASP [CDDP1-2]
1) case.cho, CHGCAR, or CHG
2) Edit > Edite Data > Volumetric Data...
3) Isosurfaces > Import... > case.cho, CHGCAR, or CHG > Replace current data, Bohr^(-3) to Angstrom^(-3) > OK > OK
4) Edit > Edite Data > Volumetric Data... > Structure parameters > Import... > case.cif > OK > OK
5) objects > Proterties > Isosurfaces... > Render from front to back > No. > Isosurface level: [Å^-3]

■ Charge Density Difference Plots [CDDP1-2]
◇ VESTA & Dysnomia (MEM)
1) case.pgrid
2) Edit > Edite Data > Volumetric Data...
3) Isosurfaces > Import... > other case.pgrid > Subtract from current data > OK > OK
4) Edit > Edite Data > Volumetric Data... > Structure parameters > Import... > case.cif > OK > OK
5) objects > Proterties > Isosurfaces... > Render from front to back > No. > Isosurface level: [Å^-3]

[CDDP1] http://vaspnotes.blogspot.jp/2013/10/charge-density-difference-plots.html
[CDDP2] http://renqinzhang.weebly.com/uploads/9/6/1/9/9619514/charge_density_difference.pdf
---------------------------------------------------------------
■ Charge Density Plots [CDDP1-3]
◇ VESTA or XcrySDen & Abinit [CDDP1-3]
1. mpirun -np 4 $HOME/abinit-6.8.2/bin/abinit < case.in
2. $HOME/abinit-6.8.2/bin/cut3d
3. What is the name of the 3D function (density, potential or wavef) file ?
  caseo.DEN
4.   Does this file contain formatted 3D ASCII data (=0)  
   or unformatted binary header + 3D data        (=1) ?
   or ETSF binary                                (=2) ?
  1
5.   What is your choice ? Type:
  9 => output .xsf file for XCrysDen
  9
6.   Enter the name of an output file:
  case_cut3d.xsf
7.  Do you want to shift the grid along the x,y or z axis (y/n)?
  n
8.   More analysis of the 3D file ? (1=default=yes,0=no)
  0
9.  case_cut3d.xsf
10. Edit > Edite Data > Volumetric Data...
11. Isosurfaces > Import... > case_cut3d.xsf, CHGCAR, or CHG > Replace current data, Bohr^(-3) to Angstrom^(-3) > OK > OK
12. Edit > Edite Data > Volumetric Data... > Structure parameters > Import... > case.cif > OK > OK
13. objects > Proterties > Isosurfaces... > Render from front to back > No. > Isosurface level: [Å^-3]
[CDDP1-3] http://www.abinit.org/documentation/helpfiles/for-v5.8/tutorial/lesson_analysis_tools.html
---------------------------------------------------------------
■ Charge Density Plots
◇ VESTA or XcrySDen & PWscf
---------------------------------------------------------------
■ Charge Density Plots [CDDP1-3]
◇ VESTA or WIEN2k [CDDP1-3]
1. run_SCF
2. cd case
3. use wien2venus.py: http://www.nims.go.jp/cmsc/staff/arai/wien/venus-j.html  (Japanease), http://www.nims.go.jp/cmsc/staff/arai/wien/venus.html  (English)
4. wien2venus.py 56 56 56  or wien2venus.py -c 56 56 56
5. mv case.cho3d case.cho
6. case.cho > VESTA
7. Edit > Edite Data > Volumetric Data...
8. Isosurfaces > Import... > case_cut3d.xsf, CHGCAR, or CHG > Replace current data, Bohr^(-3) to Angstrom^(-3) > OK > OK
9. Edit > Edite Data > Volumetric Data... > Structure parameters > Import... > case.cif > OK > OK
10. objects > Proterties > Isosurfaces... > Render from front to back > No. > Isosurface level: [Å^-3]
http://www.abinit.org/documentation/helpfiles/for-v5.8/tutorial/lesson_analysis_tools.html
---------------------------------------------------------------
■ Fourier Synthesis (Dysnomia + VESTA ver.3)
VESTA & RIETAN-FP
1) case.ins
2) Utilities > Fourier Synthesis... > Import... > case.fos > Fo-Fc > Calculate > Close
3) objects > Proterties > Isosurfaces... > Render from front to back > No. > Isosurface level:
---------------------------------------------------------------
■ Model Electron Densities
VESTA & RIETAN-FP
1) case.ins
2) Utilities > Model Electron Densities >0.042 > case.fos > Fo-Fc > Calculate > Close
3) objects > Proterties > Isosurfaces... > Render from front to back > No. > Isosurface level:
---------------------------------------------------------------
--------------------------------------------------------------------------------
■ MPF解析
1) RIETAN-FPでフィットさせる
2) case.fos を用いて Dysnomiaで解析(上記を参照)
3) case.ins で NMODE = 2 として、case.fba を用いて RIETAN-FPで計算
4) case.prf で下記のようにする。
# The coefficient, SCIO, to adjust estimated standard deviations.
_SCIO_
#E=200 300 400 などと、計算したいEの候補を複数記述
5) RIETAN_VENUS/Commands/MPF_multi.command を計算するフォルダーに入れる
 (bash.exeのファイルに関連づけられていなければ、RIETAN_VENUS/Commands/bash.exeにあるbash.exeを関連付ける)
6) MPF_multi.command をダブルクリックして実行させる
何か聞かれたら RIETAN_VENUS/Commands/bash.exe を指定する
7) case.log を見て、正しく収束していてRwpが小さいものを探す
--------------------------------------------------------------------------------
■ VESTA + Dysnomia
 ホーエンベルグ・コーンの定理により、電子密度ρ(r)とポテンシャル(エネルギー)が一対一対応する。そのため、実験的に電子密度分布が求まれば多くの情報を"理論的には"得られる。
※2014年4月時点で、VESTAに備えられているポテンシャル及びエネルギーの計算はρ(r)から計算する方法を用いている。有機分子に対しての報告ではあるが、DFTでの結果と比較的良い結果を示している(MEMはSrTiO3やBN、MgB2でも行われている)。[3-5] 一方、第一原理計算では、平面波にしてから運動エネルギーの計算を行っている(より厳密に運動エネルギーを求められるため)。

◇ Dysnomia -> VESTA で得られる情報 [1-3]
1. 電子密度 ρ(r)のラプラシアン (Laplacian of the electron density *.led): ▽2ρ(r)
2. 電子運動エネルギー密度 (Electronic kinetic-energy density *.ked) : g(r)
 λと k の2つの定数が存在する。これらの値は理論的に求められており、論文[3-5]やVESTAでは λ=1/72 と k = 1/6。
3. 電子ポテンシャルエネルギー密度 (Electronic potential-energy density *.ped) : v(r)
 正確な核の電荷Zが不明であるため、v(r) = (ℏ2 / 4m) * ▽2ρ(r) - 2*g(r) として求める。
4. 電子エネルギー密度 (Electronic energy density *.ted): he(r)
 厳密な関係式  2* g(r) + v(r) = (ℏ2 / 4m) * ▽2ρ(r) を用いて、he(r) = g(r) + v(r) から求める。
5. Natural logarithm (*.lned)
※ 「▽2ρ(r) 以外は Input data の単位に関わらず、bohr-3 で出力されるので注意」(確認中)VESTA manual [2] の 「Table 14.1: Units of converted volumetric data.」を参照。
※ 1 bohr = 0.529177 Å。1 [e/bohr3] = 6.748333 [e/Å3]。エネルギーの単位は hartree = 2 Ry = 27.211 eVである。

◇ 変換方法
1. case.pgrid を VESTAに入れる
2. Utilities -> Conversion of Electron Densities...
3. Files to be output で、欲しいデータにチェックを入れる。

◇ 結合様式の評価 (he(rb) よりも▽2ρ(rb) の使用を を推奨。rbは▽ρ(r) = 0 となる位置)
・ shared-type atomic interactions : ▽2ρ(rb) < 0 (recommend) or he(rb) < 0
 電子が両側の電子によって共有されている。共有結合性の評価
・ closed-shell interactions : ▽2ρ(rb) > 0 (recommend) or he(rb) > 0
 閉殻タイプの原子間相互作用。イオン結合性の評価
※ ファインマンによると、電子密度分布からの結合様式は、原則同じ(+ と - で働く力で分かる)で程度が違うだけという。電子密度分布が分かれば、古典的に見て説明することもできる。(文献調査中)
※上記のように定量的に比較するか、または電子密度分布から視覚的に判断する。電子密度分布から視覚的に判断する場合、「イオン結合か? 金属結合か? 共有結合か? 分子間か?」といったことを比較するのは、極端な場合を除いて、相対的な比較しかできない。

◇ Atoms in molecules (VESTAに実装されているかは不明)
 ▽ρ(r) = 0 となる点で空間を区切って、その空間内で電子密度分布を総和する。

[1] RIETAN徹底活用ガイド(3) : https://www.jstage.jst.go.jp/article/jcrsj1959/44/6/44_6_380/_pdf
[2] K. Monma and F. Izumim, "VESTA: a Three-Dimensional Visualization System
for Electronic and Structural Analysis", Oct. 22, 2013.: (14.13 Conversion of Electron Densities) in VESTAのマニュアル
[3] V. G. Tsirelson, Acta Crystallogr., Sect. B: Struct. Sci., 58 (2002) 632.: http://quant.distant.ru/files/pdf/papers/0004.PDF, http://books.google.co.jp/books?hl=ja&lr=&id=ud_DmLcBXZ4C&oi=fnd&pg=PR7&dq=Vladimir+G.+Tsirelson&ots=BcwuCtZWH1&sig=P52iWWUWVZYYOr38-dD4xPQIo-g#v=onepage&q&f=false
[4] E. A. Zhurova et al., J. AM. CHEM. SOC. 128 (2006) 14728.: http://pubs.acs.org/doi/pdf/10.1021/ja0658620
 CRSTAL98, B3LYP, 6-311G(d,p) での計算結果と比較。 対象は Pentaerythritol tetranitrate (PETN) = C5H8O12N4
[5] V. G. Tsirelson et al.,  Computational and Theoretical Chemistry 1006 (2013) 92.
[6] R. G. パールら、「原子・分子の密度汎関数法」、シュプリンガー・フェアラーク東京: 6章を参照のこと。
[7] W. Yang, Phys. Rev. A 34 (1986) 4575.
[8] W. Yang et al., Phys. Rev. A 34 (1986) 4586.
※ 素人的には、粒子密度をフーリエ変換して平面波にし、そこから運動エネルギーを計算する。ラプラシアン(Laplacian)から運動エネルギーを引けば、電子のエネルギーが求められる気がするが、そのような論文を簡単には見つけられなかった。

◇付録 g(r) の式の詳細 [6]
第1項: Thomas-Fermi の運動エネルギー(一様な電子ガスのモデル)である。
第2項: 密度勾配での補正である。
 Weizsacker の考えでは、パラメータ λは 1 であったが、その後の密度勾配展開によるやり方でλ=1/9であることが示された。
 Wizsacker補正はVESTAマニュアルに記載のあるg(r)の第2項に1/8とλが掛けられる式となっている。
 VESTAではWizsacker補正で表れている 1/8 に 1/9 を掛けて λ = 1/72 としたものと一致している。
 この密度勾配の補正(λ=1/9も含めて)は、1粒子のグリーン関数のWigner変換を半古典的にℏ展開することで得られる。その場合、高次の密度勾配補正が得られるが、原子や分子の場合には、これが正しいのは4次までである(原子の場合、6次の密度勾配補正は発散してしまうため)。
第3項: 1電子(または1粒子)のグリーン関数に対して平均経路の近似を行うとg(r)の結果が得られる。論文(Yang, 1986)[7]では1/12として記載されている(平均経路を使わない場合には1/6。1/6よりも1/12の方がHFの結果に近い事が論文[8]に記載されている)。VESTAでは1/2が既に掛けられており、k=1/6 とすることで論文[7]と一致する。

希ガスでの運動エネルギーの比較[8]
希ガス:HF値:g(r) = TF-1/9W-1/12LH の値 ( { g(r) - HF値 } / HF値 * 100 [%])
He:2.862:2.931(2.42%)
Ne:128.546:129.08 (0.41%)
Ar:526.812:528.49 (0.32%)
Kr:2752.03:2749.39 (0.10%)
Xe:7232.14:7222.61 (0.13%)

分子での運動エネルギーの比較[8]
分子:HF値:g(r) = TF-1/9W-1/12LH の値 ( { g(r) - HF値 } / HF値 * 100 [%])
H2:1.128:1.138(0.88%)
F2:198.5:199.23 (0.36%)
CO:112.7:113.51 (0.72%)
HF:100.0:100.56 (0.56%)
BF:124.2:125.07 (0.70%)
CO2:187.5:188.74 (0.66%)
C2H2:76.70:77.55 (1.11%)
※ 上記の結果から考えると、有機分子に関しては、HFと比較して運動エネルギーの誤差が3%内に収まる(分子なら約1%以内)。MEMで指定しなければならない誤差が真の値からどの程度異なっているかという問題を除けば、論文ではポジティブな結果を報告するものであるということを考えても、MEMの結果からかなりの精度で運動エネルギーや電子のエネルギーを求めることができると思われる。仕方ないのだが、HFよりはPBEで比較して欲しかった。
※ 金属に関する報告は未だ見つけていない。
--------------------------------------------------------------------------------
■ Maximum Entropy Patterson Method (MEPM): 最大エントロピーパターソン法 (ALBA)
1. case.ins (Le Bail) (change from Rietveld to Le Bail analysis) [2]
 NMODE = 4
 NPRINT = 2
 # O1/O- 1.000 ・・・・・・
 # O2/O- 1.000 ・・・・・・
 ・・・・・・
 NINT = 1
 NRANGE = 0
 LBG = 1
 NAUTO = 2
 NSFF = 0
 NOPT = 1 
 NCYCL = 50
 LQUART = 0
 NUPDT = 0
 NDA = 0
 EPSD = 0.001
2. RIETAN or case.bat file
3. 入力ファイル case.alb を作成
----------
# Title. (CHARACTER*80)
SPring-8 X-ray

# Unit cell parameters.
   6.09631  6.09631  6.09631  89.999  89.9984  90.0107

# Space-group number, setting and number of pixels.
  216    1   56   56  56

# Integrated intensities were determined from
# 0: X-ray diffraction data.
# 1: Neutron diffraction data of a compound containing no element
#    with a negative coherent-scattering length.
0

# The number of chemical species.
    3

# Real chemical species and numbers of atoms in the unit cell.
  Zr 4
 Ni 4
  Sn 4

# The name of an intensity data file including its absolute path.
# (*.ffo created by RIETAN-2000, *.hkl created by EXPO, or *.raw)
./case.ffo

# Wavelength (dummy for *.ffo files)
0.45993

# Which integrated intensities will you use?
# 0: The integrated intensity is defined as |F|^2.
# 1: The integrated intensity is defined as multiplicity X |F|^2.
0

# The initial Lagrangian multiplier.
0.005

# The coefficient, t, to adjust the Lagrangian multiplier.
0.01

# The coefficient, SCIO, to adjust estimated standard deviations.
# (dummy for *.hkl created by EXPO)
  1000

# Will you group overlapped reflections? (NGSW)
# 0: Yes (reflections overlapped heavily are grouped).
# 1: No (all the reflections are independently input).
0

# Maximum difference in d/Angstrom in grouped reflections.
# (dummy if NGSW = 1)
  0.0003

# Parameters (v, p, and b) to estimate standard deviations, s.
# s = v * sqrt[p*|F^2| + (b + |F^2|)/|F^2|].
# (dummy for *.ffo created by RIETAN-2000)
  1.000  0.0005  0.000

# The maximum number of MEP iterations.
10000
----------
4. ALBA
5. VESTA で case.pri を開く。
※ 私のやり方が悪いと思うが、ハーフホイスラーでは本来空孔サイトであるはずなのに電荷密度が充填サイトと同じように表示される。
http://fujioizumi.verse.jp/visualization/VENUS.html#ALBA
http://science.shinshu-u.ac.jp/~muki/VENUS.html: ALBAのLinux版は配られている。羨ましい・・・・・・。
--------------------------------------------------------------------------------
■ Charge Flipping Method (CFM): チャージフリッピング法
□ RIETAN & VESTA & SUPERFLIP [1]
・ superflip download & set up
1. Superflip: http://superflip.fzu.cz/
2. zipped exectutable for Windows: download at the first Table
3. rename superflip_win to SUPERFLIP
3. copy SUPERFLIP to C:\Program Files\SUPERFLIP
・ EDMA download
1. Superflip: http://superflip.fzu.cz/ 
2. zipped exectutable for Windows: download at the second Table
・ download python code and text file from [2]
1. fp2sflip.py
  change fp2sflip.py code following
  file3 = open('/Users/masami/sflip-symmetry.txt','r') > file3 = open('./sflip-symmetry.txt','r')
2. sflip-symmetry.txt
・ use
1. case.ins (Le Bail) (change from Rietveld to Le Bail analysis) [2]
 NMODE = 4
 NPRINT = 2
 # O1/O- 1.000 ・・・・・・
 # O2/O- 1.000 ・・・・・・
 ・・・・・・
 NINT = 1
 NRANGE = 0
 LBG = 1
 NAUTO = 2
 NSFF = 0
 NOPT = 1
 NCYCL = 50
 LQUART = 0
 NUPDT = 0
 NDA = 0
 EPSD = 0.001
2. RIETAN or case.bat file
3. put fp2sflip.py  and sflip-symmetry.txt in folder including case.ins
  (input files are case.ffo and case.lst )
4. fp2sflip.py case.ffo (win) or ./fp2sflip.py case.ffo (linux)
   ( option [2])
  #composition Mg8 Si4 O16 > composition your using atoms
  #histogram composition > histogram composition
5. check wheter or not the lattice constants (a, b, c) and angle ( alpha, beta, gamma) in case.inflip are matching space group. 
  e.g. 89.9900 > 90.0000
6. SUPERFLIP or superflip.exe case.inflip
EDMA
7. composition ### ### > composition your using atoms
8. EDMA case.edma
9. case.cif to VESTA
※ 私が悪いのだろうが、superflip とEDMA でハーフホイスラーでは空孔サイトも充填サイトもほぼ同じ電荷密度で出力される。
※ EDMA ではドーピングした系に対して、原子の数を小数点を入れて指定したが、占有率が1での結果が出力された。
[1] http://superflip.fzu.cz/
[2] Superflip_and_RIETAN-FP: http://www.misasa.okayama-u.ac.jp/~masami/pukiwiki/index.php?Superflip_and_RIETAN-FP
-----
□ Fox & VESTA & SUPERFLIP [3, 4]
・ Fox download
1. download: http://fox.vincefn.net/Download 
2. 解凍して .exe と記述のあるものを実行
・ fox2sflip.py
1. download: http://www.misasa.okayama-u.ac.jp/~masami/pukiwiki/index.php?Superflip_plus_Fox
・ use Fox (Le Bail analysis) (矢印をボックスなどに近づけると英文で説明が出るのでそれを参考にすると良い)
1. Object > New Crystal
 (File > open > cif を選択しても良い)
2. SpaceGroup: に空間群の記号または数値を入力する。そうすると直ぐ上の行(格子定数と角度を入力する行)が変化する。
3. Object > New PowderPattern
4. Data > Import 2Theta-Obs Pattern (RIETANでのcase.int の最初の二行を除いたものを作成して読み込ませる)
5. Pattern > Show Graph (回折パターンが表示されれば読み込みは上手くいっている) (フィッティングした結果が赤線で示される)
6. Phases > Add Crystalline Phase
7. Phases > Add Background > 50
8. Radiation > Monochromatic Wavelength
9. Max Sin (theta)/lambda: 0.400 (画面の中央のChi^2やGoF, Rwp, Rpのある行の上にある。値を大きくすると指数づけされる範囲が広がる。私の場合は高い角度まで考慮したいので 500000 にした。 )
10. Profile > Profile Fitting  + Le Bail extraction > Quick Fit (何度も押して前の画面の中央部分に記載されている Rwp が下がらなくなるまで押し続ける) > Le Bail + Fit Profile ! > File > Save
11. put fox2sflip.py and sflip-symmetry.txt in folder
12. fox2sflip.py case.xml 0.05 (win) or ./fox2sflip.py case.xml 0.05 (linux)
 最後の数値は FWHM で FWHM/2の範囲でピークを分けることを指定する。ディフォルトは 0.1 となっている。
  ( option [2])
#composition Mg8 Si4 O16 > composition your using atoms
#histogram composition > histogram composition
13. check wheter or not the lattice constants (a, b, c) and angle ( alpha, beta, gamma) in case.inflip are matching space group. 
  e.g. 89.9900 > 90.0000
14. SUPERFLIP or superflip.exe case.inflip
EDMA(未確認)
15. composition ### ### > composition your using atoms
16. EDMA case.edma
17. case.cif to VESTA
[3] http://fox.vincefn.net/
[4] Superflip_plus_Fox: http://www.misasa.okayama-u.ac.jp/~masami/pukiwiki/index.php?Superflip_plus_Fox
※ 同じ条件で何度か計算するが、フルホイスラー(FH)になったりハーフホイスラー(HH)になったりする。FHとなる確率は50%だった。
※ 高角度側を出来るだけ取り入れられないとMEMの結果には近づいていかない。
--------------------------------------------------------------------------------
■ Appendix 1 (Abinit での入力ファイル例)、γ-CuI

CuI.files
--------------
CuI.in
CuI.out
CuIi
CuIo
CuI
./I_fhi_pbe.psp
./Cu_fhi_pbe.psp
--------------

CuI.in
--------------
#************************************************************
#*      Generated by cif2cell 1.0.10 2014-04-07 21:03       *
#*                             ()                           *
#* Failed to get author information, No journal information *
#************************************************************
# Structural parameters
acell   3*11.5256284565

rprim     1.000000000000000   0.000000000000000   0.000000000000000
          0.000000000000000   1.000000000000000   0.000000000000000
          0.000000000000000   0.000000000000000   1.000000000000000
       
natom    8
ntypat   2
typat    1 1 1 1 2 2 2 2
znucl    53 29
xred     0.000000000000000   0.000000000000000   0.000000000000000
         0.000000000000000   0.500000000000000   0.500000000000000
         0.500000000000000   0.000000000000000   0.500000000000000
         0.500000000000000   0.500000000000000   0.000000000000000
         0.250000000000000   0.250000000000000   0.250000000000000
         0.750000000000000   0.750000000000000   0.250000000000000
         0.750000000000000   0.250000000000000   0.750000000000000
         0.250000000000000   0.750000000000000   0.750000000000000
#ionmov 3
#optcell 2
#ntime 10
#tolmxf 5.0d-4
prtcif 1
#nsppol 2
ecut 15
#pawecutdg 40
ixc 11
chkprim 0
#nsym 216
occopt 7
tsmear 0.01
kptopt 1
ngkpt 6 6 6
#toldfe 1.0d-5
toldff 5.0d-5
ecutsm 0.5 #
#nband 90
--------------

■ ノルム保存型擬ポテンシャルファイルの準備
計算するフォルダに I_fhi_pbe.psp と Cu_fhi_pbe.psp の名称のテキストを作成する。
Atomic data/pseudopotentials ( http://www.abinit.org/downloads/atomic-data-files ) > Norm-conserving pseudopotentials > Fritz-Haber-Institute (FHI) pseudopotentials for selected elements (Troullier-Martins scheme) > GGA (PBE) pseudopotential files from FHI code > 元素のボタンを押して、I_fhi_pbe.psp と Cu_fhi_pbe.psp にそれぞれの原子に対応するデータをコピー&ペーストする。
--------------------------------------------------------------------------------
■ Appendix 2
 神崎先生のHP(http://www.misasa.okayama-u.ac.jp/~masami/pukiwiki/にあるfp2sflip.pypythonスクリプトを書き換えてgrd2xplor.py を作成した。
 Dysnomia, etc → VESTA → *.grd → grd2xplor.py *.grd → *.xplor and *.edma
として、EDMA 用の入力ファイルを作成できる。
0. python と numpy をPCにインストールする。
1. Dysnomia などで 電子密度分布を計算する。
2. VESTA で 電子密度分布のデータ(*.pgrid など)を *.3grd に指定してSAVEする。
3. 計算したいファイルに、EDMA.exe と
sflip-symmetry.txt, grd2xplor.py(下記をコピーしてテキストなどにペーストし、名称はgrd2xplor.pyにする) を入れる。
4. コマンドプロンプト(shift + 右クリック)やTerminalなどを開いて、
 grd2xplor.py case.grd
と入力すれば、*.edma と *.xplor ファイルが生成される。
5. *.edma の composition ### ### を composition 単位胞にある原子の種類と個数を記述
 ( *.edma の終わりには、centerofcharge yes も記述してみると良いだろう)
7. EDMA *.edma で計算が実行される。
8. *.coo に charge が記述されている。そこから占有している原子を推測してみるのも良い。
※ 自作のスクリプトにおけるHFでの結果では、composition で 単位胞にある原子の種類を記述すれば、概ね想定通りの位置に原子が配置されていると結果が出力される。
※ 私の技術では、いまのところ占有率が1の結果しか得られていない。ドーピングを行った系などでは、ドーピングして入る原子の位置は正しい結果となった。しかし、ドーピングのある位置の占有率は誤って出力される。加えて、原子間にあるような占有率も上手くはいかない。
※ *.coo に記述される charge がどのような計算方法で出力されているかは、いまのところ私の勉強不足で不明だが、HFでは(理想的には)空孔サイトの位置でも charge が出力されている。

grd2xplor.py

#!/usr/bin/python
import sys,math
from numpy import *

# grd2xplor.py
# modified by **** from fp2flip.py (made by mkanzaki@me.com {Original version})
# This code reads .grd and .lst files produced by MEM analysis of
# Dysnomia (& VESTA) and RIETAN-FP, respectively,
# which containing D(I1,I2,I3) list by MEM, and
# the code then produces an input file for EDMA.
# To get cell parameters, wave length, and space group
# information, the code also reads .lst file too.
# Now you can also produce an input file for EDMA too.
#
# This is for Mac. For other OS, you may need to modify
# first line, depending on your python directory
# You may also need to change mode of the file, as follow
# > chmod 744 grd2xplor.py
#
# This code takes one parameter,
# .grd file name which Dysnomia (& VESTA) porduced.
# This code requires "sflip-symmetry.txt" in same directory.
# Original code can be downloaded from my web site.
# http://www.misasa.okayama-u.ac.jp/~masami/pukiwiki/
#
# For use, type following from terminal
# Here I assume this code and grd file are in current directory.
#  > ./grd2xplor *.grd
# You will find *.xplor and *.edma.
# You need to edit obtained .edma. Currently composition is
# commented out. So you can run EDMA,
# composition is necessary, so you need to
# edit .edma file before run EDMA.
# then run EDMA,
# > ./EDMA *.edma
#
# reading grd and lst file
if len(sys.argv) == 1:
 print 'No file name provided!'
 print 'For example, >./grd2xplor.py test.grd'
 print 'also see comment part of this code.'
 exit()
try:
 base = sys.argv[1].rstrip('.grd')
 f1 = base + '.lst'
 f2 = base + '.grd'
 file1 = open(f1,'r')
 file2 = open(f2,'r')
except IOError, (errno, msg):
 print 'lst file open error!'
 exit()

f4 = base + '.xplor'
file4 = open(f4,'w')
f5 = base + '.edma'
file5 = open(f5,'w')

# Extract title from .lst
for line in file1:
 if 'Title:' in line :
  s = line.index('Title:')
  ti = line[s+7:len(line)-3].replace(' ','')
  break
#file4.write(ti + '\n')
#file4.write('#\n')
file5.write('inputfile ' + base + '.xplor\n')

for line in file1:
 if 'Wavelength' in line :
  s = line.index('=')
  wl = line[s+1:s+10].replace(' ','')
  break

# Extract space group name from .lst
for line in file1:
 if 'Space group:' in line :
  s = line.index('Space group:')
  s2 = line.index('(')
  sg = line[s+13:s2-1].replace(' ','')
  #print sg
  break

# Finding refined cell parameters from .lst
rfine = 0
for line in file1:
 if '*** R factors,' in line :
  rfine = 1
  break
if rfine == 0 :
 print 'No refined parameters'

# Extract refined cell parameters from .lst
for line in file1:
 if 'Lattice parameter, a' in line :
  cel1 = line[16:30].replace(' ','')
  #print cel1
  a = float(cel1)
  break
for line in file1:
 if 'Lattice parameter, b' in line :
  cel2 = line[16:30].replace(' ','')
  #print cel2
  b = float(cel2)
  break
for line in file1:
 if 'Lattice parameter, c' in line :
  cel3 = line[16:30].replace(' ','')
  #print cel3
  c = float(cel3)
  break
for line in file1:
 if 'Lattice parameter, alpha' in line :
  cel4 = line[16:30].replace(' ','')
  #print cel4
  alpha = float(cel4)
  break
for line in file1:
 if 'Lattice parameter, beta' in line :
  cel5 = line[16:30].replace(' ','')
  #print cel5
  beta = float(cel5)
  break
for line in file1:
 if 'Lattice parameter, gamma' in line :
  cel6 = line[16:30].replace(' ','')
  #print cel6
  gamma = float(cel6)
  break
file1.close()

# write to these parameters to .edma
#file4.write('cell ' + cel1 + ' ' + cel2 + ' ' + cel3 + ' ' + cel4 + ' ' + cel5 + ' ' + cel6 +'\n')
file5.write('cell ' + cel1 + ' ' + cel2 + ' ' + cel3 + ' ' + cel4 + ' ' + cel5 + ' ' + cel6 +'\n')
file5.write('export ' + base + '.cif\n')
file5.write('# Chage composition below\n')
file5.write('# i.e., Mg8 Si4 O16\n')
file5.write('composition ### ###\n')
file5.write('\n')

# read sflip-symmetry.txt for space group
# change sflip-symmetry.txt directory if necessary
try:
 file3 = open('./sflip-symmetry.txt','r')
except IOError, (errno, msg):
 print 'File open error!'
 exit()
flag = 'n'
for line in file3:
 if line == '':
  print 'No space group match!'
  file3.close()
  exit()
        else:
  if sg in line :
   flag = 'y'
  if flag == 'y':
   file5.write(line.rstrip() + '\n')
          if 'endcenters' in line:
                  break
file3.close()

# read charge density list
for line in file2:
 file4.write('\n')
 file4.write('      1' + '\n')
 file4.write(line)
 break

for line in file2:
 LCA = float( line[ 1: 8].replace(' ','') )
 LCB = float( line[ 9:17].replace(' ','') )
 LCC = float( line[18:26].replace(' ','') )
 ANA = float( line[27:35].replace(' ','') )
 ANB = float( line[36:44].replace(' ','') )
 ANG = float( line[45:53].replace(' ','') )
 #print LCA, LCB, LCC, ANA, ANB, ANG
 break

for line in file2:
 PX = int( line[ 1: 5].replace(' ','') )
 PY = int( line[ 6:11].replace(' ','') )
 PZ = int( line[12:17].replace(' ','') )
 #print PX, PY, PZ
 break

file4.write("{: >8d}{: >8d}{: >8d}{: >8d}{: >8d}{: >8d}{: >8d}{: >8d}{: >8d}".format(PX, 0, (PX-1), PY, 0, (PY-1), PZ, 0, (PZ-1)) )
file4.write('\n')
file4.write("{: >12.4f}{: >12.4f}{: >12.4f}{: >12.4f}{: >12.4f}{: >12.4f}".format(LCA, LCB, LCC, ANA, ANB, ANG) )
file4.write('\n')
file4.write('ZYX')

D = zeros([PX+10, PY+10, PZ+10])
sum = 0
var = 0
esd = 0

I1 = 0
I2 = 0
I3 = 0
for line in file2:
 for I0 in range(6):
  den = line[ (1+15*I0):(15+15*I0)].replace(' ','')
  if den <> '':
   D[I1][I2][I3+I0] = float(den)
   #D[I1][I2][I3+0] = float( line[ 1:15].replace(' ','') )
   #D[I1][I2][I3+1] = float( line[16:30].replace(' ','') )
   #D[I1][I2][I3+2] = float( line[31:45].replace(' ','') )
   #D[I1][I2][I3+3] = float( line[46:60].replace(' ','') )
   #D[I1][I2][I3+4] = float( line[61:75].replace(' ','') )
   #D[I1][I2][I3+5] = float( line[76:90].replace(' ','') )
   #print D[I1][I2][I3+I0]
   #file4.write("{: >12.5e}".format( D[I1][I2][I3+I0]) )
   #file4.write('\n')
   sum += float(den)
 I3 += 6
 if I3 >= PZ :
  I3 = 0
  I2 += 1
  if I2 >= PY :
   I2 = 0
   I1 += 1
   if I1 >= PX :
    break

sum /= (PX * PY * PZ)

for I1 in range(PZ):
 file4.write( '\n' )
 file4.write("{: >7d}".format( I1 ) )
 I4 = 0
 for I2 in range(PY):
  for I3 in range(PZ):
   if ( I4 % 6 ) == 0 :
    file4.write( '\n')
   file4.write("{: >12.5e}".format( D[I1][I2][I3] ) )
   var += ( D[I1][I2][I3] - sum )**2
   I4 += 1

file4.write( '\n' + '-9999' + '\n')
var /= ((PX * PY * PZ) -1)
esd = sqrt(var)
file4.write( "{: >12.5e}{: >12.5e}".format(sum, esd) )
file4.write('\n')

#test
#print I4
#print ((PX-1) * (PY-1) * (PZ-1))

file2.close()

file4.close()
file5.close()


--------------------------------------------------------------------------------
QRコード
携帯用QRコード
アクセス数
ページビュー数