DVXa + DVME

 ここではDVXa関連の情報について纏める(この項目を立ち上げた理由は『はじめての電子状態計算』がamazonで購入できず、また中古の価格が原価の3倍程度になっているためです。利用者や注目してくれる人々を増やして改訂・増刷して頂けることを目指し、色々と情報を記述していきます)。
------------------------------------------------------------------------------
■ 基本的な情報
  書籍にあるDV-Xa法 SCAT プログラムは非相対論版である。以下のプログラムはDV-Xa研究協会や夏の学校などに参加したり、ディベロッパーにお願いして入手することができる(原子番号が50を超える原子を扱う場合は相対論版が良いだろう)。
A. 相対論版DV-Xa法プログラム(相対論版 dvxar)
B. 遷移確率計算プログラム(sxs及びsxsplt)
C. 全エネルギー計算プログラム(TESDA及びCoulomb)
D. 対称軌道作成プログラム(makef25)
E. 配置間相互作用を考慮したプログラム(DVME法)[DMEX1]
F. CVV Auger スペクトルの計算プログラム(入手可能かは不明)
※ Linux PCの方が性能が良い場合は、wineをインストールして、下記と同様にすれば良い。
wineの起動コマンドは[wine cmd]となる。$HOME/.wine/drive_c が主に使用するディレクトリになる。
------------------------------------------------------------------------------
■ セットアップ
1) 黄色い本に記載のあるようにデータを配置する。
2) http://chem.sci.hyogo-u.ac.jp/hajimete/download.html から最新のデータを得る。
3) 解凍して同名のファイルに中身をコピーする。

■ SCF計算
1) setdvxa.bat のあるファイルで、shift を押しながら右クリックして「コマンドウインドウをここで開く(W)」を選択する。
2) setdvxa
3) cd calc
4) mkdir 計算したい系の名称
5) cd 計算したい系の名称
6) notepad f01
  入力ファイル f01 を作成する。
7) makef25 (makef25 が無い場合は飛ばす)
8) makef05
9) f06zp を他のファイルからコピー&ペーストし、 ncycm の値を50 など大きな値にする。
10) dvscat
以上で基本的な使い方は終了である。後は黄色い本の裏表紙の裏に計算したい事項の記載(DOS用の入力ファイルの作り方など)や入力するコマンド(これだけでよい)が記載されている。

■ 重要な計算
□ Mulliken population analysis
1) 上記のSCF計算を終わらせる。
2) popanl
3) f08p に結果が記載されている。
※ 共有結合の強さを評価する助けとなる。

□ Bond Overlap population analysis
1) 上記のSCF計算を終わらせる。
2) bndodr
3) bn8 に結果が記載されている。
※ 軌道間の相互作用を評価する助けとなる。

□ Net charge
1) 上記のSCF計算を終わらせる。
2) netc
3) i08 に結果が記載されている。
※ イオン性と共有結合性のどちらが強いかを評価する助けとなる。

□ 分子軌道の波動関数と電子密度のプロット
1) 上記のSCF計算を終わらせる。
2) wavnum3
3) makec04
 c04 と c05 が作成される。
4) c04 を編集
  2行目:縦軸 横軸のメッシュ数
  3行目:メッシュの開始点
  4行目:横軸の間隔
  5行目:縦軸の間隔
  6行目:プロットする軌道の数 軌道の番号(プロットする軌道の数で入力した数だけ入力)
5) c05
  6行目以降(最後の行を除く):プロットする等高線の値
6) contr
7) cmap
8) dvplot で、 Wに番号が付いたものと、chg を指定する。
※  https://www.jstage.jst.go.jp/article/kenbikyo1950/35/3/35_3_221/_pdf などの例がある。
スピン密度の差を議論したい場合
9) 上向きスピンは chg をプロットする。差スピン密度は drho を用いる。
10) c04 の 6行目の最初の数値を1にする。
11) drho cXX cYY cZZ
 XX, YY, ZZ には数値が入る
12) contmap xZZ dspin
13) dvplot dspin を指定する。
※ VESTAを用いたプロット方法は http://www.chem.ous.ac.jp/~gsakane/HidemaruDV/HidemaruDV.pdf に非常に詳しく記載されている。

□ Overlap population analysis
1) 上記のSCF計算を終わらせる。
2) lb4s と lb5s を他のファイルからコピー&ペーストする
3) lb4s の編集
  1行目:原子間での結合を計算する組の数(2行目と3行目の記述を繰り返す数)
  2行目:計算する(軌道の)ペアの数 得られた結果を割る数
  3行目:上記の1項目に入力した数だけ(軌道の)ペアとなる  原子の番号 - 原子の番号 を記述する。
  4行目以降:2行目と3行目を、1行目で入力した数値 - 1 だけ繰り返す。 
4) lb5s の編集(特に編集する必要はないが一応記しておく)
  1行目:プロットのエネルギー範囲 [eV]
  2行目:x, y 軸の長さと目盛り間隔
  5行目:半値幅
  6行目:プロットのオプション(黄色い本参照)
  7行目:値の最大値
  9行目:プロットの順番
5) lvlbnds
6) dvplot で lb7s を指定する。
※ T. Mizoguchi et al., Phys. Rev. B 61 (2000) 2180.では、遷移状態計算を行った後に、Overlap population analysis を行って非占有準位でのピーク構造がどのような結合から形成されているのかを議論している。

□ 遷移状態計算[10, TSX1, TSX2, TSX3, TSX4, TSX5]
・光電子スペクトルのピークエネルギーはSlaterのtransition state methodを用いると精度良く計算することが出来、少し慎重に計算すると一般には内殻レベルに関しては1%以内、価電子領域では0.1 eV程度の正確さでピークエネルギーを再現できる[TSX2]。
・NiOの酸素K殻ELNESの理論計算を行った結果、このスペクトルについては多体効果を考慮しなくともNi同士の軌道間相互作用を考えることで解釈が可能なことが示される[TSX5]。
・下記の論文を参照するとよい。
https://www.jstage.jst.go.jp/article/kenbikyo1950/35/3/35_3_221/_pdf 
・ 遷移状態計算を行った場合、小さなクラスターでも空孔を入れた状態における波動関数が、大きなクラスターで計算したものと同様の波動関数が形成される場合は、大きなクラスターを用いなくとも実験をよく再現できる。
・ K端やL1端のようにp軌道成分だけのPDOSを用いる場合には吸収断面積(PACS; photo-absorption cross section)までの計算が必要ない場合がある。
・ 遷移エネルギーを計算する場合はスピン分極を考慮して計算を行う。
・ 多重散乱法 FFFT8と第一原理計算、実験結果との比較もされている。
計算手順
1) f01 を作成
2) makef25
  Select Symmetry では上記に表示された数値(対称性を表す記号の左の数値)を入力する。
  Select orbital では、ランタノイド系列以後の重い元素でない場合、spdである 1 を入力する。
  NEW is reset for ... と表示されるので、対称性に対応させてNEQをリセットしたい場合は Yesの 1を選択する。
3) makef05
4) dvscat
  収束するまで dvscat を繰り返す。dvscatを押すのが面倒な場合は f06zp を編集し ncycm の下の数値を200などと大きくしておく。
  後は、http://www.chem.ous.ac.jp/~gsakane/STS.html に従えばよい。
5) f08e を見て電子を0.5取り去り、0.5付けたい軌道の番号(2列目)と対称性(3列目)を決める。
  f08 の一番右側が電子の占有数である。同じ対称性の占有軌道から同じ対称性の非占有軌道をしていすると良いだろう。
  wavnum と コマンドプロンプトに入力して、出力されたファイルも参照すると良い。
6) 新しいディレクトリを作成し、先ほど計算した f17, f25, f26 を新しく作成したディレクトリに入れる。
7) f26 の6行目を編集。6行目の順番は下記のようになる。
  ノンスピンの場合:TL LSY1 REMV1 LSY2 JLV2 REMV2 LSY3 JLV3 REMV3
  スピンの場合:TL KSP1 LSY1 REMV1 KSP2 LSY2 JLV2 REMV2 KSP3 LSY3 JLV3 REMV3
  ここで、 TL:全電子数、
  KSP1,KSP2, KSP3:スピンがup なら 1、スピンがdown なら 2
  LSY1, LSY2, LSY3:f25に上から順番に対称性の記号に番号が振られていると考える(合計するとf05の中程{bohr単位で座標が書かれている部分の上に記述している}に記載されているnsymの値と一緒になることが分かる)。電子を動かす対称性の番号を記入する。
  JLV1, JLV2, JLV3:f08e に記載されている電子数を変化させたい軌道のエネルギー(6行目)と対称性(3行目)の記号が記載されている 2行目の数値を入れる。
  REMV1, REMV2, REMV3:電子の増減量
8) f26を書き換えたら、dvscat と入力する。
9) 遷移エネルギー = 遷移状態法での全エネルギー - 遷移する前での全エネルギー
  1%程度の誤差が生じる。実験スペクトルと比較する場合には、横軸(エネルギー)を、実験値のピークと合うようにシフトさせても良いだろう。
XMCD: http://tdl.libra.titech.ac.jp/hkshi/xc/contents/pdf/117076547/10

□ 紫外可視吸収スペクトル[UVX1]
A. PDOSを計算する。
B. sxs または rsxs の計算を行う。
※ H. Yoshida et al., Jpn. J. Appl. Phys. 45 (2006) 146. ではPDOSと rsxs の双方の結果が記載されている。

□ XAS&ELNESスペクトル[12-20]
  DV-Xa研究協会会報、Vol. 9, No. 2 (1996) で吸収スペクトルの報告が多い。
A. 遷移エネルギーの絶対値はSlaterの遷移状態法により求める。
  上記の遷移状態計算を参照。
B. 遷移確率<f|r|i>を求める。
  sxs または rsxs の計算を行う。rsxs の計算方法は相対論版DV-Xa計算プログラムに詳細な記載がある。下記は sxs の場合。
  1) 上記のSCF計算を終わらせる。
    余談であるが、dvscat でSCF計算して出力された f09とf39を用いる。
  2) s04 を編集する。
    1行目:IATM  JPRT (2I5の形式で入力)
       IATM は 計算の中心となるNEQの番号。つまり<f|r|i>のrの中心を示す。
       JPRT は core orbital の属する原子の番号のオプション。0が通常表示。1,2は各x, y, zでの遷移確率を出力させる。
    2行目:LNLN  E0 (I5, F15.7 の形式で入力)
       LNLN は NEQの番号 主量子数 方位量子数 0
       E0 は core orbital のLOCA+NNLZcore orbitalのエネルギー固有値(Ry単位)。小数点を付ける。
       ※ 2行目の最初は分かりにくいが、左から詰めて記述する。NEQが1桁の場合は最初は空白にし、NEQの後にはスペースを入れずに主量子数と方位量子数を書き入れる。そして、その直後に0を記入する。この後は、推定されるエネルギー準位(Ry単位)を記入する。
  3) copy d04 s05
  4) sxs
  5) sxsplot
  6) copy d04 sp4
    IGL=GAUS または LORE、 IAE=ABSO または EMIS、 NSPIN=1 または 2
    E0, E1, DX は d04 または x04 と同じ
    SIG(I), I=1, NCORE は d04 または x04 と同じ
    NPRT, IPRT(I), I=1, NCORE は d04 または x04 と同じ
    FCS(I), I=1, NCORE は d04 または x04 と同じ
    ※ 難しく書いてしまったが、要するに、sp4はd04と同じである。違うのは、1行目の第1項で、GAUS を下記に書き換える。
    発光の場合:GAUSEMIS または LOREEMIS
    吸収の場合:GAUSABSO または LOREABSO
  7) copy d05 sp5
  8) pltsxs
  9) dvplot で sp7 を指定する。
  ※ 出力されたファイルで、 s08 の一番右の数値が遷移確率を表す。ssp7 は表計算ソフトで結果を開ける。
C. X線吸収スペクトルの計算(sxsの計算の詳細は上記 遷移確率<f|r|i>を求める を参照)
  1) 遷移状態計算 を行う。
    遷移を開始する電子の占有軌道を -0.5 し、 電子が移動する先の非占有軌道を 0.5 増やす。  
  2) s04 を編集し、sp4 と sp5 を準備して、sxs と入力して計算する。
  3) sxsplt
  4) pltsxs
  5) dvplot で sp7 を指定する。
D. DVME法により求める[MEX1, MEX2]。

□ 発光(PES or XES)スペクトル
A. half-filled MOでの計算[TX1]
  1) 下方遷移することになる軌道を0.5個にして計算する。飛ばされる内殻は満たされたままにする(計算方法は上記の遷移状態計算を参照)。
