使用Multiwfn做表面活性剂aNCI分析遇到的几个小问题

牧生

看了几遍http://sobereva.com/591博文的操作，自己跃跃欲试，手痒开始了，摩拳擦掌想尝试做表面活性剂aNCI分析，请帮忙看下有无错误。

在Centos下进行的操作，使用amber力场，使用在线服务器得到表面活性剂的GAFF力场的itp。

一、向盒子中加入5个十二烷基苯磺酸钠阴离子SDBS，加满水，能量最小以后，跑npt。再冻结其中一个SDBS （resid 4），去掉控压，跑nvt，总共得到1001帧。

二、使用VMD观察轨迹运动时先载入nvt.gro，删除仅有的一帧，再载入nvt.xtc，看到resid 4没有动，其余的水分子和SDBS都在动。

三、待轨迹运动完了以后，将进度条拖回到第0帧，使用语句resid 4 or same resid as within 5 of resid 4，选择resid 4周围5 A的分子。图像如下：

第一个需要请教的：在这里，有怪异的线条，应该是水分子从盒子另一头穿过来了，这个会影响结果吗？

然后save ，选择resid 4 or same resid as within 5 of resid 4，存为aNCI.xyz

四、运行gmx editconf–f nvt.tpr –o nvt.pdb，然后vmd打开nvt.pdb，同样将resid 4 or same resid as within 5 of resid 4保存为aNCI.pdb。

     这个pdb看起来就很正常。

打开aNCI.pdb看到resid 4的原子编号为52-102

五、使用命令Multiwfn aNCI.xyz，依次选择20    3  1,1001     11    52-102     2.5 A   0.1 （希望得到的图形光滑，所以选择0.1），  计算结束以后选择6，在当前目录输出得到avgRDG.cub和avgsl2r.cub。再Multiwfn avgRDG.cub，选择13, 13, 1.5   100，  2   选择52-102，-2，重新输出avgRDG.cub。

六、将作图脚本avgRDG.vmd拷到当前目录，然后打开vmd以后，source avgRDG.vmd，得到图像如下    （我发现这里无论在当前目录，还是VMD安装目录下操作，没有差别）

需要请教的第二个问题：此前两个图中，水分子都非常靠近SDBS，但是这个图中，水分子就跑开了，这样对吗？

第三个问题，按数字0查看此时等值面内的分子，resid和residue 就都变成1了。。这个编号的变化，是正常的吗？

第四个问题，图中左边蓝色圆圈1中的浅蓝色，可以归结为氧原子与周围水的作用，但我直觉上觉得这个蓝色有点浅，。。上边圆圈2内的红色面可以归结为分子弯曲时，原子间的斥力，但是也同时出现了明显的蓝色，那么，如何解释这个蓝色呢？

我知道其中一定有不正确的操作，但是又不是太清楚错在哪
2021.9.15 编辑，重新使用了整个轨迹去计算，得到结果如下

现在看起来就比较正常了，端基磺酸根上带负电氧原子与周围的原子有强烈的作用，颜色有蓝色和红色。

但是，如果针对大体系，如何才能保存部分位置，以及特定分子及其周围的水分子

sobereva

1 没有影响。实际上你可以把体系用多帧叠加方式显示，只要SDBS附近一圈都包满了水即可

2 播放一下aNCI.xyz轨迹，看看一开始那些水是否后来已经分散开，不再紧密包裹SDBS，如果是这样，aNCI结果没法用。这种情况必须考虑所有的水，和下文的情况一样，这样在轨迹的xyz文件里SDBS附近第一配位层才始终有水
使用Multiwfn研究分子动力学中的弱相互作用
http://sobereva.com/186

3 不用管这个

4 先把前面的2解决了

牧生

sobereva 发表于 2021-9-15 04:31
1 没有影响。实际上你可以把体系用多帧叠加方式显示，只要SDBS附近一圈都包满了水即可

2 播放一下aNCI.x ...

现在情况就是，播放aNCI.xyz轨迹，水分子就会跑散开。

我在另存xyz轨迹之前已经选中了update selection every frame，播放时，SDBS周围的水分子并没有明显跑散开。但是，另存为xyz以后，播放该xyz轨迹时，水分子就散开了。

sobereva

牧生 发表于 2021-9-15 08:57
现在情况就是，播放aNCI.xyz轨迹，水分子就会跑散开。

我在另存xyz轨迹之前已经选中了update selectio ...

所以当前轨迹根本不能用

update selection every frame对于保存轨迹不生效

牧生

sobereva 发表于 2021-9-15 09:38
所以当前轨迹根本不能用

update selection every frame对于保存轨迹不生效

那么，有没有办法，在只冻住SDBS的情况下，保存SDBS及附近水分子为xyz格式的时候，保持SDBS附近一直都有水分子。

有个异想天开的办法，比如一个很大的盒子（原子数有数十万个），足够长时间的npt平衡以后，SDBS聚集形态已经稳定，只有水分子还在到处跑。。此时选择某个SDBS附近15A分子，另存为pdb格式的小盒子，冻结这个SDBS后，再对这个盒子进行nvt（原子数有几千个左右，虽然体系还是有点大，但是服务器也能算的动），结束后，对这个小盒子整体的轨迹进行分析。这样，就能实现考察动态过程中的表面活性剂-水间的相互作用。。

