计算材料学
0%

发表于 分类于
本文字数： 1.4k 阅读时长 ≈ 1 分钟 
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
from pymatgen.core.surface import Slab, SlabGenerator, generate_all_slabs, Structure, Lattice, ReconstructionGenerator, get_symmetrically_distinct_miller_indices
from pymatgen.core.structure import Structure
from pymatgen.analysis.adsorption import *
from pymatgen.ext.matproj import MPRester
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.io.vasp.inputs import Poscar
from matplotlib import pyplot as plt
import os

mpr = MPRester('Bw7HdCARiXvzEWJK')
path = os.getcwd()

# 需要修改
mp_id = "mp-153" 
# struct = Structure.from_file("Y3Mg.vasp")
struct = mpr.get_structure_by_material_id(mp_id)
struct = SpacegroupAnalyzer(struct).get_conventional_standard_structure()

formula = struct.composition.reduced_formula
print("输入的结构是:{}".format(formula))

ls = get_symmetrically_distinct_miller_indices(struct, 1, return_hkil=False)

# os.mkdir(formula)
# os.chdir(formula)
for i in ls:
    slab = SlabGenerator(struct, miller_index=i, min_slab_size=25,min_vacuum_size=15.0, lll_reduce=True, center_slab=True)
    for n, slabs in enumerate(slab.get_slabs(bonds=None, ftol=0.1, tol=0.1, max_broken_bonds=0, symmetrize=True, repair=False)):
        slabs.make_supercell([[2,0,0],[0,2,0],[0,0,1]])
    #     open(formula + mp_id[3:] + '-' + str(n+1) + '.vasp', 'w').write(str(Poscar(slabs)))
        name = str(i).split(',')[0][1]+str(i).split(',')[1][1]+str(i).split(',')[2][1]
        open(formula+'_'+name +'_' + str(n+1) + '.vasp', 'w').write(str(Poscar(slabs)))
# os.chdir(path)s

发表于 更新于 分类于
本文字数： 2.5k 阅读时长 ≈ 2 分钟 
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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
# 20200422添加mp库修正原子能量
eles = {'La':-4.9266, 'Ce':-5.9343, 'Pr':-4.7837, 'Y':-6.4674,
        'Er':-4.5674, 'Li':-1.9080, 'Sm':-4.7163, 'Ni':-5.7783,
        'Cu':-4.0993, 'Ag':-2.8262, 'Sr':-1.6841, 'Eu':-10.2906,
        'Hg':-0.3016, 'Tm':-4.4782, 'Ho':-4.5781, 'Tb':-4.6275,
        'Lu':-4.5259, 'Dy':-4.5999, 'Ba':-1.9246, 'Nd':-4.7651, 
        'Pr':-4.7837, 'F':-1.882, 'Sc' : -6.3325,'H':-3.23975,
        'He': -0.0255} 

mu_mg = -1.594

import os

ele = str(input('请输入另一种元素符号(注意大小写格式):'))
index_a =eval(input('请输入化学式中Mg的原子个数：'))
index_b =eval(input('请输入化学式中另一种元素的原子个数：'))
formation_energy=eval(input('请输入化合物每原子的生成焓：'))
mu_formula_atom=eval(input('请输入化合物每原子的DFT能量：'))
print('----------计算结果输出如下----------')

'''
total_num:单胞中总共的原子个数
mu_formula：单胞DFT修正后能量
ratio：用来判断切好的表面中元素个数的关系，要减去多余的不符合化学计量比的原子
mu_F：在富Mg情况下第二种元素的原子化学势
low：假设的一个表面能初始值，1000.用于后面找到最小的表面
'''

total_num =index_a+index_b
mu_formula =mu_formula_atom*total_num
ratio = index_b/index_a
low = 1000.0
mu_F = eles[ele]+formation_energy*total_num/index_b

