ABINIT, first (basic) lesson of the tutorial: The H2 molecule, without convergence studies.を少しずつ進めています。水素分子について,
  1. 全エネルギー【済】
  2. 結合長
  3. 電荷密度【済】
  4. 原子化エネルギー
を算出するのが目的でした。其之伍では,2. 結合長を算出します。

算出法は幾つかあります。
  1. 「結合長を与えて全エネルギーを算出」を複数の結合長について行い,全エネルギーの最小値を求める。
  2. 「結合長を与えて原子にはたらく力を算出」を複数の結合長について行い,力の最小値(0になるはずです)を求める。
  3. a.またはb.の複数の計算部分を自動で行う。(フィッティングは後で手動)
  4. 全て自動で行う。(構造最適化)
ここではまず,a.のやり方で進めます(同時に手間なくb.も行うことになります)。実はここまででまだ1度しかabinitで第一原理計算をしていません本ブログの基本方針の通り,なるべく沢山手を動かします。面倒ですが,手を動かせば動かすほど上達するはずです。

tbase1_1.inの以下が,原子群である水素分子を指定してある箇所です。28, 29行目の赤文字部分より,0.7 - (-0.7) = 1.4 Bohrの原子間距離で計算をしていたことが判ります。「これを変えて全エネルギーおよび力を算出」を繰り返します。

023: #Definition of the atoms
024: natom 2           # There are two atoms
025: typat 1 1         # They both are of type 1, that is, Hydrogen
026: xcart             # This keyword indicates that the location of the atoms
027:                   # will follow, one triplet of number for each atom
028:   -0.7 0.0 0.0    # Triplet giving the cartesian coordinates of atom 1, in Bohr
029:    0.7 0.0 0.0    # Triplet giving the cartesian coordinates of atom 2, in Bohr

ここでは,1.0 Bohr〜2.0 Bohrで,0.2 Bohr刻みで計算を行うことにします。

原子間距離1.0 Bohrのときの *.inファイルの例が以下です。28, 29行目以外は変更不要です。また,tbase1_1.filesも変更不要です。ジョブの投入の仕方を忘れていたら,其之壱を再度ご覧ください。
028:   -0.5 0.0 0.0    # Triplet giving the cartesian coordinates of atom 1, in Bohr
029:    0.5 0.0 0.0    # Triplet giving the cartesian coordinates of atom 2, in Bohr

これに対応する *.outファイルの例が以下です。全エネルギーは208行目にあります。
198:  Components of total free energy (in Hartree) :
199:
200:
201:     Kinetic energy  =  1.24897821658742E+00
202:     Hartree energy  =  8.73759293178515E-01
203:     XC energy       = -7.06145116406252E-01
204:     Ewald energy    =  4.34666042075695E-01
205:     PspCore energy  = -1.92143215271889E-05
206:     Loc. psp. energy= -2.88806161023707E+00
207:     NL   psp  energy=  0.00000000000000E+00
208:     >>>>>>>>> Etotal= -1.03682238912322E+00
209:  Other information on the energy :
210:     Total energy(eV)= -2.82133720250743E+01 ; Band energy (Ha)=  -8.1452278001E-01
また,atom1にはたらく力は175行目にあります。
158:  ----iterations are completed or convergence reached----
159:
160:  Mean square residual over all n,k,spin=   1.0119E-12; max=  1.1685E-12
161:    0.0000  0.0000  0.0000    1  1.16852E-12 kpt; spin; max resid(k); each band:
162:   8.55E-13 1.17E-12
163:  reduced coordinates (array xred) for    2 atoms
164:       -0.050000000000      0.000000000000      0.000000000000
165:        0.050000000000      0.000000000000      0.000000000000
166:  rms dE/dt=  2.1948E+00; max dE/dt=  3.8015E+00; dE/dt below (all hartree)
167:     1       3.801452853909      0.000000000000      0.000000000000
188:     2      -3.801452853909      0.000000000000      0.000000000000
169:
170:  cartesian coordinates (angstrom) at end:
171:     1     -0.26458860429500     0.00000000000000     0.00000000000000
172:     2      0.26458860429500     0.00000000000000     0.00000000000000
173:
174:  cartesian forces (hartree/bohr) at end:
175:     1     -0.38014528539087    -0.00000000000000    -0.00000000000000
176:     2      0.38014528539087    -0.00000000000000    -0.00000000000000
177:  frms,max,avg= 2.1947698E-01 3.8014529E-01   0.000E+00  0.000E+00  0.000E+00 h/b
178:
179:  cartesian forces (eV/Angstrom) at end:
180:     1    -19.54785488759569    -0.00000000000000    -0.00000000000000
181:     2     19.54785488759569    -0.00000000000000    -0.00000000000000
182:  frms,max,avg= 1.1285959E+01 1.9547855E+01   0.000E+00  0.000E+00  0.000E+00 e/A