这个，从原则上讲，是否可行？？

sobereva

牧生 发表于 2021-9-15 09:46
那么，有没有办法，在只冻住SDBS的情况下，保存SDBS及附近水分子为xyz格式的时候，保持SDBS附近一直都有 ...

没有直接办法实现。
但可以自己写脚本，循环每一帧，每一帧都用within语句选择SDBS附近一定距离的水，导出成结构文件。但同样的within距离内每一帧的水的数目可能不同，而xyz文件要求必须每一帧原子数相同，因此对每一帧还需要做个循环，逐渐改变within的阈值，当选中的水分子数和第一帧时候选中的水分子数相同时结束循环。

你说的那个异想天开的做法可以。但是你没法恰好保证小盒子里的密度合适，所以取出这么一个子体系后，应当再做一下NPT，等盒子尺寸自发稳定之后再NVT

牧生

sobereva 发表于 2021-9-16 06:39
没有直接办法实现。
但可以自己写脚本，循环每一帧，每一帧都用within语句选择SDBS附近一定距离的水，导 ...

我已经试过了，从大盒子里，取出一个小盒子，再进行NPT去调节密度的过程中，分子都会移动。即使冻结了需要分析的分子，其余的分子都会移动，造成目标分子附近的环境发生变化。如果直接NVT，分子照样会到处移动，液体形状也变化，有可能变成球形。所以，我觉得从大盒子里挖出小盒子，再想办法aNCI分析，此法行不通。再者，我现在觉得，即使强行用大体系去计算其中一个长条形分子的等值面，得到一个毛毛虫形状的等值面图，可能意义也不大。我们更多的是关注亲水头基与附近分子的作用，其余部分，有一些就行了，不必过多去计较。


照着博文http://bbs.keinsci.com/thread-9472-1-7.html 的方法，将我在1楼中，第四步中取出的pdb进行IGM分析

1. 四、运行gmx editconf–f nvt.tpr –o nvt.pdb，然后vmd打开nvt.pdb，同样将resid 4 or same resid as within 5 of resid 4保存为aNCI.pdb。

复制代码

分为了三个组，第一个SDBS，第二个SDBS，第三个为其余的水，进行IGM分析
得到的图形如下，两个SDBS分子仍然靠的近，水分子围绕在SDBS附近，没有散开。


这个看起来就比较正常了，SDBS亲水头基附近出现了明显的等值面，颜色为蓝色的部分，解释为头基的氧原子与水作用，且作用较强。SDBS碳链附近有部分离散的绿色等值面，可解释为碳链与水的作用，或者SDBS之间的弱相互作用，且该种作用较弱。SDBS附近没有毛毛虫形状的等值面，这个应该是比较简洁好看，且好分析解释。


此外，在http://bbs.keinsci.com/thread-21826-1-1.html中，sob也提到了

1. 笔者之前在《通过独立梯度模型(IGM)考察分子间弱相互作用》（http://sobereva.com/407）里还介绍过IGM分析，比起NCI分析的好处是可以自定义片段，只让片段间的相互作用显示出来，就免去了做格点数据屏蔽的步骤，而且哪怕用比aNCI明显更大的格点间距，等值面边缘的锯齿也很不明显。笔者在aNCI的思想上提出了aIGM，即平均的IGM，专用于分析动态环境中的弱相互作用，此分析将另文说明。

复制代码

那么，在大多数情况下，既然IGM的适用范围，出图质量都比NCI好，而且IGM的计算速度也比较快，是不是就应该首推用IGM，在IGM不足的时候，才选用NCI？？此外，极其期待aIGM的博文教程。

sobereva

牧生 发表于 2021-10-3 14:53
我已经试过了，从大盒子里，取出一个小盒子，再进行NPT去调节密度的过程中，分子都会移动。即使冻结了需 ...

IGM一个缺点是等值面太厚，有时候影响考察
aIGM看Multiwfn手册3.23.9节。用法很简单，选主功能20的子功能12照提示操作即可

牧生

sobereva 发表于 2021-10-3 23:01
IGM一个缺点是等值面太厚，有时候影响考察
aIGM看Multiwfn手册3.23.9节。用法很简单，选主功能20的子功 ...

已经试过了使用aIGM方法，就以教程中苯酚-水体系，只用了前10帧，大致看了下效果，那个等值面太厚了，像个馒头，能不能解决一下啊。


如果使用aIGM来分析SDBS与水作用，得到的图如下，效果似乎还没aNCI的那种薄层好看(为了省时间，只用了10帧。如果用1000帧，可能这个等值面会全部包裹住SDBS的头基，导致什么都看不出来。)


这个等值面太厚，即使在isosurface里面去改，仍然是很厚，等值面甚至比原子还大，有时候还会覆盖原子本身。

其次，有个提醒新手的地方，如果使用aIGM，要把IGM_inter.vmd脚本中的cub改成当前实际得到的cub文件名。

sobereva

牧生 发表于 2021-10-7 12:51
已经试过了使用aIGM方法，就以教程中苯酚-水体系，只用了前10帧，大致看了下效果，那个等值面太厚了，像 ...

IGM就是如此，本来就是不如NCI。这是为什么我后来特意提出了IGMH解决这个问题，不过IGMH依赖于波函数，不支持轨迹平均