timvdm

#include <iostream>

#include <vector>

#include <algorithm>

#include <cassert>

 

#include <Eigen/Core>

 

#include <openbabel/mol.h>

#include <openbabel/obconversion.h>

#include <openbabel/graphsym.h>

 

#include <graph.hh>

 

#include "permutation.h"

using namespace OpenBabel;

 

////////////////////////////////////////////////////////////////////////////////

//

//

//  O P E N B A B E L   S P E C I F I C   C O D E

//

//

////////////////////////////////////////////////////////////////////////////////

 

OpenBabel::OBMol* ReadMolecule(const std::string &filename)

{

  OpenBabel::OBConversion conv;

  OpenBabel::OBFormat *format = conv.FormatFromExt(filename.c_str());

  if (!format || !conv.SetInFormat(format)) {

    std::cout << "ERROR: Could not detect file format for: " << filename << std::endl;

    return 0;

  }

 

  std::ifstream ifs;

  ifs.open(filename.c_str());

  if (!ifs) {

    std::cout << "ERROR: Could not open file: " << filename << std::endl;

    return 0;  

  }

 

  OpenBabel::OBMol *obmol = new OpenBabel::OBMol;

  if (!conv.Read(obmol, &ifs)) {

    std::cout << "ERROR: Could not read molecule from file: " << filename << std::endl;

    delete obmol;

    ifs.close();

    return 0;

  }

 

  ifs.close();

  return obmol;

}

 

////////////////////////////////////////////////////////////////////////////////

//

//

//  U N I T   T E S T S

//

//

////////////////////////////////////////////////////////////////////////////////

 

 

 

void testPermutation()

{

  std::vector<unsigned int> e;

  e.push_back(1);

  e.push_back(2);

  e.push_back(3);

  e.push_back(4);

  Permutation p1(e), p2(e);

 

  assert(p1 == p2);

}

 

void testPermutationGroup()

{

  Permutation p;

  p.map.push_back(1);

  p.map.push_back(2);

  p.map.push_back(3);

  p.map.push_back(4);

 

  PermutationGroup group;

  assert( !group.contains(p) );

 

  group.add(p);

  assert( group.contains(p) );

}

 

void unitTests()

{

  testPermutation();

  testPermutationGroup();

 

}

 

 

 

 

void callback(void *param, unsigned int n, const unsigned int *aut)

{

  assert(param);

  /*

  fprintf((FILE*)param, "Generator: ");

  bliss::print_permutation((FILE*)param, n, aut, 0);

  fprintf((FILE*)param, "\n");

  */

 

  std::vector<unsigned int> elements;

  cout << "Generator: ";

  for (unsigned int i = 0; i < n; ++i) {

    cout << aut[i] + 1 << " ";

    elements.push_back(aut[i]+1);

  }

  cout << endl;

 

  PermutationGroup *generators = static_cast<PermutationGroup*>(param);

  generators->add(Permutation(elements));

}

 

void addInverses(PermutationGroup &G)

{

  unsigned int size = G.permutations.size();

  for (unsigned int i = 0; i < size; ++i) {

    Permutation inv_p(G.permutations.at(i).matrix().transpose());

    if (!G.contains(inv_p))

      G.add(inv_p);

  }

}

 

void addProducts(PermutationGroup &G)

{

  for (unsigned int i = 0; i < G.permutations.size(); ++i) {

    for (unsigned int j = 0; j < G.permutations.size(); ++j) {

      if (i >= j)

        continue;

 

      Permutation p = G.permutations.at(i) * G.permutations.at(j);

      if (!G.contains(p))

        G.add(p);

    }

  }

}

 

PermutationGroup findAutomorphisms(OpenBabel::OBMol *obmol, const std::vector<unsigned int> &symClasses)

{

  using OpenBabel::OBMolAtomIter;

  using OpenBabel::OBMolBondIter;

 

  // construct the bliss graph

  bliss::Graph g;

  std::map<OpenBabel::OBAtom*, unsigned int> atom2index;

  FOR_ATOMS_OF_MOL (atom, obmol) {

    atom2index[&*atom] = g.add_vertex(symClasses.at(atom->GetIndex()));

  }

  FOR_BONDS_OF_MOL (bond, obmol) {

    g.add_edge(atom2index[bond->GetBeginAtom()], atom2index[bond->GetEndAtom()]);

  }

 

  // use bliss to get the automorphism group generators

  PermutationGroup generators;

  bliss::Stats stats;

  g.find_automorphisms(stats, &callback, &generators);

 

  cout << "# Automorphisms:" << stats.group_size_approx << endl;

  unsigned long nAut = stats.group_size_approx;

 

  // construct the automorphism group

  PermutationGroup G;

 

  // add identity permutation

  std::vector<unsigned int> eElements;

  for (unsigned int i = 0; i < obmol->NumAtoms(); ++i)

    eElements.push_back(i+1);

  Permutation e(eElements);

  if (!G.contains(e))

    G.add(e);

 

  // add the generators

  for (unsigned int i = 0; i < generators.size(); ++i) {

    if (!G.contains(generators.at(i)))

      G.add(generators.at(i));

  }

 

 

  unsigned int counter = 0;

  unsigned int lastSize = 0;

  while (G.permutations.size() != lastSize) {

    lastSize = G.permutations.size();

   

    addInverses(G);

    addProducts(G);

   

    counter++;

    if (counter > 100)

      break;

  }

 

  return G;

}

 

 

 

 

 

 

int main(int argc, char **argv)

{

  unitTests();

 

  if (argc < 2) {

    std::cout << "Usage: " << argv[0] << " <filename>" << std::endl;

    return -1;

  }

 

  OpenBabel::OBMol *obmol = ReadMolecule(argv[1]);

  if (!obmol)

    return -1;

 

  std::vector<unsigned int> elements;

  for (int i = 0; i < obmol->NumAtoms(); ++i) {

    elements.push_back(i+1);

  }

 

  OpenBabel::OBGraphSym graphSym(obmol);

  std::vector<unsigned int> symClasses;

  graphSym.GetSymmetry(symClasses);

 

 

  PermutationGroup G = findAutomorphisms(obmol, symClasses);

 

  std::cout << "Automorphisms:" << std::endl;

  for (unsigned int i = 0; i < G.size(); ++i)

    G.at(i).print();

 

 

  return 0;

}