分子动力学吸附识别MATALB脚本，可以实现A物质吸附在B物质上的识别判断

%-----------------------------------------------% % Author：lammps_materials@163.com % %-----------------------------------------------% % first we need to read dump file file ="diffusion.dump"; t_start_all = cputime; try dump = fopen(file,'r'); catch error('Dumpfile not found!'); end i=1; while feof(dump) == 0 id = fgetl(dump); if (strncmpi(id,'ITEM: TIMESTEP',numel('ITEM: TIMESTEP'))) timestep(i) = str2num(fgetl(dump)); else if (strncmpi(id,'ITEM: NUMBER OF ATOMS',numel('ITEM: NUMBER OF ATOMS'))) Natoms(i) = str2num(fgetl(dump)); else if (strncmpi(id,'ITEM: BOX BOUNDS',numel('ITEM: BOX BOUNDS'))) x_bound(i,:) = str2num(fgetl(dump)); y_bound(i,:) = str2num(fgetl(dump)); z_bound(i,:) = str2num(fgetl(dump)); else if (strcmpi(id(1:11),'ITEM: ATOMS')) for j = 1 : 1: Natoms atom_data(j,:,i) = str2num(fgetl(dump)); end i=i+1; end end end end end t_start_stop = (cputime-t_start_all)/60; fprintf('Now the used time is: %.1f mins.\n',t_start_stop); disp("-------------------"); disp("----ALL DONE!!!----"); disp("-------------------"); % 这个是自动关机- % system('shutdown -s'); %---------------------------------------------------------% %{ %---------------------------------------------------------% 1.We have 4 types: 1,2,3,4. 2.We want to know how many atoms (type 4) are absorbed onto adsorbent(type 2,3) during simulation. 3.Here we define the quantitative parameter as: the number of absorbed atom/total number of type3 %---------------------------------------------------------% Here the boundary is p p p conditions in MD, We choise the cutoff is 1nm, which indicate that distance between captured atom and adsorbent is less than 1nm. %---------------------------------------------------------% %} all_frame = 100; cutoff = 1.5; %xl,yl,zl %%----------------------------------------------%% xl = x_bound(1,2)-x_bound(1,1); yl = y_bound(1,2)-y_bound(1,1); zl = z_bound(1,2)-z_bound(1,1); %%----------------------------------------------%% Nlocal = length(atom_data); % %------------------------------------------------% water_cap = zeros(all_frame,1); water_free_tag = zeros(Nlocal,all_frame); %------------------------------------------------% [fid,message] = fopen("test.txt","wt"); write_data = zeros(Nlocal,6); % loop all the frame t=cputime;disp('Calculating...'); for frame = 1:100 t_start = cputime; if(rem(frame,20)==0) disp('%------------------------------------------------%'); fprintf('Now the calculating frame is: %d.\n',frame); end % define xyz type %------------------------------------------------% now_frame = atom_data(:,:,frame); ID = now_frame(:,1); TYPE = now_frame(:,2); XYZ = now_frame(:,3:5); %------------------------------------------------% check_capture_water=0; check_cap3 = 0; in_loop_i=0; out_mark_i=0; check_in=0; check_34=0; %------------------------------------------------% % put your code here! % %------------------------------------------------% % Select start! for i = 1:Nlocal %------------------------------------------------% if(TYPE(i)==4) check_34 = check_34+1; %------------------------------------------------% for j = 1:Nlocal dx = XYZ(j,1)-XYZ(i,1); dy = XYZ(j,2)-XYZ(i,2); dz = XYZ(j,3)-XYZ(i,3); PBC(dx,dy,dz,xl,yl,zl); % pbc distance = sqrt(dx*dx+dy*dy+dz*dz); %------------------------------------------------% if(distance<=cutoff && (TYPE(j)==3)) check_in = check_in+1; in_loop_i = i; if(check_in==1) out_mark_i =i; end if(in_loop_i~=out_mark_i ) check_cap3 =check_cap3+1; check_capture_water = check_capture_water+1; water_free_tag(i,frame)=1; end out_mark_i = i; end end % loop all atom j %------------------------------------------------% end end % loop all atom i % Select finish! %------------------------------------------------% % put your code here! % %------------------------------------------------% water_cap(frame)=check_capture_water;% save %------------------------------------------------% %写出dump轨迹 write_data(:,1) = ID; write_data(:,2) = TYPE; write_data(:,3:5) = XYZ; write_data(:,6) = water_free_tag(:,frame); % %------------------------------------------------% title_first = "ITEM: TIMESTEP"; title_number = "ITEM: NUMBER OF ATOMS"; title_boundary = "ITEM: BOX BOUNDS pp pp pp"; title_all = "ITEM: ATOMS id type x y z watertag"; fprintf(fid,'%s\n',title_first); fprintf(fid,'%d\n',timestep(frame)); fprintf(fid,'%s\n',title_number); fprintf(fid,'%d\n',Natoms(1)); fprintf(fid,'%s\n',title_boundary); fprintf(fid,'%d %d\n',[x_bound(1,1),x_bound(1,2)]); fprintf(fid,'%d %d\n',[y_bound(1,1),y_bound(1,2)]); fprintf(fid,'%d %d\n',[z_bound(1,1),z_bound(1,2)]); fprintf(fid,'%s\n',title_all); [r,c]=size(write_data(:,:)); for i=1:r for j=1:c fprintf(fid,'%f\t',write_data(i,j)); end fprintf(fid,'\r\n'); end %------------------------------------------------% t_frame=cputime-t_start; if(rem(frame,20)==0) fprintf('Now the frame %d is writing in file...\n',frame); fprintf('Now the frame %d is finished.\n',frame); fprintf('Now the used time of current frame is: %8.1f s.\n',t_frame); disp('%------------------------------------------------%'); end end t_all=(cputime-t)/60; fprintf('Now the used time is: %8.1f mins.\n',t_all); disp("-------------------"); disp("----ALL DONE!!!----"); disp("-------------------");