GROMACS で計算した、メタンを水中に挿入したときのエネルギーの値が正負異常に大きい

[uran ry]

こんにちは、 お世話になります、
ERMOD のテスト計算として、　GROMACS で計算したトラジェクトリをつかってメタンを水中に入れるときの自由エネルギーを計算しようとしています

今回質問したのは、エネルギーの分布が
論文のと大きくことなり(kJ とkcal の差を含めても)
結果として溶媒和自由エネルギーの実験値の-9.2 kJ/mol か -2.21 kcal/mol　から大きくずれているので、
いろいろ試しましたがわからなかったので質問しました。


系としては以前に、
Andrey Frolov さんが投稿されたHydration free energy of amino acid analogs: methane

と同じ系です。計算条件は
GROMACSのtutorial にある BAR　の計算
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/01_theory.html
と同じものを使っています。
なおBAR では-9.3kJ /mol程度 が以下の同じ計算条件で得られたことを確認しています

version
gromacs-5.1-beta1
ermod-0.3.1     


最小化と平衡化を行ったあとに、
周期境界条件でNPT (パリネロラーマン, 300K),を行い
timestep=2fs, PME を使ってcutoff は 1.2 nm  (VDWも同じく1,2nm)
100fs ごとに 10000 このsampling を行った 溶液(水とメタン)と純溶媒 (水)に対して
gen_input を行い
それぞれsoln とrefs でermod
を行いました。
soln parameters_er
&ene_param
        boxshp = 1,
        cltype = 2,
        cmbrule = 1,
        elecut = 12,
        engdiv = 10,
        estype = 2,
        inptemp = 300,
        ljformat = 2,
        ljswitch = 1,
        lwljcut = 10,
        ms1max = 25,
        ms2max = 25,
        ms3max = 25,
        screen = 0.288243,
        sltspec = 1,
        slttype = 1,
        splodr = 6,
        upljcut = 12,
/

refs
&ene_param
        boxshp = 1,
        cltype = 2,
        cmbrule = 1,
        elecut = 12,
        engdiv = 5,
        estype = 2,
        inptemp = 300,
        ljformat = 2,
        ljswitch = 1,
        lwljcut = 10,
        ms1max = 25,
        ms2max = 25,
        ms3max = 25,
        screen = 0.288243,
        slttype = 3,
        splodr = 6,
        upljcut = 12,

/
で刻みは
&hist
      eclbin=2.0e-0, ecfbin=2.0e-2, ec0bin=2.0e-2, finfac=10.0e0,
      ecdmin=-3000.000000, ecfmns=-2.0e0, ecdcen=0.0e0, eccore=1.0e4,
      ecdmax=1.0e28, pecore=2000
ととりました。

するとそれぞれのuvrange.tt が kJ/mol で
soln
 species     minimum        maximum
    0       -5.08494        0.00001
    1      -12.22971      189.93318

refs
 species     minimum        maximum
    0       -5.01733        2.53609
    1     -843.86028       0.36868E+34

とご覧のように、水とメタンの相互作用の結果がオーダーレベルでずれています

N Matsubayashi and M. Nakahara, JCP 117 2002 では
solution で-0.2 kcal/mol から 10 kcal/mol 程度　()
refs の方は -0.2　から10e+11 程度となっています

refs に対しては、メタンを剛体で扱って計算もしましたが、そちらでもオーダーは変わりませんでした。

なお、他のパラメターは

ref directory
MDinfo
10000 1
596
3

Sltinfo
1 C -0.24 0.276144 0.35
2 H 0.06 0.12552 0.25
3 H 0.06 0.12552 0.25
4 H 0.06 0.12552 0.25
5 H 0.06 0.12552 0.25

Molprml
1 O -0.834 0.636386 0.315061
2 H 0.417 0 0
3 H 0.417 0 0

soln directory
MDinfo
10000 2
596
3

MolPrm1
1 O -0.834 0.636386 0.315061
2 H 0.417 0 0
3 H 0.417 0 0

確認のためにGROMSCS で水とメタンの相互作用のエネルギ−値をLJ と Coulomb だけ出力させますと
Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Coul-SR:SOL-CH4           0.0716258     0.0044   0.644746 -0.00518384  (kJ/mol)
LJ-SR:SOL-CH4              -13.1193      0.082    3.04473   0.171955  (kJ/mol)

最小値と最大値は
       min    max
LJ   -3      3
Coul  -19  7
でしたので、GROMACSとermod の間に大きなズレが見れます。

力場の主要部分は
topol.top

; Include forcefield parameters
#include "oplsaa.ff/forcefield.itp"

[ moleculetype ]
; Name            nrexcl
Other               3

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB      massB
; residue   1 CH4 rtp CH4  q  0.0
     1   opls_138      1    CH4      C      1      -0.24     12.011   ; qtot -0.24
     2   opls_140      1    CH4     H1      1       0.06      1.008   ; qtot -0.18
     3   opls_140      1    CH4     H2      1       0.06      1.008   ; qtot -0.12
     4   opls_140      1    CH4     H3      1       0.06      1.008   ; qtot -0.06



[ bonds ]
;  ai    aj funct            c0            c1            c2            c3
    1     2     1
    1     3     1
    1     4     1
    1     5     1

[ angles ]
;  ai    aj    ak funct            c0            c1            c2            c3
    2     1     3     1
    2     1     4     1
    2     1     5     1
    3     1     4     1
    3     1     5     1
    4     1     5     1

; Include water topology
#include "oplsaa.ff/tip3p.itp"

-------------topol.top ---
ermod とGROMACS では同じ値をつかっているように見えます

トラジェクトリに問題があるのかと、メタンのセンタリングもしましたが、
特に変化はありません。


ermod のインストールは
./configure --with-fft=fftw
make
make install
でGnu コンパイラを使っています
MAkefile の関連ありそうな部分は以下のようになっています
FCFLAGS = -DMPI -g -O2 -fdefault-real-8 -fdefault-double-8 -O3 -ffast-math -fno-finite-
math-only -march=native
BLAS_LIBS = -lblas
CC = gcc
CCDEPMODE = depmode=gcc3
CFLAGS = -g -O2
CPP = gcc -E
CPPFLAGS = -I/usr2/postdoc/rurano/apps/fftw-3.3.4/include:/usr2/postdoc/rurano/apps/boo
st_1_58_0/lib

また、添付されていました、

 ermod_GROMACS_example.tar
で計算した場合はちゃんと一致しています

自分で計算したトラジェクトリに対しては以上のようになっています。
何かわかることがあればよろしくお願いします
浦野

[Nobuyuki MATUBAYASI]

浦野さん

松林です。おっしゃる通り、おかしい結果です。

uvrange.ttのspecies 0 （周期境界条件における、メタンの自分のイメージとの相互作用）も

-5 というのは、（負に）大きすぎます。

それで、LINCSを使うと、対称性の良いメタンなんかは、構造がおかしくなることがあるようです。

soln系のMDと isolated solute のMD にて、メタンの構造を見てみてもらえますか？

さらに、もっと対称性の低い分子、

例えば、疎水性のものであれば、ブタンなんかで試してみたら、どうなるでしょうか？

[石塚良介]

浦野様:

松林研の石塚と申します。幾つか質問をさせて下さい。

1) ERmodのslvfeの溶媒和自由エネルギーの値は何kcal/molでしたか？

2) gromacs tutorialのbarは部分電荷がゼロのメタンを使用していますが、ERmodの計算には

部分電荷有りのメタンを使用されています。部分電荷有りのメタンで計算したBARの結果が有れば

教えて下さい。また、その時のcouple-lambda0とcouple-lambda1をどのように指定されているかも

教えて頂けると助かります。

ちなみに、お教え頂いたトポロジーで当方が計算したERmodの結果は、

soln uvrange.tt

 species     minimum        maximum

    0       -0.00006        0.00004

    1       -0.72955        5.42238

refs uvrange.tt  

species     minimum        maximum

    0       -0.00008        0.00005

    1       -0.85929       0.38053E+38

slvfeの結果（LJ correctionなし）で以下の通りでした。

 cumulative average & 95% error for solvation free energy

  1     2.5405

  2     2.5397     0.0016

  3     2.5406     0.0020

  4     2.5451     0.0091

  5     2.5426     0.0086

  6     2.5403     0.0084

  7     2.5398     0.0071

  8     2.5354     0.0108

  9     2.5345     0.0097

 10     2.5334     0.0090

2015年8月30日 8:00 uran ry <ryur...@gmail.com>:

--
You received this message because you are subscribed to the Google Groups "ermod-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ermod-users...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[uran ry]

松林様、石塚様
浦野です、返信ありがとう御座います

> uvrange.ttのspecies 0 （周期境界条件における、メタンの自分のイメージとの相互作用）も

> -5 というのは、（負に）大きすぎます。

なお、書き忘れましたが、
GROMACSに直接計算させた solution のトラジェクトリに対するメタンーメタンのエネルギーは次のようになっています

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------

Coul-SR:CH4-CH4          -0.0408043    3.5e-05 0.00153945 0.000218822  (kJ/mol)
LJ-SR:CH4-CH4                     0          0          0          0  (kJ/mol)

上のCoulomb の分布で　最小値と最大値は
Min.   :-0.05271 　Max.   :-0.03882
と、こちらは正しく出ているように見えます

> それで、LINCSを使うと、対称性の良いメタンなんかは、構造がおかしくなることがあるようです。

> soln系のMDと isolated solute のMD にて、メタンの構造を見てみてもらえますか？

