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
同时,绘制出来的图像如下图所示:
可以看出,拟合的还是非常不错的。
如果你有更好的方法,欢迎与我分享!
已有 4 条评论
您好,老师,为什么我按照您的方法想拟合别的晶格常数,就是更改了输入数据文件,但是通过修改N,但是始终没有拟合出好的曲线,问问什么原因。谢谢
你好,我在你的函数的基础上进行了一些优化。
我使用你的函数时发现作图有一些问题,调试后发现是最后plot时的步长问题。我使用的数据(lattice_02_processed.data)和函数(calLatticeConstant.m)已经放在github:
https://github.com/Rareson/LammpsRelated
欢迎讨论。
Domicor.
你好,能不能把你优化的函数给我一份。万分感谢!!
我的邮箱:cy2173@qq.com
非常谢谢你的反馈和改进! :)