这篇文章是我在2017年12月开的坑,但是AFK了将近一年,本质原因还是因为懒(借口)。不过在2018年9月19日开始,我正式开始了实验室生涯,做稀土光功能材料与计算方向,就不得不重新来填这个坑了2333。稀土元素独特的4f电子层轨道,让其具有独特的光、电、磁特性,使得稀土元素具有“点石成金”的能力,在先端科技产品的应用中有着不可替代的作用。然而稀土元素由于其原子半径相近,化学性质相似,使得高纯度单一稀土元素难以分离,试剂纯级的稀土试剂价格极高,作为实验耗材来说其高昂的价格往往使课题组难以承受。同时,复杂的4f轨道是人类还未真正探索的量子领域之一,探索、表征以及实现对4f电子层的量化建模,必将对化学科研和量子力学领域产生极大的推动和影响。利用机器学习方法,对4f电子层轨道进行计算拟合,会成为我接下来要做的一些微小的工作之一。本文也将成为接下来这一年的科研笔记。


种一棵树的最佳时间是十年前,其次是现在。

以记录机器学习和计算化学的单元操作coding为主,防止出现版权问题和“基础不牢地动山摇”的情况。
以记录机器学习基础与量子化学基础理论为主,纯粹作为一个技术笔记。

数学参考:

重要配置——关闭win10自动更新

首先要写的是与机器学习无关的内容,因为我所使用的设备是Windows 10系统的辣鸡电脑,而这坑爹系统总容易因为蜜汁更新导致文件缺失,所以,为了解决这个问题就一定要在进行任何软件/环境配置之前关闭自动更新!
操作方法为:计算机(此电脑)-右键-管理-服务和应用程序-服务-Windows update-右键-属性-禁用-关闭

Windows10搭建python数据分析/机器学习环境

这里同理于Linux上搭建机器学习环境,下载安装Anaconda,为了方便和现有开源包的匹配,我选择python2.7版本,对应Anaconda是2.X版本。之后像常规安装软件一样安装就行了。

在安装过程中,可能会遇到

1
failed to create anacoda menue?

解决方案有:

  • 安装目录中不得出现中文或其他奇怪的符号
  • CMD中进入安装目录,然后执行 python .\Libs_nsis.py mkmenus
  • 将JAVA的环境变量JAVA_HOME和PATH先删掉,就可以顺利安装了。

安装好之后,其中Navigator是图形界面,可以切换python版本,同时能够在右侧搜索未安装的包,然而出于众所周知的原因,Anaconda官方的源下载速度基本上你懂的,所以这个时候建议切换国内的源。
添加Anaconda的TUNA镜像

1
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/

设置搜索时显示通道地址

1
conda config --set show_channel_urls yes

对于anaconda,建议直接使用其自带的IDE——Spyder,打开后,默认就加载了ipython和python两个解释器。对于conda命令无法使用的情况,在环境变量中添加将anaconda安装目录下的Scripts。对于CMD下不能使用anaconda同理,添加安装目录到环境变量即可。

关于机器学习环境,直接借助开源包即可,这里直接pass掉某显卡包,毕竟我的电脑还要自己用!!直接使用Scikit-Learn,Numpy ,Scipy ,Matplotlib。关于这个网上有很多安装教程,我直接用conda和pip安装,之后遇到需要再拓展的包,随用随装。

1
2
3
4
5
6
7
conda install scipy
conda install theano
conda install keras
pip install tensorflow
pip install sklearn
pip install numpy
pip install matplotlib

emmmm实际上上面有点多此一举了,貌似Anaconda已经集成了这些。总之环境搭建的重心还是放在Anaconda的安装上,其他不懂的就问Google吧。

Runge-Kutta 的python实现

龙格库塔法是工程上常用的微分方程的数值计算方法,因为在冶金、化工过程当中,很多简单模型(如极其简单的那种有边界条件的一维稳态计算)都可以化简为简单的4阶龙哥库塔法模型,所以使用计算机语言掌握这种方法也是很有必要的!!

而在计算机程序运作当中,节约时间成本是非!常!必!要!的!所以这种情况下一般使用构造-调用函数的方法来完成计算。下面具体分析:
首先,根据常微分方程,将其构造为一个自定义函数,如果你是一个常微分方程的计算,就构造一个;如果是常微分方程组就构造多个。python当中使用def语法来构造,比如我将方程f=y-(2*x/y)构造为fxy,代码如下:

