<commit_before>#-*. coding: utf-8 -*-
## Copyright (c) 2008-2011, <NAME>; 2012, <NAME>; 2014, <NAME>;
## All rights reserved.
##
##  This file is part of Cinfony.
##  The contents are covered by the terms of the BSD license
##  which is included in the file LICENSE_BSD.txt.

"""
rdkit - A Cinfony module for accessing the RDKit from CPython

Global variables:
  Chem and AllChem - the underlying RDKit Python bindings
  informats - a dictionary of supported input formats
  outformats - a dictionary of supported output formats
  descs - a list of supported descriptors
  fps - a list of supported fingerprint types
  forcefields - a list of supported forcefields
"""

import os
from copy import copy

import numpy as np

from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem import Descriptors

_descDict = dict(Descriptors.descList)

import rdkit.DataStructs
import rdkit.Chem.MACCSkeys
import rdkit.Chem.AtomPairs.Pairs
import rdkit.Chem.AtomPairs.Torsions
### ODDT ###
from rdkit.Chem.Lipinski import NumRotatableBonds
from rdkit.Chem.AllChem import ComputeGasteigerCharges
from rdkit.Chem.Pharm2D import Gobbi_Pharm2D,Generate

# PIL and Tkinter
try:
    import Tkinter as tk
    import Image as PIL
    import ImageTk as PILtk
except:
    PILtk = None

# Aggdraw
try:
    import aggdraw
    from rdkit.Chem.Draw import aggCanvas
except ImportError:
    aggdraw = None

fps = ['rdkit', 'layered', 'maccs', 'atompairs', 'torsions', 'morgan']
"""A list of supported fingerprint types"""
descs = _descDict.keys()
"""A list of supported descriptors"""

_formats = {'smi': "SMILES",
            'can': "Canonical SMILES",
            'mol': "MDL MOL file",
            'mol2': "Tripos MOL2 file",
            'sdf': "MDL SDF file",
            'inchi':"InChI",
            'inchikey':"InChIKey"}
_notinformats = ['can', 'inchikey']
_notoutformats = ['mol2']
if not Chem.INCHI_AVAILABLE:
    _notinformats += ['inchi']
    _notoutformats += ['inchi', 'inchikey']

informats = dict([(_x, _formats[_x]) for _x in _formats if _x not in _notinformats])
"""A dictionary of supported input formats"""
outformats = dict([(_x, _formats[_x]) for _x in _formats if _x not in _notoutformats])
"""A dictionary of supported output formats"""

_forcefields = {'uff': AllChem.UFFOptimizeMolecule}
forcefields = _forcefields.keys()
"""A list of supported forcefields"""

