请问如何利用VMD统计水在聚合物的扩散过程中的结合水和自由水分子的个数

LLLYYHHH

如题，请问是否可以利用VMD统计水在聚合物的扩散过程中的结合水和自由水分子的个数。
如果可以，能否请教一下大致的思路，谢谢各位老师！

对抗路达摩

首先你要对结合水和自由水下一个定义，MD的分析很多都是基于轨迹的，轨迹的基本单位是帧。首先你要定义一下什么是结合水，比如距离聚合物分子的距离小于一定的cut-off并且维持一定的时间n帧，有了定义才可以对轨迹进行结合水、自由水的分析。
然后有了定义就根据具体的定义写分析脚本，如果你对于分析脚本是0基础，我建议使用python结合库MDAnalysis，python自带非常多的分析库，个人认为是最容易上手的，MDAnalysis的功能也很强大，比如按照我的定义，可以使用MDAnalysis的动态选取功能，每一帧考察聚合物周围距离多少帧内的水都是哪些，再和上一帧的进行对比，最终就能分析出结合水和自由水的个数。

sobereva

用VMD需要自己写脚本分析，至少据我所致没有现成的脚本

LLLYYHHH

对抗路达摩 发表于 2022-12-27 20:11
首先你要对结合水和自由水下一个定义，MD的分析很多都是基于轨迹的，轨迹的基本单位是帧。首先你要定义一下 ...

谢谢您的建议~

LLLYYHHH

sobereva 发表于 2022-12-28 03:18
用VMD需要自己写脚本分析，至少据我所致没有现成的脚本

谢谢sob老师~

紫丹渡津

随手写了一个Python的脚本，并行化计算结合水和自由水的分数，而且你还可以自由扩展计算结合水的状态，比如与聚合物链结合的数目等。
如何使用，脚本里关键的地方都有注释。

1. import multiprocessing as mp
   
2. import functools
   
3. 
   
4. import math
   
5. import MDAnalysis as mda
   
6. import MDAnalysis.analysis.distances as dist
   
7. import MDAnalysis.lib.distances as distances
   
8. import numpy as np
   
9. 
   
10. def assoPrint(filename, data1, data2):  # print association properties
    
11.     anaout = open(filename + '.xvg', 'w')
    
12.     print('# ' + 'condensed' + ' free', file=anaout)
    
13.     for i in range(0, len(data1)):
    
14.         print('{:10.5f} {:10.5f}'.format(data1[i], data2[i]), file=anaout)
    
15.     anaout.close()
    
16. 
    
17. 
    
18. def get_association(top, trj, ag1, ag2, cutoff, frame_id):
    
19.     uta = mda.Universe(top, trj)
    
20.     group1 = uta.select_atoms("type " + ag1)
    
21.     group2 = uta.select_atoms("type " + ag2)
    
22.     ts = uta.trajectory[frame_id]
    
23.     cell = ts.dimensions
    
24.     type1, type2 = [], []
    
25.     for atom in range(len(ag1)):
    
26.         distpair = distances.capped_distance(ag1.positions[atom], ag2.positions, cutoff, box = cell)[0]
    
27.         asso_pair.append(len(distpair))
    
28.         if len(distpair) != 0:
    
29.             type1.append(atom)
    
30.         else:
    
31.             type2.append(atom)
    
32.     # return the fractions of condensed and free water molecules
    
33.     return len(type1) / len(group1), len(type2)  / len(group1)
    
34. 
    
35. 
    
36. def main():
    
37.     top = "topol.tpr"
    
38.     trj = "traj.trr"
    
39.     atom_group1 = "AT1" # atom type for the oxygen of water
    
40.     atom_group2 = "AT2" # atom type of the given reference on your polymer
    
41.     cutoff = 5.0  # cutoff for water - polymer
    
42.     nt = 8  # number of processors
    
43.     uta = mda.Universe(top, trj)
    
44.     # generate the wanted frame_ids
    
45.     frame_ids = [ts.frame for ts in uta.trajectory]
    
46.     condensed, free = [], []
    
47.     do_analyse = functools.partial(get_association, top, trj, atom_group1, atom_group2, cutoff)
    
48.     pool = mp.Pool(processes=nt)
    
49.     data = pool.imap(do_analyse, frame_ids)
    
50.     for f1, f2 in data:
    
51.         condensed.append(f1)
    
52.         free.append(f2)
    
53.     assoPrint('water_fraction', condensed, free)
    
54. 
    
55. 
    
56. if __name__ == "__main__":
    
57.     main()

复制代码

LLLYYHHH

紫丹渡津 发表于 2022-12-29 14:21
随手写了一个Python的脚本，并行化计算结合水和自由水的分数，而且你还可以自由扩展计算结合水的状态，比如 ...

谢谢老师！我这边也写了一个tcl的脚本，您的程序我再仔细学习学习

光下澈

LLLYYHHH 发表于 2022-12-29 15:30
谢谢老师！我这边也写了一个tcl的脚本，您的程序我再仔细学习学习

你好，请教一下，有成功计算出来吗？