计算化学公社

标题: 径向分布函数求解配位数？ [打印本页]
作者
Author: Aristotler    时间: 2018-5-14 22:31
标题: 径向分布函数求解配位数？
各位老师，同学，想向大家求教一个问题：我想通过径向分布函数求解配位数，目前看到两种计算方法：第一种是利用图1的方法，那么如何求解这个积分函数呢.r^2需要积分吗？（数学渣渣）以上公式单位如何转换？

(, 下载次数 Times of downloads: 144)
第二种是如下：
(, 下载次数 Times of downloads: 147)
哪种是对的呢？第一种是怎么积分呢？感谢大家
作者
Author: wbn    时间: 2018-5-15 04:26
本帖最后由 wbn 于 2018-5-15 04:28 编辑

...你说的第二种方法用的也是第一种的公式...

我以前写过一个f90 code通过gmx的rdf输出算coordination num，你可以拿去改改。不保证对啊

  1. ! the program to calculate coordination number from g(r)
  2.       program gr2coln
  3.       implicit none
  4.       type rdf
  5.       real*8 r,gr
  6.       endtype
  7.       character rdffile*50,head
  8.       real*8 peakhead,peaktail,total,boxlength,nip
  9.       integer i
  10.       type(rdf) g_r(2000)
  11.       call getarg(1,rdffile)
  12.       write(*,*) "input where the peak begins in nm"
  13.       read(*,*) peakhead
  14.       write(*,*) "input where the peak ends in nm"
  15.       read(*,*) peaktail
  16.       write(*,*) "the box length in nm"
  17.       read(*,*) boxlength
  18.       write(*,*) "number of ion pairs in box"
  19.       read(*,*) nip
  20.       open(10,file=rdffile,status='old')
  21.       read(10,*) head
  22.       do while(head=='@'.or.head=='#')
  23.           read(10,*) head
  24.       enddo
  25.       i=1
  26.       total=0.0
  27.       read(10,*) g_r(i)%r,g_r(i)%gr
  28.       do while(g_r(i)%r<peakhead)
  29.           i=i+1
  30.           read(10,*) g_r(i)%r,g_r(i)%gr
  31.       enddo
  32.       do while(g_r(i)%r<peaktail)
  33.           i=i+1
  34.           read(10,*) g_r(i)%r,g_r(i)%gr
  35.           total=total+g_r(i)%gr*4*3.14159*g_r(i)%r*g_r(i)%r*(0.002)
  36.       enddo
  37.       total=total*nip/(boxlength*boxlength*boxlength)
  38.       write(*,*) "The coordination number is : ",total
  39.       close(10)
  40.       end
复制代码
作者
Author: sobereva    时间: 2018-5-15 06:48
稍微学点编程，用梯形积分很容易就能做
如果不编程，利用origin的曲线积分功能也完全可以
作者
Author: Aristotler    时间: 2018-5-15 07:52
本帖最后由 Aristotler 于 2018-5-15 07:58 编辑
sobereva 发表于 2018-5-15 06:48
稍微学点编程，用梯形积分很容易就能做
如果不编程，利用origin的曲线积分功能也完全可以

谢谢老师，想向您确认下，第二种方法对吗？我目前只会第二种的
作者
Author: Aristotler    时间: 2018-5-15 07:53
wbn 发表于 2018-5-15 04:26
...你说的第二种方法用的也是第一种的公式...

我以前写过一个f90 code通过gmx的rdf输出算coordination n ...

非常感谢你，想问你第二种方法对吗？
作者
Author: fhh2626    时间: 2018-5-15 11:05
第二种方法也是用第一种方法的公式积分。。。
作者
Author: wbn    时间: 2018-5-16 00:51
Aristotler 发表于 2018-5-15 07:53
非常感谢你，想问你第二种方法对吗？

根本就没有两种方法，只有这一种积分公式，你要在积分里把r^2去掉了那肯定不对啊

PS：我的code里积分直接用的累加，你需要比较精确的结果的话需要梯形积分
作者
Author: Aristotler    时间: 2018-5-16 22:53
sobereva 发表于 2018-5-15 06:48
稍微学点编程，用梯形积分很容易就能做
如果不编程，利用origin的曲线积分功能也完全可以

