处理分数占据结构

模拟方法
发布时间 :
字数:8.7k
  1. 1.pymatgen检验相同结构
    1. 1.一次性检验
    2. 2.逐个检验
    3. 3. 检验两个结构
  2. 2. pymatgen-enumeration
  3. 3. pymatgen-orderdisorder
  4. 4.计算静电能
  5. 5. sqs
    1. 1.前处理
      1. 1. 转换rndstr格式的结构
      2. 2.提交多个ip的任务
    2. 2. 后处理
      1. 1. 提取bestcorr
      2. 2.对多任务ip产生的best.out文件处理
      3. 3. 检验相同的结构
  6. 6.暴力枚举
    1. 1.枚举脚本
    2. 2. 多任务脚本 检验相同结构
    3. 3.改良版
      1. 1.限定最小距离
      2. 2.限定最小距离且保证替代发生在不同的位置
      3. 3.限定空间群
      4. 4.多任务做改良版
    4. 4.使用枚举方法处理蓝钻结构
      1. 1.完全随机的方法（完全解放双手了）
        1. 1. 提任务
        2. 2. changeatom.py
        3. 3. multitask.sh :做一些文件夹 ，提交产生结构的任务，提交计算ewald能量的任务
        4. 4. 后处理文件 get-all.sh 提取能量最低的结构
      2. 2. 替换为B_3O和B_4O构型的方法（还需要手动做一些流程）
        1. 1. 枚举氧的位置
        2. 2. 检验相同的结构
        3. 3. 替换B原子
        4. 4. 再次检验相同结构
        5. 5. 计算能量，整理结构
      3. 3. 替换为B_2O和B_3O的方法
        1. 1.枚举氧的结构和其他步骤
        2. 2-1. 替换B原子(联合取代、完全枚举)
        3. 2-2 先分类再替代
        4. 2-3 每个结构随机取样特定数目的结构
      4. 4.枚举+迭代器进行工作
      5. 5.枚举-使用静电能作为哈希值而不是fit
      6. 6. 可以分步骤来做，先替换四个，再替换两个 ！
    5. 5.检验
      1. 1. 逐个检验相似性
      2. 2. 检验两个结构
  7. 7.找出缺失的原子
  8. 7. 处理CaH6的333倍胞的问题（将324个H中的14个扣掉-可以枚举的方法

  1. 1.pymatgen检验相同结构
    1. 1.一次性检验
    2. 2.逐个检验
    3. 3. 检验两个结构
  2. 2. pymatgen-enumeration
  3. 3. pymatgen-orderdisorder
  4. 4.计算静电能
  5. 5. sqs
    1. 1.前处理
      1. 1. 转换rndstr格式的结构
      2. 2.提交多个ip的任务
    2. 2. 后处理
      1. 1. 提取bestcorr
      2. 2.对多任务ip产生的best.out文件处理
      3. 3. 检验相同的结构
  6. 6.暴力枚举
    1. 1.枚举脚本
    2. 2. 多任务脚本 检验相同结构
    3. 3.改良版
      1. 1.限定最小距离
      2. 2.限定最小距离且保证替代发生在不同的位置
      3. 3.限定空间群
      4. 4.多任务做改良版
    4. 4.使用枚举方法处理蓝钻结构
      1. 1.完全随机的方法（完全解放双手了）
        1. 1. 提任务
        2. 2. changeatom.py
        3. 3. multitask.sh :做一些文件夹 ，提交产生结构的任务，提交计算ewald能量的任务
        4. 4. 后处理文件 get-all.sh 提取能量最低的结构
      2. 2. 替换为B_3O和B_4O构型的方法（还需要手动做一些流程）
        1. 1. 枚举氧的位置
        2. 2. 检验相同的结构
        3. 3. 替换B原子
        4. 4. 再次检验相同结构
        5. 5. 计算能量，整理结构
      3. 3. 替换为B_2O和B_3O的方法
        1. 1.枚举氧的结构和其他步骤
        2. 2-1. 替换B原子(联合取代、完全枚举)
        3. 2-2 先分类再替代
        4. 2-3 每个结构随机取样特定数目的结构
      4. 4.枚举+迭代器进行工作
      5. 5.枚举-使用静电能作为哈希值而不是fit
      6. 6. 可以分步骤来做，先替换四个，再替换两个 ！
    5. 5.检验
      1. 1. 逐个检验相似性
      2. 2. 检验两个结构
  7. 7.找出缺失的原子
  8. 7. 处理CaH6的333倍胞的问题（将324个H中的14个扣掉-可以枚举的方法

1.pymatgen检验相同结构

fcc、bcc等对应着特殊空间群里特定位置的占据

1.一次性检验

import os 
from pymatgen.core import Structure
from pymatgen.io.vasp.sets import batch_write_input, MPRelaxSet
from pymatgen.io.vasp.inputs import Poscar
standard_structures= []
dir='./'
for i in range(1000):
    filename = "{}.vasp".format(i)
    if not os.path.exists(os.path.join(dir,filename)):
        continue
    print(filename)
    structure = Structure.from_file(filename)
    structure.remove_species("O")
    standard_structures.append(structure)
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
groups = matcher.group_structures([d for d in standard_structures])
print(len(groups))
groupstructures = [x[0] for x in groups]  # 将嵌套列表的每个列表的第一个元素提取出来做成一个新列表
from pymatgen.io.vasp.inputs import Poscar
n = 0
for Cry_Str in groupstructures:
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%dnew.vasp'%n)

2.逐个检验

import numpy as np
import sys,os,time
from pymatgen.core import Structure
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
all_str = []
for i in range(1,100000):
    file = f'{i}.vasp'
    if not os.path.exists(file):
       continue
    print(file)
    strua = Structure.from_file(file)
    match_found = False
    for str in all_str:
        if matcher.fit(strua,str):
           match_found = True
           break
    if not match_found:
        all_str.append(strua)
        print('%d structures' % len(all_str))
from pymatgen.io.vasp.inputs import Poscar
n = 0
for strua in all_str:
    Vasp_Str = Poscar(strua)
    n +=1
    Vasp_Str.write_file('%dnew.vasp'%n)

3. 检验两个结构

import numpy as np
import sys,os,time
from pymatgen.core import Structure
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
aa = Structure.from_file('16.vasp')
bb = Structure.from_file('4199.vasp')
matcher.fit(aa,bb)

2. pymatgen-enumeration

from pymatgen.core import Structure
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.transformations.advanced_transformations import EnumerateStructureTransformation
from pymatgen.io.vasp.sets import batch_write_input, MPRelaxSet
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
from pymatgen.analysis.bond_valence import BVAnalyzer
structure = Structure.from_file("a.cif")
data = {"Na":1, "Y":3, "Si":4,"O":-2,"Al":3,"K":1}
analyzer = SpacegroupAnalyzer(structure)
prim_cell = analyzer.find_primitive()
prim_cell.add_oxidation_state_by_element(data)
#print(prim_cell)
enum = EnumerateStructureTransformation()
enumerated = enum.apply_transformation(prim_cell, 100)
structures = [d["structure"] for d in enumerated]
print("%d structures returned." % len(structures))
from pymatgen.io.vasp.inputs import Poscar
n = 0
for Cry_Str in structures:
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%d.vasp'%n)
#from pymatgen.io.cif import CifWriter

#n = 0
#for Cry_Str in structures:
#    n = n + 1
#    cif_filename = f"{n}.cif"
#    cif_writer = CifWriter(Cry_Str)
#    cif_writer.write_file(cif_filename)   