前のAndleyさんへの投稿でそれが述べられていましたので、確認しましたが、特に問題は見られませんでした。
周期境界条件でのずれがありましたので、トラジェクトリをセンタリングしたものや、
flexible をやめてRigid 条件で計算もさせましたが、変わらずorder がずれていました

> さらに、もっと対称性の低い分子、

> 例えば、疎水性のものであれば、ブタンなんかで試してみたら、どうなるでしょうか？

試しに、親水性のEthanol 水系を試しましたが、　refsでは同じく大きなズレが見られました

水(TIP3P) -Ethanol 系でMDは上と同じ計算条件です

paramerters_er

&hist
 eclbin=2.0e-0, ecfbin=2.0e-2, ec0bin=2.0e-2, finfac=10.0e0,

      ecdmin=-2000.000000, ecfmns=-2.0e0, ecdcen=0.0e0, eccore=20.0e0,
      ecdmax=1.0e28, pecore=1000

refs uvrange.tt
species     minimum        maximum

    0       -2.06107        0.00330
    1    -1967.37450       0.47986E+28

soln uvrange.tt
 species     minimum        maximum

    0       -0.00822       -0.00348
    1       -7.21691        6.40139

比較のために、solution のtrajectory でGROMACS に計算させたエネルギーの値は

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------

Coul-14:Other-Other        -27.4339       0.21    2.97912    1.37081  (kJ/mol)
LJ-14:Other-Other          0.586377     0.0072   0.800584 -0.00278691  (kJ/mol)
最小値と最大値は
         Coul                          LJ
    Min.   :-34.03   Min.   :-0.455967 
   Max.   :-22.73   Max.   : 6.910465 
のようになっています

計算に使用しました力場は
topol.top

; Include forcefield parameters
#include "./oplsaa.ff/forcefield.itp"

[ moleculetype ]
; Name            nrexcl
Other               3

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB      massB

; residue   1 ETO rtp ETO  q  0.0
     1   opls_157      1    ETO      C      1      0.145     12.011   ; qtot 0.145
     2   opls_140      1    ETO     H1      1       0.06      1.008   ; qtot 0.205
     3   opls_140      1    ETO     H2      1       0.06      1.008   ; qtot 0.265
     4   opls_154      1    ETO     OH      1     -0.683    15.9994   ; qtot -0.418
     5   opls_155      1    ETO     HO      1      0.418      1.008   ; qtot 0
     6   opls_135      1    ETO     C1      1      -0.18     12.011   ; qtot -0.18
     7   opls_140      1    ETO     H3      1       0.06      1.008   ; qtot -0.12
     8   opls_140      1    ETO     H4      1       0.06      1.008   ; qtot -0.06
    9   opls_140      1    ETO     H5      1       0.06      1.008   ; qtot 0

[ bonds ]
;  ai    aj funct            c0            c1            c2            c3
    1     2     1
    1     3     1
    1     4     1

    1     6     1
    4     5     1
    6     7     1
    6     8     1
    6     9     1

[ pairs ]

;  ai    aj funct            c0            c1            c2            c3

    2     5     1
    2     7     1
    2     8     1
    2     9     1
   3     5     1
    3     7     1
    3     8     1
    3     9     1
    4     7     1
    4     8     1
    4     9     1
    5     6     1

[ angles ]
;  ai    aj    ak funct            c0            c1            c2            c3
    2     1     3     1
    2     1     4     1

    2     1     6     1
    3     1     4     1
    3     1     6     1
    4     1     6     1
    1     4     5     1
    1     6     7     1
    1     6     8     1
    1     6     9     1
    7     6     8     1
    7     6     9     1
    8     6     9     1

[ dihedrals ]
;  ai    aj    ak    al funct            c0            c1            c2            c3            c4            c5
    2     1     4     5     3
    3     1     4     5     3
    6     1     4     5     3
    2     1     6     7     3
    2     1     6     8     3
    2     1     6     9     3
    3     1     6     7     3
    3     1     6     8     3
    3     1     6     9     3
    4     1     6     7     3
    4     1     6     8     3
    4     1     6     9     3

> 1) ERmodのslvfeの溶媒和自由エネルギーの値は何kcal/molでしたか？

GROMACS ですので、 おそらく自由エネルギーはkJ/mol
で出ていると思っていたのですが、上のslvfe の結果は次のようになりました


  Self-energy of the solute   =        -0.2694  kcal/mol


 cumulative average & 95% error for solvation energy
  1    -3.3435
  2    -3.1208     0.4455
  3    -3.0336     0.3107
  4    -3.0673     0.2298
  5    -3.0736     0.1784
  6    -3.0237     0.1766
  7    -3.0599     0.1659
  8    -3.0405     0.1488
  9    -3.0408     0.1312
 10    -3.0577     0.1221

 group    solvation free energy     error          difference
   1             4.59696           1.03476           1.56900
   2             4.41262           1.01554           1.38465
   3             3.02797           0.69322           0.00000
   4             1.87001           0.53405          -1.15796
   5             5.79095           0.95527           2.76298
   6             0.35295           0.35901          -2.67502
   7             6.50885           0.89247           3.48088
   8            -0.68058           0.26003          -3.70855
   9             5.41319           0.50087           2.38522
  10             7.15351           0.89505           4.12555

最後に
警告として
 Warning: mesh error is    2.763 kcal/mol and is larger than the recommended value of    0.100     kcal/mol
がでています

> 2) gromacs tutorialのbarは部分電荷がゼロのメタンを使用していますが、ERmodの計算には

部分電荷有りのメタンを使用されています。部分電荷有りのメタンで計算したBARの結果が有れば

教えて下さい。また、その時のcouple-lambda0とcouple-lambda1をどのように指定されているかも

教えて頂けると助かります。

確認のために、これからやってみますが、
Tutorialによると
in good agreement with the value obtained by Shirts et al. of 2.24 ± 0.01 kcal mol^-1 (per Table II of the Supporting Information of that paper),
despite the fact that I used about half the number of λ vectors for the van der Waals transformation than the original authors did,
in addition to the other changes described previously.
上のCoulomb を無視して、Vdw のラムダの数を減らしてこのチュートリアルは実行されたが、この場合、
結果は2.24 kcal/mol のShirts の結果と十分な一致が見られる。
と言っています

私が行った結果は、今から、実行しますので、計算できましたらその部分については再度返信します

以上のようになっています
よろしくお願いします

 

[uran ry]

浦野です、
追加となります。



> 2) gromacs tutorialのbarは部分電荷がゼロのメタンを使用してい

> ますが、ERmodの計算には

> 部分電荷有りのメタンを使用されています。部分電荷有りのメタンで計算したBARの結果が有れば

> 教えて下さい。また、その時のcouple-lambda0とcouple-lambda1をどのように指定されているかも

> 教えて頂けると助かります。

につきまして、BARを用いて同じ条件で計算しましたが、
ラムダの指定は
; init_lambda_state        0    1    2    3    4    5    6    7    8    9    10   11   12   13   14   15   16   17   18   19   20
  vdw_lambdas              = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
  coul_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
で、free energy は
total      0 -     20,   DG -9.10 +/-  0.12 kJ/mol
となりました。悪くない精度です。
より上げるためには、事実上、coulomb 部分の貢献が0でしたので、vdw の lambda=0-10 を 15とか20こにすれば、BARの結果はより,2.24kcal/mol に近づくと思います

> 試しに、親水性のEthanol 水系を試しましたが、　refsでは同じく大きなズレが見られました
> 水(TIP3P) -Ethanol 系でMDは上と同じ計算条件です

また、この系の上の結果を用いた、sfvfe の結果を書いていませんでしたが、
  Number of the   1-th solvent  =         1000

  Self-energy of the solute   =        -0.0542  kcal/mol

 cumulative average & 95% error for solvation energy

  1   -20.0269
  2   -19.8783     0.2972
  3   -19.8517     0.1797
  4   -19.7853     0.1837
  5   -19.7742     0.1440
  6   -19.8197     0.1487
  7   -19.8073     0.1281
  8   -19.8529     0.1437
  9   -19.8579     0.1271
 10   -19.9255     0.1767

group    solvation free energy     error          difference

   1            -7.31934           0.08497          -0.23022
   2            -6.98349           0.06645           0.10563
   3            -7.08912           0.09454           0.00000
   4            -8.47278           0.10764          -1.38365
   5            -6.56492           0.11259           0.52420
   6            -6.39034           0.12406           0.69878
   7            -3.56741           0.07859           3.52171
   8           -10.48391           0.15241          -3.39479
   9            -2.93164           0.13438           4.15748
  10            -1.66828           0.12034           5.42085

 Warning: mesh error is    1.384 kcal/mol and is larger than the recommended value of    0.100     kcal/mol
となっています

[uran ry]

こんにちは
問題の解決はできていませんが、もう少し調べてみて、わかったことがいくつかありますので報告します、


１，　まず、soln のuvrange.tt につきましては、
用いる　.xtc をクラスター計算させる前に、元の.tpr から再作成したものを用いるとうまく行きました

水ーメタン系のsoln において
 species     minimum        maximum
    0       -0.00005        0.00000
    1       -0.68500        4.89109

となり、正常な値になりました。
水ーエタノール系 soln uvrange.tt
 species     minimum        maximum
    0       -0.00827       -0.00344
    1       -7.21683        6.40342

つまり、solution の系に関しては、おそらく計算するPC と解析するPC が異なっているために、トラジェクトリが上手に読み込めてないのが
原因であったと推察されます