老师，今天又弄了一天怎么求解配位数。还是有些问题，想求问您：我只要将g(r)用origin积分,取计算第一个峰谷的g(r)积分值就是第一壳层的配位数吗？但是配位数与径向分布函数的关系式里面又包括对r*r的积分，那么以上做法对吗？
作者
Author: sobereva    时间: 2018-5-17 09:23
Aristotler 发表于 2018-5-16 22:53
老师，今天又弄了一天怎么求解配位数。还是有些问题，想求问您：我只要将g(r)用origin积分,取计算第一个 ...

你在origin里增加一列，把g(r)乘上r^2转化成被积函数之后再用origin积分，然后再乘上前面的系数即可
作者
Author: Aristotler    时间: 2018-5-17 18:25
sobereva 发表于 2018-5-17 09:23
你在origin里增加一列，把g(r)乘上r^2转化成被积函数之后再用origin积分，然后再乘上前面的系数即可

非常感谢老师，非常
作者
Author: gpp201013    时间: 2019-7-4 21:00
你好，我目前也在计算径向分布函数，我想知道具体算法，求前辈不吝赐教，拜托啦
作者
Author: Jack    时间: 2019-7-5 07:03
楼主能分享下含有上面公式的文献么
作者
Author: Aristotler    时间: 2019-7-5 08:52
gpp201013 发表于 2019-7-4 21:00
你好，我目前也在计算径向分布函数，我想知道具体算法，求前辈不吝赐教，拜托啦

我做的就是vasp第一性原理的动力学，先用MS的AMORPHOUS CELL 模块建模，之后用vasp进行计算，之后用VMD导出给g（r）之后就是用sob老师上面提到的方法，计算配位数。希望对你有帮助
作者
Author: Aristotler    时间: 2019-7-5 08:54
Jack 发表于 2019-7-5 07:03
楼主能分享下含有上面公式的文献么

文献是在小*虫上截的
作者
Author: highlight    时间: 2019-7-5 09:01
Aristotler 发表于 2019-7-5 08:52
我做的就是vasp第一性原理的动力学，先用MS的AMORPHOUS CELL 模块建模，之后用vasp进行计算，之后用VMD导 ...

有个问题啊，你既然用的是vmd，为啥要自己算积分呢？
http://www.ks.uiuc.edu/Research/vmd/plugins/gofrgui/
The number integrals are computed directly and thus provide accurate coordination numbers.
作者
Author: Aristotler    时间: 2019-7-5 11:22
本帖最后由 Aristotler 于 2019-7-5 15:01 编辑
highlight 发表于 2019-7-5 09:01
有个问题啊，你既然用的是vmd，为啥要自己算积分呢？
http://www.ks.uiuc.edu/Research/vmd/plugins/gof ...

其实我不太会用径向分布函数求配位数的，当时还不知道怎么用VMD求配位数，你知道吗？学习一下
作者
Author: highlight    时间: 2019-7-5 11:52
Aristotler 发表于 2019-7-5 11:22
其实我不太会用径向分布函数求配位数的，但是还不知道怎么用VMD求配位数，你知道吗？学习一下

没多复杂吧，你按照 http://bbs.keinsci.com/thread-7508-1-1.html 10楼的方法做，输出就会多一列数值积分值，然后按照你的方法二找波谷读积分值就是了。
我的意思是，完全没必要自己算积分，程序已经给了。
作者
Author: Aristotler    时间: 2019-7-5 15:00
highlight 发表于 2019-7-5 11:52
没多复杂吧，你按照 http://bbs.keinsci.com/thread-7508-1-1.html 10楼的方法做，输出就会多一列数值积 ...

好的，明白了，非常感谢你
作者
Author: Aristotler    时间: 2019-12-19 20:47
wbn 发表于 2018-5-15 04:26
...你说的第二种方法用的也是第一种的公式...

我以前写过一个f90 code通过gmx的rdf输出算coordination n ...

