计算化学公社

标题: 求助：如何计算纳米孔内水分子速度的空间分布 [打印本页]
作者
Author: Heinrich    时间: 2022-2-18 05:04
标题: 求助：如何计算纳米孔内水分子速度的空间分布
在电渗流相关文献中会见到如附件所示的纳米孔道内水分子速度的空间分布曲线，这篇文献和我使用的MD软件均为GROMACS。
请问这样的数据需要如何计算得到？ (, 下载次数 Times of downloads: 63)
我的系统是一个3nm长的纳米孔，分隔两个水箱。
作者
Author: sobereva    时间: 2022-2-18 08:08
自己写脚本分析
设置nstvout让trr里记录速度信息，利用gmx traj可以从中提取特定原子的速度。利用脚本循环各帧、循环各层提取速度并统计速度。
作者
Author: Heinrich    时间: 2022-2-18 08:43
sobereva 发表于 2022-2-18 08:08
自己写脚本分析
设置nstvout让trr里记录速度信息，利用gmx traj可以从中提取特定原子的速度。利用脚本循环 ...

好的，明白了。谢谢sob老师！
作者
Author: Heinrich    时间: 2022-5-1 15:22
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)，但还是不太理解。
作者
Author: sobereva    时间: 2022-5-1 16:19
Heinrich 发表于 2022-5-1 15:22
sob老师好。我在用Gromacs做电渗流模拟的时候发现一个现象：如果我在.mdp中设置“constraints=none”的话 ...

没特殊理由绝对不要用all-bonds，论坛里我说过很多次。
作者
Author: 秃头研究生    时间: 2024-1-22 15:11
Heinrich 发表于 2022-2-18 08:43
好的，明白了。谢谢sob老师！

可以请教一下脚本怎么写吗？
作者
Author: Heinrich    时间: 2024-2-4 16:58
秃头研究生 发表于 2024-1-22 15:11
可以请教一下脚本怎么写吗？

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

老师可以求教一下你的统计代码嘛
作者
Author: Heinrich    时间: 2024-9-3 13:32
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

  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')
复制代码
作者
Author: abing    时间: 2024-9-4 10:40
Heinrich 发表于 2024-9-3 13:32
这个是我研究室的同事写的，仅供参考。当时我们刚做这个没有经验，用matlab写的。现在建议用Python写，效 ...

非常谢谢老师！！！
作者
Author: 13734436496    时间: 2024-12-26 15:45
请问，在孔隙中通过施加力的方式真的可以得到这样的水分子速度分布吗？为什么在不加力的时候速度会无限接近0m/s？我在计算时发现即使没有施加力水分子依然存在一定速度。
作者
Author: Heinrich    时间: 2025-3-7 16:19
13734436496 发表于 2024-12-26 15:45
请问，在孔隙中通过施加力的方式真的可以得到这样的水分子速度分布吗？为什么在不加力的时候速度会无限接近 ...

你好。我们的工作能做出来这样的分布。没有施加力出现速度的情况要看具体的模拟系统，比如孔内是否有电荷、是否有跨膜电位差之类的。
作者
Author: 13734436496    时间: 2025-5-20 16:29
Heinrich 发表于 2025-3-7 16:19
你好。我们的工作能做出来这样的分布。没有施加力出现速度的情况要看具体的模拟系统，比如孔内是否有电荷 ...

好的，非常感谢您的回复



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3