２．　しかし、同じく、refs の方にトラジェクトリの変換したものに対しermod を実行しても、相変わらず -1000 のminimum がspecies で出ます。
これは、水ーエタノール、水ーメタン系の両方です

一方で奇妙なことに、水ーブタン系も同じ計算条件で計算しましたが、
こちらの場合はsoln, refs をトラジェクトリを変換して、ブタンをrigide扱った場合 （エタノール、メタンはrigid でもflexible でも上の異常値を観測しています）、

refs uvrange.tt
 species     minimum        maximum

    0       -0.00006        0.00002
    1       -1.38303       0.53133E+35

soln uvrange.tt
 species     minimum        maximum

    0       -0.00012        0.00002
    1       -1.53726        4.74509

 group    solvation free energy     error          difference

   1             3.69915           0.09210           0.02887
   2             3.65826           0.09244          -0.01202
   3             3.67028           0.08940           0.00000
   4             3.72432           0.07708           0.05404
   5             3.75773           0.07532           0.08745
   6             3.74979           0.08182           0.07951
   7             3.72658           0.06236           0.05630
   8             3.70665           0.06768           0.03637
   9             3.99370           0.06455           0.32342
  10             3.67632           0.06299           0.00604

と収束の良い値が観測されました、実験値を知りませんが、おそらくこのぐらいなのではないでしょうか？

結局のところ、GROMACS のトラジェクトリの読み込みミスが原因という可能性が減少しために、
ermod のエネルギー計算が異常を起こしているようなので、どの部分か調べてみるため、いろいろ
水ーエタノール系でいろいろパラメータを変えて計算してみたところ、
静電カットオフを0.1 にした場合

refs uvrange.tt
 species     minimum        maximum

    0        0.00000        0.00000
    1       -0.56241       0.28489E+28

上の石塚さんの結果に近づきました、ということでどうやら静電計算がおかしいということが判明しました。
なので、可能性として考えられるのは
1. トラジェクトリの読み込み時の座標値がおかしくて、距離計算がおかしい
2. 読み込みのパラメータを失敗している　（GROMACS からか、パラメータファイルから?）
3. 計算条件　ボックスサイズなどがおかしい

この内,１はsoln で正しくヨメ込めているので、同じ操作をしている以上、可能性は低い。
よって、2,3 が怪しいと思われるので、この辺を調べる方法や、なにか他に調べる要素などが有りましたらよろしくお願いします

浦野

[Shun Sakuraba]

浦野様
ermodの開発チームの、東大の桜庭と申します。

> １， まず、soln のuvrange.tt につきましては、
> 用いる .xtc をクラスター計算させる前に、元の.tpr から再作成したものを用いるとうまく行きました

> ...

> 一方で奇妙なことに、水ーブタン系も同じ計算条件で計算しましたが、
> こちらの場合はsoln, refs をトラジェクトリを変換して、ブタンをrigide扱った場合 （エタノール、メタンはrigid でもflexible
> でも上の異常値を観測しています）、

.xtcファイルは、本来シミュレーションから出力されたものを直接使うことが出来るはずですが、
何らかの変換を行っているのでしょうか。
それとも、gen_input等でトラジェクトリを指定することを、「変換」と呼んでいたりしますでしょうか。
もし、ファイルフォーマットの変換等を行っている場合には、
変換に用いたコマンドの詳細を教えていただけると幸いです。

よろしくお願いします。

--
Shun SAKURABA, Ph.D.
Postdoc @ Asai Lab, Grad School of Frontier Sciences, U Tokyo

[uran ry]

桜庭様
返信ありがとうございます
変換とは
gmx trjconv -s md.trr -o traj.xtc
または
gmx trjconv -s md.xtc -o traj.xtc
をつかって、xtc を該当マシンでxtcを再出力させるか、トラじぇくとりの形式を変換する(trr から xtc)
ことです、

（予想ですが、バイナリファイルですので、
おそらく、元のファイルだとヘッダーのバイト数が異なって、解析マシンで再出力させると一致するとかになっているのだと
思います、以前にVMD でほかのソフトのMDトラジェクトリを読んだときに同じ現象に出会ったことがあります）

よろしくおねがいします
浦野

[Shun Sakuraba]

浦野様

GROMACSのトラジェクトリはポータブルなフォーマットになっており、
エンディアンの異なる機種間でも全く同一のフォーマットになります。
http://www.gromacs.org/Documentation/File_Formats/Trajectory_File
ですので、ご指摘のような問題は本来起きないはずです。

> gmx trjconv -s md.trr -o traj.xtc

こちらは、
"-f" md.trr
が正しいコマンドとなりますでしょうか？
もしその他のオプションを付けてgmx trjconvを呼んでいましたら、
正確なコマンド列をペーストしていただけると幸いです。
また、その際、「このトラジェクトリを使ったらおかしくなったが、このトラジェクトリだと正しかった」などが分かると助かります。

よろしくお願いします。

[uran ry]

桜庭様

\

GROMACSのトラジェクトリはポータブルなフォーマットになっており、
エンディアンの異なる機種間でも全く同一のフォーマットになります。
http://www.gromacs.org/Documentation/File_Formats/Trajectory_File
ですので、ご指摘のような問題は本来起きないはずです。

 > gmx trjconv -s md.trr -o traj.xtc
こちらは、
"-f" md.trr
が正しいコマンドとなりますでしょうか？
もしその他のオプションを付けてgmx trjconvを呼んでいましたら、
正確なコマンド列をペーストしていただけると幸いです。

なるほど、そのようなことはないんですか？
確かに正確にコマンドをペーストしたほうがいいですね。
私のジョブスクリプトの実行部分は

 gmx trjconv -f ../solvent/Production_MD/md.xtc -s ../solvent/Production_MD/md.tpr -pbc mol -ur compact -center -o traj_solvent.xtc < ../1.txt
 gmx trjconv -f ../solute/Production_MD/md.xtc -s ../solute/Production_MD/md.tpr -pbc mol -ur compact -center -o traj_solute.xtc  < ../1.txt

gen_input --traj traj_solvent.xtc --log ../solvent/Production_MD/md.log --rigid  ../solution/NPT/npt.gro

cp ../parameters_er_refs ./parameters_er

 rm -f progress.tt&&rm -f flcuv.tt  &&mpirun -np 4 ermod
です。solution も指定するトラジェクトリ名が変わっているだけでまったく同じコマンドを打っています

> また、その際、「このトラジェクトリを使ったらおかしくなったが、このトラジェクトリだと正しかった」などが分かると助かります。

 私の使用しているクラスターマシンにおいては、
２種類のintel prossser と　CPU Architecture　の組み合わせがあり

2 six-core 2.5 GHz Intel Xeon E5-2640

sandybridge

2 six-core 3.07 GHz Intel Xeon X5675

nehalem

上が解析時のマシン、下がMD 時のマシンです
なるべく共通部のあるエタノールー水系とブタンー水系で比較しますと
refs にの計算において
両方solvent のトラジェクトリは同じものを使っております。
また、それぞれ溶質は剛体として取り扱い、とってくる溶質は
最小化、NVT, NPT で300K, 1bar に平衡化してProduction Run をはしる前の溶質をとってきています
VMD でそれぞれチェックしていますが、特に異常な構造（結合の切れなど）は確認されませんでした。

これに対して、上記のコマンドを実行してermod を走らせると以上のようなおかしな結果になりました
平衡後の構造も溶質系の計算も上のをくぐらせれば、soln も前のとうこうのようにuvrange.ttはうまくいっていますので、
この剛体としてとってきた構造がおかしいということもすくなく思います

一応parameters を乗せておます
ethanol -水の力場は前の投稿にのせました
refs のparameters_erは

&ene_param
        boxshp = 1,
        cltype = 2,
        cmbrule = 1,

        elecut = 0.1,

        engdiv = 5,
        estype = 2,
        inptemp = 300,
        ljformat = 2,
        ljswitch = 1,
        lwljcut = 10,

        ms1max = 28,
        ms2max = 28,
        ms3max = 28,

        screen = 0.288243,
        slttype = 3,
        splodr = 6,
        upljcut = 12,
/

&hist
 eclbin=2.0e-0, ecfbin=2.0e-2, ec0bin=2.0e-2, finfac=10.0e0,

      ecdmin=-5000.000000, ecfmns=-2.0e0, ecdcen=0.0e0, eccore=20.0e0,
      ecdmax=1.0e28, pecore=1000

MolPrml

1 O -0.834 0.636386 0.315061
2 H 0.417 0 0
3 H 0.417 0 0

MDinfo
10000 1
1000
3

SltInfo
1 C 0.145 0.276144 0.35 25.36 12.4 3.75
2 H 0.06 0.12552 0.25 26.39 12.13 3.51
3 H 0.06 0.12552 0.25 25.07 13.06 2.93
4 O -0.683 0.71128 0.312 25.28 13.07 4.94
5 H 0.418 0 0 24.64 13.77 4.94
6 C -0.18 0.276144 0.35 24.46 11.2 3.79
7 H 0.06 0.12552 0.25 23.5 11.39 4.26
8 H 0.06 0.12552 0.25 24.4 10.65 2.85
9 H 0.06 0.12552 0.25 24.91 10.52 4.51

水-Buthane 系では
refs
parameter_er

&ene_param
        boxshp = 1,
        cltype = 2,
        cmbrule = 1,
        elecut = 12,
        engdiv = 5,
        estype = 2,
        inptemp = 300,
        ljformat = 2,
        ljswitch = 1,
        lwljcut = 10,

        ms1max = 28,
        ms2max = 28,
        ms3max = 28,
        screen = 0.288243,
        slttype = 2,

        splodr = 6,
        upljcut = 12,
