[urano]

お世話になります、名古屋大学の浦野です

現在CGmodel をversion 1.4 をもとにしたgenesis に入れようとしているのですが、

gromacs のinput 情報(.top, .gro)をもとに、祖視化したwater モデルで計算をさせようとすると、

エラーを起こします。

調べると、どうも水分子の原子数が3 (TIP3P,SPC)の場合や、4の場合(TIP4P) は対応しているが、祖視化のモデルで１原子１分子の水モデルは対応してないようです。

そこで、水分子を０にして、祖視化モデルの水分子を水分子でないものとして、扱わせると（系の中に水が０になる）、そのような状況を想定していないのか、配列の0添字にアクセスするというエラーを起こします。

今は、ただのエネルギー関数の確認なので、ひとつだけ、擬似的にTIP3P 水分子を入れて、揚げて、他のはCG model のままにしようとすると、

forrtl: 致命的なエラー (408): fort: (2): 添字 #2 (配列 NUM_NONB_EXCL) に値 23087 (9477 の上限より上) が含まれています。

なるエラーで止まります。

 

ですので、質問としては、

上のような、１つの水分子が１個の原子からなるシステムを実行させるにはどうすればいいでしょうか？

（なお、このinputはGROMACS では正確に動きますので、そこにエラーがある心配はないです）

 

 

 

 

 

[urano]

すいません、genesis のinput を付け忘れました。
coding　の関係でENERGY forcefield としては、CHARMM を利用したいです。

[INPUT]
grotopfile = topol.top
grocrdfile = snap.gro

[OUTPUT]

[ENERGY]
forcefield = CHARMM
electrostatic = PME
switchdist = 8.0
cutoffdist = 10.0
pairlistdist = 12.0
dielec_const = 1.0
pme_alpha = 0.34
pme_ngrid_x = 64
pme_ngrid_y = 64
pme_ngrid_z = 64
pme_nspline = 4
water_model = SOL
vdw_force_switch = NO
output_style = GENESIS

[DYNAMICS]
integrator = LEAP
nsteps = 20
timestep = 0.001
eneout_period = 2
nbupdate_period = 5
iseed = 314159
verbose = yes

[CONSTRAINTS]
rigid_bond = NO
shake_iteration = 500
shake_tolerance = 1.0e-10

[ENSEMBLE]
ensemble = NVE
tpcontrol = NO
temperature = 0

[BOUNDARY]
type = PBC
box_size_x = 61.9
box_size_y = 61.9
box_size_z = 61.9

[SELECTION]
group1 = rnam:MIZ

[RESTRAINTS]

[ckobayashi]

GENESIS 開発チームです。

1点目のご質問はspdynに対しての事かと思われます。現在、spdynでは一部の力場以外では水分子のないシステムは想定していないため、エラーになると思います。現在準備中のupdate版ではその部分は修正される予定です。

2点目ですが、現在のGENESISではGROMACSのファイルフォーマットを使って、CHARMM力場を利用することはできません。マニュアルの4.1をご覧になっていただければと思います。

[urano]

早速の返答、ありがとうございます。

1点目については、修正版では解決するとのことですが、
spdyn だけでなく、atdyn においても同様のエラーを確認しています。
上の水分子のない場合を想定していないはatdyn でも共通ではないのでしょうか？それともinput の条件がありますか次の、力場とinput の組み合わせのせいでしょうか)。

なお、atdyn の場合のエラーメッセージは次です

forrtl: 致命的なエラー (408): fort: (7): ポインター LIST をポインターがターゲットと関連していないときに使用しようとしました。

また、次のように、gromacs と charmm の組み合わせが許されていないとのことですが、その際にでるエラーメッセージはどのようなものでしょうか？何か特定のメッセージはでるでしょうか？それとも単純にseg fault で終わりますか？

2点目について、GROMACS力場を使わない場合に、
GENESIS　において、CHARMM 力場のinput を使う場合に、LJパラメータを 原子種単位でなくて、原子ペア単位で(combniation rule で生み出さない)与える方法は存在するでしょうか？
それとも自分でコードを追加する必要がありますか？

後、CHARMM力場について、確認したいのですが、input にexculuded を与える部分がないのですが、
CHARMM force field を使う場合に、nonbonded のexcluded list (1-2,1-3,1-4) は正確に対応していますでしょうか？ CHARMM 力場の中に、NBXmod (1-5) で使われる部分です。

質問多くてすみませんが、返答よろしくお願いします。

[urano]

度々すみません、
まとめると、２点めについては、CHARMM 力場のNBFIX option とNBXMod option が機能確認済みかということです。

[ckobayashi]

御連絡ありがとうございます。

atdyn/spdyn 双方でエラーを出力するとのこと、エラーログから推測するにatdyn/spdynでは別の理由であることと、力場とファイルフォーマットの組み合わせの問題かと推測しています。

