MATLAB计算平衡晶格常数

材料的平衡晶格常数是比较常见和基础的性质,在开始入门计算材料时一般都会接触到。

一般的方法是计算晶体在不同晶格常数下的结合能,然后通过对结合能进行拟合,找到能量的最小值,对应的就是平衡晶格常数。

这里忽略具体如何获取这些“晶格常数——结合能”数值对,而介绍如何对这些数值对进行处理,获取平衡晶格常数。

比较常见的方法是使用origin的多项式拟合功能,获得多项式系数;然后再对对多项式进行处理,获得多项式的极值点对应的晶格常数就是平衡晶格常数。我个人感觉,还是有点小麻烦,所以就写了下面的代码,可以直接处理得到平衡晶格常数和相应的结合能。具体代码如下:

function cal_lat_const(N,inFileName)
% calcuate the lattice constant according to "lat_const cohesive_energy"
% Input:
%   N: order of the polynomial fitting.
%   inFileNmae: name of the file storing "lat_const cohesive_energy"
% Example:
% 	cal_lat_const(8,'Au')
%	Here, 'Au' is a file stores "lat_const cohesive_energy"
% Powered by Xianbao Duan
% Email: xianbao.d@gmail.com
% Website: http://www.52souji.net/

% read in data from the file
data = load(inFileName,'-ascii');
x = data(:,1);
y = data(:,2);

Ecoh = polyfit(x,y,N);    	% polynomial fitting
dEcoh = polyder(Ecoh);  	% derivation of the polynomial equation
zero_points = roots(dEcoh);    	% solve the zero points

% calculate the lattice constant and coresponding cohesive energy
for i = 1: length(zero_points)
    if isreal(zero_points(i))
        if zero_points(i) > x(1)
            if zero_points(i) < x(end)
                lat_const = zero_points(i)
                coh_energy = spline(x,y,lat_const)
            end
        end
    end
end

% plot the curve
xx = x(1):0.1:x(end);
yy = polyval(Ecoh,xx);
figure;
plot(x,y,'o',xx,yy);
xlabel('lattice constant / A');
ylabel('cohesive energy / eV');
legend('Original points','Fitted curve');
说明
1. 第20行是对拟合的多项式求零点。一般会有多个零点,所以zero_points为一数组;
2. 第27,28行会直接把平衡晶格常数和相应的结合能输出,但因为可能会有多个,所以需要你根据平衡晶格常数的大致范围,自己做出选择;
3. 输入参数N为多项式拟合的阶数,一般取 6~8。需要根据具体情况进行调整,并不是越大越好;
4. 第35-41行是绘制拟合的曲线。

下面以Fe为例具体说明。
我先用LAMMPS算了不同晶格常数下,Fe的结合能,存为文件 fe.txt,数据如下:

2.1	18.26577944
2.2	6.940458875
2.3	0.09168609
2.4	-2.348192601
2.5	-3.180605085
2.6	-3.665433307
2.7	-3.9828838
2.8	-4.133455473
2.9	-4.141436809
3	-4.031262224
3.1	-3.870968947
3.2	-3.657669589
3.3	-3.415075988
3.4	-3.163401596
3.5	-2.89362391
3.6	-2.602716638
3.7	-2.302031334
3.8	-2.008237816
3.9	-1.742084327
4	-1.506097434
4.1	-1.295680398
4.2	-1.104512187
4.3	-0.923554689
4.4	-0.748623841
4.5	-0.583055534
4.6	-0.437166093
4.7	-0.259909843

在MATLAB命令行输入:

>> cal_lat_const(8,'fe.txt')

得到如下结果(但一般情况下会输出多组):

lat_const =
    2.8962
coh_energy =
   -4.1436

同时,绘制出来的图像如下图所示:

lattice_constant_vs_cohesive_energy_of_fe

可以看出,拟合的还是非常不错的。

如果你有更好的方法,欢迎与我分享!

标签: matlab

相关文章推荐

添加新评论 (无需注册,可直接评论)

已有 4 条评论

  1. LSW

    您好,老师,为什么我按照您的方法想拟合别的晶格常数,就是更改了输入数据文件,但是通过修改N,但是始终没有拟合出好的曲线,问问什么原因。谢谢

  2. Domicor

    你好,我在你的函数的基础上进行了一些优化。
    我使用你的函数时发现作图有一些问题,调试后发现是最后plot时的步长问题。我使用的数据(lattice_02_processed.data)和函数(calLatticeConstant.m)已经放在github:
    https://github.com/Rareson/LammpsRelated
    欢迎讨论。

    Domicor.

    1. 陈尊友

      你好,能不能把你优化的函数给我一份。万分感谢!!
      我的邮箱:cy2173@qq.com

    2. 我爱搜集网博主

      非常谢谢你的反馈和改进! :)