/

&hist
      eclbin=5.0e-2, ecfbin=2.0e-3, ec0bin=2.0e-4, finfac=10.0e0,
      ecdmin=-20.000000, ecfmns=-0.20e0, ecdcen=0.0e0, eccore=20.0e0,
      ecdmax=1.0e11, pecore=200
/

MDinfo

10000 1
1000
3

SltInfo

1 C -0.18 0.276144 0.35 23.26 16.07 27.06
2 H 0.06 0.12552 0.25 22.79 15.09 26.92
3 H 0.06 0.12552 0.25 24.18 16.22 26.5
4 H 0.06 0.12552 0.25 23.64 16.03 28.07
5 C -0.12 0.276144 0.35 22.23 17.22 26.91
6 H 0.06 0.12552 0.25 21.32 16.86 27.38
7 H 0.06 0.12552 0.25 22.06 17.38 25.85
8 C -0.12 0.276144 0.35 22.52 18.55 27.55
9 H 0.06 0.12552 0.25 22.46 18.44 28.63
10 H 0.06 0.12552 0.25 21.78 19.28 27.22
11 C -0.18 0.276144 0.35 23.87 19.04 27.09
12 H 0.06 0.12552 0.25 23.93 20.13 27.01
13 H 0.06 0.12552 0.25 24.16 18.48 26.21
14 H 0.06 0.12552 0.25 24.66 18.81 27.8

1 O -0.834 0.636386 0.315061
2 H 0.417 0 0
3 H 0.417 0 0

topol.top

; Include forcefield parameters
#include "./oplsaa.ff/forcefield.itp"

[ moleculetype ]
; Name            nrexcl
Other               3

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB      massB

; residue   1 CH3 rtp CH3  q  0.0
     1   opls_135      1    CH3      C      1      -0.18     12.011   ; qtot -0.18
     2   opls_140      1    CH3     H1      1       0.06      1.008   ; qtot -0.12
     3   opls_140      1    CH3     H2      1       0.06      1.008   ; qtot -0.06
     4   opls_140      1    CH3     H3      1       0.06      1.008   ; qtot 0
; residue   2 CH2 rtp CH2  q  0.0
     5   opls_136      2    CH2      C      2      -0.12     12.011   ; qtot -0.12
     6   opls_140      2    CH2     H1      2       0.06      1.008   ; qtot -0.06
     7   opls_140      2    CH2     H2      2       0.06      1.008   ; qtot 0
; residue   3 CH2 rtp CH2  q  0.0
     8   opls_136      3    CH2      C      3      -0.12     12.011   ; qtot -0.12
     9   opls_140      3    CH2     H1      3       0.06      1.008   ; qtot -0.06
    10   opls_140      3    CH2     H2      3       0.06      1.008   ; qtot 0
; residue   4 CH3 rtp CH3  q  0.0
    11   opls_135      4    CH3      C      4      -0.18     12.011   ; qtot -0.18
    12   opls_140      4    CH3     H1      4       0.06      1.008   ; qtot -0.12
    13   opls_140      4    CH3     H2      4       0.06      1.008   ; qtot -0.06
    14   opls_140      4    CH3     H3      4       0.06      1.008   ; qtot 0

[ bonds ]
;  ai    aj funct            c0            c1            c2            c3
    1     2     1
    1     3     1
    1     4     1

    1     5     1

    5     6     1
    5     7     1

    5     8     1
    8     9     1
    8    10     1
    8    11     1
   11    12     1
   11    13     1
   11    14     1

[ pairs ]
;  ai    aj funct            c0            c1            c2            c3

    1     9     1
    1    10     1
    1    11     1
    2     6     1

    2     7     1
    2     8     1

    3     6     1
    3     7     1
    3     8     1
    4     6     1

    4     7     1
    4     8     1

    5    12     1
    5    13     1
    5    14     1
    6     9     1
    6    10     1
    6    11     1
    7     9     1
    7    10     1
    7    11     1
    9    12     1
    9    13     1
    9    14     1
   10    12     1
   10    13     1
   10    14     1

[ angles ]
;  ai    aj    ak funct            c0            c1            c2            c3
    2     1     3     1
    2     1     4     1

    2     1     5     1
    3     1     4     1

    3     1     5     1
    4     1     5     1

    1     5     6     1
    1     5     7     1
    1     5     8     1
    6     5     7     1
    6     5     8     1
    7     5     8     1
    5     8     9     1
    5     8    10     1
    5     8    11     1
    9     8    10     1
    9     8    11     1
   10     8    11     1
    8    11    12     1
    8    11    13     1
    8    11    14     1
   12    11    13     1
   12    11    14     1
   13    11    14     1

[ dihedrals ]
;  ai    aj    ak    al funct            c0            c1            c2            c3            c4            c5

    2     1     5     6     3
    2     1     5     7     3
    2     1     5     8     3
    3     1     5     6     3
    3     1     5     7     3
    3     1     5     8     3
    4     1     5     6     3
    4     1     5     7     3
    4     1     5     8     3
    1     5     8     9     3
    1     5     8    10     3
    1     5     8    11     3
    6     5     8     9     3
    6     5     8    10     3
    6     5     8    11     3
    7     5     8     9     3
    7     5     8    10     3
    7     5     8    11     3
    5     8    11    12     3
    5     8    11    13     3
    5     8    11    14     3
    9     8    11    12     3
    9     8    11    13     3
    9     8    11    14     3
   10     8    11    12     3
   10     8    11    13     3
   10     8    11    14     3

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include water topology
#include "./oplsaa.ff/tip3p.itp"


[ molecules ]
; Compound        #mols
Other               1
SOL              1000



のようになっています
よろしくおねがいします

[Shun Sakuraba]

浦野様

返信ありがとうございます。確認したい点が2点あるのですが、

1.:

> gmx trjconv -f ../solvent/Production_MD/md.xtc -s
> ../solvent/Production_MD/md.tpr -pbc mol -ur compact -center -o
> traj_solvent.xtc < ../1.txt

> ..

> これに対して、上記のコマンドを実行してermod を走らせると以上のようなおかしな結果になりました
> 平衡後の構造も溶質系の計算も上のをくぐらせれば、soln も前のとうこうのようにuvrange.ttはうまくいっていますので、
> この剛体としてとってきた構造がおかしいということもすくなく思います

・「平衡後の構造」とは、純溶媒系の計算結果のトラジェクトリの事を指している、という理解で良いでしょうか？
・-pbc mol -ur compact -center を付けてトラジェクトリを変換すると、
　溶液系(soln)の場合では、付けなかった場合に比較して、「結果が改善された」
　という事で良いでしょうか？

2.:

> gmx trjconv -f ../solute/Production_MD/md.xtc -s
> ../solute/Production_MD/md.tpr -pbc mol -ur compact -center -o
> traj_solute.xtc < ../1.txt

Soluteは周期境界系では無いトラジェクトリであり、box情報を持たないはずだと思うのですが、
-pbc molで周期境界処理を行えていることから、異なる条件でシミュレーションを行っていることが予想されます。
Soluteを適切に（真空、非周期境界条件が正しくなります）シミュレーションしているかどうかの確認をお願いできますでしょうか。

よろしくお願いします。

> <http://ark.intel.com/products/64591/Intel-Xeon-Processor-E5-2640-15M-Cache-2_50-GHz-7_20-GTs-Intel-QPI>

>
>
> sandybridge
> 2 six-core 3.07 GHz Intel Xeon X5675

> <http://ark.intel.com/products/52577/Intel-Xeon-Processor-X5675-12M-Cache-3_06-GHz-6_40-GTs-Intel-QPI?q=x5675>
>
>
> nehalem上が解析時のマシン、下がMD 時のマシンです

[Nobuyuki MATUBAYASI]

> Soluteは周期境界系では無いトラジェクトリであり、box情報を持たないはずだと思うのですが、-pbc molで周期境界処理を行えていることから、異なる条件でシミュレーションを行っていることが予想されます。 
> Soluteを適切に（真空、非周期境界条件が正しくなります）シミュレーションしているかどうかの確認をお願いできますでしょうか。

松林です。少しだけ補足です。

isolated soluteを周期境界でやっても構いません。

isolated soluteは、分子内ポテンシャルだけでサンプリングしたいので、

PMEのように、unit cellの存在を前提とし、neutralizing backgroundとの相互作用が入るものはNGですが、

coulombtype = Cut-off として、cutoff lengthを事実上の無限大にすれば問題はありません。

isolated soluteの場合、最初の速度を０にしていなければ、周期境界を使わないと、

座標値がどんどん大きくなって、桁がおかしくなる可能性があります。

形式的に周期境界を使えば、この危険性は避けられます

[uran ry]

浦野です、返答ありがとうございます



 >   gmx trjconv -f ../solvent/Production_MD/md.

xtc -s
 > ../solvent/Production_MD/md.tpr -pbc mol -ur compact -center -o
 > traj_solvent.xtc < ../1.txt
 > ..
 > これに対して、上記のコマンドを実行してermod を走らせると以上のようなおかしな結果になりました
 > 平衡後の構造も溶質系の計算も上のをくぐらせれば、soln も前のとうこうのようにuvrange.ttはうまくいっていますので、
 > この剛体としてとってきた構造がおかしいということもすくなく思います
・「平衡後の構造」とは、純溶媒系の計算結果のトラジェクトリの事を指している、という理解で良いでしょうか？