def readfile(format, filename, *args, **kwargs):
    """Iterate over the molecules in a file.

    Required parameters:
       format - see the informats variable for a list of available
                input formats
       filename

    You can access the first molecule in a file using the next() method
    of the iterator:
        mol = readfile("smi", "myfile.smi").next()

    You can make a list of the molecules in a file using:
        mols = list(readfile("smi", "myfile.smi"))

    You can iterate over the molecules in a file as shown in the
    following code snippet:
    >>> atomtotal = 0
    >>> for mol in readfile("sdf", "head.sdf"):
    ...     atomtotal += len(mol.atoms)
    ...
    >>> print atomtotal
    43
    """
    if not os.path.isfile(filename):
        raise IOError, "No such file: '%s'" % filename
    format = format.lower()
    # Eagerly evaluate the supplier functions in order to report
    # errors in the format and errors in opening the file.
    # Then switch to an iterator...
    if format=="sdf":
        def filereader_sdf():
            block = ''
            n = 0
            f = open(filename)
            for line in f:
                block += line
                if line[:4] == '$$$$':
                    yield Molecule(source={'fmt': format, 'string': block, 'n': n, 'filename': filename})
                    n += 1
                    block = ''
            f.close()
        return filereader_sdf()
    elif format=="mol":
        def mol_reader():
            yield Molecule(Chem.MolFromMolFile(filename, *args, **kwargs))
        return mol_reader()
    elif format=="pdb":
        def mol_reader():
            yield Molecule(Chem.MolFromPDBFile(filename, *args, **kwargs))
        return mol_reader()
    elif format=="mol2":
        def filereader_mol2():
            block = ''
            data = ''
            n = 0
            f = open(filename)
            for line in f:
                if line[:1] == '#':
                    data += line
                elif line[:17] == '@<TRIPOS>MOLECULE':
                    if n>0: #skip `zero` molecule (any preciding comments and spaces)
                        yield Molecule(source={'fmt': format, 'string': block, 'n': n, 'filename': filename})
                    n += 1
                    block = data
                    data = ''
                block += line
            f.close()
            
        return filereader_mol2()
    
    elif format=="smi":
        iterator = Chem.SmilesMolSupplier(filename, delimiter=" \t",
                                          titleLine=False, *args, **kwargs)
        def smi_reader():
            for mol in iterator:
                yield Molecule(mol)
        return smi_reader()
    elif format=='inchi' and Chem.INCHI_AVAILABLE:
        def inchi_reader():
            for line in open(filename, 'r'):
                mol = Chem.inchi.MolFromInchi(line.strip(), *args, **kwargs)
                yield Molecule(mol)
        return inchi_reader()
    else:
        raise ValueError, "%s is not a recognised RDKit format" % format

def readstring(format, string):
    """Read in a molecule from a string.

    Required parameters:
       format - see the informats variable for a list of available
                input formats
       string

    Example:
    >>> input = "C1=CC=CS1"
    >>> mymol = readstring("smi", input)
    >>> len(mymol.atoms)
    5
    """
    format = format.lower()
    if format=="mol" or format=="sdf":
        mol = Chem.MolFromMolBlock(string)
    elif format=="mol2":
        mol = Chem.MolFromMol2Block(string)
    elif format=="pdb":
        mol = Chem.MolFromPDBBlock(string)
    elif format=="smi":
        mol = Chem.MolFromSmiles(string)
    elif format=='inchi' and Chem.INCHI_AVAILABLE:
        mol = Chem.inchi.MolFromInchi(string)
    else:
        raise ValueError,"%s is not a recognised RDKit format" % format
    if mol:
        return Molecule(mol)
    else:
        raise IOError, "Failed to convert '%s' to format '%s'" % (
            string, format)

class Outputfile(object):
    """Represent a file to which *output* is to be sent.

    Required parameters:
       format - see the outformats variable for a list of available
                output formats
       filename

    Optional parameters:
       overwite -- if the output file already exists, should it
                   be overwritten? (default is False)

    Methods:
       write(molecule)
       close()
    """
    def __init__(self, format, filename, overwrite=False):
        self.format = format
        self.filename = filename
        if not overwrite and os.path.isfile(self.filename):
            raise IOError, "%s already exists. Use 'overwrite=True' to overwrite it." % self.filename
        if format=="sdf":
            self._writer = Chem.SDWriter(self.filename)
        elif format=="smi":
            self._writer = Chem.SmilesWriter(self.filename, isomericSmiles=True)
        elif format in ('inchi', 'inchikey') and Chem.INCHI_AVAILABLE:
            self._writer= open(filename, 'w')
        else:
            raise ValueError,"%s is not a recognised RDKit format" % format
        self.total = 0 # The total number of molecules written to the file

    def write(self, molecule):
        """Write a molecule to the output file.

        Required parameters:
           molecule
        """
        if not self.filename:
            raise IOError, "Outputfile instance is closed."
        if self.format in ('inchi', 'inchikey'):
            self._writer.write(molecule.write(self.format) +'\n')
        else:
            self._writer.write(molecule.Mol)
        self.total += 1

    def close(self):
        """Close the Outputfile to further writing."""
        self.filename = None
        self._writer.flush()
        del self._writer

