计算化学公社

标题: 计算rdf和配位数和snapshot看起来不一样 [打印本页]

作者
Author: www1 时间: 2023-3-10 09:55
标题: 计算rdf和配位数和snapshot看起来不一样
本帖最后由 www1 于 2023-3-27 12:48 编辑

请教各位老师同学，

我跑完gromacs的md模拟K离子在水溶液中的行为（主要是要看K的水合离子的周围是有几个水分子）之后，使用vmd进行rdf的统计和配位数的计算。rdf的峰明显低很多，配位数算出来只有2。同样的力场文献里的配位数是4到5。
我check了轨迹之后，发现看起来轨迹没有什么问题，确实是跟文献里报道的基本一致，但是无论我怎么样统计出来的结果都还是有问题的。所以请大家帮我看看，到底是哪里的问题，感谢。

下面附上一些截图。

1. 我用graphical representation 里用selected atom 检查了我的选择原子的语法，是没问题的。蓝色是K 离子，红色是H2O里的O。也能大致看得出配位数肯定不止2.

2. 下面是我进行rdf calculation的时候的选择。

selection 1： name K, selection2 name OW.
3. 计算出来的rdf跟文献里的趋势是完全一致的，只是高度不一样。计算配位数是用自己写的脚本计算的，也test过没有任何问题了。

------------------------------------------3-27 update

脚本里的公式进行了更正，但是计算出来的rdf和轨迹中观察到的，还是不一致~请问是什么问题~

------------------------------------------

希望各位老师同学可以帮忙找到问题所在，非常感谢！！

下面附上我根据rdf计算配位数的脚本，有需要的同学可以拿去用。

import numpy as np

rdf_file = "rdf.dat"
rdf_data = np.loadtxt(rad_file)

r_min =2
r_max = 3.5

r_idx = np.where((rdf_data[:,0] > r_min) & (rdf_data[:,0] < r_max))[0]
#coord_number = np.trapz(rdf_data[r_idx,1], rdf_data[r_idx,0])
coord_num = 4 * np.pi * density * np.trapz(r[idx_min:idx_max]**2 * g_r[idx_min:idx_max], r[idx_min:idx_max])

print("Coordination number:")

作者
Author: lyj714 时间: 2023-3-10 11:34
配位数公式你用的不对啊，是4pi rho r**2 gr dr

作者
Author: www1 时间: 2023-3-11 11:32

lyj714 发表于 2023-3-10 11:34
配位数公式你用的不对啊，是4pi rho r**2 gr dr

啊…… 是这样。谢谢

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