#获取当前目录
path = os.getcwd()
dirs = os.listdir(path)
#进入到表面目录进行循环
for dir in dirs:
    if dir[0]=='0' or dir[0]=='1':
        os.chdir(path+'/'+dir)
        temp_poscar = open("POSCAR", "r", encoding="utf-8").readlines()
        poscar = []
        for i in temp_poscar:
            poscar.append(i.strip())
        oszicar = open("OSZICAR", "r", encoding="utf-8").readlines()
        x=poscar[5][0:-1].split(" ")[0]
        y=poscar[5][0:-1].split(" ")[-1]
        #获取元素的系数，x是Mg
        if x[0] == 'M':
            num_x=eval(poscar[6].strip().split()[0])
            num_y=eval(poscar[6].strip().split()[1])
        else:
            num_x=eval(poscar[6].strip().split()[1])
            num_y=eval(poscar[6].strip().split()[0])
        #获取表面积
        temp_a = list(filter(None, poscar[2].split(" ")))
        temp_b = list(filter(None, poscar[3].split(" ")))
        area = eval(temp_a[0])*eval(temp_b[1])
        #获取收敛能量，但目前无法判断结构是否是收敛
        try:
            energy = eval(oszicar[-1].strip().split()[4])
        except:
            energy = eval(oszicar[-1].split(" ")[7])
        n = 0
        m = 0
        #判断系数关系
        if (num_y/num_x) == ratio:
            surf = (energy - num_x/index_a*mu_formula)/2/area
        elif (num_y/num_x) < ratio:
            while (num_y/num_x) < ratio:
                n = n+1
                num_x = num_x-1
            surf = (energy - num_x/index_a*mu_formula - n*mu_mg)/2/area
        else:
            while (num_y/num_x) > ratio:
                num_y = num_y-1
                m = m+1
            surf = (energy - num_x/index_a*mu_formula - m*mu_F)/2/area
        #输出表面能
        convert_surf = surf*16000
#         print(os.getcwd())
#         print(energy)
        print("{}面,{:.3f} mj/m2".format(os.getcwd().split('/')[-1],convert_surf))
        if low > convert_surf:
            low = convert_surf
        else:
            continue
    else:
        continue
# -0.018
# -2.4328
print('----------计算结果为----------')
print('表面能计算结束,最稳定表面为：'+str(round(low,3)))

发表于 分类于
本文字数： 183 阅读时长 ≈ 1 分钟

注意正负数要多加验证，比如-0.5取整到-1还是0？

1
2
3
4
5
6
7
8
9
10
11
3**2 = 9 #次方运算
1/2 = 0.5 #正常除法
3%5 = 3 #取余数
int(11/5) = 2 #取整数部分，向0取整
round(3.5) = 4 # 四舍五入
11//3=3 #取整数 

#涉及到负数时，用以下方式
import math 
math.ceil(3.25)=4 #向上取整
math.floor() #向下取整

发表于 分类于
本文字数： 2.1k 阅读时长 ≈ 2 分钟 
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
44
45
46
47
48
49
# -*- coding:utf-8 -*-
import json
from pymatgen.apps.borg.hive import VaspToComputedEntryDrone
from pymatgen.apps.borg.queen import BorgQueen
from pymatgen.io.vasp.outputs import Outcar
from pymatgen.analysis.phase_diagram import PhaseDiagram
from pymatgen.analysis.phase_diagram import PDPlotter
from pymatgen.entries.compatibility import MaterialsProjectCompatibility
from pymatgen.ext.matproj import MPRester