ここでの平衡後の構造というのは、solution 系、solvent系に対して、最小化を行ったあと、NVT, NPT (200ps) を行って、300K, 1気圧に平衡化した状態のsolution の中のの溶質の構造を指しています。
計算に用いる構造はそれぞれその後の、Production_run の計算結果(1ns )を用いています。
rigid では１点構造を用いていますので、最小化と,NPT への平衡化が終わった状態の１点構造を用いたということです

 

・-pbc mol -ur compact -center を付けてトラジェクトリを変換すると、
　溶液系(soln)の場合では、付けなかった場合に比較して、「結果が改善された」
　という事で良いでしょうか？

はい、そのとおりです、水メタンにおいては結果は改善されました。


すいませんが、よく調べると、水ブタン系ですと solnにおいては
上のコマンドは不要でした。
trjconv 非実行、trjconv 実行後
どちらも同じ値が得られました
 species     minimum        maximum
    0       -0.00011        0.00001
    1       -1.53815        4.83950

水ーエタノール系 soln uvrange.tt
 species     minimum        maximum

    0       -0.00822       -0.00348
    1       -7.21691        6.40139

ただし、refs　においては水エタノール系、水メタンともにtrjconv を実行しても１桁ほど最小値が減少しましたが、
以前-2000 と異常な負の最小値をとっております

2.:
 >   gmx trjconv -f ../solute/Production_MD/md.xtc -s
 > ../solute/Production_MD/md.tpr -pbc mol -ur compact -center -o
 > traj_solute.xtc  < ../1.txt

Soluteは周期境界系では無いトラジェクトリであり、box情報を持たないはずだと思うのですが、
-pbc molで周期境界処理を行えていることから、異なる条件でシミュレーションを行っていることが予想されます。
Soluteを適切に（真空、非周期境界条件が正しくなります）シミュレーションしているかどうかの確認をお願いできますでしょうか。

VMDで見て確認したのと、上にも出しましたがGROMACS で出力させたエネルギーの値（全ポテンシャル、や溶質ー溶質相互作用、溶媒溶媒相互作用、溶質ー溶媒相互作用）を
確認した範囲では妥当な値が出力されています。
トラジェクトリを読み込んでermod で計算させると発生する問題のようです。

一つ質問なのですが、
さらに、念の為、見逃していてもいいように --flexible ではなく --rigid で計算をさせていますが、私の認識ではflexible では挿入する溶質の構造は solute のトラジェクトリから取ってくるのにたいし
rigid では 指定したgro のなから溶質の構造を撮ってきて、それだけをつかう。
という認識なのですが、この場合一点構造は最小化を行っていれば、万が一溶質の計算が失敗していても最小化の構造なので、精度は落ちてもermodの計算に異常な構造が使われることが
ないと思っていたのですが、なにか思い違いはあるのでしょうか？

> Soluteは周期境界系では無いトラジェクトリであり、box情報を持たないはずだと思うのですが、-pbc molで周期境界処理を行えていることから、異なる条件でシミュレーションを行っていることが予想されます。 

> Soluteを適切に（真空、非周期境界条件が正しくなります）シミュレーションしているかどうかの確認をお願いできますでしょうか。
parameters_er で指定する場所が無いようなのですが、GROMACS と同じBox size　などが使われているかなどを確認することはできるでしょうか？

> isolated soluteを周期境界でやっても構いません。

> isolated soluteは、分子内ポテンシャルだけでサンプリングしたいので、

> PMEのように、unit cellの存在を前提とし、neutralizing backgroundとの相互作用が入るものはNGですが、

> coulombtype = Cut-off として、cutoff lengthを事実上の無限大にすれば問題はありません。

> isolated soluteの場合、最初の速度を０にしていなければ、周期境界を使わないと、

> 座標値がどんどん大きくなって、桁がおかしくなる可能性があります。

> 形式的に周期境界を使えば、この危険性は避けられます

確かに溶質計算にPME が入っていますので、一度なくしてカットオフ最大で計算してみます

以上です
よろしくお願いします

2015年9月1日火曜日 6時49分20秒 UTC-4 Nobuyuki MATUBAYASI:

[uran ry]

浦野です、

水ーメタン系において
溶質単体の計算を周期境界条件をつかい、
rlist                    = 1.2
; Electrostatics
coulombtype              = Cut-off
rcoulomb                 = 1.2

そのトラジェクトリに対して
 gmx trjconv -f ../solute/Production_MD/md.trr -s ../solute/Production_MD/md

.tpr -pbc mol -ur compact -center -o traj_solute.xtc  < ../1.txt

と
gmx trjconv -f ../solvent/Production_MD/md.trr -s ../solvent/Production_MD/md.tpr -pbc mol -ur compact -center -o traj_solvent.xtc < ../1.txt
をしたあとのxtc ファイルに対して
gen_input --traj traj_solvent.xtc --log ../solvent/Production_MD/md.log --fl
exible  traj_solute.xtc
をしましたが
結果はやはり
refs uvrange.tt
 species     minimum        maximum
    0       -0.00001        0.00001
    1     -104.00627       0.27605E+28

最小値が-104 と相変わらず大きく出ています

それから、

もうひとつ気になった事項がありますので、
伝えます

メタン水のsoln の計算において
trjconv をかけずに直接Production Run のトラジェクトリからermodを計算する場合、

  energy of    2008.     for   1-th species
 The eccore parameter is too small; there should be no distribution at energy coordinate larger than eccore in the solution system

のようにeccore が小さすぎるというエラーがよく発生します。eccore をとても大きくとれば (1.0e+8 など)　計算はできるのですが、
不思議なのですが、eccore のあとでは locarithm meshの計算があるはずだから
eccore が大きなエラーがでても問題ないと思われるのですが、

これもなにか問題を示唆しているのでしょうか。

以上のようになっています
よろしくお願いします

浦野

[uran ry]

浦野です

もう少し調べてみたのですが、
GROMACS ４からGROMACS ５になった時に
すこし xtc ファイルのフォーマットが変わっているようなのですが、
どちらで計算されているのでしょうか？
私はGROMACS 5 で計算しているのですが、そちらもGROMACS5で計算されているのでしょうか

GROMACS Release note によると http://www.gromacs.org/About_Gromacs/Release_Notes/Versions_5.0.x
Fixed selection of xtc groups not starting at index 0 #1561.

とのことです。
確認よろしくお願いします

[Nobuyuki MATUBAYASI]

浦野さん

今回の問題の解決には至らない部分ですが、いくつかの点について回答します

> 一方で奇妙なことに、水ーブタン系も同じ計算条件で計算しましたが、

> こちらの場合はsoln, refs をトラジェクトリを変換して、ブタンをrigide扱った場合 

> （エタノール、メタンはrigid でもflexible でも上の異常値を観測しています）、

・・・

> と収束の良い値が観測されました、実験値を知りませんが、おそらくこのぐらいなのではないでしょうか？

ブタンに関する計算の点、返事が遅れましたが、その通りです。

実験値もそのようなものですし、我々の方での計算値も同じような値です

> rigid では 指定したgro のなから溶質の構造を撮ってきて、それだけをつかう。

> という認識なのですが、この場合一点構造は最小化を行っていれば、

> 万が一溶質の計算が失敗していても最小化の構造なので、精度は落ちてもermodの計算に

> 異常な構造が使われることがないと思っていたのですが、なにか思い違いはあるのでしょうか？

はい、それで結構と思います。

実際には、壊れた構造とかでなければ、最小化も不要です。

計算するのは、溶質ー溶媒相互作用だけです。

溶質の構造がおかしければ、溶質の分子内ポテンシャルがおかしなことになるのですが、

そのことは、溶質ー溶媒相互作用がヘンになることを意味するわけでは無いので

（なので、遷移状態にある溶質でもOKです）

> trjconv をかけずに直接Production Run のトラジェクトリからermodを計算する場合、

>  energy of    2008.     for   1-th species

> The eccore parameter is too small; there should be no distribution at energy coordinate larger than eccore in the solution system

> のようにeccore が小さすぎるというエラーがよく発生します。

> eccore をとても大きくとれば (1.0e+8 など)　計算はできるのですが、

> 不思議なのですが、eccore のあとでは locarithm meshの計算があるはずだから

> eccore が大きなエラーがでても問題ないと思われるのですが、これもなにか問題を示唆しているのでしょうか。

エネルギーがeccoreより大きいと、soln系では分布関数（engsln.XX）が０になるように、

eccoreを設定しています。つまり、eccoreより大きいところが、coreの領域というわけです。

eccoreの値は、ある程度以上大きければなんでも良いのですが、

溶質ー溶媒のペアエネルギー値は、soln系では、普通は、20 kcal/mol程度以下のはずなので、

さすがに、2000 kcal/molというのは、何かおかしいということになります

> GROMACS ４からGROMACS ５になった時に

> すこし xtc ファイルのフォーマットが変わっているようなのですが、どちらで計算されているのでしょうか？

> 私はGROMACS 5 で計算しているのですが、そちらもGROMACS5で計算されているのでしょうか

こちらでは、gromacs 4の人もいるし、5の人もいて、両方が同じERmodを使っていますが、

特に問題は出ていません

松林

[Shun Sakuraba]

浦野様

桜庭です。
現象と、各種パラメータの値から考えますと、
現在の系のbox sizeはz方向がcutoffの2倍ぎりぎりの値となっており、
これによって何らかの不具合を起こしている可能性が高いのでは無いかと考えております。
(GROMACS4→5の当該チケットでは、xtcファイルについての特定の処理のバグを修正しただけで、ファイルフォーマットの変更はありません)

