// -*- C++ -*-

/* 
 * GChemPaint library
 * molecule.cc 
 *
 * Copyright (C) 2001-2007 Jean Bréfort <jean.brefort@normalesup.org>
 *
 * This program is free software; you can redistribute it and/or 
 * modify it under the terms of the GNU General Public License as 
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
 * USA
 */

#include "gchempaint-config.h"
#include "molecule.h"
#include "document.h"
#include "view.h"
#include "application.h"
#include "tool.h"
#include "stringdlg.h"
#include <glib/gi18n-lib.h>
#include <unistd.h>
#include <openbabel/obconversion.h>
#include <cmath>
#include <clocale>
#include <cstring>

static void do_export_to_ghemical (gcpMolecule* pMol)
{
	pMol->ExportToGhemical ();
}

static void do_select_alignment (GObject *action, gcpMolecule* pMol)
{
	Object *object = (Object*) g_object_get_data (action, "item");
	pMol->SelectAlignmentItem (object);
}

static void do_build_inchi (gcpMolecule* pMol)
{
	pMol->ShowInChI ();
}

static void do_build_smiles (gcpMolecule* pMol)
{
	pMol->BuildSmiles ();
}

static void do_show_webbook (gcpMolecule* pMol)
{
	pMol->ShowWebBase ("http://webbook.nist.gov/cgi/cbook.cgi?Name=", "&Units=SI");
}

static void do_show_pubchem (gcpMolecule* pMol)
{
	pMol->ShowWebBase ("http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?CMD=search&DB=pccompound&term=\"", "\"");
}

static void do_open_in_calc (gcpMolecule* pMol)
{
	pMol->OpenCalc ();
}

gcpMolecule::gcpMolecule(TypeId Type): Object(Type)
{
	m_Alignment = NULL;
	m_Changed = true;
}

gcpMolecule::gcpMolecule(gcpAtom* pAtom): Object(MoleculeType)
{
	AddAtom(pAtom);
	gcpChain* pChain = new gcpChain(this, pAtom); //will find the cycles
	delete pChain;
	m_Alignment = NULL;
	m_Changed = true;
}

gcpMolecule::~gcpMolecule()
{
	std::list<gcpBond*>::iterator n;
	for (n = m_Bonds.begin(); n != m_Bonds.end(); n++) (*n)->RemoveAllCycles();
	while (!m_Cycles.empty())
	{
		delete m_Cycles.front();
		m_Cycles.pop_front();
	}
	while (!m_Chains.empty())
	{
		delete m_Chains.front();
		m_Chains.pop_front();
	}
}

void gcpMolecule::AddAtom (gcpAtom* pAtom)
{
	m_Atoms.remove (pAtom); // avoid duplicates
	m_Atoms.push_back (pAtom);
	AddChild (pAtom);
}

void gcpMolecule::AddFragment(gcpFragment* pFragment)
{
	m_Fragments.push_back(pFragment);
	AddChild(pFragment);
}

void gcpMolecule::AddBond(gcpBond* pBond)
{
	if (pBond->GetAtom (0) && pBond->GetAtom (1))
		CheckCrossings (pBond);
	m_Bonds.push_back(pBond);
	AddChild(pBond);
	EmitSignal (OnChangedSignal);
}

void gcpMolecule::Remove(Object* pObject)
{
	if (pObject == m_Alignment)
		m_Alignment = NULL;
	switch (pObject->GetType())
	{
		case AtomType:
			m_Atoms.remove((gcpAtom*)pObject);
			break;
		case BondType:
			m_Bonds.remove((gcpBond*)pObject);
			break;
		case FragmentType:
			m_Fragments.remove((gcpFragment*)pObject);
			break;
	}
	pObject->SetParent(GetParent());
}

void gcpMolecule::UpdateCycles(gcpBond* pBond)	//FIXME: function not totally implemented
{
	gcpChain* pChain = new gcpChain(this, pBond); //will find the cycles
	delete pChain;
}

