计算化学公社

标题: 脚本分享_取动力学轨迹中某分子一定距离内完整分子的cluster [打印本页]

---

作者
Author: 风起~    时间: 2025-3-7 20:55
标题: 脚本分享_取动力学轨迹中某分子一定距离内完整分子的cluster


写在前面：当时有一个工作需要考虑动力学静态无序，就写了这样一个脚本。

结果发现这个还需要考虑动力学轨迹周期性的问题，还是有一丢丢的麻烦的。这里需要感谢钟老师对代码的技术支持
想着论坛对于MDAnalysis的代码还是比较少的，现在分享出来，大家可以参考一下。

pbc_cluster('A-toluene','A-toluene-npt.gro','A-toluene-npt.xtc',[num_frame-500,num_frame],[16,17],6,4)
指的是，对于A-toluene-npt.xtc轨迹中，选取最后500帧的第17个分子的{质心6埃内，任一原子4埃内}的完整分子的cluster

1. import os
   
2. import numpy as np
   
3. import MDAnalysis as mda
   
4. import MDAnalysis.analysis.distances as mda_distances
   
5. 
   
6. def remove_duplicates(lst):
   
7.     return [lst[i] for i in range(len(lst)) if lst[i] not in lst[:i]]
   
8.    
   
9. def transcell(atom_group, cell_dim):
   
10.     coord = atom_group.positions.copy()
    
11.     x_coord = []
    
12.     y_coord = []
    
13.     z_coord = []
    
14.     x_width = max(np.float64(coord[:,0])) - min(np.float64(coord[:,0]))
    
15.     y_width = max(np.float64(coord[:,1])) - min(np.float64(coord[:,1]))
    
16.     z_width = max(np.float64(coord[:,2])) - min(np.float64(coord[:,2]))
    
17.    
    
18.     if x_width > cell_dim[0]/2:
    
19.         for a in coord[:,0]:
    
20.             a=np.float64(a)
    
21.             if a < cell_dim[0] /2:
    
22.                 x_coord.append(a + cell_dim[0])
    
23.             else:
    
24.                 x_coord.append(a)
    
25.     else:
    
26.         x_coord = coord[:,0]
    
27.         
    
28.     if y_width > cell_dim[1]/2:
    
29.         for a in coord[:,1]:
    
30.             a=np.float64(a)  
    
31.             if a < cell_dim[1] /2:
    
32.                 y_coord.append(a + cell_dim[1])
    
33.             else:
    
34.                 y_coord.append(a)
    
35.     else:
    
36.         y_coord = coord[:,1]   
    
37.         
    
38.     if z_width > cell_dim[2]/2:
    
39.         for a in coord[:,2]:
    
40.             a=np.float64(a)  
    
41.             if a < cell_dim[2] /2:
    
42.                 z_coord.append(a + cell_dim[2])
    
43.             else:
    
44.                 z_coord.append(a)
    
45.     else:
    
46.         z_coord = coord[:,2]           
    
47.         
    
48.     new_positions = np.stack((x_coord,y_coord,z_coord), axis=-1)
    
49.     atom_group.positions = new_positions
    
50.     return atom_group   
    
51. 
    
52. def pbc_cluster(file_name,coord_file,traj_file,frame,mol_idx,cluster_central,cluster_range):
    
53.     u = mda.Universe(coord_file,traj_file)
    
54.     for i in range(frame[0],frame[1]):
    
55.         u.trajectory[i]
    
56.         cell_dims = u.trajectory[i].dimensions[:3]
    
57. #       coordinates = u.select_atoms('resid '+str(mol_idx[0])+':'+str(mol_idx[1])).positions
    
58.         for j in range(mol_idx[0],mol_idx[1]):
    
59.             #选中溶剂质心范围
    
60.             selection = u.select_atoms('resid '+str(j)) + u.select_atoms('around '+str(cluster_central)+' resid '+str(j))
    
61.             unique_residues  = set(atom.resid for atom in selection)
    
62.             unique_selection = u.select_atoms(f"resid {' '.join(map(str, unique_residues))}")
    
63.             center = u.select_atoms('resid '+str(j)).center_of_mass()                                                                                                                                                                                                                                   
    
64.             # Translate the coordinates                                                                  
    
65.             new_position = transcell(unique_selection, cell_dims)                 
    
66.             #处理质心问题
    
67.             selection_guest = u.select_atoms('resid '+str(j))
    
68.             nearby_molecules=[] #溶剂质心与guest的任一原子的距离小于5埃
    
69.             for molecule in new_position.split('residue'):         
    
70.                 solvent_centroid = molecule.center_of_mass()
    
71.                 distances = mda_distances.distance_array(u.select_atoms('resid '+str(j)).positions, solvent_centroid)
    
72.                 if any([i <= cluster_central for i in distances]):
    
73.                     nearby_molecules.append(molecule)
    
74. 
    
75.             #选中溶剂原子相距范围   
    
76.             selection_2 = u.select_atoms('resid '+str(j)) + u.select_atoms('around '+str(cluster_range) +' resid '+str(j))
    
77.             unique_residues_2  = set(atom.resid for atom in selection_2)
    
78.             unique_selection_2 = u.select_atoms(f"resid {' '.join(map(str, unique_residues_2))}")                                                                                                                                                                                                                                 
    
79.             # Translate the coordinates                                                                  
    
80.             new_position_2 = transcell(unique_selection_2, cell_dims)  
    
81.             for molecule in new_position_2.split('residue'):
    
82.                 nearby_molecules.append(molecule)
    
83.                                                    
    
84.             #选中的分子的信息      
    
85.             atoms_name = []
    
86.             coord = []
    
87.             for q in nearby_molecules:
    
88.                 for t in q.names:
    
89.                     atoms_name.append(t)
    
90.                 for p in q.positions:
    
91.                     coord.append(p)
    
92. 
    
93.             atom_names_2d = [[name] for name in atoms_name]
    
94.             merged_array = np.hstack((atom_names_2d, coord))
    
95.             merged_array = [list(i) for i in merged_array]
    
96.             merged_array = remove_duplicates(merged_array)
    
97. #            merged_array = transcell(np.array(merged_array), cell_dims)
    
98.             np.savetxt(str(file_name)+'_'+str(i).zfill(5)+'_'+str(j).zfill(3)+'.xyz',np.array(merged_array),\
    
99.                        fmt='%-6s %12s %12s %12s',header=str(len(merged_array))+'\n',comments='')
    
100. 
     
101. u1 = mda.Universe('A-toluene-npt.gro','A-toluene-npt.xtc')
     
102. num_frame = len(u1.trajectory)
     
103. pbc_cluster('A-toluene','A-toluene-npt.gro','A-toluene-npt.xtc',[num_frame-500,num_frame],[16,17],6,4)                                          
     

复制代码

---

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

Powered by Discuz! X3.3