#energy= [d["energy"] for d in enumerated]
#fl = open('energy.dat','w')
#fl.writelines(str(energy))
#data = {"K":1, "O":-2, "Fe":2,"Ti":4.25,"Mn":3.75}
#abc=[]
#for i in range(1,101):
#    m = Structure.from_file('%d.vasp'%i)
#    m.add_oxidation_state_by_element(data)
#    energy = EwaldElectrostaticModel().get_energy(m)
#    abc.append(energy)
#f2 = open('energy2.dat','w')
#f2.writelines(str(abc))
with open('energy.dat', 'w') as f:
    for item in abc:
        f.write(str(item) + '\n')

#strcturelist2.sort(key = EwaldElectrostaticModel().get_energy)

3. pymatgen-orderdisorder

import os
from pymatgen.io.cif import CifParser, CifWriter
from pymatgen.transformations.standard_transformations import SubstitutionTransformation, OrderDisorderedStructureTransformation, SupercellTransformation
from pymatgen.core import Structure
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
from pymatgen.transformations.advanced_transformations import EnumerateStructureTransformation
from pymatgen.io.vasp.sets import batch_write_input, MPRelaxSet
from pymatgen.io.vasp.inputs import Poscar
structure = Structure.from_file("a.cif")
data = {"Ni":2, "Co":3, "S":-2,"Se":-2}
analyzer = SpacegroupAnalyzer(structure)
prim_cell = analyzer.find_primitive()
prim_cell.add_oxidation_state_by_element(data)
print(prim_cell)
order = OrderDisorderedStructureTransformation()
standard_structures=order.apply_transformation(prim_cell,return_ranked_list=1000)
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
groups = matcher.group_structures([d["structure"] for d in standard_structures]) #这应该是个字典，把字典中的structure关键字提出来组>成列表,groups里只剩下结构了

#structures = [d["structure"] for d in standard_structures]
#print("%d structures returned." % len(structures))

#from pymatgen.io.vasp.inputs import Poscar
#n = 0
#for Cry_Str in structures:
    #print(Cry_Str)
#    Vasp_Str = Poscar(Cry_Str)
#    n = n+1
#    Vasp_Str.write_file('%d.vasp'%n)
groupstructures = [x[0] for x in groups]  # 将嵌套列表的每个列表的第一个元素提取出来做成一个新列表
from pymatgen.io.vasp.inputs import Poscar
n = 0
for Cry_Str in groupstructures:
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%d.vasp'%n)
#产生等差数列 for i in range(len(groups))
#提取出列举时计算的能量 energy = [x('energy') for x in standard_structures]
#验证同一类群下的结构静电能相等
#suoyoua = [x for x in groups[0]]
#for cry in suoyoua:
#   energy = EwaldElectrostaticModel().get_energy(cry)
#  print(energy)
#验证输出的顺序是按照能量升序进行的
#suoyou1 = [x[0] for x in groups]
#for cry in suoyou1:
#    energy = EwaldElectrostaticModel().get_energy(cry)
#    print(energy)

4.计算静电能

import os
from pymatgen.core import Structure
from pymatgen.io.vasp.inputs import Poscar
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
abc=[]
data = {"Ca":-6, "H":1}
for i in range(1,5034):
    m = Structure.from_file('%d.vasp'%i)
    m.add_oxidation_state_by_element(data)
    energy = EwaldElectrostaticModel().get_energy(m)
    abc.append(energy)
with open('energy.dat', 'w') as f:
    for item in abc:
        f.write(str(item) + '\n')

5. sqs

1.前处理

1. 转换rndstr格式的结构

from pymatgen.core import Structure
structure = Structure.from_file("a.cif")
structure.to(filename='rndstr.in')

2.提交多个ip的任务

for i in $(seq 1 1 20)
do
cat > sqs_$i.slurm << liz
#!/bin/sh
#SBATCH -J analyze
#SBATCH -p xieyu
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -w c0406


source activate liz

mcsqs -n=36 -ip=$i
liz
sbatch sqs_$i.slurm
done
#grep Objective_function bestcorr*.out|sort -n

2. 后处理

1. 提取bestcorr

for i in $(seq 1 1 12)
do
tail -n 1 bestcorr.out >> corr.dat
done

2.对多任务ip产生的best.out文件处理

==注意，数字前的空格会对结构在vesta中的显示产生影响==

sqs2poscar文件在 同步空间/计算脚本/sqs处

for i in $(seq 1 1 12)
do 
./sqs2poscar bestsqs$i.out
mv bestsqs$i.out-POSCAR $i.vasp
sed -i '2s/.*/1/' $i.vasp
sed -i '3s/-10\.57524/ -10.57524/' $i.vasp
done

3. 检验相同的结构

import os
from pymatgen.core import Structure
from pymatgen.io.vasp.sets import batch_write_input, MPRelaxSet
from pymatgen.io.vasp.inputs import Poscar
standard_structures= []
for i in range(1,13):
    filename = "{}.vasp".format(i)
    print(filename)
    structure = Structure.from_file(filename)
    standard_structures.append(structure)
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
groups = matcher.group_structures([d for d in standard_structures])  #这应该是个字典，把字典中的structure关键字提出来组成列表,groups里只剩下结构了

#structures = [d["structure"] for d in standard_structures]
#print("%d structures returned." % len(structures))

#from pymatgen.io.vasp.inputs import Poscar
#n = 0
#for Cry_Str in structures:
    #print(Cry_Str)
#    Vasp_Str = Poscar(Cry_Str)
#    n = n+1
#    Vasp_Str.write_file('%d.vasp'%n)
print(len(groups))
groupstructures = [x[0] for x in groups]  # 将嵌套列表的每个列表的第一个元素提取出来做成一个新列表
from pymatgen.io.vasp.inputs import Poscar
n = 0
for Cry_Str in groupstructures:
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%dnew.vasp'%n)
#产生等差数列 for i in range(len(groups))
#提取出列举时计算的能量 energy = [x('energy') for x in standard_structures]
#验证同一类群下的结构静电能相等
#suoyoua = [x for x in groups[0]]
#for cry in suoyoua:
#   energy = EwaldElectrostaticModel().get_energy(cry)
#  print(energy)
#验证输出的顺序是按照能量升序进行的
#suoyou1 = [x[0] for x in groups]
#for cry in suoyou1:
#    energy = EwaldElectrostaticModel().get_energy(cry)
#    print(energy)

6.暴力枚举

1.枚举脚本

在替换多种元素时，可以只替换一种元素，然后在生成的poscar中直接修改元素种类和个数

#!/usr/bin/python2
#Change the format of crystal structure
#python Str_Format.py infile outfile
try:
    from pymatgen.core.structure          import IStructure, Structure
    from pymatgen.io.vasp.inputs          import Poscar
except:
    print ('You should install pymatgen first.')
    exit(0)

from monty.io import zopen
import random
import sys,os,time
import numpy as np

def readvasp(infile):
    with zopen(infile, "rt") as f:
        contents = f.read()

    Cry_Str = IStructure.from_str(contents, fmt="poscar")
    lattice = Cry_Str.lattice
    species = [ x.species.elements[0] for x in Cry_Str.sites ]
    coords  = [ x.coords for x in Cry_Str.sites ]
    Cry_Str = Structure(lattice,species,coords,coords_are_cartesian=True)
    return Cry_Str

def replace(Cry_Str, replace_list, outfile):
    natom = len(Cry_Str.sites)
    for ele in replace_list:
        for i in ele[0]:
            Cry_Str.replace(i,ele[1])

    Cry_Str.sort()
    Vasp_Str = Poscar(Cry_Str)
    Vasp_Str.write_file(outfile)
    return