ですが、これ以上は原因がはっきりしませんので、
問題の追跡のため、以下3点をお願いできますでしょうか。
1. parameters_erの &ene_param 内に、一行
block_threshold = 1.0,
と加えた場合、問題は改善しますでしょうか。
2. GROMACS 5.1 betaを利用中とのことですが、なにぶんβ版ですし、GROMACSはバージョンアップ時に良く問題を起こします。
　ついては、5.0系でももし計算可能でしたら、大変お手数ですが、5.0系で問題が起きないかどうかを確認願えませんでしょうか。
3. 上記で解決が困難な場合、
solution, reference系のmdrun時の入力ファイルを桜庭まで送付いただけませんでしょうか。
　.tprファイルのみを送っていいただければ充分です。

よろしくお願いします。

[uran ry]

桜庭様、松林様、石塚様
返信ありがとうございます、

石塚さんから送られてきたinput ファイル群を元に（こちらではuvrange.ttで正しそうな

結果を再現しました、ただし自由エネルギーは負の値が出ていますので、メタン系としては
おかしいと思えます、slvfe の自由エネルギーの出力値を教えていだたけますか？）
条件を変えて、順番に試していったところ、uvrange.tt に異常値が発生する条件らしきものがわかりました


まず試しにGROMACS4 で同じ条件で計算したところ、同じ異常な結果が発生しましたので、ご指摘のようにGROMACSのバージョンが原因ではありませんでした。

水ーメタン系において
1. 石塚さんから送られてきたファイルではuvrange.ttに異常がないことを確認
soln uvrange.tt
 species     minimum        maximum
    0       -0.00007        0.00004
    1       -0.67294        4.44217

slvfe.log
refs
 species     minimum        maximum
    0       -0.00006        0.00005
    1       -0.56027       0.66506E+40


上のにたいするslfve の結果は

  Number of the   1-th solvent  =         1000

  Self-energy of the solute   =        -0.0000  kcal/mol

 cumulative average & 95% error for solvation energy

  1    -3.1653
  2    -3.1680     0.0053
  3    -3.1846     0.0333
  4    -3.1705     0.0366
  5    -3.1662     0.0296
  6    -3.1667     0.0242
  7    -3.1736     0.0247
  8    -3.1678     0.0243
  9    -3.1654     0.0220
 10    -3.1653     0.0197

 group    solvation free energy     error          difference

   1           -51.47256           1.98368          -3.80705
   2           -44.85675           1.84581           2.80876
   3           -47.66552           1.99794           0.00000
   4           -43.84219           1.78652           3.82332
   5           -45.02757           1.84800           2.63794
   6           -44.02303           1.79894           3.64248
   7           -19.97499           1.22779          27.69053
   8           -43.54243           1.77218           4.12309
   9           -12.79311           0.98585          34.87241
  10           -43.47141           1.76887           4.19411

自由エネルギーが-3.16 kcal/mol と 2.21 kcal /mol から全く異なる値を得ていますが　（--rigid にしたせいでしょうか？）
とりあえず、この結果を参照として、徐々に元の計算条件に戻していきました

2. 最初に力場ファイルを送られてきたものから、自身の使用していたものに変えたところ、参照系に
近い結果でした

refs uvrange.tt
 species     minimum        maximum

    0       -0.00006        0.00005
    1       -0.56027       0.66506E+40

slvfe.log

  Self-energy of the solute   =        -0.0000  kcal/mol

 cumulative average & 95% error for solvation energy

  1    -3.0873
  2    -3.1222     0.0698
  3    -3.1354     0.0482
  4    -3.1300     0.0357
  5    -3.1237     0.0304
  6    -3.1283     0.0265
  7    -3.1248     0.0234
  8    -3.1289     0.0218
  9    -3.1348     0.0226
 10    -3.1340     0.0203

 group    solvation free energy     error          difference

   1           -52.55578           3.17061          -3.72545
   2           -45.93310           2.94517           2.89723
   3           -48.83033           3.16837           0.00000
   4           -44.84112           2.83706           3.98921
   5           -46.06060           2.95555           2.76973
   6           -45.07355           2.87650           3.75678
   7           -20.67760           1.85427          28.15273
   8           -44.54125           2.83038           4.28908
   9           -13.30199           1.47440          35.52834
  10           -44.45226           2.81440           4.37806


 Warning: mesh error is    3.989 kcal/mol and is larger than the recommended value of    0.100     kcal/mol

3.  ２からさらに　input ファイル(水1000 個)を私がつかっていた元のpdbファイルに戻しました(水　600)。すると異常な結果が発生しました
結果として溶媒のサイズが小さすぎるときに()、周期境界条件をつかってERMOD を計算すると起きるようです(PME のon off は関係なかったです)
ですので
> 現象と、各種パラメータの値から考えますと、
> 現在の系のbox sizeはz方向がcutoffの2倍ぎりぎりの値となって
> これによって何らかの不具合を起こしている可能性が高いのでは無いかと考えております。
が可能性が高いようです。
さらに、このプロダクションRUN のトラジェクトリにおいて、ボックスサイズがかなり大きく変化しているので、
こちらが原因として有り得そうですが、ボックスサイズの可変性などは倍などに変化しても大丈夫なのでしょうか？

それから、 uvrange.tt に関しては解決したようですが、結果として出力される溶媒和自由エネルギーは実験値の2.21kcal/mol から大きく外れているようですが、
こちらは何を試せばいいでしょうか？
現在のparameter_fe は
&fevars
inptemp=300.000000
/
のみです

この実験値からのズレは再計算した水エタノール系(水1000 こ)では確認されませんでした、パラメータは上と同じです

soln uvrange.tt

 species     minimum        maximum
    0       -0.00866       -0.00356
    1       -7.01317        5.79841

aveuv.tt
    1      -20.05247
    2      -19.86766
    3      -19.86083
    4      -19.82894
    5      -19.77203
    6      -19.72598
    7      -20.06918
    8      -20.18950
    9      -19.66349
   10      -19.66214

refs uvrange.tt
 species     minimum        maximum

    0       -0.00743       -0.00469
    1       -9.08044       0.25874E+38


slvfe.log

  Number of the   1-th solvent  =         1000

  Self-energy of the solute   =        -0.0060  kcal/mol

 cumulative average & 95% error for solvation energy

  1   -20.0585
  2   -19.9661     0.1848
  3   -19.9330     0.1255
  4   -19.9085     0.1014
  5   -19.8824     0.0943
  6   -19.8573     0.0919
  7   -19.8885     0.0995
  8   -19.9269     0.1154
  9   -19.8983     0.1168
 10   -19.8753     0.1141

group    solvation free energy     error          difference

   1            -4.72146           0.05502           0.01251
   2            -4.72760           0.04875           0.00637
   3            -4.73396           0.05096           0.00000
   4            -4.72334           0.04967           0.01062
   5            -4.73544           0.05247          -0.00148
   6            -4.71643           0.04903           0.01754
   7            -4.72962           0.05305           0.00435
   8            -4.75476           0.05038          -0.02080
   9            -4.71478           0.05018           0.01918
  10            -4.73107           0.05254           0.00290

実験値は-4.9 kcal/mol ですので十分近い一致が見れました

> 問題の追跡のため、以下3点をお願いできますでしょうか。
> 1. parameters_erの &ene_param 内に、一行
> block_threshold = 1.0,
> と加えた場合、問題は改善しますでしょうか。

試してみましたところ
segmentation fault
がでて止まりました

[:84251] *** Process received signal ***
[:84251] Signal: Segmentation fault (11)
[:84251] Signal code: Address not mapped (1)
[:84251] Failing at address: 0x19c6b90
[:84248] [ 0] /lib64/libpthread.so.0(+0xf710) [0x2b8037bde710]
84248] [ 1] ermod(__engproc_MOD_enginit+0xa7a) [0x414b5a]
:84248] [ 2] ermod(enganal_init_+0x12) [0x4077b2]
:84248] [ 3] ermod(MAIN__+0x59) [0x4073c9]
:84248] [ 4] ermod(main+0x2a) [0x42881a]
:84248] [ 5] /lib64/libc.so.6(__libc_start_main+0xfd) [0x2b8037e0ad5d]
:84248] [ 6] ermod() [0x4072a9]

> 2. GROMACS 5.1 betaを利用中とのことですが、なにぶんβ版ですし、GROMACSはバージョンアップ時に良く問題を起こします。
>　ついては、5.0系でももし計算可能でしたら、大変お手数ですが、5.0系で問題が起きないかどうかを確認願えませんでしょうか。

4.6系で起きたので、5.0系でも起きるとおもいます


最後に水メタン系においても、上ではsolute をrigid に扱っていたため、負の自由エネルギーがでていましたが、
GAS フェイズのトラジェクトリをつかって、flexible にして、 ermod, slvfe をしたところ

  Number of the   1-th solvent  =         1000

  Self-energy of the solute   =        -0.0000  kcal/mol

 cumulative average & 95% error for solvation energy

  1    -3.1654
  2    -3.1680     0.0053
  3    -3.1846     0.0333
  4    -3.1705     0.0366
  5    -3.1662     0.0296
  6    -3.1667     0.0242
  7    -3.1736     0.0247
  8    -3.1678     0.0243
  9    -3.1654     0.0220
 10    -3.1653     0.0197

 group    solvation free energy     error          difference

   1             2.55686           0.02803           0.01670
   2             2.53450           0.02631          -0.00566
   3             2.54016           0.02632           0.00000
   4             2.55182           0.02340           0.01166
   5             2.52452           0.02505          -0.01564
   6             2.58178           0.02361           0.04161
   7             2.51477           0.02419          -0.02539
   8             2.50745           0.02312          -0.03271
   9             2.62903           0.02009           0.08887
  10             2.53969           0.02304          -0.00047


