Changeset f5ea10

Timestamp:

Jun 22, 2018, 7:26:09 AM (7 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

Candidate_v1.6.1, ChemicalSpaceEvaluator

Children:

99c705

Parents:

8d56a6

git-author:

Frederik Heber <frederik.heber@…> (09/26/17 22:30:01)

git-committer:

Frederik Heber <frederik.heber@…> (06/22/18 07:26:09)

Message:

Added Graph6Reader, extended BoostGraphCreator, added ChemicalSpaceEvaluatorAction.

• added visible generateAllInducedSubgraphs to Extractors.
• TESTS: due to new option "graph6" containing a digit we needed to modify moltest_check.py to also scan for digits and not just letters.
• DOCU: Added evaluate-chemical-space to userguide.

Files:

• doc/userguide/userguide.xml (modified) (1 diff)
• src/Actions/GlobalListOfActions.hpp (modified) (1 diff)
• src/Actions/GraphAction/ChemicalSpaceEvaluatorAction.cpp (added)
• src/Actions/GraphAction/ChemicalSpaceEvaluatorAction.def (added)
• src/Actions/GraphAction/ChemicalSpaceEvaluatorAction.hpp (added)
• src/Actions/Makefile.am (modified) (3 diffs)
• src/FunctionApproximation/Extractors.cpp (modified) (7 diffs)
• src/FunctionApproximation/Extractors.hpp (modified) (2 diffs)
• src/Graph/BoostGraphCreator.cpp (modified) (2 diffs)
• src/Graph/BoostGraphCreator.hpp (modified) (1 diff)
• src/Graph/Graph6Reader.cpp (added)
• src/Graph/Graph6Reader.hpp (added)
• src/Graph/Makefile.am (modified) (2 diffs)
• tests/Python/AllActions/moltest_check.py (modified) (1 diff)
• tests/Python/AllActions/options.dat (modified) (1 diff)

doc/userguide/userguide.xml

@@ -1268 +1268 @@
           to both other atoms serves as rotation joint).</para>
         </section>
+        <section xml:id="bond.evaluate-chemical-space">
+          <title xml:id="bond.evaluate-chemical-space.title">Evaluate the Chemical Space</title>
+          <para>Imagine that we are given a graph consisting of nodes and edges.
+          As we have been speaking extensively of the adjacency graph before, then
+          any graph where nodes are assigned a chemical element and edges are given
+          a degree, effectively represent a molecule of covalently bonded atoms.</para>
+          <para>This is the notion behind constructing parts of the chemical space that
+          encompasses all molecules that either already exist or could be devised,
+          i.e. that are stable.</para>
+          <para>One method of creating points in this chemical space, i.e. stable
+          molecules, is presented in the article by Hamaekers et al. 2017, where
+          they follow the well-known <emphasis>octet rule</emphasis>.</para>
+          <para>In this action essentially the same method is implemented. Given
+          an arbitrary graph, encoded as Graph6 string, such as produced by
+          the <link xlink:href="http://users.cecs.anu.edu.au/~bdm/nauty/">nauty</link> 
+          toolset (Brendan McKay and Adolfo Piperno), and a set of chemical
+          elements, the action creates every possible molecule using the graph
+          and the given elements for each of its nodes, changing bond degrees
+          as long as to fulfil the octet rule.</para>
+          <programlisting>--evaluate-chemical-space \
+  --graph6 "B`" --elements C C</programlisting>
+        <para>The graph6 string "B`" represents the simplest graph consisting
+        of just two nodes, connected by a single edge. Designating each node
+        to be a carbon atom, this action will then produce the carbonhydrates
+        C2H6, C2H4, and C2H2.</para>
+        <para>At the moment, the molecules are not created but their signatures
+        are used in looking up the fragments in the internal homology container,
+        see <link linkend="fragmentation">Fragmentation </link> and 
+        <link linkend="homology">Homologies</link>. Instead of creating atomic
+        coordinates, approximate energies are given that allow to estimate
+        whether the molecule candidate is stable or not. Note that this requires
+        that a suitable homology container file has been loaded first that contains
+        all necessary fragments.</para>
+        <note>The graph is automatically saturated with hydrogens in obeying
+        the octet rule, i.e. hydrogen should not be given in the list of elements.</note>
+        <para>Using nauty that produces all possible graphs with a fixed number
+        of nodes, these can be simply fed into <productname>MoleCuilder </productname>
+        to produce all associated molecules.</para>
+        </section>
       </section>
       <section xml:id="molecule">