if __name__ == '__main__':
    instr = 'a.vasp'

    Cry_Str = readvasp(instr)
    #modify
    Cry_Str.make_supercell([1,1,1])
    for i in range(10):

        Cry_tp = Cry_Str.copy()
        re_index = np.arange(28,82)
        np.random.shuffle(re_index)
        re_index1 = re_index[:3]
        re_index2 = re_index[7:9]



        replace_list = [ [re_index1, 'B'],[re_index2,'O'] ]

        tmp = '%d.vasp'%i
        replace(Cry_tp, replace_list, tmp)

2. 多任务脚本 检验相同结构

for i in $(seq 0 1000 25000)
do 
cat > job_$i.sh << liz
#!/bin/sh
#SBATCH -J liz
#SBATCH -p xieyu
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -w c0607
source activate liz

python jianyan_$i.py

liz
cat > jianyan_$i.py << liz
import os
from pymatgen.core import Structure
from pymatgen.io.vasp.sets import batch_write_input, MPRelaxSet
from pymatgen.io.vasp.inputs import Poscar
standard_structures= []
dir='/work/home/liz/workspace/D-LuNH-sqs/2-Lu-25电子/1-暴力破解/3-zuhe'
for i in range($i,$((i+1000))):
    filename = "{}.vasp".format(i)
    if not os.path.exists(os.path.join(dir,filename)):
        continue
    print(filename)
    structure = Structure.from_file(filename)
    standard_structures.append(structure)
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
groups = matcher.group_structures([d for d in standard_structures])
print(len(groups))
groupstructures = [x[0] for x in groups]  # 将嵌套列表的每个列表的第一个元素提取出来做成一个新列表
from pymatgen.io.vasp.inputs import Poscar
n = $i
for Cry_Str in groupstructures:
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%dnew.vasp'%n)
liz
sbatch job_$i.sh
done

3.改良版

1.限定最小距离

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure

strua = Structure.from_file("a.vasp")
#print(strua.sites[re_index1[1]])
#print(strua.get_distance(re_index1[0],re_index1[1]))

#a_frac_coords=strua.sites[re_index1[0]].frac_coords
#b_frac_coords=strua.sites[re_index1[1]].frac_coords
#c_frac_coords=strua.sites[re_index1[2]].frac_coords
#arr1=a_frac_coords-b_frac_coords
#arr1[arr1 >= 0.5] = 1
#arr1[arr1 <= -0.5] = -1
#arr1[np.logical_and(arr1 > -0.5, arr1 < 0.5)] = 0
#print(arr1)
#strua.sites[re_index1[1]].frac_coords +=arr1
#print(strua.sites[re_index1[1]])
#print(strua.get_distance(re_index1[0],re_index1[1]))
#diff1=np.sqrt(np.dot(strua.sites[re_index1[0]].coords-strua.sites[re_index1[1]].coords,strua.sites[re_index1[0]].coords-strua.sites[re_index1[1]].coords))
#if diff1<3
#print(strua.sites[re_index1[1]].frac_coords)
#h2=a_frac_coords-c_frac_coords
#h3=b_frac_coords-c_frac_coords
#lattice = strua.lattice
#lattice_matrix =np.array(lattice.matrix)
#coords  = [ x.coords for x in strua.sites ]
#矩阵乘法 np.dot(A,B)
#矩阵改变维数 A.reshape(1,3)
#sites从0开始
#np.dot(a,a)求平方和
#strua.sites[0].coords 输出笛卡尔坐标 .frac_coords分数坐标
#print(np.sqrt(np.dot(strua.sites[70].coords-strua.sites[71].coords,strua.sites[70].coords-strua.sites[71].coords)))
#求两个原子的距离
for i in range(20):
    Cry_tp = strua.copy()
    re_index = np.arange(27,82)
    np.random.shuffle(re_index)
    re_index1 = re_index[:3]
    replace_list = [ [re_index1, 'N'] ]
    if Cry_tp.get_distance(re_index1[0],re_index1[1])<3.6 or Cry_tp.get_distance(re_index1[0],re_index1[2])<3.6 or Cry_tp.get_distance(re_index1[1],re_index1[2])<3.6:
        continue
    for ele in replace_list:
        for j in ele[0]:
            Cry_tp.replace(j,ele[1])
    Cry_tp.sort()
    Vasp_Str = Poscar(Cry_tp)
    Vasp_Str.write_file('%d.vasp'%i)

2.限定最小距离且保证替代发生在不同的位置

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure

strua = Structure.from_file("a.vasp")
#print(strua.sites[re_index1[1]])
#print(strua.get_distance(re_index1[0],re_index1[1]))

#a_frac_coords=strua.sites[re_index1[0]].frac_coords
#b_frac_coords=strua.sites[re_index1[1]].frac_coords
#c_frac_coords=strua.sites[re_index1[2]].frac_coords
#arr1=a_frac_coords-b_frac_coords
#arr1[arr1 >= 0.5] = 1
#arr1[arr1 <= -0.5] = -1
#arr1[np.logical_and(arr1 > -0.5, arr1 < 0.5)] = 0
#print(arr1)
#strua.sites[re_index1[1]].frac_coords +=arr1
#print(strua.sites[re_index1[1]])
#print(strua.get_distance(re_index1[0],re_index1[1]))
#diff1=np.sqrt(np.dot(strua.sites[re_index1[0]].coords-strua.sites[re_index1[1]].coords,strua.sites[re_index1[0]].coords-strua.sites[re_index1[1]].coords))
#if diff1<3
#print(strua.sites[re_index1[1]].frac_coords)
#h2=a_frac_coords-c_frac_coords
#h3=b_frac_coords-c_frac_coords
#lattice = strua.lattice
#lattice_matrix =np.array(lattice.matrix)
#coords  = [ x.coords for x in strua.sites ]
#矩阵乘法 np.dot(A,B)
#矩阵改变维数 A.reshape(1,3)
#sites从0开始
#np.dot(a,a)求平方和
#strua.sites[0].coords 输出笛卡尔坐标 .frac_coords分数坐标
#print(np.sqrt(np.dot(strua.sites[70].coords-strua.sites[71].coords,strua.sites[70].coords-strua.sites[71].coords)))
#求两个原子的距离
for i in range(20):
    Cry_tp = strua.copy()
    re_index = np.arange(27,108)
    np.random.shuffle(re_index)
    re_index1 = re_index[:3]
    if all(27<= num <=53 for num in re_index1):
        continue
    if all(54<= num <=107 for num in re_index1):
        continue
    replace_list = [ [re_index1, 'N'] ]
    if Cry_tp.get_distance(re_index1[0],re_index1[1])<3.6 or Cry_tp.get_distance(re_index1[0],re_index1[2])<3.6 or Cry_tp.get_distance(re_index1[1],re_index1[2])<3.6:
        continue
    for ele in replace_list:
        for j in ele[0]:
            Cry_tp.replace(j,ele[1])
    Cry_tp.sort()
    Vasp_Str = Poscar(Cry_tp)
    Vasp_Str.write_file('%d.vasp'%i)

3.限定空间群

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.analysis.structure_matcher import StructureMatcher
strua = Structure.from_file("a.vasp")
#print(strua.sites[re_index1[1]])
#print(strua.get_distance(re_index1[0],re_index1[1]))