B. SXS-TS (蛍光X線スペクトル)(sxsの計算の詳細は上記 遷移確率<f|r|i>を求める を参照)
  1) 遷移状態計算 を行う。
     この場合に、遷移すると考える占有軌道を二つとも -0.5 とする。
  2) s04 を編集し、sp4 と sp5 を準備して、sxs と入力して計算する。
  3) sxsplt
  4) pltsxs
  5) dvplot で sp7 を指定する。

□ Auger electron energy [AX1, AX2]
A. 計算方法[AX2]
  EAuger = |Ei - Ef| ≒ |E(A) - E(B) - E(C) - EC(BC) - EE(BC)|
  EAuger: Auger electron energy (オージェ電子エネルギー)
  Ei: Total energy of initial state(始状態の全エネルギー)
  Ef: Total energy of final state(終状態の全エネルギー)
  E(A): Binding energy of A orbital in ground state(軌道Aの基底状態エネルギー。軌道Aは電子線やX線によって励起されて最初に正孔が出来る軌道)
  E(B): Binding energy of B orbital in ground state(軌道Bの基底状態のエネルギー。軌道Bは電子が軌道B→軌道Aとなる軌道B)
  E(C): Binding energy of C orbital in ground state(軌道Cの定常状態のエネルギー。軌道Cは電子が軌道B→軌道Aへ落ちる時にエネルギーを得る軌道)
  EC(BC): Coulomb integral between B and C orbital
  EE(BC): Exchange integral between B and C orbital
  文献によっては EC(BC)及びEE(BC)を終状態での正孔同士の反発又は相関エネルギーとして一括して記述している文献もある[TX1]。

B. Augerスペクトルの計算法[AX3, TX1]
a. Core-Core-Valence (CCV) Auger 遷移[TX1]
  計算式: ACCV(E) = Cs ρ's(E) + Cp ρ'p(E)
  ここで、ρ'l は正孔のある終状態の局所状態密度分布、Clは原子のAuger遷移行列要素。lはsまたはpを示す。
  DV-Xa法SCATでの計算方法[AX3]
  1) SCF計算を行う
  2) maked04 でd04とd05を作成し、x04とx05にそれぞれ名称を変更する
  3) HOMO を 0 eV にする(2行目の右端を1にする)
  4) エネルギー範囲と半値幅は実験スペクトルと比較できるように適切な値を記入する
  5) Augerに関係する原子のみを選択し、考えられるAuger遷移に関連するValenceの原子に(DOSをプロットするときのように0以外の)数値を記入する。それ以外のものは0へ。
  6)  倍率を記入するところには、考えられるAuger遷移の比率をデータテーブル[DT]を参考に記入する。
  7) xps のプログラムを動かす
  8) 出力されたx07をdvplotで、xx7をigorやsma4winでプロットする
  9) 実験結果ではプラズモン励起によるエネルギー損失サテライト線を差し引いて、計算結果と比較してみるとよい
  [DT] M. H. Chem et al., Atomic Data Nucl. Data Tables 45 (1990) 1.
b. Core-Valence-Valence(CVV) Auger 遷移[AX4]
  計算式:ACVV(E) = Css R2s ρs(E)∗ρs(E) + Csp Rs Rp ρs(E)∗ρp(E) + Cpp R2p ρp(E)∗ρp(E)
  ACVV(E) 〜( ρ∗ρ'(E) ) / ( [1 - ΔU ∫( ρ∗ρ(E) / (E-ε) ) dε]2 + [ΔUπρ∗ρ'(E)]2 )
  ここでの ∗ は畳み込み積分を示す。 
  1) f08またはf08pから手作業で計算
c. half-filled MOsでの計算[TX1]
  1) Auger遷移に関係する軌道全てを0.5個ずつにして計算する(計算方法は上記の 遷移状態計算を参照)。

□ その他
■ 鏡面対称性の作成の仕方1
1. displatでZ軸が画面に対して垂直であるようにしてf01を出力する
2. symOrbでsを押して探索させCsが表示される
3. f01を出力する
4. f01のNEWを見やすいように並べて、同じNEWのXYZの値を揃える
5. makef25
6. makef05
7. dvscat

■ 鏡面対称性の作成の仕方2(鏡面対称性のあるf01がmakef25で自動認識されない場合など)
1. displatでZ軸が画面に対して垂直であるようにしてf01を出力する
 (displatで鏡面対称性を作ると入力ファイルに対して少し安心できる)
(ここでNEWが順番になるようにコピー&ペーストしておくとf25を比較しやすくなるが、f25.inでのチェックは難しくなる)
2. Excelなどでf01を読み取り、angstromの場合は1.889725989倍してbohr(atomic)単位にする。桁数は揃えておきたい(Excelのボタンにある)。
3. f25.inの名称のファイルを作る(makef25sでzに対して+と-の座標ができる)。Excelのデータをコピー&ペーストする。置換でTabを空白に変えておく。
Cs
 考慮に入れる軌道(第四と五周期の原子があるならばd、第六周期があるならばf)
bohr(atomic)単位で座標を入力する(鏡面に対する片方のみを記述する)
4. makef25sを入力し、作成されたf25を確認、f25の形式と同じ順番となるようにf01を確認する
5. makef05
6. dvscat
------------------------------------------------------------------------------
■ プログラムの入手法[DX00]
□ DV-Xa分子軌道法プログラム本体(次世代DVSCAT実行ファイル群、DISPLAT一式)
  http://chem.sci.hyogo-u.ac.jp/hajimete/download.html
□ 教育用分子軌道計算システム eduDV(Fortranプログラム群、バッチファイル群)
  http://www.chem.ous.ac.jp/%7Egsakane/fun/index.html#edudv
□ 結晶構造及び電子・核密度等の三次元可視化プログラムVESTA
  http://jp-minerals.org/vesta/jp/
□ 構造ファイル変換ユーティリティ Open Babel
  http://openbabel.org/wiki/Main_Page
□ 秀丸エディタ(シェアウェア)
  http://hide.maruo.co.jp/software/hidemaru.html
□ DV-Xa法計算支援環境(秀丸エディタマクロ集)
  http://www.chem.ous.ac.jp/~gsakane/dvxa_assistance_environment_7.html
□ DV-Xa法計算支援環境利用の手引き(PDF文書)
  http://www.chem.ous.ac.jp/~gsakane/dvxa_assistance_environment_1.html
□ DV-Xa関連のページ
  http://imac.eng.kagawa-u.ac.jp/dvxa/dvxa.html 
------------------------------------------------------------------------------
■ upgrade(詳細は文献[4]を参照)
1) CDからデータをインストール
2) 「はじめての電子状態計算ホームページ」での「ダウンロードページ」から dvxa_v1_04 [DX0]をダウンロードし置き換える。displat_fullやsymchkなども置き換えると良い。
※ dvscat2を利用したい場合は、windows7 64 bitを利用している場合は DV-Xa法に関するブログ にて dvscat2 を書き換える。
------------------------------------------------------------------------------
■ Linux上でのセットアップ(編集中)
------------------------------------------------------------------------------
■ 入力ファイルの変換(詳細は文献[4]を参照)
  「秀丸エディタ」を有していれば文献[4]によりVESTAからDV-Xaの入力ファイルのf01を簡単に得ることが出来る。ここでは「秀丸エディタ」を用いずに無料で入力ファイルを変換する方法を解説する(予定)。
------------------------------------------------------------------------------
DV-Xa法による新材料設計[APX1]
□ 基本的な材料設計指針:
  伝導性や磁性、光応答性など物性測定結果をもとに、分子軌道計算と分子動力学計算を行い、材料合成にフィードバックする。

□ 伝導性の良い材料を設計するためには:
  エネルギー準位計算を行い、易動度(mobility)を高くするための指針(電子相関)、HOMO-LUMOギャップの減少(半導体の場合)、キャリア濃度の増大(金属の場合)など。

□ 新しい超伝導体を得るためには:
  これまでのクーパー対とは異なる考え方。互いに反発する電子間に引力を持たせるためにスピン(反強磁性)を用いる。d-π相互作用を用いた超伝導材料の設計。

□ 光応答性の高い材料を設計するには:
  DOSの議論(光の吸収・発光)。金属錯体における色のメカニズム(d-d, d-π, π-π*, M→L遷移など)。遷移状態法によるエネルギー計算。

□ 熱電性能の高い材料を設計するためには:
  熱電性能指数Zを大きくする。具体的にはゼーベック係数Sを大きくし、電気伝導度σを大きくし、熱伝導を小さくする。そのためには半導体が良い。格子振動の寄与も重要。

□ 材料としての金属錯体の特徴:
  高い設計性と高い機能発現性を同時に兼ね備える。特殊な酸化状態・電子状態を有し、低次元物性材料としての応用も可能。

□ 強い磁性材料を設計するためには:
  スピン間に働く相互作用の調整と、バルクにスピンを同じ向きに合わせるための単結晶合成技術が必要。さらにVESTAを用いて波動関数を三次元可視化することにより、視覚的にスピン間の相互作用を予測することが可能。
------------------------------------------------------------------------------
■ 便利な情報
A. 『はじめての電子状態計算』の裏表紙の裏には便利なコマンド表が記述されている。
B. RIETAN-PFとDV-Xa法SCATは「秀丸エディタ」との併用が便利である。
C. 収束条件を0.01 eとした報告もある[12]。(計算が収束しなければ検討してみても良いだろう)
------------------------------------------------------------------------------
Q&A [QX1, QX2]
A. クラスターのとり方について
  共有結合が強く、group unitがハッキリしている場合は、それを途切れないようにすることが重要であろう。従って、中心原子から一定距離で切った球形のクラスターよりもgroup unitを数個つないだクラスターが良い結果を生むようである。(下記の「C. クラスター内の電子数の決め方」も参照)

B. クラスターの終端の取り扱い
  水素終端(クラスターの外側に水素を配置する)では、水素化合物の電子状態を計算していることになっており、水素のレベルが価電子帯あるいは伝導帯に混成してしまう。半導体のバンドギャップの計算にはほとんど問題はないが、X線蛍光あるいは吸収スペクトルの計算では、注意しないと水素と混成したレベルが重畳してしまう。

C. クラスター内の電子数の決め方
  全ての化合物のクラスター計算において、電子数は化学ポテンシャルを表現しているので、原則として、中性のクラスターをとるのではなく、適当な電子数のクラスター(イオン)で計算すべきである。しかし、多くの金属性金属間化合物のように、構成原子間での電荷移動量が少なく、電子数の微少な変化によって全体の電子状態が大きく変化しない場合は、あまりこの電子数に神経質になる必要はないだろう。また大きなクラスターをとると、総電子数を中性にとっていても、クラスター表面を一皮むいたクラスターの内側での電子数が、実際の物質の状態をよく反映している場合もある。(計算する系の価数をよく考慮して電子数の取り方を決める。同じクラスターを用いても、電子数を変化させることで、電気伝導性のあるTiO{クラスターでは(TiO6)10-あるいは(TiO6)10+}と絶縁体TiO2{クラスターでは(TiO6)8-}を表現できる)(私はあまり頭が良くなかったので、30〜100原子程度の大規模なクラスターを構築して計算し、中心部分の電子状態をバルクの電子状態とした)

D. 不純物を挿入した場合、クラスターの大きさは影響するのか?
  浅い準位を発生するドナーやアクセプター不純物を挿入した場合、不純物準位の波動関数が大きく広がっているため、クラスターの大きさは電子状態に影響を与える。深い準位を発生するような不純物を挿入した場合、不純物準位の波動関数は比較的局在しているためクラスターをあまり大きくとらなくても電子状態への影響は小さいと思われる。しかし、シリコン母体そのものの波動関数が非局在化しており、クラスターのサイズによって影響を受けていることに変わりはない。不純物原子の周囲で格子歪が生じている可能性が高いため、このことも考慮する必要がある。

E. SXSの精度はサンプル点数にどう依存するか?
  配布版SXSでは、軽原子1個あたり、約6000点のサンプル点を用いて、積分値が10〜20%ばらつく。

F. 高いAOのpupulation が負になる場合があるが、それで良いのか?
  負のorbital populationは、その原子軌道が分子軌道のなかで、おもに反結合性の寄与をしていることを示していて、とくに気にする必要はない。結晶を調べるため、負に帯電したクラスターを考える場合は現れやすい。解析上扱いにくい場合は、それら軌道を除くのも方法である。