1
2
3
def fxy(x,y):  
f = y - (2 * x / y)
return f

接下来就是构造一个迭代器,用于解常微分方程(组)。龙格库塔法的数学表达并不难,这里必须强调的几个点是:

  • x0,y0均为初值条件
  • h为步长,即微积分当中所谓的极小长度dx或dy
  • N为迭代数,迭代次数为N-n,一般来说迭代次数越多,输出值约精确。这是一种迭代逼近。

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
def rk4(x0,y0,h,N):  
n = 1
while(n != N):
x1 = x0 + h
k1 = fxy(x0, y0)
k2 = fxy(x0+h/2, y0+h*k1/2)
k3 = fxy(x0+h/2, y0+h*k2/2)
k4 = fxy(x1, y0+h*k3)
y1 = y0 + h * (k1 + 2* k2 + 2 * k3 + k4) / 6
print("%.2f, %.6f" %(x1, y1))
n = n + 1
x0 = x1
y0 = y1

接下来就是主要计算的部分了,为了方便设定,初值条件、步长、迭代数实现自设定,直接人工输入,之后调用上面写好的函数即可,这里你也可以把主要计算构造成一个函数,然后执行即可。代码如下:

1
2
3
4
5
6
7
def main():  
x0=float(input("x0= "))
y0=float(input("y0= "))
h=float(input("h= "))
N=float(input("N= "))
rk4(x0,y0,h,N)
main()

接下来可以举个例子来计算常微分方程组:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def fxyzw(x, y, z, w):
f = z * 0.1
return f
def gxyzw(x, y, z, w):
g = w * 0.1
return g
def pxyzw(x, y, z, w):
p = (x **2 + y ** 2 + z ** 2 + w ** 2) ** (1 / 2) * 0.1
return p
def rg4(x0,y0,z0,w0,h,N):
n = 1
while(n != N):
x1 = x0 + h
k1 = fxyzw(x0, y0, z0, w0)
l1 = gxyzw(x0, y0, z0, w0)
m1 = pxyzw(x0, y0, z0, w0)
k2 = fxyzw(x0+h/2, y0+h*k1/2, z0+h*l1/2, w0+h*m1/2)
l2 = gxyzw(x0+h/2, y0+h*k1/2, z0+h*l1/2, w0+h*m1/2)
m2 = pxyzw(x0+h/2, y0+h*k1/2, z0+h*l1/2, w0+h*m1/2)
k3 = fxyzw(x0+h/2, y0+h*k2/2, z0+h*l2/2, w0+h*m2/2)
l3 = gxyzw(x0+h/2, y0+h*k2/2, z0+h*l2/2, w0+h*m2/2)
m3 = pxyzw(x0+h/2, y0+h*k2/2, z0+h*l2/2, w0+h*m2/2)
k4 = fxyzw(x1, y0+h*k3, z0+h*l3, w0+h*m3)
l4 = gxyzw(x1, y0+h*k3, z0+h*l3, w0+h*m3)
m4 = pxyzw(x1, y0+h*k3, z0+h*l3, w0+h*m3)
y1 = y0 + h * (k1 + 2* k2 + 2 * k3 + k4) / 6
z1 = z0 + h * (l1 + 2* l2 + 2 * l3 + l4) / 6
w1 = w0 + h * (m1 + 2* m2 + 2 * m3 + m4) / 6
print("%.2f, %.6f, %.6f, %.6f" %(x1, y1, z1, w1))
n = n + 1
x0 = x1
y0 = y1
z0 = z1
w0 = w1
def main():
x0=float(input("x0= "))
y0=float(input("y0= "))
z0=float(input("z0= "))
w0=float(input("w0= "))
h=float(input("h= "))
N=float(input("N= "))
rg4(x0,y0,z0,w0,h,N)
main()

这里注意2点:

  • Python中有单目运算符+(正号)、-(负号)。双目运算符:+、 - 、 *、 /、 %、**、//,分别表示加、减、乘、除法、取余、求幂、整除。所以要用对符号!!
  • 类似一个多元常微分方程组,记得把程序里的变量数和元数匹配!!

经典机器学习方法

1.数据降维与回归方法

2.数据聚类与分类方法

争对上述2类方法,一般在实际数据处理时要进行补值(暂且认为在已做好降维的情况下),不同的方法有与之对应不同的补值方法。而特征选择方法一般有以下2种:

  • 集成式特征选择方法(排序过滤算法方法部分):BLogReg、SBMLR、T-test、Relief、Fisher、ChiSquare、InforGain、Gini
  • RPFS方法:RSM、PCA、Fisher得分、SFS