typedef struct {
	gcpAtom *pAtom;
	int sb; // number of shared bonds
	int so0, so1; // total bond order of shared bonds for each atom
} MergedAtom;

bool gcpMolecule::Merge(gcpMolecule* pMolecule, bool RemoveDuplicates)
{
	gcpAtom* pAtom;
	gcpFragment* pFragment;
	gcpBond* pBond;
	gcpChain* pChain;
	gcpCycle* pCycle;
	if (RemoveDuplicates) {
		std::list<gcpAtom*>::iterator i = m_Atoms.begin (), end = m_Atoms.end ();
		double x, y, x0, y0, x1, y1;
		MergedAtom *ma;
		gcpBond *b0, *b1;
		bool DoMerge;
		int n;
		map<gcpAtom*, MergedAtom*> AtomMap;
		map<gcpAtom*, MergedAtom*>::iterator j, endj; 
		map<gcpBond*, gcpBond*> BondMap;
		map<gcpBond*, gcpBond*>::iterator b, endb; 
		for (; i != end; i++) {
			(*i)->GetCoords (&x, &y);
			pAtom = (gcpAtom*) pMolecule->GetAtomAt (x, y);
			if (pAtom) {
				if ((*i)->GetZ()!= pAtom->GetZ())
				{
					// Cannot merge atoms which are not of the same element.
					endj = AtomMap.end ();
					for (j = AtomMap.begin (); j != endj; j++)
						delete (*j).second;
					return false;
				}
				ma = new MergedAtom;
				ma->sb = ma->so0 = ma->so1 = 0;
				ma->pAtom = pAtom;
				endj = AtomMap.end ();
				for (j = AtomMap.begin (); j != endj; j++) {
					if ((b1 = (gcpBond*) pAtom->GetBond ((*j).second->pAtom))) {
						b0 = (gcpBond*) (*i)->GetBond ((*j).first);
						if (b0) {
							ma->sb++;
							ma->so0 += b0->GetOrder ();
							ma->so1 += b1->GetOrder ();
							(*j).second->sb++;
							(*j).second->so0 += b0->GetOrder ();
							(*j).second->so1 += b1->GetOrder ();
							BondMap[b0] = b1;
						}
					}
				}
				AtomMap[*i] = ma;
			}
		}
		// Now check if merging is possible for each shared atom.
		DoMerge = AtomMap.size() != 0;
		if (DoMerge){
			x = y = 0.;
			endj = AtomMap.end ();
			for (j = AtomMap.begin (); j != endj; j++) {
				ma = (*j).second;
				n = ma->pAtom->GetTotalBondsNumber () - ma->so0 - ma->so1 + ma->sb;
				if (!(*j).first->AcceptNewBonds (n)) {
					DoMerge = false;
					break;
				}
				n = (*j).first->GetTotalBondsNumber () - ma->so0 - ma->so1 + ma->sb;
				if (!ma->pAtom->AcceptNewBonds (n)) {
					DoMerge = false;
					break;
				}
				(*j).first->GetCoords (&x0, &y0);
				ma->pAtom->GetCoords (&x1, &y1);
				x += x1 - x0;
				y += y1 - y0;
			}
		}
		if (DoMerge) {
			//First align molecules
			x /= 2.* AtomMap.size();
			y /= 2.* AtomMap.size();
			Move (x, y);
			pMolecule->Move (-x, -y);

			//Then align each atom individually
			endj = AtomMap.end ();
			for (j = AtomMap.begin (); j != endj; j++) {
				(*j).first->GetCoords (&x0, &y0);
				(*j).second->pAtom->GetCoords (&x1, &y1);
				(*j).first->Move ((x1 - x0) / 2.,(y1 - y0) / 2.); 
			}
			gcpView *pView = ((gcpDocument*) GetDocument ())->GetView ();

			/* Treat shared bonds (set order to 1, store max order in b1 and remove the bond
			 * from pMolecule. */
			endb = BondMap.end ();
			for (b = BondMap.begin (); b != endb; b++) {
				b1 = (*b).second;
				n = (*b).first->GetOrder ();
				pView->Remove (b1);
				pMolecule->Remove (b1);
				(*b).first->SetOrder (1);
				pAtom = (gcpAtom*) b1->GetAtom (0);
				pAtom->RemoveBond(b1);
				b1->ReplaceAtom (pAtom, NULL);
				pAtom = (gcpAtom*) b1->GetAtom (1);
				pAtom->RemoveBond(b1);
				b1->ReplaceAtom (pAtom, NULL);
				if (n > b1->GetOrder ())
					b1->SetOrder (n);
			}

			// Treat shared atoms and delete from pMolecule
			map< Atom *, Bond * >::iterator ai;  	 
			endj = AtomMap.end ();
			for (j = AtomMap.begin (); j != endj; j++) {
				b0 = (gcpBond*) (*j).second->pAtom->GetFirstBond (ai);
				while (b0) {
					b0->ReplaceAtom ((*j).second->pAtom, (*j).first);
					(*j).first->AddBond (b0);
					b0 = (gcpBond*) (*j).second->pAtom->GetNextBond (ai);
				}
				pMolecule->Remove ((*j).second->pAtom);
				pView->Remove ((*j).second->pAtom);
				delete (*j).second->pAtom;
			}			

			// Try to restore max bond order for shared bonds and destroy old bonds
			endb = BondMap.end ();
			for (b = BondMap.begin (); b != endb; b++) {
				n = (*b).second->GetOrder () - 1;
				b0 = (*b).first;
				while ((n > 0) &&
						(!((gcpAtom*) b0->GetAtom(0))->AcceptNewBonds (n) ||
						!((gcpAtom*) b0->GetAtom(1))->AcceptNewBonds (n)))
					n--;
				if (n > 0)
					b0->SetOrder (n + 1);
				delete (*b).second;
			}
		}

		// Clean memory
		endj = AtomMap.end ();
		for (j = AtomMap.begin (); j != endj; j++)
			delete (*j).second;
		//return if merging is not possible
		if (!DoMerge) return false;
	}
	while (!pMolecule->m_Atoms.empty())
	{
		pAtom = pMolecule->m_Atoms.front();
		AddAtom(pAtom);
		pMolecule->m_Atoms.pop_front();
	}
	while (!pMolecule->m_Fragments.empty())
	{
		pFragment = pMolecule->m_Fragments.front();
		AddFragment(pFragment);
		pMolecule->m_Fragments.pop_front();
	}
	while (!pMolecule->m_Bonds.empty())
	{
		pBond = pMolecule->m_Bonds.front();
		AddBond(pBond);
		pMolecule->m_Bonds.pop_front();
	}
	while (!pMolecule->m_Chains.empty())	//FIXME: Chains should change
	{
		pChain = pMolecule->m_Chains.front();
		m_Chains.push_back(pChain);
		pMolecule->m_Chains.pop_front();
	}
	while (!pMolecule->m_Cycles.empty())
	{
		pCycle = pMolecule->m_Cycles.front();
		m_Cycles.push_back(pCycle);
		pMolecule->m_Cycles.pop_front();
	}
	Object* pObj = pMolecule->GetParent ();
	delete pMolecule;
	pObj->EmitSignal (OnChangedSignal);
	if (RemoveDuplicates)
		UpdateCycles ();
	EmitSignal (OnChangedSignal);
	return true;
}

