#+TITLE: dump2sq
#+AUTHOR: 大窪 貴洋
#+LANGUAGE: en
#+EMAIL: ohkubo.takahiro@faculty.chiba-u.jp
#+OPTIONS: toc:2 num:nil author:t creator:nil LaTeX:nil \n:nil ^:{}
#+SETUPFILE: https://amorphous.tf.chiba-u.jp/org-html-themes/org/theme-readtheorg.setup

#+MACRO: color @@html:<font color="$1">$2</font>@@
#+MACRO: kakomi @@html:<div style="margin: 1em auto; padding: 1em; width: 100%; border: 2px dotted #ff0000;background-color: #ff;border-radius: 1px">$1</div>@@

\\

# * Author
# [[https://amorphous.tf.chiba-u.jp/][大窪 貴洋]] (千葉院工) 

* はじめに
複数の時間ステップからなる時間ステップの構造データから、動径分布関数
$g(r)$ 、積算配位数 $CN(r)$ 、構造因子 $S(Q)$ を計算するプログラムです。
Triclinicなセルでも計算できます。原子ペア毎の$S(Q)$ を計算した後、原子
散乱因子または中性子散乱長を使って、X-rayと中性子回折実験的で観測され
る $S(Q)$ を計算します。$g(r)$ は、セルサイズ $L/2$ を超えても計算でき
ます。$L/2\sqrt{3}$ まで確保しても良さそうです。産総研の土田さんにアド
バイスいただきました。

https://www.cmu.edu/biolphys/deserno/pdf/gr_periodic.pdf

1ステップ毎に$g(r)$ と$S(Q)$ を計算して最後に時間平均しています。よっ
て、NPT計算でセルサイズが変わったとしても問題なく計算できます。計算コ
ストのほとんどは、原子間距離の計算です。

* 入力ファイル
  拡張子で入力ファイルのformatを判断しています。
** LAMMPS
   dumpコマンドで生成したファイルで、拡張子は.trajectryです。1ス
   テップ目で原子のならび決めて、2ステップ目以降は、同じ原子の並びとし
   て計算しています。時間ステップで原子の並びが変わると、正しい結果に
   ならないので注意してください。データの並びを決め打ちしています。そ
   のため、以下のようにdumpとdump_modifyで元素記号の出力し、idでsortし
   ておく必要があります。
   #+begin_example
   dump              1 all custom 10 hoge.lammpstrj id type element x y z
   dump_modify       1 element Si O sort id
   #+end_example
** XYZ
   xyz形式で、拡張子は.xyzまたは.10です。す。1行目に原子数、2行目に
   Step数と、格子定数を記述しておきます。alpha, beta, gammaがない場合
   は、90, 90, 90とします。3行目以降にデータがならびます。2ステップ目
   以降も同じ形式でデータを並べます。
   #+begin_example
   420
   Step      225837   18.4116484 18.4116484 18.4116484 90 90 90
   Si    10.560073  17.438902  6.614775
   O     17.812448   0.110128  0.3462461
   .
   .
   .
   #+end_example

  
* 出力ファイル
  4つのファイルに拡張子を追加して出力します。basenameは、指定しなけれ
  ば、インプットファイルと同じになります。
** 積算配位数
   拡張子はcnです。原子ペア毎の $CN(r)$ が保存されています。
** 動径分布関数
   拡張子はgrです。原子ペア毎の $g(r)$ と全原子 $G(r)$ が保存されています。
** 構造因子
   拡張子はsqです。原子ペア毎の $S(Q)$ とX線と中性子で観測される $S(Q)$ が保存されています。
