计算化学公社

 找回密码 Forget password
 注册 Register
Views: 389|回复 Reply: 0
打印 Print 上一主题 Last thread 下一主题 Next thread

[程序/脚本开发] 脚本分享_取动力学轨迹中某分子一定距离内完整分子的cluster

[复制链接 Copy URL]

125

帖子

0

威望

2200

eV
积分
2325

Level 5 (御坂)

本帖最后由 风起~ 于 2025-3-7 20:58 编辑

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

结果发现这个还需要考虑动力学轨迹周期性的问题,还是有一丢丢的麻烦的。这里需要感谢钟老师对代码的技术支持
想着论坛对于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. def remove_duplicates(lst):
  6.     return [lst[i] for i in range(len(lst)) if lst[i] not in lst[:i]]
  7.    
  8. def transcell(atom_group, cell_dim):
  9.     coord = atom_group.positions.copy()
  10.     x_coord = []
  11.     y_coord = []
  12.     z_coord = []
  13.     x_width = max(np.float64(coord[:,0])) - min(np.float64(coord[:,0]))
  14.     y_width = max(np.float64(coord[:,1])) - min(np.float64(coord[:,1]))
  15.     z_width = max(np.float64(coord[:,2])) - min(np.float64(coord[:,2]))
  16.    
  17.     if x_width > cell_dim[0]/2:
  18.         for a in coord[:,0]:
  19.             a=np.float64(a)
  20.             if a < cell_dim[0] /2:
  21.                 x_coord.append(a + cell_dim[0])
  22.             else:
  23.                 x_coord.append(a)
  24.     else:
  25.         x_coord = coord[:,0]
  26.         
  27.     if y_width > cell_dim[1]/2:
  28.         for a in coord[:,1]:
  29.             a=np.float64(a)  
  30.             if a < cell_dim[1] /2:
  31.                 y_coord.append(a + cell_dim[1])
  32.             else:
  33.                 y_coord.append(a)
  34.     else:
  35.         y_coord = coord[:,1]   
  36.         
  37.     if z_width > cell_dim[2]/2:
  38.         for a in coord[:,2]:
  39.             a=np.float64(a)  
  40.             if a < cell_dim[2] /2:
  41.                 z_coord.append(a + cell_dim[2])
  42.             else:
  43.                 z_coord.append(a)
  44.     else:
  45.         z_coord = coord[:,2]           
  46.         
  47.     new_positions = np.stack((x_coord,y_coord,z_coord), axis=-1)
  48.     atom_group.positions = new_positions
  49.     return atom_group   

  50. def pbc_cluster(file_name,coord_file,traj_file,frame,mol_idx,cluster_central,cluster_range):
  51.     u = mda.Universe(coord_file,traj_file)
  52.     for i in range(frame[0],frame[1]):
  53.         u.trajectory[i]
  54.         cell_dims = u.trajectory[i].dimensions[:3]
  55. #       coordinates = u.select_atoms('resid '+str(mol_idx[0])+':'+str(mol_idx[1])).positions
  56.         for j in range(mol_idx[0],mol_idx[1]):
  57.             #选中溶剂质心范围
  58.             selection = u.select_atoms('resid '+str(j)) + u.select_atoms('around '+str(cluster_central)+' resid '+str(j))
  59.             unique_residues  = set(atom.resid for atom in selection)
  60.             unique_selection = u.select_atoms(f"resid {' '.join(map(str, unique_residues))}")
  61.             center = u.select_atoms('resid '+str(j)).center_of_mass()                                                                                                                                                                                                                                   
  62.             # Translate the coordinates                                                                  
  63.             new_position = transcell(unique_selection, cell_dims)                 
  64.             #处理质心问题
  65.             selection_guest = u.select_atoms('resid '+str(j))
  66.             nearby_molecules=[] #溶剂质心与guest的任一原子的距离小于5埃
  67.             for molecule in new_position.split('residue'):         
  68.                 solvent_centroid = molecule.center_of_mass()
  69.                 distances = mda_distances.distance_array(u.select_atoms('resid '+str(j)).positions, solvent_centroid)
  70.                 if any([i <= cluster_central for i in distances]):
  71.                     nearby_molecules.append(molecule)

  72.             #选中溶剂原子相距范围   
  73.             selection_2 = u.select_atoms('resid '+str(j)) + u.select_atoms('around '+str(cluster_range) +' resid '+str(j))
  74.             unique_residues_2  = set(atom.resid for atom in selection_2)
  75.             unique_selection_2 = u.select_atoms(f"resid {' '.join(map(str, unique_residues_2))}")                                                                                                                                                                                                                                 
  76.             # Translate the coordinates                                                                  
  77.             new_position_2 = transcell(unique_selection_2, cell_dims)  
  78.             for molecule in new_position_2.split('residue'):
  79.                 nearby_molecules.append(molecule)
  80.                                                    
  81.             #选中的分子的信息      
  82.             atoms_name = []
  83.             coord = []
  84.             for q in nearby_molecules:
  85.                 for t in q.names:
  86.                     atoms_name.append(t)
  87.                 for p in q.positions:
  88.                     coord.append(p)

  89.             atom_names_2d = [[name] for name in atoms_name]
  90.             merged_array = np.hstack((atom_names_2d, coord))
  91.             merged_array = [list(i) for i in merged_array]
  92.             merged_array = remove_duplicates(merged_array)
  93. #            merged_array = transcell(np.array(merged_array), cell_dims)
  94.             np.savetxt(str(file_name)+'_'+str(i).zfill(5)+'_'+str(j).zfill(3)+'.xyz',np.array(merged_array),\
  95.                        fmt='%-6s %12s %12s %12s',header=str(len(merged_array))+'\n',comments='')

  96. u1 = mda.Universe('A-toluene-npt.gro','A-toluene-npt.xtc')
  97. num_frame = len(u1.trajectory)
  98. pbc_cluster('A-toluene','A-toluene-npt.gro','A-toluene-npt.xtc',[num_frame-500,num_frame],[16,17],6,4)                                          
复制代码



评分 Rate

参与人数
Participants 1
eV +5 收起 理由
Reason
sobereva + 5

查看全部评分 View all ratings

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2025-8-18 05:19 , Processed in 0.629874 second(s), 21 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list