研: phreeqc初心者のおべんきょう(1)

Fe2+のFe(III) oxideへの吸着を計算しなくてはいけないことになって、共同研究者から

「PHREEQCやMINTEQLなら異なるpH条件で様々なFe(III) oxideに対するFe2+の吸着量を計算できるから便利だよ」

と助言を受け、なるほどそうか、じゃあ良く聞くPHREEQC（ちなみに読み方はフリクシー、かわいい）を使ってみるか…と思い立ち、とりあえず覚えたことをメモ。
間違いなど教えていただけると嬉しいです。


1. とりあえずこちらからMacOS(Intel) batch versionをダウンロード


2. dmgをマウントするとphreeqc-2.x（わたしがダウンロードしたバージョンはphreeqc-2.18.3）フォルダがあるので、すきなところに保存


3. phreeqc-2.x直下のアップルスクリプトもしくはbin/phreeqcよりターミナルで起動。以下ターミナルでの利用。binにパスを通す。


4. とりあえずテスト。作業用フォルダを好きなところに作成して、作業用フォルダにphreeqc-2.x/examplesからex1とdatabaseからphreeqc.datをコピペ、作業用フォルダに移動後ターミナルよりphreeqcで起動

Name of input file?
ex1
Input file: ex1

Name of output file?
Default: ex1.out

Output file: ex1.out

Name of database file?
Default: phreeqc.dat

Database file: phreeqc.dat

と、ダイアログに沿って入力すると、ex1.outに結果が出力される。

ex1は海水中のUのスペシエーションの例で、まだ完全に理解していないですが


SOLUTION 1  SEAWATER FROM NORDSTROM ET AL. (1979)
# 以下は濃度などの設定
        units   ppm
        pH      8.22
        pe      8.451
        density 1.023
        temp    25.0
        redox   O(0)/O(-2)
        Ca              412.3
        Mg              1291.8
        Na              10768.0
        K               399.1
        Fe              0.002
        Mn              0.0002  pe
        Si              4.28
        Cl              19353.0
        Alkalinity      141.682 as HCO3
        S(6)            2712.0
        N(5)            0.29    gfw   62.0
        N(-3)           0.03    as    NH4
        U               3.3     ppb   N(5)/N(-3)
        O(0)            1.0     O2(g) -0.7
# ここまで濃度などの設定

# SOLUTION_MASTER_SPECIES: 元素名、主要な種の荷電状態を含む化学式、アルカリ度への寄与、化学式のグラム重量、元素のグラム重量

SOLUTION_MASTER_SPECIES
        U       U+4     0.0     238.0290     238.0290
        U(4)    U+4     0.0     238.0290
        U(5)    UO2+    0.0     238.0290
        U(6)    UO2+2   0.0     238.0290

# SOLUTION_SPECIES: 会合反応の反応式、平衡定数をlog_k、25度のときの反応エンタルピーをdelta_hで与える

SOLUTION_SPECIES
        #primary master species for U
        #is also secondary master species for U(4)
        U+4 = U+4
                log_k          0.0
        U+4 + 4 H2O = U(OH)4 + 4 H+
                log_k          -8.538
                delta_h        24.760 kcal
        U+4 + 5 H2O = U(OH)5- + 5 H+
                log_k          -13.147
                delta_h        27.580 kcal
        #secondary master species for U(5)
        U+4 + 2 H2O = UO2+ + 4 H+ + e-
                log_k          -6.432
                delta_h        31.130 kcal
        #secondary master species for U(6)
        U+4 + 2 H2O = UO2+2 + 4 H+ + 2 e-
                log_k          -9.217
                delta_h        34.430 kcal
        UO2+2 + H2O = UO2OH+ + H+
                log_k          -5.782
                delta_h        11.015 kcal
        2UO2+2 + 2H2O = (UO2)2(OH)2+2 + 2H+
                log_k          -5.626
                delta_h        -36.04 kcal
        3UO2+2 + 5H2O = (UO2)3(OH)5+ + 5H+
                log_k          -15.641
                delta_h        -44.27 kcal
        UO2+2 + CO3-2 = UO2CO3
                log_k          10.064
                delta_h        0.84 kcal
        UO2+2 + 2CO3-2 = UO2(CO3)2-2
                log_k          16.977
                delta_h        3.48 kcal
        UO2+2 + 3CO3-2 = UO2(CO3)3-4
                log_k          21.397
                delta_h        -8.78 kcal