もしお手数でなければ、デバッグオプション”–enable-debug=3″をつけてコンパイルしていただき、そのログをいただければ、もう少しちゃんと回答できるのではと思います。

% ./configure –enable-debug=3

% make clean

% make -j4

% make install

gromacs/charmmの組み合わせ等、ファイルフォーマットと力場についての禁忌となる組み合わせは、現時点の公開版ではエラーが出力されていません。修正版ではgromacs/charmmだけでなく他の組み合わせに対してもエラーストップするようになっています。

CHARMM力場はPSF/PARファイルからの入力になります。こちらは、GROMACSファイルフォーマットとかなり異なるやり方を取っています。

LJのパラメータに関しては、基本的にパラメータ(PAR)ファイルから各原子種のパラメータを読み取
り、combination ruleで作るか、NBFIXで2つの原子種間の例外的なLJパラメータを決めるの2パター
ンしか基本ないのではと思います。

現行のGENESISですと、GromacsファイルフォーマットはGOモデル等に利用するのですが、それらのモデルはGromacsフォーマットの[pairs]を読みますので、setup_enefunc_contact(spdyn, sp_enefunc_gromacs.fpp)/setup_enefunc_nonb(atdyn, at_enefunc_gromacs.fpp)を読む部分を参考にされる
と良いかと思います。

NBFIX optionは対応済みです。

NBXMOD optioinは、NBXMOD=+5対応と考えています。

[urano]

返答ありがとうございます。
NBFIX option に対応されているならば、CHARMM力場に変えても特に問題ないので、そちらに変えることにします。

debug option 入で、ログファイルは、両方共用意したのですが、
こちらのシステム,imageとURL のリンク方法等は見つかるのですが、ファイルの添付が見当たらないのですが、どのようにお送りすればいいのでしょうか？

[ckobayashi]

GENESISのメールアドレス genesis”at”riken.jp (atは@に変更)に送っていただければと思います。

[urano]

メールでお送りしました。何かわかりましたらよろしくお願いします

[ckobayashi]

ログと入力ファイルの送信ありがとうございました。

1. ログに出力されているエラーを起こしているライン数に対応するソースコードがない
2. GENESISの公開版をそのまま動かすと以下のエラーでストップする
   Setup_Solute_And_Water> # of water is incorrect. rank_no = 0

上の理由から、他にユーザが修正された跡が見えますので、ユーザ修正されたコードと解釈しました。ユーザの方が修正されたコードについては正直、対応できかねますが、行数から類推するにsp_domain.fpp内のsubroutine setup_atom_by_table内の水分子内のdomainの設定(original codeの2850/2858/2866行目のいずれか)から呼ばれた、subroutine molecule_to_domain内でのエラーではないかと思います。(original codeで3952行目)

私の推測が正しい場合は、水分子のリストが正しく作られていないのが原因だと思いますので、enefunc%table%water_listを作成したあたりを確認されたらいかがでしょうか？

[urano]

返信ありがとうございます。

1については、エラーの場所（サブルーチン）を確認するために、write文を入れているのと水分子1原子で動かせないと試行錯誤したのがコメントアウトされただけで残っているので、公開版とずれています。現在input 周りはwrite(MSG) 以外は変更していないので、その辺を修正前の公開版に戻したので、log を再度送り直します。

2. については、genesis　の同じところで止まります。推測については、まさにそこが最初の質問に関係するのですが、水分子が3原子や4原子の場合は動くのですが、祖視化したモデルのような水分子1こが１つの原子になる場合は動かないことが問題になっています。

一時的に0 stepだけでも動かそうとして、ダミーの水分子を１つだけ作成し、CGの水分子を水分子以外として取扱うようにしたがエラーを吐いたのが、新しく送信しなおしたlog_charmmです。

エネルギーのdebug をしたいので、0 stepだけでも動かしてエネルギーを出力することはできませんか？それともたくさんの水分子がないと、genesis は動かないのでしょうか(おそらく各ドメインに1つ以上の水分子)。

logとinput は再度メールでお送りします

よろしくお願いします

 

 

 

[ckobayashi]

ログを送付いただいて再度確認しました。

そこで分かったことは、水分子がないことが問題ではなく、CG分子の分子内相互作用(dihedral angle)の1-4の原子間距離がpairlist距離の半分より長い事が原因でした。(10Åを越える物が複数見られました）

pairlist距離を20Å程度にして、MPI数を1にするといただいたインプットでも動作すること自体は確認しました。

[urano]

parilist を20にした場合動きました。

とりあえず、１つだけでも、DUMMYの水分子を入れれば動くということですね。とりあえず、これで試してみることにします。

ありがとうございます。