from rdkit import Chem
from rdkit.Chem import rdMolDescriptors

# Function to compute LabuteASA for each atom
def compute_labute_asa(mol, includeHs=True):
    # Compute LabuteASA for the molecule
    labute_asa = rdMolDescriptors.CalcLabuteASA(mol, includeHs=includeHs)
    # Get the contribution of each atom
    (atom_contribs, H_contrib) = rdMolDescriptors._CalcLabuteASAContribs(mol)
    # Print the contribution of each atom
    for i, contrib in enumerate(atom_contribs):
        print(f"Atom {i} ({mol.GetAtomWithIdx(i).GetSymbol()}): {contrib:.2f}")
    print('Hydrogen contribution', H_contrib)
    return labute_asa

smiles = "CCO"
# Without explicit H
mol = Chem.MolFromSmiles(smiles)
print("Without explicit H:")
print("-------------------")
total_asa = compute_labute_asa(mol)
print(f"Total LabuteASA: {total_asa:.2f}")
print("===")
mol2 = Chem.AddHs(mol)
print("With explicit H:")
print("-------------------")
total_asa = compute_labute_asa(mol2, includeHs=False)
print(f"Total LabuteASA: {total_asa:.2f}")`

Without explicit H:
-------------------
Atom 0 (C): 6.92
Atom 1 (C): 6.61
Atom 2 (O): 5.11
Hydrogen contribution 1.2612803353779227
Total LabuteASA: 19.90
===
With explicit H:
-------------------
Atom 0 (C): 6.88
Atom 1 (C): 6.58
Atom 2 (O): 5.11
Atom 3 (H): 1.31
Atom 4 (H): 1.31
Atom 5 (H): 1.31
Atom 6 (H): 1.31
Atom 7 (H): 1.31

Vi[idx1] += (Rj * Rj - (Ri - dij)**2) / dij;
Vi[idx2] += (Ri * Ri - (Rj - dij)**2) / dij;

Vi[idx1] += Rj * Rj - (Ri - dij)**2 / dij;
Vi[idx2] += Ri * Ri - (Rj - dij)**2 / dij;