# PHASES: 解離反応の化学式、平衡定数をlog_k、25度のときの反応エンタルピーをdelta_hで与える

PHASES
        Uraninite
        UO2 + 4 H+ = U+4 + 2 H2O
        log_k          -3.490
        delta_h        -18.630 kcal
END



Uの反応などはデータベース（phreeqc.dat）に含まれていないので、追加している。

あれ…これはもしや、attachment constantとかは自分で入力しなければいけないのだろうか…?それならパラメータがあれば自分でプログラム組む方が速いような...

とりあえず、ぼちぼちと試してみます。

本：思考の整理学

先週「つん読」していた本を読み終わったので、昨日は新しい書籍を購入しに大阪のジュンク堂へ。

数冊を補充し、昨日は

思考の整理学　外山滋比古著　ちくま文庫

を読了。東大・京大で５年間販売冊数第１位という帯に惹かれました。

しまった、これはもっと早く読んでおくべきだった。
そして大学生の皆さん、特にこれから卒論を書くことになる皆さん、是非読んでください。
この本には、卒論を書くとはどういうことか、卒論を書くためにどうすべきかの示唆がふんだんにちりばめられています。

以下わたしの心に止まったことを箇条書きにて

1 グライダーにエンジンを持たせる教育とは？
2 快活な頭の状態を意識的に作る（忘却することの大事さを含む）
3 １つのアイディアにとらわれすぎないこと
4 整理する、寝かせる、書く、話す
5 創造の芽を育てる

2-4は近年その重要性を痛感しており、自分なりに気をつけていたけれど、1と5は未だ難しく。

わたしの場合、4の「整理」のためのアイディアメモにはWunderlistを利用しており、「寝かせる・書く」のためには研究用雑記ノートと研究テーマ（プロジェクト）に特化したノートを作成、結果のまとめはLaTeXか場合によってはWordで、「話す」は近くに研究の話ができる方が少ないため（ぐすん）メールとEvernoteで行うことが多いです。

ノートに書くことが未だに好きなのは、「考える」「寝かせる」が一緒にできるからなのではないかと思う次第。

そういう理由から、わたしは講義中に板書を撮影することは不可としています。
何が大事か考えながら、自分なりにまとめてメモをすることのほうが重要ですよ。

四方山：清須会議

下田の出張が延期になったので（Iさんとそのご家族が早く元気になりますように！）、昨日は来週日曜日の女子中高生のための関西科学塾の講義資料の作成。中高生の皆さんを対象に二酸化炭素濃度上昇に伴う海洋酸性化のシミュレーションをします！

今日は放射性物質関連のお仕事を再開。そして久々に定時に大学を出て、清須会議を観てきてしまいました。今回は三谷さんがメガホンを取る初の時代物、さらには豪華な顔ぶれが出演とのことで楽しみにしていました。実際に楽しめました。三谷さんの作品はいつも人物間の駆け引きが絶妙ですが、今回も健在。役者の皆様も本当に素晴らしかったですが、個人的に小日向文世さん演じる丹羽長秀が一番素敵でした...！

三谷さんの映画作品に連なるテーマとして夢と挫折、仕事と責任が描かれているように感じており、活力をいただいております。

四方山：下田に出張してきます

昨日は鉄還元細菌による鉄還元動態のモデルのアイディアをまとめ、簡単なシミュレーション結果も添えて共同研究者に送るところまでなんとか終了（今週の目標にしていたので良かった良かった）。

忙しいときもすぐに返事をくれる人で、いつもポジティブなコメントをくれるのが嬉しい（今回の場合はvery nice!でした）。自分もそうありたいと思います...いや、思っているだけではいかんですね、そうせねば。

