Learn VASP The Hard Way (Ex35): Birch-Murnaghan方程拟合脚本-2-脚本精解篇

2017-08-30


  Ex35 Birch-Murnaghan方程拟合脚本-2


脚本精解篇


阅读本文前,为避免看不懂,请先浏览前面两节的内容!

 

Ex34 Birch-Murnaghan方程拟合脚本-1


http://www.bigbrosci.com/newsitem/278011645



Ex33 晶格参数的确定


http://www.bigbrosci.com/newsitem/278011633


通过上一节(Ex34)的脚本,我们可以拟合出稳定晶体结构的晶格参数,但是,这个脚本是怎么运行的呢?大家在浏览之后,有没有什么疑问?本节我们重点介绍一下:这个脚本的运行原理。

 



1  数据

 

首先我们看下读取的data文件:



A) 共两列,两列之间用tab分开;

B) 第一列为缩放系数,第二列为对应结构的能量。

C) 在本文例子中,晶格参数等于缩放系数 * 2.8664 Å (记住这个公式)

 



2  脚本内容: 浏览一遍即可,详解在下面: 





详解:


2.1  1-15行:准备部分



A) 1-4 行为脚本的一些基本说明,为养成一个良好的写脚本习惯,以及避免后期不必要的麻烦,第一行必须要放进去。如果你用bash写的话,把python替换为bash即可,其他部分不用动,例如: #!/usr/bin/env bash

B2-3行注释,解释了一下脚本干什么用的,谁写的,怎么用。

C)  6-7 行调用python库中的math  numpy 模块,见下面两个链接:

 

https://docs.python.org/2/library/math.html

http://www.numpy.org/

 

D9-15行:通过使用 print 命令输出一个警告,因为本文用的是缩放系数,很多人测试的时候保持缩放系数不变,直接采用不同的晶格参数。

 



2.2  读取data信息,准备拟合



A)  16-18行为脚本的注释说明,

B)  20行中,我们使用 np.loadtxt读取data文件的内容,

注意:

1)  data 用单引号括起来; 

2)   usecols=(0,1) 代表我们要读取第一列和第二列的内容。Python中的第一个用0表示,所以这里我们用的(0,1)来表示读取第1和第2列的内容;

3)   delimiter='\t',意思是第一列和第二列之间用tab'\t')分开;

4)    unpack=True 意思是把读取的数据生成一个数列,并分配给 E,具体见链接:

 

https://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html

  

通过第20行,我们可以获得两个数列,第一个是 a 包含了data文件中的缩放系数, 第二个是E,包含了data文件中的能量信息

C) 22-28行,把缩放系数转换为晶格常数(a*2.8664),然后平方取倒数(**(-2)),用来得到我们前面讲解的x。这是因为:当我们把BM方程写成:y(x) = c0 + c1x+ c2x2+c3x3 的形式时,前提是令 (1/a)2 = x。不懂的具体见Ex33中的内容。

 

http://www.bigbrosci.com/newsitem/278011633

 

截止到现在,我们有了三个数列: a x

 



2.3 29-40行:拟合多项式方程: 



A) 最小二乘法拟合y(x) = c0 + c1x+ c2x2+c3x3这个方程,得到常数系列c,这里我们调用了numpy polyfit这个模块来进行拟合。具体见:

 

https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html

 

B) 31行:本脚本的核心: p = np.polyfit(x, E, 3)

 

1) np.polyfit 读取x E的数值,拟合一个y(x) = c0 + c1x+ c2x2+c3x3的方程,

 

注意:这里拟合的时候,用的数值是xE(1/a)2 = x 所以,在进行拟合之前,务必确保你已经得到了正确的x数列。具体根据data文件的内容进行修改设置24行的公式

 

2) 括号中的3代表的是多项式方程中最大的那个次数。

扩展:如果你想拟合一个 y = c4x4+ …的方程,把3改成即可。

3)拟合的结果输出到数列 p 中,数列p中共有4个常数(p[0],  p[1]p[2]p[3]),按照次数从高到底的顺序排列,对应着我们Ex33中的(c3c2, c1,  c0 

4) 33-36我们将中的数值分别分配给c3c2, c1,  c0

 



2.4   38-40行,获取晶格参数:

 

到目前位置,多项式方程的常数我们已经获得了,下面就是求解问题了。


A)  42行为下图中公式的代码描述,获取方程的解x,因为前面我们已经有了一个x变量了,这里方程的解用x1表示,

B) 43行获取晶格常数,公式为: (1/a)2= x

C) 45行,输出晶格常数的数值。

 



3 扩展练习:

 

3.1 根据解释,重新浏览脚本,知道每一行代表的含义。

3.2 本脚本使用的时候,需要自己根据data数列中的内容,设置更改一下获取 x 数列的公式(24);

3.3 阅读Ex32-本节的所有内容,整个流程不要再存在什么疑问。

 



 4 总结

 

本节我们详细讲解了BM方程拟合脚本,相信大家应该了解的差不多了,拟合的办法有很多,但BM方程只有一个。请务必掌握至少一个拟合的方法。此外,本人马上就要毕业答辩了,后面更新会有所放缓,努力争取保证每周一篇的速度。如果你有什么科研心得,欢迎分享。

 


如有错误,疑问,或者建议,请发邮件联系大师兄: qli@bigbrosci.com

给大师兄留言:QQ号(2674006510)  微信: BigBroSci

加入大师兄QQ群:遇见大师兄 217821116 (先根据群公告修改群名片,若不修改,一经发现,立即踢出。

如果你有自己的科研经验和心得,也欢迎分享给大家!

此外,QQ群专注于科研思维的碰撞与科研生活的分享,本书中已经详细解释或者指明的易出错部分,不建议在群中继续咨询,请大家认真学习并主动积极地去思考和练习。


如果喜欢大师兄的文章,欢迎关注我们,转载,转发。

打赏一下,鼓励大师兄们写出更好的文章!

本网站由阿里云提供云计算及安全服务