本帖最后由 Graphite 于 2025-9-24 20:36 编辑
这个算法是捣鼓reax_tools时某天看b站无意刷到的群论原理得到的灵感。
简单地说,我们给分子中每一个中心原子赋予一个大质数,例如C=701,对键连的原子赋予一个小质数,例如H=3,O=23,一个连接了H,H,O的原子C的值是701 * 3 * 3 * 23,由于乘法具有交换律,所以无论是C-H-H-O还是C-O-H-O结果都是一样的,结果称为中心原子的hash。
再将所有中心原子的hash相加,得到初步的分子hash,加法也是符合交换律的。最后再用几个数字搓几遍hash,把结果打散。
这个hash是一个usigned int,也就是说到达42.9亿后会安全溢出回到0,每个分子会有一个唯一的hash,并且结构微微改变一个H,hash就会大不一样,理想情况下,能够均匀分布在0-42.9亿的范围内。我们知道std::unordered_map和std::unordered_set利用的是哈希桶,我们使键为这个hash,值为分子对象,可以实现O(1)复杂度的分子查找和去重。
要确定体系中是否有一个分子,只要查有没有hash就可以,要维护体系中分子的模版和建立反应网络,只需要有所有的hash,和一份分子对象池就可以。
于是,在reax_tools中,最后收集和处理所有相关反应时,总时间复杂度是O(N)。N为体系原子数或分子数或轨迹文件的大小。
不过,这个算法只能处理构造异构体之类,旋光异构体不行。hash碰撞的概率不高,至少到达出现碰撞的时候,对于总体的数值结果影响已经很小了。
inline void update_topology() { // Formula std::string mol_formula; for (auto& mol_atom : mol_atoms) { types_nums[mol_atom->type_name]++; }
for (auto& pair : types_nums) { mol_formula += fmt::format("{}{}", pair.first,pair.second); } formula = rename_formula(mol_formula, ELEMENT_DISPLAY_ORDER);
// Hash unsigned int mol_hash = 1; for (const auto& atom : mol_atoms) { unsigned int atom_hash = BIGGER_PRIME_NUMBERS[atom->type_id]; unsigned int bonded_hash = 1;
for (const auto& bonded : atom->bonded_atoms){ bonded_hash *= PRIME_NUMBERS[bonded->type_id]; // multiply is commutative } // Combine atom and its bonded environment, order independent atom_hash *= bonded_hash; mol_hash += atom_hash; // Addition is commutative } // Final mixing for better distribution mol_hash ^= (mol_hash >> 16); mol_hash *= 0x85ebca6b; mol_hash ^= (mol_hash >> 13); hash = mol_hash;
for (auto& atom : mol_atoms) { atom->belong_molecule = this; } };
|