想问您一下 peakhead, peaktail, total, boxlength, nip对应的是什么意思呢？
作者
Author: Aristotler    时间: 2019-12-19 21:33
wbn 发表于 2018-5-15 04:26
...你说的第二种方法用的也是第一种的公式...

我以前写过一个f90 code通过gmx的rdf输出算coordination n ...

      open(10,file=rdffile,status='old') 请问这个是什么呢？我运行了一下报错
作者
Author: Aristotler    时间: 2019-12-21 09:56
highlight 发表于 2019-7-5 09:01
有个问题啊，你既然用的是vmd，为啥要自己算积分呢？
http://www.ks.uiuc.edu/Research/vmd/plugins/gof ...

计算质心间的配位数VMD实现不了
作者
Author: 星空DMU    时间: 2020-1-7 17:31
sobereva 发表于 2018-5-17 09:23
你在origin里增加一列，把g(r)乘上r^2转化成被积函数之后再用origin积分，然后再乘上前面的系数即可

请问sob老师，怎么转化成被积函数。
作者
Author: sobereva    时间: 2020-1-7 18:39
星空DMU 发表于 2020-1-7 17:31
请问sob老师，怎么转化成被积函数。

你设定列的表达式，设成g(r)那一列的数据与相应的r那一列数据平方的乘积，就成了被积函数
作者
Author: 星空DMU    时间: 2020-1-7 19:42
sobereva 发表于 2020-1-7 18:39
你设定列的表达式，设成g(r)那一列的数据与相应的r那一列数据平方的乘积，就成了被积函数

谢谢sob老师，还有一个问题想向您请教，我计算A原子和B原子之间的配位数，那数密度怎么算？是A和B之间的相对密度还是A和B单独的体系密度。
作者
Author: sobereva    时间: 2020-1-7 19:55
星空DMU 发表于 2020-1-7 19:42
谢谢sob老师，还有一个问题想向您请教，我计算A原子和B原子之间的配位数，那数密度怎么算？是A和B之间的 ...

还是对rdf积分得到
差异在于算rdf的时候谁是参考组是谁被计算的组
作者
Author: 星空DMU    时间: 2020-1-7 21:23
sobereva 发表于 2020-1-7 19:55
还是对rdf积分得到
差异在于算rdf的时候谁是参考组是谁被计算的组

我看书上写的是系统密度。您说的对RDF积分，origin里面结果是面积值；谁是参考组，谁是被计算的组意思是以谁为中心原子，谁是外层分布原子吗？
作者
Author: sobereva    时间: 2020-1-8 21:24
星空DMU 发表于 2020-1-7 21:23
我看书上写的是系统密度。您说的对RDF积分，origin里面结果是面积值；谁是参考组，谁是被计算的组意思是 ...

是
参考组是-ref，被算的组是-sel

那应当叫“体系的密度”而不是“系统的密度”，那本书质量很差，大部分都是抄来的。
rho对应的是被计算的组在盒子中的密度

作者
Author: 星空DMU    时间: 2020-1-9 08:31
sobereva 发表于 2020-1-8 21:24
是
参考组是-ref，被算的组是-sel

感谢sob老师，那被计算的组在盒子中的密度是用被计算的组的个数除以体系的体积吗，密度单位是什么啊。
作者
Author: sobereva    时间: 2020-1-10 01:52
星空DMU 发表于 2020-1-9 08:31
感谢sob老师，那被计算的组在盒子中的密度是用被计算的组的个数除以体系的体积吗，密度单位是什么啊。

是
密度显而易见
作者
Author: 星空DMU    时间: 2020-1-10 08:31
sobereva 发表于 2020-1-10 01:52
是
密度显而易见

谢谢老师。
作者
Author: Edward    时间: 2020-10-26 19:44
楼主你好 请问利用N/V计算密度的话 V是指构建晶体体积还是rdf中以r值为半径的球的体积 谢谢
作者
Author: 顺利毕业发paper    时间: 2021-5-28 09:05
Edward 发表于 2020-10-26 19:44
**** 作者被禁止或删除 内容自动屏蔽 ****

应该是球的体积，因为r^2是在积分里边的，所以体积是个变量



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3