第一原理分子動力学プログラム STATE Senri Wiki

エチレン分子（振動解析） †

この例ではエチレン分子を例として振動解析を行います。まずC2H4ディレクトリに移動します。

$ cd C2H4

構造最適化 †

振動解析を行うため、先ず分子構造を最適化します。この例ではGDIISを用いて構造最適化を行います。

• 入力ファイル (nfinp_gdiis)

  #
  # Ethylene molecule in a box: geometry optimization with the GDIIS method
  #
  WF_OPT DAV
  GEO_OPT GDIIS
  NTYP   2
  NATM   6
  TYPE   0
  GMAX    5.00
  GMAXP  15.00
  MIX_ALPHA 0.8
  WIDTH   0.0010
  EDELTA  0.1000D-08
  NEG     10
  FMAX    0.5000D-03
  CELL   12.00  12.00  12.00  90.00  90.00  90.00
  &ATOMIC_SPECIES
   C  12.0107  pot.C_pbe3
   H   1.0079  pot.H_lda3
  &END
  &ATOMIC_COORDINATES CARTESIAN
        1.262722983300      0.000000000000      0.000000000000    1    1    1
        2.348328846800      1.753458668500      0.000000000000    1    1    2
        2.348328846800     -1.753458668500      0.000000000000    1    1    2
       -1.262722983300      0.000000000000      0.000000000000    1    1    1
       -2.348328846800      1.753458668500      0.000000000000    1    1    2
       -2.348328846800     -1.753458668500      0.000000000000    1    1    2
  &END

振動解析 †

構造最適化が終了したら得られた構造を元に新しい入力ファイルを作成します。 ここではgeom2nfinpを使って作成します。

$ geom2nfinp -i nfinp_gdiis -g GEOMETRY -o nfinp_vib

これにより生成されたnfinp_vibを振動解析用に修正します。

• 入力ファイル (nfinp_vib)

  #
  # Ethylene molecule in a box: geometry optimization with the GDIIS method
  #
  TASK   VIB
  WF_OPT DAV
  NTYP   2
  NATM   6
  TYPE   0
  GMAX    5.00
  GMAXP  15.00
  MIX_ALPHA 0.8
  WIDTH   0.0010
  EDELTA  0.1000D-08
  NEG     10
  FMAX    0.5000D-03
  CELL   12.00  12.00  12.00  90.00  90.00  90.00
  &ATOMIC_SPECIES
   C  12.0107  pot.C_pbe3
   H   1.0079  pot.H_lda3
  &END
  &ATOMIC_COORDINATES CARTESIAN
        1.260767348060     -0.000000889176      0.000000061206    1    1    1
        2.337934105040      1.755199776368      0.000000035554    1    1    2
        2.337933682371     -1.755198581491      0.000000037135    1    1    2
       -1.260766004354     -0.000000071340      0.000000050715    1    1    1
       -2.337933757669      1.755199342527      0.000000064907    1    1    2
       -2.337933482763     -1.755199042963      0.000000067944    1    1    2
  &END

  ここで振動計算を行うためには

  TASK   VIB

  TASKにVIBというキーワードを指定します。 古いバージョンのSTATEでは振動計算を行う際にnfvibrate.dataというファイルを作成する必要がありましたが、現在のバージョンではnfvibrate.dataが存在しない場合、各原子についての変位（振幅は0.01 Bohr）を記述したnfvibrate.dataが自動で生成されます。

この入力ファイルを用いて計算を実行すると

nfforce.data

というファイルが生成されます。このファイルには原子変位を与えた際に各原子にかかる力が記述されています。 このデータを元に、以下のようにgifというプログラムを使って振動数を計算します。

$ gif -f nfforce.data

振動数は以下のように出力されます。

             =========             
              SUMMARY              
             =========             

 MODE  WR       : NU(meV)  NU(cm-1)
    1 -0.42D-03 :   12.97    104.63
    2 -0.19D-03 :    8.75     70.59
    3 -0.62D-04 :    4.98     40.19
    4 -0.18D-04 :    2.67     21.54
    5  0.30D-04 :    3.46     27.93
    6  0.28D-03 :   10.70     86.33
    7  0.25D-01 :  100.48    810.43
    8  0.32D-01 :  114.17    920.87
    9  0.34D-01 :  116.25    937.60
   10  0.41D-01 :  128.26   1034.48
   11  0.55D-01 :  148.39   1196.82
   12  0.68D-01 :  165.42   1334.18
   13  0.76D-01 :  175.51   1415.54
   14  0.10D+00 :  201.49   1625.12
   15  0.36D+00 :  379.55   3061.29
   16  0.36D+00 :  381.80   3079.41
   17  0.37D+00 :  388.22   3131.17
   18  0.38D+00 :  393.55   3174.18

振動モードはvib.dataというファイル (gifファイル) に出力されています。 では次に振動モードを可視化してみましょう。可視化を行うためのプログラムとしてgif2xsfが用意されています。gif2xsfを実行するにはvib.dataに加えてXSFフォーマットの原子座標が必要です。 特に分子の場合は非周期のXSFフォーマットの方が扱いやすいので、以下のようにchkinpfを実行します。

$ chkinpf --atom nfinp_vib

これでC2H4.xsfというファイルが生成されます。 次に以下を実行します。

$ gif2xsf -s --xsf C2H4.xsf --gif vib.data --prefix vib

ここで'-s'というオプションは各振動モード毎にファイルを作る（split）オプションで、今の場合prefixとしてvibという名前を指定しているので、vib_[mode#].xsfというファイルが複数生成されます。 '-s'を指定しない場合はaxsf形式ですべての振動モードが一つのファイルに出力されます。 コマンドを実行するとオプションに応じて出力ファイル名あるいはprefixを聞かれるので、それに従って入力します。