class Molecule(object):
    """Represent an rdkit Molecule.

    Required parameter:
       Mol -- an RDKit Mol or any type of cinfony Molecule

    Attributes:
       atoms, data, formula, molwt, title

    Methods:
       addh(), calcfp(), calcdesc(), draw(), localopt(), make3D(), removeh(),
       write()

    The underlying RDKit Mol can be accessed using the attribute:
       Mol
    """
    _cinfony = True

    def __init__(self, Mol = None, source = None, protein = False):
        if hasattr(Mol, "_cinfony"):
            a, b = Mol._exchange
            if a == 0:
                molecule = readstring("smi", b)
            else:
                molecule = readstring("mol", b)
            Mol = molecule.Mol

        self.Mol = Mol
        ### ODDT ###
        self.protein = protein
        self._atom_dict = None
        self._res_dict = None
        self._ring_dict = None
        self._coords = None
        self._charges = None
        # lazy
        self._source = source # dict with keys: n, fmt, string, filename
    
    # lazy Molecule parsing requires masked Mol
    @property
    def Mol(self):
        if not self._Mol and self._source:
            self._Mol = readstring(self._source['fmt'], self._source['string']).Mol
            self._source = None
        return self._Mol
        
    @Mol.setter
    def Mol(self, value):
        self._Mol = value
    
    @property
    def atoms(self): return [Atom(rdkatom) for rdkatom in self.Mol.GetAtoms()]
    @property
    def data(self): return MoleculeData(self.Mol)
    @property
    def molwt(self): return Descriptors.MolWt(self.Mol)
    @property
    def formula(self): return Descriptors.MolecularFormula(self.Mol)
    def _gettitle(self):
        # Note to self: maybe should implement the get() method for self.data
        if "_Name" in self.data:
            return self.data["_Name"]
        else:
            return ""
    def _settitle(self, val): self.Mol.SetProp("_Name", val)
    title = property(_gettitle, _settitle)
    @property
    def _exchange(self):
        if self.Mol.GetNumConformers() == 0:
            return (0, self.write("smi"))
        else:
            return (1, self.write("mol"))
    
    # cache frequently used properties and cache them in prefixed [_] variables
    @property
    def coords(self):
        if self._coords is None:
            self._coords = np.array([atom.coords for atom in self.atoms])
        return self._coords
    
    @property
    def charges(self):
        if self._charges is None:
            self._charges = np.array([atom.partialcharge for atom in self.atoms])
        return self._charges
    
    #### Custom ODDT properties ####
    @property
    def sssr(self):
        return [list(path) for path in list(Chem.GetSymmSSSR(self.Mol))]
    @property
    def num_rotors(self):
        return NumRotatableBonds(self.Mol)
    
    @property
    def atom_dict(self):
        # check cache and generate dicts
        if self._atom_dict is None:
            self._dicts()
        return self._atom_dict
        
    @property
    def res_dict(self):
        # check cache and generate dicts
        if self._res_dict is None:
            self._dicts()
        return self._res_dict
        
    @property
    def ring_dict(self):
        # check cache and generate dicts
        if self._ring_dict is None:
            self._dicts()
        return self._ring_dict
    
    @property
    def clone(self):
        return Molecule(copy(self.Mol))
    
    def clone_coords(self, source):
        self.Mol.RemoveAllConformers()
        for conf in source.Mol.GetConformers():
            self.Mol.AddConformer(conf)
        return self
    
    def _dicts(self):
        # Atoms
        atom_dtype = [('id', 'int16'),
                 # atom info
                 ('coords', 'float16', 3),
                 ('charge', 'float16'),
                 ('atomicnum', 'int8'),
                 ('atomtype','a4'),
                 ('hybridization', 'int8'),
                 ('neighbors', 'float16', (4,3)), # non-H neighbors coordinates for angles (max of 6 neighbors should be enough)
                 # residue info
                 ('resid', 'int16'),
                 ('resname', 'a3'),
                 ('isbackbone', 'bool'),
                 # atom properties
                 ('isacceptor', 'bool'),
                 ('isdonor', 'bool'),
                 ('isdonorh', 'bool'),
                 ('ismetal', 'bool'),
                 ('ishydrophobe', 'bool'),
                 ('isaromatic', 'bool'),
                 ('isminus', 'bool'),
                 ('isplus', 'bool'),
                 ('ishalogen', 'bool'),
                 # secondary structure
                 ('isalpha', 'bool'),
                 ('isbeta', 'bool')
                 ]

        a = []
        atom_dict = np.empty(self.Mol.GetNumAtoms(), dtype=atom_dtype)
        metals = [3,4,11,12,13,19,20,21,22,23,24,25,26,27,28,29,30,31,37,38,39,40,41,42,43,44,45,46,47,48,49,50,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,87,88,89,90,91,
    92,93,94,95,96,97,98,99,100,101,102,103]
        for i, atom in enumerate(self.atoms):
            
            atomicnum = atom.atomicnum
            partialcharge = atom.partialcharge
            coords = atom.coords
            atomtype = atom.Atom.GetProp("_TriposAtomType") if atom.Atom.HasProp("_TriposAtomType") else atom.Atom.GetSymbol()
