在gromacs里应该如何计算体系的结构因子？

diaok

在gromacs里应该如何计算体系的结构因子？

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

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


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

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


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

wbn

我算过很多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

diaok

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）应该没太大贡献吧？

wbn

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方法，氢的贡献太大了。

diaok

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值的归一化吧？

wbn

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的话理论和实验峰高差个一两倍也是有可能的...

diaok

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

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


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

wbn

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

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

yuyu116

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求积分，就不太懂了

少年爱吃地三鲜

请问楼主，可否知道在分子动力学模拟中，静态结构因子S(q) 和散射强度I(q)直接按照什么换算公式计算？

weaey

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

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

wbn

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

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

wbn

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

有什么问题直接问吧

weaey

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进行套嵌吗。

wbn

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个