|
|
这个是我研究室的同事写的，仅供参考。当时我们刚做这个没有经验，用matlab写的。现在建议用Python写，效率高一点。
- %% from trajectory file to plot the z-velocity profile in y-direction
- % Che Wei Ou in Daiguji & Hsu lab, The University of Tokyo
- % gromacs traj output the y-coordinate and z-velocity
- % 20230123 Fix the bug so now can set the range over 4~8 in y direction
- %% Reset
- clear variables
- clc
- format longG
- v_file = 'veloc-h2o-md2-10ns-20ns.xvg';
- c_file = 'coordinate-h2o-md2-10ns-20ns.xvg';
- %% Open file
- fid_c = fopen(c_file); fid_v = fopen(v_file);
- pat = ['#','@']; TF_v = 1;TF_c =1;
- %% Skip headline until the first useful data is read
- while TF_v == 1
- tline_v = fgetl(fid_v);
- TF_v = (contains(tline_v,pat(1))+contains(tline_v,pat(2)));
- end
- while TF_c == 1
- tline_c = fgetl(fid_c);
- TF_c = (contains(tline_c,pat(1))+contains(tline_c,pat(2)));
- end
- %% Accumulate the average velocity in each section
- edges = 3.5:0.04:8.5; % y-direction
- j= 1;
- tic
- while ischar(tline_c)
- temp_c = str2num(tline_c);
- temp_v = str2num(tline_v);
- disp(num2str(temp_c(1)))
- temp_c(1) = [];temp_v(1) = []; % delete the frame number
- [~, ~, bins] = histcounts(temp_c, edges);
- valididx = reshape(find(bins), [], 1);
- B = accumarray( reshape(bins(valididx), [], 1), valididx, [], @(V){V.'});
- for i= 1:length(edges)-1
- if i <= length(B)
- % sum_veloc(j,i)=sum(temp_v(B{i}));
- distr_veloc(j,i)=mean(temp_v(B{i}));
- end
- end
- % number_ions(j,:)= bins;
- clear temp_c temp_v bins valididx B
- tline_v = fgetl(fid_v);
- tline_c = fgetl(fid_c);
- j = j+1;
-
- end
- toc
- fclose(fid_v); fclose(fid_c);% close file
- distr_veloc(find(isnan(distr_veloc)==1)) = 0;
- %% plot
- DL = length(abs(mean(distr_veloc)))+1;
- plot((edges(1:DL-1)+edges(2:DL))/2,abs(mean(distr_veloc))*1000,'o--')
- xlabel('y [nm]')
- ylabel('Absolute Velocity [m/s]') % 10^-9 10^-12
- title('Average Velocity Profile')
- save('velocity.mat',"distr_veloc","edges")
- % for i = 1:length(edges)-1
- % each_slice_ions(i) = sum(sum(number_ions == i));
- % end
- % hold on
- % plot((edges(1:end-1)+edges(2:end))/2,sum(sum_veloc)./each_slice_ions,'r')
|
|
发表于 Post on 2022-2-18 05:04:34
|
只看该作者 Only view this author
|只看大图 View big picture
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
发表于 Post on 2022-2-18 08:08:00
楼主 Author