bool gcpMolecule::Load (xmlNodePtr node)
{
	char* buf;
	xmlNodePtr child;
	Object* pObject;
	gcpDocument* pDoc = (gcpDocument*) GetDocument ();

	buf = (char*) xmlGetProp (node, (xmlChar*) "id");
	if (buf) {
		SetId (buf);
		xmlFree (buf);
	}
	child = GetNodeByName (node, (char*) "atom");
	while (child) {
		pObject = new gcpAtom ();
		if (pDoc)
			AddChild (pObject);
		if (!pObject->Load (child)) {
			delete pObject;
			return false;
		}
		if (pDoc)
			pDoc->AddAtom ((gcpAtom*) pObject);
		AddAtom ((gcpAtom*) pObject);
		child = GetNextNodeByName (child->next, (char*) "atom");
	}
	
	child = GetNodeByName (node, (char*) "fragment");
	while (child) {
		pObject = new gcpFragment ();
		AddChild (pObject);
		if (!pObject->Load (child))  {
			delete pObject;
			return false;
		}
		if (pDoc)
			pDoc->AddFragment ((gcpFragment*) pObject);
		AddFragment ((gcpFragment*) pObject);
		child = GetNextNodeByName (child->next, (char*) "fragment");
	}

	child = GetNodeByName (node, (char*) "bond");
	while (child) {
		pObject = new gcpBond ();
		AddBond ((gcpBond*) pObject);
		if (!pObject->Load (child)) {
			delete pObject;
			m_Bonds.remove ((gcpBond*) pObject);
			return false;
		}
		if (pDoc)
			pDoc->AddBond ((gcpBond*) pObject);
		child = GetNextNodeByName (child->next, (char*) "bond");
		CheckCrossings ((gcpBond*) pObject);
	}
	if (!m_Atoms.empty ()) {
		gcpAtom* pAtom = m_Atoms.front ();
		std::list<gcpAtom*>::iterator i = m_Atoms.begin ();
		i++;
		for (; i != m_Atoms.end (); i++)
			(*i)->SetParent (NULL);
		gcpChain* pChain = new gcpChain (this, pAtom); //will find the cycles
		delete pChain;
	}
	buf = (char*) xmlGetProp (node, (const xmlChar*) "valign");
	if (buf) {
		m_Alignment = GetDescendant (buf);
		xmlFree (buf);
		if (!m_Alignment)
			return false;
	}
	m_Changed = true;
	return true;
}