src/Actions/GlobalListOfActions.hpp

@@ -75 +75 @@
   (GeometryPositionToVector) \
   (GeometryRemove) \
+  (GraphChemicalSpaceEvaluator) \
   (GraphUpdateMolecules) \
   (GraphCorrectBondDegree) \

src/Actions/Makefile.am

@@ -312 +312 @@
 
 GRAPHACTIONSOURCE = \
+  Actions/GraphAction/ChemicalSpaceEvaluatorAction.cpp \
   Actions/GraphAction/CorrectBondDegreeAction.cpp \
   Actions/GraphAction/CreateAdjacencyAction.cpp \
@@ -319 +320 @@
   Actions/GraphAction/UpdateMoleculesAction.cpp
 GRAPHACTIONHEADER = \
+  Actions/GraphAction/ChemicalSpaceEvaluatorAction.hpp \
   Actions/GraphAction/CorrectBondDegreeAction.hpp \
   Actions/GraphAction/CreateAdjacencyAction.hpp \
@@ -326 +328 @@
   Actions/GraphAction/UpdateMoleculesAction.hpp
 GRAPHACTIONDEFS = \
+  Actions/GraphAction/ChemicalSpaceEvaluatorAction.def \
   Actions/GraphAction/CorrectBondDegreeAction.def \
   Actions/GraphAction/CreateAdjacencyAction.def \

src/FunctionApproximation/Extractors.cpp

@@ -34 +34 @@
 #include <config.h>
 #endif
-
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/breadth_first_search.hpp>
-#include <boost/graph/subgraph.hpp>
 
 //#include "CodePatterns/MemDebug.hpp"
@@ -48 +44 @@
 #include <vector>
 #include <boost/assign.hpp>
-#include <boost/bimap.hpp>
-#include <boost/bimap/set_of.hpp>
-#include <boost/bimap/multiset_of.hpp>
 #include <boost/bind.hpp>
 #include <boost/foreach.hpp>
@@ -67 +60 @@
 
 using namespace boost::assign;
+
+using namespace Extractors;
 
 FunctionModel::arguments_t
@@ -186 +181 @@
 }
 
