最大エントロピー法[MEM] PRIMA について

改訂 2010-11-23 by H.M.
補遺 2013-7-24/8-7 by H.M.
MEMプログラムPRIMAは 泉 富士夫氏のホームページ
homepage.mac.com/fujioizumi/index.html
からVENUSのパッケージプログラムの1つとして提供されている。以前のPRIMA Version 3.7(2009/11/5)は不安定であった。
2010年11月現在のPRIMA Ver.3.7 (2010/9)は、一定以上のCPUおよびメモリーを搭載しておれば、Windows XPで安定して走る。
 ここではVer.2.0.5(一部に不具合)と、検証済みのVer.3.5.6(2006-3-23)について記す(2010/4/12)が、Ver.3.7でもたいていは当てはまる。
 PRIMAの使用法の概略
   PRIMAのWINDOWS版実行形式 by Izumi  (注:zipファイル;Ver.3.5.6および3.7)
   PRIMAのテストデータ1(亜硝酸ナトリウムの強誘電構造のテキスト形式のファイル)
   PRIMAのテストデータ2(塩化ナトリウムのテキスト形式のファイル)

 PRIMA(Version 3.5.6および3.7)に関する補遺(2013-7-24)   
  • KClをInitial lambda given by MEM file:0 で実行する MEMファイル0 入力パラメータファイル0   
  • KClをInitial lambda calculated by PRIMA:1 で実行する MEMファイル1 入力パラメータファイル1
      つまり、1 で実行するときは、X線であっても、MEMファイルのtotal b-の箇所に 0.0 と書き込んでおかないと、うまく走らない! (Version 3.7でも走るようになった、と言うべきか)
     正確に言うと、KClのこのデータの場合はnon-linear method計算で初期ラグランジュ係数をPRIMA計算という組み合わせでは、サイクル1回目から発散してNANが出て計算は出来なかったが、ダミーでtotal b-を0.0 として与えると、何故か計算がうまく行くのであった。
     ちなみに、total b-を入れずに、計算がうまく行くのは、
    (i)non-linear methodで初期ラグランジュ係数をMEMファイルに書き込んだ0.001という値を使うと479サイクルで計算は収束する。
    (ii)0-th order approximation methodで初期ラグランジュ係数をMEMファイルの値から出発すると、上記(i)よりもさらに収束は遅いが1000サイクル程度で計算は収束する。
     これに対して、0-th order approximationで初期ラグランジュ係数をPRIMAで計算させると、coefficient to adjust lambdaを大きくとれは初期係数すら求まらず、小さくとれは計算は進むがほとんど収束を水に数百サイクルで計算が自動的に打ちきりとなる。
      初期λの与え方で、以後の計算アルゴリズムに差があり、結果は微妙に異なる。
      いずれにせよ、コマンドプロンプト(CMD)のコマンドラインで
        c:\user\work>PRIMA.exe  < primainp0.txt
      のようにすれば実行できる。ただし、workというフォルダには
        PRIMA.exe,  spgra.dat,   KC+0.mem,  primainp0.txt
      が存在し、かつ、CMDの作業フォルダがカレントフォルダとなっていること(CMDのプロパティの設定)
    
    PRIMAにおけるsigma(F)の重要性についての追記(2013-8-7)
      この点は、つとに言われていることである。しかし、元の回折データを測定・処理する過程において、
     小さすぎるsigma(F)/Fが得られることがある。このような場合には、PRIMAはうまく走らなかったり、
     ゴーストが出たり、原子密度がイビツになったり、丸いはずの原子がバラバラになったような像を与えることがある。
      上記KClで、数割小さいsigma(F)で求めた電子密度図
      本来的には元の測定データを適切に取り直すべきであるが、経験的には、次のようにして「もっともらしい」
     原子密度を表示することもできる。それは、下記のプログラムmegsgcを使い、sigma(F)を一律に何倍かにし、
     PRIMAが収束した時点で得られるR値が、最小二乗法のR値と同程度にすることである。
      上記KClで、sigma(F)を何倍かして求めた電子密度図
      なお、低角の強度が強い反射が消衰効果や強度飽和で失われていると、PRIMAが推論したピーク強度や
     ピークの積分値には不確かさが生ずることに注意しなければならない。
      上記KClで、欠落していた内角の反射を補ってから求めた電子密度図
         原点のK+のピーク強度は体心のCl-のピーク強度より強くなっている。
     もう一つの例は2010年に解析したタウリンの室温構造P21/cである。
     RIGAKU-AFC5Rで測定したリガク提供の標準試料の最小二乗法の結果は1286の独立な反射に対して水素原子を含めてR=3.453%, S=0.598で収束した。
    このFobsとsigmaFを用いてPRIMA(Ver=3.7)を走らせると初期の段階でNANがでて、計算が進行しない。そこで、sigmaFを5倍したmemファイルで行うと、
    129サイクルで収束し、RF=0.027274となった。少し、MEMのRFの方が下がり気味ではあるが、VESTAで表示した電子密度(Isosurface=0.9)では、
    丸い原子密度+水素の膨らみを表現できている。
     なお、この図での水素位置はおなじFobsを用いて最小二乗法で得られた座標で示している。(プロトン位置よりもNやCに寄っている)
     このタウリンの例でも、sigmaFが小さすぎて評価されているとも言える (注)このタウリンの場合、sigmaFを一律大きくしなくても、PRIMA(Ver=3.5.6)で0th order single-pixcel approximationを選択し   CONSTRが7.6程度という収束に達していない段階(サイクル数200)で計算を打ち切れば、同様な電子密度を得られることを   念のためコメントしておく。2010年のページに記載の標準データの結果はそのようにして得られた。

    MEMにおけるsigma(F)と計算の収束(サイクル数、原子のピーク強度)
     極低温において、中性子による核密度を求めるときの留意すべき問題
     ・ NaNO2におけるsigma(F)の影響の例 2010/3/18 by H.M.
     ・ NaNO2におけるサイクル数の原子密度への影響の例 2010/4/2 by H.M.
     ・ PRIMA ver3.5.6によるsigmaFの原子密度への影響例 2010/4/7 by H.M.

    FONDERの室温データでも、sigma(F)に注意が必要である。  2010/11/23 by H.M.

    FONDERデータのsigFを変更するプログラム 2010/10/11 zipに書かれているメモを参照されたし

    cont3d  : 3次元原子密度から2次元データを切り出す   2010/1/26
      MEMの出力(例えば、DEMO0125.pri)をPRIMAで3次元テキストファイルに変換(DEMO0125.3ed)後、 さらに2次元データに変換する。
      概要
      実行形式、ソース、実行例 zip形式 234KB

    memsgc  : prima用入力ファイル *.mem を読み込んで、sigmaFに指定した係数を掛けたものを出力する 2010/4/6
      概要
      実行形式、ソースプログラム zip形式 157KB
    conmem06  : 3次元原子密度から2次元データを切り出す   2006/5/2
      MEMの出力(例えば、DEMO0125.den)を2次元データ(SYFRDAT形式;"DEMO0125.rho")に変換する。
    denmom  : 2次元原子密度(SYFRDAT形式;"DEMO0125.rho")から単位胞中の全原子密度を計算したり、重心(モーメント)を 計算する   2005/12/26
    denmom2  : 計算時に2次関数で内挿して精度を上げる。(大差無いが)  2006/1/12

    もどる
    (c)2010/1/9 by H.M.
    Rev. 2010/11/23;2013/7/24;2013/8/7