如何使用matlab绘制晶胞结构示意图
今天分享的内容是如何运用matlab编程实现晶胞结构示意图绘制。
一、前期准备工作（晶体基础知识+晶胞结构参数获取）
- 1 什么是晶体？
晶体（crystal）即是物质的质点（分子、原子、离子）在三维空间作有规律的周期性重复排列所形成的物质。
- 2 晶体长什么样？
- 3 晶体晶胞结构
晶胞分为元胞和超胞，元胞是指晶体结构最小重复单元，超胞是原胞的扩展，比如3*3*3，2*2*2等，扩展而成新的重复单元。
- 4 晶体结构参数
元胞是晶体结构最小重复单元，如果确定描述元胞的所有特性参数，我们就可以准确的绘制元胞，再通过元胞复制扩展，从而实现任意尺寸晶胞的绘制。
那么，元胞完整特性参数包括如下图所示，一个完整晶胞包括两部分：晶胞骨架+晶胞中分布原子（分子或离子）。
晶胞骨架可以通过6个晶胞参数进行完整描述，分别是晶胞三个维度方向长度a，b，c，以及相邻维度方向之间夹角alpha，beta，gamma。
晶胞中原子位置由晶体点群对称性决定，例如NaCl型晶胞属于Fm3m点群，Cl－构成面心立方，Na+填充在所有的正八面体空隙，在氯化钠晶体中,每个氯离子的周围都有6个钠离子,每个钠离子的周围也有6个氯离子。对于多数晶体，准确描述晶体中每一个原子的具体位置是一件很复杂事情。
- 5 晶体结构信息获取
为了方便获取晶体的结构信息，在这里给大家分享两个重要网站。这里Crystallographic Information File，简称CIF。
a、 CCDC-剑桥晶体数据库（CSD）
剑桥晶体数据中心（The Cambridge Crystallographic Data Centre ，CCDC)编制、发行的The Cambridge Structural Database (CSD)是一个管理有效、内容丰富的小分子有机物和金属有机化合物晶体结构的世界仓储，目前包含从X-射线到中子衍射分析的近百万个条目。
CCDC-剑桥晶体数据库（CSD）网址 https://www.ccdc.cam.ac.uk/
b、 ICSD-无机晶体数据库（ICSD）
ICSD（Inorganic Crystal Structure Database）是世界上最大的结构完全确定的无机晶体结构数据库，由德国FIZ Karlsruhe出版发行，每年春、秋季各更新一次。
ICSD-无机晶体数据库（ICSD）网址 https://icsd.fiz-karlsruhe.de/
- 6 CIF文件读取
简单认识一下CIF文件中内容，CIF文件可以直接用TXT记事本打开，CIF文件中记录晶体结构信息主要包括以下几个部分：晶胞参数，对称等效位置，简并原子位置。这里需要注意一点，需要将简并原子位置参数代入到对称等效位置，获得所有原子位置，再根据对应晶体点群类型，做进一步位置变换，删除其中超出元胞的原子位置，就可以获得晶胞中所有原子的位置参数，结合晶胞结构参数，就可以绘制晶体结构示意图了。
二、晶胞绘制数学模型搭建
Matlab绘制晶胞结构示意图的整体逻辑框架如下图所示。
在coding之前，我们需要解决两个数学问题：
- 如何基于晶胞参数确定晶胞骨架的8个顶点坐标？
- 如何通过元胞扩展实现超晶胞的绘制？
对于问题1，实际上我们需要解决一个如下数学问题，在一个三维直角坐标系下，已知三个相交于一点的三个线段，以及两两线段之间夹角，如何计算其中一条线段与另两条线段组成平面之间夹角。
问题可以简化成下图几何模型，构建三维直角坐标系，使线段a在x轴上，线段b在xy平面上，a与b夹角为alpha，线段c在xy平面以外，分别与a，b夹角为beta，gamma。基于此求红框所示六面体8个顶点1-2-3-4-4‘-3’-2‘-1’坐标。
首先，可以显然写出1-2-3-4四个顶点坐标，
顶点1坐标：[0, 0, 0];
顶点2坐标：[a, 0, 0];
顶点3坐标：[b*cos(alpha), b*sin(alpha), 0];
顶点4坐标：[a+b*cos(alpha), b*sin(alpha), 0];
只需求解1‘坐标，我们就可以基于对称平移思路，计算2’-3‘-4’坐标。
1‘坐标可以建立以下方程进行求解，该方程可以通过matlab求解实现。
假设1‘坐标为[x,y,z],则有
[x y z] *[a 0 0]/(a*c)=cos(beta);
[x y z]*[b*cos(alpha),b*sin(alpha),0]/(b*c)=cos(gamma);
x^2+y^2+z^2=c^2;
考虑到1‘在z轴上方，增加限制条件z>0;
解决了问题1，对于问题2，就相对简单，对于任意坐标沿a,b,c三个方向平移单位晶胞长度，其坐标变换分别为[a, 0 ,0],[0 b*cos(beta),0],[0 0 z]。
三，matlab函数代码
第一步，撰写CIF文件读取函数代码；函数名cifread，核心思路是基于特色字符识别关键信息读取所需参数。
function [cellpara,position,site,spacegroup]=cifread(filename)
fidin=fopen(filename);
fidout(1)=fopen('position.txt','w'); % 创建MKMATLAB.txt文件
fidout(2)=fopen('cellpara.txt','w'); % 创建MKMATLAB.txt文件
fidout(3)=fopen('site.txt','w'); % 创建MKMATLAB.txt文件
linecount=0;headerline=0;
while ~feof(fidin)
tline=fgetl(fidin); % 从文件读行
linecount=linecount+1;
while isempty(tline)
% tline=fgetl(fidin); % 从文件读行
% linecount=linecount+1;
break;
end
if contains(tline,'x')&&contains(tline,'y')&&contains(tline,'z')
k=1;
fprintf(fidout(k),'%sn',tline); %写入文本;
end
%%读取cellpara
if contains(tline,'_cell_length_a')
cellpara(1,1)=str2double(replace(tline,['_cell_length_a','(',")",' '],''));
k=2;
fprintf(fidout(k),'%sn',tline); %写入文本;
elseif contains(tline,'_cell_length_b')
cellpara(2,1)=str2double(replace(tline,['_cell_length_b','(',")",' '],''));
k=2;
fprintf(fidout(k),'%sn',tline); %写入文本;
elseif contains(tline,'_cell_length_c')
cellpara(3,1)=str2double(replace(tline,['_cell_length_c','(',")",' '],''));
k=2;
fprintf(fidout(k),'%sn',tline); %写入文本;
elseif contains(tline,'_cell_angle_alpha')
cellpara(4,1)=str2double(replace(tline,['_cell_angle_alpha','(',")",' '],''));
k=2;
fprintf(fidout(k),'%sn',tline); %写入文本;
elseif contains(tline,'_cell_angle_beta')
cellpara(5,1)=str2double(replace(tline,['_cell_angle_beta','(',")",' '],''));
k=2;
fprintf(fidout(k),'%sn',tline); %写入文本;
elseif contains(tline,'_cell_angle_gamma')
cellpara(6,1)=str2double(replace(tline,['_cell_angle_gamma','(',")",' '],''));
k=2;
fprintf(fidout(k),'%sn',tline); %写入文本;
end
%% 读取siteposition
if contains(tline,'_symmetry_space_group_name_H-M')
spacegroup=replace(tline,{'_symmetry_space_group_name_H-M',' ','''',},'');
end
if contains(tline,'_atom_site_')&&~contains(tline,'_geom_')
headerline=1;
elseif contains(tline,'loop_')
headerline=0;
end
if headerline~=0
k=3;
fprintf(fidout(k),'%sn',tline); %写入文本;
end
end
% cellpara1=importdata('cellpara.txt');% 将生成的data.txt文件导入工作空间，变量名为data，实际上它不显示出来
position1=importdata('position.txt');
site1=importdata('site.txt');site=struct;
pos=position1(2:end);pos1=pos;position=cell(length(pos(:,1)),3);
for i=1:length(pos(:,1))
if contains(pos(i,1),'''')
[startindex,endindex] = regexpi(pos(i,1),'''');
if startindex{1,1}(1)+1<=endindex{1,1}(end)-1
pos1{i,1}=pos{i,1}(startindex{1,1}(1)+1:endindex{1,1}(end)-1);
end
end
if contains(pos1(i,1),',')
position(i,:)=strsplit(pos1{i,1},',');
end
end
if isstruct(site1)
if length(site1.data(1,:))>=3
site.data=site1.data(:,1:3);
else
site.data=str2double(site1.textdata(:,3:5));
end
j=1;
for i=1:length(site1.textdata(:,1))
if ~isempty(site1.textdata{i,2})
site.textdata(j,1)=site1.textdata(i,2);
site.textdata{j,1}=replace(site.textdata{j,1},{' ','1','2','3','4','5','6','7','8','+','-'},'');
j=j+1;
end
end
site.data(isnan(site.data(:,1)),:)=[];
else
j=0;
for i=1:length(site1(:,1))
if ~contains(site1{i,1},'_atom_site_')
j=j+1;
txt=strsplit(site1{i,1});
site.textdata{j,1}=txt{1,2};
site.textdata{j,1}=replace(site.textdata{j,1},{' ','1','2','3','4','5','6','7','8','+','-'},'');
txt=replace(txt,{'(',')'},'');
site.data{j,1}=[str2double(txt{1,5}),str2double(txt{1,6}),str2double(txt{1,7})];
end
end
site.data=cell2mat(site.data);
end
fclose(fidin);第二步，根据晶胞参数计算晶胞骨架顶点；函数名zpoint，主要通过matlab求解三元二次方程组。
function [solx,soly,solz]=zpoint(cellpara)
syms x y z
eqns=[x^2+y^2+z^2==cellpara(3)^2,...
[x y z]*[cellpara(1) 0 0]'/cellpara(3)*cellpara(1)==cos(cellpara(5)/180*pi),...
[x y z]*[cellpara(2)*cos(cellpara(4)/180*pi) cellpara(2)*sin(cellpara(4)/180*pi) 0]'...
/cellpara(3)*cellpara(2)==cos(cellpara(6)/180*pi),z>0];
vars=[x y z];
[solx,soly,solz] = solve(eqns, vars);
end
第三步，基于晶胞对称等效位置和简并位置参数，计算晶胞中每个原子位置。函数名siteposition，主要通过晶体结构参数，对称等效位置及简并原子位置确定元胞中每个原子的位置坐标。
function siteposition = sitepositionfun(cellpara,position,site,spacegroup)
%UNTITLED 此处显示有关此函数的摘要
% 此处显示详细说明
sitepos=cell(length(site.data(:,1)),2);
for k=1:length(site.data(:,1))
x=site.data(k,1);y=site.data(k,2);z=site.data(k,3);
[l,d]=size(position);
% posvalue=zeros(l,d);
for i=1:l
for j=1:d
sitepos{k,2}(i,j)=eval(position{i,j});
end
end
sitepos(k,1)=site.textdata(k,1);
sitepos{k,2}=unique(sitepos{k,2},'row');
end
for k=1:length(site.data(:,1))
switch spacegroup(1)
case 'C'
sitepos{k,2}=sitepos{k,2}(:,:)+[0.5 0.5 0];
case 'P'
case 'A'
sitepos{k,2}(:,:)=sitepos{k,2}(:,:)+[0 0.5 0.5];
case 'I'
sitepos{k,2}(:,:)=sitepos{k,2}(:,:)+[0 0.5 0.5];
case 'R'
sitepos{k,2}(1:length(sitepos{k,2}(:,1)),:)=sitepos{k,2}(:,:)+[2/3 1/3 1/3];
sitepos{k,2}(length(sitepos{k,2}(:,1))+1:2*length(sitepos{k,2}(:,1)),:)=sitepos{k,2}(:,:)+[1/3 2/3 2/3];
end
[xx,yy,zz]=zpoint(cellpara);
deltax=1;
deltay=sin(cellpara(4)/180*pi);
deltaz=eval(zz)/cellpara(3);
xyz(1,:)=[eval(xx)/cellpara(1),eval(xx)/cellpara(1)+1];
xyz(2,:)=[eval(yy)/cellpara(2),eval(yy)/cellpara(2)+1];
xyz(3,:)=[0,1];
[m,n]=size(sitepos{k,2}(:,:));
for i=1:m
for j=1:n
if sitepos{k,2}(i,j)>xyz(j,2)||sitepos{k,2}(i,j)<xyz(j,1)
sitepos{k,2}(:,:)=sitepos{k,2}(:,:)-floor(sitepos{k,2}(:,:));
end
end
if isequal(sitepos{k,2}(i,:),[0 0 0])
sitepos{k,2}=[sitepos{k,2};[deltax deltay deltaz];[deltax deltay sitepos{k,2}(i,3)];[deltax sitepos{k,2}(i,2) deltaz];[sitepos{k,2}(i,1) deltay deltaz];...
[deltax sitepos{k,2}(i,2) sitepos{k,2}(i,3)];[sitepos{k,2}(i,1) deltay sitepos{k,2}(i,3)];[sitepos{k,2}(i,1) sitepos{k,2}(i,2) deltaz]];
elseif sitepos{k,2}(i,1)==0&&all(sitepos{k,2}(i,[2 3])~=0)
sitepos{k,2}=[sitepos{k,2};[deltax sitepos{k,2}(i,2) sitepos{k,2}(i,3)]];
elseif sitepos{k,2}(i,2)==0&&all(sitepos{k,2}(i,[1 3])~=0)
sitepos{k,2}=[sitepos{k,2};[sitepos{k,2}(i,1) deltay sitepos{k,2}(i,3)]];
elseif sitepos{k,2}(i,3)==0&&all(sitepos{k,2}(i,[1 2])~=0)
sitepos{k,2}=[sitepos{k,2};[sitepos{k,2}(i,1) sitepos{k,2}(i,2) deltaz]];
elseif isequal(sitepos{k,2}(i,[1 2]),[0 0])&&sitepos{k,2}(i,3)~=0
sitepos{k,2}=[sitepos{k,2};[deltax sitepos{k,2}(i,2) sitepos{k,2}(i,3)];[sitepos{k,2}(i,1) deltay sitepos{k,2}(i,3)];[deltax deltay sitepos{k,2}(i,3)]];
elseif isequal(sitepos{k,2}(i,[1 3]),[0 0])&&sitepos{k,2}(i,2)~=0
sitepos{k,2}=[sitepos{k,2};[deltax sitepos{k,2}(i,2) sitepos{k,2}(i,3)];[sitepos{k,2}(i,1) sitepos{k,2}(i,2) deltaz];[deltax sitepos{k,2}(i,2) deltaz]];
elseif isequal(sitepos{k,2}(i,[2 3]),[0 0])&&sitepos{k,2}(i,1)~=0
sitepos{k,2}=[sitepos{k,2};[sitepos{k,2}(i,2) deltay sitepos{k,2}(i,3)];[sitepos{k,2}(i,1) sitepos{k,2}(i,2) deltaz];[sitepos{k,2}(i,1) deltay deltaz]];
end
end
sitepos{k,2}=uniquetol(sitepos{k,2},'Byrows',true);
[m,~]=size(sitepos{k,2}(:,:));
for i=1:m
sitepos{k,2}(i,1)= sitepos{k,2}(i,1)*cellpara(1)+xx*sitepos{k,2}(i,3);
sitepos{k,2}(i,2)= sitepos{k,2}(i,2)*cellpara(2)*sin(cellpara(4)/180*pi)+yy*sitepos{k,2}(i,3);
sitepos{k,2}(i,3)= sitepos{k,2}(i,3)*zz;
end
end
siteposition=sitepos;
end
第四步，绘制元胞，根据8个顶点绘制晶胞12条框线，根据元胞中所有原子位置绘制原子。函数名plotunitcell，主要应用matlab中line函数和sphere函数。
function plotunitcell(siteposition,cellpara)
%PLOTATOM 此处显示有关此函数的摘要
% 此处显示详细说明
r=[0.4 0.5 0.6 0.6];color=['b','g','r','r'];
atoms=length(siteposition(:,1));
for i=1:atoms
for j=1:length(siteposition{i,2}(:,1))
plotsphere(siteposition{i,2}(j,:),r(i),color(i));
hold on
end
end
light;
set(gca,'xtick',[],'xticklabel',[]);
set(gca,'ytick',[],'yticklabel',[]);
set(gca,'ztick',[],'zticklabel',[]);
axis equal;
axis off
[zx,zy,zz]=zpoint(cellpara);
%%晶体元胞角点
point=zeros(8,3);
point(1,:)=[0 0 0];
point(2,:)=[cellpara(1),0,0];
point(3,:)=[cellpara(2)*cos(cellpara(4)/180*pi) cellpara(2)*sin(cellpara(4)/180*pi) 0];
point(4,:)=[cellpara(1)+cellpara(2)*cos(cellpara(4)/180*pi) cellpara(2)*sin(cellpara(4)/180*pi) 0];
point(5,:)=[zx zy zz];
point(6,:)=point(2,:)+point(5,:);
point(7,:)=point(3,:)+point(5,:);
point(8,:)=point(4,:)+point(5,:);
line(point([1,2],1),point([1,2],2),point([1,2],3),'LineWidth',2,'Color','k');
hold on
line(point([1,3],1),point([1,3],2),point([1,3],3),'LineWidth',2,'Color','k');
hold on
line(point([1,5],1),point([1,5],2),point([1,5],3),'LineWidth',2,'Color','k');
hold on
line(point([5,6],1),point([5,6],2),point([5,6],3),'LineWidth',2,'Color','k');
hold on
line(point([5,7],1),point([5,7],2),point([5,7],3),'LineWidth',2,'Color','k');
hold on
line(point([3,4],1),point([3,4],2),point([3,4],3),'LineWidth',2,'Color','k');
hold on
line(point([3,7],1),point([3,7],2),point([3,7],3),'LineWidth',2,'Color','k');
hold on
line(point([2,4],1),point([2,4],2),point([2,4],3),'LineWidth',2,'Color','k');
hold on
line(point([2,6],1),point([2,6],2),point([2,6],3),'LineWidth',2,'Color','k');
hold on
line(point([6,8],1),point([6,8],2),point([6,8],3),'LineWidth',2,'Color','k');
hold on
line(point([7,8],1),point([7,8],2),point([7,8],3),'LineWidth',2,'Color','k');
hold on
line(point([4,8],1),point([4,8],2),point([4,8],3),'LineWidth',2,'Color','k');
hold on
function plotsphere(pos,r,color)
[u,v,w]=sphere(100);
x=pos(1)+r.*u;
y=pos(2)+r.*v;
z=pos(3)+r.*w;
surf(x,y,z,'FaceColor',color,'EdgeColor','none','FaceAlpha',1);
colormap hot;
grid off;
end
end
第五步，绘制超晶胞，计算超晶胞顶点，计算超晶胞所有原子位置，重复第四步绘制超晶胞。
function plotsupercell(supercell,siteposition,cellpara,down,up)
%UNTITLED4 此处显示有关此函数的摘要
% 此处显示详细说明
% Check number of input arguments
global linecolor linewidth
if isempty(linecolor)
linecolor='k';
end
if isempty(linewidth)
linewidth=2;
end
try
narginchk(3,5);
catch ME
error(message('plotsupercell:numberOfInputs', ME.message));
end
if nargin < 5, up = [0 0 0];
if nargin < 4, down = up + supercell;
end
end
down=[min(down(1),up(1)),min(down(2),up(2)),min(down(3),up(3))];
up=down+supercell;
[dx,dy,deltaz]=zpoint(cellpara);
deltax=cellpara(1);
deltay=cellpara(2)*sin(cellpara(4)*pi/180);
pos=siteposition;
for i0=1:length(siteposition(:,1))
posnum=length(siteposition{i0,2}(:,1));plus=0;
for i1=1:supercell(1)
for i2=1:supercell(2)
for i3=1:supercell(3)
pos{i0,2}(plus*posnum+1:plus*posnum+posnum,:)=siteposition{i0,2}+[(down(1)+i1-1)*deltax+...
(down(2)+i2-1)*cellpara(2)*cos(cellpara(4)/180*pi)+(down(3)+i3-1)*dx,...
(down(2)+i2-1)*deltay+(down(3)+i3-1)*dy,(down(3)+i3-1)*deltaz];
plus=plus+1;
end
end
end
pos{i0,2}=uniquetol(pos{i0,2},'Byrows',true);
end
r=[0.4 0.5 0.6];color=['b','g','r'];%%原子半径及颜色
atoms=length(pos(:,1));
for i=1:atoms
for j=1:length(pos{i,2}(:,1))
plotsphere(pos{i,2}(j,:),r(i),color(i));
hold on
end
end
light;
set(gca,'xtick',[],'xticklabel',[]);
set(gca,'ytick',[],'yticklabel',[]);
set(gca,'ztick',[],'zticklabel',[]);
axis equal;
axis off
% supercellpara=cellpara;
% supercellpara(1:3)=cellpara(1:3)'.*supercell;
% [zx,zy,zz]=zpoint(supercellpara);
%%晶体元胞角点
point=zeros(8,3);
point(1,:)=[cellpara(1)*down(1)+cellpara(2)*cos(cellpara(4)/180*pi)*down(2) cellpara(2)*sin(cellpara(4)/180*pi)*down(2) 0]...
+[dx,dy,deltaz].*down(3);
point(2,:)=point(1,:)+[cellpara(1)*supercell(1),0,0];
point(3,:)=point(1,:)+[cellpara(2)*cos(cellpara(4)/180*pi)*supercell(2) cellpara(2)*sin(cellpara(4)/180*pi)*supercell(2) 0];
point(4,:)=point(1,:)+[cellpara(1)*supercell(1)+cellpara(2)*cos(cellpara(4)/180*pi)*supercell(2) cellpara(2)*sin(cellpara(4)/180*pi)*supercell(2) 0];
point(5,:)=point(1,:)+[dx dy deltaz]*supercell(3);
point(6,:)=point(2,:)+point(5,:)-point(1,:);
point(7,:)=point(3,:)+point(5,:)-point(1,:);
point(8,:)=point(4,:)+point(5,:)-point(1,:);
line(point([1,2],1),point([1,2],2),point([1,2],3),'LineWidth',linewidth,'Color',linecolor);
hold on
line(point([1,3],1),point([1,3],2),point([1,3],3),'LineWidth',linewidth,'Color',linecolor);
hold on
line(point([1,5],1),point([1,5],2),point([1,5],3),'LineWidth',linewidth,'Color',linecolor);
hold on
line(point([5,6],1),point([5,6],2),point([5,6],3),'LineWidth',linewidth,'Color',linecolor);
hold on
line(point([5,7],1),point([5,7],2),point([5,7],3),'LineWidth',linewidth,'Color',linecolor);
hold on
line(point([3,4],1),point([3,4],2),point([3,4],3),'LineWidth',linewidth,'Color',linecolor);
hold on
line(point([3,7],1),point([3,7],2),point([3,7],3),'LineWidth',linewidth,'Color',linecolor);
hold on
line(point([2,4],1),point([2,4],2),point([2,4],3),'LineWidth',linewidth,'Color',linecolor);
hold on
line(point([2,6],1),point([2,6],2),point([2,6],3),'LineWidth',linewidth,'Color',linecolor);
hold on
line(point([6,8],1),point([6,8],2),point([6,8],3),'LineWidth',linewidth,'Color',linecolor);
hold on
line(point([7,8],1),point([7,8],2),point([7,8],3),'LineWidth',linewidth,'Color',linecolor);
hold on
line(point([4,8],1),point([4,8],2),point([4,8],3),'LineWidth',linewidth,'Color',linecolor);
hold on
function plotsphere(pos,r,color)
[u,v,w]=sphere(100);
x=pos(1)+r.*u;
y=pos(2)+r.*v;
z=pos(3)+r.*w;
surf(x,y,z,'FaceColor',color,'EdgeColor','none','FaceAlpha',1);
colormap hot;
grid off;
end
end
四，代码测试及演示
最后，我们通过调用各个函数实现晶胞结构示意图绘制。
clc;clear;
global linewidth linecolor;
linewidth=2;linecolor=[0 0 0];
[cellpara,position,site,spacegroup]=cifread('[muchong.com]LiCoO2.cif');
siteposition = sitepositionfun(cellpara,position,site,spacegroup);
figure(1)
plotunitcell(siteposition,cellpara);
figure(2)
supercell=[2 2 2];
plotsupercell(supercell,siteposition,cellpara);最后展示一下，锂离子电池常见层状LiCoO2晶体结构示意图。
本图文内容来源于网友网络收集整理提供，作为学习参考使用，版权属于原作者。
THE END
二维码