#a_frac_coords=strua.sites[re_index1[0]].frac_coords
#b_frac_coords=strua.sites[re_index1[1]].frac_coords
#c_frac_coords=strua.sites[re_index1[2]].frac_coords
#arr1=a_frac_coords-b_frac_coords
#arr1[arr1 >= 0.5] = 1
#arr1[arr1 <= -0.5] = -1
#arr1[np.logical_and(arr1 > -0.5, arr1 < 0.5)] = 0
#print(arr1)
#strua.sites[re_index1[1]].frac_coords +=arr1
#print(strua.sites[re_index1[1]])
#print(strua.get_distance(re_index1[0],re_index1[1]))
#diff1=np.sqrt(np.dot(strua.sites[re_index1[0]].coords-strua.sites[re_index1[1]].coords,strua.sites[re_index1[0]].coords-strua.sites[re_index1[1]].coords))
#if diff1<3
#print(strua.sites[re_index1[1]].frac_coords)
#h2=a_frac_coords-c_frac_coords
#h3=b_frac_coords-c_frac_coords
#lattice = strua.lattice
#lattice_matrix =np.array(lattice.matrix)
#coords  = [ x.coords for x in strua.sites ]
#矩阵乘法 np.dot(A,B)
#矩阵改变维数 A.reshape(1,3)
#sites从0开始
#np.dot(a,a)求平方和
#strua.sites[0].coords 输出笛卡尔坐标 .frac_coords分数坐标
#print(np.sqrt(np.dot(strua.sites[70].coords-strua.sites[71].coords,strua.sites[70].coords-strua.sites[71].coords)))
#求两个原子的距离
total_structure=[]
for i in range(200000):
    Cry_tp = strua.copy()
    re_index = np.arange(0,96)
    np.random.shuffle(re_index)
    re_index1 = re_index[:4]
    Cry_tp.remove_sites(re_index1)
    a=SpacegroupAnalyzer(Cry_tp).get_space_group_number()
    if a<= 140:
      continue
    print(i)
    print(a)
    Vasp_Str = Poscar(Cry_tp)
    Vasp_Str.write_file('%d.vasp'%i)
    total_structure.append(Cry_tp)
matcher = StructureMatcher()
groups = matcher.group_structures([d for d in total_structure])
groupstructures = [x[0] for x in groups]
from pymatgen.io.vasp.inputs import Poscar
n = 10001
for Cry_Str in groupstructures:
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%dnew.vasp'%n)

4.多任务做改良版

for i in $(seq 1 100 10100)
do
cat > baoli_$i.py << liz
from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.analysis.structure_matcher import StructureMatcher
strua = Structure.from_file("a.vasp")
#print(strua.sites[re_index1[1]])
#print(strua.get_distance(re_index1[0],re_index1[1]))

#a_frac_coords=strua.sites[re_index1[0]].frac_coords
#b_frac_coords=strua.sites[re_index1[1]].frac_coords
#c_frac_coords=strua.sites[re_index1[2]].frac_coords
#arr1=a_frac_coords-b_frac_coords
#arr1[arr1 >= 0.5] = 1
#arr1[arr1 <= -0.5] = -1
#arr1[np.logical_and(arr1 > -0.5, arr1 < 0.5)] = 0
#print(arr1)
#strua.sites[re_index1[1]].frac_coords +=arr1
#print(strua.sites[re_index1[1]])
#print(strua.get_distance(re_index1[0],re_index1[1]))
#diff1=np.sqrt(np.dot(strua.sites[re_index1[0]].coords-strua.sites[re_index1[1]].coords,strua.sites[re_index1[0]].coords-strua.sites[re_index1[1]].coords))
#if diff1<3
#print(strua.sites[re_index1[1]].frac_coords)
#h2=a_frac_coords-c_frac_coords
#h3=b_frac_coords-c_frac_coords
#lattice = strua.lattice
#lattice_matrix =np.array(lattice.matrix)
#coords  = [ x.coords for x in strua.sites ]
#矩阵乘法 np.dot(A,B)
#矩阵改变维数 A.reshape(1,3)
#sites从0开始
#np.dot(a,a)求平方和
#strua.sites[0].coords 输出笛卡尔坐标 .frac_coords分数坐标
#print(np.sqrt(np.dot(strua.sites[70].coords-strua.sites[71].coords,strua.sites[70].coords-strua.sites[71].coords)))
#求两个原子的距离
total_structure=[]
for i in range(20000):
    Cry_tp = strua.copy()
    re_index = np.arange(0,96)
    np.random.shuffle(re_index)
    re_index1 = re_index[:4]
    Cry_tp.remove_sites(re_index2)
    a=SpacegroupAnalyzer(Cry_tp).get_space_group_number()
    if a<= 140:
      continue
    print(i)
    print(a)
    Vasp_Str = Poscar(Cry_tp)
    Vasp_Str.write_file('%d.vasp'%i)
    total_structure.append(Cry_tp)
matcher = StructureMatcher()
groups = matcher.group_structures([d for d in total_structure])
groupstructures = [x[0] for x in groups]
from pymatgen.io.vasp.inputs import Poscar
n = $((i))
for Cry_Str in groupstructures:
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%dnew.vasp'%n)
liz
done

4.使用枚举方法处理蓝钻结构

  1. 这是第一次系统使用暴力枚举法产生结构，其中，由于pymatgen用于计算ewald加和的程序太慢，使用了组内ARES中用于计算ewald能量的脚本ewalds.x
  2. 一共有两种方式产生，一种是完全随机的，另一种相当于对整个分子的替代（使用了pymatgen中找原子附近其他原子的功能
  3. 脚本、流程都经过了反复的改进

1.完全随机的方法（完全解放双手了）

  1. 一共需要四个文件

a.vasp

multitask.sh

ewalds.x

changeatom.py:用于产生随机的结构

以及一个后处理文件

get-all.sh

2.分为三级别目录

一级： 1-21

二级： 1-48

三级： 1-random 2-low ；其中 1-random是随机产生的结构，并计算了静电能，2-low是能量最低的前30个结构以及对应的能量

  1. 通过设置不同的价态，可以得到不同种原子靠近或者远离的结构