para = ["is_hubbard", "hubbards","potcar_symbols", "run_type","structures","efermi","is_spin","potcar_spec","atomic_symbols"]#,"get_computed_entry"]#"potcar_spec"
dat = ["final_energy"]
etry = VaspToComputedEntryDrone(parameters = para,data = dat)
#entry.assimilate(a)#.to_json())
path = "/lustre/home/acct-umjzhh/umjzhh-2/zhenming/e-hull/"
queen = BorgQueen(etry,rootpath=path,number_of_drones=1)
#queen.serial_assimilate(path)
entries = queen.get_data()
#save_data = queen.save_data('entry_save.txt')
#entries = queen.load_data('entry_save.json')
#print("entries:",entries)
#print("enrty0:",entries[0])
#print("entries:",dir(entries[0]))
a=MPRester("mqzy3IaxSNegF941")
compat = MaterialsProjectCompatibility()
computed_entries = compat.process_entries(entries)
print("len_computed_entries:",len(computed_entries))
elements = [ ]
for entry in computed_entries:
    #print("entry.as_dict():",entry.as_dict()[u'parameters'][u'atomic_symbols'])
    #print("set(entry.as_dict():",set(entry.as_dict()[u'parameters'][u'atomic_symbols']))
    element = list(set(entry.as_dict()[u'parameters'][u'atomic_symbols']))
    #print("list(set(entry.as_dict():)",element)
    s = [ ]
    for i in range(len(element)):
        s.append(str(element[i]))
    elements.append(s)#json.dumps(#.get_reduced_composition_and_factor#[u'atomic_symbols']#composition.reduced_composition#.atomic_symbols
#print(elements)
ss=[]  
for element in elements:
    for i in range(len(element)):
        ss.append(element[i])
ss = list(set(ss))
mp_entries = a.get_entries_in_chemsys(ss)
total_entries = computed_entries+mp_entries
pd = PhaseDiagram(total_entries)
for i in range(len(computed_entries)):
    decomp,e_above_hull = pd.get_decomp_and_e_above_hull(computed_entries[i])
    print("name and ebh:",(computed_entries[i].name,e_above_hull))

发表于 分类于
本文字数： 1.2k 阅读时长 ≈ 1 分钟

简介

奇异值分解（singular Value Decomposition），简称SVD，线性代数中矩阵分解的方法。假如有一个矩阵A，对它进行奇异值分解，可以得到三个矩阵：

20141228153754275

这三个矩阵的大小：

20141228153852359

矩阵sigma(即上图U和V中间的矩阵)除了对角元素不为0，其他元素都为0，并且对角元素是从大到小排列的，前面的元素比较大，后面的很多元素接近0。这些对角元素就是奇异值。

sigma中有n个奇异值，但是由于排在后面的很多接近0，所以我们可以仅保留比较大的r个奇异值：

20141228155108176

实际应用中，我们仅需保留着三个比较小的矩阵，就能表示A，不仅节省存储量，在计算的时候更是减少了计算量。SVD在信息检索（隐性语义索引）、图像压缩、推荐系统、金融等领域都有应用。

在主成分分析中，我是通过特征值分解的方法来实现PCA的，除了特征值分解，还可以用奇异值分解来实现PCA。特征值和奇异值二者之间是有关系的：上面我们由矩阵A获得了奇异值sigma(i)，假如方阵A*A’的特征值为lamda(i)，则：sigma(i)^2=lamda(i)。可以发现，求特征值必须要求矩阵是方阵，而求奇异值对任意矩阵都可以，因此PCA的实现其实用SVD的更多，在scikit-learn中，PCA算法其实也是通过SVD来实现的。

python中实现SVD

numpy中的linalg已经实现了SVD，可以直接调用

1
2
3
4
5
6
7
8
9
10
11
12
>>> A=mat([[1,2,3],[4,5,6]])
>>> from numpy import linalg as la
>>> U,sigma,VT=la.svd(A)
>>> U
matrix([[-0.3863177 , -0.92236578],
        [-0.92236578,  0.3863177 ]])
>>> sigma
array([ 9.508032  ,  0.77286964])
>>> VT
matrix([[-0.42866713, -0.56630692, -0.7039467 ],
        [ 0.80596391,  0.11238241, -0.58119908],
        [ 0.40824829, -0.81649658,  0.40824829]])

有一点需要注意，sigma本来应该跟A矩阵的大小2*3一样，但linalg.svd()只返回了一个行向量的sigma，并且只有2个奇异值（本来应该有3个），这是因为第三个奇异值为0，舍弃掉了。之所以这样做，是因为当A是非常大的矩阵时，只返回奇异值可以节省很大的存储空间。当然，如果我们要重构A，就必须先将sigma转化为矩阵。