これら一連の計算は,本家チュートリアル1 〜其之壱〜のときと同じ作業ディレクトリD:¥abinit¥L01で行えば良いのですが,1つ注意があります。以下は,D:¥abinit¥L01で1.0 → 1.2 → 1.4 → 1.6 → 1.8 → 2.0 Bohrの順に6回の計算(1.4 Bohrは再度計算)を行った後に,dirコマンドでファイル一覧表示をしたものです。
D:¥abinit¥L01>dir 
 ドライブ D のボリューム ラベルは XXXXX です
 ボリューム シリアル番号は XXXX-XXXX です

 D:¥abinit¥L01 のディレクトリ

2014/11/15  01:12    <DIR>          .
2014/11/15  01:12    <DIR>          ..
2014/11/15  01:15               378 01h.pspgth
2014/11/20  18:22            28,131 log.txt
2014/11/15  01:17                65 tbase1_1.files
2014/11/20  18:21             2,933 tbase1_1.in
2014/11/15  01:21            17,450 tbase1_1.out
2014/11/20  18:16            17,451 tbase1_1.outA
2014/11/20  18:17            17,451 tbase1_1.outB
2014/11/20  18:19            17,451 tbase1_1.outC
2014/11/20  18:20            17,451 tbase1_1.outD
2014/11/20  18:21            17,369 tbase1_1.outE
2014/11/20  18:22            17,451 tbase1_1.outF
2014/11/15  02:20           666,179 tbase1_1.xsf
2014/11/20  18:22             5,141 tbase1_1o_DDB
2014/11/20  18:22           217,634 tbase1_1o_DEN
2014/11/20  18:22               236 tbase1_1o_EIG
2014/11/20  18:22               548 tbase1_1o_EIG.nc
2014/11/20  18:22             4,512 tbase1_1o_GSR
2014/11/20  18:22             2,160 tbase1_1o_OUT.nc
2014/11/20  18:22            34,798 tbase1_1o_WFK
              19 個のファイル           1,084,789 バイト
               2 個のディレクトリ  XX,XXX,XXX,XXX バイトの空き容量

D:¥abinit¥L01>
*.outファイルの末尾にA〜Fが付いて,メイン出力ファイルが複数(計7つ)あります。
abinitでは,メイン出力ファイルが既に存在するとそのファイルを保持したまま,新たに末尾にA〜Zを付け加えた *.outファイルが出力されます。ここで注意したいのは,最新の計算結果がどの *.outファイルに出力されるか?です。ファイルのタイムスタンプを見れば判断できるので,日時分を見ます。この例では,本家チュートリアル1 〜其之壱〜で得たメイン出力ファイルは *.outです。その後1.0 Bohrについて得たメイン出力ファイルは *.outA,1.2 Bohrが *.outB,その次が *.outC...という具合になります。なお,メイン出力ファイル以外の出力ファイル(*o_DENなど)は最後のジョブのものに更新されます。

さて。原子間距離を1.0 Bohr〜2.0 Bohrの範囲で0.2 Bohr刻みで計算した結果が以下の表になります。

原子間距離 (Bohr) 全エネルギー (Ha) 力 (Ha/Bohr)
1.0 -1.03682238912322E+00 -0.38014528539087
1.2 -1.08658147848712E+00 -0.14716643363971
1.4 -1.10372242131361E+00 -0.03740558871227
1.6 -1.10511251078528E+00 0.01763045171617
1.8 -1.09828719403481E+00 0.04785128411887
2.0 -1.08679728680317E+00 0.06544241709374

これを図示したのが下の図です。
ByHand
全エネルギーには最小値,にはゼロの点が存在することが分かります。表および図より,水素分子の再安定な構造は,原子間距離が1.5 Bohr少々(例えば,1.52 Bohr)と予想できます。

ただ,データ6点ではフィッティングするには少々雑です。もっと多くの計算データ点が欲しいところです。刻み幅を0.1 Bohrや0.01 Bohrに小さくして,データ点数を増やしさえすれば良いのですが,人力ではなくコンピュータ(abinit)に任せることにします。次回この方法を示します。基本方針に反するようですが,基本方針は目的達成のための手段ですので。

以上で,結合長の算出(人力)が完了です。

おまけ
*.outの末尾追加文字のA〜Zの全てが使われてもなおジョブを投入した場合,エラーとなり計算されません。不要ファイルを削除・整頓するよう促されます。つまり,1つのディレクトリに27つまでのメイン出力ファイルが出力可能です。