计算化学公社

标题: 在gromacs里应该如何计算体系的结构因子？ [打印本页]

作者
Author: diaok 时间: 2017-3-22 21:54
标题: 在gromacs里应该如何计算体系的结构因子？
在gromacs里应该如何计算体系的结构因子？

gmx sans和saxs对大的体系完全算不动

saxs得到的结果和盒子有关，但是原理是hf wave拟合
sans的原理看起来像散射公式，但是精度却和盒子大小无关
这两个命令分别适用于什么体系？

然后gmx4.6.7版的rdf可以-sq得结构因子，但是后面的版本都没有了
不知道是什么原因。

实验上一般sans和saxs结果峰位置类似

这些计算方法有什么适用范围？如何取舍？
另外x射线散射原理中散射强度还和电子数目有关，上面的计算里如何引入这个贡献？

作者
Author: wbn 时间: 2017-3-23 03:24
我算过很多X-ray structure factor，还用这个灌了不少水 http://dx.doi.org/10.1063/1.4972410 doi: http://dx.doi.org/10.1063/1.4962257
这种方法通过g_rdf 得到的rdf 来进行变换，就可以得到S(q)了。具体公式见两篇文章的公式（1）。用这种方法，如果模拟力场合格的话算出的S(q)和实验基本是重合的，见文章中的图。

具体的Fourier变换过程我写了段fortran代码来实现，如果老板同意我可以把代码发给你，不过我写的代码很蹩脚，如果你试着用matlab来做可以更加简单。

等你做出来了拜托引一下我的文章吧... 虽然是灌水，但一个引用也没有太尴尬了...

我没有试过gromacs里自带的方法，但是看起来似乎不是很靠谱, manual 里也没有说他的算法。Cromer's method 是一种计算form factor的方法，具体怎么算S(q)他也没说。算S(q)一定要给所有原子assign element，我不知道在gmx里怎么实现这一点。

原子散射X-ray的能力大致与电子数成正比，但是 Cromer 的贡献就是用HF方法算出来各个原子的atomic form factor，然后用公式拟合F(q)与q的关系。具体数字在他196X年的一篇JCP和international tables of crystalline volume C 都可以查到。我的文章中就是用的他的数值。

我还没有做过neutron scattering, 但是用MD轨迹生成S(q)的方式在neutron scattering分析里似乎不大常见，更常见的方法是先做neutron scattering，用得到的衍射图样拟合Monte Carlo 模拟的结构，具体搜索 Empirical Potential Structure Refinement

作者
Author: diaok 时间: 2017-3-23 14:02
本帖最后由 diaok 于 2017-3-23 14:22 编辑

wbn 发表于 2017-3-23 03:24
我算过很多X-ray structure factor，还用这个灌了不少水 http://dx.doi.org/10.1063/1.4972410 doi: http: ...

首先感谢指导

gmx saxs的时候有assign element
就是太慢不知道和rdf的比起来怎么样
现在确定了g_rdf的-sq和saxs完全一样也需要指定元素
但是这个fft过程为什么用到了元素，和我们用的好像不同

rdf到s（q）的matlab傅里叶变换我这里已经有了，之前粗粒化体系之间是用这个算了体系里的亲水组分，完全没考虑其他组分和交换作用等

1）你的文章中提到的公式里面是否用到了cromer的系数？
rdf用基团团质心时直接算基团的总电子数目吗？
2）我的体系还是粗粒化的，完全还原回全原子存在困难，目前只实现了部分还原
也不知道是否需要还原
氢原子对体系的s（q）应该没太大贡献吧？

作者
Author: wbn 时间: 2017-3-23 23:06
本帖最后由 wbn 于 2017-3-24 00:39 编辑

diaok 发表于 2017-3-23 14:02
首先感谢指导

gmx saxs的时候有assign element

（1）我用了Cromer的系数来算F(q)，电子数目是对atomic form factor F(q) 在q=0时的估计，最精确的方法是把所有原子对的rdf都算出来，然后把每个原子对的S(q)分别乘以两种原子的F(q)，再累加起来。

（2）氢原子对X-ray scattering 贡献有限，但是最好不要忽略掉。但是neutron scattering 他的贡献可是比谁都大。

如果你的力场没法做到all atom的话，我建议的方法是用我文章中的公式（1），把你所有的粒子都当作一种原子

(1) 算出所有粒子两两之间的rdf，做傅立叶变换

(2) 将变换结果乘以两种粒子的F(q)，得到两种粒子间的S(q)。一个粒子的F(q)是里面所有原子F(q)的累加。原子F(q)用Cromer系数算

(3)把所有partial S(q)全加起来得到总S(q)，当然要分别乘以两种粒子的fraction