Matlab实现SVD

1
[u,s,v] = svd(f);

注意：Python中的SVD和matlab中的SVD中最后一个矩阵v互为转置关系

发表于 分类于
本文字数： 163 阅读时长 ≈ 1 分钟

读取

1
2
3
4
5
6
#适用于文本但不适用于数字
with open ('test.dat', 'r+') as f:
  f.read()

#pandas读取方法
df = pd.read_csv('test.dat',encoding = 'GBK')

存储

1
2
data1 = pd.DataFrame(arr1)
data1.to_csv('data1.csv')

发表于 更新于 分类于
本文字数： 744 阅读时长 ≈ 1 分钟 
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
# 逻辑回归 平方差画图
import numpy as np 
# mpl_toolkits是matplotlib官方的工具包 mplot3d是用来画三维图像的工具包
from mpl_toolkits.mplot3d import Axes3D
# pyplot 是一个有命令风格的的函数集合，与matlab相似。
from matplotlib import pyplot as plt
import math

# 创建一个图像窗口
fig = plt.figure()
# 在图像窗口添加3d坐标轴
ax = Axes3D(fig)

# 使用np.linspace定义 x:范围(-10,10);个数为100
x = np.linspace(-10,10,100)
# 定义 y:范围(-3,3);个数为50
y = np.linspace(-10,10,100)
# 创建x-y平面网络
x,y = np.meshgrid(x,y)
# 定义函数 r=1/2*(x-y)^2
# r = 1/2*np.square(x-y)
r = x**2+y**2-math.sqrt(2)*x*y

# 将函数显示为3d  rstride 和 cstride 代表 row(行)和column(列)的跨度 get_cmap为色图分类 
ax.plot_surface(x,y,r,rstride = 1, cstride = 1,cmap='rainbow')

# 投影
# ax.contourf(x,y,r,zdir= 'x', offset = (-2),cmap=plt.get_cmap('rainbow')) 

# 显示创建的图像
plt.savefig('minus.tif', dpi=200)
plt.show()

实际效果如图：

minus.tif

发表于 更新于 分类于
本文字数： 1k 阅读时长 ≈ 1 分钟 
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
44
45
46
47
48
49
50
51
52
53
54
55
# fix_atom_v1.0
# 仅针对于pymatgen生成的标准POSCAR，若其他格式POSCAR导致坐标位数不同需要更改脚本

import os

#Need Modifycation!!!
#对称表面上下可以动的层数，读取当前文件下POSCAR内容
# layers = eval(input("请输入slab上下要松弛的层数"))
layers = 3
f = open('POSCAR', 'r')
lines = f.readlines()
f.close()
new = []
for line in lines:
    new.append(line.strip())
os.rename('POSCAR','POSCAR_unfixed_atoms')

#加上Selective Dynamics变成可改变模式的
new.insert(7,'Selective Dynamics')

#去掉最后的字母元素符号
for i in range(9,len(new)):
    if new[i][-2] == ' ':
        new[i] = new[i][0:-2]
    else:
        new[i] = new[i][0:-3]

#获得slab层数及坐标
ls = []
for i in range(9,len(new)):
#     num = 1
#     while new[i][-num] != ' ':
#         num += 1
#     print(num-1)
#     print(i,new[i][1-num:])
    ls.append(eval(new[i][-8:]))
ls = list((set(ls)))
ls.sort()
print("此slab模型共有{}层, 松弛上下{}层".format(len(ls),layers))

#加上FT模式
for i in range(9,len(new)):
    if eval(new[i][-8:]) in ls[0:layers]:
        new[i] = new[i]+"  T T T"
    elif eval(new[i][-8:]) in ls[-layers:]:
        new[i] = new[i]+"  T T T"
    else:
        new[i] = new[i]+"  F F F"

#保存文件
f = open('POSCAR', 'w')
for para in new:
    f.write(para)
    f.write('\n')
f.close()