常用软件与数据库

对合成材料进行表征,之后一般需要计算机绘图,比如画一些二维的分子式、反应过程,或者三维的分子式。还有就是数据图的精修,比如XRD精修一类的。目前在做的主要有生物分子和无机材料2个方向,当然这里也有一部分是和有机材料通用的。

  • 生物分子

关于生物分子,一般是做分子对接,其步骤如下:
从数据库下载分子结构文件,然后用chemdraw软件绘制分子式,之后保存为pdb格式。使用网页程序Z-Dock进行结合位点的选择和绑定点残留的选择,并利用分子对接软件Accelrys Discovery Studio 3.5采用CDOCKER模块进行分子对接模拟。

常用数据库:RCSB Protein Data Bank

网页程序:Z-Dock

  • 无机材料

对于无机材料主要是绘制3D分子结构图,虽然chemdraw也可以画,但是crystalmaker更方便一些,而只要有了无机晶体的一些晶体参数即可画晶体结构,这类参数主要有:晶体常数、对称性以及位置坐标、原子半径,定义化学键等。可以使用数据编辑器来编辑它们。

而一般晶体文件可以从数据库里找,主要是cif文件。

常用的cif数据库有:The Cambridge Crystallographic Data CentreIUCR

常用软件:mercury
(实际上crystalmaker和chemdraw更好用)
还有一些文件是mol文件,这类文件的格式是:cas号.mol

基本上应用第一广泛的分子结构文件。因为它就是一段标准文本,可以用文本编辑器看。几乎所有的可视化分子结构软件都能读取mol文件并渲染其描述的结构。

常用mol数据库:chemexperchemicalbook

  • 有机材料

至于绘画什么的都按照分子式来,具体的和上面那些相同,可以下载现成的mol/cif文件,或者用chemdraw自己画。

最后,关于有机/生物材料,最大的问题是结构式/分子式很难获得,即使得到cif文件也难以去进行后期处理,这里推荐Jmol这个软件,获取mol或获取pdb都不错,而且可读DDL2 的CIF文件。

至于软件库,可以看看这里IUCr

pdb数据库与pymol的安装(for Windows)、蛋白质在线计算库

由于项目涉及到做蛋白质分子对接这一块儿,所以必须自己去搞相关工作,然而之前用jmol提取了cif文件的三维结构,不过显示的是球棒模型,而非蛋白质一般的带条模型,本来去问老板联系的另一个老师,结果她啥都不会,所以只能是我自己去搞了。

蛋白质分子对接的思路有重要的一点就是由cif文件提取pdb文件,之后才能做下一步工作:计算位点、对接。

因此,pdb数据库有时候也是急需的,可以节省大量时间。这里从小木虫上找了几个:

另外就是需要比较好的编辑工具,比较同源建模获得的pdb文件是需要优化一下H的,使得能量最低。这里所找到的最合适也最有普遍意义的就是pymol了,然而pymol在官网上找不到免费版本!!而它是开源软件!原则上是免费的!!因此只能找国外大学的源,这里就要用UCI的源了(国内需要梯子访问!!):

从这里找到pymol即可,可惜只有Windows版本的。由于我之前在电脑上装过Anconada2,所以是集成了python2.7环境的。就下载python2.7版本64位的。之后,在dos下输入

1
python

验证是打开了Anconada2中集成的python2.7,随后输入

1
quite()

退回到dos之后,输入

1
pip

显示已经装好pip了。由于下载的文件是whl文件,需要用到wheel来安装,所以

1
pip install wheel

安装之后,安装pymol,我下载的是pymol-2.1.0-cp27-cp27m-win_amd64.whl

1
pip install pymol-2.1.0-cp27-cp27m-win_amd64.whl

到这里就安装好了,接下来就能愉快的使用pymol了,输入

1
pymol

即可打开,可视化界面就不多说了,关于教程网上有很多,这里也不废话了。但是个人觉得如果单纯是用球棍模型,jmol就完全够用了。

而最后,由于蛋白质分子领域的计算很麻烦,装很多软件的操作实在是坑爹,这里参考这篇文章找到了一些在线计算的网站:

蛋白质理化性质分析:

点击Resources A…Z,从右侧选择要检测的项目即可,基本上需要的检测项目都有。