void gcpMolecule::UpdateCycles()
{
	Lock (true);
	std::list<gcpBond*>::iterator n;
	for (n = m_Bonds.begin(); n != m_Bonds.end(); n++) (*n)->RemoveAllCycles();
	while (!m_Cycles.empty())
	{
		delete m_Cycles.front();
		m_Cycles.pop_front();
	}
	if (!m_Atoms.empty())
	{
		std::list<gcpAtom*>::iterator i = m_Atoms.begin();
		i++;
		for (; i != m_Atoms.end(); i++) (*i)->SetParent(NULL);
		gcpChain* pChain = new gcpChain(this, *(m_Atoms.begin())); //will find the cycles
		delete pChain;
	}
	Lock (false);
}

void gcpMolecule::Clear()
{
	m_Bonds.clear();
	m_Atoms.clear();
	m_Fragments.clear();
}

void gcpMolecule::SetSelected (GtkWidget* w, int state)
{
	map<string, Object*>::iterator i;
	Object *child = GetFirstChild (i);
	while (child)
	{
		child->SetSelected (w, state);
		child = GetNextChild (i);
	}
}

void gcpMolecule::Transform2D (Matrix2D& m, double x, double y)
{
	Object::Transform2D (m, x, y);
	std::list<gcpAtom*>::iterator i = m_Atoms.begin ();
	for (; i != m_Atoms.end (); i++)
	if (((*i)->GetZ () != 6) && (*i)->GetAttachedHydrogens () &&
		(*i)->GetBondsNumber ()) (*i)->Update ();
}

Object* gcpMolecule::GetAtomAt(double x, double y, double z)
{
	// Make use of gcpBond::GetAtomAt
	std::list<gcpBond*>::iterator n, end = m_Bonds.end();
	Object* pObj = NULL;
	for (n = m_Bonds.begin(); n != end; n++)
		if ((pObj = (*n)->GetAtomAt (x, y)))
			break;
	return pObj;
}
	