週末から火曜日まで共同研究の打ち合わせで筑波大学下田臨海実験センターに滞在してきます。ついにマンガン、ヒ素に関わる研究にも乗り出せるか？！（MnはCe anomalyの研究の時に少しだけ関わったけれど）あと、DOC関連のことも話し合えたらな...楽しみ！

四方山：今日は臨時休業日でした

...ということに大学に着いてから気がつきました。
学園祭中はサイエンスオープンラボにたくさんの皆様にお越しいただき、誠に有り難うございました！
情報科学科の催しへの来場者数は計679名でした。

今日はひたすら鉄還元細菌についての論文を読んでおりました。Fe(II)による還元の阻害とか、electron shuttle（日本語でなんと訳すのかしら？）とか、どう考えるべきかな...！

そして16時半からは高知大学の加藤元海さんによるセミナー、その後懇親会でした。

MathematicaでのStepwise変数選択（変数減少法のみ）

MathematicaでStepwise変数選択（変数減少法のみ）のプログラム作成してみました。

改良などご助言ありましたらご連絡いただけると嬉しいです。

AICやBICなどによる選択法のプログラムも作ったのですが、また次回にでも。

(* テスト用のデータ作成 *)

tmp = {RandomReal[{0, 20}, 40], RandomReal[{0, 20}, 40] , RandomReal[{0, 20}, 40],
   Table[0.5*x + RandomVariate[NormalDistribution[0, 1]], {x, 1, 40, 1}],  Table[1*x + RandomVariate[NormalDistribution[0, 1]], {x, 1, 40, 1}]};
data = {"x1", "x2", "x3", "x4", "y"};
data = Join[{data}, Transpose[tmp]];
data // TableForm

(* dataの最後のリストを従属変数として、他の変数を独立変数として重回帰し、結果をreslistに保存。推定されたパラメタに対するt検定のp値が全てcritを下回ったら終了（誤差項は除いています）。そうでなければ最大のp値が返されたパラメタを持つ独立変数を除き再度重回帰。独立変数の数が0になったら変数選択を終了。*)

crit = 0.2;
label = data[[1]];
df = Transpose[Drop[data, {1}]];
xlist = Table[Subscript[x, i] , {i, 1, Length[df] - 1}];
reslist = {}; tmplist = {};
While[Length[df] > 1,
 {
  data = Transpose[df];
  lm = LinearModelFit[data, xlist, xlist];
  AppendTo[reslist, lm];
  tmp = Drop[lm["ParameterTable"][[1]][[1]][[All, 5]], 2];
  AppendTo[tmplist, tmp];
  Print[label];
  If[Length[Select[tmp, # > crit &]] == 0, Break[]];
  remove = 1; maxP = -1;
  For[i = 1, i <= Length[tmp], i++,
   If[tmp[[i]] > maxP, {remove = i, maxP = tmp[[i]]}]];
  df = Delete[df, {remove}];
  label = Delete[label, {remove}];
  xlist = Delete[xlist, {remove}];
  }
 ]

(* ベストモデルの抽出 *)

Normal[Last[reslist]]
Last[reslist]["ParameterTable"]

論：微生物群衆の情報は窒素循環を知るためには役に立たない？

---

Do we need to understand microbial communities to predict ecosystem function? A comparison of statistical models of nitrogen cycling processes
Graham et al. (2014) Soil Biol. Biochem., 68, 279-282

---

微生物群衆の遺伝子存在量データを含む土壌多変量データに対して、硝化速度やN2Oフラックスを従属変数としてステップワイズ変数選択、重回帰分析を用い、微生物群衆データがない方がモデルの説明力が高いことを示した論文。

それで微生物群衆データの情報が窒素循環を知るためには有益でないとするのはちと乱暴すぎるのでは？
まず線形和だけで従属変数を説明できるとする考え方は仮定に過ぎないので、その仮定に基づいたら微生物群衆データ情報は役に立たないという結論になった、というだけの話だと思うのです。
非線形関係を考慮したらどうなるのか？という議論も必要。

統計モデルはプロセスを駆動するメカニズムを知るための指標を提示するものであって、統計モデルを構築することでメカニズムを理解することにはならないと考えます。

---