implement generic inchitomol() with stereo support (#138)

* implement generic inchitomol() with stereo support

* fix stereo handling in inchitomol for double bonds + code styling

* add/modify tests

* more code cleaning

* inchitosdf: convert via MolGraph in order to support stereo

---------

Co-authored-by: hhaensel <[email protected]>

Author: hhaensel

src/inchi.jl

@@ -2,6 +2,38 @@
 # This file is a part of MolecularGraph.jl
 # Licensed under the MIT License http://opensource.org/licenses/MIT
 #
+const INCHI_StereoType_None       = 0
+const INCHI_StereoType_DoubleBond = 1
+const INCHI_StereoType_Tetrahedral= 2
+const INCHI_StereoType_Allene     = 3
+
+const INCHI_PARITY_NONE      = 0
+const INCHI_PARITY_ODD       = 1   # 'o'
+const INCHI_PARITY_EVEN      = 2   # 'e'
+const INCHI_PARITY_UNKNOWN   = 3   # 'u'
+const INCHI_PARITY_UNDEFINED = 4   # '?', rarely used
+
+const NO_ATOM = -1             # from header: NO_ATOM == -1
+
+const MAXVAL = 20
+const ATOM_EL_LEN = 6
+const NUM_H_ISOTOPES = 3
+const AT_NUM = Int16
+const S_CHAR = Int8
+struct inchi_Atom
+    x::Cdouble
+    y::Cdouble
+    z::Cdouble
+    neighbor::NTuple{MAXVAL,AT_NUM}
+    bond_type::NTuple{MAXVAL,S_CHAR}
+    bond_stereo::NTuple{MAXVAL,S_CHAR}
+    elname::NTuple{ATOM_EL_LEN,Cchar}
+    num_bonds::AT_NUM
+    num_iso_H::NTuple{NUM_H_ISOTOPES+1,S_CHAR}
+    isotopic_mass::AT_NUM
+    radical::S_CHAR
+    charge::S_CHAR
+end
 
 mutable struct inchi_Output
     szInChI::Cstring
@@ -18,10 +50,17 @@ mutable struct inchi_InputINCHI
                            # '/' or '-' depending on OS and compiler
 end
 
+struct inchi_Stereo0D
+    neighbor::NTuple{4, AT_NUM} # always 4 atoms
+    central_atom::AT_NUM        # central atom index (or NO_ATOM)
+    type::S_CHAR                # inchi_StereoType0D
+    parity::S_CHAR              # inchi_StereoParity0D (possibly combined)
+end
+
 mutable struct inchi_Input
     # the caller is responsible for the data allocation and deallocation
-    atom::Ptr{Cvoid}            # array of num_atoms elements
-    stereo0D::Ptr{Cvoid}        # array of num_stereo0D 0D stereo elements or NULL
+    atom::Ptr{inchi_Atom}            # array of num_atoms elements
+    stereo0D::Ptr{inchi_Stereo0D}        # array of num_stereo0D 0D stereo elements or NULL
     szOptions::Cstring          # InChI options: space-delimited; each is preceded by
                                 # '/' or '-' depending on OS and compiler
     num_atoms::Cshort           # number of atoms in the structure < MAX_ATOMS
@@ -30,20 +69,20 @@ end
 
 mutable struct inchi_InputEx
     # the caller is responsible for the data allocation and deallocation
-    atom::Ptr{Cvoid}            # array of num_atoms elements
-    stereo0D::Ptr{Cvoid}        # array of num_stereo0D 0D stereo elements or NULL
-    szOptions::Cstring          # InChI options: space-delimited; each is preceded by
-                                # '/' or '-' depending on OS and compiler
-    num_atoms::Cshort           # number of atoms in the structure < MAX_ATOMS
-    num_stereo0D::Cshort        # number of 0D stereo elements
+    atom::Ptr{inchi_Atom}          # array of num_atoms elements
+    stereo0D::Ptr{inchi_Stereo0D}  # array of num_stereo0D 0D stereo elements or NULL
+    szOptions::Cstring             # InChI options: space-delimited; each is preceded by
+                                   # '/' or '-' depending on OS and compiler
+    num_atoms::Cshort              # number of atoms in the structure < MAX_ATOMS
+    num_stereo0D::Cshort           # number of 0D stereo elements
     polymer::Ptr{Cvoid}
     v3000::Ptr{Cvoid}
 end
 
 mutable struct inchi_OutputStruct
     # the caller is responsible for the data allocation and deallocation
-    atom::Ptr{Cvoid}        # array of num_atoms elements
-    stereo0D::Ptr{Cvoid}    # array of num_stereo0D 0D stereo elements or NULL
+    atom::Ptr{inchi_Atom}        # array of num_atoms elements
+    stereo0D::Ptr{inchi_Stereo0D}    # array of num_stereo0D 0D stereo elements or NULL
     num_atoms::Cshort;      # number of atoms in the structure < MAX_ATOMS
     num_stereo0D::Cshort    # number of 0D stereo elements
     szMessage::Cstring      # Error/warning ASCIIZ message
@@ -62,13 +101,13 @@ end
 
 mutable struct inchi_OutputStructEx
     # the caller is responsible for the data allocation and deallocation
-    atom::Ptr{Cvoid}        # array of num_atoms elements
-    stereo0D::Ptr{Cvoid}    # array of num_stereo0D 0D stereo elements or NULL
-    num_atoms::Cshort;      # number of atoms in the structure < MAX_ATOMS
-    num_stereo0D::Cshort    # number of 0D stereo elements
-    szMessage::Cstring      # Error/warning ASCIIZ message
-    szLog::Cstring          # log-file ASCIIZ string, contains a human-readable list
-                            # of recognized options and possibly an Error/warn message
+    atom::Ptr{inchi_Atom}          # array of num_atoms elements
+    stereo0D::Ptr{inchi_Stereo0D}  # array of num_stereo0D 0D stereo elements or NULL
+    num_atoms::Cshort;             # number of atoms in the structure < MAX_ATOMS
+    num_stereo0D::Cshort           # number of 0D stereo elements
+    szMessage::Cstring             # Error/warning ASCIIZ message
+    szLog::Cstring                 # log-file ASCIIZ string, contains a human-readable list
+                                   # of recognized options and possibly an Error/warn message
     warningflags::NTuple{2, NTuple{2, Clong}} # warnings, see INCHIDIFF in inchicmp.h                        */
                             # [x][y]:
                             # x=0 => Reconnected if present in InChI
@@ -160,52 +199,235 @@ end
 inchikey(mol::SimpleMolGraph) = inchikey(inchi(mol))
 
 """
-    inchitosdf(inchi::String; options::String = "") -> Union{String,Nothing}
+    inchitosdf(inchi::AbstractString; options::String = "", config::Union{Nothing,Dict{Symbol,Any}})
 
-Generate sdf string from inchi string, `options` are specified in https://github.com/mojaie/libinchi/blob/master/INCHI_BASE/src/inchi_api.h
+Generate sdf string from inchi string.
+This new version goes via parsing to a MolGraph and exporting as sdf.
+Parsing `options` are specified in https://github.com/mojaie/libinchi/blob/master/INCHI_BASE/src/inchi_api.h
+Coordinates are generated via coordgen (Schrödinger coordgenlibs).
 """
-function inchitosdf(inchi::String; options::String = "", verbose::Bool = false)
-    # support the correct options format depending on OS
-    opts = opt_array(options)
-    # add a timeout of 60s per molecule
-    any(occursin.(r"^Wm?\d+$", opts)) || push!(opts, "W60")
-    # switch output to sdf format
-    "OutputSDF" ∈ opts || push!(opts, "OutputSDF")
-    options = opt_string(opts)
-
-    structure = inchi_OutputStructEx()
-    inchi_input = inchi_InputINCHI(Base.unsafe_convert(Cstring, inchi), Base.unsafe_convert(Cstring, options))
-    @ccall libinchi.GetStructFromINCHIEx(
-        inchi_input::Ref{inchi_InputINCHI}, structure::Ref{inchi_OutputStructEx})::Int32
-
-    input = inchi_InputEx(
-        structure.atom,
-        structure.stereo0D,
-        Base.unsafe_convert(Cstring, options),
-        structure.num_atoms,
-        structure.num_stereo0D,
-        structure.polymer,
-        structure.v3000
-    )
-    output = inchi_Output()
-    @ccall libinchi.GetINCHIEx(
-        input::Ref{inchi_InputEx}, output::Ref{inchi_Output})::Cint
-    report_output(output, verbose)
+function inchitosdf(inchi::AbstractString; options::String = "", config::Union{Nothing,Dict{Symbol,Any}} = nothing)
+    printv2mol(inchitomol(inchi; options, config))
+end
 
-    res = output.szInChI == C_NULL ? nothing : unsafe_string(output.szInChI)
+# decode parity byte (connected in low 3 bits, disconnected in bits 3..5)
+# returns the parity value to use (0..4)
+decode_parity(p::Integer) = begin
+    p_raw  = UInt8(p)                    # make sure unsigned before bit ops
+    p_conn = Int(p_raw & 0x07)           # low 3 bits
+    p_disc = Int((p_raw >> 3) & 0x07)    # next 3 bits
+    p_conn != 0 ? p_conn : p_disc        # prefer connected parity if present
+end
 
-    # Free buffers allocated by GetStructFromINCHIEx and GetINCHI
-    @ccall libinchi.FreeStructFromINCHIEx(structure::Ref{inchi_OutputStructEx})::Cvoid
-    @ccall libinchi.FreeINCHI(output::Ref{inchi_Output})::Cvoid
+function process_inchi_stereo!(g::T, structure::inchi_OutputStructEx) where T <: MolGraph
+    ET = edgetype(T)
+    
+    stereocenters = Dict{Int64, Tuple{Int64,Int64,Int64,Bool}}()   # center => (looking_from,v1,v2,is_clockwise)
+    stereobonds   = Dict{Edge{Int64}, Tuple{Int64,Int64,Bool}}()   # edge => (a1,a2,is_cis)
+
+    atoms = unsafe_wrap(Vector{inchi_Atom}, structure.atom, structure.num_atoms)
+    stereos = unsafe_wrap(Vector{inchi_Stereo0D}, structure.stereo0D, structure.num_stereo0D)
+    n_heavy = length(atoms)  # number of heavy atoms (without isotopic H/D/T)
+    is_heavy(idx) = 0 < idx <= n_heavy # for 1-based indices
+
+    for s in stereos
+        parity = decode_parity(s.parity)
+        parity == INCHI_PARITY_NONE || parity == INCHI_PARITY_UNKNOWN && continue  # nothing to set
+        parity == INCHI_PARITY_UNDEFINED && continue  # explicit undefined -> skip
+
+        # convert to 1-based indexing
+        nbrs = s.neighbor .+ 1 
+        center = s.central_atom + 1
+
+        if s.type == INCHI_StereoType_Tetrahedral
+            center == 0 && continue  # sanity check
+
+            # InChI neighbors order corresponds to W,X,Y,Z (see InChI comments).
+            w, x, y, z = nbrs   # these are already expanded indices (or NO_ATOM)
+            is_clockwise = (parity == INCHI_PARITY_EVEN)   # 'e' == clockwise per docs
+            stereocenters[center] = (w, x, y, is_clockwise)
+        elseif s.type == INCHI_StereoType_DoubleBond
+            all(is_heavy.(nbrs)) || continue  # skip disconnected or missing substituents
+
+            # the double bond is between the middle two atoms
+            # by InChI convention, neighbor[1,2] on one end, neighbor[3,4] on the other
+            # determine the actual bond by looking up the shared bond in the graph
+            bond_candidates = [u_edge(g, nbrs[i], nbrs[j]) for i in 1:2, j in 3:4]
+            # pick the one that actually exists in your bond list
+            edges = collect(g.eprops)  # materialize iterator
+            bond_key_index = findfirst(e -> e[1] in bond_candidates, edges)
+            bond_key_index === nothing && continue
+            bond_key = edges[bond_key_index][1]
+
+            is_trans = if parity == INCHI_PARITY_ODD
+                true      # InChI uses ODD -> trans
+            elseif parity == INCHI_PARITY_EVEN
+                false     # EVEN -> cis
+            else
+                # we currently don't get here, because we filtered out NONE/UNKNOWN/UNDEFINED above
+                # once we can handle undetermined stereo, we need to implement that logic here
+                continue   # unknown
+            end
+            
+            # store stereobond
+            stereobonds[bond_key] = (nbrs[1], nbrs[3], is_trans)
+        elseif s.type == INCHI_StereoType_Allene
+            @info "Allene stereo, not yet tested"
+            # Tried testing with `mol = smilestomol("C[C@]=C=C(C)F")`,
+            # but the InChI generated from it does not contain stereo information.
+            # Allene (axial) handling
+            # InChI gives neighbors as W,X,Y,Z where W/X are substituents on one end
+            # and Y/Z on the other. central_atom is the central cumulated C.
+            # We need to pick a viewing atom and two reference atoms to encode axial chirality.
+            center == 0 && continue
+
+            # For an allene to be chiral, both ends must have two substituents (i.e. W/X and Y/Z present)
+            all(is_heavy.(nbrs)) || continue
+            w, x, y, z = nbrs
+
+            is_clockwise = (parity == INCHI_PARITY_EVEN)
+
+            # store into same stereocenters dict so coordgen!/rendering code can use it
+            stereocenters[center] = (w, x, y, is_clockwise)
+        else
+            @warn "Unknown stereo type $(s.type), skipping"
+        end
+    end
 
-    return res
+    # merge into mol.gprops the same way you already did:
+    merge!(g.gprops.stereocenter, stereocenters)
+    merge!(g.gprops.stereobond, stereobonds)
+    return g
 end
 
 """
-    function inchitomol(inchi::String; options = "", verbose = false)
+    function inchitomol(inchi::String; options = "", config::Union{Nothing,Dict{Symbol,Any}} = nothing)
 
 Generate molecule from inchi string, `options` are specified in https://github.com/mojaie/libinchi/blob/master/INCHI_BASE/src/inchi_api.h
+`config` is for internal or advanced use only. Maybe removed in a future release.
 """
-function inchitomol(inchi::String; options = "", verbose = false)
-    inchitosdf(inchi; options, verbose) |> IOBuffer |> sdftomol
+function inchitomol(::Type{T}, inchi::AbstractString;
+    options::String = "",
+    config::Union{Nothing,Dict{Symbol,Any}} = nothing
+) where T <: MolGraph
+    inchi isa String || (inchi = "$inchi")
+
+    structure = inchi_OutputStructEx()
+    inchi_input = inchi_InputINCHI(Base.unsafe_convert(Cstring, inchi),
+                                   Base.unsafe_convert(Cstring, options))
+
+    ret = @ccall libinchi.GetStructFromINCHIEx(
+        inchi_input::Ref{inchi_InputINCHI},
+        structure::Ref{inchi_OutputStructEx}
+    )::Int32
+
+    if ret < 0
+        error("libinchi failed with code $ret")
+    end
+    
+    config = if config === nothing && T === SMILESMolGraph
+        Dict{Symbol, Any}(:on_init => smiles_on_init!, :on_update => smiles_on_update!)
+    elseif config === nothing && T === SDFMolGraph
+        Dict{Symbol, Any}(:on_init => sdf_on_init!, :on_update => sdf_on_update!)
+    else
+        config
+    end
+    
+    # determine the MolecularGraph atom/bond types and constructors
+    V = vproptype(T) === AbstractAtom ? SDFAtom : vproptype(T)
+    E = eproptype(T) === AbstractBond ? SDFBond : eproptype(T)
+    ET = edgetype(T)
+
+    atoms = unsafe_wrap(Vector{inchi_Atom}, structure.atom, structure.num_atoms)
+
+    # heavy atoms
+    expanded_syms = String[
+        String(reinterpret(UInt8, collect(Iterators.takewhile(!=(0x00), a.elname))))
+    for a in atoms]
+
+    # isotopic H/D/T
+    for a in atoms
+        # a.num_iso_H is a NTuple{4,S_CHAR}, where index 1 = non-specified, 2 = H, 3 = D, 4 = T
+        # let's add symbols only for explicit isotopes [1H], [D] and [T]
+        # as NUM_H_ISOTOPES = 3, we only need to check indices 2, 3, and 4
+        for iso_index in 1:NUM_H_ISOTOPES
+            n = a.num_iso_H[iso_index + 1]
+            iso_sym = ["H", "D", "T"][iso_index]
+            for _ in 1:n
+                push!(expanded_syms, iso_sym)
+            end
+        end
+    end
+
+    # build bonds to isotopic H/D/T
+    bonds = Tuple{Int,Int,Int}[]
+    first_isotope_index = length(atoms) + 1
+    next_iso_idx = first_isotope_index
+    for (parent_idx, a) in enumerate(atoms)
+        for iso_index in 1:NUM_H_ISOTOPES
+            n = a.num_iso_H[iso_index + 1]
+            for _ in 1:n
+                push!(bonds, (parent_idx, next_iso_idx, 1))
+                next_iso_idx += 1
+            end
+        end
+    end
+
+    # heavy atom bonds
+    for (i, a) in enumerate(atoms)
+        for k in 1:a.num_bonds
+            j = a.neighbor[k] + 1
+            j > i && push!(bonds, (i, j, Int(a.bond_type[k])))
+        end
+    end
+
+    N = length(expanded_syms)
+    vprops = V[]
+    edges = ET[]
+    eprops = E[]
+
+    for idx in 1:N
+        if idx <= length(atoms)
+            a = atoms[idx]
+            symbol = expanded_syms[idx]
+            d = Dict{String,Any}(
+                "symbol" => symbol,
+                "charge" => Int(a.charge),
+                "multiplicity" => (Int(a.radical) == 0 ? 1 : Int(a.radical)),
+                "isotope" => Int(a.isotopic_mass),
+            )
+            :coords in fieldnames(V) && push!(d, "coords" => [Float64(a.x), Float64(a.y), Float64(a.z)])
+            push!(vprops, V(d))
+        else
+            symbol = expanded_syms[idx]
+            d = Dict(
+                "symbol" => symbol,
+                "charge" => 0,
+                "multiplicity" => 1,
+                "isotope" => 0
+            )
+            :coords in fieldnames(V) && push!(d, "coords" => [0.0, 0.0, 0.0])
+            push!(vprops, V(d))
+        end
+    end
+
+    for (i, j, order) in bonds
+        ekey = i < j ? ET(i, j) : ET(j, i)
+        push!(edges, ekey)
+        push!(eprops, E(; order))
+    end
+
+    g = T(edges, vprops, eprops; NamedTuple((k, v) for (k, v) in config)...)
+
+    process_inchi_stereo!(g, structure)
+    
+    @ccall libinchi.FreeStructFromINCHIEx(structure::Ref{inchi_OutputStructEx})::Cvoid
+
+    coordgen!(g)
+    remove_hydrogens!(g)
+    return g
+end
+
+function inchitomol(inchi::String; options::String = "", config::Union{Nothing,Dict{Symbol,Any}} = nothing)
+    inchitomol(SDFMolGraph, inchi; options, config)
 end