BARで得た、2.2 kcal/mol からは 0.3kcal /mol ほど離れていますが、温度や密度を考慮すれば良い一致が見られました。
ですので、溶媒和の自由エネルギー計算としては無事成功したようにおもいます
皆様、ありがとうございました

最後に、２つ疑問なのですが、
1. 上の計算の失敗に関して、周期境界条件を使い、セルサイズが小さいときにおきているようなのですが、
もし、圧力が高かったり、低かったりして、セルサイズが変更しやすい場合も上と同じような計算の失敗が起こる可能性があるのでしょうか？

2. 上で見たようにメタン系において、rigid と flexible の場合に slvfe の結果に大きな違いが見られたのですが、
メタン程度のサイズが自由エネルギーの正負がひっくり返るような大きな
違いが起きるのはよくあることなのでしょうか？　（rigid の場合に、自由エネルギーの精度が悪くなるのは予想していましたが、
ここまで大きな違いが出るとはおもっていなかったのですが）

よろしくおねがいします

浦野

[Nobuyuki MATUBAYASI]

浦野さん

松林です。

> 水ーメタン系において
> 1. 石塚さんから送られてきたファイルではuvrange.ttに異常がないことを確認

> 自由エネルギーが-3.16 kcal/mol と 2.21 kcal /mol から全く異なる値を得ていますが　

> （--rigid にしたせいでしょうか？）

--rigidとはおそらく関係はないと思います。

> メタン系としておかしいと思えます、

> slvfe の自由エネルギーの出力値を教えていだたけますか？

おっしゃる通り、おかしい結果ではあります。

問題がなければ、2.5 kcal/molぐらいになるはずです

> 3.  ２からさらに　input ファイル(水1000 個)を私がつかっていた元のpdbファイルに

> 戻しました(水　600)。すると異常な結果が発生しました結果として

> 溶媒のサイズが小さすぎるときに()、周期境界条件をつかってERMOD を計算すると

> 起きるようです

実は、そうでもなくて、例えば、私は、メタンやエタノールで、テストの目的で、

周期境界条件で、水の分子数を１００まで減らしたことがありませんが、問題はありませんでした。

ちなみに、周期境界を使わずに、

dropletでやった場合は、溶媒の個数を１にしても、まともな数字が得られます

> このプロダクションRUN のトラジェクトリにおいて、

> ボックスサイズがかなり大きく変化しているので、こちらが原因として有り得そうですが、

> ボックスサイズの可変性などは倍などに変化しても大丈夫なのでしょうか？

ボックスサイズが倍も変わるのは、MD自体の問題だと思います。

系の体積をVとすると、体積変化の相対幅は１／sqrt(V)の程度であるべきで、

ある程度の大きさの系であれば、Vの最初の桁の話でありません

> 1. 上の計算の失敗に関して、周期境界条件を使い、

> セルサイズが小さいときにおきているようなのですが、

> もし、圧力が高かったり、低かったりして、セルサイズが変更しやすい場合も

> 上と同じような計算の失敗が起こる可能性があるのでしょうか？

計算が熱力学の範疇に入って、体積変化の相対幅が１／sqrt(V)の程度であれば、

大丈夫と思います。

ちなみに、臨界点の近くにきて、cell sizeが大いに揺らぐ場合でも、

体積変化の相対幅が１／sqrt(V)となるように、Vをちゃんと大きく取れば問題はでない、

というのがこれまでの経験です

> 2. 上で見たようにメタン系において、rigid と flexible の場合に s

> lvfe の結果に大きな違いが見られたのですが、
> メタン程度のサイズが自由エネルギーの正負がひっくり返るような大きな違いが

> 起きるのはよくあることなのでしょうか？

今回が初めてです。

メタン程度の小さな分子であれば、rigid でも flexible でも、

溶媒和自由エネルギーの値は、それほど変わりません

[uran ry]

松林様
返答ありがとうございます、

2015年9月3日木曜日 17時00分30秒 UTC-4 Nobuyuki MATUBAYASI:

浦野さん

松林です。

> 水ーメタン系において
> 1. 石塚さんから送られてきたファイルではuvrange.ttに異常がないことを確認

> 自由エネルギーが-3.16 kcal/mol と 2.21 kcal /mol から全く異なる値を得ていますが　

> （--rigid にしたせいでしょうか？）

--rigidとはおそらく関係はないと思います。

> メタン系としておかしいと思えます、

> slvfe の自由エネルギーの出力値を教えていだたけますか？

おっしゃる通り、おかしい結果ではあります。

問題がなければ、2.5 kcal/molぐらいになるはずです

> 3.  ２からさらに　input ファイル(水1000 個)を私がつかっていた元のpdbファイルに

> 戻しました(水　600)。すると異常な結果が発生しました結果として

> 溶媒のサイズが小さすぎるときに()、周期境界条件をつかってERMOD を計算すると

> 起きるようです

実は、そうでもなくて、例えば、私は、メタンやエタノールで、テストの目的で、

周期境界条件で、水の分子数を１００まで減らしたことがありませんが、問題はありませんでした。

ちなみに、周期境界を使わずに、

dropletでやった場合は、溶媒の個数を１にしても、まともな数字が得られます

確認のために、
ermod の計算時に周期境界条件を使用しないとして
boxshp=0
を用いて、ermod を実行することができませんでしたので、私のほうでは確認できませんでした。
となるとやはり、トラジェクトリの間にボックスサイズの可変が大きいというのが原因として高いそうです。
 

> このプロダクションRUN のトラジェクトリにおいて、

> ボックスサイズがかなり大きく変化しているので、こちらが原因として有り得そうですが、

> ボックスサイズの可変性などは倍などに変化しても大丈夫なのでしょうか？

ボックスサイズが倍も変わるのは、MD自体の問題だと思います。

系の体積をVとすると、体積変化の相対幅は１／sqrt(V)の程度であるべきで、

ある程度の大きさの系であれば、Vの最初の桁の話でありません

すみません、この点に関しては、私の聞き方が悪かったようです。
私の質問の意図としましては
MDで体積サイズが非常に大きく変化するのは確かに計算条件としてよろしくないのですが、
計算時にボックスサイズの変化が起きても、その情報は
GROMACS のトラジェクトリ情報に含まれています（実際GROMACS自体でエネルギー計算をする場合には
それをつかってまともなエネルギー値を出しています）。
その一方で、ERNOD をつかって周期境界条件の計算をするときには、
そのボックスサイズの変更の情報もGROMACSと同じように使えるはずなのですが、
そのようなボックスサイズの変更が行われていないために、本来のカットされるはずのイメージセルのぶんなどが余計に計算されて
uvrange.tt で以上に低い値を出しているのではないでしょうか？
とうことで、ボックスサイズが大きく変化しても大丈夫なのでしょうか？と尋ねました
 

> 1. 上の計算の失敗に関して、周期境界条件を使い、

> セルサイズが小さいときにおきているようなのですが、

> もし、圧力が高かったり、低かったりして、セルサイズが変更しやすい場合も

> 上と同じような計算の失敗が起こる可能性があるのでしょうか？

計算が熱力学の範疇に入って、体積変化の相対幅が１／sqrt(V)の程度であれば、

大丈夫と思います。

ちなみに、臨界点の近くにきて、cell sizeが大いに揺らぐ場合でも、

体積変化の相対幅が１／sqrt(V)となるように、Vをちゃんと大きく取れば問題はでない、

というのがこれまでの経験です

なるほど、Vをおおきくとる場合におきる密度の違いなどはどうなるのだろうか、
とか少し疑問に思いますが、今のところ私はこれに関係する計算は行わないので
とりあえず大丈夫であると理解します。
 

> 2. 上で見たようにメタン系において、rigid と flexible の場合に s

> lvfe の結果に大きな違いが見られたのですが、
> メタン程度のサイズが自由エネルギーの正負がひっくり返るような大きな違いが

> 起きるのはよくあることなのでしょうか？

今回が初めてです。

メタン程度の小さな分子であれば、rigid でも flexible でも、

溶媒和自由エネルギーの値は、それほど変わりません

そうですよね、uvrange.tt はsoln , refs ともに妥当な値を出しており、
refs のuvrange.tt はrigid , solid でもほとんど変わらなかったのに
上記のような結果になったのはやはりおかしいですよね。そもそも大きな負の値(-0.1 以下)が殆ど無いにもかかわらず
積分の結果が大きな負になるというのが不自然ですよね。
mesh error がやたら大きいのが原因なのだろうかと思いましたが、それも変ですし。
あとで、もう一度計算しなおして見ることにします

返信ありがとうございました
浦野

[石塚良介]

浦野様:

石塚です。１点だけ確認させて下さい。

私がお送りしたmdpファイル、トポロジーファイルを、

何の変更もせずに計算した結果は2.5 kcal/molを

再現できなかったということでしょうか？

2015年9月4日 8:30 uran ry <ryur...@gmail.com>:

[Nobuyuki MATUBAYASI]

浦野さん

松林です。一点だけ

> MDで体積サイズが非常に大きく変化するのは確かに計算条件としてよろしくないのですが、

> 計算時にボックスサイズの変化が起きても、その情報はGROMACS のトラジェクトリ情報に