------------------------------------------------------------------------------
DV-Xaの新展開[DX1, DX2, DX3]
------------------------------------------------------------------------------
■ 理論的背景[7, 8, TX1, TX2]
A. 一電子近似(または平均場近似)
  多体波動関数 ψ(r1, r2, r3 ......, rN) を ψ1(r1) × ψ2(r2) × ...... × ψN(rN) とした一電子の波動関数の掛け算に近似する。当然ながら両者は「イコール」では結べないので、この近似に対処することが必要となる(この対処方法が後に述べる「密度汎関数理論」である)。

B. 時間に依存しないSchrödinger方程式(原子単位系)
  ハミルトニアン H = Σm( -1/2▽2m - ΣN ZN/riN) + Σi≠j 1/rij + ΣN≠l ZNZl/rNl
  エネルギー E = <ψ|H|ψ>
ここで、mは全粒子(電子と原子核)に対する添字、iとjは電子の添字、kとlは原子核の添字、ZNは原子Nの原子番号、riNとrijとrNlは二個の粒子間の距離である。最初の項が運動エネルギー項、二番目の項が電子-核間の引力、三番目の項が電子間反発、四番目が核間反発に相当する。

C. Born-Oppenheimer(断熱)近似
  陽子は電子の約1860倍の質量を持つために、電子から見れば陽子はほとんど動いていないように考えることが出来る。電子と核を分離すれば下記のようになる。
  ハミルトニアン H = Σi ( -1/2▽2i - ΣN ZN/riN) + Σi≠j 1/rij + ΣN (-1/2▽2N)  + ΣN≠l ZNZl/rNl
 ここで、iは電子に対する添字。 四番目の核の運動エネルギー項と五番目の核間反発は核の位置を固定したので、繰り返し計算(SCF計算) のどこか(最初や最後)で計算してしまえば定数として扱うことが可能となり、計算コストを削減することができる。
  もし、陽子の位置を固定していなければ計算時間はとても掛かるのだけれど、ハミルトニアンHから、核の振動状態、回転状態、並進状態に関する情報を得ることが出来る。つまり分子の振動スペクトルを予測することが可能となることは覚えておかれたい。

D. Hartree-Fock方程式
  Schrödinger方程式に変分法を適用し、ラグランジュの未定乗数法を用いると下記のHartree-Fock方程式を得ることが出来る。
  hi ψi(ri) = εi ψi(ri)
  ハミルトニアン hi = -1/2▽2 - ΣNZN/riN + Σj ∫ (ψ*j(rjj(rj)/rij ) drj + Vxc(ri)
ここで、最後の項は記述が複雑なので省略する。この最後の項は次に述べる「密度汎関数理論」によって、様々に記述される。

E. 密度汎関数理論
  ホーエンベルグとコーンにより、背理法を用いて「電荷密度の分布とエネルギー」が1対1の関係を持つことが証明された。これにより、一電子近似(または平均場近似)やBorn-Oppenheimer近似などの近似を用いても、厳密解と同じ「電荷密度の分布とエネルギーの関係」を与えることができれば良いことが分かった。そこで、複雑な計算となる交換・相関項を簡単な形に書き換えることが行われている。例えば、交換相関項は波動関数では無く電荷密度分布を関数としたものにするなどである。[XCX1, XCX2]

F. 数値計算
  分子軌道を原子軌道の線形結合として下記のように表してみる。
  ψi(ri) = Σm Cmi Φm(ri)
  これをHartree-Fock方程式やHartree-Fock-Slater方程式などに代入して展開し、左から原子軌道を掛けて積分すると、原子軌道は直交基底なので、次の永年方程式を得ることが出来る。
  (H - εS)C = 0
  ここで、Hは共鳴積分、εは固有値、Sは重なり積分、Cは固有ベクトルである。この永年方程式を手作業では数学の授業で習った固有値を解く手法で求めることができるが、計算機ではそれに適した固有値と固有ベクトルを求めるアルゴリズムを用いて解く。また、HやSでは積分を行う必要があるので、Σで代用してある程度の実空間の点で済ませる。この点の取り方がDV-Xa法ではDVサンプルポイントの取り方となる。この部分は並列化が可能なため、誰かが並列化用に書き換えてくれると嬉しい。
※ 余談であるが、並列化の例であれば、DV-Xa法SCATをFortran77からC++に書き換えてMPI並列した報告がある[MX1]。
------------------------------------------------------------------------------
■ 重要なプログラム・アルゴリズム解説[7]
------------------------------------------------------------------------------
■ プログラミングでの重要な点
A. DV-Xa法 SCAT は Fortran77形式で記述されている。
B. Fortran77形式ではcommon文を用いてルーチン間のデータをやり取りすることが主流であった。
C. common文では旧帝などでの大型計算機の場合などを除き、多くの場合、メモリーは2GB程度までしか利用出来ない。
D. どのプログラムでも、配列やcommon文でのデータの利用範囲には注意が必要となる。また、データの書き出し、読み出しにも注意が必要となる(小規模の系で計算が可能な場合はそのバグの多くがこの問題であることが多い)。
E. 上記Dの問題に対処するには、「徐々に系の規模を大きくしてどの時点でエラーが表示されるかを知る」及び「その時点でのデータのチェック」、「初期値の設定忘れが無いかチェック」、「オプションを付けてコンパイルを行う」ことである。
F. コンパイラでのオプションは左欄にある「コンパイラの設定」を参照するとよいだろう。
------------------------------------------------------------------------------
■ サブルーチン[7]
1) MAINS
  このプログラムのメインルーチンである。計算を実際に行う下記の2)〜6)までの各サブルーチンを呼び出している。

2) DINTAL
  □ 計算に必要なデータの設定
  □ MAINSから呼ばれる。計算を行おうとする系についての各種のデータを設定するルーチン。計算に必要なデータは(一部は別のサブルーチンで読まれる)ここで読み込まれたり、計算を行って設定される。サブルーチンDINTALは主としていくつかのread文とサブルーチンのcall文、その他から成っている。

3) ATCALC
  □ 原子軌道の計算
  □ MAINSから呼ばれる。クラスターを構成する原子のうち種類が(原子の種類のみならず、空間的性質についても)同じ原子、すなわちreference原子についての諸量を計算するルーチンである。結果はcommon文 PRCORE、PROND、PRMESHでやり取りされる。

4) CATOM
  □ 分子軌道の繰り返し計算
  □ SCF分子軌道の繰り返し計算の部分である。サブルーチンの概要は以下の通りである。
  ポテンシャル用原子を加える場合、その加えるポテンシャル用原子のクーロンポテンシャルを集計しておいてからSCF計算に必要な必要な諸量をread文で読み、基底関数をセットして(SSOGT)、Fock行列及び重なり行列を形成し(FOCKM)、対角化問題を解き(EMERGE)、MullkenのPopulation Analysisに従って密度解析を行い(CHRANL)、SCF収束判定をして(JUDGES)、次のステップのデータを整える(CHMIXS)、といった手順で繰り返し計算を進める。

5) POPANL
  □ Mullikenの電荷密度解析
  □ 分子軌道が求められた後、密度解析を行うルーチンである。このルーチンを通って全体の計算は終了するので、次に引き渡す結果となるものはない。また、必要なデータはcommon文から得られる。密度解析及び結果の印字は各対称ブロックごとに行われる。実際の計算は、ループ内でサブルーチンを呼んで、そこで行われている。

6) SFSTOR
  □ ファイル操作
※ 上記の2)〜6)までの各サブルーチンでは更に下位のサブルーチンを多数呼び出している。
------------------------------------------------------------------------------
■ DINTALとそれから主に呼ばれる下位のサブルーチンの詳細
□ DINTAL
  read文で読み込んでいる内容は以下の通りである。
    iout: 出力ファイル場号。標準出力は6番であるが、それ以外(尚且つ入力等の既に使われている番号以外)を指定することによって、計算結果をファイルに落とすことができる。
    mmol: クラスターの名前
    rido: 原子軌道の計算において加える井戸型ポテンシャルの半径(a.u. 単位)
    vido: 原子軌道の計算において加える井戸型ポテンシャルの深さ(a.u. 単位)
    idoa: 原子軌道の計算を行う場合に井戸型ポテンシャルを加えるかどうかのパラメータ。0のときは加えない。1, 2はそれぞれ、タイプ1、タイプ2である。タイプ1は距離が遠くなるにつれ値が減少していくポテンシャル。タイプ2は距離が遠くなるについて値が増加していくポテンシャルである。
    macf:
    ngrd:
    ncycle: SCF分子軌道計算をする際の繰り返しの上限
    ntabpt: 原子軌道計算をする際の原子のmesh数
    nnst:
    ltime: CPUの上限(msec単位)
    ntape:
    jprt: 出力のオプション。-4から5までの値を判定に用いている。
    tolm: SCF分子軌道計算でのthreshold
    nat: クラスターを構成する原子の個数
    ndat: クラスターを構成する原子のうち原子軌道計算をする原子、すなわちreference原子の個数
    ncore: reference原子の原子軌道の和
    natoms: ポテンシャルを含めたクラスターを構成する原子の個数
    ntypes: ポテンシャルを含めたクラスターを構成する原子のうちのreference原子の個数
    mq: Madelungポテンシャルを加えるかどうかのパラメータ。0で加えない、1で加える。
    jptg: DVサンプルポイントを計算するかどうかのパラメータ。0でない値が入っている場合、DVサンプルポイントの計算は行われない(サンプルポイント及び球面調和関数は既に分かってファイルにある場合)
   neq(i): クラスターを構成する原子のreference原子番号
   z(i): クラスターを構成する原子の原子番号
   tl: クラスターの全電子数
   lsy1, jlv1, remv1: 遷移状態計算を行う場合の電子の数を変化させる分子軌道の属する対称ブロックと分子軌道のレベル及び変化させる電子数である。remv1は-0.5で0.5個電子を放出し、+0.5で0.5個電子を受け取ることになる。
   lys2, jlv2, remv2: クラスターのイオン化エネルギーの計算の場合、電子の数が変化する分子軌道は1つであるが、1電子励起状態の場合は2つになる。
   lsy, jlv3, remv3:
   itime: 前回の計算はどのサブルーチンで終了したか? 今回はどこから計算するかを示すパラメータ(CPUを制限しているために終了した場合)
   jout(i,j): 原子軌道に関するデータ。一般に0である。1の場合、原子軌道のbaseから除く。2の場合、最後の結果の出力に際してこの原子軌道については出力しない。

  DINTALからはいつかのサブルーチンを呼んでいる。呼んでいるサブルーチンとそのサブルーチンで行われていることの概略は以下の通りと成る。
  □ FREAD 入力データをフォーマット等は気にせずに単なる文字列としてそのまま別のファイルにコピーする
  □ CONSNT 定数の設定
  □ FILGEN ファイル番号の設定
  □ INITIA 配列の初期化
  □ DINPT 各reference原子についての変数の読み込みと設定
  □ RHOB 原子の電荷密度の初期値の読み込み又は計算
  □ DIMCHK 各データの次元のチェックとファイルへの書き込み
  □ CORES 原子軌道についてのいくつかの変数の設定
  □ SYMTRY 対称軌道を用いる場合はその構築又は原子軌道に関する情報のまとめ
  □ ATMPOS 原子の位置や係数について
  □ MDINPT mq=1 の場合呼ばれ、Madelungポテンシャルについてのデータを読み込む
  □ PTGET DVサンプルポイントと各点における各原子軌道の動径関数を計算する

  □ FREAD DINTAL及びCATOMにより呼ばれる。入力データをファイルに書き込むルーチンである。DINTALで呼ばれる場合、引数のiparaは1であり、行番号30で文番号600に飛ぶ。標準入力データファイルを行うごとに最後の行まで呼んでいく。ifl7ファイルに内容をそのまま書き込んでいる。CATOMで呼ばれる場合は引数のiparaは0である。ifl7ファイルの内容をiwrt、即ちjfl3ファイルにそのまま書き込んでいる。ただし、各原子軌道に割り振られた電子数は新しい値としている。
  □ CONSNT 定数を設定するルーチンである。引数msは単位系のオプションで、msが3以下の場合Rydberg単位、4以上の場合au単位の値になる。結果はcommon文CNSTで引き渡される
  □ FILGEN いろいろなファイル番号を設定している。結果はcommon文IOCNTLで引き渡される。
    iegi: 分子軌道係数coeを保存しておくファイル
    itp: DVサンプルポイントに関するデータであるサンプルポイント密度、各原子との距離、mesh点を保存しておくファイル
    iyl: 使われていない
    isav: 使われていない
    ipun: 計算の途中で得られた原子の電荷密度を保存しておくファイル
    ifl1: 使われていない
    ifl2: 対称軌道に関するデータを読み込むファイル。
    ifl3: 原子の電子密度を読み込むファイル(初期データを作成する場合もブランクデータをreference原子の数だけ入れておく必要がある)。
    ifl4: 使われていない
    ifl5: 結果の出力ファイルである。一般にioutに等しい。
    ifl6: 結果の出力ファイルである。一般にioutに等しい。
    jfl1: 使われていない
    jfl2: Madelungポテンシャルに関する入力データのファイル
    jfl3: FREADにおいて入力データを書き込むファイル
    jfl4: 密度解析の結果である出力することにした軌道についての情報や軌道エネルギー、overlap population を保存しておくファイル
    jfl5: 使われていない
    jfl6: 使われていない
    jfl7: 使われていない
    jfl8: 使われていない
  □ INITIA k個の要素を持つ一次元配列に0.0を代入、即ち初期化するルーチン
  □ RNUCLR 原子の半径を計算するルーチン。引数のうちzは原子番号、rnucは原子半径、anucは核の質量数である。データbnucには各原子の核の質量数が格納されている。
  □ DINPT 各reference原子(ne)ごとに呼ばれる。引数のうちne、ms以外は全てここでの結果である。read文で読み込んでいる内容は以下の通りである。
    name: reference原子の名前
    loco: この原子の順番。大抵は順番に読み込んでいく。
    n: mesh数
    j: 原子の原子軌道の数
    ja: 
    zn: 原子の原子番号
    xalph: 原子のXa値
    rn: 最大mesh点距離
    h: mesh間の増分
    xion: 
    anuc: 原子の質量数
    xni: 原子軌道のn値(主量子数)
    xli: 原子軌道のl値(方位量子数) 
    xei: 原子軌道の固有値
    xzi: 原子軌道の占有電子数
    xz2i:
  □ RHOB 各reference原子(ne)ごとに呼ばれる。aがrh, iatcの結果、bがworking file、それ以外は入力データである。番号ifl3だるファイルに何か値が入っているかどうか最初の5点をreadしてみる。もし何か値が入っていた場合は、iatcは0となり残りのデータを読み込む(readする)。値が入っていなかった場合はiatcは1となり、このサブルーチンで初期データを作成する。電荷密度はρ(r)であるが、ここでは-4πr2ρ(r)の値を各mesh点上の値としてau単位で計算する。
  □ DIMCHK ここまでのデータをもとに次元のチェックを行い、入力データとして出力ファイルに書き込む。データのやり取りは全てcommon文を通して行われる。
  □ CORES reference原子の原子軌道についての変数を設定する。mstar、nco、indxが結果である。reference原子を順に並べて、その原子軌道について順に見ていったときのi番目のreference原子の原子軌道の始まる数がmstar(i)、i番目のreference原子のもつ原子軌道の数がnco(i)である。indx(i)はi番目の原子軌道の属するreference原子(何番目のreference原子か?)である。
  □ SYMTRY DV-Xa法では分子軌道は原子軌道(SCF計算毎に得られる電荷密度分布を用いてreference原子毎に計算した原子軌道)の線形結合で表されるが、計算及び定性的解釈の容易さから対称軌道の線形結合で表すことが多い。原子軌道及び対称軌道についての諸量を整えているのがこのサブルーチンである。引数のうちkkk2、nsym、naos、nmos、maos、misos、isyml、jsyml、lorb、mtitl、lmaxz、neq、nco、mstar、locor、indx、eignl、lval、loca、nbpot、kat、mcopt、mval、npart、mval1、ichi、kido、nchi、mcor、mn、natom、cnが結果である。対称軌道を用いない場合はこのサブルーチン内でASYMTYを読んでおり、そこで諸量を整える。このサブルーチンの詳細は以下の通りである。
  分子の軌道関数は分子軌道の1つのSlater行列式で表される。分子軌道は原子軌道の線形径都合で表される。原子軌道は分子を構成する原子の原子軌道とそのものを用いても良いが、対称性のある分子の場合、対称軌道を用いた方が重なり積分行列やFock行列がblock outあするので、計算時間の短縮及び定性的解釈の容易さが期待される。
  サブルーチンSYMTRYでは、対称軌道の読み込みと各種データの構築を行っている。また、対称軌道を用いない場合も、このサブルーチンの中で、サブルーチンASYMTYを呼び出し、データを整えている。
    mtitle: 対称ブロックの各名称及びコメントに関するデータ
    npat(i): 対称軌道を構成するbasisを線形結合の順に並べた際のbasisの番号
    eignl(i): 対称軌道を構成するbasisの規格化された線形結合を、順番に書き下したもの。npart(i)の順に規格化した値が入る。
    kido(i): 対称軌道を構成するbasisの(n * 100 + l)値
    kat(i): 対称軌道を構成するbasisの個数
    mval1(i): 対称軌道を構成するbasisの線形結合のうち最後のbasisのm値
    ichi(i): 対称軌道を構成するbasisのうち初頭のbasisの属する原子番号
    nchi(i): 対称軌道を構成するbasisのうち初頭のbasisの属するreference原子番号
    mcor(i): 対称軌道を構成するbasisの原子軌道番号
    loca(i): 対称軌道を構成するbasisの属する原子番号
    nbpt(i): lcoa(i) のbasisの属するreference原子番号
    mcopt(i): joca(i) のbasisの属するreference原子の原子軌道番号
    mval(i): lcoa(i)のbasisのm値
    lval(i): loca(i)のbasisのl値
    lmaxz(i): 分子を構成する各原子の最大値
    norb(i): 各対称ブロックを構成する対称軌道の数。入力のndimの値とは異なる場合がある。
    lorb(i): その対称ブロックの前のブロックまでに用いられた対称軌道の数
    配列になっていない下記の変数の意味は次の通りである。
    naos: basisの数
    nmos: 対称軌道の数(=naos)
    maos: npart(i)の要素の数、すなわち、basisの総数
    misos: 構成する対称軌道の数が最大のブロックのその対称軌道の数
    ksiz: 一次元化したときの行列要素の総数
    対称性を利用しない場合は、サブルーチン SYMTRYの中で、サブルーチンASYMTYが呼ばれ、そこで各データが構築される。SYMTRYで入力されるデータはnsymとisyml(i)、jsyml(i)であり、それぞれ、nsym=1、isyml(1)=0、jsyml(1)=0である。各変数は1つのbasisより成る対称軌道と考えられる。
  □ ASYMTY(SYMTRYより呼ばれる) 分子軌道がクラスターを構成する原子の原子軌道そのものの線形結合で表されるとした場合の原子軌道に関する諸量を整えている。引数のうちnaos、mm、ibs、ii、mval1、ichi、nchi、kido、mcor、mcopt、mval、lval、loca、nbpt、kat、indx、lmaxz、npart、eign1が結果である。
  ASYMTYが呼ばれる場合、SYMTRYで入力されるデータはnsym、isyml(i)、jsyml(i)についてそれぞれ、nsym=1、isyml(1)=0、jsyml(1)=0である。これは分子によらず対称性を用いない場合は常に同じである。各変数は1つのbasisより成る対称軌道であると考えると対称軌道を用いる場合との対称を付けやすい。
  □ ATMPOS クラスターを構成する原子のx, y, z座標及びDVサンプルポイントについての諸量を整えている。引数のうち結果は neqp、kplv、exalf、exafc、xfac、xvp、yvp、zvp、vechas、tick、alf、r0、beta、gamma、wltabである。read文がいつかあるが、それらの内容は以下の通りである。
    xvp, yvp, zvp: クラスターを構成する原子のx, y, z座標及び点電荷を加える場合のそのx, y, z座標
    neqp: その原子の種類
    exafl: 各原子の勢力範囲外でのXa値
    vecbas:  DVサンプル点をランダムに発生させるもとになる無理数
    tick: クラスターを構成する原子に割り当てるDVサンプル点の割合
    alf:
    r0: クラスターを構成する原子の勢力範囲の半径(原子半径程度の値 in au)
  □ MDINPT Madelungポテンシャルを加えた計算を行う場合に必須な各種データを読み込み、計算を行い、点電荷を置くx, y, z座標と点電荷の大きさについての情報を得るサブルーチンである。引数のうち結果は、axx、ayy、azz、chz、kmadであり、xx、yy、zz、neqv、xot、yot、zot、ezはこのサブルーチン内でread文で読み込まれる。読み込まれる変数の意味をまとめると以下の通りである。
    katoms: 単位格子内の原子の数(Madelungポテンシャルを発生させるための原子数)
    ntype: イオン種の数
    nout: defectとして除く原子数
    ez: イオンの価数
    ax, ay, az: 単位格子の格子軸の長さ(原子単位)格子定数
    xx, yy, zz: 格子軸の長さを単位とする格子内原子のx, y, z座標
    neqv: 格子内原子のイオン種の番号
    mx, my, mz: 計算上あらかじめ発生させるMadelungポテンシャルに含める点電荷の範囲(格子定数単位)
    xot, yot, zot: defectとして除く原子の座標(原子単位)
  このサブルーチンにおいて実際になされているのは、必要なデータのread文による読み込みと結果の書きだしであり、点電荷についての諸量の計算はLTTCEいて行われている。
  □ LTTCE (MDINTPにより呼ばれる) Madelungポテンシャルとして考える点電荷の位置や大きさを計算するルーチンである。引数のうち結果はkk、axx、azz、chz、kmadである。Madelungポテンシャルについての入力データは単位格子に関してであるので、ここでそれらをもとに全点電荷の位置と大きさを求める。
  □ POTM(PTGETより呼ばれる) Madelungポテンシャルの値を計算するルーチンである。あるDVサンプルポイントごとに呼ばれる。
  □ PTGET DVサンプルポイント全てのx, y, z座標と各点におけるクラスターを構成する原子の原子軌道の角度成分の値を計算し、結果をファイルに書き出しておく。引数の大半は入力データで、結果はkpl1、yl1、ylである。
    call DVPT DVサンプルポイントのx, y, z座標を求める。
    r2(i): あるDVサンプルポイントとi番目の原子との距離
    call KPLGET あるDVサンプルポイントと原子iとの間の距離r2(i)に相当するその原子が参照しているreference原子のmesh点( kpl1(i) )を求める。
    den: DVサンプルポイント密度である。原子に関して和を取って得られる。
    call YLMGET, HARM 原子軌道の角度部分、すなわち球面調和関数の値を求め順番に格納する。
  □ DVPT あるDVサンプルポイントのx, y, z座標を計算する。引数のうちxp, yp, zpが結果、すなわちnn番目のDVサンプルポイントのx, y, z座標である。
  □ KPLGET あるDVサンプルポイントとi番目の原子との距離(rnu)に相当するその原子が参照しているreference原子のmesh点(kpl1)を求める。kpl1が原子中心から近すぎたり(原子中心付近)、遠すぎたり(最遠mesh点距離)、原子半径程度の値(r0)付近であった場合は適宜補正を加える。
  □ YLMGET あるDVサンプルポイントにおけるi番目の原子の原子軌道の角度部分を求める。引数のうち結果はylである。DVサンプルポイントのx, y, z座標(xz, yz, zz)及び考えている原子との距離(rnu)から実際の計算に便利なように諸量を求め、ALM1を呼び出して球面調和関数を求めている。
  □ ALM1(YLMGETより呼ばれる) 実数型球面調和関数を逐次的に求める。ylが結果である。
  □ HARM(PTGETより呼ばれる) YLMGET及びALM1で計算された実数型球面調和関数は2次元配列の形で結果が得られているので、このルーチンで一次元配列に直している。yl1がその結果になる。
------------------------------------------------------------------------------
■ ATCALCから主に呼ばれる下位のサブルーチンの詳細
  □ ATNSCL reference原子の原子軌道計算を行って原子の電荷密度を求める。引数のうち、結果はrh、xz、xe、iido、xn、xl、rbar1、vbar1であり、a、c、d、y、 vrはworking fileである。
    jj: 原子を構成する原子軌道の数
    nc1: 原子軌道計算における繰り返しの上限
    delvr: 原子軌道計算におけるthreashold
    rbar1: 原子の井戸型ポテンシャルの半径
    vbar1: 原子の井戸型ポテンシャルの深さ
    ph: 繰り返しの際に取り込む前回の結果の割合
    phi: 繰り返しの際に取り込む前回の結果の割合を調整するパラメータ
    xn(i): 原子軌道のn値
    xl(i): 原子軌道のl値
    xe(i): 原子軌道の軌道エネルギー
    xz(i): 原子軌道の占有電子数
  □ CLPOT(ATNSCL、CATOM及びSSOGTより呼ばれる) 原子のクーロンポテンシャルをシンプソン則で求めるルーチンである。引数のうち結果はクーロンポテンシャルvhであれい、c、dはworking fileである。
    p4: 1/4π
  □ HFCAL(ATNSCLから呼ばれる) 原子の原子軌道の動径分布をSchrödinger方程式を数値的に解くことで得、原子の電荷密度を計算するルーチンである。引数のうち結果はrhである。
    pi4: 4π
    third: 1/3
  □ WELPOT(HFSCAL及びSSOWVより呼ばれる) r > r0 での井戸型ポテンシャルの計算を行うルーチンである。引数のうち結果はvrである。他のポテンシャルが計算上 rV(r)の形をしているのと同様に井戸型ポテンシャルも rV(r) の形で求める。
    タイプ1 井戸型ポテンシャルに加えて、r > r0 では -1/r に比例するポテンシャルにする
    タイプ2 井戸型ポテンシャルに加えて、r > r0 で直線的に大きくなるポテンシャルにする
  □ SINPT(HFSCAL及びSSOWVより呼ばれる) Schrödinger方程式の動径部分を数値的に求めていく際の初期データ、すなわちHamming法における積分間隔の設定とポテンシャル rV(r) の3字のべき級数展開における展開係数の計算を行うルーチンである。引数のうち結果は、eps、rbar、vbar、h3、ha1、ha2、ha3、ha4、ha5、ha6、ha7、ha8、ha9、fdll、voc、ndebuであるが、このうち実際には、eps、rbar、vbar、ha1、ha2、fdll、voc、ndebuしか今後の計算で用いられていない。
  □ HUMMG(SHEQより呼ばれる) Hamming法による微分方程式の解を求めるルーチンである。引数のうち結果はda、db、aである。このうちda、dbは次の初期値であり、aが求めたい波動関数である。求めようとしている方程式は2次微分方程式である。波動関数の2次微分量から1次微分量を求め、1次微分量から波動関数を求めるといった2段階の手法をとっている。変数のうちaで始まるものが1次微分量から波動関数を得る際に関係するもの、bで始まるものが2次微分量から1次微分量に関するものである。
    fmodc: 112/121
    fcorc: 9/121
  □ DIFSTR(SHEQより呼ばれる) Hamming法による微分方程式の数値解ルーチンにおいて必要となる初期値の計算を行うルーチンである。引数のうち結果はa、da、dbである。
  □ DIADL(SHEQより呼ばれる) 波動関数の規格化の計算をシンプソン則を用いて行うfunction文である。
  □ SHEQ(HFSCAL及びSSOWVより呼ばれる) Schrödinger方程式を数値的に解くルーチンである。引数のうち結果は波動関数の動径部分であるaとその軌道エネルギーであるeの2つである。計算は主として7つの部分より成る。
------------------------------------------------------------------------------
■ CATOMから主に呼ばれる下位のサブルーチンの詳細
  read文で読み込んでいる内容は以下の通りである。
    nscf: SCF法のタイプをしていするパラメータ(SCCFなど)及び収束の如何にかかわらず、繰り返しを打ち切るかどうかのパラメータ。wellで井戸型ポテンシャルを考慮する。1111またはaaaaで打ち切りと見なす。入力データの最後には必ず必要である。
    nscc: これといって使われていない
    mmod: これといって使われていない
    mfix: nfixでitape=1、ireeでitape=0
    mflz: nflzでntape=1、nfulでntape=0
    nmax: DVサンプルポイント数。0でディフォルト値 500が代入される。プログラムの初めにreadされた値より大きい場合は、初めの値になる。
    phm: 繰り返しの際に取り込む前回の結果の割合。1.0d-7より小さいとき、ディフォルト値になり、0.3が代入される。
    del: baseの計算をする際のthreasholdの元となるパラメータ。使われていない
    shikii: 原子の電界に関するしきい値
    jprt: 出力のオプション
    iwre: 出力のオプション
    ndwav: 出力のオプション
    ndvr: 出力のオプション
    ja: SUBWTRで使われている。結果の出力を何番目から始めるか。0でディフォルト値1が代入され、最初から出力を始めることになる。
    ls: SUBWTRで使われる。結果の出力をどこまで行うか。
  □ SSOGT DV-Xa法では分子軌道は原子軌道の線形結合で表されるLCAO近似を採用している。原子軌道、すなわち基底関数は解析的な形ではなく数値的な形で得られるが、それを求めるのが本サブルーチンである。最終的な結果はeco、vht、snloであり、common文でやり取りされる。
  □ BASCRT ある原子の原子軌道、すなわち基底関数を算出するルーチンである。引数のうち結果はvht、snlo、ecoであり、c、dはworking fileである。
    call SSOPOT 周囲の原子の影響を考慮した全ポテンシャルの算出を行う。ndvrが5より大きい場合、全ポテンシャルを出力する。
    call SSOWV 原子軌道の動径部分を求める。
  □ SSOPOT(BASCRTより呼ばれる) 周囲の原子の影響を考慮した全ポテンシャルを算出するルーチンである。引数のうち結果は全ポテンシャルvhtであるサブルーチンBASCRTは各reference原子毎に呼ばれるものであった。何番目のreference原子かを示す引数がndであり、naはndをreference原子をとする原子の中の最初の原子を示している。
    call SUMCH 周囲の原子の影響を考慮したポテンシャルを算出する
  □ SUMCH(SSOPOTより呼ばれる) マフィンテイン型のポテンシャル及び電荷密度を得るルーチンである。引数のうち結果はchsmである。入力のchatがポテンシャルか又は電荷密度であるので、結果もそれに従ってどちらかの値になる。   
  □ DINTX(SUMCHより呼ばれる) マフィンティン型積分を数値積分で行うルーチンである。引数のうち結果はsumである。rhが積分計算をする関数である。
  □ SSOWV(BASCRTより呼ばれる) 原子の波動関数の動径部分を求めるルーチンである。ATCALCブロックにおけるサブルーチンHFSCALとほぼ同じルーチンである。引数のうち結果はxx、aである。
    call SINPT Schrödinger方程式の数値解を求めるルーチンSHEQに必要な諸量を計算する。
    call SINPT 全ポテンシャルに井戸型ポテンシャルを加えた場合のSHEQの初期データを構築する
  □ FOCKM Fock行列及び重なり行列を求めるルーチンである。結果smx1、smx2、fmxはcommon文でやり取りされる。各DVサンプルポイントについてループを取り幾つかのサブルーチンを呼んでFock行列と重なり行列を算出している。このサブルーチン内でcallしているサブルーチンとその内容の概略は以下の通りである。
  ・CPOTG 全ポテンシャルについてそのDVサンプルポイントでの値を補間して得る
  ・CHARB 電荷密度についてそのDVサンプルポイントでの値を補間して得る
  ・XPOTG 交換ポテンシャルについてそのDVサンプルポイントでの値を補間して得る
  ・BASGET 波動関数についてそのDVサンプルポイントでの値を補間して得る
  ・SOGET 対称軌道を考慮したDVサンプルポイントでの波動関数を得る
  ・HSMATX Fock行列、重なり行例を足しあげることで求める
  ・DEGEN DVサンプルポイントの重みについて考慮する
  サブルーチンFOCKMで呼び出すcall については以下の通りである。
    call BASGET h(rk) Φv(rk)とΦv(rk)を求める
    call SOGET BASGETで求められた波動関数についての量はbaseに関するものであるので、ここで対称性も考慮したh(rkv(rk)とΦv(rk)を求める
    call HSMATX Fock行列要素、重なり行列要素を得る。ループの途中ではどちらも途中結果であり、ループが全て終了した時点でそれぞれ意味ある値になる。
  □ CPOTG(FOCKMより呼ばれる) あるDVサンプルポイントにおける全ポテンシャルの値を、補間により求めるルーチンである。引数のうち結果はpotである。全ポテンシャルはreference原子の各mesh点上の値として得られている。DVサンプルポイントはランダムに決められているので、当然mesh点上に来るとは限らない。そこで、そのサンプルポイントに最近のmesh点の値を用いて補間することでポイント上の値を得る。
  □ CHARB(FOCKMより呼ばれる) あるDVサンプルポイントにおける電荷密度の値を補間により求めるサブルーチンである。引数のうち結果はcharである。これまでで得られた電荷密度はreference原子についての各mesh点上での値であり、連続的に得られているのではなかった。サンプルポイント上での値は、そのmesh点上での値を用いて補間することで得られる。
  □ XPOTG(FOCKMより呼ばれる) あるDVサンプルポイントにおける交換ポテンシャルの値を求めるサブルーチンである。引数のうち結果はxpotである。交換ポテンシャルは-6a{3/8π * ρ(r)}^1/3の形をしている。ここでaはXa値であり、各原子によって異なった値を取りうる量である。原子の計算ではXa値は一義的に決まるが、分子では構成原子から均一に寄与を受けた値と考えて計算を進める。
  □ BASGET(FOCKMより呼ばれる) あるDVサンプルポイントにおける波動関数の値を補間により求めるサブルーチンである。引数のうち結果はbas1とtbas1である。Schrödinger方程式を数値的に解くルーチンSHEQで得られた波動関数は、各mesh点上での値としてであった。一方、角度部分である球面調和関数は、各DVサンプルポイント上での値が得られていた。このサブルーチンでは、DVサンプルポイント上で補間した動径部分と球面調和関数により、そのDVサンプルポイントでの原子軌道を求める。
  □ SOGET(FOCKMより呼ばれる) 対称軌道は原子軌道をbasisとして構築される。あるDVサンプルポイントでの原子軌道の値はBASGETで得られるので、ここではそれらを用いて対称軌道を作成し、軌道wav11(i)及び運動エネルギーが作用したtwav11(i)を得る。引数のうち結果はwav11, twav11である。
  □ HSMAX(FOCKMより呼ばれる) Fock行列及び重なり行列を求めるルーチンである。引数のうち結果はfmx1、smx1である。Fock行列及び重なり行列は実対称行列になるので、下三角部分のみ計算し、一次元化して配列に記録する。また、対称軌道を用いる場合は、各対称性に関して行列がブロックアウトする。そこで、その各ブロックの下三角部分を一次元化していく。対称性に関して縮退がある場合は、縮退度分足し込んで各行列要素が得られる。
  □ DEGEN(FOCKMより呼ばれる) DV数値積分は各サンプルポイントでの重みをかけて和をとるものであるが、サブルーチンHSMATXではサンプルポイント密度の逆数を掛けて行列要素を得ており、サンプルポイントの総数分大きい値になっている。HSMATXはサンプルポイントの数だけ呼ばれるので、この中でサンプルポイントに依存しない量の演算を行うのは、CPU的に無駄である。そこで、このルーチンで正しいDV数値積分値を求めている。HSMATXでは対称性に関して縮退がある場合はその縮退度分足し込んであった。ここでは、そこで、その分余計に割らなければならない。引数のうち結果は、smx1、fmx1である。a1は重み因子と言えるが、大体が1.0が代入される。
  □ INTERP(CPOTG、CHARB、BASGETより呼ばれる) Lagrange補間多項式による関数の補間を行うルーチンである。引数のうちrは距離、pは補間すべき関数、rsは補間すべき点、psは補間された結果である。
  □ EMERGE(CATOMより呼ばれる) これまででFock行列及び重なり行列が求められている。ここでは、分子軌道法の主たる計算部分である対角化を行っている。今得られている永年方程式はHC=ESCといった形をしており、一般的な対角化問題になっていない。また、各行列はブロック対角化されているので、各ブロックごとに対角化問題を解けばよいことになる。結果はcommon文を通してやり取りされている。結果は、eet、ititl、ntitl、coeである。
    call INDXG 下三角行列を一次元化する場合の区切りを示すパラメータであるndexを得ている
    call HAUSU 一般対角か問題への帰着及び対角化はサブルーチンHAUSUで行う。
    call ORDRNG 各ブロックごとに固有値が得られているが、ここでORDRNGを呼んで全体として小さい順に入れ替えを行っている(実際には入れ替えているのではなく順番をファイルにとっておく)。
  □ HAUSU(EMERGEより呼ばれる) HC=ESCといった形を持つ対角化問題をHouseholder法で解くサブルーチンである。結果は固有値がd、固有ベクトルがcoeである。ルーチンはいくつかのサブルーチンを呼んでいる。以下、内容を纏める。
  ・ HREDUO HC=ECといった一般対角化問題に帰着させる
  ・ HTRIDI 三重対角化を行う
  ・ UNITM 固有ベクトルに初期値を代入する
  ・ TQL2 QL分解を行って固有値、固有ベクトルを求める
  ・ HTRIBK 得られた固有ベクトルを2次元化する
  ・ HREBAK 得られた固有ベクトルを元の値に還元する
  □ ORDRNG(EMERGEより呼ばれる) 対角化により得られた固有値を小さい順に並べ替え、固有ベクトルをそれに対応させるルーチンである。引数のうち結果は、ititl、ntitl、eet、ispinである。入力データはブロックアウトした各ブロックの固有値smxwである。これを小さい順にして全固有値が入るべくeetに入れ直す。ititlにはどの対称ブロックに属しているかの情報が、ntitlには小さい順に並べ替えたときの固有値が各対称ブロックにおいて元々何番目にあったかについての情報が入れられる。ispinはスピンの別であるが、閉殻系の場合、常に1となる。また、対称軌道を用いない場合、ititlも常に1が入ることになる(対称軌道を用いないということは、1種類の対称ブロックであると考えると明白である)。縮退している時にはその縮退度の数だけ呼ばれ、固有値が整えられる。
  □ CHRANL(CATOMより呼ばれる) これまで得られた分子軌道からMullikenの密度解析を行って電荷の再分配、すなわち原子軌道に入る電子数の計算を行うルーチンである。結果はcommon文でやり取りされ、xam、clno、niiである。ss、rhoi、csijはworking filesである。
    call OCCUPS 電子を下から順に詰めていく。遷移状態の計算もここで行われる。
    call NIIG 軌道係数はORDRNGにおいて変化しないので、全固有値の小さい順ではなく、対角化されて得られた順になっている。ここで、その固有ベクトルが何番目のレベルのものであるかの対応づけを行っている。
  □ OCCUPS(CHRANLより呼ばれる) これまでで分子軌道及び軌道エネルギーが得られているので、ここでは、各軌道にレベルの低い方から順に電子を詰めていく。一般的閉殻系の場合、電子は下から2個ずつ入っていく。遷移状態の計算を行う場合は、電子数を操作して軌道に詰めていく。引数のうち結果はelnoであり、各軌道に入る電子数が得られる。
  □ NIIG(CHRANLより呼ばれる) 分子軌道の軌道係数、すなわち対角化で得られた固有ベクトルは得られた順にファイルに保存されている。一方、軌道エネルギーはORDRNGにおいて小さい順に並べ変えられている。ここでは、その軌道が何番目のレベルのものに相当するかをniiとして求めている。
    call CIJS ある分子軌道のoverlap populationを得る
    call ORBPOP orbital population のもとを得る
    call GRSPOP 原子軌道に電子を割り振り各軌道に属する電子数xzmを得る
  □ CIJS(MULCHS、SLVPOP及びDENMTXより呼ばれる) 二次元配列を一次元化した場合の掛け算。あるlのもので、c(i,l)*c(l,j)**(ij)を一次元化して得る。軌道lのoverlap populationである。nは行列cの大きさ、cr、srはそれぞれ一次元化された軌道係数と重なり行列、csijが一次元化された結果である。   
  □ ORBPOP(MULCHS、SLVPOPより呼ばれる) あるlのもとで、Σjc(i,l)*c(j,l)*s(i,j) = Σjcs(i,j) = ss(i)の計算を行っている。nは行列ssの大きさ、csijは一次元化されたcs(i,j)である。ssがこのルーチンの結果になる。
  □ GRSPOP(MULCHSより呼ばれる) 各原子軌道に電子の再分配を行う。mcorは、対称軌道を構成するbasisである原子軌道のうち初頭のbasisの属するreference原子軌道についてのデータである。それが同じ物について和を取っていく。
  □ JUDGES(CATOMより呼ばれる) 分子軌道計算が収束したかどうかを判定するルーチンである。引数のうち結果はlitern、xzxである。収束している場合liternは1、していない場合は0の値が代入される。
  □ CHMIXS(CATOMより呼ばれる) 繰り返しの次の初期データを作成するルーチンである。引数のうち結果はrhoである。これまでで、分子軌道とMullikenの密度解説により再分配された、原子軌道の電子の占有数が得られている。こから新たに電荷密度を計算して求め、混合比を考慮して繰り返しの次の初期データとしている。

------------------------------------------------------------------------------
■ POPANLから主に呼ばれる下位のサブルーチンの詳細
  □ CNUM 密度行列の出力に要する諸量の設定
  □ SLVPOP orbital populationの計算
  □ LVLPOP 結果の出力
  □ DENMTX 密度行列の計算
  □ KAKZO 密度解析の結果の出力

  □ CNUM 各対称ブロックごとに呼ばれ、密度解析の結果の出力に必要な諸量を整える。引数のうち結果はkct、numb、num、knum、kkである。
  □ SLVPOP 各対称ブロックごとに呼ばれ、ここで、その対称ブロックの軌道のoverlap population の計算を行っている(分子軌道係数から得られる量であり、まだ軌道の占有数は作用されていない)。引数のうち結果はdである。また、結果は軌道エネルギーと共にファイル番号jfl4に保存される。
  □ LVLPRT 各対称ブロックごとに呼ばれ、軌道エネルギー、軌道のoverlap populationについて出力を行うルーチンである。クラスターが大きくなるにつれ結果も増加するので、CNUMのインデックスをもとに、出力の制限を行っている。
  □ DENMTX 各対称ブロックごとに呼ばれる。あるレベルを構成する各軌道に電子を振り分けたときの値、すなわちorbital populationであるeiejを得る。
  □ KAKZO 各対称ブロックごとに呼ばれる。niは出力の桁数を指定するパラメーターであり、実際の出力はni-1桁になる。
  □ UNITM(HAUSUより呼ばれる) ある行列を単位行列とするルーチンである(対角項は1.0、非対角項は0.0をい代入する)。
  □ ACOPY(EMERGE、FOCKMより呼ばれる) 引数のうち結果はf2であり、f1をf2にコピーするルーチンである。
  □ RHSAVE(ATCALC、CATOMより呼ばれる) 原子の電荷密度について得られた結果をファイルに書き出しておくルーチンである。4πr2ρ(r)の形にしてファイル番号ipunに書き込んでいる。読み込むときは通常のread文で行われる。
  □ SNMES(DINTXより呼ばれる) シンプソン則による積分計算のfunction文である。
------------------------------------------------------------------------------
重要な文献&HP
[1] はじめての電子状態計算ホームページ: http://chem.sci.hyogo-u.ac.jp/hajimete/
[2] DV-Xa法に関するブログ: http://chem.sci.hyogo-u.ac.jp/dvxa_blog/ 
[3] ベータ版ダウンロードページ: http://www.dvxa.org/ 
[4] 坂根弦太・DV-Xa法計算支援環境を構築するためのリンク集01: http://www.chem.ous.ac.jp/~gsakane/dvxa_assistance_environment_1.html から 「DV-Xα法計算支援環境利用の手引き(pdf文書, 2010/11/27; 146頁; 11.0 MB)」をダウンロードする
[5] DV-Xaホームページ: http://www.dvxa.org/
[6] 楽しい電子状態計算: http://www.chem.ous.ac.jp/~gsakane/fun/
[7] 岩沢美佐子ら、DV‐Xα法による電子状態計算―そのプログラムと解説、三共出版
[8] 『電子構造論による化学の探究』
[9] 第一原理計算 〜構造最適化に向けた 材料・デバイス別 事例集〜、情報機構
[10] Slater Transition Stateの計算やイオン化エネルギーの計算、1電子励起状態の計算などを行う際の電子の移動方法: http://www.chem.ous.ac.jp/~gsakane/STS.html
[TX1] H. Adachi et al., "Hartree-Fock-Slater Method for Materials Science", Springer, 2005., google books; http://books.google.co.jp/books?id=A0-sWjTpi0gC&pg=PA338&lpg=PA338&dq=Auger+CVV+DV-Xa&source=bl&ots=yUeiIwVzgl&sig=bXu1vKtuyEPNnSKB_vxlXVf70Sw&hl=ja&sa=X&ei=dwEiUaqAGYbymAX2-YGICA&ved=0CDMQ6AEwAQ#v=onepage&q=Auger%20CVV%20DV-Xa&f=false
------------------------------------------------------------------------------
DV-Xa研究協会会報
[TSX1] 足立裕彦、DVXa研究協会会報、8 (1995) 94.: ピーク強度に関しては、ScofieldがHartree-Fock-Slater modelで計算した原子のionization cross-sectionなどを用いて、それにatomic orbital populationをかけたものの和を取って求めている。この論文ではDV-Xa法を用いてイオン化確率を求めるプログラムを開発している。(FeO68-は非常に良い一致を示している)
[TSX2] 宇田英一郎ら、DVXa研究協会会報、5 (1992) 25.: 硫黄化合物の蛍光X線(Kβ線)スペクトルの計算が行われている。その結果では遷移状態計算よるよりも通常の基底状態の計算が良い一致を示す結果となった。これは「S6のような価電子帯が込み入った多数の分子軌道からなる系で内殻と外殻に0.5個ずつ空孔を配置した遷移状態のSCF計算をすると空孔の影響で価電子帯中の分子軌道の順位が入れ替わり空孔が別の分子軌道に移動してしまうことが多い。(中略)目的とする分子軌道に空孔を置いて計算を収束させたかどうか疑わしいというところにある」としている。計算値のブロードニングは、ローレンツ曲線で行い、半値幅は気体で1 eV, 分子性結晶の固体で1.5 eV, イオン性結晶の固体で2 eVとしている。
[TSX3] 関根理香、DVXa研究協会会報、9 (1996) 41.:
[TSX4] 向山毅、DVXa研究協会会報、8 (1995) 310.: COでのO 1s 及び C 1s空孔に対するK X線の放出率を、Frozen-orbital(FR)近似、遷移状態(Transition state, TR)法、そして、始状態では内殻に正孔があり、終状態では正孔が外殻に移った後のポテンシャルで計算した波動関数により遷移確率を計算したRelaxed-orbital approximation(RX)法で比較している。これらの計算結果の値は異なっているが、相対放出率はどのような近似で計算してもあまり差のないことが示されている。
[TSX5] 神田英之ら、DVXa研究協会会報、10 (1997) 41.: NiOの酸素K殻ELNESの計算を行っている。[Ni19O44]50-クラスターで、遷移状態法を用いて、酸素2p非占有軌道の部分状態密度(PDOS)と実験結果と比較している。ピークの数と相対的な位置は、実験と比較的良い一致を示す。ピーク強度比は525 eV近傍のピークのみ大きい。
[UVX1] 丸山睦弘ら、DVXa研究協会会報、8 (1995) 229.: M(bpy)2(CN)(M=Fe, Ru, Os)錯体の電子状態を計算している。遷移状態法を用いて紫外可視吸収スペクトルを計算。実験結果と比較的良い一致。
[12] 脇田久伸ら、DVXa研究協会会報、6 (1993) 87.: 窒素を配位原子とする銅(U)錯体のXANESスペクトルを計算。モデルは結晶の構造解析結果をもとに構築。遷移状態法を用いたのかなどの詳細は未記載。(私の私見では、ブロードニング幅を調節すれば実験値と比較的良い一致を示しているように思える)
[13] 田中功ら、DVXa研究協会会報、7 (1994) 59.:(Mg13O14)2-と(Si5O16)12-クラスターでELNESを計算。遷移エネルギーの絶対値をSlaterの遷移状態法により求めた。遷移エネルギーの実験値との誤差は0.5%であった。(ピーク強度は合わないがピークの帰属は可能に見える)
[14] 松尾修司ら、DVXa研究協会会報、7 (1994) 145.: 窒素を配位原子とする銅(U)錯体のXANESスペクトルを計算。原子間距離をXANESスペクトルからDV-Xa法を用いて見積もっている。
[15] 吉矢真人ら、DVXa研究協会会報、8 (1995) 102.: a-Quartz([Si5O15]12-)のX線吸収スペクトルの計算を行っている。ELNESスペクトルのSi L2,3 edgeとO K edgeの実験値と、Si-3s/3dとO-2pの基底状態と遷移状態での部分状態密度(PDOS)との比較を行っている。また、光学遷移確率(PACS{Photoabsorption Cross Section})をそれぞれ基底状態と遷移状態について計算している。結果は、Si L23 edgeとの比較において基底状態と遷移状態のPDOSには大きな違いが認められない(どちらかというと基底状態のPDOSの方が実験値によく対応している)。一方、O K edge は遷移状態でのPDOSの方が一致が良い。PDOSとPACSの差は非常に小さい。
[16] 川崎晋司ら、DVXa研究協会会報、8 (1995) 237.: C60F46クラスターを用いて、遷移状態計算法を用いて1つのC 1s に対して計算を行い、C 2p 非占有準位と実験値と比較。比較的良い一致。
[17] 溝口照康ら、DVXa研究協会会報、13 (2000) 175.: 吸収断面積(PACS)を計算。(この他にも、この年の会報では参考になる吸収スペクトルの報告がある)
[18] 四元健司ら、DVXa研究協会会報、13 (2000) 178.: sxs計算によりFe K edge XANESのスペクトル計算をしている。
[19] 松尾修司ら、DVXa研究協会会報、16 (2003) 182.: sxs計算によりBr L3 edge XANESのスペクトル計算をしている。ブロードニングの半値幅はNaBrで1.5 eV、KBrで1.25 eVとしている。L2,3吸収近傍スペクトルの形状に目立った違いが見られなければ、非相対論のコードでもスペクトル解析に有効であることを示している。
[20] 越野雅至ら、DVXa研究協会会報、17 (2004) 60.:  F4TCNQスペクトルをDV-Xa法で再現するには、1s準位のエネルギーシフトが大きいため、core-hole効果と遷移確率を取り入れた遷移状態計算が不可欠であるとしている。
[MEX1] 宮前亨ら、DVXa研究協会会報、15 (2002) 109.
[MEX2] 渡邊真太、DVXa研究協会会報、21 (2008) 18.
[AX1] 田中康三ら、DVXa研究協会会報、5 (1992) 17.: Zn4Oクラスターに対するAugerスペクトルの計算。遷移状態法によって酸素の部分状態密度を求め、さらにそれのSelf-convolutionによりAuger電子スペクトルを計算している。計算結果は実験結果と比較的良い一致を示している。
[AX2] T. Yamamoto and M. Uda, DVXa研究協会会報、7 (1994) 154.:
[AX3] 山本知之、DVXa研究協会会報、12 (1999) 2.: http://www.process.mtl.kyoto-u.ac.jp/DV-XaPDF/DVXa12-2-1999.pdf
[AX4] 茂木昌都ら、DVXa研究協会会報、16 (2003) 260.: Al化合物のKLVオージェスペクトル
[APX1] 石井知彦、DVXa研究協会会報、17 (2004) 335.
[QX1] DVXa研究協会会報、7 (1994):
[QX2] DVXa研究協会会報、7 (1994) 269:
[DX00] 水野正隆、DVXa研究協会会報、21 (2008) 13:
[DX0] 水野正隆、DVXa研究協会会報、20 (2007) 242:
[DX1] 足立裕彦、DVXa研究協会会報、9 (1996) 180:
[DX2] DVXa研究協会会報、9 (1996) 112:
[DX3]DVXa研究協会会報、16 (2003) 84:
[TX1] 森永正彦、DVXa研究協会会報、6 (1993) 167.:
[TX2] 早藤貴範、DVXa研究協会会報、7 (1994) 134.:
[XCX1] 岩沢美佐子ら、DVXa研究協会会報、7 (1994) 197.:
[XCX2] 小野倫也ら、DVXa研究協会会報、10 (1997) 40.: LDAに拡張する方法がよく分かる。GGAは実空間差分法が良いとの結果。
[MX1] 野村重孝ら、DVXa研究協会会報、16 (2003) 106.:
[DMEX1] 吉田尚史、小笠原一禎、DVXa研究協会会報、20 (2007) 235.:
------------------------------------------------------------------------------

(下記は編集中)
-------------------------------------------------------------------------------
VESTAで出力可能なXYZファイルから SCAT F01のひな形を作成する Fortran90 プログラム

下記はVESTAで出力可能なXYZファイルからSACAT f01を作るプログラムです。
 
xyz2f01.f 90のコンパイル
1. テキストファイルに下記のプログラムをコピー&ペーストし、xyz2f01.f 90という名称でSAVEします。
2. Terminal を開いて、cd で xyz2f01.f90   のあるディレクトリをカレントディレクトにします。
3. g95 xyz2f01.f90  -o xyz2f01  などとしてコンパイルします(コンパイラによっては、g95 が gfortran や ifort になる)。そうすると、実行ファイルが作成されます。

使い方
1. VESTAでxyz形式で出力します。
2. xyz2f01.exe と VESTAでxyz形式で出力させたファイル を同じフォルダに入れます。
3. コマンドプロンプト上で、xyz2f01.exe VESTAでxyz形式で出力させた名称 を入力します。
4. f01のひな形が作られます(マーデルングの指定はご自身でして下さい)。

xyz2f01.f90
-------------------------------------------------------------------------------
      program xyz2f01
*     XMol XYZ output by VESTA --> SCAT F01 (template)

      character file1*512, line1*100, line2*100, symbol*2, elem(111)*2
     
data (elem(j), j = 1, 111)/ ' H','He','Li','Be',' B',' C',' N',
     &' O',' F','Ne','Na','Mg','Al','Si',' P',' S','Cl','Ar',' K','Ca',
     &'Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As',
     &'Se','Br','Kr','Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd',
     &'Ag','Cd','In','Sn','Sb','Te',' I','Xe','Cs','Ba','La','Ce','Pr',
     &'Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf',
     &'Ta',' W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At',
     &'Rn','Fr','Ra','Ac','Th','Pa',' U','Np','Pu','Am','Cm','Bk','Cf',
     &'Es','Fm','Md','No','Lr','Rf','Db','Sg','Bh','Hs','Mt','Ds','Rg'/

      open (unit  = 2, file =  'F01', form = 'formatted',
     &status = 'unknown', access =  'sequential')
      write(2, '(a)') '| Z ||NEQ||   X    ||   Y    ||   Z    |'

      call getarg(1, file1)
      open (unit = 1, file = file1, form = 'formatted', status = 'old',
     &  access = 'sequential')
      read(1, *) natom
      read(1,'()') ! Skip a line
      do i = 1, natom
         read(1, '(a)') line1
         do j = 1, 100
            if (line1(j:j) .ne. ' ') exit
         end do
         ibegin = j
         do j = ibegin + 1, 100
            if (line1(j:j) .eq. ' ') exit
         end do
         iend = j - 1
         len_symbol = iend - ibegin + 1 ! length of the chemical symbol
         select case (len_symbol)
         case (1)
            line2 = ''' '//line1(ibegin:iend)//''''//
     &      trim(line1(iend+1:))
         case (2)
            line2 = ''''//line1(ibegin:iend)//''''//
     &      trim(line1(iend+1:))
         case default
            write(6,'(a)') ' Too long symbol: '//line1(ibegin:iend)
            stop
         end select
        
         read(line2,*) symbol, x, y, z
         do j = 1, 111
            if (symbol .eq. elem(j)) go to 1
         end do
         write(6, '(2a)') ' A bad chemical symbol in *.xyz: ', symbol
         stop
    1    write(2, '(2i5, 3(f10.5))') j, i, x, y, z
      end do
      close(unit = 1)

      write(2,'(a)') '---------------------------------------------'
      write(2,'(a)') '|NEQ||  CHG   ||U/D||   RD   ||   VD   |    1'
      write(2,'(a)') '---------------------------------------------'
      write(2,'(a)') '    0     Unit     (0: Angstrom  1: Atomic)'
      write(2,'(a)') '    0     Spin     (0: Non-spin  1: Spin  )'
      write(2,'(a)') '    0     M.P.     (0: No        1: Yes   )'
      write(2,'(a)') '    0     Sample point (<100000, = 0 Autoset)'
      close(unit = 2)
      end
-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
f03ひな形作成 Fortran77 プログラム
-------------------------------------------------------------------------------
下記はf01からf03のひな形を作るプログラムです(編集中)
A. makef03.bat の作成
1. テキストフォルダを開いて 下記のmakef03.bat をコピー&ペースとして、makef03.bat という名称でSAVEする。
2. 作成した makef03.bat を exec ファイルへ入れる。

B. makef03.f のコンパイル
1. テキストファイルに下記のプログラムをコピー&ペーストし、makef03.f という名称でSAVEします。
2. 各行の位置はFortran77の形式に揃えてください。
3. Terminal を開いて、cd で makef03.f  のあるディレクトリをカレントディレクトにします。
4. gfortran makef03.f  -o makef03 と入力します。そうすると、実行ファイルが作成されます。
5. 作成したmakef03.exe をobject ファイルへ入れる。

使い方
1. f01に通常と同じようにマーデルングポテンシャルの記述をします。
2. f01に下記のオプション記述を追加します。
(編集中)
3. makef03と入力すると、f03が形成されます。

makef05をベースに作成(エンベッド法を念頭において作っていたので、Embedded Method とありますが、気にしないで下さい)

f01の下の部分「Deviation of Length from Center  :  1 %」から下記
のように記入します(コピーアンドペーストで対処できます)。

記入例:f01

| Z ||NEQ||   X    ||   Y    ||   Z    |
   13    1   0.00000   0.00000   0.00000
   26    2   1.44025   1.44025   1.44025
------------------------------------------------------------------
|NEQ||  CHG   ||U/D||   RD   ||   VD   |    1
   1        -3.0
   2         3.0
------------------------------------------------------------------
    0     Unit     (0:angstrom  1:atomic)
    0     Spin     (0:non-spin  1:spin  )
    1     M.P.     (0:No        1:Yes   )
    0     Sample Point ( <100000, =0 autoset )

Deviation of Coordination Number :  1 %
Deviation of Length from Center  :  1 %

------------------------------------------------------------------
   Embedded Method   |      XD     ||      YD     ||      ZD     |
LATTICE CONSTANT    | 200.0000000000 200.0000000000   4.0000000000 (real   /angstrom)
TRANSITION EXPANSION|              0              0            999 (integer/LATTICE CONSTANT)
SHIFT QUANTITY      |   0.0000000000   0.0000000000   0.0000000000 (real   /LATTICE CONSTANT)
CLACULATION RANGE ±|   0.5000000000   0.5000000000 999.0000000000 (real   /LATTICE CONSTANT)
------------------------------------------------------------------

「LATTICE CONSTANT」の行は格子定数をÅ単位で記入してください。
オプションは、*の後に入れます。
1:x軸ならx軸の座標で最も数値の大きい値をもつ原子を省く。y,zも同様。
そのため、最外部の原子が省かれ重なりを気にしなくてよくなる。ただし、斜めの構造を持つ結晶には対応していない。
2:手動で重なる部分を省いて、それに対応する格子定数を入れる場合に用いる。
以下は、格子定数を単位として何倍の範囲を考慮したいかを記入していきます。
「TRANSITION EXPANSION」の行は並進拡張(より広げる数)の記入場所です。
「SHIFT QUANTITY」の行は原点からどれだけ電荷を移動させるかを入力します。
「CLACULATION RANGE ±」の行は±の範囲で実際に計算する範囲を指定します。

*注1:LATTICE CONSTANTの値は、1次元周期の場合、他の軸で電荷の影響が殆んど無いように距離を多く記入してください。2次元の場合も同様です。

makef03.bat
------------------------------------------------------------------------------
@echo off
if not "%1"=="" cd %1
if not exist f01 goto quit
%dvdir%\object\makef03.exe
goto end
:quit
echo f01 not exist
:end
------------------------------------------------------------------------------


makef03.f プログラム
------------------------------------------------------------------------------
c=======================================================================
c  makef03.f (from makef05.f)
c   make f03 from f01 
c   2007. 9.26  Version 1.0   written by ******
c   2007. 9.26  Version 1.1   auto lattice parameter option
c   *計算機センター等で処理を行なわせる場合は、この日本語の文は省いて
c   ください。*を先頭につけておきますので、検索して取り除いてください。
c=======================================================================
      IMPLICIT REAL*8(A-H,O-Z)
      character*6 chstr,dum
      character*4 chstr2
      character*2 dum2,chstr3
      character*1 uord(100)
      DIMENSION NEQ(1000),Z(1000),NAME(10)
     & ,XN(100,100),XL(100,100),XE(100,100),XZ(100,100)
     & ,KCORE(100),ISYML(21),JSYML(21),IAT(100),JJ(100)
     & ,RR(10),TIC(2100)
      dimension xv(1000),yv(1000),zv(1000),nz(1000),ch(1000)
      dimension rda(200),vda(200)
      dimension xz2(100,100)
      dimension nzs(200)
c ---------------------------------------------------------------------
c additon 2007.9.26
c =====================================================================
      real*8 lcx,lcy,lcz
      integer*4 ntex,ntey,ntez
      real*8 sqx,sqy,sqz
      real*8 crx,cry,crz
      integer*4 zero
      real*8 lcux,lcuy,lcuz
      real*8 xvut,yvut,zvut
      character*1 odlc
      real*8 mxv,myv,mzv
      integer*4 nc
      integer*4 ATOMS

      DATA zero/0/
      DATA mxv,myv,mzv/0,0,0/
      DATA nc/0/
c =====================================================================
c  *新たに追加したデータの定義です
c ---------------------------------------------------------------------

c ---------------------------------------------------------------------
c  chstr= DATA chstr/'------'/
c  dum=*f01の読み取り
c  chstr2= DATA chstr2/'Nsym'/
c  du2=
c  chstr3= DATA chstr3/'XD'/
c  uord=

c  NEQ=*クラスターを構成する原子が何番目の原子にあたるか
c  Z=*クラスターを構成する原子の原子番号
c  NAME=*対称ブロックの名前
c  XN=*原子軌道のn値
c  XL=*原子軌道のl値
c  XE=*原子軌道の軌道エネルギー
c  XZ=*原子軌道の電子(の占有)数
c  KCORE=*その軌道について出力するかどうかの情報が入れられている。
c   joutが通し番号で入っている。
c   jout=原子軌道に関するデータ。一般に0である。
c       =1 原子軌道のbaseから除く
c       =2 最後の出力onlyのルーチンで使われる。
c  ISYML=*対称ブロックを扱うかのパラメータ、通常1
c  JSYML=*対称ブロックの縮重度 1:縮重なし 2:2重縮重 3:3重縮重
c  IAT=
c  JJ=*原子を構成する原子軌道の数
c  RR= DATA RR/1.0,2.0,2.0,2.5,2.5,3.0,3.0,3.0,3.0,3.0/
c  TIC=

c  xv,yv,zv=原子のx,y,z座標
c  nz=原子番号
c  ch=
c  rda=
c  vda=
c  xz2=
c  nzs=
c --------------------------------------------------------------------
      DATA RR/1.0,2.0,2.0,2.5,2.5,3.0,3.0,3.0,3.0,3.0/
      DATA chstr/'------'/
      DATA chstr2/'Nsym'/
      DATA chstr3/'XD'/
      IOUT=8
      INPUT=3
c      OPEN(3,FILE='nonrel')
c      OPEN(4,FILE='radii')
      OPEN(5,FILE='f01')
c      OPEN(8,FILE='f05')
c      OPEN(25,FILE='f25')
c ---------------------------------------------------------------------
      OPEN(13,FILE='f03')
c ---------------------------------------------------------------------
      NCYCLE=10
      N999=99999
      NTABPT=300
      NMAXK=5000
      NNST=0
      LT=999
c ---------------------------------------------------------------------
c  NCYCLE=*SCF軌道計算をする際の繰り返しの上限
c  N999=*
c  NTABPT=*原子軌道計算をする際のメッシュ数。
c   0でデフォルト値の選択になり、300が入る
c  NMAXK=*DVサンプルポイントの数
c  NNST=*
c  LT=*
c ---------------------------------------------------------------------
      ATOA=0.529177000000D0
      UT=1/ATOA
      NCORE=0
      NDAT=0
      TL=0
      extch=0.0
      nspat=500
c ---------------------------------------------------------------------
c  ATOA=*Åから原子単位に変換する場合に割る値
c  UT=1/ATOA *こちらを掛けると計算が速くなる(演算は掛け算の方が速い)
c  NCORE=*reference原子の原子軌道の総数
c  NDAT=*クラスターを構成する原子のうち原子軌道計算をする原子、即ち
c   reference原子の総数
c  TL=*全電子数
c  extch=*
c  nspat=*
c ---------------------------------------------------------------------

      MDBL=0
c ---------------------------------------------------------------------
c MDBL=*f01のx,y,z座標がdoubleタイプであるかを判定。0:普通 1:double型
c ---------------------------------------------------------------------
c =====================================================================
      read(5,740) dum2
      if (dum2.eq.chstr3) MDBL=1
      rewind 5
  740 format(17x,a2)
c ---------------------------------------------------------------------
c *f05の一行目を読み込んで、doubleタイプであれば MDBL=1 とする。
c =====================================================================

      jsep=0
      do i=1,1000
      read(5,750) dum
       if (dum.eq.chstr) then
        if (jsep.eq.0) then
         isep1=i
         jsep=jsep+1
         NAT=i-2
         NATOMS=NAT
c   -------------------------------------------------------------------
c   *dumでf01ファイルの内容を文字として読み込んでいって、
c   行を---が見つかるまで数えて上下2つを引いて原子の個数とする
c   -------------------------------------------------------------------
        else if (jsep.eq.1) then
         rewind 5
         isep2=i
         goto 70
        end if
       end if
      end do
  750 format(a6)
c ---------------------------------------------------------------------
c  jsep=*原子の位置座標(0)のデータと電荷のデータ(1)の切り替えに用いる
c  isep1=*---から---までを含めた原子の位置座標の行数を数える
c  NAT=*ポテンシャル用原子も含めたクラスターを構成する原子の個数
c  NATOMS=*ポテンシャル用原子も含めたクラスターを構成する原子の個数
c  isep2=*原子の位置座標からデータと電荷のデータの---まで行数を数える
c ---------------------------------------------------------------------

   70 continue

      KSP=0
      REMV=0.0
      RD=2.5
      VD=-2.0
      IDOA=-1
      TOLM=0.1D-05
      JPRT=-3
c --------------------------------------------------------------------
c  KSP=*
c  REMV=*
c  RD=*
c  VD=*
c  IDOA=*井戸型ポテンシャルを加えるかのパラメータ
c   0:加えない 1:タイプ1 2:タイプ2
c  TOLM=*分子軌道の繰り返し計算におけるthreshold
c  JPRT=*出力のオプション
c --------------------------------------------------------------------
    
      read(5,*)
      do i=1,NAT
       if (MDBL.eq.0) then
        read(5,710) nz(i),NEQ(i),xv(i),yv(i),zv(i)
       else
        read(5,715) nz(i),NEQ(i),xv(i),yv(i),zv(i)
       end if
       nzs(neq(i))=nz(i)
       Z(i)=real(nz(i))
       if (NEQ(i).gt.NDAT) then
        NDAT=NEQ(i)
        NTYPES=NDAT
       end if
      end do
  710 format(2i5,3f10.5)
  715 format(2i5,3f15.10)
c ----------------------------------------------------------------------
c  NTYPES=*ポテンシャル用原子も含めたクラスターを構成する原子のうちの
c   reference原子の和
c ----------------------------------------------------------------------
     
c----- Error check in f01 ----------------------------------------------
      do i=1,NAT
       if (nzs(neq(i)).ne.nz(i)) then
        write(*,713) neq(i)
c        goto 999
       end if
      end do
  713 format('Error -- duplicated Z in NEQ(',i2,')')

      do i=1,NDAT
       do j=1,NAT
        if (neq(j).eq.i) goto 80
       end do
        write(*,714) i
c        goto 999
   80  continue
      end do
  714 format('Error -- missing NEQ(',i2,')')
c-----------------------------------------------------------------------

      do i=1,NDAT
       ch(i)=0
       vda(i)=0
       rda(i)=0
      end do

      read(5,*)
      read(5,711) IDOA
  711 format(40x,i5)

      do i=1,isep2-isep1-2
       read(5,712) ii,ch(ii),uord(ii),rdi,vdi
       if (ii.le.NDAT) then
        if (rdi.ne.0) rda(ii)=rdi
        if (vdi.ne.0) vda(ii)=vdi
        if (ii.eq.0.and.ch(ii).ne.0) then
         extch=extch+ch(ii)
        end if
       end if
      end do
      IF (IDOA.eq.-1) IDOA=0
  712 format(i5,f10.5,2x,a,2x,2f10.5)
c -----------------------------------------------------------------------
c  *|NEQ||  CHG   ||U/D||   RD   ||   VD   |    1 の欄の読み取り
c  井戸型ポテンシャルの1タイプを加える
c -----------------------------------------------------------------------

       icalt=1
      do i=2,NAT
       if (nz(i).ne.nz(1)) icalt=2
      end do
c -----------------------------------------------------------------------
c  icalt=*
c  icalt2=*
c -----------------------------------------------------------------------
c  *
c -----------------------------------------------------------------------
     
      do i=1,NDAT
       if (rda(i).eq.0.or.vda(i).eq.0) then
c        call setwell(rdi,vdi,IDOA,nzs(i),ch(i),icalt)
        rda(i)=rdi
        vda(i)=vdi
        if (rda(i).eq.0) rda(i)=rd
        if (vda(i).eq.0) vda(i)=vd
       end if
      end do
c -----------------------------------------------------------------------
c  *
c -----------------------------------------------------------------------

      read(5,*)
      read(5,700) jut
       if (jut.eq.1) UT=1.0
      read(5,700) ISN
      read(5,700) MQ
  700 format(i5)
c -----------------------------------------------------------------------
c  *UnitからM.P.までのオプションの設定を読み込む
c -----------------------------------------------------------------------

      read(5,720) NMAXK
       if (NMAXK.eq.0) NMAXK=nspat*NAT
       if (NMAXK.gt.99999) then
        NNST=log10(real(NMAXK))-4
        NMAXK=NMAXK/10**NNST
       end if
  720 format(i5)
c     read(5,730) icalc
c      if (icalc.eq.0.or.icalc.gt.10) icalc=10
c 730 format(i5)
      icalc=1
c -----------------------------------------------------------------------
c  *DVサンプルポイントの数(NMASK)でSample Pointのオプションを読み込む
c -----------------------------------------------------------------------

c -----------------------------------------------------------------------
c additon 2007.9.26
c =======================================================================
      read(5,721) dum
      read(5,721) dum
      read(5,721) dum
      read(5,721) dum
      read(5,721) dum
      read(5,721) dum
      read(5,731) odlc,lcx,lcy,lcz
      read(5,741) ntex,ntey,ntez
      read(5,751) sqx,sqy,sqz
      read(5,751) crx,cry,crz
     
  721 format(A6)
  731 format(18x,A1,3x,3F15.10)
  741 format(22x,3I15)
  751 format(22x,3F15.10)
     
c =======================================================================
c  *追加のデータを読み込ませる
c -----------------------------------------------------------------------
     
c      do i=1,NAT
c       TL=TL+Z(i)-ch(neq(i))
c      end do
c       TL=TL-extch
c -----------------------------------------------------------------------
c  *
c -----------------------------------------------------------------------
c -----------------------------------------------------------------------
c additon 2007.9.26
c =======================================================================
c -----[pre1]
     
      do i=1,NATOMS
       if(odlc.eq.'1')then
        if(xv(i).ge.mxv)then
         mxv=xv(i)
        endif
        if(yv(i).ge.myv)then
         myv=yv(i)
        endif
        if(zv(i).ge.mzv)then
         mzv=zv(i)
        endif
       endif
      end do
c ----------
c  *最外にある原子の座標を探す
c ----------
c -----[pre2]
      do i=1,NATOMS
       if(odlc.eq.'1')then
        if(xv(i).eq.mxv .or. yv(i).eq.myv .or. zv(i).eq.mzv)then
         nc=nc+1
        endif
       endif
      end do
c ----------
c  *最外にある原子の数をカウント
c ------[pre3]
      if(odlc.eq.'1')then
       ATOMS=NATOMS-nc
      else
       ATOMS=NATOMS
      endif
c -----
c  *基本的な周期を作るのに必要な原子の数を示すようにさせる
c -----
      WRITE(13,801) ATOMS,NTYPES,zero
c -----
       WRITE(13,811) (real(-nzs(i)),i=1,NTYPES)
c -----LATTICE CONSTANT / atomic-----
      lcux=lcx*UT
      lcuy=lcy*UT
      lcuz=lcz*UT
      WRITE(13,831) lcux,lcuy,lcuz
c -----charge coordinate-----
      do i=1,NATOMS
       if(odlc.eq.'1')then
        if(xv(i).ne.mxv .and. yv(i).ne.myv .and. zv(i).ne.mzv)then
         xvut=xv(i)*UT/lcux
         yvut=yv(i)*UT/lcuy
         zvut=zv(i)*UT/lcuz
         WRITE(13,821) xvut,yvut,zvut,NEQ(i)
        endif
       else
         xvut=xv(i)*UT/lcux
         yvut=yv(i)*UT/lcuy
         zvut=zv(i)*UT/lcuz
         WRITE(13,821) xvut,yvut,zvut,NEQ(i)
       endif
      end do
c -----TRANSITION EXPANSION-----
      WRITE(13,801) ntex,ntey,ntez
c -----SHIFT QUANTITY -----
      WRITE(13,831) sqx,sqy,sqz
c -----CLACULATION RANGE *±----
      WRITE(13,831) crx,cry,crz
c -----------------------------------------------------------------------
c define format
  801 format(3I5)
  811 format(8F10.5)
  821 format(3F20.10,I5)
  831 format(3F20.10)
c -----------------------------------------------------------------------
      END
c =======================================================================

-----------------------------------------------------------------------------
 
QRコード
携帯用QRコード
アクセス数
ページビュー数
[無料でホームページを作成] [通報・削除依頼]