理论上分子模拟得到的S(q)中q最小只能到2pi/L, L是box length，而你又没有用AA力场，那么 q > 2 A^{-1} 的部分肯定也不对了。所以除非你的盒子很大，否则理论上你算出来的S(q)能和实验能吻合的部分就非常有限。

不要用gmx saxs算，他应该是没考虑你这种情况。

neutron scattering建议用全原子Monte Carlo方法，氢的贡献太大了。

作者
Author: diaok 时间: 2017-3-24 20:08
本帖最后由 diaok 于 2017-3-24 21:28 编辑

wbn 发表于 2017-3-23 23:06
（1）我用了Cromer的系数来算F(q)，电子数目是对atomic form factor F(q) 在q=0时的估计，最精确的方法 ...

公式大概理解了
以前用的是另一个公式的结构因子
结果里没有用过F（q），只对密度归一化了
不归一化的时候只是多了个常数1

公式1里面的洛伦兹W(r)是必须的吗？
在matlab中加入这个完全变成震荡了
我直接用的累加得到积分已解决乘除号错了。。。这个一般要用吗？

另外还有rdf的归一化问题：
先得到所有归一化的rdf，求出所有fft结果
然后求得各组成的f（q），最终统一归一化
这个流程正确吗？归一化是针对曲线里横坐标上每个固定q值的归一化吧？

作者
Author: wbn 时间: 2017-3-25 02:24
本帖最后由 wbn 于 2017-3-25 02:48 编辑

diaok 发表于 2017-3-24 20:08
公式大概理解了
以前用的是另一个公式的结构因子
结果里没有用过F（q），只对密度归一化了

我做的的是小分子均相体系WAXS，你如果做的是SAXS的话可以不用W(r)。但是要保证你的盒子够大，平衡够充分，不然low q 部分理论和实验对应很难的。

但是加入这个变成震荡不太科学，他的效果是在low q的时候使曲线变平滑才对，你是不是 theta 定义错了？

两种粒子间的partial rdf 就直接用g_rdf的结果，就是收敛到1的那个，不要做额外处理。整个过程不需要归一化，严格按照公式(1)把所有参数都设置对的话最后的结果是直接可以和实验比较的，最多可能需要在结果加上1，因为有的程序提取的S(q)是收敛到1的，而公式(1)收敛到0.

当然你如果只做了SAXS的话有可能计算和实验峰高偏差比较大...因为从实验intensity还原出S(q)还是比较困难的，很多参数靠猜。因为我们用的软件（pdfgetX2）从实验提取S(q)是最后收敛到1的，我们一般测到20 A-1，通过调节参数使high q的部分收敛到1，这样low q的部分就会和计算吻合。如果只做SAXS的话理论和实验峰高差个一两倍也是有可能的...

作者
Author: diaok 时间: 2017-4-15 16:42

wbn 发表于 2017-3-25 02:24
我做的的是小分子均相体系WAXS，你如果做的是SAXS的话可以不用W(r)。但是要保证你的盒子够大，平衡够充分 ...

尝试了下用Cromer的系数换算，结果得到的峰非常不明显。。。。
按照这个系数换算后峰位置有可能移动吗？

也可能是我模型的问题，
体系里面6种物质
得到了所有物质的各种交叉的RDF
6
5
4
3
2
1
分子 = sum（1-6本身的项 * 系数^2）+2*（各种交叉项*系数*系数）
下面这部分的原理应该没错吧？

作者
Author: wbn 时间: 2017-4-15 22:48

diaok 发表于 2017-4-15 16:42
尝试了下用Cromer的系数换算，结果得到的峰非常不明显。。。。
按照这个系数换算后峰位置有可能移动吗？ ...

看不懂你问的是什么哎，如果你不清楚做的对不对的话请详细描述一下你的算法，说的时候请用术语，partial radial distribution gA-B(r), atomic form factor F(q), partial structure factor SA-B(q) 之类的语言，你说的什么项，系数之类的东西我不明白在说什么

作者
Author: yuyu116 时间: 2020-8-20 09:30

wbn 发表于 2017-3-23 23:06
（1）我用了Cromer的系数来算F(q)，电子数目是对atomic form factor F(q) 在q=0时的估计，最精确的方法 ...

老师您好，文献里有这样的傅里叶公式 S(q)=1+4*Pi*P 积分号 sin（qr）/qr * g（r）*r*r*dr 。散射因子q最小值取2*π/L，然后依次取n*2*π/L，n=1,2,3，...。不知我的理解是否正确，q和r有什么关系吗,如果没关系，有两个变量q和r，还是对r求积分，就不太懂了

作者
Author: 少年爱吃地三鲜 时间: 2021-1-2 17:07
请问楼主，可否知道在分子动力学模拟中，静态结构因子S(q) 和散射强度I(q)直接按照什么换算公式计算？

作者
Author: weaey 时间: 2021-5-11 21:02