流程(产生1000万个随机结构）：

1. 提任务

for i in $(seq 1 1 10);do mkdir $i;cp a.vasp ewalds.x multitask.sh changeatom.py $i;cd $i; ./multitask.sh ;cd ..;done

2. changeatom.py

需要更改的地方是产生结构的数量、以及替换的种类和位置

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
strua = Structure.from_file("a.vasp")
#需要考虑
for i in range(10000):
    print(i)
    Cry_tp = strua.copy()
    re_index = np.arange(0,512)
    np.random.shuffle(re_index)
    re_index1 = re_index[:7]
    re_index2 = re_index[7:9]
    replace_list = [ [re_index1, 'B'],[re_index2,'O'] ]
    for ele in replace_list:
        for j in ele[0]:
            Cry_tp.replace(j,ele[1])
    Cry_tp.sort()
    Vasp_Str = Poscar(Cry_tp)
    Vasp_Str.write_file('%d.vasp'%i)

3. multitask.sh :做一些文件夹 ，提交产生结构的任务，提交计算ewald能量的任务

需要注意的是，使用ewalds.x需要先source环境

for i in $(seq 1 1 48)
do
mkdir $i
cd $i
mkdir 1-random 2-low
cd ..
cp changeatom.py $i/1-random
cp ewalds.x $i/1-random
cp a.vasp $i/1-random
cd $i/1-random
cat > job.slurm << liz
#!/bin/sh
#SBATCH  --job-name=opt 
#SBATCH  --output=log.out.%j
#SBATCH  --error=log.err.%j
#SBATCH  --partition=xieyu
#SBATCH  --nodes=1
#SBATCH  --ntasks=1
#SBATCH  --cpus-per-task=1
source activate liz

python changeatom.py > log 2>&1

liz

sbatch job.slurm
cd ..
cd ..
done



sleep 1m
for i in $(seq 1 1 48)
do
cd $i/1-random
cat > ewaldenergy.slurm <<liz2
#!/bin/sh
#SBATCH  --job-name=opt 
#SBATCH  --output=log.out.%j
#SBATCH  --error=log.err.%j
#SBATCH  --partition=xieyu
#SBATCH  --nodes=1
#SBATCH  --ntasks=1
#SBATCH  --cpus-per-task=1
source /work/env/intel2020

for j in \$(seq 0 1 9999);
do
echo -n \$j.vasp >> energy.dat
./ewalds.x \$j.vasp 3 5 -2 >> energy.dat
done
sort -k3n energy.dat > ewaldenergy.dat
head -n 30 ewaldenergy.dat > ../2-low/ewaldenergy.dat
head -n 30 ewaldenergy.dat | awk '{print \$1}' | while read -r line; do
cp "\$line" ../2-low/
done
liz2
sbatch ewaldenergy.slurm
cd ..
cd ..
done

4. 后处理文件 get-all.sh 提取能量最低的结构

rm all-sorted.dat
for i in $(seq 1 1 21)
do
cd $i
for j in $(seq 1 1 48)
do
cd $j
cd 2-low
#awk -v var="$j" '{ $1 = var " " $1 } 1' ewaldenergy.dat > data.dat
awk -v var="$i" -v var2="$j" '{ $1 = var " "  var2 " " $1 } 1' ewaldenergy.dat > data.dat
cat data.dat >> ../../data.dat
rm data.dat
cd ..
cd ..
done
cd ..
done
for h in $(seq 1 1 21)
do
cat $h/data.dat >> ./all.dat
rm $h/data.dat
done
sort -k5n all.dat  > all-sorted.dat
rm all.dat

echo 'sucess1'
# 设置计数器
count=1
rm -r all
mkdir all
# 读取包含三列信息的文件
head -n 100 all-sorted.dat | while read -r line; do
  # 提取大目录、小目录和文件名
  directory=$(echo "$line" | awk '{print $1}')
  subdirectory=$(echo "$line" | awk '{print $2}')
  filename=$(echo "$line" | awk '{print $3}')

  # 构建目标文件名
  target_filename="${count}_${directory}_${subdirectory}_${filename}.vasp"

  # 复制文件到当前目录并重命名
  cp "${directory}/${subdirectory}/2-low/${filename}" "./all/${target_filename}"

  # 增加计数器
  ((count++))
done
echo 'sucess2'

2. 替换为B_3O和B_4O构型的方法（还需要手动做一些流程）

  1. 替换为特定构型的方法利用了pymatgen中寻找相邻原子的方法，基本的流程是

先随机找两个可能替换氧原子的位置

检验氧原子的距离是否符合要求（在本例中设定两个氧原子不要靠的太近，以免产生O-B-O键

查找两个位置附近的原子的位置

替换附近原子为B

检验相同结构

计算静电能

  1. 另一套流程是：(下文的流程)

先随机替换氧的位置

检验氧原子距离

检验氧的相同结构

替换B

检验相同结构

计算静电能

这样流程的好处是可以大大减少产生结构的数量，缺点是不能把算能量的步骤集成

  1. 值得注意的一点是，在这次的产生结构的脚本中，是可以枚举O的所有可能情况的，所以最终结果是枚举的

image-20230612211658194

1. 枚举氧的位置

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
from itertools import combinations
strua = Structure.from_file("a.vasp")
hh=np.arange(512)
combinations_list1 = list(combinations(hh, 2))
print(len(combinations_list1))
combinations_array1 = np.array(combinations_list1)
numberbb = 0
numberaa = 0
for ar in combinations_array1:
    ara=strua.copy()
    numberbb +=1
    print(numberbb)
    if ara.get_distance(ar[0],ar[1])>11 or ara.get_distance(ar[0],ar[1])<3:
        continue
    numberaa +=1
    replace_list = [[ar,'O'] ]
    for ele in replace_list:
        for j in ele[0]:
            ara.replace(j,ele[1])
    ara.sort()
    Vasp_Str = Poscar(ara)
    Vasp_Str.write_file('%d.vasp'% numberaa)

2. 检验相同的结构

需要多次分批次检验，否则太慢了，有一个多任务脚本和一个py脚本

在不同文件夹中多次进行多脚本计算，会比较方便

image-20230612172111471

cp ../*new.vasp ./1-first

第一次检验

for i in $(seq 0 1000 130000)
do
cat > job_$i.sh << liz
#!/bin/sh
#SBATCH  --job-name=opt 
#SBATCH  --output=log.out.%j
#SBATCH  --error=log.err.%j
#SBATCH  --partition=xieyu
#SBATCH  --nodes=1
#SBATCH  --ntasks=1
#SBATCH  --cpus-per-task=1
source activate liz

python jianyan_$i.py >> jianyan_$i.dat 2>&1

liz
cat > jianyan_$i.py << liz
import os
from pymatgen.core import Structure
from pymatgen.io.vasp.sets import batch_write_input, MPRelaxSet
from pymatgen.io.vasp.inputs import Poscar
standard_structures= []
dir='/work/home/liz/workspace/3-temporary-duty/2-blue/4-zhengshi/3-B3O-B4O/4-xinzuo/1-替换O/1-str'
for j in range($i,$((i+1000))):
    filename = "{}.vasp".format(j)
    filepath = os.path.join(dir, filename)
    if not os.path.exists(filepath):
        continue
    print(filename)
    structure = Structure.from_file(filepath)
    standard_structures.append(structure)
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
groups = matcher.group_structures([d for d in standard_structures])
print(len(groups))
groupstructures = [x[0] for x in groups]  # 将嵌套列表的每个列表的第一个元素提取出来做成一个新列表
from pymatgen.io.vasp.inputs import Poscar
n = $i
for Cry_Str in groupstructures:
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%dnew.vasp'%n)
liz
sbatch job_$i.sh
done

第二次检验

for i in $(seq 0 10000 130000)
do
cat > job_$i.sh << liz
#!/bin/sh
#SBATCH  --job-name=opt 
#SBATCH  --output=log.out.%j
#SBATCH  --error=log.err.%j
#SBATCH  --partition=xieyu
#SBATCH  --nodes=1
#SBATCH  --ntasks=1
#SBATCH  --cpus-per-task=1
source activate liz

python jianyan_$i.py >> jianyan_$i.dat 2>&1

liz
cat > jianyan_$i.py << liz
import os
from pymatgen.core import Structure
from pymatgen.io.vasp.sets import batch_write_input, MPRelaxSet
from pymatgen.io.vasp.inputs import Poscar
standard_structures= []
dir='/work/home/liz/workspace/3-temporary-duty/2-blue/4-zhengshi/3-B3O-B4O/4-xinzuo/2-去除相同结构/1-first'
for j in range($i,$((i+11000))):
    filename = "{}new.vasp".format(j)
    filepath = os.path.join(dir, filename)
    if not os.path.exists(filepath):
        continue
    print(filename)
    structure = Structure.from_file(filepath)
    standard_structures.append(structure)
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
groups = matcher.group_structures([d for d in standard_structures])
print(len(groups))
groupstructures = [x[0] for x in groups]  # 将嵌套列表的每个列表的第一个元素提取出来做成一个新列表
from pymatgen.io.vasp.inputs import Poscar
n = $i
for Cry_Str in groupstructures:
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%dnew.vasp'%n)
liz
sbatch job_$i.sh
done

最后一步的检验

import os
from pymatgen.core import Structure
from pymatgen.io.vasp.sets import batch_write_input, MPRelaxSet
from pymatgen.io.vasp.inputs import Poscar
standard_structures= []
dir='/work/home/liz/workspace/3-temporary-duty/2-blue/4-zhengshi/3-B3O-B4O/4-xinzuo/2-去除相同结构/2-new'
for j in range(124000):
    filename = "{}new.vasp".format(j)
    filepath = os.path.join(dir, filename)
    if not os.path.exists(filepath):
        continue
    print(filename)
    structure = Structure.from_file(filepath)
    standard_structures.append(structure)
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
groups = matcher.group_structures([d for d in standard_structures])
print(len(groups))
groupstructures = [x[0] for x in groups]  # 将嵌套列表的每个列表的第一个元素提取出来做成一个新列表
from pymatgen.io.vasp.inputs import Poscar
n = 0
for Cry_Str in groupstructures:
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%dnew.vasp'%n)

3. 替换B原子

先替换B4，再替换B3

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
from itertools import combinations
dir='/work/home/liz/workspace/3-temporary-duty/2-blue/4-zhengshi/3-B3O-B4O/4-xinzuo/2-去除相同结构/2-new'
counter = 0
for j in range(124000):
    filename = "{}new.vasp".format(j)
    filepath = os.path.join(dir, filename)
    if not os.path.exists(filepath):
        continue
    print(filename)
    structure = Structure.from_file(filepath)
    o_indices = [i for i, site in enumerate(structure) if site.specie.name == "O"]
    o_array = np.array(o_indices)
    ara=structure.copy()
    list1=ara.get_neighbors(ara[o_array[0]], r=2.0)
    list2=ara.get_neighbors(ara[o_array[1]], r=2.0)
    list1index=np.array([site.index for site in list1])
    list2index=np.array([site.index for site in list2])

    replace_list = [ [list1index, 'B'] ]
    for ele in replace_list:
        for j in ele[0]:
            ara.replace(j,ele[1])
    combinations_list = list(combinations(list2index, 3))
    combinations_array = np.array(combinations_list)
    for br in combinations_array:
        counter += 1
        brb=ara.copy()
        replace_listara = [ [br, 'B'] ]
        for ele in replace_listara:
           for j in ele[0]:
                 brb.replace(j,ele[1])
        brb.sort()
        Vasp_Str = Poscar(brb)
        Vasp_Str.write_file('%d.vasp'% counter)
        print(counter)

4. 再次检验相同结构

import os
from pymatgen.core import Structure
from pymatgen.io.vasp.sets import batch_write_input, MPRelaxSet
from pymatgen.io.vasp.inputs import Poscar
standard_structures= []
dir='/work/home/liz/workspace/3-temporary-duty/2-blue/4-zhengshi/3-B3O-B4O/4-xinzuo/3-替换B'
for j in range(500):
    filename = "{}.vasp".format(j)
    filepath = os.path.join(dir, filename)
    if not os.path.exists(filepath):
        continue
    print(filename)
    structure = Structure.from_file(filepath)
    standard_structures.append(structure)
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
groups = matcher.group_structures([d for d in standard_structures])
print(len(groups))
groupstructures = [x[0] for x in groups]  # 将嵌套列表的每个列表的第一个元素提取出来做成一个新列表
from pymatgen.io.vasp.inputs import Poscar
n = 0
for Cry_Str in groupstructures:
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%dnew.vasp'%n)

5. 计算能量，整理结构

rm energy.dat
for i in $(seq 1 1 99);
do 
echo -n ${i}new.vasp >> energy.dat
./ewalds.x ../4-去除相同结构/${i}new.vasp -1 0 2 >> energy.dat
done
sort -k3n energy.dat  > sorted.dat
mkdir str
cd str
count=1
head -n 100 ../sorted.dat |while read -r line;do
filename=$(echo "$line" | awk '{print $1}')
target_filename=${count}.vasp
cp ../../4-去除相同结构/${filename} ./${target_filename}
count=$((count+ 1))
done
cd ..

3. 替换为B_2O和B_3O的方法

1.枚举氧的结构和其他步骤

2-1. 替换B原子(联合取代、完全枚举)

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
from itertools import combinations
dir='/work/home/liz/workspace/3-temporary-duty/2-blue/4-zhengshi/1-512-B-7-O-2/3-B2O-B3O/1-O-结构'
counter = 0
filename = "1new.vasp"
filepath = os.path.join(dir, filename)
print(filename)
structure = Structure.from_file(filepath)
o_indices = [i for i, site in enumerate(structure) if site.specie.name == "O"]
o_array = np.array(o_indices)
c_indices = [i for i, site in enumerate(structure) if site.specie.name == "C"]
c_array = np.array(c_indices)
ara=structure.copy()
list1=ara.get_neighbors(ara[o_array[0]], r=2.0)
list2=ara.get_neighbors(ara[o_array[1]], r=2.0)
list1index=np.array([site.index for site in list1])
list2index=np.array([site.index for site in list2])
combined_array = np.concatenate((list1index, list2index))
c_array2 = np.setdiff1d(c_array, combined_array)
combinations_list1 = list(combinations(list1index, 2))
combinations_list2 = list(combinations(list2index, 3))
combinations_list3 = list(combinations(c_array2, 2))
array1 = np.array(combinations_list1)
array2 = np.array(combinations_list2)
array3 = np.array(combinations_list3)
print(array3.shape)
print(len(array1))
print(len(array2))
print(len(array3))
for i in range(len(array1)):
    for j in range(len(array2)):
        for k in range(len(array3)):
            brb = ara.copy()
            counter += 1
           # print(array1[i])
           # print(array2[j])
           # print(array3[k])
            combined_array33=np.concatenate((array1[i], array2[j], array3[k]))
            replace_listara = [ [combined_array33, 'B'] ]
            for ele in replace_listara:
                for m in ele[0]:
                   brb.replace(m,ele[1])
            brb.sort()
            Vasp_Str = Poscar(brb)
            Vasp_Str.write_file('%d.vasp'% counter)
            print(counter)

2-2 先分类再替代

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
from itertools import combinations
dir='/work/home/liz/workspace/3-temporary-duty/2-blue/4-zhengshi/1-512-B-7-O-2/3-B2O-B3O/1-O-结构'
counter = 0
filename = "1new.vasp"
filepath = os.path.join(dir, filename)
print(filename)
structure = Structure.from_file(filepath)
o_indices = [i for i, site in enumerate(structure) if site.specie.name == "O"]
o_array = np.array(o_indices)
c_indices = [i for i, site in enumerate(structure) if site.specie.name == "C"]
c_array = np.array(c_indices)
ara=structure.copy()
list1=ara.get_neighbors(ara[o_array[0]], r=2.0)
list2=ara.get_neighbors(ara[o_array[1]], r=2.0)
list1index=np.array([site.index for site in list1])
list2index=np.array([site.index for site in list2])
combined_array = np.concatenate((list1index, list2index))
print(combined_array)
c_array2 = np.setdiff1d(c_array, combined_array)
combinations_list1 = list(combinations(list1index, 2))
combinations_list2 = list(combinations(list2index, 3))
combinations_list3 = list(combinations(c_array2, 2))
array1 = np.array(combinations_list1)
array2 = np.array(combinations_list2)
array3 = np.array(combinations_list3)
print(array3.shape)
print(len(array1))
print(len(array2))
print(len(array3))
standard_structures=[]
for i in range(len(array1)):
    for j in range(len(array2)):
#        for k in range(len(array3)):
            brb = ara.copy()
           # print(array1[i])
           # print(array2[j])
           # print(array3[k])
            combined_array33=np.concatenate((array1[i], array2[j]))
#            print(combined_array33)
            replace_listara = [ [combined_array33, 'B'] ]
            for ele in replace_listara:
                for m in ele[0]:
                   brb.replace(m,ele[1])
            standard_structures.append(brb)
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
groups = matcher.group_structures([d for d in standard_structures])
print(len(groups))
groupstructures = [x[0] for x in groups]
for straaa in groupstructures:
    for k in range(len(array3)):
        copy2=straaa.copy()
        counter += 1
        replace_listara = [ [array3[k], 'B'] ]
        for ele in replace_listara:
            for m in ele[0]:
                copy2.replace(m,ele[1])
        copy2.sort()
        Vasp_Str = Poscar(copy2)

        Vasp_Str.write_file('%d.vasp'% counter)
        print(counter)

2-3 每个结构随机取样特定数目的结构

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
from itertools import combinations
dir='/work/home/liz/workspace/3-temporary-duty/2-blue/4-zhengshi/1-512-B-7-O-2/3-B2O-B3O/1-O-结构'
counter = 0
for j in range(100):
    filename = "{}new.vasp".format(j)
    filepath = os.path.join(dir, filename)
    if not os.path.exists(filepath):
        continue
    print(filename)
    structure = Structure.from_file(filepath)
    o_indices = [i for i, site in enumerate(structure) if site.specie.name == "O"]
    o_array = np.array(o_indices)
    c_indices = [i for i, site in enumerate(structure) if site.specie.name == "C"]
    c_array = np.array(c_indices)
    ara=structure.copy()
    list1=ara.get_neighbors(ara[o_array[0]], r=2.0)
    list2=ara.get_neighbors(ara[o_array[1]], r=2.0)
    list1index=np.array([site.index for site in list1])
    list2index=np.array([site.index for site in list2])
    combined_array = np.concatenate((list1index, list2index))
    c_array2 = np.setdiff1d(c_array, combined_array)
    combinations_list1 = list(combinations(list1index, 2))
    combinations_list2 = list(combinations(list2index, 3))
    combinations_list3 = list(combinations(c_array2, 2))
    array1 = np.array(combinations_list1)
    array2 = np.array(combinations_list2)
    array3 = np.array(combinations_list3)
    arr1=np.arange(len(array1))
    arr2=np.arange(len(array2))
    arr3=np.arange(len(array3))
    for ww in range(5000):
        np.random.shuffle(arr1)
        np.random.shuffle(arr2)
        np.random.shuffle(arr3)

        brb = ara.copy()
        counter += 1
        combined_array33=np.concatenate((array1[arr1[1]], array2[arr2[1]], array3[arr3[1]]))
        replace_listara = [ [combined_array33, 'B'] ]
        for ele in replace_listara:
            for m in ele[0]:
               brb.replace(m,ele[1])
        brb.sort()
        Vasp_Str = Poscar(brb)
        Vasp_Str.write_file('%d.vasp'% counter)
        print(counter)

4.枚举+迭代器进行工作

可以是可以，但是现在太慢了，问题出在挨个去对比结构相似性，一个可行的方案是不断地细分到不同的结构特征集合，用结构特征做一个初筛

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
from itertools import combinations
import itertools
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
strua = Structure.from_file("a.vasp")
hh=np.arange(16,112)
combinations_generator = itertools.combinations(hh, 6)
numberbb = 0
numberaa = 0
all_structure=[]
for ar in combinations_generator:
    ara=strua.copy()
    numberbb +=1
    print(numberbb)
    remove_list = list(ar)
    ara.remove_sites(remove_list)
    ara.sort()
    match_found =False
    for aa in all_structure:
        if matcher.fit(ara,aa):
            match_found = True
            break
    if not match_found:
        all_structure.append(ara)
        print("%d structures " %len(all_structure))
        
bb=0
for aa in all_structure:    
    bb +=1
    Vasp_Str = Poscar(aa)
    Vasp_Str.write_file('%d.vasp'% bb)

5.枚举-使用静电能作为哈希值而不是fit

这一方法验证错误，因为两个不同的结构可能具有相同的静电能

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
from itertools import combinations
import itertools
from pymatgen.analysis.structure_matcher import StructureMatcher

strua = Structure.from_file("a.vasp")
hh = np.arange(16, 112)
combinations_generator = itertools.combinations(hh, 6)
data = {"Ca":6, "H":-1}
structure_dict = {}  # 用于存储已处理的结构
strua.add_oxidation_state_by_element(data)
print(strua)
aaa= 0
for ar in combinations_generator:
    aaa +=1
    print(aaa)
    ara = strua.copy()
    remove_list = list(ar)
    ara.remove_sites(remove_list)
    ara.sort()
    
    # 计算结构的哈希值作为字典的键
    ara_hash = EwaldElectrostaticModel().get_energy(ara)
    
    # 使用字典进行匹配检查
    if ara_hash not in structure_dict:
        structure_dict[ara_hash] = ara
        print("New structure added! Total structures: %d" % len(structure_dict))


        

for Cry_Str in structure_dict.values():
    #print(Cry_Str)
    Vasp_Str = Poscar(Cry_Str)
    n = n+1
    Vasp_Str.write_file('%d.vasp'%n)

6. 可以分步骤来做，先替换四个，再替换两个 ！

先替换1个，再替换一个，再替换一个，再，，，，这样一定是可以做的

从上往下试，直到能够使用程序做了就停止（95–92），如果从下往上试（92—95）会面临文件数太多的困难

5.检验

1. 逐个检验相似性

import numpy as np
import sys,os,time
from pymatgen.core import Structure
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
all_str = []
for i in range(1,100000):
    file = f'{i}.vasp'
    if not os.path.exists(file):
       continue
    print(file)
    strua = Structure.from_file(file)
    match_found = False
    for str in all_str:
        if matcher.fit(strua,str):
           match_found = True
           break
    if not match_found:
        all_str.append(strua)
        print('%d structures' % len(all_str))
from pymatgen.io.vasp.inputs import Poscar
n = 0
for strua in all_str:
    Vasp_Str = Poscar(strua)
    n +=1
    Vasp_Str.write_file('%dnew.vasp'%n)

2. 检验两个结构

import numpy as np
import sys,os,time
from pymatgen.core import Structure
from pymatgen.analysis.structure_matcher import StructureMatcher
matcher = StructureMatcher()
aa = Structure.from_file('16.vasp')
bb = Structure.from_file('4199.vasp')
matcher.fit(aa,bb)

7.找出缺失的原子

from pymatgen.core.structure import IStructure, Structure
from itertools import combinations
stru_b=Structure.from_file(filename='supercell.vasp')
frac_coords_b=stru_b.frac_coords.tolist()
normalized_coords_b = [[x % 1 for x in coord] for coord in frac_coords_b]
set2 = set(map(tuple, normalized_coords_b))
for i in range(1, 1000001):
    file = f'{i}.vasp'
    stru_a=Structure.from_file(file)
    # 检查文件是否存在
    if not os.path.exists(file):
        continue
    frac_coords_a=stru_a.frac_coords.tolist()
    normalized_coords_a = [[x % 1 for x in coord] for coord in frac_coords_a]
    set1 = set(map(tuple, normalized_coords_a))
    result= [list(i) for i in set2.difference(set1)]
    print(result)
    atom_indices=[]
    for coord in result:
        try:
            index = normalized_coords_b.index(coord)
            atom_indices.append(index)
        except ValueError:
          # 如果找不到匹配的坐标，可以添加一个标记，如-1，以表示未匹配的坐标
            atom_indices.append(-1)
    close_atom_pairs = []
    combinations_of_missing_atoms = combinations(atom_indices, 2)
    threshold_distance = 7
    for pair in combinations_of_missing_atoms:
        distance = stru_b.get_distance(pair[0], pair[1])
        if distance < threshold_distance:
            close_atom_pairs.append(pair)
    # 如果 close_atom_pairs 列表不为空，则表示存在距离小于阈值的原子对
    if close_atom_pairs:
        print("存在距离小于阈值的原子对：")
        print(len(close_atom_pairs))
    else:
        print("不存在距离小于阈值的原子对。")

7. 处理CaH6的333倍胞的问题（将324个H中的14个扣掉-可以枚举的方法

==/work/home/liz/workspace/3-temporary-duty/3-zwb/6-333/1-kou3/2/1-2/ee.py==

from pymatgen.core.structure import IStructure, Structure
from pymatgen.io.vasp.inputs import Poscar
from monty.io import zopen
import random
import sys,os,time
import numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.energy_models import EwaldElectrostaticModel
from itertools import combinations
strua = Structure.from_file("a.vasp")

H_indice = [i for i, site in enumerate(strua) if site.species_string == "H"]
H_indices=np.array(H_indice)

distanceall=strua.distance_matrix
distance=distanceall[H_indices][:,H_indices]
print(distance.shape)
print(distance)

N = distance.shape[0]
two_dim_matrix = np.zeros((N, N), dtype=int)
for i in range(N):
    for j in range(N):
        dist1 = distance[i][j]
        if dist1 <= 3.6:
            two_dim_matrix[i][j] = 0
        else:
            two_dim_matrix[i][j] = 1
allnumber = 0
matching_columns = np.where(two_dim_matrix[0, 1:] == 1)[0]
matching_columns = matching_columns  + 1
print(matching_columns)
for j in matching_columns:
    print(j)
    matching_columns2=np.where(two_dim_matrix[j, j+1:] == 1)[0]
    matching_columns2 = matching_columns2 + j + 1
    matching_columns2 = np.intersect1d(matching_columns2,matching_columns)
   # print(len(matching_columns2))
    for jj in matching_columns2:
        matching_columns3=np.where(two_dim_matrix[jj, jj+1:] == 1)[0]
        matching_columns3 = matching_columns3 + jj + 1
        matching_columns3 = np.intersect1d(matching_columns2,matching_columns3)
      #  print(len(matching_columns3))
        for jjj in matching_columns3:
            matching_columns4=np.where(two_dim_matrix[jjj, jjj+1:] == 1)[0]
            matching_columns4 = matching_columns4 + jjj + 1
            matching_columns4 = np.intersect1d(matching_columns3,matching_columns4)
           # print(len(matching_columns4))
            for jjjj in matching_columns4:
                matching_columnsddd=np.where(two_dim_matrix[jjjj, jjjj+1:] == 1)[0]
                matching_columnsddd = matching_columnsddd + jjjj + 1
                matching_columnsddd = np.intersect1d(matching_columns4,matching_columnsddd)
               # print(len(matching_columnsddd))
                for gg in matching_columnsddd:
                    matching_columns5=np.where(two_dim_matrix[gg, gg+1:] == 1)[0]
                    matching_columns5 = matching_columns5 + gg + 1
                    matching_columns5 = np.intersect1d(matching_columnsddd,matching_columns5)
                   # print(len(matching_columns5))
                    for ggg in matching_columns5:
                        matching_columns6=np.where(two_dim_matrix[ggg, ggg+1:] == 1)[0]
                        matching_columns6 = matching_columns6 + ggg + 1
                        matching_columns6 = np.intersect1d(matching_columns5,matching_columns6)
                       # print(len(matching_columns6))
                        for gggg in matching_columns6:
                            matching_columns7=np.where(two_dim_matrix[gggg, gggg+1:] == 1)[0]
                            matching_columns7 = matching_columns7 + gggg + 1
                            matching_columns7 = np.intersect1d(matching_columns6,matching_columns7)
                           # print(len(matching_columns7))
                            for dd in matching_columns7:
                                matching_columns8=np.where(two_dim_matrix[dd, dd+1:] == 1)[0]
                                matching_columns8 = matching_columns8 + dd + 1
                                matching_columns8 = np.intersect1d(matching_columns7,matching_columns8)
                               # print(len(matching_columns8))
                                for ddd in matching_columns8:
                                    matching_columns9=np.where(two_dim_matrix[ddd, ddd+1:] == 1)[0]
                                    matching_columns9 = matching_columns9 + ddd + 1
                                    matching_columns9 = np.intersect1d(matching_columns8,matching_columns9)
                                   # print(len(matching_columns9))
                                    for dddd in matching_columns9:
                                        matching_columns10=np.where(two_dim_matrix[dddd, dddd+1:] == 1)[0]
                                        matching_columns10 = matching_columns10 + dddd + 1
                                        matching_columns10 = np.intersect1d(matching_columns9,matching_columns10)
                                       # print(len(matching_columns10))
                                        for ww in matching_columns10:
                                            matching_columns11=np.where(two_dim_matrix[ww, ww+1:] == 1)[0]
                                            matching_columns11 = matching_columns11 + ww + 1
                                            matching_columns11 = np.intersect1d(matching_columns10,matching_columns11)
                                           # print(len(matching_columns11))
                                            for www in matching_columns11:
                                                matching_columns12=np.where(two_dim_matrix[www, www+1:] == 1)[0]
                                                matching_columns12 = matching_columns12 + www + 1
                                                matching_columns12 = np.intersect1d(matching_columns11,matching_columns12)
                                               # print(len(matching_columns12))
                                                for wwww in matching_columns12:
                                                    matching_columns13=np.where(two_dim_matrix[wwww, wwww+1:] == 1)[0]
                                                    matching_columns13 = matching_columns13 + wwww + 1
                                                    listaaa=[0,j,jj,jjj,jjjj,gg,ggg,gggg,dd,ddd,dddd,ww,www,wwww]
                                                   # print(listaaa)
                                                    allnumber += 1
                                                    #print(allnumber)
                                                   # print(listaaa)
                                                    if allnumber % 100000 == 0:
                                                       print(allnumber)
                                                       print(listaaa)
                                                       listbbb= H_indices[listaaa]
                                                       print(listbbb)
                                                       ara=strua.copy()
                                                       for ele in listbbb:
                                                           ara.replace(ele,'O')
                                                       ara.sort()
                                                       Vasp_Str = Poscar(ara)
                                                       Vasp_Str.write_file('%d.vasp'% allnumber)
                                                       
 #                                                   for wwwww in matching_columns14:
#                                                        matching_columns15=np.where(two_dim_matrix[wwwww, wwwww+1:] == 1)[0]
 #                                                       matching_columns15 = matching_columns15 + wwwww + 1

image-20231028104117829

image-20231028104128748

image-20231028104147185

image-20231028113559666

转载请注明来源 有问题可通过github提交issue

©2024 xaz

Built with Hexo and 3-hexo theme