feram: 強誘電体の高速分子動力学シミュレーション
feramは強誘電体薄膜の高速分子動力学シミュレーターです．
Table of Contents:
feramは強誘電体薄膜の高速分子動力学シミュレーターです．
feram用のパラメータファイルはコメント行と 'tag = value(s)' というタグ行 からなるテキストファイルです．パラメータファイルのファイル名を foo123 または foo123.feram とした場合，出力ファイルは foo123.log, foo123.0000000000.cood, ... などとなります． それぞれのパラメータをどう決定するかは parameters/parameters.html に説明があります．
'#' で始まる行はコメント行で，無視されます． 空行も無視されます．
# This is a comment line. # Here are two more # comment lines.
.feram入力ファイルの１行目にシバンを書いておいてスクリプトのように実行することができます．
タグ行には前後がスペースで挟まれた等号 ' = ' が入っていなくてはいけません． 次のとおり，その左にタグを右に設定値を置きます：
tag = 1.0 tag = -2.0 -3.0 -4.0 tag = 5.0 6.0 7.0
以下に現在feramで利用可能なタグの一覧を示します．
.logファイル（詳細は「出力ファイル」の項を参照）に出力するメッセージの量を指定します．
verbose = 0
0のとき，.logファイルは作られません（未実装）．
verbose = 1
1（デフォルト値）のとき，.logファイルには最初と最後のイテレーションのエネルギーの値のみがレポートされます．
verbose = 2
2のとき，.logファイルにはすべてのイテレーション毎にエネルギーの値がレポートされます．.logファイルのサイズは巨大になります．
verbose = 3
3のとき，.logファイルにはcallされたサブルーチンがメッセージとして出力されます．
verbose = 4
4のとき，.logファイルにはさらに詳細なメッセージが出力されます．
シミュレーションの方法を指定します．
method = 'md'
'md'を指定すると能勢−ポアンカレのサーモスタットを使ったカノニカルアンサンブルでの分子動力学シミュレーションを行います．
method = 'vs'
'vs'を指定すると速度スケーリング法を使ったカノニカルアンサンブルでの分子動力学シミュレーションを行います．
method = 'lf'
'lf'を指定すると蛙飛び (leapfrog) 法を使ったミクロカノニカルアンサンブルでの分子動力学シミュレーションを行います．
method = 'mc'
'mc'を指定するとモンテカルロシミュレーションを行います．まだインプリメントされていません．
method = 'hl'
'hl'の機能はちかぢか削除されます． ヒステリシスループ (hysteresis loop) のシミュレーションには n_E_wave_periodとE_wave_typeとを使って下さい．
圧力をGPa単位で．過小評価されるTcを負圧をかけることにより補正することがある．
GPa = -5.0
また，３つの値を指定することにより垂直応力をかけることもできる．
GPa = 0.0 0.0 -3.0
温度．
kelvin = 100
Rにあるunit cellの双極子の長さu(R)の有効質量． 単位はatomic mass unit．
もしacoustic_mass_amuが設定されないと，音響変位，すなわち局所歪みは{u(R)}を元に最適化される． acoustic_mass_amuとして音響変位の有効質量が
acoustic_mass_amu = 46.44
と与えられると，音響変位についても分子動力学計算が行われる． 単位はatomic mass unit．一般的に全５原子の平均質量とする． 今のところ，method = 'md'では音響変位の分子動力学計算はできません． method = 'vs'もしくはmethod = 'lf'で使って下さい．
能勢のQ．MD計算に必要なパラメータで，熱欲と系とのカップリングの強さを表す． 大きなsupercell，高温には大きな値を， 小さなsupercell，低温には小さな値を使う．
Q_Nose = 14.4
システムの構造を指定します． 'bulk' は無限に周期的なバルクです． 'film' は「自由な」薄膜キャパシタです．電極は短絡されています． 'epit' はエピタキシャル歪み（epi_strainで指定）を伴った薄膜キャパシタです．これも電極は短絡されています．
bulk_or_film = 'epit'
System size, L_x, L_y and L_z. They must be equal or larger than 4.
L = 32 32 4
計算の高速化のためのy方向のパディングの値． 最適なパディングの値はシステムやコンパイラに依存します． 64x64x64 と 96x96x96 のシステムサイズなら
padding_y = 5
128x128x128 なら
padding_y = 1
など，いちどベンチマークを取って最適値を出しておくとよいでしょう． padding_yの既定値は3です． 詳しくは src/coord_module.F を参照してください．
0: dead layer はなし． 1: 片側だけに dead layer がある． 2: 両側に dead layers がある．
gap_id = 1
dead layer の厚さは常に 1 unit cell．
dead layerの分極．単位はu(R)と同じÅ．
gap_dipole_u = 0.0 0.0 0.16
ペロブスカイト強誘電体の立方晶の格子定数をÅで指定します．
a0 = 3.99
エピタキシャル歪み．
epi_strain = -0.01
epi_strainはbulk_or_film = 'epit'のときにだけ有効であることに注意して下さい． 詳しくは src/optimize-homo-strain.F を見てください．
タイムステップの長さをピコ秒で．
dt = 0.002
thermalizeにかけるステップ数，物理量の平均を取るのにかける ステップ数，.coordファイルを書き出す頻度，ヒステリシスループの 計算のときに.hlファイルに書き出す頻度．
n_thermalize = 40000 n_average = 10000 n_coord_freq = 50000 n_hl_freq = 5000
external_E_fieldの振幅を持つ外部電場を時間により振動させるときの周期． E_wave_typeに'triangular_sin'か'triangular_cos'を指定します． src/27example-BaTiO3-new-param-E_wave/ 以下に例があります． 'triangular_sin'か'triangular_cos'の場合 n_E_wave_periodは４で割り切れる正の整数でなくてはなりません． 'ramping_off'と'ramping_on'も使うことができます． 先にE_wave_typeに'ramping_off'または'ramping_on'を指定しておけば n_E_wave_periodはどんな正の整数でもかまいません．
ヒステリシスループ (hysteresis loop) のシミュレーションを n_hysteresis_loop_continue回目のイテレーションから 再開することができます．
n_hysteresis_loop_continue = 20000
このタグを
coord_directory = 'coord_files'
と設定すると，feramは.coordファイルをそのディレクトリに書き出します． この例では coord_files/foo.0000020000.coord などとなります． 巨大な.coordファイルをNFSのホームディレクトリでなくローカルの ハードディスクに書き出したい場合に便利です．特に指定がなければカレント ディレクトリに.coordファイルは書き出されます．
このタグを 'never' に設定すれば，.coordファイルをけっして書き出しません．
coord_directory = 'never'
このタグを
distribution_directory = 'distribution_files'
と設定すると，foo123.distribution と foo123.distribution3d とは distribution_files/ ディレクトリに書き出されます．とくに後者は ファイルサイズが約180MBと大きいので，ローカルのハードディスクに 書き出したい場合があるかもしれません．特に指定がなければカレント ディレクトリに２つのファイルは書き出されます．
もし .distribution と .distribution3d とが必要なければ
distribution_directory = 'never'
としてください．こうすることで計算時間も大幅に節約できます．
外部電場Eを[V/Angstrom]の単位で． methodが'md', 'lf', 'mc'のときはこの電場で一定． methodが'hl'のヒステリシスの計算ではn_averageステップの間に電場が+Eから-Eまで変化する．
external_E_field = 0.00 0.00 -0.01
また，E_wave_typeを指定すると，Eはn_E_wave_periodの周期で振動する．
4次の多項式の係数．
P_kappa2 = 5.502 [eV/Angstrom^2] # P_4(u) = kappa2*u^2 + alpha*u^4 P_alpha = 110.4 [eV/Angstrom^4] # + gamma*(u_y*u_z+u_z*u_x+u_x*u_y), P_gamma = -163.1 [eV/Angstrom^4] # where u^2 = u_x^2 + u_y^2 + u_z^2
6次の項の係数．
u^8の係数．
短距離相互作用の係数．
j = -2.648 3.894 0.898 -0.789 0.562 0.358 0.179 j(i) [eV/Angstrom^2]
弾性定数 B11=C11*a0^3 [eV]．
B11 = 126. B12 = 44.9 B44 = 50.3 [eV]
歪みとuとのカップリング係数．
B1xx = -211. [eV/Angstrom^2] B1yy = -19.3 [eV/Angstrom^2] B4yz = -7.75 [eV/Angstrom^2]
疑似乱数生成器 (RNG) のための２つの整数の種を指定します． Marsaglia-Tsang 64-bit uiniversal RNG が使われています．
seed = 143456789 987654321
デフォルトの値は123456789 987654321です． 次のようにjot(1)を使って種を生成するのもよいでしょう．
$ jot -r 2 1 2147483648 639779840 1297280417
bash(1)の$((RANDOM))機能を使うのもよいでしょう．
$ echo $((RANDOM*65536+RANDOM)) $((RANDOM*65536+RANDOM)) 2146723680 1863789493
dipoleの初期状態の平均値とばらつき具合． ２つの値が
init_dipo_avg = 0.10 0.00 0.00 [Angstrom] # Average of initial dipole displacements init_dipo_dev = 0.02 0.02 0.02 [Angstrom] # Deviation of initial dipole displacements,
のとき，４つの双極子があったら
0.09 0.01 0.02 0.11 -0.02 0.01 0.08 0.02 -0.03 0.12 -0.01 0.00
などとなるだろう．.restartファイルが見つからない場合に限り init_dipo_avgタグとinit_dipo_devタグとで指定される初期状態が 疑似乱数を使って作られます．
unit cellあたりの有効電荷．単位はe．
Z_star = 9.956
光学誘電率ε_∞．
epsilon_inf = 5.24
この値は自動的に最適な値に設定されるので，ユーザが指定する必要はない．デバッグのために用いる． src/dipole-dipole-long-range.Fで双極子のエワルト和に使われる任意のκである．
kappa = 0.15
もし
plot_dispersion = .true.
が指定されたなら，分散関係のプロットのためのデータファイルが書き出されます． デフォルトは.false.です．
plot_dispersion = .true.のときは L_x, L_y, L_z はすべて同じ偶数でなくてはいけません．
L = 32 32 32
がよいでしょう．
例題の src/08example-BaTiO3-dispersion/ を見ていただくか、たとえば， 次のようにテストスクリプト dispersion_check.sh を実行すると分散関係を描くことができます．
$ tar xf SOMEWHERE/feram-X.YY.ZZ.tar.xz
$ cd feram-X.YY.ZZ
$ mkdir builddirectory
$ cd builddirectory
$ ../configure && make -j4
$ make check TESTS='dispersion_check.sh dispersion_32x32x32.inhomo-K.gp'
$ less dispersion_32x32x32.feram # the input file
$ less dispersion_32x32x32.gp # the GNUPLOT script
$ gv src/dispersion_32x32x32.long+short.energy.eps &
$ gv src/dispersion_32x32x32.long+short.eps &
$ gv src/dispersion_32x32x32.long+short.interaction.eps &
$ gv src/dispersion_32x32x32.inhomo-K.eps &
# Previews of plots. You may also use evince(1) instead of gv(1).
[Takeshi Nishimatsu el al.: Phys. Rev. B 82, 134106 (2010), http://dx.doi.org/10.1103/PhysRevB.82.134106 ] の中の Fig.3 や parameters/parameters.html の中の Fig.5 と比較してください．
feramが
$ feram foo.feram bar.feram
と実行されて，bar.feramの中で
continue = .true.
が指定されたなら，bar.feramはfoo.feramの結果を初期状態として計算を始めます． .restartファイルは使われません． デフォルトは.false.です．
feramを次のようにパラメータ・ファイル "foo123.feram" を引数として
$ feram foo123.feram
と実行すると，feramは "foo123.restart" ファイルを探して， 見つかればそれを初期状態として用います．.restartファイルの書式は .coordファイルと同じです．.restartファイルが見つからない場合は init_dipo_avgタグとinit_dipo_devタグとで指定される初期状態が 疑似乱数を使って作られます．
２次の原子間力定数行列 (qudratic inter-atomic force constant (IFC) matrices) を.quadraticファイルから読み込ませることができます． 3x3の行列が(Lx/2+1)*Ly*Lz個ひつようです．.quadraticファイルを使う場合， j1-j7 をすべてゼロにして下さい．また，P_kappa2を *elevation* のために使って下さい．src/25example-BaTiO3-read-quadratic/ にその 使用例があります．
feramを次のようにパラメータ・ファイル "foo123.feram" を引数として
$ feram foo123.feram
と実行すると，feramは "foo123.defects" ファイルを探して， 見つかればそこに指定された双極子u(R)を固定します． このファイルが見つからない場合，上記のような「格子欠陥」は導入されません．
.defects 入力ファイルの書式は
0 4 0 0.1 0.1 0.3 4 6 13 0.0 0.0 0.3 : ix iy iz u_x u_z u_z :
と固定するu(R)の位置とベクトルを順に１行ずつ書きます．
src/feram_defects_maker.rb はこの .defects 入力ファイルを作るのに便利です．
feramを次のようにパラメータ・ファイル "foo123.feram" を引数として
$ feram foo123.feram
と実行すると，feramは "foo123.localfield" ファイルを探して， 見つかればそのファイルの中で指定されたサイトに指定された局所場をかけます． .localfieldファイルの書式は
4 6 13 0.00 0.00 0.10 : ix iy iz Ex Ey Ez :
局所場の単位は V/Angstrom です．このファイルが見つからない場合， 局所場はもちろんどこのサイトでも０になります．
各サイトにランダムな局所場をかける場合は付属のferam_random_field_generatorを使って.localfieldファイルを作るとよいでしょう．
$ ./feram_random_field_generator Lx Ly Lz mu_x mu_y mu_z sigma_x sigma_y sigma_z seed1 seed2 $ ./feram_random_field_generator 32 32 243 0.0 0.0 0.0 0.02 0.02 0.02 123456789 987654321
各サイトの acoustic modulation a(R) を指定します． acoustic modulation に依存するエネルギは
E_acoustic_modulation = modulation_constant Sigma_R_alpha eta_alpha(R) a(R)
となります．
必須ではありませんが，FFTW用のwisdomファイルを用意すると， より最適なFFTのサブルーチンが選ばれるので， 全体の計算時間を短縮できることがあります． 逆に全体の計算時間が長くなってしまうこともあります． トライ・アンド・エラーで最適な計算条件を探してください．
wisdomファイルはferamと同時にsrc/内にビルドされるferam_fftw_wisdomコマンドで生成できます． たとえば，src/18example-benchmark/ ディレクトリ内で forward.feram のシステムサイズ Lx*Ly*Lz=32*32*243 の計算をするなら，
$ w # Before executing feram_fftw_wisdom, confirm that load average is zero and there is no background job.
$ OMP_NUM_THREADS=6 ../feram_fftw_wisdom 1000 32 32 243 3
:
It takes a few minutes or a few hours depending on the system size.
:
$ mv wisdom_new wisdom
$ OMP_NUM_THREADS=6 ../feram forward.feram
$ grep FFTW_WISDOM forward.log
feram_common.F: 47: FFTW_WISDOM: Successfully imported FFTW wisdom in current directory.
$ grep TIMING_REPORT forward.log
とします．
feram_fftw_wisdomに与えるコマンドライン引数は
feram_fftw_wisdom <N_TIMES> <Lx> <Ly> <Lz> <padding_y>
です．
FFTW wisdomについての詳細は http://www.fftw.org/doc/Other-Important-Topics.html をご覧下さい．もしMKLの中のFFTW3ラッパーを使っている場合は wisdomファイルを準備する必要はありません．
ソフトウエアパッケージの次のディレクトリに入力ファイル（と結果）が入っています．
feram-X.YY.ZZ/src/03example-BaTiO3-bulk-phase-transition/ feram-X.YY.ZZ/src/09example-BaTiO3-bulk-leap-frog/ feram-X.YY.ZZ/src/17example-PbTiO3-100-900K/ feram-X.YY.ZZ/src/20example-BaTiO3-new-param-pyro/ feram-X.YY.ZZ/src/21example-KNbO3/ feram-X.YY.ZZ/src/22example-BaTiO3-new-param-hl/ feram-X.YY.ZZ/src/23example-BaTiO3-new-param-hl-film-acoustic-mass/ feram-X.YY.ZZ/src/24example-BaTiO3-new-param-heating-cooling/ feram-X.YY.ZZ/src/25example-BaTiO3-read-quadratic/ feram-X.YY.ZZ/src/26example-BaTiO3-acoustic-MD/ feram-X.YY.ZZ/src/27example-BaTiO3-new-param-E_wave/ feram-X.YY.ZZ/src/28example-PbTiO3-elastocaloric-770K/ feram-X.YY.ZZ/src/29example-BaTiO3-new-param-defects/ feram-X.YY.ZZ/src/30example-BaTiO3-new-param-ramping-off/ feram-X.YY.ZZ/src/31example-BaTiO3-new-param-thermal-conductivity/ feram-X.YY.ZZ/src/32example-BaTiO3-old-64x64x1024-benchmark/ feram-X.YY.ZZ/src/34example-BST/
ご支援に感謝します：
http://loto.sourceforge.net/feram/ is hosted by .