#            if self.protein:
#                residue = pybel.Residue(atom.OBAtom.GetResidue())
#            else:
#                residue = False
            residue = None
            
            # get neighbors, but only for those atoms which realy need them
            neighbors = np.empty(4, dtype=[('coords', 'float16', 3),('atomicnum', 'int8')])
            neighbors.fill(np.nan)
            for n, nbr_atom in enumerate(atom.neighbors):
                neighbors[n] = (nbr_atom.coords, nbr_atom.atomicnum)
            atom_dict[i] = (atom.idx,
                      coords,
                      partialcharge,
                      atomicnum,
                      atomtype,
                      np.clip(atom.Atom.GetHybridization()-1, 0, 3),
                      neighbors['coords'],
                      # residue info
                      residue.idx if residue else 0,
                      residue.name if residue else '',
                      False, #residue.OBResidue.GetAtomProperty(atom.OBAtom, 2) if residue else False, # is backbone
                      # atom properties
                      False, #atom.OBAtom.IsHbondAcceptor(),
                      False, #atom.OBAtom.IsHbondDonor(),
                      False, #atom.OBAtom.IsHbondDonorH(),
                      atomicnum in metals,
                      atomicnum == 6 and not (np.in1d(neighbors['atomicnum'], [6,1])).any(), #hydrophobe #doble negation, since nan gives False
                      atom.Atom.GetIsAromatic(),
                      atomtype in ['O3-', '02-' 'O-'], # is charged (minus)
                      atomtype in ['N3+', 'N2+', 'Ng+'], # is charged (plus)
                      atomicnum in [9,17,35,53], # is halogen?
                      False, # alpha
                      False # beta
                      )
        
#        if self.protein:
#            # Protein Residues (alpha helix and beta sheet)
#            res_dtype = [('id', 'int16'),
#                         ('resname', 'a3'),
#                         ('N', 'float16', 3),
#                         ('CA', 'float16', 3),
#                         ('C', 'float16', 3),
#                         ('isalpha', 'bool'),
#                         ('isbeta', 'bool')
#                         ] # N, CA, C

#            b = []
#            for residue in self.residues:<commit_msg>Reimplement RDKit features for interactions<commit_after>#-*. coding: utf-8 -*-
## Copyright (c) 2008-2011, <NAME>; 2012, <NAME>; 2014, <NAME>;
## All rights reserved.
##
##  This file is part of Cinfony.
##  The contents are covered by the terms of the BSD license
##  which is included in the file LICENSE_BSD.txt.

"""
rdkit - A Cinfony module for accessing the RDKit from CPython

Global variables:
  Chem and AllChem - the underlying RDKit Python bindings
  informats - a dictionary of supported input formats
  outformats - a dictionary of supported output formats
  descs - a list of supported descriptors
  fps - a list of supported fingerprint types
  forcefields - a list of supported forcefields
"""

import os
from copy import copy

import numpy as np

