# Tested with RDKit 2026.03.6. Run: python compare_enantiomers.py enantiomers_pubchem.json
import json, csv, sys
from pathlib import Path
import rdkit
from rdkit import Chem, DataStructs
from rdkit.Chem import Descriptors, rdMolDescriptors, rdFingerprintGenerator
source = Path(sys.argv[1]) if len(sys.argv)>1 else Path('enantiomers_pubchem.json')
raw=json.loads(source.read_text())
records=raw['records']
assert len(records)==40
gens = {
 'morgan_off':rdFingerprintGenerator.GetMorganGenerator(radius=2,fpSize=2048,includeChirality=False),
 'morgan_on':rdFingerprintGenerator.GetMorganGenerator(radius=2,fpSize=2048,includeChirality=True),
 'atom_pair_off':rdFingerprintGenerator.GetAtomPairGenerator(fpSize=2048,includeChirality=False,countSimulation=False),
 'atom_pair_on':rdFingerprintGenerator.GetAtomPairGenerator(fpSize=2048,includeChirality=True,countSimulation=False)
}
rows=[]
for i in range(0,40,2):
 a,b=records[i:i+2]
 pa,pb=a['result']['properties'][0],b['result']['properties'][0]
 ma,mb=[Chem.MolFromSmiles(p['SMILES']) for p in (pa,pb)]
 assert ma is not None and mb is not None
 assert pa['Charge']==pb['Charge']==0
 for m in (ma,mb): Chem.AssignStereochemistry(m,cleanIt=True,force=True)
 ca,cb=[Chem.FindMolChiralCenters(m,includeUnassigned=True,useLegacyImplementation=False) for m in (ma,mb)]
 iso=lambda m:Chem.MolToSmiles(m,isomericSmiles=True)
 plain=lambda m:Chem.MolToSmiles(m,isomericSmiles=False)
 mirror=Chem.Mol(ma)
 for atom in mirror.GetAtoms():
  if atom.GetChiralTag() in (Chem.ChiralType.CHI_TETRAHEDRAL_CW,Chem.ChiralType.CHI_TETRAHEDRAL_CCW): atom.InvertChirality()
 Chem.AssignStereochemistry(mirror,cleanIt=True,force=True)
 checked=(plain(ma)==plain(mb) and iso(ma)!=iso(mb) and iso(mirror)==iso(mb) and len(ca)==len(cb)==1 and all(s!='?' for _,s in ca+cb))
 assert checked and ca[0][1] != cb[0][1], (a['query'],ca,cb)
 row=dict(pair=a['query'].removeprefix('L-').removeprefix('(S)-'),name_A=a['query'],name_B=b['query'],
 cid_A=pa['CID'],cid_B=pb['CID'],url_A=f"https://pubchem.ncbi.nlm.nih.gov/compound/{pa['CID']}",url_B=f"https://pubchem.ncbi.nlm.nih.gov/compound/{pb['CID']}",
 smiles_A=pa['SMILES'],smiles_B=pb['SMILES'],cip_A=ca[0][1],cip_B=cb[0][1],
 stereo_centers_A=str(ca),stereo_centers_B=str(cb),enantiomer_check=checked,
 formula_A=pa['MolecularFormula'],formula_B=pb['MolecularFormula'],
 pubchem_weight_A=pa['MolecularWeight'],pubchem_weight_B=pb['MolecularWeight'],
 rdkit_formula_A=rdMolDescriptors.CalcMolFormula(ma),rdkit_formula_B=rdMolDescriptors.CalcMolFormula(mb),
 rdkit_weight_A=Descriptors.MolWt(ma),rdkit_weight_B=Descriptors.MolWt(mb),
 connectivity_equal=plain(ma)==plain(mb),isomeric_smiles_equal=iso(ma)==iso(mb))
 assert row['formula_A']==row['formula_B']==row['rdkit_formula_A']==row['rdkit_formula_B']
 assert float(row['pubchem_weight_A'])==float(row['pubchem_weight_B'])
 assert row['rdkit_weight_A']==row['rdkit_weight_B']
 for name,gen in gens.items():
  fa,fb=[gen.GetFingerprint(m) for m in (ma,mb)]
  row[name+'_tanimoto']=DataStructs.TanimotoSimilarity(fa,fb)
  row[name+'_identical']=fa.ToBitString()==fb.ToBitString()
  row[name+'_bits_A']=fa.GetNumOnBits()
  row[name+'_bits_B']=fb.GetNumOnBits()
  row[name+'_shared_bits']=(fa&fb).GetNumOnBits()
 fa,fb=[Chem.RDKFingerprint(m,fpSize=2048) for m in (ma,mb)]
 row['rdk_path_tanimoto']=DataStructs.TanimotoSimilarity(fa,fb)
 fa,fb=[gens['morgan_on'].GetSparseCountFingerprint(m) for m in (ma,mb)]
 row['morgan_on_sparse_counts_equal']=fa.GetNonzeroElements()==fb.GetNonzeroElements()
 rows.append(row)
with open('enantiomer_comparison.csv','w') as f:
 w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
print('RDKit',rdkit.__version__,'verified pairs',len(rows))
for r in rows: print(r['pair'],r['cid_A'],r['cid_B'],r['formula_A'],r['pubchem_weight_A'],r['cip_A']+'/'+r['cip_B'],*[round(r[k+'_tanimoto'],6) for k in gens],r['rdk_path_tanimoto'])
print('identical pair counts:',{k:sum(r[k+'_identical'] for r in rows) for k in gens})
