330 lines
12 KiB
Python
330 lines
12 KiB
Python
from fontTools.misc.plistlib import end_string
|
||
from pymatgen.core import Structure
|
||
import spglib
|
||
from fractions import Fraction
|
||
import random
|
||
from pymatgen.core.sites import PeriodicSite, Species,Element,Lattice
|
||
import numpy as np
|
||
import re
|
||
def typejudge(number):
|
||
if number in [1, 2]:
|
||
return "Triclinic"
|
||
elif 3 <= number <= 15:
|
||
return "Monoclinic"
|
||
elif 16 <= number <= 74:
|
||
return "Orthorhombic"
|
||
elif 75 <= number <= 142:
|
||
return "Tetragonal"
|
||
elif 143 <= number <= 167:
|
||
return "Trigonal"
|
||
elif 168 <= number <= 194:
|
||
return "Hexagonal"
|
||
elif 195 <= number <= 230:
|
||
return "Cubic"
|
||
else:
|
||
return "Unknown"
|
||
|
||
|
||
def extract_oxi_state(species_string):
|
||
"""
|
||
从 species_string（如 "Li+:0.645"）提取氧化态
|
||
|
||
Args:
|
||
species_string: 包含元素和氧化态的字符串（如 "Li+:0.645"、"Fe2-:0.5"）
|
||
|
||
Returns:
|
||
int: 氧化态数值（如 "+" -> 1, "2-" -> -2），默认返回 0 如果解析失败
|
||
"""
|
||
# 使用正则表达式匹配化学符号和氧化态（如 Li+, Fe2-）
|
||
match = re.search(r"([A-Za-z]+)([+-]?\d*[+-])", species_string)
|
||
if not match:
|
||
return 0 # 默认中性
|
||
|
||
# 提取氧化态部分（如 "+", "2-", "3+"）
|
||
oxi_str = match.group(2)
|
||
|
||
# 处理单符号情况（如 "+" -> 1, "-" -> -1）
|
||
if oxi_str in ("+", "-"):
|
||
return 1 if oxi_str == "+" else -1
|
||
|
||
# 处理多数字情况（如 "2+" -> 2, "3-" -> -3）
|
||
try:
|
||
num = int(oxi_str[:-1]) # 提取数字部分
|
||
sign = 1 if oxi_str[-1] == "+" else -1
|
||
return num * sign
|
||
except (ValueError, IndexError):
|
||
return 0 # 解析失败时返回中性
|
||
|
||
def process_cif_file(cif_file_path,explict_element):
|
||
structure = Structure.from_file(cif_file_path)
|
||
result_list = []
|
||
for index, site in enumerate(structure.sites, start=1):
|
||
occupancy = site.species.element_composition.num_atoms
|
||
species = site.species.chemical_system
|
||
if species in explict_element:
|
||
break
|
||
if occupancy < 1.0:
|
||
entry = next((r for r in result_list if
|
||
r["species"] == species and r["occupancy"] == occupancy), None)
|
||
if entry:
|
||
entry["number"].append(index)
|
||
else:
|
||
result_list.append({
|
||
"species": species,
|
||
"number": [index],
|
||
"occupancy": occupancy
|
||
})
|
||
return result_list
|
||
|
||
|
||
def factorize_to_three_factors(n,type_sym=None):
|
||
factors = []
|
||
|
||
# 遍历可能的x值
|
||
if type_sym == None:
|
||
for x in range(1, n + 1):
|
||
if n % x == 0:
|
||
remaining_n = n // x
|
||
# 遍历可能的y值
|
||
for y in range(1, remaining_n + 1):
|
||
if remaining_n % y == 0:
|
||
z = remaining_n // y
|
||
factors.append({'x': x, 'y': y, 'z': z})
|
||
if type_sym == "xyz":
|
||
for x in range(1, n + 1):
|
||
if n % x == 0:
|
||
remaining_n = n // x
|
||
# 遍历可能的y值
|
||
for y in range(1, remaining_n + 1):
|
||
if remaining_n % y == 0 and y <= x:
|
||
z = remaining_n // y
|
||
if z <= y:
|
||
factors.append({'x': x, 'y': y, 'z': z})
|
||
|
||
def sum_score(factor):
|
||
x, y, z = factor['x'], factor['y'], factor['z']
|
||
return x + y + z
|
||
|
||
# 按照sum_score从小到大排序
|
||
sorted_factors = sorted(factors, key=sum_score)
|
||
return sorted_factors
|
||
|
||
|
||
|
||
|
||
|
||
def calculate_supercell_factor(occupancy):
|
||
# 将浮点数转换为分数形式
|
||
fraction = Fraction(occupancy).limit_denominator()
|
||
|
||
# 获取分子和分母
|
||
numerator = fraction.numerator
|
||
denominator = fraction.denominator
|
||
|
||
return numerator,denominator
|
||
|
||
|
||
def mark_atoms_randomly(factors, atom_number):
|
||
"""
|
||
根据扩胞因子和占据数量生成随机占据字典
|
||
|
||
Args:
|
||
factors: 扩胞因子字典 {'x': int, 'y': int, 'z': int}
|
||
atom_number: 需要占据的副本数量
|
||
|
||
Returns:
|
||
字典 {0: 1或0, 1: 1或0, ..., total_copies-1: 1或0}
|
||
"""
|
||
x, y, z = factors['x'], factors['y'], factors['z']
|
||
total_copies = x * y * z
|
||
|
||
if atom_number > total_copies:
|
||
raise ValueError(f"atom_number ({atom_number}) 不能超过扩胞总数 (x*y*z = {total_copies})")
|
||
|
||
# 生成所有副本索引 [0, 1, 2, ..., total_copies-1]
|
||
atom_dice = list(range(total_copies))
|
||
|
||
# 随机选择 atom_number 个副本占据
|
||
selected_atoms = random.sample(atom_dice, atom_number)
|
||
|
||
# 创建结果字典 {0: 1或0, 1: 1或0, ...}
|
||
result = {atom: 1 if atom in selected_atoms else 0 for atom in atom_dice}
|
||
|
||
return result
|
||
|
||
|
||
def generate_random_list(total_elements, atom_number):
|
||
# 确保 atom_number 不超过 total_elements
|
||
if atom_number > total_elements:
|
||
raise ValueError("atom_number cannot be greater than the total number of elements (x * y * z)")
|
||
|
||
# 创建一个全0的列表
|
||
result = [0] * total_elements
|
||
|
||
# 随机选择 atom_number 个位置，并将这些位置的值设为1
|
||
indices = random.sample(range(total_elements), atom_number)
|
||
for index in indices:
|
||
result[index] = 1
|
||
|
||
return result
|
||
|
||
|
||
def merge_structures(struct_copies, factors):
|
||
"""
|
||
将多个副本结构按三维顺序合并为扩胞后的结构
|
||
|
||
Args:
|
||
struct_copies: 副本结构列表（长度 = x*y*z）
|
||
factors: 扩胞因子字典 {"x": int, "y": int, "z": int}
|
||
|
||
Returns:
|
||
合并后的扩胞结构
|
||
"""
|
||
x, y, z = factors["x"], factors["y"], factors["z"]
|
||
total_copies = x * y * z
|
||
|
||
if len(struct_copies) != total_copies:
|
||
raise ValueError("副本数量与扩胞因子不匹配")
|
||
|
||
# 获取原结构的晶格
|
||
original_lattice = struct_copies[0].lattice
|
||
|
||
# 创建扩胞后的新晶格（直接按倍数缩放）
|
||
new_lattice_matrix = np.dot(original_lattice.matrix, np.diag([x, y, z]))
|
||
new_lattice = Lattice(new_lattice_matrix)
|
||
|
||
# 初始化合并后的结构
|
||
merged_structure = Structure(
|
||
lattice=new_lattice,
|
||
species=[],
|
||
coords=[],
|
||
coords_are_cartesian=False
|
||
)
|
||
|
||
# 按三维顺序填充每个副本的原子
|
||
for copy_idx in range(total_copies):
|
||
# 计算当前副本的分数坐标偏移量
|
||
offset = np.array([
|
||
copy_idx // (y * z), # x方向偏移
|
||
(copy_idx % (y * z)) // z, # y方向偏移
|
||
copy_idx % z # z方向偏移
|
||
])
|
||
|
||
# 将当前副本的原子添加到合并结构中（考虑偏移）
|
||
for site in struct_copies[copy_idx]:
|
||
if site.species: # 跳过空位
|
||
merged_structure.append(
|
||
species=site.species,
|
||
coords=site.frac_coords + offset,
|
||
coords_are_cartesian=False,
|
||
properties=site.properties
|
||
)
|
||
|
||
return merged_structure
|
||
|
||
def expand_structure(structure, factors, atom_indices, atom_number):
|
||
# 参数检查
|
||
x, y, z = factors['x'], factors['y'], factors['z']
|
||
total_copies = x * y * z
|
||
if not all(1 <= idx <= len(structure.sites) for idx in atom_indices):
|
||
raise ValueError("atom_indices包含无效原子索引")
|
||
|
||
# 生成独立副本
|
||
struct_copies = [structure.copy() for _ in range(total_copies)]
|
||
atom_dice = list(range(total_copies)) # 所有副本索引 [0,1,2,...]
|
||
|
||
# 处理每个目标原子
|
||
for atom_idx in atom_indices:
|
||
original_site = structure.sites[atom_idx - 1]
|
||
element = original_site.species.chemical_system
|
||
|
||
|
||
# 生成当前原子的占据字典（如{0:1, 1:0, 2:1,...}）
|
||
occupancy_dict = mark_atoms_randomly(factors,atom_number)
|
||
# 修改每个副本
|
||
for copy_idx, occupy in occupancy_dict.items():
|
||
# 或者方法2：使用remove/insert
|
||
struct_copies[copy_idx].remove_sites([atom_idx - 1])
|
||
oxi_state = extract_oxi_state(original_site.species_string)
|
||
if occupy:
|
||
new_site = PeriodicSite(
|
||
species=Species(element, oxi_state),
|
||
coords=original_site.frac_coords,
|
||
lattice=struct_copies[copy_idx].lattice,
|
||
to_unit_cell=True,
|
||
label=original_site.label
|
||
)
|
||
struct_copies[copy_idx].sites.insert(atom_idx - 1, new_site)
|
||
else:
|
||
species_dict = {Species(element, oxi_state): 0.0}
|
||
new_site = PeriodicSite(
|
||
species = species_dict,
|
||
coords=original_site.frac_coords,
|
||
lattice=struct_copies[copy_idx].lattice,
|
||
to_unit_cell=True,
|
||
label=original_site.label
|
||
)
|
||
struct_copies[copy_idx].sites.insert(atom_idx - 1, new_site)
|
||
|
||
# 合并副本
|
||
expanded_structure = Structure(
|
||
lattice=np.dot(structure.lattice.matrix, np.diag([x, y, z])),
|
||
species=[],
|
||
coords=[],
|
||
coords_are_cartesian=False
|
||
)
|
||
|
||
for copy_idx in range(total_copies):
|
||
offset = np.array([
|
||
copy_idx // (y * z),
|
||
(copy_idx % (y * z)) // z,
|
||
copy_idx % z
|
||
])
|
||
for site in struct_copies[copy_idx]:
|
||
if site.species: # 只添加非空位
|
||
expanded_structure.append(
|
||
species=site.species,
|
||
coords=site.frac_coords + offset,
|
||
coords_are_cartesian=False,
|
||
properties=site.properties
|
||
)
|
||
expanded_structure = merge_structures(struct_copies,factors)
|
||
return expanded_structure
|
||
|
||
|
||
def process_occupation(input_file,output_file,explict_element = ["Li"],expect_cifnumber = 10,random_time=1):
|
||
struct = Structure.from_file(input_file)
|
||
space_group_info = struct.get_space_group_info()
|
||
space_group_symbol = space_group_info[0]
|
||
all_spacegroup_symbols = [spglib.get_spacegroup_type(i) for i in range(1, 531)]
|
||
symbol = all_spacegroup_symbols[0]
|
||
for symbol_i in all_spacegroup_symbols:
|
||
if space_group_symbol == symbol_i.international_short:
|
||
symbol = symbol_i
|
||
space_type = typejudge(symbol.number)
|
||
print(f"当前空间群符号为{space_group_symbol},序号为{symbol.number},对应的晶体体系为{space_type}")
|
||
occupation_list = process_cif_file(cif_file_path=input_file,explict_element=explict_element)
|
||
print(occupation_list)
|
||
for occupation in occupation_list:
|
||
atom_number, target_multiplier=calculate_supercell_factor(occupation["occupancy"])
|
||
divides = []
|
||
if space_type == "Hexagonal":
|
||
print('当前为六方晶系，暂不处理')
|
||
if space_type == "Cubic":
|
||
print("当前为立方晶体，三个方向同步")
|
||
divides = factorize_to_three_factors(target_multiplier,"xyz")
|
||
else:
|
||
print("为其他晶系，假设三个方向不同")
|
||
divides = factorize_to_three_factors(target_multiplier,)
|
||
print(divides)
|
||
for it in divides:
|
||
end_str = f'x{it["x"]}y{it["y"]}z{it["z"]}'
|
||
for i in range(random_time):
|
||
expand_struct=expand_structure(struct,it,occupation["number"],atom_number)
|
||
expand_struct.to_file(output_file)
|
||
print(it)
|
||
else:
|
||
print(f"不存在除{explict_element}以外的共占位原子")
|
||
process_occupation("../data/input_pre/ICSD_1234.cif", "haha.cif", explict_element=[], expect_cifnumber=1, random_time=1)
|
||
|