from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem import Descriptors
from rdkit import RDConfig

_descDict = dict(Descriptors.descList)

import rdkit.DataStructs
import rdkit.Chem.MACCSkeys
import rdkit.Chem.AtomPairs.Pairs
import rdkit.Chem.AtomPairs.Torsions
### ODDT ###
from rdkit.Chem.Lipinski import NumRotatableBonds
from rdkit.Chem.AllChem import ComputeGasteigerCharges
from rdkit.Chem.Pharm2D import Gobbi_Pharm2D,Generate

# PIL and Tkinter
try:
    import Tkinter as tk
    import Image as PIL
    import ImageTk as PILtk
except:
    PILtk = None

# Aggdraw
try:
    import aggdraw
    from rdkit.Chem.Draw import aggCanvas
except ImportError:
    aggdraw = None

fps = ['rdkit', 'layered', 'maccs', 'atompairs', 'torsions', 'morgan']
"""A list of supported fingerprint types"""
descs = _descDict.keys()
"""A list of supported descriptors"""

_formats = {'smi': "SMILES",
            'can': "Canonical SMILES",
            'mol': "MDL MOL file",
            'mol2': "Tripos MOL2 file",
            'sdf': "MDL SDF file",
            'inchi':"InChI",
            'inchikey':"InChIKey"}
_notinformats = ['can', 'inchikey']
_notoutformats = ['mol2']
if not Chem.INCHI_AVAILABLE:
    _notinformats += ['inchi']
    _notoutformats += ['inchi', 'inchikey']

informats = dict([(_x, _formats[_x]) for _x in _formats if _x not in _notinformats])
"""A dictionary of supported input formats"""
outformats = dict([(_x, _formats[_x]) for _x in _formats if _x not in _notoutformats])
"""A dictionary of supported output formats"""