-typedef size_t level_t;
-typedef size_t node_t;
-typedef std::multimap< level_t, node_t > nodes_per_level_t;
-typedef std::set<node_t> nodes_t;
-typedef std::set<nodes_t> set_of_nodes_t;
-typedef boost::property_map < boost::adjacency_list <>, boost::vertex_index_t >::type index_map_t;
-
-typedef boost::bimap<
-    boost::bimaps::set_of< size_t >,
-    boost::bimaps::multiset_of< Extractors::ParticleType_t >
-> type_index_lookup_t;
-
-typedef std::set<node_t> set_type;
-typedef std::set<set_type> powerset_type;
-
-typedef boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS,
-    boost::no_property, boost::no_property > UndirectedGraph;
-typedef boost::subgraph< UndirectedGraph > UndirectedSubgraph;
-
-typedef std::map< node_t, std::pair<Extractors::ParticleType_t, size_t> > node_FragmentNode_map_t;
-
-typedef std::map< argument_t::indices_t, size_t> argument_placement_map_t;
-
-typedef std::map<size_t, size_t> argindex_to_nodeindex_t;
-
 void insertIntoNodeFragmentMap(
     node_FragmentNode_map_t &_node_FragmentNode_map,
@@ -334 +304 @@
  * \param index_map with indices per \a graph' vertex
  */
-static void generateAllInducedConnectedSubgraphs(
+void Extractors::generateAllInducedConnectedSubgraphs(
     const size_t N,
     const level_t level,
@@ -438 +408 @@
 }
 
-static HomologyGraph createHomologyGraphFromNodes(
+HomologyGraph Extractors::createHomologyGraphFromNodes(
     const nodes_t &nodes,
     const type_index_lookup_t &type_index_lookup,
@@ -501 +471 @@
 
   return HomologyGraph(graph_nodes, graph_edges);
-}
-
-/**
- * I have no idea why this is so complicated with BGL ...
- *
- * This is taken from the book "The Boost Graph Library: User Guide and Reference Manual, Portable Documents",
- * chapter "Basic Graph Algorithms", example on calculating the bacon number.
- */
-template <typename DistanceMap>
-class distance_recorder : public boost::default_bfs_visitor
-{
-public:
-  distance_recorder(DistanceMap dist) : d(dist) {}
-
-  template <typename Edge, typename Graph>
-  void tree_edge(Edge e, const Graph &g) const {
-    typename boost::graph_traits<Graph>::vertex_descriptor u = source(e,g), v = target(e,g);
-    d[v] = d[u] + 1;
-  }
-
-private:
-  DistanceMap d;
-};
-
-template <typename DistanceMap>
-distance_recorder<DistanceMap> record_distance(DistanceMap d)
-{
-  return distance_recorder<DistanceMap>(d);
 }
 

src/FunctionApproximation/Extractors.hpp

@@ -14 +14 @@
 #endif
 
+#include <boost/bimap.hpp>
+#include <boost/bimap/set_of.hpp>
+#include <boost/bimap/multiset_of.hpp>
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/breadth_first_search.hpp>
+#include <boost/graph/subgraph.hpp>
 #include <boost/function.hpp>
+
+#include <map>
+#include <set>
 
 #include "Fragmentation/EdgesPerFragment.hpp"
@@ -51 +60 @@
   typedef std::vector<ParticleType_t> ParticleTypes_t;
 
+  typedef size_t level_t;
+  typedef size_t node_t;
+  typedef std::multimap< level_t, node_t > nodes_per_level_t;
+  typedef std::set<node_t> nodes_t;
+  typedef std::set<nodes_t> set_of_nodes_t;
+
+  typedef boost::bimap<
+      boost::bimaps::set_of< size_t >,
+      boost::bimaps::multiset_of< Extractors::ParticleType_t >
+  > type_index_lookup_t;
+
+  typedef std::set<node_t> set_type;
+  typedef std::set<set_type> powerset_type;
+
+  typedef boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS,
+      boost::property<boost::vertex_name_t, atomId_t>,
+      boost::property<boost::vertex_color_t, boost::default_color_type> /* needed for limited-depth DFS,
+      otherwise the property_map gets full size of graph */
+      > UndirectedGraph;
+  typedef boost::subgraph< UndirectedGraph > UndirectedSubgraph;
+
+  typedef boost::property_map < UndirectedGraph, boost::vertex_index_t >::type index_map_t;
+
+  typedef std::map< node_t, std::pair<Extractors::ParticleType_t, size_t> > node_FragmentNode_map_t;
+
+  typedef std::map< argument_t::indices_t, size_t> argument_placement_map_t;
+
+  typedef std::map<size_t, size_t> argindex_to_nodeindex_t;
+
+  /**
+   * I have no idea why this is so complicated with BGL ...
+   *
+   * This is taken from the book "The Boost Graph Library: User Guide and Reference Manual, Portable Documents",
+   * chapter "Basic Graph Algorithms", example on calculating the bacon number.
+   */
+  template <typename DistanceMap>
+  class distance_recorder : public boost::default_bfs_visitor
+  {
+  public:
+    distance_recorder(DistanceMap dist) : d(dist) {}
+
+    template <typename Edge, typename Graph>
+    void tree_edge(Edge e, const Graph &g) const {
+      typename boost::graph_traits<Graph>::vertex_descriptor u = source(e,g), v = target(e,g);
+      d[v] = d[u] + 1;
+    }
+
+  private:
+    DistanceMap d;
+  };
+
+  template <typename DistanceMap>
+  distance_recorder<DistanceMap> record_distance(DistanceMap d)
+  {
+    return distance_recorder<DistanceMap>(d);
+  }
+
+  HomologyGraph createHomologyGraphFromNodes(
+      const nodes_t &nodes,
+      const type_index_lookup_t &type_index_lookup,
+      const UndirectedGraph &graph,
+      const index_map_t &index_map
+      );
+
+  void generateAllInducedConnectedSubgraphs(
+      const size_t N,
+      const level_t level,
+      const nodes_t &nodes,
+      set_of_nodes_t &set_of_nodes,
+      const nodes_per_level_t &nodes_per_level,
+      const UndirectedGraph &graph,
+      const std::vector<size_t> &_distance,
+      const index_map_t &index_map);
+
   /** Namespace for some internal helper functions.
    *

src/Graph/BoostGraphCreator.cpp

@@ -43 +43 @@
 
 #include "Atom/atom.hpp"
+#include "Graph/Graph6Reader.hpp"
 #include "molecule.hpp"
 
@@ -94 +95 @@
 }
 
+void BoostGraphCreator::createFromGraph6String(
+    const std::string &_graph_string)
+{
+  Graph6Reader reader;
+  {
+    std::stringstream inputstream(_graph_string);
+    reader(inputstream);
+  }
+
+  graph = UndirectedGraph();
+
+  // add nodes
+  for(int numbers = 0; numbers < reader.get_num_nodes(); ++numbers) {
+    const atomId_t atomid = numbers;
+    Vertex v = boost::add_vertex(atomid, graph);
+    const atomId_t vertexname = boost::get(boost::get(boost::vertex_name, graph), v);
+    const nodeId_t vertexindex = boost::get(boost::get(boost::vertex_index, graph), v);
+    LOG(2, "DEBUG: Adding node " << vertexindex << " associated to atom #" << vertexname);
+    ASSERT( vertexname == atomid,
+        "BoostGraphCreator::createFromRange() - atomid "+toString(atomid)
+        +" is not name of vertex "+toString(vertexname)+".");
+    atomids_nodeids.insert( std::make_pair(vertexname, vertexindex) );
+  }
+
+  // add edges
+  const Graph6Reader::edges_t &edges = reader.get_edges();
+  for(Graph6Reader::edges_t::const_iterator iter = edges.begin();
+      iter != edges.end(); ++iter) {
+    // graph6 contains only upper triangle of adjacency matrix, hence add only once
+    const nodeId_t leftnodeid = iter->first;
+    const nodeId_t rightnodeid = iter->second;
+    boost::add_edge(leftnodeid, rightnodeid, graph);
+  }
+}
+
 BoostGraphCreator::nodeId_t BoostGraphCreator::getNodeId(
     const atomId_t &_atomid) const

src/Graph/BoostGraphCreator.hpp

@@ -90 +90 @@
       const predicate_t &_pred);
 
+  /** Creates the boost::graph from the defining graph6 string where the atom
+   * nodes map is simply the identity.
+   *
+   * \param _graph_string graph6 string defining the graph
+   */
+  void createFromGraph6String(
+      const std::string &_graph_string);
+
   /** Getter for the created graph.
    *

src/Graph/Makefile.am

@@ -10 +10 @@
         Graph/ConnectedSubgraph.cpp \
         Graph/CyclicStructureAnalysis.cpp \
-        Graph/DepthFirstSearchAnalysis.cpp
-                                  
+        Graph/DepthFirstSearchAnalysis.cpp \
+        Graph/Graph6Reader.cpp
+
 GRAPHHEADER = \
         Graph/AdjacencyList.hpp \
@@ -23 +24 @@
         Graph/CyclicStructureAnalysis.hpp \
         Graph/DepthFirstSearchAnalysis.hpp \
+        Graph/Graph6Reader.hpp \
         Graph/ListOfLocalAtoms.hpp
 

tests/Python/AllActions/moltest_check.py

@@ -24 +24 @@
 def CheckParameters(docstring):
     result = 0
-    params = re.findall(r'\(str\)([-_a-zA-Z]*)', docstring)
+    params = re.findall(r'\(str\)([-_a-zA-Z0-9]*)', docstring)
 
     for param in params:

tests/Python/AllActions/options.dat

@@ -92 +92 @@
 fragment_path   "test/"
 fragment_prefix "BondFragment"
+graph6  "B`"
 grid_level      "5"
 help    "help"