wbn 发表于 2017-3-23 03:24
我算过很多X-ray structure factor，还用这个灌了不少水 http://dx.doi.org/10.1063/1.4972410 doi: http: ...

您好，您的文章我看了，关于结构因子的公式有点问题想要请假，可以加qq聊一下吗，我的qq是645344968

作者
Author: wbn 时间: 2021-5-12 10:18

yuyu116 发表于 2020-8-20 09:30
老师您好，文献里有这样的傅里叶公式 S(q)=1+4*Pi*P 积分号 sin（qr）/qr * g（r）*r*r*dr 。散射因子q ...

这就是空间变换啊，了解一下傅立叶变换的定义

作者
Author: wbn 时间: 2021-5-12 10:19

weaey 发表于 2021-5-11 21:02
您好，您的文章我看了，关于结构因子的公式有点问题想要请假，可以加qq聊一下吗，我的qq是645344968

有什么问题直接问吧

作者
Author: weaey 时间: 2021-5-12 17:31

wbn 发表于 2021-5-12 10:19
有什么问题直接问吧

1、这个原子散射因子是一个常数还是随q变化，我看了您引用的文章，散射因子是随角度变化的，如果在公式中随q变化的话，应该如何转化。
2、您在文章中偏结构因子的分母是xifi(q)求和的平方，可以使用xjfj(q)吗，使用xifi(q)的意义是以i原子为中心寻找j原子吗。
3、S(q)分子上积分前的两个求和的意义是什么，如果我的i原子有56个，j原子有24个，写fortran时，是用两个do循环i=1,56和j=1,24进行套嵌吗。

作者
Author: wbn 时间: 2021-5-12 17:55

weaey 发表于 2021-5-12 17:31
1、这个原子散射因子是一个常数还是随q变化，我看了您引用的文章，散射因子是随角度变化的，如果在公式中 ...

1. 随q变化。q(r)如何获得见International Tables for Crystallography: Volume C, 如果你要把它当作一个常数的话就用原子序数，虽然q(r)随r变化很大，但是不同元素的比值随r变化不大
2.分母就是把它归一化一下，不弄的话分子分母量纲都不统一
3.求和是对所有的原子类型求和，i和j都应该有n个

作者
Author: weaey 时间: 2021-5-12 19:14
本帖最后由 weaey 于 2021-5-12 20:30 编辑

wbn 发表于 2021-5-12 17:55
1. 随q变化。q(r)如何获得见International Tables for Crystallography: Volume C, 如果你要把它当作一个 ...

十分感谢，差不多看明白了

作者
Author: Nieqixia 时间: 2021-12-1 14:14

diaok 发表于 2017-3-23 14:02
首先感谢指导

gmx saxs的时候有assign element

您好，请问fdf到s（q）的matlab傅立叶变换可以分享吗

作者
Author: will维维 时间: 2021-12-7 10:45
本帖最后由 will维维于 2021-12-7 11:34 编辑

wbn 发表于 2017-3-25 02:24
我做的的是小分子均相体系WAXS，你如果做的是SAXS的话可以不用W(r)。但是要保证你的盒子够大，平衡够充分 ...

老师您好，计算gA-B(r)的时候，如果原子类型A、B在同一个分子内（尤其是成键的两原子），那么就会出现一个极高的分子内峰，对S(q)的计算结果影响很大。请问我计算S(q)的时候需要排除gA-B(r)的分子内峰吗？
另外，实验上做saxs需要扣掉纯水的基底，不知道用文献里的eq(1)计算的时候, 是不是也要另外做一个纯水的模拟？

作者
Author: StormSpirts 时间: 2022-3-22 14:41
本帖最后由 StormSpirts 于 2022-4-10 16:06 编辑

搜集了一些关于GROMACS的Cromer and Mann method的资料，发上来给大家作参考：

GROMACS的Cromer and Mann method的源代码：https://fossies.org/dox/gromacs- ... _source.html#l00075

Cromer and Mann method原文：http://scripts.iucr.org/cgi-bin/paper?S0567739468000550

一篇使用Cromer and Mann method处理X射线数据，对比MD结构参数的文章，原理写得挺详细：https://journals.iucr.org/s/issues/2021/05/00/xh5057/index.html
文中列出了atomic form factor和structure factor的公式，两者形式非常相似。又因为atomic form factor可以用Cromer and Mann method方法近似得到，所以GROMACS用Cromer and Mann method方法近似计算结构因子可能是有合理性的（但我对这方面不懂，不敢乱说）。

作者
Author: hanshuai222 时间: 2024-10-24 22:03

wbn 发表于 2017-3-23 03:24
我算过很多X-ray structure factor，还用这个灌了不少水 http://dx.doi.org/10.1063/1.4972410 doi: http: ...

请问代码能分享下吗

欢迎光临计算化学公社 (http://bbs.keinsci.com/)