** 原子散乱因子と中性子散乱長
   拡張子はcoeffです。中性子散乱長と原子散乱因子を出力します。原子散乱
   因子は、[[https://www.szfki.hu/~nphys/rmc++/downloads.html][xcoeff.f]]から抽出しています。中性子散乱長は、[[https://www.ncnr.nist.gov/resources/n-lengths/][NIST]]のデータを
   参照しています。
   
* インストール
  プログラムはCで作成しています。Cのコンパイラを準備してmakeしてください。
  ダウンロードして単にmakeしてください。

  [[file:dump2sq-1.0.tar.gz][dump2sq-1.0.tar.gz]]

  #+begin_example
  make
  #+end_example

  # iccでOMT_NUM_THREAD > 1にすると、S(Q)がおかしくなる。g(r)は正しいのに…。
  # 原因不明。Threadは1にしとく。
  # icc用のMakefileも準備しています。
  # #+begin_example
  # make -f Makefile.icc
  # #+end_example
  # これで、dump2sqが生成されます。
   
* 使い方
  並列計算のために、環境変数をセットしてスレッド数を指定しておく。
  #+begin_example
  export OMP_NUM_THREADS=4
  #+end_example
  example以下のglass.lammpstrjを引数で渡して実行してみます。
  #+begin_example
  dump2sq example/glass.lammpstrj
  #+end_example
  デフォルトで、$g(r)$ を計算する$r$ の最大値は、もっとも小さいセル長
  さの半分としています。プログラムが実行されると、次のように出力されて、
  4つのファイルが出力されます。
#+begin_example
basename: example/glass
Maximum cell length = 17.284100
Rmax = 8.642050
Start Gr and Sq calculation...
Calculation 0 steps...
Calculation 10 steps...
Calculation 20 steps...
Calculation 30 steps...
Calculation 40 steps...
Calculation 50 steps...
Calculation 60 steps...
Calculation 70 steps...
Calculation 80 steps...
Calculation 90 steps...
Total: 99 steps

example/glass.coeff was generated.
example/glass.cn was generated.
example/glass.gr was generated.
example/glass.sq was generated.
  #+end_example
  
  $g(r)$ の$r$ ステップや最大値、$Q$ の範囲をコマンドラインで指定できます。詳細
  は、ヘルプを参照してください。ヘルプで表示される[ ]の値はデフォルト
  値です。GRmaxは、セルサイズを超えても計算できますが、物理的な意味は
  ないでしょう。
  #+begin_example
  ./dump2sq -h
  usage: dump2sq [options] inputfile
  inputfile format
   LAMMPS:  *.lammpstrj
   FEMTECK: *.10
  
  options:
   -h  --help          show this help message
       --SQmin=%f      Q minimum [0.2]
       --SQmax=%f      Q maximum [30]
       --SQdel=%d      Q delta [0.02]
       --GRmin=%f      r minimum [0.025]
   -r  --GRmax=%f      r maximum [L(minium cell length)*0.5]
   -s  --GRdel=%d      r delta [0.05]
   -o  --output=%s     Output basename [basename]
   #+end_example
    
* tools
** coeffploter.py
   元素を指定して、原子散乱因子と中性子散乱長をプロットする。引数に元
   素名を複数していできる。

** gr2sq.py
   Qの刻みを変えて再計算するために、dump2sqで生成したgrファイルから
   S(Q)を計算しなおす。結晶のパターンを計算するために、$2\theta$ で$Q$
   を指定して計算することもできる。$\theta$ と$Q$ の関係は以下のとおり。
   defaultで波長は、CuのK alpha線の1.54183 Åとしている。
   $$
   Q = \frac{4\pi}{\lambda}\sin(\theta)
   $$
   計算した$S(Q)$はpngとxlsxを生成する

** cn2Iq.c
  dump2sqで生成された*.cnと*.coeffから$I(q)$を計算するCプログラム。
  openmpでスレッド並列する。f*fは足していない。sincのwindow関数をかけ
  る。コンパイルは、
  #+begin_example
  gcc -fopenmp -O3 -Wall cn2Iq.c  -o cn2Iq
  #+end_example
  実行は、*.cnと*.coeffを渡す。I(q)は標準出力に出力されうｒ．
  #+begin_example
   ./cn2Iq fort.cn fort.cn.coeff
  #+end_example
  複数のcnファイルやxslxに保存したい場合は、cn2Iq.pyを使う。

** cn2Iq.py
  subprocessでCのcn2Iqを呼び出して、I(Q)のpngとxlsxを生成する。Qの範囲
  またはthetaの範囲を指定できる。