Heinrich 发表于 2025-3-7 16:19
你好。我们的工作能做出来这样的分布。没有施加力出现速度的情况要看具体的模拟系统，比如孔内是否有电荷 ...

好的，非常感谢您的回复

13734436496 发表于 2024-12-26 15:45
请问，在孔隙中通过施加力的方式真的可以得到这样的水分子速度分布吗？为什么在不加力的时候速度会无限接近 ...

你好。我们的工作能做出来这样的分布。没有施加力出现速度的情况要看具体的模拟系统，比如孔内是否有电荷、是否有跨膜电位差之类的。

请问，在孔隙中通过施加力的方式真的可以得到这样的水分子速度分布吗？为什么在不加力的时候速度会无限接近0m/s？我在计算时发现即使没有施加力水分子依然存在一定速度。

Heinrich 发表于 2024-9-3 13:32
这个是我研究室的同事写的，仅供参考。当时我们刚做这个没有经验，用matlab写的。现在建议用Python写，效 ...

非常谢谢老师！！！

abing 发表于 2024-7-25 21:51
老师可以求教一下你的统计代码嘛

这个是我研究室的同事写的，仅供参考。当时我们刚做这个没有经验，用matlab写的。现在建议用Python写，效率高一点。

1. %% from trajectory file to plot the z-velocity profile in y-direction
   
2. % Che Wei Ou in Daiguji & Hsu lab, The University of Tokyo
   
3. % gromacs traj output the y-coordinate and z-velocity
   
4. % 20230123 Fix the bug so now can set the range over 4~8 in y direction
   
5. %% Reset
   
6. clear variables
   
7. clc
   
8. format longG
   
9. v_file = 'veloc-h2o-md2-10ns-20ns.xvg';
   
10. c_file = 'coordinate-h2o-md2-10ns-20ns.xvg';
    
11. %% Open file
    
12. fid_c = fopen(c_file); fid_v = fopen(v_file);
    
13. pat = ['#','@']; TF_v = 1;TF_c =1;
    
14. %% Skip headline until the first useful data is read
    
15. while TF_v == 1
    
16.     tline_v = fgetl(fid_v);
    
17.     TF_v = (contains(tline_v,pat(1))+contains(tline_v,pat(2)));
    
18. end
    
19. while TF_c == 1
    
20.     tline_c = fgetl(fid_c);
    
21.     TF_c = (contains(tline_c,pat(1))+contains(tline_c,pat(2)));
    
22. end
    
23. 
    
24. %% Accumulate the average velocity in each section
    
25. edges = 3.5:0.04:8.5; % y-direction
    
26. j= 1;
    
27. tic
    
28. while ischar(tline_c)
    
29.     temp_c = str2num(tline_c);
    
30.     temp_v = str2num(tline_v);
    
31.     disp(num2str(temp_c(1)))
    
32.     temp_c(1) = [];temp_v(1) = []; % delete the frame number
    
33.     [~, ~, bins] = histcounts(temp_c, edges);
    
34.     valididx = reshape(find(bins), [], 1);
    
35.     B = accumarray( reshape(bins(valididx), [], 1), valididx, [], @(V){V.'});
    
36.     for i= 1:length(edges)-1
    
37.         if i <= length(B)
    
38.     %         sum_veloc(j,i)=sum(temp_v(B{i}));
    
39.             distr_veloc(j,i)=mean(temp_v(B{i}));
    
40.         end
    
41.     end
    
42. %     number_ions(j,:)= bins;
    
43.     clear temp_c temp_v bins valididx B
    
44.     tline_v = fgetl(fid_v);
    
45.     tline_c = fgetl(fid_c);
    
46.     j = j+1;
    
47.    
    
48. end
    
49. toc
    
50. fclose(fid_v); fclose(fid_c);% close file
    
51. distr_veloc(find(isnan(distr_veloc)==1)) = 0;
    
52. %% plot
    
53. DL = length(abs(mean(distr_veloc)))+1;
    
54. plot((edges(1:DL-1)+edges(2:DL))/2,abs(mean(distr_veloc))*1000,'o--')
    
55. xlabel('y [nm]')
    
56. ylabel('Absolute Velocity [m/s]') % 10^-9 10^-12
    
57. title('Average Velocity Profile')
    
58. save('velocity.mat',"distr_veloc","edges")
    
59. %  for i = 1:length(edges)-1
    
60. %      each_slice_ions(i) = sum(sum(number_ions == i));
    
61. %  end
    
62. %  hold on
    
63. %  plot((edges(1:end-1)+edges(2:end))/2,sum(sum_veloc)./each_slice_ions,'r')
    

复制代码

Heinrich 发表于 2024-2-4 16:58
我用的工具是matlab。目前有两种计算方法，一种是用gmx traj从.trr中导出水分子的坐标和速度，用matlab一 ...

老师可以求教一下你的统计代码嘛

秃头研究生 发表于 2024-1-22 15:11
可以请教一下脚本怎么写吗？

我用的工具是matlab。目前有两种计算方法，一种是用gmx traj从.trr中导出水分子的坐标和速度，用matlab一行一行读取输出的数据；然后将y坐标切成很多个slides，根据水分子的坐标把这些分子放进不同的slides里，然后分别计算每个slide中水分子的平均速度。第二个方法是只导出水分子的坐标，然后通过计算两帧之间的分子位移确定速度，其他部分和第一种方法相同。

Heinrich 发表于 2022-2-18 08:43
好的，明白了。谢谢sob老师！

可以请教一下脚本怎么写吗？

Heinrich 发表于 2022-5-1 15:22
sob老师好。我在用Gromacs做电渗流模拟的时候发现一个现象：如果我在.mdp中设置“constraints=none”的话 ...

没特殊理由绝对不要用all-bonds，论坛里我说过很多次。

sobereva 发表于 2022-2-18 08:08
自己写脚本分析
设置nstvout让trr里记录速度信息，利用gmx traj可以从中提取特定原子的速度。利用脚本循环 ...

sob老师好。我在用Gromacs做电渗流模拟的时候发现一个现象：如果我在.mdp中设置“constraints=none”的话就不会出现电渗流，如果设置“constraints=all-bonds”就会出现，请问应该如何理解？我读了您之前写的《解析gromacs的restraint、constraint和freeze》(http://sobereva.com/10)，但还是不太理解。

sobereva 发表于 2022-2-18 08:08
自己写脚本分析
设置nstvout让trr里记录速度信息，利用gmx traj可以从中提取特定原子的速度。利用脚本循环 ...

好的，明白了。谢谢sob老师！

自己写脚本分析
设置nstvout让trr里记录速度信息，利用gmx traj可以从中提取特定原子的速度。利用脚本循环各帧、循环各层提取速度并统计速度。