#!/usr/bin/env python
"""
Protonates the input molecules...

@author: JP

"""

from rdkit import Chem
from rdkit.Chem import AllChem 


molstr = """MolPort-006-830-028
     RDKit          2D

 26 30  0  0  0  0  0  0  0  0999 V2000
   10.4974    7.8006    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    9.7943    6.5019    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
    9.0042    7.8254    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    9.0166    6.9238    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    6.8866    6.5681    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
    9.7198    8.2515    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   10.5098    6.9196    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    7.5566    6.9362    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    6.2413    6.9320    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    8.3177    6.5019    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    5.6003    5.8070    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
    6.8990    5.8277    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    9.7859    5.6747    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    6.2537    5.4390    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    8.2639    8.2060    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   11.2129    8.2266    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    7.5566    7.7841    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    5.9808    6.1751    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    5.6003    6.5598    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    9.3682    4.9633    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   10.1913    4.9633    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    9.7404    9.0828    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
   11.2046    9.1283    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
    6.8452    8.2597    0.0000 F   0  0  0  0  0  0  0  0  0  0  0  0
   11.9409    7.7717    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
    4.9633    5.4306    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
  2  7  1  0
  3  6  1  0
  4  2  1  0
  5  8  1  0
  6  1  1  0
  7  1  2  0
  8 10  1  0
  9  5  1  0
 10  4  2  0
 14 11  1  0
 12  5  1  0
 13  2  1  0
 14 12  1  0
 15  3  2  0
 16  1  1  0
 17 15  1  0
  9 18  1  0
  9 19  1  0
 20 13  1  0
 21 13  1  0
 22  6  2  0
 23 16  2  0
 24 17  1  0
 25 16  1  0
 26 11  1  0
  4  3  1  0
 20 21  1  0
 17  8  2  0
 19 11  1  0
 14 18  1  0
M  END
$$$$
"""


rules=(
# Carboxylic Acid ionization
   ('[O-]','[$([OH]C(=O))]'),

# Now the bases

# Amines 3 
   ('[NH1+]','[$([NX3;H0]([CX4])([CX4])[CX4])]')
)


# Load my charming molecule
mol = Chem.MolFromMolBlock(molstr)
# sanitize
Chem.SanitizeMol(mol)

orig = Chem.MolToSmiles(mol)
for repl, patt in rules:
    repl = Chem.MolFromSmiles(repl, sanitize=False)
    patt = Chem.MolFromSmarts(patt)
    ps = AllChem.ReplaceSubstructs(mol, patt, repl, replaceAll=True)
    ps[0].UpdatePropertyCache(False)
    Chem.GetSymmSSSR(ps[0])                
    new = Chem.MolToSmiles(ps[0])
    modified = Chem.MolToSmiles(mol) != Chem.MolToSmiles(ps[0])
    if modified:
        # emebed these as
        print "Replaced" 
        #AllChem.Compute2DCoords(ps[0])
        orig = Chem.MolToSmiles(mol)
        mol = ps[0]
## Finally, clean up and write smiles and sdf versions
Chem.SanitizeMol(mol)

print Chem.MolToMolBlock(mol)