"""Usage: python v113-rdkit-reference.py > v113-rdkit-reference.json"""
import json
import math
import rdkit
from rdkit import Chem
from rdkit.Chem import rdMolDescriptors
from rdkit.Geometry import Point3D

assert rdkit.__version__ == "2026.09.1"
Chem.SetAllowNontetrahedralChirality(True)
a = 1.01 / math.sqrt(3)
cases = [
    ("[OH-]", [[0, 0, 0], [0.97, 0, 0]]),
    ("[NH4+]", [[0, 0, 0], [a, a, a], [a, -a, -a], [-a, a, -a], [-a, -a, a]]),
]
eem = []
for smiles, coords in cases:
    mol = Chem.AddHs(Chem.MolFromSmiles(smiles))
    assert mol.GetNumAtoms() == len(coords)
    conformer = Chem.Conformer(mol.GetNumAtoms())
    for index, xyz in enumerate(coords):
        conformer.SetAtomPosition(index, Point3D(*xyz))
    mol.AddConformer(conformer)
    charges = list(rdMolDescriptors.CalcEEMcharges(mol))
    eem.append({"smiles": smiles, "formal_charge": Chem.GetFormalCharge(mol),
                "coordinates_angstrom": coords, "charges": charges,
                "total_charge": sum(charges)})

allene = []
for smiles in ["CC=[C@AL1]=CCC", "CC=[C@AL2]=CCC", "CC=[C@]=CCC", "CC=[C@@]=CCC"]:
    mol = Chem.MolFromSmiles(smiles)
    allene.append({"input": smiles,
                   "canonical": Chem.MolToSmiles(mol, canonical=True, isomericSmiles=True)})

print(json.dumps({"rdkit_version": rdkit.__version__,
                  "allow_non_tetrahedral_stereo": Chem.GetAllowNontetrahedralChirality(),
                  "eem": eem, "allene": allene}, indent=2))