double gcpMolecule::GetYAlign ()
{
	if (m_Alignment)
		return m_Alignment->GetYAlign ();
	double y, maxy = - DBL_MAX, miny = DBL_MAX;
	std::list<gcpAtom*>::iterator i = m_Atoms.begin(), end = m_Atoms.end();
	for (; i != end; i++) {
		y = (*i)->GetYAlign ();
		if (y < miny)
			miny = y;
		if (y > maxy)
			maxy = y;
	}
	std::list<gcpFragment*>::iterator ig = m_Fragments.begin(), endg = m_Fragments.end();
	for (; ig != endg; ig++) {
		y = (*ig)->GetYAlign ();
		if (y < miny)
			miny = y;
		if (y > maxy)
			maxy = y;
	}
	return (miny + maxy) / 2.0;
}

bool gcpMolecule::BuildContextualMenu (GtkUIManager *UIManager, Object *object, double x, double y)
{
	bool result = false;
	GtkActionGroup *group = gtk_action_group_new ("molecule");
	GtkAction *action;
	action = gtk_action_new ("Molecule", _("Molecule"), NULL, NULL);
	gtk_action_group_add_action (group, action);
	if (!m_Fragments.size ()) {
		if (((gcpDocument*) GetDocument ())->GetApplication ()->HaveGhemical ()) {
			action = gtk_action_new ("ghemical", _("Export molecule to Ghemical"), NULL, NULL);
			g_signal_connect_swapped (action, "activate", G_CALLBACK (do_export_to_ghemical), this);
			gtk_action_group_add_action (group, action);
			g_object_unref (action);
			gtk_ui_manager_add_ui_from_string (UIManager, "<ui><popup><menu action='Molecule'><menuitem action='ghemical'/></menu></popup></ui>", -1, NULL);
		}
		if (((gcpDocument*) GetDocument ())->GetApplication ()->HaveInChI ()) {
			action = gtk_action_new ("inchi", _("Generate InChI"), NULL, NULL);
			g_signal_connect_swapped (action, "activate", G_CALLBACK (do_build_inchi), this);
			gtk_action_group_add_action (group, action);
			g_object_unref (action);
			gtk_ui_manager_add_ui_from_string (UIManager, "<ui><popup><menu action='Molecule'><menuitem action='inchi'/></menu></popup></ui>", -1, NULL);
			action = gtk_action_new ("webbook", _("NIST WebBook page for this molecule"), NULL, NULL);
			g_signal_connect_swapped (action, "activate", G_CALLBACK (do_show_webbook), this);
			gtk_action_group_add_action (group, action);
			g_object_unref (action);
			gtk_ui_manager_add_ui_from_string (UIManager, "<ui><popup><menu action='Molecule'><menuitem action='webbook'/></menu></popup></ui>", -1, NULL);
			action = gtk_action_new ("pubchem", _("PubChem page for this molecule"), NULL, NULL);
			g_signal_connect_swapped (action, "activate", G_CALLBACK (do_show_pubchem), this);
			gtk_action_group_add_action (group, action);
			g_object_unref (action);
			gtk_ui_manager_add_ui_from_string (UIManager, "<ui><popup><menu action='Molecule'><menuitem action='pubchem'/></menu></popup></ui>", -1, NULL);
		}
		action = gtk_action_new ("smiles", _("Generate Smiles"), NULL, NULL);
		g_signal_connect_swapped (action, "activate", G_CALLBACK (do_build_smiles), this);
		gtk_action_group_add_action (group, action);
		g_object_unref (action);
		gtk_ui_manager_add_ui_from_string (UIManager, "<ui><popup><menu action='Molecule'><menuitem action='smiles'/></menu></popup></ui>", -1, NULL);
		action = gtk_action_new ("calc", _("Open in Calculator"), NULL, NULL);
		g_signal_connect_swapped (action, "activate", G_CALLBACK (do_open_in_calc), this);
		gtk_action_group_add_action (group, action);
		g_object_unref (action);
		gtk_ui_manager_add_ui_from_string (UIManager, "<ui><popup><menu action='Molecule'><menuitem action='calc'/></menu></popup></ui>", -1, NULL);
		result = true;
	}
	if (m_Bonds.size ()) {
		action = gtk_action_new ("select-align", _("Select alignment item"), NULL, NULL);
		g_signal_connect (action, "activate", G_CALLBACK (do_select_alignment), this);
		g_object_set_data (G_OBJECT (action), "item", object);
		gtk_action_group_add_action (group, action);
		g_object_unref (action);
		gtk_ui_manager_add_ui_from_string (UIManager, "<ui><popup><menu action='Molecule'><menuitem action='select-align'/></menu></popup></ui>", -1, NULL);
		result = true;
	}
	if (result)
		gtk_ui_manager_insert_action_group (UIManager, group, 0);
	g_object_unref (group);
	return result | GetParent ()->BuildContextualMenu (UIManager, object, x, y);
}

