计算化学公社

标题: 求助cp2k计算得到的轨迹文件统计Ca-Oc和Ca-Ow键个数 [打印本页]
作者
Author: 李思玫    时间: 2025-3-2 18:16
标题: 求助cp2k计算得到的轨迹文件统计Ca-Oc和Ca-Ow键个数
采用cp2k软件计算分子动力学，得到轨迹文件。已知一个体系中含有Ca,Si,H,O和C，其中O分为H2O中的Ow和碳酸分子中的Oc,以及Si-(OH)4四面体项链的Osi，现在通过cp2k软件得到一个.xyz的轨迹文件，轨迹文件中没有区分Ow，Oc和Osi，怎么利用python脚本统计Ca-Ow，Ca-Oc和Ca-Osi化学键数量分被是多少呢
作者
Author: 丁越    时间: 2025-3-3 08:57
本帖最后由 丁越 于 2025-3-3 09:22 编辑

Oc, Osi, Ow在xyz文件中原子序号是固定不变的，直接通过原子序号查找某类O原子不就行了，再结合VMD的原子选择语法比如index 0 1 2 3选择Ow。之后就是用TCL语法遍历计算选中的原子与体系所有原子的键长，你自己设定一个判据标准，小于多少则视为成键。但缺点是TCL脚本运行起来挺慢。

第二种方法是通过ASE，里面的几何模块中有计算原子间的距离矩阵的函数，下面给你一段代码你自己参考一下：
  1. import numpy as np
  2. from ase import io
  3. from ase.data import covalent_radii
  4. from ase.geometry import get_distances



  8. def coordinated_bonds_num(atoms, index=None, scale=1.15):
  9.     """
  10.     calculate the number of coordinated bonds of each atom, simply based on geometry criteria.

  12.     parameters:
  13.         atoms: an ASE Atoms object.
  14.         index: list.
  15.             atomic index to be calculated.
  16.         scale: float.
  17.             the scale factor of covalent radii.
  18.     return: list.
  19.     """

  21.     _, distances = get_distances(atoms.positions, pbc=atoms.pbc, cell=atoms.cell)
  22.     radiis = (covalent_radii[atoms.get_atomic_numbers()] +
  23.               covalent_radii[atoms.get_atomic_numbers()].reshape(-1, 1))
  24.     distances = np.round(distances, decimals=1)
  25.     radiis = np.round(radiis, decimals=1)
  26.     is_bonded = np.where(distances <= scale * radiis, 1, 0) - np.eye(len(atoms))
  27.     bond_num = np.sum(is_bonded, axis=1)
  28.     if index is not None:
  29.         bond_num = bond_num[index]

  31.     return bond_num
复制代码
最后注意一点是CP2K的MD轨迹是不考虑PBC的，所以不论是拿VMD还是ASE统计都需要把原子卷回到盒子内。比如用ASE，参考下面代码：
  1. def wraped_coord(atoms, cell=None):
  2.     """
  3.     wrap atomic coordinates into unit cell.

  5.     parameters:
  6.         atoms: An Atoms object of ASE.
  7.         cell: 3 * 3 nparray. cell parameters in vector formate.
  8.     return:
  9.         an Atoms object with wrapped coordinates.
  10.     """
  11.     if not all(atoms.pbc):
  12.         if cell is None:
  13.             raise ValueError("Cell parameters must be defined when atoms are wrapped.")
  14.         else:
  15.             atoms.set_cell(cell)
  16.             atoms.pbc = (True, True, True)
  17.             atoms.center(vacuum=7.5, axis=2)
  18.             atoms.wrap(pretty_translation=True)
  19.     return atoms
复制代码



作者
Author: 李思玫    时间: 2025-3-3 10:29
丁越 发表于 2025-3-3 08:57
Oc, Osi, Ow在xyz文件中原子序号是固定不变的，直接通过原子序号查找某类O原子不就行了，再结合VMD的原子选 ...

感谢楼主，这个O原子类型他的原子序号很乱，也不是那种连续的，如果用VMD是挨个确定好原子序号编进去吗
作者
Author: 丁越    时间: 2025-3-3 10:45
李思玫 发表于 2025-3-3 10:29
感谢楼主，这个O原子类型他的原子序号很乱，也不是那种连续的，如果用VMD是挨个确定好原子序号编进去吗

轨迹文件中原子序号又不随时间变化。gview也行，自己分类一下O原子类型。



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