CSM及TET,CS

This commit is contained in:
2025-12-07 20:08:19 +08:00
parent 08f5a51fc4
commit 3d44b31194
8 changed files with 1087 additions and 0 deletions

356
py/utils/CS_analyse.py Normal file
View File

@@ -0,0 +1,356 @@
from typing import List, Dict
from pymatgen.core.structure import Structure
from pymatgen.analysis.local_env import VoronoiNN
import numpy as np
def check_real(nearest):
real_nearest = []
for site in nearest:
if np.all((site.frac_coords >= 0) & (site.frac_coords <= 1)):
real_nearest.append(site)
return real_nearest
def special_check_for_3(site, nearest):
real_nearest = []
distances = []
for site2 in nearest:
distance = np.linalg.norm(np.array(site.frac_coords) - np.array(site2.frac_coords))
distances.append(distance)
sorted_indices = np.argsort(distances)
for index in sorted_indices[:3]:
real_nearest.append(nearest[index])
return real_nearest
def CS_catulate(
struct,
sp: str = 'Li',
anion: List[str] = ['O'],
tol: float = 0,
cutoff: float = 3.0,
notice: bool = False
) -> Dict[str, Dict[str, int]]:
"""
计算结构中不同类型阳离子多面体之间的共享关系(角、边、面共享)。
该函数会分别计算以下三种情况的共享数量:
1. 目标原子 vs 目标原子 (e.g., Li-Li)
2. 目标原子 vs 其他阳离子 (e.g., Li-X)
3. 其他阳离子 vs 其他阳离子 (e.g., X-Y)
参数:
struct (Structure): 输入的pymatgen结构对象。
sp (str): 目标元素符号,默认为 'Li'
anion (list): 阴离子元素符号列表,默认为 ['O']。
tol (float): VoronoiNN 的容差。对于Li通常设为0。
cutoff (float): VoronoiNN 的截断距离。对于Li通常设为3.0。
notice (bool): 是否打印详细的共享信息。
返回:
dict: 一个字典,包含三类共享关系的统计结果。
"sp_vs_sp", "sp_vs_other", "other_vs_other" 分别对应上述三种情况。
每个键的值是另一个字典统计了共享2个(边)、3个(面)等情况的数量。
例如: {'sp_vs_sp': {'1': 10, '2': 4}, 'sp_vs_other': ...}
共享1个阴离子为角共享2个为边共享3个为面共享。
"""
# 初始化 VoronoiNN 对象
voro_nn = VoronoiNN(tol=tol, cutoff=cutoff)
# 1. 分类存储所有阳离子的近邻阴离子信息
target_sites_info = []
other_cation_sites_info = []
for index, site in enumerate(struct.sites):
# 跳过阴离子本身
if site.species.chemical_system in anion:
continue
# 获取当前位点的近邻阴离子
try:
# 使用 get_nn_info 更直接
nn_info = voro_nn.get_nn_info(struct, index)
nearest_anions = [
nn["site"] for nn in nn_info
if nn["site"].species.chemical_system in anion
]
except Exception as e:
print(f"Warning: Could not get neighbors for site {index} ({site.species_string}): {e}")
continue
if not nearest_anions:
continue
# 整理信息
site_info = {
'index': index,
'element': site.species.chemical_system,
'nearest_anion_indices': {nn.index for nn in nearest_anions}
}
# 根据是否为目标原子进行分类
if site.species.chemical_system == sp:
target_sites_info.append(site_info)
else:
other_cation_sites_info.append(site_info)
# 2. 初始化结果字典
# 共享数量key: 1-角, 2-边, 3-面
results = {
"sp_vs_sp": {"1": 0, "2": 0, "3": 0, "4": 0},
"sp_vs_other": {"1": 0, "2": 0, "3": 0, "4": 0},
"other_vs_other": {"1": 0, "2": 0, "3": 0, "4": 0},
}
# 3. 计算不同类别之间的共享关系
# 3.1 目标原子 vs 目标原子 (sp_vs_sp)
for i in range(len(target_sites_info)):
for j in range(i + 1, len(target_sites_info)):
atom_i = target_sites_info[i]
atom_j = target_sites_info[j]
shared_anions = atom_i['nearest_anion_indices'].intersection(atom_j['nearest_anion_indices'])
shared_count = len(shared_anions)
if shared_count > 0 and str(shared_count) in results["sp_vs_sp"]:
results["sp_vs_sp"][str(shared_count)] += 1
if notice:
print(
f"[Li-Li] Atom {atom_i['index']} and {atom_j['index']} share {shared_count} anions: {shared_anions}")
# 3.2 目标原子 vs 其他阳离子 (sp_vs_other)
for atom_sp in target_sites_info:
for atom_other in other_cation_sites_info:
shared_anions = atom_sp['nearest_anion_indices'].intersection(atom_other['nearest_anion_indices'])
shared_count = len(shared_anions)
if shared_count > 0 and str(shared_count) in results["sp_vs_other"]:
results["sp_vs_other"][str(shared_count)] += 1
if notice:
print(
f"[Li-Other] Atom {atom_sp['index']} and {atom_other['index']} share {shared_count} anions: {shared_anions}")
# 3.3 其他阳离子 vs 其他阳离子 (other_vs_other)
for i in range(len(other_cation_sites_info)):
for j in range(i + 1, len(other_cation_sites_info)):
atom_i = other_cation_sites_info[i]
atom_j = other_cation_sites_info[j]
shared_anions = atom_i['nearest_anion_indices'].intersection(atom_j['nearest_anion_indices'])
shared_count = len(shared_anions)
if shared_count > 0 and str(shared_count) in results["other_vs_other"]:
results["other_vs_other"][str(shared_count)] += 1
if notice:
print(
f"[Other-Other] Atom {atom_i['index']} and {atom_j['index']} share {shared_count} anions: {shared_anions}")
return results
def CS_catulate_old(struct, sp='Li', anion=['O'], tol=0, cutoff=3.0,notice=False,ID=None):
"""
计算结构中目标元素与最近阴离子的共享关系。
参数:
struct (Structure): 输入结构。
sp (str): 目标元素符号,默认为 'Li'
anion (list): 阴离子列表,默认为 ['O']。
tol (float): VoronoiNN 的容差,默认为 0。
cutoff (float): VoronoiNN 的截断距离,默认为 3.0。
返回:
list: 包含每个目标位点及其最近阴离子索引的列表。
"""
# 初始化 VoronoiNN 对象
if sp=='Li':
tol = 0
cutoff = 3.0
voro_nn = VoronoiNN(tol=tol, cutoff=cutoff)
# 初始化字典,用于统计共享关系
shared_count = {"2": 0, "3": 0,"4":0,"5":0,"6":0}
# 存储结果的列表
atom_dice = []
# 遍历结构中的每个位点
for index,site in enumerate(struct.sites):
# 跳过阴离子位点
if site.species.chemical_system in anion:
continue
# 跳过Li原子
if site.species.chemical_system == sp:
continue
# 获取 Voronoi 多面体信息
voro_info = voro_nn.get_voronoi_polyhedra(struct, index)
# 找到最近的阴离子位点
nearest_anions = [
nn_info["site"] for nn_info in voro_info.values()
if nn_info["site"].species.chemical_system in anion
]
# 如果没有找到最近的阴离子,跳过
if not nearest_anions:
print(f"No nearest anions found for {ID} site {index}.")
continue
if site.species.chemical_system == 'B' or site.species.chemical_system == 'N':
nearest_anions = special_check_for_3(site,nearest_anions)
nearest_anions = check_real(nearest_anions)
# 将结果添加到 atom_dice 列表中
atom_dice.append({
'index': index,
'nearest_index': [nn.index for nn in nearest_anions]
})
# 枚举 atom_dice 中的所有原子对
for i, atom_i in enumerate(atom_dice):
for j, atom_j in enumerate(atom_dice[i + 1:], start=i + 1):
# 获取两个原子的最近阴离子索引
nearest_i = set(atom_i['nearest_index'])
nearest_j = set(atom_j['nearest_index'])
# 比较最近阴离子的交集大小
shared_count_key = str(len(nearest_i & nearest_j))
# 更新字典中的计数
if shared_count_key in shared_count:
shared_count[shared_count_key] += 1
if notice:
if shared_count_key=='2':
print(f"{atom_j['index']}{atom_i['index']}之间存在共线")
print(f"共线的阴离子为{nearest_i & nearest_j}")
if shared_count_key=='3':
print(f"{atom_j['index']}{atom_i['index']}之间存在共面")
print(f"共面的阴离子为{nearest_i & nearest_j}")
# # 最后将字典中的值除以 2因为每个共享关系被计算了两次
# for key in shared_count.keys():
# shared_count[key] //= 2
return shared_count
def CS_count(struct, sharing_results: Dict[str, Dict[str, int]], sp: str = 'Li') -> float:
"""
分析多面体共享结果,计算平均每个目标原子参与的共享阴离子数。
这个函数是 calculate_polyhedra_sharing 的配套函数。
参数:
struct (Structure): 输入的pymatgen结构对象用于统计目标原子总数。
sharing_results (dict): 来自 calculate_polyhedra_sharing 函数的输出结果。
sp (str): 目标元素符号,默认为 'Li'
返回:
float: 平均每个目标原子sp参与的共享阴离子数量。
例如结果为2.5意味着平均每个Li原子通过共享与其他阳离子
包括Li和其他阳离子连接了2.5个阴离子。
"""
# 1. 统计结构中目标原子的总数
target_atom_count = 0
for site in struct.sites:
if site.species.chemical_system == sp:
target_atom_count += 1
# 如果结构中没有目标原子直接返回0避免除以零错误
if target_atom_count == 0:
return 0.0
# 2. 计算加权的共享阴离子总数
total_shared_anions = 0
# 处理 sp_vs_sp (例如 Li-Li) 的共享
# 每个共享关系涉及两个目标原子,所以权重需要乘以 2
if "sp_vs_sp" in sharing_results:
sp_vs_sp_counts = sharing_results["sp_vs_sp"]
for num_shared_str, count in sp_vs_sp_counts.items():
num_shared = int(num_shared_str)
# 权重 = 共享阴离子数 * 涉及的目标原子数 (2) * 出现次数
total_shared_anions += num_shared * 2 * count
# 处理 sp_vs_other (例如 Li-X) 的共享
# 每个共享关系涉及一个目标原子,所以权重乘以 1
if "sp_vs_other" in sharing_results:
sp_vs_other_counts = sharing_results["sp_vs_other"]
for num_shared_str, count in sp_vs_other_counts.items():
num_shared = int(num_shared_str)
# 权重 = 共享阴离子数 * 涉及的目标原子数 (1) * 出现次数
total_shared_anions += num_shared * 1 * count
# 3. 计算平均值
# 平均每个目标原子参与的共享阴离子数 = 总的加权共享数 / 目标原子总数
average_sharing_per_atom = total_shared_anions / target_atom_count
return average_sharing_per_atom
def CS_count_old(struct, shared_count, sp='Li'):
count = 0
for site in struct.sites:
if site.species.chemical_system == sp:
count += 1 # 累加符合条件的原子数量
CS_count = 0
for i in range(2, 7): # 遍历范围 [2, 3, 4, 5]
if str(i) in shared_count: # 检查键是否存在
CS_count += shared_count[str(i)] * i # 累加计算结果
if count > 0: # 防止除以零
CS_count /= count # 平均化结果
else:
CS_count = 0 # 如果 count 为 0直接返回 0
return CS_count
def check_only_corner_sharing(sharing_results: Dict[str, Dict[str, int]]) -> int:
"""
检查目标原子(sp)是否只参与了角共享共享1个阴离子
该函数是 calculate_polyhedra_sharing 的配套函数。
参数:
sharing_results (dict): 来自 calculate_polyhedra_sharing 函数的输出结果。
返回:
int:
- 1: 如果 sp 的共享关系中,边共享(2)、面共享(3)等数量均为0
并且至少存在一个角共享(1)。
- 0: 如果 sp 存在任何边、面等共享,或者没有任何共享关系。
"""
# 提取与目标原子 sp 相关的共享数据
sp_vs_sp_counts = sharing_results.get("sp_vs_sp", {})
sp_vs_other_counts = sharing_results.get("sp_vs_other", {})
# 1. 检查是否存在任何边共享、面共享等 (共享数 > 1)
# 检查 sp-sp 的共享
for num_shared_str, count in sp_vs_sp_counts.items():
if int(num_shared_str) > 1 and count > 0:
return 0 # 发现了边/面共享,立即返回 0
# 检查 sp-other 的共享
for num_shared_str, count in sp_vs_other_counts.items():
if int(num_shared_str) > 1 and count > 0:
return 0 # 发现了边/面共享,立即返回 0
# 2. 检查是否存在至少一个角共享 (共享数 == 1)
# 运行到这里,说明已经没有任何边/面共享了。
# 现在需要确认是否真的存在角共享,而不是完全没有共享。
corner_share_sp_sp = sp_vs_sp_counts.get("1", 0) > 0
corner_share_sp_other = sp_vs_other_counts.get("1", 0) > 0
if corner_share_sp_sp or corner_share_sp_other:
return 1 # 确认只存在角共享
else:
return 0 # 没有任何共享关系,也返回 0
# structure = Structure.from_file("../raw/0921/wjy_001.cif")
# a = CS_catulate(structure,notice=True)
# b = CS_count(structure,a)
# print(f"{a}\n{b}")
# print(check_only_corner_sharing(a))