void gcpMolecule::ExportToGhemical ()
{
	OBMol Mol;
	OBConversion Conv;
	OBFormat* pOutFormat = Conv.FindFormat ("gpr");
	Conv.SetInAndOutFormats (pOutFormat, pOutFormat);
	BuildOBMol (Mol);
	char *tmpname = g_strdup ("/tmp/gcp2gprXXXXXX");
	int f = g_mkstemp (tmpname);
	gchar *old_num_locale;
	close (f);
	ofstream ofs;
	ofs.open(tmpname);
	if (!ofs) throw (int) 1;
	old_num_locale = g_strdup(setlocale(LC_NUMERIC, NULL));
	setlocale(LC_NUMERIC, "C");
	Conv.Write (&Mol, &ofs);
	setlocale(LC_NUMERIC, old_num_locale);
	g_free(old_num_locale);
	ofs.close();
	char *command_line = g_strconcat ("ghemical -f ", tmpname, NULL);	
	g_free (tmpname);
	g_spawn_command_line_async (command_line, NULL);
	g_free (command_line);
}

void gcpMolecule::SelectAlignmentItem (Object *child)
{
	m_Alignment = (m_Alignment != child)? child: NULL;
	EmitSignal (OnChangedSignal);
}

xmlNodePtr gcpMolecule::Save (xmlDocPtr xml)
{
	xmlNodePtr node = Object::Save (xml);
	if (!node)
		return NULL;
	if (m_Alignment)
		xmlNewProp (node, (const xmlChar*) "valign", (const xmlChar*) m_Alignment->GetId ());
	return node;
}

void gcpMolecule::BuildOBMol (OBMol &Mol)
{
	double xav = 0., yav = 0., zf;
	unsigned n = m_Atoms.size();
	map<string, unsigned> AtomTable;
	list<gcpAtom*>::iterator ia, enda = m_Atoms.end ();
	list<gcpBond*> BondList;
	double x, y, z;
	for (ia = m_Atoms.begin(); ia != enda; ia++) {
		(*ia)->GetCoords(&x, &y, &z);
		xav += x;
		yav += y;
	}
	xav /= n;
	yav /= n;
	gcpAtom* pgAtom;
	OBAtom obAtom;
	unsigned index = 1;
	map<Atom*, Bond*>::iterator i;
	gcpBond *pBond;
	Mol.BeginModify();
	Mol.ReserveAtoms(n);
	for (ia = m_Atoms.begin(); ia != enda; ia++) {
		pgAtom = *ia;
		AtomTable [pgAtom->GetId ()] = index;
		obAtom.SetIdx(index++);
		obAtom.SetAtomicNum(pgAtom->GetZ());
		pgAtom->GetCoords(&x, &y, &z);
		// Scans the atom bonds and change z to try conservation of stereochemistry
		pBond = (gcpBond*) pgAtom->GetFirstBond (i);
		while (pBond) {
			zf = (pBond->GetAtom (0) == pgAtom)? 1.: -1.;
			switch (pBond->GetType ()) {
			case UpBondType:
				z += zf * 50.;
				break;
			case DownBondType:
				z -= zf * 50.;
				break;
			default:
				break;
			}
			pBond = (gcpBond*) pgAtom->GetNextBond (i);
		}
		obAtom.SetVector((xav - x )/ 100, (yav - y) / 100, z / 100);
		Mol.AddAtom(obAtom);
		obAtom.Clear();
	}
	list<gcpBond*>::iterator j, endb = m_Bonds.end ();
	int start, end, order;
	for (j = m_Bonds.begin(); j != endb; j++)
	{
		order = (*j)->GetOrder();
		start = AtomTable[(*j)->GetAtom(0)->GetId()];
		end = AtomTable[(*j)->GetAtom(1)->GetId()];
		Mol.AddBond(start, end, order, 0);
	}
	Mol.EndModify();
}

