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

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

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

感谢楼主，这个O原子类型他的原子序号很乱，也不是那种连续的，如果用VMD是挨个确定好原子序号编进去吗

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
   
5. 
   
6. 
   
7. 
   
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.
    
11. 
    
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.     """
    
20. 
    
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]
    
30. 
    
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.
   
4. 
   
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

复制代码