base_feature_factory = AllChem.BuildFeatureFactory(os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef'))
""" Global feature factory based on BaseFeatures.fdef """

_forcefields = {'uff': AllChem.UFFOptimizeMolecule}
forcefields = _forcefields.keys()
"""A list of supported forcefields"""

def readfile(format, filename, *args, **kwargs):
    """Iterate over the molecules in a file.

    Required parameters:
       format - see the informats variable for a list of available
                input formats
       filename

    You can access the first molecule in a file using the next() method
    of the iterator:
        mol = readfile("smi", "myfile.smi").next()

    You can make a list of the molecules in a file using:
        mols = list(readfile("smi", "myfile.smi"))

    You can iterate over the molecules in a file as shown in the
    following code snippet:
    >>> atomtotal = 0
    >>> for mol in readfile("sdf", "head.sdf"):
    ...     atomtotal += len(mol.atoms)
    ...
    >>> print atomtotal
    43
    """
    if not os.path.isfile(filename):
        raise IOError, "No such file: '%s'" % filename
    format = format.lower()
    # Eagerly evaluate the supplier functions in order to report
    # errors in the format and errors in opening the file.
    # Then switch to an iterator...
    if format=="sdf":
        def filereader_sdf():
            block = ''
            n = 0
            f = open(filename)
            for line in f:
                block += line
                if line[:4] == '$$$$':
                    yield Molecule(source={'fmt': format, 'string': block, 'n': n, 'filename': filename})
                    n += 1
                    block = ''
            f.close()
        return filereader_sdf()
    elif format=="mol":
        def mol_reader():
            yield Molecule(Chem.MolFromMolFile(filename, *args, **kwargs))
        return mol_reader()
    elif format=="pdb":
        def mol_reader():
            yield Molecule(Chem.MolFromPDBFile(filename, *args, **kwargs))
        return mol_reader()
    elif format=="mol2":
        def filereader_mol2():
            block = ''
            data = ''
            n = 0
            f = open(filename)
            for line in f:
                if line[:1] == '#':
                    data += line
                elif line[:17] == '@<TRIPOS>MOLECULE':
                    if n>0: #skip `zero` molecule (any preciding comments and spaces)
                        yield Molecule(source={'fmt': format, 'string': block, 'n': n, 'filename': filename})
                    n += 1
                    block = data
                    data = ''
                block += line
            f.close()
            
        return filereader_mol2()
    
    elif format=="smi":
        iterator = Chem.SmilesMolSupplier(filename, delimiter=" \t",
                                          titleLine=False, *args, **kwargs)
        def smi_reader():
            for mol in iterator:
                yield Molecule(mol)
        return smi_reader()
    elif format=='inchi' and Chem.INCHI_AVAILABLE:
        def inchi_reader():
            for line in open(filename, 'r'):
                mol = Chem.inchi.MolFromInchi(line.strip(), *args, **kwargs)
                yield Molecule(mol)
        return inchi_reader()
    else:
        raise ValueError, "%s is not a recognised RDKit format" % format

def readstring(format, string):
    """Read in a molecule from a string.

    Required parameters:
       format - see the informats variable for a list of available
                input formats
       string

    Example:
    >>> input = "C1=CC=CS1"
    >>> mymol = readstring("smi", input)
    >>> len(mymol.atoms)
    5
    """
    format = format.lower()
    if format=="mol" or format=="sdf":
        mol = Chem.MolFromMolBlock(string)
    elif format=="mol2":
        mol = Chem.MolFromMol2Block(string)
    elif format=="pdb":
        mol = Chem.MolFromPDBBlock(string)
    elif format=="smi":
        mol = Chem.MolFromSmiles(string)
    elif format=='inchi' and Chem.INCHI_AVAILABLE:
        mol = Chem.inchi.MolFromInchi(string)
    else:
        raise ValueError,"%s is not a recognised RDKit format" % format
    if mol:
        return Molecule(mol)
    else:
        raise IOError, "Failed to convert '%s' to format '%s'" % (
            string, format)

class Outputfile(object):
    """Represent a file to which *output* is to be sent.

    Required parameters:
       format - see the outformats variable for a list of available
                output formats
       filename

    Optional parameters:
       overwite -- if the output file already exists, should it
                   be overwritten? (default is False)

    Methods:
       write(molecule)
       close()
    """
    def __init__(self, format, filename, overwrite=False):
        self.format = format
        self.filename = filename
        if not overwrite and os.path.isfile(self.filename):
            raise IOError, "%s already exists. Use 'overwrite=True' to overwrite it." % self.filename
        if format=="sdf":
            self._writer = Chem.SDWriter(self.filename)
        elif format=="smi":
            self._writer = Chem.SmilesWriter(self.filename, isomericSmiles=True)
        elif format in ('inchi', 'inchikey') and Chem.INCHI_AVAILABLE:
            self._writer= open(filename, 'w')
        else:
            raise ValueError,"%s is not a recognised RDKit format" % format
        self.total = 0 # The total number of molecules written to the file

    def write(self, molecule):
        """Write a molecule to the output file.

        Required parameters:
           molecule
        """
        if not self.filename:
            raise IOError, "Outputfile instance is closed."
        if self.format in ('inchi', 'inchikey'):
            self._writer.write(molecule.write(self.format) +'\n')
        else:
            self._writer.write(molecule.Mol)
        self.total += 1

    def close(self):
        """Close the Outputfile to further writing."""
        self.filename = None
        self._writer.flush()
        del self._writer

class Molecule(object):
    """Represent an rdkit Molecule.

    Required parameter:
       Mol -- an RDKit Mol or any type of cinfony Molecule

    Attributes:
       atoms, data, formula, molwt, title

    Methods:
       addh(), calcfp(), calcdesc(), draw(), localopt(), make3D(), removeh(),
       write()

    The underlying RDKit Mol can be accessed using the attribute:
       Mol
    """
    _cinfony = True

    def __init__(self, Mol = None, source = None, protein = False):
        if hasattr(Mol, "_cinfony"):
            a, b = Mol._exchange
            if a == 0:
                molecule = readstring("smi", b)
            else:
                molecule = readstring("mol", b)
            Mol = molecule.Mol

        self.Mol = Mol
        ### ODDT ###
        self.protein = protein
        self._atom_dict = None
        self._res_dict = None
        self._ring_dict = None
        self._coords = None
        self._charges = None
        # lazy
        self._source = source # dict with keys: n, fmt, string, filename
    
    # lazy Molecule parsing requires masked Mol
    @property
    def Mol(self):
        if not self._Mol and self._source:
            self._Mol = readstring(self._source['fmt'], self._source['string']).Mol
            self._source = None
        return self._Mol
        
    @Mol.setter
    def Mol(self, value):
        self._Mol = value
    
    @property
    def atoms(self): return [Atom(rdkatom) for rdkatom in self.Mol.GetAtoms()]
    @property
    def data(self): return MoleculeData(self.Mol)
    @property
    def molwt(self): return Descriptors.MolWt(self.Mol)
    @property
    def formula(self): return Descriptors.MolecularFormula(self.Mol)
    def _gettitle(self):
        # Note to self: maybe should implement the get() method for self.data
        if "_Name" in self.data:
            return self.data["_Name"]
        else:
            return ""
    def _settitle(self, val): self.Mol.SetProp("_Name", val)
    title = property(_gettitle, _settitle)
    @property
    def _exchange(self):
        if self.Mol.GetNumConformers() == 0:
            return (0, self.write("smi"))
        else:
            return (1, self.write("mol"))
    
    # cache frequently used properties and cache them in prefixed [_] variables
    @property
    def coords(self):
        if self._coords is None:
            self._coords = np.array([atom.coords for atom in self.atoms])
        return self._coords
    
    @property
    def charges(self):
        if self._charges is None:
            self._charges = np.array([atom.partialcharge for atom in self.atoms])
        return self._charges
    
    #### Custom ODDT properties ####
    @property
    def sssr(self):
        return [list(path) for path in list(Chem.GetSymmSSSR(self.Mol))]
    @property
    def num_rotors(self):
        return NumRotatableBonds(self.Mol)
    
    @property
    def atom_dict(self):
        # check cache and generate dicts
        if self._atom_dict is None:
            self._dicts()
        return self._atom_dict
        
    @property
    def res_dict(self):
        # check cache and generate dicts
        if self._res_dict is None:
            self._dicts()
        return self._res_dict
        
    @property
    def ring_dict(self):
        # check cache and generate dicts
        if self._ring_dict is None:
            self._dicts()
        return self._ring_dict
    
    @property
    def clone(self):
        return Molecule(copy(self.Mol))
    
    def clone_coords(self, source):
        self.Mol.RemoveAllConformers()
        for conf in source.Mol.GetConformers():
            self.Mol.AddConformer(conf)
        return self
    
    def _dicts(self):
        # Atoms
        atom_dtype = [('id', 'int16'),
                 # atom info
                 ('coords', 'float16', 3),
                 ('charge', 'float16'),
                 ('atomicnum', 'int8'),
                 ('atomtype','a4'),
                 ('hybridization', 'int8'),
                 ('neighbors', 'float16', (4,3)), # non-H neighbors coordinates for angles (max of 6 neighbors should be enough)
                 # residue info
                 ('resid', 'int16'),
                 ('resname', 'a3'),
                 ('isbackbone', 'bool'),
                 # atom properties
                 ('isacceptor', 'bool'),
                 ('isdonor', 'bool'),
                 ('isdonorh', 'bool'),
                 ('ismetal', 'bool'),
                 ('ishydrophobe', 'bool'),
                 ('isaromatic', 'bool'),
                 ('isminus', 'bool'),
                 ('isplus', 'bool'),
                 ('ishalogen', 'bool'),
                 # secondary structure
                 ('isalpha', 'bool'),
                 ('isbeta', 'bool')
                 ]

        a = []
        atom_dict = np.empty(self.Mol.GetNumAtoms(), dtype=atom_dtype)
        metals = [3,4,11,12,13,19,20,21,22,23,24,25,26,27,28,29,30,31,37,38,39,40,41,42,43,44,45,46,47,48,49,50,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,87,88,89,90,91,
    92,93,94,95,96,97,98,99,100,101,102,103]
        for i, atom in enumerate(self.atoms):
            
            atomicnum = atom.atomicnum
            partialcharge = atom.partialcharge
            coords = atom.coords
            atomtype = atom.Atom.GetProp("_TriposAtomType") if atom.Atom.HasProp("_TriposAtomType") else atom.Atom.GetSymbol()
#            if self.protein:
#                residue = pybel.Residue(atom.OBAtom.GetResidue())
#            else:
#                residue = False
            residue = None
            
            # get neighbors, but only for those atoms which realy need them
            neighbors = np.empty(4, dtype=[('coords', 'float16', 3),('atomicnum', 'int8')])
            neighbors.fill(np.nan)
            for n, nbr_atom in enumerate(atom.neighbors):
                neighbors[n] = (nbr_atom.coords, nbr_atom.atomicnum)
            atom_dict[i] = (atom.idx,
                      coords,
                      partialcharge,
                      atomicnum,
                      atomtype,
                      np.clip(atom.Atom.GetHybridization()-1, 0, 3),
                      neighbors['coords'],
                      # residue info
                      residue.idx if residue else 0,
                      residue.name if residue else '',
                      False, #residue.OBResidue.GetAtomProperty(atom.OBAtom, 2) if residue else False, # is backbone
                      # atom properties
                      False, #atom.OBAtom.IsHbondAcceptor(),
                      False, #atom.OBAtom.IsHbondDonor(),
                      False, #atom.OBAtom.IsHbondDonorH(),
                      atomicnum in metals,
                      False, #atomicnum == 6 and not (np.in1d(neighbors['atomicnum'], [6,1])).any(), #hydrophobe #doble negation, since nan gives False
                      atom.Atom.GetIsAromatic(),
                      atomtype in ['O3-', '02-' 'O-'], # is charged (minus)
                      atomtype in ['N3+', 'N2+', 'Ng+'], # is charged (plus)
                      atomicnum in [9,17,35,53], # is halogen?
                      False, # alpha
                      False # beta
                      )
        
        # Match features and mark them in atom_dict
        feats = base_feature_factory.GetFeaturesForMol(self.Mol)
        translate_feats = {'Donor':'isdonor',
                   'Acceptor':'isacceptor',
                   'NegIonizable':'isminus',
                   'PosIonizable':'isplus',
                   'Aromatic':'isaromatic',
                   'Hydrophobe':'ishydrophobe'
                   }
        feat_atom_ids = {}
        for f in feats:
            if f.GetFamily() in translate_feats:
                if not translate_feats[f.GetFamily()] in feat_atom_ids:
                    feat_atom_ids[translate_feats[f.GetFamily()]] = tuple()
                feat_atom_ids[translate_feats[f.GetFamily()]] += f.GetAtomIds()
        for key, aids in feat_atom_ids.items():
            atom_dict[key][np.array(aids)] = True
        
        ### FIX: remove acidic carbons from isminus group (they are part of smarts)
        atom_dict['isminus'][atom_dict['isminus'] & (atom_dict['atomicnum'] == 6)] = False
        
#        if self.protein:
#            # Protein Residues (alpha helix and beta sheet)
#            res_dtype = [('id', 'int16'),
#                         ('resname', 'a3'),
#                         ('N', 'float16', 3),
#                         ('CA', 'float16', 3),
#                         ('C', 'float16', 3),
#                         ('isalpha', 'bool'),
#                         ('isbeta', 'bool')
#                         ] # N, CA, C

#            b = []
#            for residue in self.residues:
{'fraction_non_alphanumeric': 0.11489936917993392, 'fraction_numerical': 0.015830579753679783, 'mean_word_length': 1.90421355666056, 'pattern_counts': {'":': 22, '<': 11, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3996756', 'n_tokens_mistral': 10325, 'n_tokens_neox': 9658, 'n_words': 2741}