void gcpMolecule::BuildOBMol2D (OBMol &Mol)
{
	double xav = 0., yav = 0.;
	unsigned n = m_Atoms.size();
	map<string, unsigned> AtomTable;
	list<gcpAtom*>::iterator ia, enda = m_Atoms.end ();
	list<gcpBond*> BondList;
	double x, y, z;
	for (ia = m_Atoms.begin(); ia != enda; ia++) {
		(*ia)->GetCoords(&x, &y, &z);
		xav += x;
		yav += y;
	}
	xav /= n;
	yav /= n;
	gcpAtom* pgAtom;
	OBAtom obAtom;
	unsigned index = 1;
	map<Atom*, Bond*>::iterator i;
	Mol.BeginModify();
	Mol.ReserveAtoms(n);
	Mol.SetDimension (2);
	for (ia = m_Atoms.begin(); ia != enda; ia++) {
		pgAtom = *ia;
		AtomTable [pgAtom->GetId ()] = index;
		obAtom.SetIdx(index++);
		obAtom.SetAtomicNum(pgAtom->GetZ());
		pgAtom->GetCoords(&x, &y, &z);
		// Scans the atom bonds and change z to try conservation of stereochemistry
		obAtom.SetVector((xav - x )/ 100, (yav - y) / 100, 0.);
		Mol.AddAtom(obAtom);
		obAtom.Clear();
	}
	list<gcpBond*>::iterator j, endb = m_Bonds.end ();
	int start, end, order, flag;
	for (j = m_Bonds.begin(); j != endb; j++)
	{
		order = (*j)->GetOrder();
		start = AtomTable[(*j)->GetAtom(0)->GetId()];
		end = AtomTable[(*j)->GetAtom(1)->GetId()];
		switch((*j)->GetType())
		{
			case UpBondType:
				flag = OB_WEDGE_BOND;
				break;
			case DownBondType:
				flag = OB_HASH_BOND;
				break;
			default:
				flag = 0;
		}
		Mol.AddBond(start, end, order, flag);
	}
	Mol.EndModify();
}

void gcpMolecule::BuildInChI ()
{
	OBMol Mol;
	OBConversion Conv;
	BuildOBMol2D (Mol);
	OBFormat *pInChIFormat = Conv.FindFormat ("inchi"), *pMolFormat = Conv.FindFormat ("mol");
	if (pInChIFormat) {
		Conv.SetInAndOutFormats (pMolFormat, pInChIFormat);
		Conv.SetOptions ("xt", OpenBabel::OBConversion::OUTOPTIONS);
		ostringstream ofs;
		char *old_num_locale = g_strdup (setlocale (LC_NUMERIC, NULL));
		setlocale (LC_NUMERIC, "C");
		Conv.Write (&Mol, &ofs);
		setlocale (LC_NUMERIC, old_num_locale);
		g_free (old_num_locale);
		// remove "INCHI=" and the new line char at the end
		m_InChI = ofs.str ().substr (0, ofs.str ().length () - 2);
	} else {
		Conv.SetInAndOutFormats (pMolFormat, pMolFormat);
		char *tmpname = g_strdup ("/tmp/inchiXXXXXX");
		int f = g_mkstemp (tmpname);
		close (f);
		ofstream ofs;
		ofs.open (tmpname);
		char *old_num_locale = g_strdup (setlocale (LC_NUMERIC, NULL));
		setlocale (LC_NUMERIC, "C");
		Conv.Write (&Mol, &ofs);
		setlocale (LC_NUMERIC, old_num_locale);
		ofs.close ();
		// calling main_inchi -STDIO -AuxNone -NoLabels
		char *command = g_strdup_printf ("main_inchi %s -STDIO -AuxNone -NoLabels", tmpname);
		char *output, *errors;
		g_spawn_command_line_sync (command, &output, &errors, NULL, NULL);
		if (output) {
			// remove the new line char at the end
			output[strlen (output) - 1] = 0;
			m_InChI = output + 6;
			g_free (output);
		}
		if (errors)
			g_free (errors);
		g_free (command);
		g_free (old_num_locale);
		remove (tmpname);
		g_free (tmpname);
	}
	m_Changed = false;
}