> 含まれています（実際GROMACS自体でエネルギー計算をする場合には
> それをつかってまともなエネルギー値を出しています）。
> その一方で、ERNOD をつかって周期境界条件の計算をするときには、
> そのボックスサイズの変更の情報もGROMACSと同じように使えるはずなのですが、
> そのようなボックスサイズの変更が行われていないために、

> 本来のカットされるはずのイメージセルのぶんなどが余計に計算されて
> uvrange.tt で以上に低い値を出しているのではないでしょうか？
> とうことで、ボックスサイズが大きく変化しても大丈夫なのでしょうか？と尋ねました

そういうことでしたか。

ERmodの方でも問題は無いと思っておりましたが、

cutoff長がcell lengthをはみ出る場合は、ちょっと良く分かりません。

桜庭さんが

>>>> 現在の系のbox sizeはz方向がcutoffの2倍ぎりぎりの値となっており、 
>>>> これによって何らかの不具合を起こしている可能性が高いのでは無いかと考えております。

と言っていましたが、この点と関係するかもしれませんね。

ともあれ、上の件は誤解していましたが、現状では、良く分からないです。すみません
 

[uran ry]

浦野です

> 石塚です。１点だけ確認させて下さい。

> 私がお送りしたmdpファイル、トポロジーファイルを、

> 何の変更もせずに計算した結果は2.5 kcal/molを

> 再現できなかったということでしょうか？

solution, solvent, solute をすべて計算した場合は再現しましたが、
solution, solventは同じトラジェクトリをつかって
solute の計算をせずに、
--rigid system.gro をつかって計算した場合に
上のような結果になりました

それから、
2015年9月4日金曜日 6時30分15秒 UTC-4 Nobuyuki MATUBAYASI:

わかりました、丁寧にありがとうございます。
もしかして、エネルギーの増幅具合も20倍程度でしたので、ちょうどイメージセルの数ぐらいでしたので
トラジェクトリのボックスサイズの変更をermod においては行っていないのかなと思いましたが、
私のようにテスト計算だから小さな系をつかって計算してみよう、などとならない限り
本計算においては、多分問題にならないのでしょう。


浦野
 

[Nobuyuki MATUBAYASI]

浦野さん

松林です。しつこくて、すみません

> トラジェクトリのボックスサイズの変更をermod においては

> 行っていないのかなと思いましたが、

> 私のようにテスト計算だから小さな系をつかって計算してみよう、などとならない限り
> 本計算においては、多分問題にならないのでしょう

ボックスサイズの変更は、もちろん、ermodの中で行っております

（トラジェクトリファイルの中の、cell sizeの情報を読込んでいる）。

また、今回の浦野さんの計算は、水分子の個数が６００という記憶ですが、

メタンやエタノールであれば、全く小さい系ではなく、十分な大きさと思います。

これまでの経験では、もっと小さくても（１００や５０個の水）でも大丈夫でした。

一応、蛇足まで

[石塚良介]

浦野様:

石塚です。

>solution, solvent, solute をすべて計算した場合は再現しましたが、
>solution, solventは同じトラジェクトリをつかって
>solute の計算をせずに、
>--rigid system.gro をつかって計算した場合に
>上のような結果になりました

solutionは同じものを使わない方が良いと思います。

freezegrpsでCH4を固定したsolution MDと、

--rigidを指定したrefsから計算したslvfeは

以下の結果になりました。

ちなみに、--rigidではメタン分子１個だけが

記載されているgroファイルを指定しました。

 cumulative average & 95% error for solvation free energy

  1     2.4891

  2     2.4856     0.0070

  3     2.4965     0.0222

  4     2.4931     0.0171

  5     2.5003     0.0195

  6     2.4954     0.0187

  7     2.4967     0.0160

  8     2.4925     0.0161

  9     2.4900     0.0151

 10     2.4863     0.0154

確か、--rigidは結構注意が必要なオプションだった気がします。

2015年9月5日 17:37 Nobuyuki MATUBAYASI <nobu...@scl.kyoto-u.ac.jp>:

--

[uran ry]

松林様、石塚様
浦野です、返信ありがとうございます、

> ボックスサイズの変更は、もちろん、ermodの中で行っております

>（トラジェクトリファイルの中の、cell sizeの情報を読込んでいる）。

> また、今回の浦野さんの計算は、水分子の個数が６００という記憶ですが、

> メタンやエタノールであれば、全く小さい系ではなく、十分な大きさと思います。

> これまでの経験では、もっと小さくても（１００や５０個の水）でも大丈夫でした。

では、カットオフとボックスサイズがほとんど同じというのが問題だったのかもしれないですね。
本計算ではそうおきない気がしますが、注意することにします

> solutionは同じものを使わない方が良いと思います。

> freezegrpsでCH4を固定したsolution MDと、

> --rigidを指定したrefsから計算したslvfeは

> 以下の結果になりました。

> ちなみに、--rigidではメタン分子１個だけが

> 記載されているgroファイルを指定しました。

> 確か、--rigidは結構注意が必要なオプションだった気がします。
サイズの変更が行われており、上のことより指定rigid を変更しましたら
wiki のexmaple において、solution をrigid の対象に選んでいましたので、同じようにしました。
一応、CH4.gro を --rigid の対象に選んでも計算しましたが値は上と変わりませんでした(slvfe ~ -50kcal/mol)。
soln の結果は同じものをつかっているのですが.。
rigidの場合だと、なにかが起こっているのはまちがいないのですが、やはりよくわからないですね。
uvrange.tt, weight_refs の結果にはそこまで大きな違いがないのですが。
とりあえず, ethanol では問題がなかったですので、メタンのrigid だと注意が必要ということなのでしょう。

rigid オプションを使う場合は結果の確認に気をつけることにします

ありがとうございました
浦野

[Shun Sakuraba]

浦野様

もしよければ、問題の起こるケースでの.tprファイル並びに.top/.itpファイルを桜庭まで送付していただいて、問題の再現に協力していただけると幸いです。
よろしくお願いします。

[石塚良介]

浦野様、櫻庭さん:

石塚です。

--rigidでrefsを計算し、soln MDをフレキシブルメタン
（結合長のみLINCSで固定）で計算した場合は、

 cumulative average & 95% error for solvation free energy

  1   -49.0573
  2   -48.8966     0.3214
  3   -48.6574     0.5131
  4   -48.6506     0.3631
  5   -48.3512     0.6616
  6   -48.1842     0.6352
  7   -48.3685     0.6512
  8   -48.2921     0.5842
  9   -48.0761     0.6724
 10   -48.1832     0.6385

でしたので、浦野さんの結果（slvfe ~ -50kcal/mol）を再現しています。

freezegrpsでメタンを固定したMDを用いて計算したsolnを用いると、
前回報告した2.5 kcal/molになります。

2015年9月8日 12:47 Shun Sakuraba <shun.s...@gmail.com>:

[Shun Sakuraba]

浦野様

> では、カットオフとボックスサイズがほとんど同じというのが問題だったのかもしれないですね。
> 本計算ではそうおきない気がしますが、注意することにします
本件の問題はこの点ではなく、"--rigid"について誤った使い方を行っているのが原因と思われます。

rigid を使う場合、チュートリアルにも明示してありますが、
Solution系のシミュレーションを、分子が剛体であるように扱って実行する必要があります。
（通常、GROMACSではこれは[ constraints ]またはfreezegrpsを用いて実行します）
総合しますと、今回はそのようなシミュレーションを行っていないとの事ですので、これが原因と思われます。
溶液系と、参照溶媒系が一致しない場合に当たりますので、異常な数字が出る事になります。

もし--flexible条件下でも同様の問題が発生した場合は、.tprファイルならびに.top/.itpファイルを桜庭までお願いします。

よろしくお願いします。

On 2015/09/08 4:28, uran ry wrote:

[uran ry]

石塚様
確認ありがとうございます。


桜庭様
2015年9月8日火曜日 3時42分29秒 UTC-4 Shun:

浦野様

 > では、カットオフとボックスサイズがほとんど同じというのが問題だったのかもしれないですね。
 > 本計算ではそうおきない気がしますが、注意することにします
本件の問題はこの点ではなく、"--rigid"について誤った使い方を行っているのが原因と思われます。

rigid を使う場合、チュートリアルにも明示してありますが、
Solution系のシミュレーションを、分子が剛体であるように扱って実行する必要があります。
（通常、GROMACSではこれは[ constraints ]またはfreezegrpsを用いて実行します）
総合しますと、今回はそのようなシミュレーションを行っていないとの事ですので、これが原因と思われます。
溶液系と、参照溶媒系が一致しない場合に当たりますので、異常な数字が出る事になります。

soln 系のための計算も、剛体で取り扱う必要があるんですね。なるほど、すみませんがそれは見落としていました。
flexible の方では一致していますので問題ないと思います。
refsではただ純溶媒系にテストパーティクルを放り込んで計算しているだけと思っていましたので、多少の上下がでるだけで、
オーダーのズレが出るほどとは思っていませんでした、今後注意します

ありがとうございました

[uran ry]

他の閲覧者用のNoteです

本件につきまして
別トピックに有りますように、一部は
「dodecahedron、octahedronを使用した際の不具合」
であることがわかりました。

私がtrjconv を行う際に溶質をセンタリングしていましたので、dodecahedron をつかった私の計算において
trjconv を行ったあとのsoln 系では
正しい計算が行われるようになったようです。

ermodの最新版では,対応しているようですので、同じ問題に出会った方は
最新版をご利用ください。

浦野

2015年9月8日火曜日 12時23分29秒 UTC-4 uran ry: