|
本帖最后由 风起~ 于 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
- import os
- import numpy as np
- import MDAnalysis as mda
- import MDAnalysis.analysis.distances as mda_distances
- def remove_duplicates(lst):
- return [lst[i] for i in range(len(lst)) if lst[i] not in lst[:i]]
-
- def transcell(atom_group, cell_dim):
- coord = atom_group.positions.copy()
- x_coord = []
- y_coord = []
- z_coord = []
- x_width = max(np.float64(coord[:,0])) - min(np.float64(coord[:,0]))
- y_width = max(np.float64(coord[:,1])) - min(np.float64(coord[:,1]))
- z_width = max(np.float64(coord[:,2])) - min(np.float64(coord[:,2]))
-
- if x_width > cell_dim[0]/2:
- for a in coord[:,0]:
- a=np.float64(a)
- if a < cell_dim[0] /2:
- x_coord.append(a + cell_dim[0])
- else:
- x_coord.append(a)
- else:
- x_coord = coord[:,0]
-
- if y_width > cell_dim[1]/2:
- for a in coord[:,1]:
- a=np.float64(a)
- if a < cell_dim[1] /2:
- y_coord.append(a + cell_dim[1])
- else:
- y_coord.append(a)
- else:
- y_coord = coord[:,1]
-
- if z_width > cell_dim[2]/2:
- for a in coord[:,2]:
- a=np.float64(a)
- if a < cell_dim[2] /2:
- z_coord.append(a + cell_dim[2])
- else:
- z_coord.append(a)
- else:
- z_coord = coord[:,2]
-
- new_positions = np.stack((x_coord,y_coord,z_coord), axis=-1)
- atom_group.positions = new_positions
- return atom_group
- def pbc_cluster(file_name,coord_file,traj_file,frame,mol_idx,cluster_central,cluster_range):
- u = mda.Universe(coord_file,traj_file)
- for i in range(frame[0],frame[1]):
- u.trajectory[i]
- cell_dims = u.trajectory[i].dimensions[:3]
- # coordinates = u.select_atoms('resid '+str(mol_idx[0])+':'+str(mol_idx[1])).positions
- for j in range(mol_idx[0],mol_idx[1]):
- #选中溶剂质心范围
- selection = u.select_atoms('resid '+str(j)) + u.select_atoms('around '+str(cluster_central)+' resid '+str(j))
- unique_residues = set(atom.resid for atom in selection)
- unique_selection = u.select_atoms(f"resid {' '.join(map(str, unique_residues))}")
- center = u.select_atoms('resid '+str(j)).center_of_mass()
- # Translate the coordinates
- new_position = transcell(unique_selection, cell_dims)
- #处理质心问题
- selection_guest = u.select_atoms('resid '+str(j))
- nearby_molecules=[] #溶剂质心与guest的任一原子的距离小于5埃
- for molecule in new_position.split('residue'):
- solvent_centroid = molecule.center_of_mass()
- distances = mda_distances.distance_array(u.select_atoms('resid '+str(j)).positions, solvent_centroid)
- if any([i <= cluster_central for i in distances]):
- nearby_molecules.append(molecule)
- #选中溶剂原子相距范围
- selection_2 = u.select_atoms('resid '+str(j)) + u.select_atoms('around '+str(cluster_range) +' resid '+str(j))
- unique_residues_2 = set(atom.resid for atom in selection_2)
- unique_selection_2 = u.select_atoms(f"resid {' '.join(map(str, unique_residues_2))}")
- # Translate the coordinates
- new_position_2 = transcell(unique_selection_2, cell_dims)
- for molecule in new_position_2.split('residue'):
- nearby_molecules.append(molecule)
-
- #选中的分子的信息
- atoms_name = []
- coord = []
- for q in nearby_molecules:
- for t in q.names:
- atoms_name.append(t)
- for p in q.positions:
- coord.append(p)
- atom_names_2d = [[name] for name in atoms_name]
- merged_array = np.hstack((atom_names_2d, coord))
- merged_array = [list(i) for i in merged_array]
- merged_array = remove_duplicates(merged_array)
- # merged_array = transcell(np.array(merged_array), cell_dims)
- np.savetxt(str(file_name)+'_'+str(i).zfill(5)+'_'+str(j).zfill(3)+'.xyz',np.array(merged_array),\
- fmt='%-6s %12s %12s %12s',header=str(len(merged_array))+'\n',comments='')
- u1 = mda.Universe('A-toluene-npt.gro','A-toluene-npt.xtc')
- num_frame = len(u1.trajectory)
- pbc_cluster('A-toluene','A-toluene-npt.gro','A-toluene-npt.xtc',[num_frame-500,num_frame],[16,17],6,4)
复制代码
|
评分 Rate
-
查看全部评分 View all ratings
|