void gcpMolecule::BuildSmiles ()
{
	OBMol Mol;
	OBConversion Conv;
	OBFormat* pOutFormat = Conv.FindFormat ("smi");
	Conv.SetInAndOutFormats (pOutFormat, pOutFormat);
	BuildOBMol2D (Mol);
	ostringstream ofs;
	char *old_num_locale = g_strdup (setlocale (LC_NUMERIC, NULL));
	setlocale (LC_NUMERIC, "C");
	Conv.Write (&Mol, &ofs);
	setlocale (LC_NUMERIC, old_num_locale);
	//TODO: do something with the string
	g_free (old_num_locale);
	string str = ofs.str ().substr (0, ofs.str ().length () - 2);
	new gcpStringDlg (reinterpret_cast<gcpDocument *>(GetDocument ()), str, gcpStringDlg::SMILES);
}

void gcpMolecule::ShowInChI ()
{
	if (m_Changed)
		BuildInChI ();
	new gcpStringDlg (reinterpret_cast<gcpDocument *>(GetDocument ()), m_InChI, gcpStringDlg::INCHI);
}

void gcpMolecule::ShowWebBase (char const* uri_start, char const *uri_end)
{
	if (m_Changed)
		BuildInChI ();
	if (m_InChI.length () == 0)
		return; //should emit at least a warning
	string::size_type t;
	while ((t = m_InChI.find ('+')) != string::npos)
		m_InChI.replace (t, 1, "%2b");
	string uri = string (uri_start) + m_InChI + uri_end;
	((gcpDocument*) GetDocument ())->GetApplication ()->ShowURI (uri);
}

bool gcpMolecule::OnSignal (SignalId Signal, Object *Child)
{
	m_Changed = true;
	return true;
}

void gcpMolecule::OpenCalc ()
{
	list<gcpAtom*>::iterator ia, enda = m_Atoms.end ();
	ostringstream ofs;
	int nH;
	ofs << "gchemcalc ";
	for (ia = m_Atoms.begin(); ia != enda; ia++) {
		ofs << (*ia)->GetSymbol();
		nH = (*ia)->GetAttachedHydrogens ();
		if (nH > 0) {
			ofs << "H" << nH;
		}
	}
	g_spawn_command_line_async (ofs.str ().c_str (), NULL);
}

void gcpMolecule::CheckCrossings (gcpBond *pBond)
{
	gcpView *pView = reinterpret_cast<gcpDocument*> (GetDocument ())->GetView ();
	list<gcpBond*>::iterator i, iend = m_Bonds.end ();
	for (i = m_Bonds.begin (); i != iend; i++)
		if (((*i) != pBond) && (*i)->IsCrossing (pBond)) {
			pView->Update (pBond);
			pView->Update (*i);
		}
}

void gcpMolecule::Add (GtkWidget* w)
{
	/* Now that Object::Add adds children, we need to override this method
	because a molecule children (atoms, bonds, and fragments) are added to
	the widget by the document. FIXME! FIXME! FIXME! CHANGE THIS OLD WEIRD CODE!
	*/
}