表面吸着分子の振動モード計算を効率よく進める方法につ
ついて解説します。


PHASEのPhonon計算では、力計算と、dynamical matrixの対
角化の計算を別々に行うことが出来ます。この機能を利用
することによって一部の原子のみ考慮した振動計算を行う
ことが可能になります。

手順としては以下のようになります。
１．dynamical matrix用に力計算を行う（必要な原子数分）
２．force.dataエディットする
３．dynamical matrixの対角化を行う

それぞれについて以下に説明します。


１．dynamical matrix用に力計算を行う（必要な原子数分）

Phonon計算では、１原子についてｘ、ｙ、ｚの３方向それぞれ
に対して＋、−の２方向の変位をさせて力計算を行います。つ
まり原子数の６倍の回数の力計算が必要になります。大きな系
を計算する場合にはその回数が膨大になるために、PHASEでは
それらを分けて実行する機能があります（現バージョンのPHASE
のマニュアルにはまだ記載されていません）。

そのためにはPHASEの以下のようなforce_calcタグを使用します。
force_calc以下のついている３行です。

Phonon{
        sw_phonon = on
        sw_calc_force = on
        force_calc{
               start = 1, end = 1
        }
        displacement = 0.05
        norder = 1
        sw_polynomial_fit = off
        sw_vibrational_modes = on
}

start、endは何回目の計算を実行するかを示します。上記の例
では原子リストの１番目の原子の−ｘ方向に変位させた場合の
計算のみ行いますし、

        force_calc{
               start = 13, end = 18
        }

とするば、３番目の原子の−ｘ方向から＋ｚ方向までの６通り
を計算します。

この機能を使用して、表面吸着子の振動計算に必要な分だけの
力計算を実行します。例えば、吸着子の原子数が２で、それの
みの計算を実行する場合には（入力ファイルの１，２番目の原
子が吸着子と仮定しています）

Phonon{
        sw_phonon = on
        sw_calc_force = on
        force_calc{
               start = 1, end = 12
        }
        displacement = 0.05
        norder = 1
        sw_polynomial_fit = off
        sw_vibrational_modes = off
}

とします。ここでは吸着子のみ考慮しまましたが、一般的には
基板側の原子もある程度考慮する必要があるために、この数よ
りは大きく取ります。ここでは力計算しかしない（対角化をし
ない）ため、sw_vibrational_modes = off　とします。


２．force.dataエディットする

その後、force.dataをエディットします。このファイルには基
板側の原子も全て含んだ力が出力されていますので、上記で考
慮した原子以外のデータを消去します。書式についてはマニュ
アルに掲載されています。force.dataの行数は1+n*6*(1+n)　
ですので（nは原子数）、エディット後は37行となります。

参考までに、force.dataの原子数を減らすためのプログラムを
添付します（大したものではないですが）。あまり少なくする
と、計算がおかしくなるようなので、基板側もある程度考慮し
たほうがいいと思われます。

このプログラムは対称性については全く配慮していないので、
対称性のある系を計算するときには注意してください。


３．dynamical matrixの対角化を行う

その後、対角化を行います。Phononタグは以下のようにします。

Phonon{
        sw_phonon = on
        sw_calc_force = off
        force_calc{
               start = 1, end = 12
        }
        displacement = 0.05
        norder = 1
        sw_polynomial_fit = off
        sw_vibrational_modes = on
}

この計算では力計算は行わないため、sw_calc_force = off　
とします。その場合、force_calcタグは意味がなくなります。
対角化を行うので、sw_vibrational_modes = onとします。

原子数を変えているので、入力ファイルの原子リストも変更す
る必要があります。この場合では、吸着子の２原子のみにしま
す。基板側も考慮する場合にはその数になります。

この計算が終了すれば、mode.dataに振動モードが出力されます。
出力内容についてはマニュアルにあります。


＊注意点＊

上記にも書きましたが、基板側の原子数が極端に少ない場合、
吸着子の振動に強く関係する原子も省いてしまうことになり、
振動数の結果が本来の結果（全原子を考慮した場合のこと）と
大きく異なることになってしまいます。

ですので、実際の計算においては、上記２．及び３．の手続き
について原子数を変えながら色々と試して、十分収束した値を
採用することをおすすめします。