matlab dump,matlab 读取lammps dump 文件

发一个我修改过的m文件吧：

% MATLAB后处理 LAMMPS数据：求解管道流动，各量沿着管道半径方向的分布：

function [varargout] = scandump(varargin) % Function to scan LAMMPS dump file

clc;

clear all

try

dump = fopen('water_pos+vel.txt','r');

catch

error('Dumpfile not found!');

end

%定义初始的一些量：

R=30;     %定义管道的半径

long=75;   %管道的长度

h=input('请输入每一层的厚度：h=');    %每一层的厚度

nlayer=floor(R/h+0.5); %沿着管道半径方向，所分的层数

i=1;

out.position(i,1) = 0;

fid = fopen('okokokokokokokok.txt','w'); % 把处理好的数据存储到density_ok.txt

fprintf(fid,'# r  density/number ave_vel\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%% 读取并处理LAMMPS输出文件部分%%%%%%%%%%%%%%%%%%%%%%%%

while ~feof(dump);   %判断是否到了文件的结尾处，如果没有到文件的结尾处则：

id = fgetl(dump);  %从文件中读取一行，赋值给id

for k=1:nlayer;

r(k)=(k-0.5)*h;

v(k)=pi*h*long*(2*r(k)+h)/1000; % r(k)是沿着管道半径方向，第k层的半径值

count(k)=0;

sumv(k)=0;

end

switch id

case 'ITEM: TIMESTEP'

out.timestep(i,1) = str2num(fgetl(dump)) %这种情况下，下一行的值就是timestep

fprintf(fid,'TIMESTEP          Layer\n');    %把这些信息内容输出到处理之后的文件里：

fprintf(fid,'%10d    %10d \n',out.timestep(i,1),nlayer);

case 'ITEM:        NUMBER        OF        ATOMS                OF        ATOMS'

out.Natoms(i,1) = str2num(fgetl(dump));  %这种情况下，下一行的值就是原子的个数，也就是说原子信息部分的列数

case 'ITEM:        BOX        BOUNDS                        BOUNDS'

out.boxbound(1,:,i) = str2num(fgetl(dump)); %下一行是x的范围

out.boxbound(2,:,i) = str2num(fgetl(dump)); %y的范围

out.boxbound(3,:,i) = str2num(fgetl(dump)); %z的范围

case 'ITEM:        ATOMS        id        type        x        id        type        x'

for j = 1

ut.Natoms(i,1);

A=str2num(fgetl(dump)); % 读取所有原子信息的部分，并且做相关的统计处理

x(j)=A(:,3);           % LAMMPS输出的dump文件里面，第七列是每个原子的r值

y(j)=A(:,4);

vz1(j)=A(:,8);         %vz的值，vz是沿着通道方向的速度

vz(j)=vz1(j)*20.45482706*100; %把速度的值转换成 nm/ns

v_r(j)=sqrt(x(j)*x(j)+y(j)*y(j));

if (v_r(j)<=30.5)

jlayer=floor(v_r(j)/h); %jlayer 就是第j层

if (jlayer <= 0)

jlayer = 1;

end

if (jlayer > nlayer)

jlayer = nlayer;

end

% 设置每一层的计算器初值

layer(j)=jlayer;

count(jlayer)= count(jlayer)+1;

sumv(jlayer)=sumv(jlayer)+vz(jlayer)

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%% 后处理的结果存入最终的数据文件中 %%%%%%%%%%%%%%%%%%%%

for k=1:nlayer;

density_number(k)=count(k)/v(k);

ave_vel(k)=sumv(k)/count(k);

fprintf(fid,'%5d %10f %10d %10f \n ',r(k),density_number(k),ave_vel(k));

end

end

end

fclose(fid); %存储结束之后关闭文件