210
py/utils/analyze_env_st.py Normal file
View File

@@ -0,0 +1,210 @@
#!/usr/bin/env python
# This code extracts the lithium environment of all of lithium sites provided in a structure file.
import os, sys
import numpy as np
import scipy
import argparse
from scipy.spatial import ConvexHull
from itertools import permutations
from pymatgen.core.structure import Structure
from pymatgen.core.periodic_table import *
from pymatgen.core.composition import *
from pymatgen.ext.matproj import MPRester
from pymatgen.io.vasp.outputs import *
from pymatgen.analysis.chemenv.coordination_environments.coordination_geometry_finder import LocalGeometryFinder
from pymatgen.analysis.chemenv.coordination_environments.structure_environments import LightStructureEnvironments
from pymatgen.analysis.chemenv.coordination_environments.chemenv_strategies import SimplestChemenvStrategy
from pymatgen.analysis.chemenv.coordination_environments.coordination_geometries import *
__author__ = "KyuJung Jun"
__version__ = "0.1"
__maintainer__ = "KyuJung Jun"
__email__ = "kjun@berkeley.edu"
__status__ = "Development"
'''
Input for the script : path to the structure file supported by Pymatgen
Structures with partial occupancy should be ordered or modified to full occupancy by Pymatgen.
'''
parser = argparse.ArgumentParser()
parser.add_argument('structure', help='path to the structure file supported by Pymatgen', nargs='?')
parser.add_argument('envtype', help='both, tet, oct, choosing which perfect environment to reference to', nargs='?')
args = parser.parse_args()
class HiddenPrints:
'''
class to reduce the output lines
'''
def __enter__(self):
self._original_stdout = sys.stdout
sys.stdout = open(os.devnull, 'w')
def __exit__(self, exc_type, exc_val, exc_tb):
sys.stdout.close()
sys.stdout = self._original_stdout
def non_elements(struct, sp='Li'):
"""
struct : 必须是一个有序结构
sp : the mobile specie
returns a new structure containing only the framework anions (O, S, N).
"""
anions_to_keep = {"O", "S", "N"}
stripped = struct.copy()
species_to_remove = [el.symbol for el in stripped.composition.elements
if el.symbol not in anions_to_keep]
if species_to_remove:
stripped.remove_species(species_to_remove)
return stripped
def site_env(coord, struct, sp="Li", envtype='both'):
'''
coord : Fractional coordinate of the target atom
struct : structure object from Pymatgen
sp : the mobile specie
envtype : This sets the reference perfect structure. 'both' compares CSM_tet and CSM_oct and assigns to the lower one.
'tet' refers to the perfect tetrahedron and 'oct' refers to the perfect octahedron
result : a dictionary of environment information
'''
stripped = non_elements(struct)
with_li = stripped.copy()
with_li.append(sp, coord, coords_are_cartesian=False, validate_proximity=False)
with_li = with_li.get_sorted_structure()
tet_oct_competition = []
if envtype == 'both' or envtype == 'tet':
for dist in np.linspace(1, 4, 601):
neigh = with_li.get_neighbors(with_li.sites[0], dist)
if len(neigh) < 4:
continue
elif len(neigh) > 4:
break
neigh_coords = [i.coords for i in neigh]
with HiddenPrints():
lgf = LocalGeometryFinder(only_symbols=["T:4"])
lgf.setup_structure(structure=with_li)
lgf.setup_local_geometry(isite=0, coords=neigh_coords)
try:
site_volume = ConvexHull(neigh_coords).volume
tet_env_list = []
for i in range(20):
tet_env = {'csm': lgf.get_coordination_symmetry_measures()['T:4']['csm'], 'vol': site_volume,
'type': 'tet'}
tet_env_list.append(tet_env)
tet_env = min(tet_env_list, key=lambda x: x['csm'])
tet_oct_competition.append(tet_env)
except Exception as e:
print(e)
print("This site cannot be recognized as tetrahedral site")
if len(neigh) == 4:
break
if envtype == 'both' or envtype == 'oct':
for dist in np.linspace(1, 4, 601):
neigh = with_li.get_neighbors(with_li.sites[0], dist)
if len(neigh) < 6:
continue
elif len(neigh) > 6:
break
neigh_coords = [i.coords for i in neigh]
with HiddenPrints():
lgf = LocalGeometryFinder(only_symbols=["O:6"], permutations_safe_override=False)
lgf.setup_structure(structure=with_li)
lgf.setup_local_geometry(isite=0, coords=neigh_coords)
try:
site_volume = ConvexHull(neigh_coords).volume
oct_env_list = []
for i in range(20):
'''
20 times sampled in case of the algorithm "APPROXIMATE_FALLBACK" is used. Large number of permutations
are performed, but the default value in the function "coordination_geometry_symmetry_measures_fallback_random"
(NRANDOM=10) is often too small. This is not a problem if algorithm of "SEPARATION_PLANE" is used.
'''
oct_env = {'csm': lgf.get_coordination_symmetry_measures()['O:6']['csm'], 'vol': site_volume,
'type': 'oct'}
oct_env_list.append(oct_env)
oct_env = min(oct_env_list, key=lambda x: x['csm'])
tet_oct_competition.append(oct_env)
except Exception as e:
print(e)
print("This site cannot be recognized as octahedral site")
if len(neigh) == 6:
break
if len(tet_oct_competition) == 0:
return {'csm': np.nan, 'vol': np.nan, 'type': 'Non_' + envtype}
elif len(tet_oct_competition) == 1:
return tet_oct_competition[0]
elif len(tet_oct_competition) == 2:
csm1 = tet_oct_competition[0]
csm2 = tet_oct_competition[1]
if csm1['csm'] > csm2['csm']:
return csm2
else:
return csm1
def extract_sites(struct, sp="Li", envtype='both'):
"""
struct : structure object from Pymatgen
envtype : 'tet', 'oct', or 'both'
sp : target element to analyze environment
"""
envlist = []
# --- 关键修改:直接遍历原始结构,即使它是无序的 ---
# 我们不再调用 get_sorted_structure()
# 我们只关心那些含有目标元素 sp 的位点
# 遍历每一个位点 (site)
for i, site in enumerate(struct):
# 检查当前位点的组分(site.species)中是否包含我们感兴趣的元素(sp)
# site.species.elements 返回该位点上的元素列表,例如 [Element Li, Element Fe]
# [el.symbol for el in site.species.elements] 将其转换为符号列表 ['Li', 'Fe']
site_elements = [el.symbol for el in site.species.elements]
if sp in site_elements:
# 如果找到了Li我们就对这个位点进行环境分析
# 注意:我们将原始的、可能无序的 struct 传递给 site_env
# 因为 site_env 内部的函数 (如 LocalGeometryFinder) 知道如何处理它
# 为了让下游函数(特别是 non_elements能够工作
# 我们在这里创建一个一次性的、临时的有序结构副本给它
# 这可以避免我们之前遇到的所有 'ordered structures only' 错误
temp_ordered_struct = struct.get_sorted_structure()
singleenv = site_env(site.frac_coords, temp_ordered_struct, sp, envtype)
envlist.append({'frac_coords': site.frac_coords, 'type': singleenv['type'], 'csm': singleenv['csm'],
'volume': singleenv['vol']})
if not envlist:
print(f"警告: 在结构中未找到元素 {sp} 的占位。")
return envlist
def export_envs(envlist, sp='Li', envtype='both', fname=None):
'''
envlist : list of dictionaries of environment information
fname : Output file name
'''
if not fname:
fname = "extracted_environment_info" + "_" + sp + "_" + envtype + ".dat"
with open(fname, 'w') as f:
f.write('List of environment information\n')
f.write('Species : ' + sp + "\n")
f.write('Envtype : ' + envtype + "\n")
for index, i in enumerate(envlist):
f.write("Site index " + str(index) + ": " + str(i) + '\n')
# struct = Structure.from_file("../raw/0921/wjy_475.cif")
# site_info = extract_sites(struct, envtype="both")
# export_envs(site_info, sp="Li", envtype="both")