ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/openbabel/mol.cpp

Comparing trunk/src/openbabel/mol.cpp (file contents):
Revision 741 by tim, Wed Nov 16 19:42:11 2005 UTC vs.
Revision 1081 by gezelter, Thu Oct 19 20:49:05 2006 UTC

#	Line 36 | Line 36 | namespace OpenBabel
36		namespace OpenBabel
37		{
38
39	<	extern bool SwabInt;
40	<	extern OBPhModel        phmodel;
41	<	extern OBAromaticTyper  aromtyper;
42	<	extern OBAtomTyper      atomtyper;
43	<	extern OBBondTyper      bondtyper;
39	>	extern bool SwabInt;
40	>	extern OBPhModel      phmodel;
41	>	extern OBAromaticTyper  aromtyper;
42	>	extern OBAtomTyper      atomtyper;
43	>	extern OBBondTyper      bondtyper;
44
45
46	<	//
47	<	// OBMol member functions
48	<	//
49	<	void  OBMol::SetTitle(const char *title)
50	<	{
51	<	_title = title;
52	<	Trim(_title);
53	<	}
46	>	/** \class OBMol
47	>	\brief Molecule Class
48	>
49	>	The most important class in Open Babel is OBMol, or the molecule class.
50	>	The OBMol class is designed to store all the basic information
51	>	associated with a molecule, to make manipulations on the connection
52	>	table of a molecule facile, and to provide member functions which
53	>	automatically perceive information about a molecule. A guided tour
54	>	of the OBMol class is a good place to start.
55	>
56	>	An OBMol class can be declared:
57	>	\code
58	>	OBMol mol;
59	>	\endcode
60
61	<	void  OBMol::SetTitle(std::string &title)
62	<	{
63	<	_title = title;
64	<	Trim(_title);
65	<	}
61	>	For example:
62	>	\code
63	>	#include <iostream.h>
64	>	#include "mol.hpp"
65	>	#include "obconversion.hpp"
66	>	int main(int argc,char **argv)
67	>	{
68	>	OBConversion conv(&cin,&cout);
69	>	if(conv.SetInAndOutFormats("SDF","MOL2"))
70	>	{
71	>	OBMol mol;
72	>	if(conv.Read(&mol))
73	>	...manipulate molecule
74	>
75	>	conv->Write(&mol);
76	>	}
77	>	return(1);
78	>	}
79	>	\endcode
80	>
81	>	will read in a molecule in SD file format from stdin
82	>	(or the C++ equivalent cin) and write a MOL2 format file out
83	>	to standard out. Additionally, The input and output formats can
84	>	be altered using the OBConversion class
85
86	<	bool SortVVInt(const vector<int> &a,const vector<int> &b)
87	<	{
86	>	Once a molecule has been read into an OBMol (or created via other methods)
87	>	the atoms and bonds
88	>	can be accessed by the following methods:
89	>	\code
90	>	OBAtom *atom;
91	>	atom = mol.GetAtom(5); //random access of an atom
92	>	\endcode
93	>	or
94	>	\code
95	>	OBBond *bond;
96	>	bond = mol.GetBond(14); //random access of a bond
97	>	\endcode
98	>	or
99	>	\code
100	>	FOR_ATOMS_IN_MOL(atom, mol) // iterator access (see OBMolAtomIter)
101	>	\endcode
102	>	or
103	>	\code
104	>	FOR_BONDS_IN_MOL(bond, mol) // iterator access (see OBMolBondIter)
105	>	\endcode
106	>	It is important to note that atom arrays currently begin at 1 and bond arrays
107	>	begin at 0. Requesting atom 0 (\code
108	>	OBAtom *atom = mol.GetAtom(0); \endcode
109	>	will result in an error, but
110	>	\code
111	>	OBBond *bond = mol.GetBond(0);
112	>	\endcode
113	>	is perfectly valid.
114	>	Note that this is expected to change in the near future to simplify coding
115	>	and improve efficiency.
116	>
117	>	The ambiguity of numbering issues and off-by-one errors led to the use
118	>	of iterators in Open Babel. An iterator is essentially just a pointer, but
119	>	when used in conjunction with Standard Template Library (STL) vectors
120	>	it provides an unambiguous way to loop over arrays. OBMols store their
121	>	atom and bond information in STL vectors. Since vectors are template
122	>	based, a vector of any user defined type can be declared. OBMols declare
123	>	vector<OBNodeBase*> and vector<OBEdgeBase*> to store atom and bond information.
124	>	Iterators are then a natural way to loop over the vectors of atoms and bonds.
125	>
126	>	A variety of predefined iterators have been created to simplify
127	>	common looping requests (e.g., looping over all atoms in a molecule,
128	>	bonds to a given atom, etc.)
129	>
130	>	\code
131	>	#include "obiter.hpp"
132	>	...
133	>	#define FOR_ATOMS_OF_MOL(a,m)     for( OBMolAtomIter     a(m); a; a++ )
134	>	#define FOR_BONDS_OF_MOL(b,m)     for( OBMolBondIter     b(m); b; b++ )
135	>	#define FOR_NBORS_OF_ATOM(a,p)    for( OBAtomAtomIter    a(p); a; a++ )
136	>	#define FOR_BONDS_OF_ATOM(b,p)    for( OBAtomBondIter    b(p); b; b++ )
137	>	#define FOR_RESIDUES_OF_MOL(r,m)  for( OBResidueIter     r(m); r; r++ )
138	>	#define FOR_ATOMS_OF_RESIDUE(a,r) for( OBResidueAtomIter a(r); a; a++ )
139	>	...
140	>	\endcode
141	>
142	>	These convenience functions can be used like so:
143	>	\code
144	>	#include "obiter.hpp"
145	>	#include "mol.hpp"
146	>
147	>	OBMol mol;
148	>	double exactMass = 0.0f;
149	>	FOR_ATOMS_OF_MOL(a, mol)
150	>	{
151	>	exactMass +=  a->GetExactMass();
152	>	}
153	>	\endcode
154	>
155	>	Note that with these convenience macros, the iterator "a" (or
156	>	whichever name you pick) is declared for you -- you do not need to
157	>	do it beforehand.
158	>	*/
159	>
160	>	//
161	>	// OBMol member functions
162	>	//
163	>	void  OBMol::SetTitle(const char *title)
164	>	{
165	>	_title = title;
166	>	Trim(_title);
167	>	}
168	>
169	>	void  OBMol::SetTitle(std::string &title)
170	>	{
171	>	_title = title;
172	>	Trim(_title);
173	>	}
174	>
175	>	bool SortVVInt(const vector<int> &a,const vector<int> &b)
176	>	{
177		return(a.size() > b.size());
178	<	}
178	>	}
179
180	<	bool SortAtomZ(const pair<OBAtom*,double> &a, const pair<OBAtom*,double> &b)
181	<	{
180	>	bool SortAtomZ(const pair<OBAtom*,double> &a, const pair<OBAtom*,double> &b)
181	>	{
182		return (a.second < b.second);
183	<	}
183	>	}
184
185	<	double OBMol::GetTorsion(int a,int b,int c,int d)
186	<	{
185	>	double OBMol::GetTorsion(int a,int b,int c,int d)
186	>	{
187		return(CalcTorsionAngle(((OBAtom*)_vatom[a-1])->GetVector(),
188		((OBAtom*)_vatom[b-1])->GetVector(),
189		((OBAtom*)_vatom[c-1])->GetVector(),
190		((OBAtom*)_vatom[d-1])->GetVector()));
191	<	}
191	>	}
192
193	<	void OBMol::SetTorsion(OBAtom *a,OBAtom *b,OBAtom *c, OBAtom *d, double ang)
194	<	{
193	>	void OBMol::SetTorsion(OBAtom *a,OBAtom *b,OBAtom *c, OBAtom *d, double ang)
194	>	{
195		vector<int> tor;
196		vector<int> atoms;
197
198	<	obErrorLog.ThrowError(__FUNCTION__,
198	>	obErrorLog.ThrowError(__func__,
199		"Ran OpenBabel::SetTorsion", obAuditMsg);
200
201		tor.push_back(a->GetCIdx());
#	Line 92 | Line 206 | void OBMol::SetTorsion(OBAtom *a,OBAtom *b,OBAtom *c,
206		FindChildren(atoms, b->GetIdx(), c->GetIdx());
207		int j;
208		for (j = 0 ; (unsigned)j < atoms.size() ; j++ )
209	<	atoms[j] = (atoms[j] - 1) * 3;
209	>	atoms[j] = (atoms[j] - 1) * 3;
210
211		double v1x,v1y,v1z,v2x,v2y,v2z,v3x,v3y,v3z;
212		double c1x,c1y,c1z,c2x,c2y,c2z,c3x,c3y,c3z;
#	Line 124 | Line 238 | void OBMol::SetTorsion(OBAtom *a,OBAtom *b,OBAtom *c,
238		c1mag = SQUARE(c1x)+SQUARE(c1y)+SQUARE(c1z);
239		c2mag = SQUARE(c2x)+SQUARE(c2y)+SQUARE(c2z);
240		if (c1mag*c2mag < 0.01)
241	<	costheta = 1.0; //avoid div by zero error
241	>	costheta = 1.0; //avoid div by zero error
242		else
243	<	costheta = (c1x*c2x + c1y*c2y + c1z*c2z)/(sqrt(c1mag*c2mag));
243	>	costheta = (c1x*c2x + c1y*c2y + c1z*c2z)/(sqrt(c1mag*c2mag));
244
245		if (costheta < -0.999999)
246	<	costheta = -0.999999;
246	>	costheta = -0.999999;
247		if (costheta >  0.999999)
248	<	costheta =  0.999999;
248	>	costheta =  0.999999;
249
250		if ((v2x*c3x + v2y*c3y + v2z*c3z) > 0.0)
251	<	radang = -acos(costheta);
251	>	radang = -acos(costheta);
252		else
253	<	radang = acos(costheta);
253	>	radang = acos(costheta);
254
255		//
256		// now we have the torsion angle (radang) - set up the rot matrix
#	Line 173 | Line 287 | void OBMol::SetTorsion(OBAtom *a,OBAtom *b,OBAtom *c,
287		tz = _c[tor[1]+2];
288		vector<int>::iterator i;
289		for (i = atoms.begin(),j=*i;i != atoms.end();i++,j=*i)
290	<	{
290	>	{
291		_c[j] -= tx;
292		_c[j+1] -= ty;
293		_c[j+2]-= tz;
#	Line 186 | Line 300 | void OBMol::SetTorsion(OBAtom *a,OBAtom *b,OBAtom *c,
300		_c[j] += tx;
301		_c[j+1] += ty;
302		_c[j+2] += tz;
303	<	}
304	<	}
303	>	}
304	>	}
305
306
307	<	double OBMol::GetTorsion(OBAtom *a,OBAtom *b,OBAtom *c,OBAtom *d)
308	<	{
307	>	double OBMol::GetTorsion(OBAtom *a,OBAtom *b,OBAtom *c,OBAtom *d)
308	>	{
309		return(CalcTorsionAngle(a->GetVector(),
310		b->GetVector(),
311		c->GetVector(),
312		d->GetVector()));
313	<	}
313	>	}
314
315	<	void OBMol::ContigFragList(std::vector<std::vector<int> >&cfl)
316	<	{
315	>	void OBMol::ContigFragList(std::vector<std::vector<int> >&cfl)
316	>	{
317		int j;
318		OBAtom *atom;
319		OBBond *bond;
#	Line 214 | Line 328 | void OBMol::ContigFragList(std::vector<std::vector<int
328		frag.Resize(NumAtoms()+1);
329
330		while ((unsigned)used.CountBits() < NumAtoms())
331	<	{
331	>	{
332		curr.Clear();
333		frag.Clear();
334		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
335	<	if (!used.BitIsOn(atom->GetIdx()))
335	>	if (!used.BitIsOn(atom->GetIdx()))
336		{
337	<	curr.SetBitOn(atom->GetIdx());
338	<	break;
337	>	curr.SetBitOn(atom->GetIdx());
338	>	break;
339		}
340
341		frag |= curr;
342		while (!curr.IsEmpty())
343	<	{
343	>	{
344		next.Clear();
345		for (j = curr.NextBit(-1);j != curr.EndBit();j = curr.NextBit(j))
346	<	{
346	>	{
347		atom = GetAtom(j);
348		for (bond = atom->BeginBond(k);bond;bond = atom->NextBond(k))
349	<	if (!used.BitIsOn(bond->GetNbrAtomIdx(atom)))
350	<	next.SetBitOn(bond->GetNbrAtomIdx(atom));
351	<	}
349	>	if (!used.BitIsOn(bond->GetNbrAtomIdx(atom)))
350	>	next.SetBitOn(bond->GetNbrAtomIdx(atom));
351	>	}
352
353		used |= curr;
354		used |= next;
355		frag |= next;
356		curr = next;
357	<	}
357	>	}
358
359		tmp.clear();
360		frag.ToVecInt(tmp);
361		cfl.push_back(tmp);
362	<	}
362	>	}
363
364		sort(cfl.begin(),cfl.end(),SortVVInt);
365	<	}
365	>	}
366
367	<	/*!
368	<	**\brief Fills the OBGeneric OBTorsionData with torsions from the mol
369	<	*/
370	<	void OBMol::FindTorsions()
371	<	{
367	>	/*!
368	>	**\brief Fills the OBGeneric OBTorsionData with torsions from the mol
369	>	*/
370	>	void OBMol::FindTorsions()
371	>	{
372		//if already has data return
373		if(HasData(OBGenericDataType::TorsionData))
374	<	return;
374	>	return;
375
376		//get new data and attach it to molecule
377		OBTorsionData *torsions = new OBTorsionData;
#	Line 270 | Line 384 | void OBMol::FindTorsions()
384
385		//loop through all bonds generating torsions
386		for(bond = BeginBond(bi1);bond;bond = NextBond(bi1))
387	<	{
387	>	{
388		b = bond->GetBeginAtom();
389		c = bond->GetEndAtom();
390		if(b->IsHydrogen() || c->IsHydrogen())
391	<	continue;
391	>	continue;
392
393		for(a = b->BeginNbrAtom(bi2);a;a = b->NextNbrAtom(bi2))
394	<	{
394	>	{
395		if(a == c)
396	<	continue;
396	>	continue;
397
398		for(d = c->BeginNbrAtom(bi3);d;d = c->NextNbrAtom(bi3))
399	<	{
399	>	{
400		if(d == b)
401	<	continue;
401	>	continue;
402		torsion.AddTorsion(a,b,c,d);
403	<	}
404	<	}
403	>	}
404	>	}
405		//add torsion to torsionData
406		if(torsion.GetSize())
407	<	torsions->SetData(torsion);
407	>	torsions->SetData(torsion);
408		torsion.Clear();
409	<	}
409	>	}
410
411		return;
412	<	}
412	>	}
413
414	<	void OBMol::FindLargestFragment(OBBitVec &lf)
415	<	{
414	>	void OBMol::FindLargestFragment(OBBitVec &lf)
415	>	{
416		int j;
417		OBAtom *atom;
418		OBBond *bond;
#	Line 308 | Line 422 | void OBMol::FindLargestFragment(OBBitVec &lf)
422
423		lf.Clear();
424		while ((unsigned)used.CountBits() < NumAtoms())
425	<	{
425	>	{
426		curr.Clear();
427		frag.Clear();
428		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
429	<	if (!used.BitIsOn(atom->GetIdx()))
429	>	if (!used.BitIsOn(atom->GetIdx()))
430		{
431	<	curr.SetBitOn(atom->GetIdx());
432	<	break;
431	>	curr.SetBitOn(atom->GetIdx());
432	>	break;
433		}
434
435		frag |= curr;
436		while (!curr.IsEmpty())
437	<	{
437	>	{
438		next.Clear();
439		for (j = curr.NextBit(-1);j != curr.EndBit();j = curr.NextBit(j))
440	<	{
440	>	{
441		atom = GetAtom(j);
442		for (bond = atom->BeginBond(k);bond;bond = atom->NextBond(k))
443	<	if (!used.BitIsOn(bond->GetNbrAtomIdx(atom)))
444	<	next.SetBitOn(bond->GetNbrAtomIdx(atom));
445	<	}
443	>	if (!used.BitIsOn(bond->GetNbrAtomIdx(atom)))
444	>	next.SetBitOn(bond->GetNbrAtomIdx(atom));
445	>	}
446
447		used |= curr;
448		used |= next;
449		frag |= next;
450		curr = next;
451	<	}
451	>	}
452
453		if (lf.Empty() || lf.CountBits() < frag.CountBits())
454	<	lf = frag;
455	<	}
456	<	}
454	>	lf = frag;
455	>	}
456	>	}
457
458	<	//! locates all atoms for which there exists a path to 'end'
459	<	//! without going through 'bgn'
460	<	//! children must not include 'end'
461	<	void OBMol::FindChildren(vector<OBAtom*> &children,OBAtom *bgn,OBAtom *end)
462	<	{
458	>	//! locates all atoms for which there exists a path to 'end'
459	>	//! without going through 'bgn'
460	>	//! children must not include 'end'
461	>	void OBMol::FindChildren(vector<OBAtom*> &children,OBAtom *bgn,OBAtom *end)
462	>	{
463		OBBitVec used,curr,next;
464
465		used |= bgn->GetIdx();
#	Line 358 | Line 472 | void OBMol::FindChildren(vector<OBAtom*> &children,OBA
472		vector<OBEdgeBase*>::iterator j;
473
474		for (;;)
475	<	{
475	>	{
476		next.Clear();
477		for (i = curr.NextBit(-1);i != curr.EndBit();i = curr.NextBit(i))
478	<	{
478	>	{
479		atom = GetAtom(i);
480		for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
481	<	if (!used[nbr->GetIdx()])
481	>	if (!used[nbr->GetIdx()])
482		{
483	<	children.push_back(nbr);
484	<	next |= nbr->GetIdx();
485	<	used |= nbr->GetIdx();
483	>	children.push_back(nbr);
484	>	next |= nbr->GetIdx();
485	>	used |= nbr->GetIdx();
486		}
487	<	}
487	>	}
488		if (next.Empty())
489	<	break;
489	>	break;
490		curr = next;
491	<	}
492	<	}
491	>	}
492	>	}
493
494	<	//! locates all atoms for which there exists a path to 'second'
495	<	//! without going through 'first'
496	<	//! children must not include 'second'
497	<	void OBMol::FindChildren(vector<int> &children,int first,int second)
498	<	{
494	>	//! locates all atoms for which there exists a path to 'second'
495	>	//! without going through 'first'
496	>	//! children must not include 'second'
497	>	void OBMol::FindChildren(vector<int> &children,int first,int second)
498	>	{
499		int i;
500		OBBitVec used,curr,next;
501
#	Line 394 | Line 508 | void OBMol::FindChildren(vector<int> &children,int fir
508		vector<OBEdgeBase*>::iterator j;
509
510		while (!curr.IsEmpty())
511	<	{
511	>	{
512		next.Clear();
513		for (i = curr.NextBit(-1);i != curr.EndBit();i = curr.NextBit(i))
514	<	{
514	>	{
515		atom = GetAtom(i);
516		for (j = atom->BeginBonds(),bond=(OBBond *)*j;
517	<	j != atom->EndBonds();j++,bond=(OBBond *)*j)
518	<	if (!used.BitIsOn(bond->GetNbrAtomIdx(atom)))
519	<	next.SetBitOn(bond->GetNbrAtomIdx(atom));
520	<	}
517	>	j != atom->EndBonds();j++,bond=(OBBond *)*j)
518	>	if (!used.BitIsOn(bond->GetNbrAtomIdx(atom)))
519	>	next.SetBitOn(bond->GetNbrAtomIdx(atom));
520	>	}
521
522		used |= next;
523		curr = next;
524	<	}
524	>	}
525
526		used.SetBitOff(first);
527		used.SetBitOff(second);
528		used.ToVecInt(children);
529	<	}
529	>	}
530
531	<	/*!
532	<	**\brief Calculates the graph theoretical distance of each atom.
533	<	** Vector is indexed from zero
534	<	*/
535	<	bool OBMol::GetGTDVector(vector<int> &gtd)
536	<	//calculates the graph theoretical distance for every atom
537	<	//and puts it into gtd
538	<	{
531	>	/*!
532	>	**\brief Calculates the graph theoretical distance of each atom.
533	>	** Vector is indexed from zero
534	>	*/
535	>	bool OBMol::GetGTDVector(vector<int> &gtd)
536	>	//calculates the graph theoretical distance for every atom
537	>	//and puts it into gtd
538	>	{
539		gtd.clear();
540		gtd.resize(NumAtoms());
541
#	Line 435 | Line 549 | bool OBMol::GetGTDVector(vector<int> &gtd)
549		next.Clear();
550
551		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
552	<	{
552	>	{
553		gtdcount = 0;
554		used.Clear();
555		curr.Clear();
#	Line 443 | Line 557 | bool OBMol::GetGTDVector(vector<int> &gtd)
557		curr.SetBitOn(atom->GetIdx());
558
559		while (!curr.IsEmpty())
560	<	{
560	>	{
561		next.Clear();
562		for (natom = curr.NextBit(-1);natom != curr.EndBit();natom = curr.NextBit(natom))
563	<	{
563	>	{
564		atom1 = GetAtom(natom);
565		for (bond = atom1->BeginBond(j);bond;bond = atom1->NextBond(j))
566	<	if (!used.BitIsOn(bond->GetNbrAtomIdx(atom1)) && !curr.BitIsOn(bond->GetNbrAtomIdx(atom1)))
567	<	if (!(bond->GetNbrAtom(atom1))->IsHydrogen())
568	<	next.SetBitOn(bond->GetNbrAtomIdx(atom1));
569	<	}
566	>	if (!used.BitIsOn(bond->GetNbrAtomIdx(atom1)) && !curr.BitIsOn(bond->GetNbrAtomIdx(atom1)))
567	>	if (!(bond->GetNbrAtom(atom1))->IsHydrogen())
568	>	next.SetBitOn(bond->GetNbrAtomIdx(atom1));
569	>	}
570
571		used |= next;
572		curr = next;
573		gtdcount++;
574	<	}
574	>	}
575		gtd[atom->GetIdx()-1] = gtdcount;
576	<	}
576	>	}
577		return(true);
578	<	}
578	>	}
579
580	<	/*!
581	<	**\brief Calculates a set of graph invariant indexes using
582	<	** the graph theoretical distance, number of connected heavy atoms,
583	<	** aromatic boolean, ring boolean, atomic number, and
584	<	** summation of bond orders connected to the atom.
585	<	** Vector is indexed from zero
586	<	*/
587	<	void OBMol::GetGIVector(vector<unsigned int> &vid)
588	<	{
580	>	/*!
581	>	**\brief Calculates a set of graph invariant indexes using
582	>	** the graph theoretical distance, number of connected heavy atoms,
583	>	** aromatic boolean, ring boolean, atomic number, and
584	>	** summation of bond orders connected to the atom.
585	>	** Vector is indexed from zero
586	>	*/
587	>	void OBMol::GetGIVector(vector<unsigned int> &vid)
588	>	{
589		vid.clear();
590		vid.resize(NumAtoms()+1);
591
#	Line 482 | Line 596 | void OBMol::GetGIVector(vector<unsigned int> &vid)
596		OBAtom *atom;
597		vector<OBNodeBase*>::iterator j;
598		for (i=0,atom = BeginAtom(j);atom;atom = NextAtom(j),i++)
599	<	{
599	>	{
600		vid[i] =  (unsigned int)v[i];
601		vid[i] += (unsigned int)(atom->GetHvyValence()*100);
602		vid[i] += (unsigned int)(((atom->IsAromatic()) ? 1 : 0)*1000);
603		vid[i] += (unsigned int)(((atom->IsInRing()) ? 1 : 0)*10000);
604		vid[i] += (unsigned int)(atom->GetAtomicNum()*100000);
605		vid[i] += (unsigned int)(atom->GetImplicitValence()*10000000);
606	<	}
607	<	}
606	>	}
607	>	}
608
609	<	static bool OBComparePairSecond(const pair<OBAtom*,unsigned int> &a,const pair<OBAtom*,unsigned int> &b)
610	<	{
609	>	static bool OBComparePairSecond(const pair<OBAtom*,unsigned int> &a,const pair<OBAtom*,unsigned int> &b)
610	>	{
611		return(a.second < b.second);
612	<	}
612	>	}
613
614	<	static bool OBComparePairFirst(const pair<OBAtom*,unsigned int> &a,const pair<OBAtom*,unsigned int> &b)
615	<	{
614	>	static bool OBComparePairFirst(const pair<OBAtom*,unsigned int> &a,const pair<OBAtom*,unsigned int> &b)
615	>	{
616		return(a.first->GetIdx() < b.first->GetIdx());
617	<	}
617	>	}
618
619	<	//! counts the number of unique symmetry classes in a list
620	<	static void ClassCount(vector<pair<OBAtom*,unsigned int> > &vp,unsigned int &count)
621	<	{
619	>	//! counts the number of unique symmetry classes in a list
620	>	static void ClassCount(vector<pair<OBAtom*,unsigned int> > &vp,unsigned int &count)
621	>	{
622		count = 0;
623		vector<pair<OBAtom*,unsigned int> >::iterator k;
624		sort(vp.begin(),vp.end(),OBComparePairSecond);
625	+	#if 0 // original version
626
627	+	unsigned int id=0; // [ejk] appease gcc's bogus "might be undef'd" warning
628	+	for (k = vp.begin();k != vp.end();k++)
629	+	{
630	+	if (k == vp.begin())
631	+	{
632	+	id = k->second;
633	+	k->second = count = 0;
634	+	}
635	+	else
636	+	if (k->second != id)
637	+	{
638	+	id = k->second;
639	+	k->second = ++count;
640	+	}
641	+	else
642	+	k->second = count;
643	+	}
644	+	count++;
645	+	#else // get rid of warning, moves test out of loop, returns 0 for empty input
646	+
647		k = vp.begin();
648		if (k != vp.end())
649	<	{
649	>	{
650		unsigned int id = k->second;
651		k->second = 0;
652		++k;
653		for (;k != vp.end(); ++k)
654	<	{
654	>	{
655		if (k->second != id)
656	<	{
656	>	{
657		id = k->second;
658		k->second = ++count;
659	<	}
659	>	}
660		else
661	<	k->second = count;
662	<	}
661	>	k->second = count;
662	>	}
663		++count;
664	<	}
664	>	}
665		else
666	<	{
666	>	{
667		// [ejk] thinks count=0 might be OK for an empty list, but orig code did
668		//++count;
669	<	}
670	<	}
669	>	}
670	>	#endif
671	>	}
672
673	<	//! creates a new vector of symmetry classes base on an existing vector
674	<	//! helper routine to GetGIDVector
675	<	static void     CreateNewClassVector(vector<pair<OBAtom*,unsigned int> > &vp1,vector<pair<OBAtom*,unsigned int> > &vp2)
676	<	{
673	>	//! creates a new vector of symmetry classes base on an existing vector
674	>	//! helper routine to GetGIDVector
675	>	static void   CreateNewClassVector(vector<pair<OBAtom*,unsigned int> > &vp1,vector<pair<OBAtom*,unsigned int> > &vp2)
676	>	{
677		int m,id;
678		OBAtom *nbr;
679		vector<OBEdgeBase*>::iterator j;
#	Line 546 | Line 682 | static void    CreateNewClassVector(vector<pair<OBAtom*,u
682		sort(vp1.begin(),vp1.end(),OBComparePairFirst);
683		vp2.clear();
684		for (i = vp1.begin();i != vp1.end();i++)
685	<	{
685	>	{
686		vector<unsigned int> vtmp;
687		for (nbr = i->first->BeginNbrAtom(j);nbr;nbr = i->first->NextNbrAtom(j))
688	<	vtmp.push_back(vp1[nbr->GetIdx()-1].second);
688	>	vtmp.push_back(vp1[nbr->GetIdx()-1].second);
689		sort(vtmp.begin(),vtmp.end(),OBCompareUnsigned);
690		for (id=i->second,m=100,k = vtmp.begin();k != vtmp.end();k++,m*=100)
691	<	id += *k * m;
691	>	id += *k * m;
692
693		vp2.push_back(pair<OBAtom*,unsigned int> (i->first,id));
694	<	}
695	<	}
694	>	}
695	>	}
696
697	<	/*!
698	<	**\brief Calculates a set of symmetry identifiers for a molecule.
699	<	** Atoms with the same symmetry ID are symmetrically equivalent.
700	<	** Vector is indexed from zero
701	<	*/
702	<	void OBMol::GetGIDVector(vector<unsigned int> &vgid)
703	<	{
697	>	/*!
698	>	**\brief Calculates a set of symmetry identifiers for a molecule.
699	>	** Atoms with the same symmetry ID are symmetrically equivalent.
700	>	** Vector is indexed from zero
701	>	*/
702	>	void OBMol::GetGIDVector(vector<unsigned int> &vgid)
703	>	{
704		vector<unsigned int> vgi;
705		GetGIVector(vgi);  //get vector of graph invariants
706
#	Line 573 | Line 709 | void OBMol::GetGIDVector(vector<unsigned int> &vgid)
709		vector<OBNodeBase*>::iterator j;
710		vector<pair<OBAtom*,unsigned int> > vp1,vp2;
711		for (i=0,atom = BeginAtom(j);atom;atom = NextAtom(j),i++)
712	<	vp1.push_back(pair<OBAtom*,unsigned int> (atom,vgi[i]));
712	>	vp1.push_back(pair<OBAtom*,unsigned int> (atom,vgi[i]));
713
714		unsigned int nclass1,nclass2; //number of classes
715		ClassCount(vp1,nclass1);
716
717		if (nclass1 < NumAtoms())
718	<	{
718	>	{
719		for (i = 0;i < 100;i++) //sanity check - shouldn't ever hit this number
720	<	{
720	>	{
721		CreateNewClassVector(vp1,vp2);
722		ClassCount(vp2,nclass2);
723		vp1 = vp2;
724		if (nclass1 == nclass2)
725	<	break;
725	>	break;
726		nclass1 = nclass2;
727	<	}
728	<	}
727	>	}
728	>	}
729
730		vgid.clear();
731		sort(vp1.begin(),vp1.end(),OBComparePairFirst);
732		vector<pair<OBAtom*,unsigned int> >::iterator k;
733		for (k = vp1.begin();k != vp1.end();k++)
734	<	vgid.push_back(k->second);
735	<	}
734	>	vgid.push_back(k->second);
735	>	}
736
737	<	unsigned int OBMol::NumHvyAtoms()
738	<	{
737	>	unsigned int OBMol::NumHvyAtoms()
738	>	{
739		OBAtom *atom;
740		vector<OBNodeBase*>::iterator(i);
741		unsigned int count = 0;
742
743		for(atom = this->BeginAtom(i);atom;atom = this->NextAtom(i))
744	<	{
744	>	{
745		if(!atom->IsHydrogen())
746	<	count++;
747	<	}
746	>	count++;
747	>	}
748
749		return(count);
750	<	}
750	>	}
751
752	<	unsigned int OBMol::NumRotors()
753	<	{
752	>	unsigned int OBMol::NumRotors()
753	>	{
754		OBBond *bond;
755		vector<OBEdgeBase*>::iterator i;
756
757		unsigned int count = 0;
758		for (bond = BeginBond(i);bond;bond = NextBond(i))
759	<	if (bond->IsRotor())
760	<	count++;
759	>	if (bond->IsRotor())
760	>	count++;
761
762		return(count);
763	<	}
763	>	}
764
765	<	//! Returns a pointer to the atom after a safety check
766	<	//! 0 < idx <= NumAtoms
767	<	OBAtom *OBMol::GetAtom(int idx)
768	<	{
765	>	//! Returns a pointer to the atom after a safety check
766	>	//! 0 < idx <= NumAtoms
767	>	OBAtom *OBMol::GetAtom(int idx)
768	>	{
769		if ((unsigned)idx < 1 || (unsigned)idx > NumAtoms())
770	<	{
771	<	obErrorLog.ThrowError(__FUNCTION__, "Requested Atom Out of Range", obDebug);
770	>	{
771	>	obErrorLog.ThrowError(__func__, "Requested Atom Out of Range", obDebug);
772		return((OBAtom*)NULL);
773	<	}
773	>	}
774
775		return((OBAtom*)_vatom[idx-1]);
776	<	}
776	>	}
777
778	<	OBAtom *OBMol::GetFirstAtom()
779	<	{
778	>	OBAtom *OBMol::GetFirstAtom()
779	>	{
780		return((_vatom.empty()) ? (OBAtom*)NULL : (OBAtom*)_vatom[0]);
781	<	}
781	>	}
782
783	<	//! Returns a pointer to the bond after a safety check
784	<	//! 0 <= idx < NumBonds
785	<	OBBond *OBMol::GetBond(int idx)
786	<	{
783	>	//! Returns a pointer to the bond after a safety check
784	>	//! 0 <= idx < NumBonds
785	>	OBBond *OBMol::GetBond(int idx)
786	>	{
787		if (idx < 0 || (unsigned)idx >= NumBonds())
788	<	{
789	<	obErrorLog.ThrowError(__FUNCTION__, "Requested Bond Out of Range", obDebug);
788	>	{
789	>	obErrorLog.ThrowError(__func__, "Requested Bond Out of Range", obDebug);
790		return((OBBond*)NULL);
791	<	}
791	>	}
792
793		return((OBBond*)_vbond[idx]);
794	<	}
794	>	}
795
796	<	OBBond *OBMol::GetBond(int bgn, int end)
797	<	{
796	>	OBBond *OBMol::GetBond(int bgn, int end)
797	>	{
798		return(GetBond(GetAtom(bgn),GetAtom(end)));
799	<	}
799	>	}
800
801	<	OBBond *OBMol::GetBond(OBAtom *bgn,OBAtom *end)
802	<	{
801	>	OBBond *OBMol::GetBond(OBAtom *bgn,OBAtom *end)
802	>	{
803		OBAtom *nbr;
804		vector<OBEdgeBase*>::iterator i;
805
806		for (nbr = bgn->BeginNbrAtom(i);nbr;nbr = bgn->NextNbrAtom(i))
807	<	if (nbr == end)
808	<	return((OBBond *)*i);
807	>	if (nbr == end)
808	>	return((OBBond *)*i);
809
810		return(NULL); //just to keep the SGI compiler happy
811	<	}
811	>	}
812
813	<	OBResidue *OBMol::GetResidue(int idx)
814	<	{
813	>	OBResidue *OBMol::GetResidue(int idx)
814	>	{
815		if (idx < 0 || (unsigned)idx >= NumResidues())
816	<	{
817	<	obErrorLog.ThrowError(__FUNCTION__, "Requested Residue Out of Range", obDebug);
816	>	{
817	>	obErrorLog.ThrowError(__func__, "Requested Residue Out of Range", obDebug);
818		return((OBResidue*)NULL);
819	<	}
819	>	}
820
821		return (_residue[idx]);
822	<	}
822	>	}
823
824	<	std::vector<OBInternalCoord*> OBMol::GetInternalCoord()
825	<	{
824	>	std::vector<OBInternalCoord*> OBMol::GetInternalCoord()
825	>	{
826		if (_internals.empty())
827	<	{
827	>	{
828		_internals.push_back((OBInternalCoord*)NULL);
829		for(unsigned int i = 1; i <= NumAtoms(); i++)
830	<	{
830	>	{
831		_internals.push_back(new OBInternalCoord);
832	<	}
832	>	}
833		CartesianToInternal(_internals, *this);
834	<	}
834	>	}
835		return _internals;
836	<	}
836	>	}
837
838	<	//! Implements <a href="http://qsar.sourceforge.net/dicts/blue-obelisk/index.xhtml#findSmallestSetOfSmallestRings">blue-obelisk:findSmallestSetOfSmallestRings</a>.
839	<	vector<OBRing*> &OBMol::GetSSSR()
840	<	{
838	>	//! Implements <a href="http://qsar.sourceforge.net/dicts/blue-obelisk/index.xhtml#findSmallestSetOfSmallestRings">blue-obelisk:findSmallestSetOfSmallestRings</a>.
839	>	vector<OBRing*> &OBMol::GetSSSR()
840	>	{
841		if (!HasSSSRPerceived())
842	<	FindSSSR();
842	>	FindSSSR();
843
844		if (!HasData(OBGenericDataType::RingData))
845	<	SetData(new OBRingData);
845	>	SetData(new OBRingData);
846
847		OBRingData *rd = (OBRingData *) GetData(OBGenericDataType::RingData);
848		return(rd->GetData());
849	<	}
849	>	}
850
851	<	double OBMol::GetMolWt()
852	<	{
853	<	double molwt=0.0;
851	>	double OBMol::GetMolWt()
852	>	{
853	>	double mass = 0.0;
854		OBAtom *atom;
855		vector<OBNodeBase*>::iterator i;
856
857	+	bool UseImplicitH = NumHvyAtoms() && (NumBonds()!=0 || NumAtoms()==1);
858		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
859	<	molwt += atom->GetAtomicMass();
859	>	{
860	>	if(UseImplicitH)
861	>	{
862	>	if (!atom->IsHydrogen())
863	>	mass += isotab.GetExactMass(1,1) * atom->ImplicitHydrogenCount();
864	>	}
865	>	mass += atom->GetAtomicMass();
866	>	}
867	>	return(mass);
868	>	}
869
870	<	return(molwt);
871	<	}
726	<
727	<	double OBMol::GetExactMass()
728	<	{
870	>	double OBMol::GetExactMass()
871	>	{
872		double mass=0.0;
873		OBAtom *atom;
874		vector<OBNodeBase*>::iterator i;
875
876	+	bool UseImplicitH = NumHvyAtoms() && (NumBonds()!=0 || NumAtoms()==1);
877		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
878	<	mass += atom->GetExactMass();
879	<
878	>	{
879	>	if(UseImplicitH)
880	>	{
881	>	if (!atom->IsHydrogen())
882	>	mass += isotab.GetExactMass(1,1) * atom->ImplicitHydrogenCount();
883	>	}
884	>	mass += atom->GetExactMass();
885	>	}
886		return(mass);
887	<	}
887	>	}
888
889	<	//! Stochoimetric formula (e.g., C4H6O).
890	<	//!   This is either set by OBMol::SetFormula() or generated on-the-fly
891	<	//!   using the "Hill order" -- i.e., C first if present, then H if present
892	<	//!   all other elements in alphabetical order.
893	<	string OBMol::GetFormula()
894	<	{
895	<	string attr = "Formula";
896	<	OBPairData *dp = (OBPairData *) GetData(attr);
889	>	//! Stochoimetric formula (e.g., C4H6O).
890	>	//!   This is either set by OBMol::SetFormula() or generated on-the-fly
891	>	//!   using the "Hill order" -- i.e., C first if present, then H if present
892	>	//!   all other elements in alphabetical order.
893	>	string OBMol::GetFormula()
894	>	{
895	>	string attr = "Formula";
896	>	OBPairData *dp = (OBPairData *) GetData(attr);
897
898	<	if (dp != NULL) // we already set the formula
899	<	return dp->GetValue();
898	>	if (dp != NULL) // we already set the formula
899	>	return dp->GetValue();
900
901	<	obErrorLog.ThrowError(__FUNCTION__,
902	<	"Ran OpenBabel::SetFormula -- Hill order formula",
903	<	obAuditMsg);
901	>	obErrorLog.ThrowError(__func__,
902	>	"Ran OpenBabel::SetFormula -- Hill order formula",
903	>	obAuditMsg);
904
905	<	// OK, now let's generate the formula and store it for future use.
906	<	// These are the atomic numbers of the elements in alphabetical order.
907	<	const int NumElements = 110;
908	<	const int alphabetical[NumElements] = {
909	<	89, 47, 13, 95, 18, 33, 85, 79, 5, 56, 4, 107, 83, 97, 35, 6, 20, 48,
910	<	58, 98, 17, 96, 27, 24, 55, 29, 105, 66, 68, 99, 63, 9, 26, 100, 87, 31,
911	<	64, 32, 1, 2, 72, 80, 67, 108, 53, 49, 77, 19, 36, 57, 3, 103, 71, 101,
912	<	12, 25, 42, 109, 7, 11, 41, 60, 10, 28, 102, 93, 8, 76, 15, 91, 82, 46,
913	<	61, 84, 59, 78, 94, 88, 37, 75, 104, 45, 86, 44, 16, 51, 21, 34, 106, 14,
914	<	62, 50, 38, 73, 65, 43, 52, 90, 22, 81, 69, 92, 110, 23, 74, 54, 39, 70,
915	<	30, 40 };
905	>	// OK, now let's generate the formula and store it for future use.
906	>	// These are the atomic numbers of the elements in alphabetical order.
907	>	const int NumElements = 110;
908	>	const int alphabetical[NumElements] = {
909	>	89, 47, 13, 95, 18, 33, 85, 79, 5, 56, 4, 107, 83, 97, 35, 6, 20, 48,
910	>	58, 98, 17, 96, 27, 24, 55, 29, 105, 66, 68, 99, 63, 9, 26, 100, 87, 31,
911	>	64, 32, 1, 2, 72, 80, 67, 108, 53, 49, 77, 19, 36, 57, 3, 103, 71, 101,
912	>	12, 25, 42, 109, 7, 11, 41, 60, 10, 28, 102, 93, 8, 76, 15, 91, 82, 46,
913	>	61, 84, 59, 78, 94, 88, 37, 75, 104, 45, 86, 44, 16, 51, 21, 34, 106, 14,
914	>	62, 50, 38, 73, 65, 43, 52, 90, 22, 81, 69, 92, 110, 23, 74, 54, 39, 70,
915	>	30, 40 };
916
917	<	int atomicCount[NumElements];
918	<	//  int index;
917	>	int atomicCount[NumElements];
918	>	//  int index;
919		#ifdef HAVE_SSTREAM
920	<	stringstream formula;
920	>	stringstream formula;
921		#else
922	<	strstream formula;
922	>	strstream formula;
923		#endif
924
925	<	for (int i = 0; i < NumElements; i++)
926	<	atomicCount[i] = 0;
925	>	for (int i = 0; i < NumElements; i++)
926	>	atomicCount[i] = 0;
927
928	<	FOR_ATOMS_OF_MOL(a, *this)
929	<	atomicCount[a->GetAtomicNum() - 1]++;
928	>	FOR_ATOMS_OF_MOL(a, *this)
929	>	{
930	>	int anum = a->GetAtomicNum();
931	>	if (anum == 1) continue;    // skip explicit hydrogens
932	>	atomicCount[anum - 1]++;
933	>	atomicCount[0] += a->ImplicitHydrogenCount() + a->ExplicitHydrogenCount();
934	>	}
935
936	<	if (atomicCount[5] != 0) // Carbon (i.e. 6 - 1 = 5)
937	<	{
938	<	if (atomicCount[5] > 1)
939	<	formula << "C" << atomicCount[5];
940	<	else if (atomicCount[5] == 1)
941	<	formula << "C";
936	>	if (atomicCount[5] != 0) // Carbon (i.e. 6 - 1 = 5)
937	>	{
938	>	if (atomicCount[5] > 1)
939	>	formula << "C" << atomicCount[5];
940	>	else if (atomicCount[5] == 1)
941	>	formula << "C";
942
943	<	atomicCount[5] = 0; // So we don't output C twice
943	>	atomicCount[5] = 0; // So we don't output C twice
944
945	<	// only output H if there's also carbon -- otherwise do it alphabetical
946	<	if (atomicCount[0] != 0) // Hydrogen (i.e., 1 - 1 = 0)
947	<	{
948	<	if (atomicCount[0] > 1)
949	<	formula << "H" << atomicCount[0];
950	<	else if (atomicCount[0] == 1)
951	<	formula << "H";
945	>	// only output H if there's also carbon -- otherwise do it alphabetical
946	>	if (atomicCount[0] != 0) // Hydrogen (i.e., 1 - 1 = 0)
947	>	{
948	>	if (atomicCount[0] > 1)
949	>	formula << "H" << atomicCount[0];
950	>	else if (atomicCount[0] == 1)
951	>	formula << "H";
952
953	<	atomicCount[0] = 0;
954	<	}
955	<	}
953	>	atomicCount[0] = 0;
954	>	}
955	>	}
956
957	<	for (int j = 0; j < NumElements; j++)
958	<	{
959	<	if (atomicCount[ alphabetical[j]-1 ] > 1)
960	<	formula << etab.GetSymbol(alphabetical[j])
961	<	<< atomicCount[ alphabetical[j]-1 ];
962	<	else if (atomicCount[ alphabetical[j]-1 ] == 1)
963	<	formula << etab.GetSymbol( alphabetical[j] );
964	<	}
957	>	for (int j = 0; j < NumElements; j++)
958	>	{
959	>	if (atomicCount[ alphabetical[j]-1 ] > 1)
960	>	formula << etab.GetSymbol(alphabetical[j])
961	>	<< atomicCount[ alphabetical[j]-1 ];
962	>	else if (atomicCount[ alphabetical[j]-1 ] == 1)
963	>	formula << etab.GetSymbol( alphabetical[j] );
964	>	}
965
966	<	dp = new OBPairData;
967	<	dp->SetAttribute(attr);
968	<	dp->SetValue( formula.str() );
969	<	SetData(dp);
966	>	dp = new OBPairData;
967	>	dp->SetAttribute(attr);
968	>	dp->SetValue( formula.str() );
969	>	SetData(dp);
970
971	<	return (formula.str());
972	<	}
971	>	return (formula.str());
972	>	}
973
974	<	void OBMol::SetFormula(string molFormula)
975	<	{
976	<	string attr = "Formula";
977	<	OBPairData *dp = (OBPairData *) GetData(attr);
974	>	void OBMol::SetFormula(string molFormula)
975	>	{
976	>	string attr = "Formula";
977	>	OBPairData *dp = (OBPairData *) GetData(attr);
978
979	<	if (dp == NULL)
980	<	{
981	<	dp = new OBPairData;
982	<	dp->SetAttribute(attr);
983	<	}
984	<	dp->SetValue(molFormula);
979	>	if (dp == NULL)
980	>	{
981	>	dp = new OBPairData;
982	>	dp->SetAttribute(attr);
983	>	}
984	>	dp->SetValue(molFormula);
985
986	<	SetData(dp);
987	<	}
986	>	SetData(dp);
987	>	}
988
989	<	void OBMol::SetTotalCharge(int charge)
990	<	{
989	>	void OBMol::SetTotalCharge(int charge)
990	>	{
991		SetFlag(OB_TCHARGE_MOL);
992		_totalCharge = charge;
993	<	}
993	>	}
994
995	<	//! Returns the total molecular charge -- if it has not previously been set
996	<	//!  it is calculated from the atomic formal charge information.
997	<	//!  (This may or may not be correct!)
998	<	//!  If you set atomic charges with OBAtom::SetFormalCharge()
999	<	//!   you really should set the molecular charge with OBMol::SetTotalCharge()
1000	<	int OBMol::GetTotalCharge()
1001	<	{
995	>	//! Returns the total molecular charge -- if it has not previously been set
996	>	//!  it is calculated from the atomic formal charge information.
997	>	//!  (This may or may not be correct!)
998	>	//!  If you set atomic charges with OBAtom::SetFormalCharge()
999	>	//!   you really should set the molecular charge with OBMol::SetTotalCharge()
1000	>	int OBMol::GetTotalCharge()
1001	>	{
1002		if(HasFlag(OB_TCHARGE_MOL))
1003	<	return(_totalCharge);
1003	>	return(_totalCharge);
1004		else // calculate from atomic formal charges (seems the best default)
1005	<	{
1006	<	obErrorLog.ThrowError(__FUNCTION__,
1007	<	"Ran OpenBabel::GetTotalCharge -- calculated from formal charges",
1008	<	obAuditMsg);
1005	>	{
1006	>	obErrorLog.ThrowError(__func__,
1007	>	"Ran OpenBabel::GetTotalCharge -- calculated from formal charges",
1008	>	obAuditMsg);
1009
1010		OBAtom *atom;
1011		vector<OBNodeBase*>::iterator i;
1012		int chg = 0;
1013
1014		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
1015	<	chg += atom->GetFormalCharge();
1015	>	chg += atom->GetFormalCharge();
1016		return (chg);
1017	<	}
1018	<	}
1017	>	}
1018	>	}
1019
1020	<	void   OBMol::SetTotalSpinMultiplicity(unsigned int spin)
1021	<	{
1020	>	void   OBMol::SetTotalSpinMultiplicity(unsigned int spin)
1021	>	{
1022		SetFlag(OB_TSPIN_MOL);
1023		_totalSpin = spin;
1024	<	}
1024	>	}
1025
1026	<	//! Returns the total spin multiplicity -- if it has not previously been set
1027	<	//!  it is calculated from the atomic spin multiplicity information
1028	<	//!  assuming the high-spin case (i.e. it simply sums the atomic spins,
1029	<	//!  making no attempt to pair spins).
1030	<	//!  However, if you set atomic spins with OBAtom::SetSpinMultiplicity()
1031	<	//!   you really should set the molecular spin with
1032	<	//!   OBMol::SetTotalSpinMultiplicity()
1033	<	unsigned int OBMol::GetTotalSpinMultiplicity()
1034	<	{
1026	>	//! Returns the total spin multiplicity -- if it has not previously been set
1027	>	//!  it is calculated from the atomic spin multiplicity information
1028	>	//!  assuming the high-spin case (i.e. it simply sums the atomic spins,
1029	>	//!  making no attempt to pair spins).
1030	>	//!  However, if you set atomic spins with OBAtom::SetSpinMultiplicity()
1031	>	//!   you really should set the molecular spin with
1032	>	//!   OBMol::SetTotalSpinMultiplicity()
1033	>	unsigned int OBMol::GetTotalSpinMultiplicity()
1034	>	{
1035		if (HasFlag(OB_TSPIN_MOL))
1036	<	return(_totalSpin);
1036	>	return(_totalSpin);
1037		else // calculate from atomic spin information (assuming high-spin case)
1038	<	{
1039	<	obErrorLog.ThrowError(__FUNCTION__,
1038	>	{
1039	>	obErrorLog.ThrowError(__func__,
1040		"Ran OpenBabel::GetTotalSpinMultiplicity -- calculating from atomic spins and high spin",
1041		obAuditMsg);
1042
#	Line 890 | Line 1045 | unsigned int OBMol::GetTotalSpinMultiplicity()
1045		unsigned int spin = 1;
1046
1047		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
1048	<	{
1048	>	{
1049		if (atom->GetSpinMultiplicity() > 1)
1050	<	spin += atom->GetSpinMultiplicity() - 1;
1051	<	}
1050	>	spin += atom->GetSpinMultiplicity() - 1;
1051	>	}
1052		return (spin);
1053	<	}
1054	<	}
1053	>	}
1054	>	}
1055
1056	<	OBMol &OBMol::operator=(const OBMol &source)
1057	<	//only atom and bond info is copied from src to dest
1058	<	//Conformers are now copied also, MM 2/7/01
1059	<	//Rotamers and residue information are copied, MM 4-27-01
1060	<	{
1056	>	OBMol &OBMol::operator=(const OBMol &source)
1057	>	//only atom and bond info is copied from src to dest
1058	>	//Conformers are now copied also, MM 2/7/01
1059	>	//Rotamers and residue information are copied, MM 4-27-01
1060	>	{
1061		OBMol &src = (OBMol &)source;
1062		vector<OBNodeBase*>::iterator i;
1063		vector<OBEdgeBase*>::iterator j;
#	Line 916 | Line 1071 | OBMol &OBMol::operator=(const OBMol &source)
1071		_vbond.reserve(src.NumBonds());
1072
1073		for (atom = src.BeginAtom(i);atom;atom = src.NextAtom(i))
1074	<	AddAtom(*atom);
1074	>	AddAtom(*atom);
1075		for (bond = src.BeginBond(j);bond;bond = src.NextBond(j))
1076	<	AddBond(*bond);
1076	>	AddBond(*bond);
1077
1078		this->_title  = src.GetTitle();
1079		this->_energy = src.GetEnergy();
#	Line 928 | Line 1083 | OBMol &OBMol::operator=(const OBMol &source)
1083		//Copy Residue information
1084		unsigned int NumRes = src.NumResidues();
1085		if (NumRes)
1086	<	{
1086	>	{
1087		unsigned int k;
1088		OBResidue *src_res=NULL;
1089		OBResidue *res=NULL;
#	Line 936 | Line 1091 | OBMol &OBMol::operator=(const OBMol &source)
1091		OBAtom *atom=NULL;
1092		vector<OBAtom*>::iterator ii;
1093		for (k=0 ; k<NumRes ; k++)
1094	<	{
1094	>	{
1095		res = NewResidue();
1096		src_res = src.GetResidue(k);
1097		res->SetName(src_res->GetName());
#	Line 944 | Line 1099 | OBMol &OBMol::operator=(const OBMol &source)
1099		res->SetChain(src_res->GetChain());
1100		res->SetChainNum(src_res->GetChainNum());
1101		for (src_atom=src_res->BeginAtom(ii) ; src_atom ; src_atom=src_res->NextAtom(ii))
1102	<	{
1102	>	{
1103		atom = GetAtom(src_atom->GetIdx());
1104		res->AddAtom(atom);
1105		res->SetAtomID(atom,src_res->GetAtomID(src_atom));
1106		res->SetHetAtom(atom,src_res->IsHetAtom(src_atom));
1107		res->SetSerialNum(atom,src_res->GetSerialNum(src_atom));
1108	<	}
1109	<	}
1110	<	}
1108	>	}
1109	>	}
1110	>	}
1111
1112		//Copy conformer information
1113		if (src.NumConformers() > 1)
1114	<	{
1114	>	{
1115		int k,l;
1116		vector<double*> conf;
1117		double* xyz = NULL;
1118		for (k=0 ; k<src.NumConformers() ; k++)
1119	<	{
1119	>	{
1120		xyz = new double [3*src.NumAtoms()];
1121		for (l=0 ; l<(int) (3*src.NumAtoms()) ; l++)
1122	<	xyz[l] = src.GetConformer(k)[l];
1122	>	xyz[l] = src.GetConformer(k)[l];
1123		conf.push_back(xyz);
1124	<	}
1124	>	}
1125		SetConformers(conf);
1126	<	}
1126	>	}
1127
1128		//Copy rotamer list
1129		OBRotamerList *rml = (OBRotamerList *)src.GetData(OBGenericDataType::RotamerList);
1130		if (rml && rml->NumAtoms() == src.NumAtoms())
1131	<	{
1131	>	{
1132		//Destroy old rotamer list if necessary
1133		if ((OBRotamerList *)GetData(OBGenericDataType::RotamerList))
1134	<	{
1134	>	{
1135		DeleteData(OBGenericDataType::RotamerList);
1136	<	}
1136	>	}
1137
1138		//Set base coordinates
1139		OBRotamerList *cp_rml = new OBRotamerList;
#	Line 987 | Line 1142 | OBMol &OBMol::operator=(const OBMol &source)
1142		double *c=NULL;
1143		double *cc=NULL;
1144		for (k=0 ; k<rml->NumBaseCoordinateSets() ; k++)
1145	<	{
1145	>	{
1146		c = new double [3*rml->NumAtoms()];
1147		cc = rml->GetBaseCoordinateSet(k);
1148		for (l=0 ; l<3*rml->NumAtoms() ; l++)
1149	<	c[l] = cc[l];
1149	>	c[l] = cc[l];
1150		bc.push_back(c);
1151	<	}
1151	>	}
1152		if (rml->NumBaseCoordinateSets())
1153	<	cp_rml->SetBaseCoordinateSets(bc,rml->NumAtoms());
1153	>	cp_rml->SetBaseCoordinateSets(bc,rml->NumAtoms());
1154
1155		//Set reference array
1156		unsigned char *ref = new unsigned char [rml->NumRotors()*4];
1157		if (ref)
1158	<	{
1158	>	{
1159		rml->GetReferenceArray(ref);
1160		cp_rml->Setup((*this),ref,rml->NumRotors());
1161		delete [] ref;
1162	<	}
1162	>	}
1163
1164		//Set Rotamers
1165		unsigned char *rotamers = new unsigned char [(rml->NumRotors()+1)*rml->NumRotamers()];
1166		if (rotamers)
1167	<	{
1167	>	{
1168		vector<unsigned char*>::iterator kk;
1169		unsigned int idx=0;
1170		for (kk = rml->BeginRotamer();kk != rml->EndRotamer();kk++)
1171	<	{
1171	>	{
1172		memcpy(&rotamers[idx],(const unsigned char*)*kk,sizeof(unsigned char)*(rml->NumRotors()+1));
1173		idx += sizeof(unsigned char)*(rml->NumRotors()+1);
1174	<	}
1174	>	}
1175		cp_rml->AddRotamers(rotamers,rml->NumRotamers());
1176		delete [] rotamers;
1177	<	}
1177	>	}
1178		SetData(cp_rml);
1179	<	}
1179	>	}
1180
1181		return(*this);
1182	<	}
1182	>	}
1183
1184	<	OBMol &OBMol::operator+=(const OBMol &source)
1185	<	{
1184	>	OBMol &OBMol::operator+=(const OBMol &source)
1185	>	{
1186		OBMol &src = (OBMol &)source;
1187		vector<OBNodeBase*>::iterator i;
1188		vector<OBEdgeBase*>::iterator j;
#	Line 1041 | Line 1196 | OBMol &OBMol::operator+=(const OBMol &source)
1196		_title += "_" + string(src.GetTitle());
1197
1198		for (atom = src.BeginAtom(i) ; atom ; atom = src.NextAtom(i))
1199	<	AddAtom(*atom);
1199	>	AddAtom(*atom);
1200		for (bond = src.BeginBond(j) ; bond ; bond = src.NextBond(j))
1201	<	AddBond(bond->GetBeginAtomIdx() + prevatms, bond->GetEndAtomIdx() + prevatms, bond->GetBO(), bond->GetFlags());
1201	>	AddBond(bond->GetBeginAtomIdx() + prevatms, bond->GetEndAtomIdx() + prevatms, bond->GetBO(), bond->GetFlags());
1202
1203		EndModify();
1204
1205		return(*this);
1206	<	}
1206	>	}
1207
1208	<	bool OBMol::Clear()
1209	<	{
1210	<	obErrorLog.ThrowError(__FUNCTION__,
1211	<	"Ran OpenBabel::Clear Molecule", obAuditMsg);
1208	>	bool OBMol::Clear()
1209	>	{
1210	>	obErrorLog.ThrowError(__func__,
1211	>	"Ran OpenBabel::Clear Molecule", obAuditMsg);
1212
1213		vector<OBNodeBase*>::iterator i;
1214		vector<OBEdgeBase*>::iterator j;
1215		for (i = _vatom.begin();i != _vatom.end();i++)
1216	<	{
1216	>	{
1217		DestroyAtom(*i);
1218		*i = NULL;
1219	<	}
1219	>	}
1220		for (j = _vbond.begin();j != _vbond.end();j++)
1221	<	{
1221	>	{
1222		DestroyBond(*j);
1223		*j = NULL;
1224	<	}
1224	>	}
1225
1226		_natoms = _nbonds = 0;
1227
1228		//Delete residues
1229		unsigned int ii;
1230		for (ii=0 ; ii<_residue.size() ; ii++)
1231	<	{
1231	>	{
1232		delete _residue[ii];
1233	<	}
1233	>	}
1234		_residue.clear();
1235
1236		//clear out the multiconformer data
1237		vector<double*>::iterator k;
1238		for (k = _vconf.begin();k != _vconf.end();k++)
1239	<	delete [] *k;
1239	>	delete [] *k;
1240		_vconf.clear();
1241
1242		if (!_vdata.empty()) //clean up generic data
1243	<	{
1243	>	{
1244		vector<OBGenericData*>::iterator m;
1245		for (m = _vdata.begin();m != _vdata.end();m++)
1246	<	delete *m;
1246	>	delete *m;
1247		_vdata.clear();
1248	<	}
1248	>	}
1249
1250		_c = (double*) NULL;
1251		_flags = 0;
1252		_mod = 0;
1253
1254		return(true);
1255	<	}
1255	>	}
1256
1257	<	void OBMol::BeginModify()
1258	<	{
1257	>	void OBMol::BeginModify()
1258	>	{
1259		//suck coordinates from _c into _v for each atom
1260		if (!_mod && !Empty())
1261	<	{
1261	>	{
1262		OBAtom *atom;
1263		vector<OBNodeBase*>::iterator i;
1264		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
1265	<	{
1265	>	{
1266		atom->SetVector();
1267		atom->ClearCoordPtr();
1268	<	}
1268	>	}
1269
1270		vector<double*>::iterator j;
1271		for (j = _vconf.begin();j != _vconf.end();j++)
1272	<	delete [] *j;
1272	>	delete [] *j;
1273
1274		_c = NULL;
1275		_vconf.clear();
1276
1277		//Destroy rotamer list if necessary
1278		if ((OBRotamerList *)GetData(OBGenericDataType::RotamerList))
1279	<	{
1279	>	{
1280		delete (OBRotamerList *)GetData(OBGenericDataType::RotamerList);
1281		DeleteData(OBGenericDataType::RotamerList);
1282	<	}
1283	<	}
1282	>	}
1283	>	}
1284
1285		_mod++;
1286	<	}
1286	>	}
1287
1288	<	void OBMol::EndModify(bool nukePerceivedData)
1289	<	{
1288	>	void OBMol::EndModify(bool nukePerceivedData)
1289	>	{
1290		if (_mod == 0)
1291	<	{
1292	<	obErrorLog.ThrowError(__FUNCTION__, "_mod is negative - EndModify() called too many times", obDebug);
1291	>	{
1292	>	obErrorLog.ThrowError(__func__, "_mod is negative - EndModify() called too many times", obDebug);
1293		return;
1294	<	}
1294	>	}
1295
1296		_mod--;
1297
1298		if (_mod)
1299	<	return;
1299	>	return;
1300
1301		if (nukePerceivedData)
1302	<	_flags = 0;
1302	>	_flags = 0;
1303		_c = NULL;
1304
1305		/*
1306		leave generic data alone for now - just nuke it on clear()
1307	<	if (HasData("Comment")) delete [] (char*)GetData("Comment");
1308	<	_vdata.clear();
1307	>	if (HasData("Comment")) delete [] (char*)GetData("Comment");
1308	>	_vdata.clear();
1309		*/
1310
1311		if (Empty())
1312	<	return;
1312	>	return;
1313
1314		//if atoms present convert coords into array
1315		double *c = new double [NumAtoms()*3];
#	Line 1164 | Line 1319 | void OBMol::EndModify(bool nukePerceivedData)
1319		OBAtom *atom;
1320		vector<OBNodeBase*>::iterator j;
1321		for (idx=0,atom = BeginAtom(j);atom;atom = NextAtom(j),idx++)
1322	<	{
1322	>	{
1323		atom->SetIdx(idx+1);
1324		(atom->GetVector()).Get(&_c[idx*3]);
1325		atom->SetCoordPtr(&_c);
1326	<	}
1326	>	}
1327		_vconf.push_back(c);
1328
1329	<	//kekulize structure
1329	>	//kekulize structure
1330		SetAromaticPerceived();
1331		Kekulize();
1332		//kekulize();
#	Line 1186 | Line 1341 | void OBMol::EndModify(bool nukePerceivedData)
1341		//      bond->UnsetAromatic();
1342
1343		UnsetImplicitValencePerceived();
1344	<	}
1344	>	}
1345
1346	<	OBAtom *OBMol::CreateAtom(void)
1347	<	{
1346	>	OBAtom *OBMol::CreateAtom(void)
1347	>	{
1348		return new OBAtom;
1349	<	}
1349	>	}
1350
1351	<	OBBond *OBMol::CreateBond(void)
1352	<	{
1351	>	OBBond *OBMol::CreateBond(void)
1352	>	{
1353		return new OBBond;
1354	<	}
1354	>	}
1355
1356	<	void OBMol::DestroyAtom(OBNodeBase *atom)
1357	<	{
1356	>	void OBMol::DestroyAtom(OBNodeBase *atom)
1357	>	{
1358		if (atom)
1359	<	{
1359	>	{
1360		delete atom;
1361		atom = NULL;
1362	<	}
1363	<	}
1362	>	}
1363	>	}
1364
1365	<	void OBMol::DestroyBond(OBEdgeBase *bond)
1366	<	{
1365	>	void OBMol::DestroyBond(OBEdgeBase *bond)
1366	>	{
1367		if (bond)
1368	<	{
1368	>	{
1369		delete bond;
1370		bond = NULL;
1371	<	}
1372	<	}
1371	>	}
1372	>	}
1373
1374	<	//! \brief Get a new atom to add to a molecule
1375	<	//!
1376	<	//! Also checks bond_queue for any bonds that should be made to the new atom
1377	<	OBAtom *OBMol::NewAtom()
1378	<	{
1374	>	//! \brief Get a new atom to add to a molecule
1375	>	//!
1376	>	//! Also checks bond_queue for any bonds that should be made to the new atom
1377	>	OBAtom *OBMol::NewAtom()
1378	>	{
1379		BeginModify();
1380
1381		OBAtom *obatom = CreateAtom();
#	Line 1231 | Line 1386 | OBAtom *OBMol::NewAtom()
1386		#define OBAtomIncrement 100
1387
1388		if (_vatom.empty() || _natoms+1 >= (signed)_vatom.size())
1389	<	{
1389	>	{
1390		_vatom.resize(_natoms+OBAtomIncrement);
1391		vector<OBNodeBase*>::iterator j;
1392		for (j = _vatom.begin(),j+=(_natoms+1);j != _vatom.end();j++)
1393	<	*j = (OBNodeBase*)NULL;
1394	<	}
1393	>	*j = (OBNodeBase*)NULL;
1394	>	}
1395		#undef OBAtomIncrement
1396
1397		_vatom[_natoms] = obatom;
1398		_natoms++;
1399
1400		if (HasData(OBGenericDataType::VirtualBondData))
1401	<	{
1401	>	{
1402		/*add bonds that have been queued*/
1403		OBVirtualBond *vb;
1404		vector<OBGenericData*> verase;
1405		vector<OBGenericData*>::iterator i;
1406		for (i = BeginData();i != EndData();i++)
1407	<	if ((*i)->GetDataType() == OBGenericDataType::VirtualBondData)
1407	>	if ((*i)->GetDataType() == OBGenericDataType::VirtualBondData)
1408		{
1409	<	vb = (OBVirtualBond*)*i;
1410	<	if (vb->GetBgn() > _natoms || vb->GetEnd() > _natoms)
1411	<	continue;
1412	<	if (obatom->GetIdx() == static_cast<unsigned int>(vb->GetBgn())
1413	<	|| obatom->GetIdx() == static_cast<unsigned int>(vb->GetEnd()))
1409	>	vb = (OBVirtualBond*)*i;
1410	>	if (vb->GetBgn() > _natoms || vb->GetEnd() > _natoms)
1411	>	continue;
1412	>	if (obatom->GetIdx() == static_cast<unsigned int>(vb->GetBgn())
1413	>	|| obatom->GetIdx() == static_cast<unsigned int>(vb->GetEnd()))
1414		{
1415	<	AddBond(vb->GetBgn(),vb->GetEnd(),vb->GetOrder());
1416	<	verase.push_back(*i);
1415	>	AddBond(vb->GetBgn(),vb->GetEnd(),vb->GetOrder());
1416	>	verase.push_back(*i);
1417		}
1418		}
1419
1420		if (!verase.empty())
1421	<	DeleteData(verase);
1422	<	}
1421	>	DeleteData(verase);
1422	>	}
1423
1424		EndModify();
1425
1426		return(obatom);
1427	<	}
1427	>	}
1428
1429	<	OBResidue *OBMol::NewResidue()
1430	<	{
1429	>	OBResidue *OBMol::NewResidue()
1430	>	{
1431		OBResidue *obresidue = new OBResidue;
1432		obresidue->SetIdx(_residue.size());
1433		_residue.push_back(obresidue);
1434		return(obresidue);
1435	<	}
1435	>	}
1436
1437	<	//! \brief Add an atom to a molecule
1438	<	//!
1439	<	//! Also checks bond_queue for any bonds that should be made to the new atom
1440	<	bool OBMol::AddAtom(OBAtom &atom)
1441	<	{
1437	>	//! \brief Add an atom to a molecule
1438	>	//!
1439	>	//! Also checks bond_queue for any bonds that should be made to the new atom
1440	>	bool OBMol::AddAtom(OBAtom &atom)
1441	>	{
1442		BeginModify();
1443
1444		OBAtom *obatom = CreateAtom();
#	Line 1295 | Line 1450 | bool OBMol::AddAtom(OBAtom &atom)
1450		#define OBAtomIncrement 100
1451
1452		if (_vatom.empty() || _natoms+1 >= (signed)_vatom.size())
1453	<	{
1453	>	{
1454		_vatom.resize(_natoms+OBAtomIncrement);
1455		vector<OBNodeBase*>::iterator j;
1456		for (j = _vatom.begin(),j+=(_natoms+1);j != _vatom.end();j++)
1457	<	*j = (OBNodeBase*)NULL;
1458	<	}
1457	>	*j = (OBNodeBase*)NULL;
1458	>	}
1459		#undef OBAtomIncrement
1460
1461		_vatom[_natoms] = (OBNodeBase*)obatom;
1462		_natoms++;
1463
1464		if (HasData(OBGenericDataType::VirtualBondData))
1465	<	{
1465	>	{
1466		/*add bonds that have been queued*/
1467		OBVirtualBond *vb;
1468		vector<OBGenericData*> verase;
1469		vector<OBGenericData*>::iterator i;
1470		for (i = BeginData();i != EndData();i++)
1471	<	if ((*i)->GetDataType() == OBGenericDataType::VirtualBondData)
1471	>	if ((*i)->GetDataType() == OBGenericDataType::VirtualBondData)
1472		{
1473	<	vb = (OBVirtualBond*)*i;
1474	<	if (vb->GetBgn() > _natoms || vb->GetEnd() > _natoms)
1475	<	continue;
1476	<	if (obatom->GetIdx() == static_cast<unsigned int>(vb->GetBgn())
1477	<	|| obatom->GetIdx() == static_cast<unsigned int>(vb->GetEnd()))
1473	>	vb = (OBVirtualBond*)*i;
1474	>	if (vb->GetBgn() > _natoms || vb->GetEnd() > _natoms)
1475	>	continue;
1476	>	if (obatom->GetIdx() == static_cast<unsigned int>(vb->GetBgn())
1477	>	|| obatom->GetIdx() == static_cast<unsigned int>(vb->GetEnd()))
1478		{
1479	<	AddBond(vb->GetBgn(),vb->GetEnd(),vb->GetOrder());
1480	<	verase.push_back(*i);
1479	>	AddBond(vb->GetBgn(),vb->GetEnd(),vb->GetOrder());
1480	>	verase.push_back(*i);
1481		}
1482		}
1483
1484		if (!verase.empty())
1485	<	DeleteData(verase);
1486	<	}
1485	>	DeleteData(verase);
1486	>	}
1487
1488		EndModify();
1489
1490		return(true);
1491	<	}
1491	>	}
1492
1493	<	bool OBMol::InsertAtom(OBAtom &atom)
1494	<	{
1493	>	bool OBMol::InsertAtom(OBAtom &atom)
1494	>	{
1495		BeginModify();
1496
1497		OBAtom *obatom = CreateAtom();
#	Line 1348 | Line 1503 | bool OBMol::InsertAtom(OBAtom &atom)
1503		#define OBAtomIncrement 100
1504
1505		if (_vatom.empty() || _natoms+1 >= (signed)_vatom.size())
1506	<	{
1506	>	{
1507		_vatom.resize(_natoms+OBAtomIncrement);
1508		vector<OBNodeBase*>::iterator j;
1509		for (j = _vatom.begin(),j+=(_natoms+1);j != _vatom.end();j++)
1510	<	*j = (OBNodeBase*)NULL;
1511	<	}
1510	>	*j = (OBNodeBase*)NULL;
1511	>	}
1512		#undef OBAtomIncrement
1513
1514		_vatom[_natoms] = (OBNodeBase*)obatom;
1515		_natoms++;
1516
1517		if (HasData(OBGenericDataType::VirtualBondData))
1518	<	{
1518	>	{
1519		/*add bonds that have been queued*/
1520		OBVirtualBond *vb;
1521		vector<OBGenericData*> verase;
1522		vector<OBGenericData*>::iterator i;
1523		for (i = BeginData();i != EndData();i++)
1524	<	if ((*i)->GetDataType() == OBGenericDataType::VirtualBondData)
1524	>	if ((*i)->GetDataType() == OBGenericDataType::VirtualBondData)
1525		{
1526	<	vb = (OBVirtualBond*)*i;
1527	<	if (vb->GetBgn() > _natoms || vb->GetEnd() > _natoms)
1528	<	continue;
1529	<	if (obatom->GetIdx() == static_cast<unsigned int>(vb->GetBgn())
1530	<	|| obatom->GetIdx() == static_cast<unsigned int>(vb->GetEnd()))
1526	>	vb = (OBVirtualBond*)*i;
1527	>	if (vb->GetBgn() > _natoms || vb->GetEnd() > _natoms)
1528	>	continue;
1529	>	if (obatom->GetIdx() == static_cast<unsigned int>(vb->GetBgn())
1530	>	|| obatom->GetIdx() == static_cast<unsigned int>(vb->GetEnd()))
1531		{
1532	<	AddBond(vb->GetBgn(),vb->GetEnd(),vb->GetOrder());
1533	<	verase.push_back(*i);
1532	>	AddBond(vb->GetBgn(),vb->GetEnd(),vb->GetOrder());
1533	>	verase.push_back(*i);
1534		}
1535		}
1536
1537		if (!verase.empty())
1538	<	DeleteData(verase);
1539	<	}
1538	>	DeleteData(verase);
1539	>	}
1540
1541		EndModify();
1542
1543		return(true);
1544	<	}
1544	>	}
1545
1546	<	bool OBMol::AddResidue(OBResidue &residue)
1547	<	{
1546	>	bool OBMol::AddResidue(OBResidue &residue)
1547	>	{
1548		BeginModify();
1549
1550		OBResidue *obresidue = new OBResidue;
#	Line 1402 | Line 1557 | bool OBMol::AddResidue(OBResidue &residue)
1557		EndModify();
1558
1559		return(true);
1560	<	}
1560	>	}
1561
1562	<	bool OBMol::StripSalts()
1563	<	{
1562	>	bool OBMol::StripSalts()
1563	>	{
1564		vector<vector<int> > cfl;
1565		vector<vector<int> >::iterator i,max;
1566
1567		ContigFragList(cfl);
1568		if (cfl.empty() || cfl.size() == 1)
1569	<	return(false);
1569	>	return(false);
1570
1571	<	obErrorLog.ThrowError(__FUNCTION__,
1571	>	obErrorLog.ThrowError(__func__,
1572		"Ran OpenBabel::StripSalts", obAuditMsg);
1573
1574		max = cfl.begin();
1575		for (i = cfl.begin();i != cfl.end();i++)
1576	<	if ((*max).size() < (*i).size())
1577	<	max = i;
1576	>	if ((*max).size() < (*i).size())
1577	>	max = i;
1578
1579		vector<int>::iterator j;
1580		vector<OBNodeBase*> delatoms;
1581		for (i = cfl.begin();i != cfl.end();i++)
1582	<	if (i != max)
1583	<	for (j = (*i).begin();j != (*i).end();j++)
1584	<	delatoms.push_back(GetAtom(*j));
1582	>	if (i != max)
1583	>	for (j = (*i).begin();j != (*i).end();j++)
1584	>	delatoms.push_back(GetAtom(*j));
1585
1586		if (!delatoms.empty())
1587	<	{
1587	>	{
1588		int tmpflags = _flags & (~(OB_SSSR_MOL));
1589		BeginModify();
1590		vector<OBNodeBase*>::iterator k;
1591		for (k = delatoms.begin();k != delatoms.end();k++)
1592	<	DeleteAtom((OBAtom*)*k);
1592	>	DeleteAtom((OBAtom*)*k);
1593		EndModify();
1594		_flags = tmpflags;
1595	<	}
1595	>	}
1596
1597		return(true);
1598	<	}
1598	>	}
1599
1600	<	bool OBMol::DeleteNonPolarHydrogens()
1601	<	{
1600	>	bool OBMol::DeleteNonPolarHydrogens()
1601	>	{
1602		OBAtom *atom;
1603		vector<OBNodeBase*>::iterator i;
1604		vector<OBNodeBase*> delatoms;
1605
1606	<	obErrorLog.ThrowError(__FUNCTION__,
1606	>	obErrorLog.ThrowError(__func__,
1607		"Ran OpenBabel::DeleteHydrogens -- nonpolar",
1608		obAuditMsg);
1609
1610		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
1611	<	if (atom->IsNonPolarHydrogen())
1612	<	delatoms.push_back(atom);
1611	>	if (atom->IsNonPolarHydrogen())
1612	>	delatoms.push_back(atom);
1613
1614		if (delatoms.empty())
1615	<	return(true);
1615	>	return(true);
1616
1462	–	/*
1463	–
1464	–
1465	–	int idx1,idx2;
1466	–	vector<double*>::iterator j;
1467	–	for (idx1=0,idx2=0,atom = BeginAtom(i);atom;atom = NextAtom(i),idx1++)
1468	–	if (!atom->IsHydrogen())
1469	–	{
1470	–	for (j = _vconf.begin();j != _vconf.end();j++)
1471	–	memcpy((char*)&((*j)[idx2*3]),(char*)&((*j)[idx1*3]),sizeof(double)*3);
1472	–	idx2++;
1473	–	}
1474	–	*/
1475	–
1617		IncrementMod();
1618
1619		for (i = delatoms.begin();i != delatoms.end();i++)
1620	<	DeleteAtom((OBAtom *)*i);
1620	>	DeleteAtom((OBAtom *)*i);
1621
1622		DecrementMod();
1623
1624		return(true);
1625	<	}
1625	>	}
1626
1627	<	bool OBMol::DeleteHydrogens()
1628	<	{
1627	>	bool OBMol::DeleteHydrogens()
1628	>	{
1629		OBAtom *atom,*nbr;
1630		vector<OBNodeBase*>::iterator i;
1631		vector<OBNodeBase*> delatoms,va;
1632
1633	<	obErrorLog.ThrowError(__FUNCTION__,
1633	>	obErrorLog.ThrowError(__func__,
1634		"Ran OpenBabel::DeleteHydrogens", obAuditMsg);
1635
1636		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
1637	<	if (atom->IsHydrogen())
1638	<	delatoms.push_back(atom);
1637	>	if (atom->IsHydrogen())
1638	>	delatoms.push_back(atom);
1639
1640		if (delatoms.empty())
1641	<	return(true);
1641	>	return(true);
1642
1643		/* decide whether these flags need to be reset
1644	<	_flags &= (~(OB_ATOMTYPES_MOL));
1645	<	_flags &= (~(OB_HYBRID_MOL));
1646	<	_flags &= (~(OB_PCHARGE_MOL));
1647	<	_flags &= (~(OB_IMPVAL_MOL));
1644	>	_flags &= (~(OB_ATOMTYPES_MOL));
1645	>	_flags &= (~(OB_HYBRID_MOL));
1646	>	_flags &= (~(OB_PCHARGE_MOL));
1647	>	_flags &= (~(OB_IMPVAL_MOL));
1648		*/
1649
1650		//find bonds to delete
1651		vector<OBEdgeBase*> vdb;
1652		vector<OBEdgeBase*>::iterator j;
1653		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
1654	<	if (!atom->IsHydrogen())
1655	<	for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
1656	<	if (nbr->IsHydrogen())
1657	<	vdb.push_back(*j);
1654	>	if (!atom->IsHydrogen())
1655	>	for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
1656	>	if (nbr->IsHydrogen())
1657	>	vdb.push_back(*j);
1658
1659		IncrementMod();
1660		for (j = vdb.begin();j != vdb.end();j++)
1661	<	DeleteBond((OBBond *)*j); //delete bonds
1661	>	DeleteBond((OBBond *)*j); //delete bonds
1662		DecrementMod();
1663
1664		int idx1,idx2;
1665		vector<double*>::iterator k;
1666		for (idx1=0,idx2=0,atom = BeginAtom(i);atom;atom = NextAtom(i),idx1++)
1667	<	if (!atom->IsHydrogen())
1667	>	if (!atom->IsHydrogen())
1668		{
1669	<	//Update conformer coordinates
1670	<	for (k = _vconf.begin();k != _vconf.end();k++)
1671	<	memcpy((char*)&((*k)[idx2*3]),(char*)&((*k)[idx1*3]),sizeof(double)*3);
1669	>	//Update conformer coordinates
1670	>	for (k = _vconf.begin();k != _vconf.end();k++)
1671	>	memcpy((char*)&((*k)[idx2*3]),(char*)&((*k)[idx1*3]),sizeof(double)*3);
1672
1673	<	idx2++;
1674	<	va.push_back(atom);
1673	>	idx2++;
1674	>	va.push_back(atom);
1675		}
1676
1677		for (i = delatoms.begin();i != delatoms.end();i++)
1678	<	{
1678	>	{
1679		DestroyAtom(*i);
1680		_natoms--;
1681	<	}
1681	>	}
1682
1683		_vatom.clear();
1684		for (i = va.begin();i != va.end();i++)
1685	<	_vatom.push_back((OBNodeBase*)*i);
1685	>	_vatom.push_back((OBNodeBase*)*i);
1686
1687		//_atom = va;
1688		//_atom.resize(_atom.size()+1);
#	Line 1550 | Line 1691 | bool OBMol::DeleteHydrogens()
1691
1692		//reset all the indices to the atoms
1693		for (idx1=1,atom = BeginAtom(i);atom;atom = NextAtom(i),idx1++)
1694	<	atom->SetIdx(idx1);
1694	>	atom->SetIdx(idx1);
1695
1696		return(true);
1697	<	}
1697	>	}
1698
1699	<	bool OBMol::DeleteHydrogens(OBAtom *atom)
1700	<	//deletes all hydrogens attached to the atom passed to the function
1701	<	{
1699	>	bool OBMol::DeleteHydrogens(OBAtom *atom)
1700	>	//deletes all hydrogens attached to the atom passed to the function
1701	>	{
1702		OBAtom *nbr;
1703		vector<OBNodeBase*>::iterator i;
1704		vector<OBEdgeBase*>::iterator k;
1705		vector<OBNodeBase*> delatoms;
1706
1707		for (nbr = atom->BeginNbrAtom(k);nbr;nbr = atom->NextNbrAtom(k))
1708	<	if (nbr->IsHydrogen())
1709	<	delatoms.push_back(nbr);
1708	>	if (nbr->IsHydrogen())
1709	>	delatoms.push_back(nbr);
1710
1711		if (delatoms.empty())
1712	<	return(true);
1712	>	return(true);
1713
1714		IncrementMod();
1715		for (i = delatoms.begin();i != delatoms.end();i++)
1716	<	DeleteHydrogen((OBAtom*)*i);
1716	>	DeleteHydrogen((OBAtom*)*i);
1717		DecrementMod();
1718
1719		return(true);
1720	<	}
1720	>	}
1721
1722
1723	<	bool OBMol::DeleteHydrogen(OBAtom *atom)
1724	<	//deletes the hydrogen atom passed to the function
1725	<	{
1723	>	bool OBMol::DeleteHydrogen(OBAtom *atom)
1724	>	//deletes the hydrogen atom passed to the function
1725	>	{
1726		//find bonds to delete
1727		OBAtom *nbr;
1728		vector<OBEdgeBase*> vdb;
1729		vector<OBEdgeBase*>::iterator j;
1730		for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
1731	<	vdb.push_back(*j);
1731	>	vdb.push_back(*j);
1732
1733		IncrementMod();
1734		for (j = vdb.begin();j != vdb.end();j++)
1735	<	DeleteBond((OBBond*)*j); //delete bonds
1735	>	DeleteBond((OBBond*)*j); //delete bonds
1736		DecrementMod();
1737
1738		int idx;
1739		if (atom->GetIdx() != NumAtoms())
1740	<	{
1740	>	{
1741		idx = atom->GetCIdx();
1742		int size = NumAtoms()-atom->GetIdx();
1743		vector<double*>::iterator k;
1744		for (k = _vconf.begin();k != _vconf.end();k++)
1745	<	memmove((char*)&(*k)[idx],(char*)&(*k)[idx+3],sizeof(double)*3*size);
1745	>	memmove((char*)&(*k)[idx],(char*)&(*k)[idx+3],sizeof(double)*3*size);
1746
1747	<	}
1747	>	}
1748
1749		_vatom.erase(_vatom.begin()+(atom->GetIdx()-1));
1750		DestroyAtom(atom);
#	Line 1612 | Line 1753 | bool OBMol::DeleteHydrogen(OBAtom *atom)
1753		//reset all the indices to the atoms
1754		vector<OBNodeBase*>::iterator i;
1755		for (idx=1,atom = BeginAtom(i);atom;atom = NextAtom(i),idx++)
1756	<	atom->SetIdx(idx);
1756	>	atom->SetIdx(idx);
1757
1758		return(true);
1759	<	}
1759	>	}
1760
1761	<	bool OBMol::AddHydrogens(bool polaronly,bool correctForPH)
1762	<	{
1761	>	bool OBMol::AddHydrogens(bool polaronly,bool correctForPH)
1762	>	{
1763		if (!IsCorrectedForPH() && correctForPH)
1764	<	CorrectForPH();
1764	>	CorrectForPH();
1765
1766		if (HasHydrogensAdded())
1767	<	return(true);
1767	>	return(true);
1768		SetHydrogensAdded();
1769
1770		if (!polaronly)
1771	<	obErrorLog.ThrowError(__FUNCTION__,
1771	>	obErrorLog.ThrowError(__func__,
1772		"Ran OpenBabel::AddHydrogens", obAuditMsg);
1773		else
1774	<	obErrorLog.ThrowError(__FUNCTION__,
1775	<	"Ran OpenBabel::AddHydrogens -- polar only", obAuditMsg);
1774	>	obErrorLog.ThrowError(__func__,
1775	>	"Ran OpenBabel::AddHydrogens -- polar only", obAuditMsg);
1776
1777		//count up number of hydrogens to add
1778		OBAtom *atom,*h;
#	Line 1639 | Line 1780 | bool OBMol::AddHydrogens(bool polaronly,bool correctFo
1780		vector<pair<OBAtom*,int> > vhadd;
1781		vector<OBNodeBase*>::iterator i;
1782		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
1783	<	{
1783	>	{
1784		if (polaronly && !(atom->IsNitrogen() || atom->IsOxygen() ||
1785		atom->IsSulfur() || atom->IsPhosphorus()))
1786	<	continue;
1786	>	continue;
1787
1788		hcount = atom->GetImplicitValence() - atom->GetValence();
1789
#	Line 1654 | Line 1795 | bool OBMol::AddHydrogens(bool polaronly,bool correctFo
1795		hcount-=2;
1796
1797		if (hcount < 0)
1798	<	hcount = 0;
1798	>	hcount = 0;
1799		if (hcount)
1800	<	{
1800	>	{
1801		vhadd.push_back(pair<OBAtom*,int>(atom,hcount));
1802		count += hcount;
1803	<	}
1804	<	}
1803	>	}
1804	>	}
1805
1806		if (count == 0)
1807	<	return(true);
1807	>	return(true);
1808		bool hasCoords = HasNonZeroCoords();
1809
1810		//realloc memory in coordinate arrays for new hydrogens
1811		double *tmpf;
1812		vector<double*>::iterator j;
1813		for (j = _vconf.begin();j != _vconf.end();j++)
1814	<	{
1814	>	{
1815		tmpf = new double [(NumAtoms()+count)*3];
1816		memset(tmpf,'\0',sizeof(double)*(NumAtoms()+count)*3);
1817		if (hasCoords)
1818	<	memcpy(tmpf,(*j),sizeof(double)*NumAtoms()*3);
1818	>	memcpy(tmpf,(*j),sizeof(double)*NumAtoms()*3);
1819		delete []*j;
1820		*j = tmpf;
1821	<	}
1821	>	}
1822
1823		IncrementMod();
1824
#	Line 1688 | Line 1829 | bool OBMol::AddHydrogens(bool polaronly,bool correctFo
1829
1830
1831		for (k = vhadd.begin();k != vhadd.end();k++)
1832	<	{
1832	>	{
1833		atom = k->first;
1834		double bondlen = hbrad+etab.CorrectedBondRad(atom->GetAtomicNum(),atom->GetHyb());
1835		for (m = 0;m < k->second;m++)
1836	<	{
1836	>	{
1837		for (n = 0;n < NumConformers();n++)
1838	<	{
1838	>	{
1839		SetConformer(n);
1840		if (hasCoords)
1841	<	{
1841	>	{
1842		atom->GetNewBondVector(v,bondlen);
1843		_c[(NumAtoms())*3]   = v.x();
1844		_c[(NumAtoms())*3+1] = v.y();
1845		_c[(NumAtoms())*3+2] = v.z();
1846	<	}
1846	>	}
1847		else
1848	<	memset((char*)&_c[NumAtoms()*3],'\0',sizeof(double)*3);
1849	<	}
1848	>	memset((char*)&_c[NumAtoms()*3],'\0',sizeof(double)*3);
1849	>	}
1850		h = NewAtom();
1851		h->SetType("H");
1852		h->SetAtomicNum(1);
#	Line 1713 | Line 1854 | bool OBMol::AddHydrogens(bool polaronly,bool correctFo
1854		// copy parent atom residue to added hydrogen     REG 6/30/02
1855
1856		if (atom->HasResidue())
1857	<	{
1857	>	{
1858
1859		string aname;
1860
#	Line 1727 | Line 1868 | bool OBMol::AddHydrogens(bool polaronly,bool correctFo
1868
1869		atom->GetResidue()->SetAtomID(h,aname);
1870
1871	<	}
1871	>	}
1872
1873		AddBond(atom->GetIdx(),h->GetIdx(),1);
1874		h->SetCoordPtr(&_c);
1875	<	}
1876	<	}
1875	>	}
1876	>	}
1877
1878		DecrementMod();
1879		SetConformer(0);
#	Line 1741 | Line 1882 | bool OBMol::AddHydrogens(bool polaronly,bool correctFo
1882		_flags &= (~(OB_PCHARGE_MOL|OB_ATOMTYPES_MOL));
1883
1884		return(true);
1885	<	}
1885	>	}
1886
1887	<	bool OBMol::AddPolarHydrogens()
1888	<	{
1887	>	bool OBMol::AddPolarHydrogens()
1888	>	{
1889		return(AddHydrogens(true));
1890	<	}
1890	>	}
1891
1892	<	bool OBMol::AddHydrogens(OBAtom *atom)
1893	<	{
1892	>	bool OBMol::AddHydrogens(OBAtom *atom)
1893	>	{
1894		OBAtom *h;
1895
1896		//count up number of hydrogens to add
#	Line 1758 | Line 1899 | bool OBMol::AddHydrogens(OBAtom *atom)
1899
1900		hcount = atom->GetImplicitValence() - atom->GetValence();
1901
1902	<	//Jan 05 Implicit valency now left alone; use spin multiplicity for implicit Hs
1903	<	int mult = atom->GetSpinMultiplicity();
1904	<	if(mult==2) //radical
1905	<	hcount-=1;
1906	<	else if(mult==1 || mult==3) //carbene
1907	<	hcount-=2;
1902	>	//Jan 05 Implicit valency now left alone; use spin multiplicity for implicit Hs
1903	>	int mult = atom->GetSpinMultiplicity();
1904	>	if(mult==2) //radical
1905	>	hcount-=1;
1906	>	else if(mult==1 || mult==3) //carbene
1907	>	hcount-=2;
1908
1909	<	if (hcount < 0)
1910	<	hcount = 0;
1909	>	if (hcount < 0)
1910	>	hcount = 0;
1911		if (hcount)
1912	<	{
1912	>	{
1913		vhadd.push_back(pair<OBAtom*,int>(atom,hcount));
1914		count += hcount;
1915	<	}
1915	>	}
1916
1917		if (count == 0)
1918	<	return(true);
1918	>	return(true);
1919
1920		//realloc memory in coordinate arrays for new hydroges
1921		double *tmpf;
1922		vector<double*>::iterator j;
1923		for (j = _vconf.begin();j != _vconf.end();j++)
1924	<	{
1924	>	{
1925		tmpf = new double [(NumAtoms()+count)*3+10];
1926		memcpy(tmpf,(*j),sizeof(double)*NumAtoms()*3);
1927		delete []*j;
1928		*j = tmpf;
1929	<	}
1929	>	}
1930
1931		IncrementMod();
1932
#	Line 1795 | Line 1936 | bool OBMol::AddHydrogens(OBAtom *atom)
1936		double hbrad = etab.CorrectedBondRad(1,0);
1937
1938		for (k = vhadd.begin();k != vhadd.end();k++)
1939	<	{
1939	>	{
1940		atom = k->first;
1941		double bondlen = hbrad+etab.CorrectedBondRad(atom->GetAtomicNum(),atom->GetHyb());
1942		for (m = 0;m < k->second;m++)
1943	<	{
1943	>	{
1944		for (n = 0;n < NumConformers();n++)
1945	<	{
1945	>	{
1946		SetConformer(n);
1947		atom->GetNewBondVector(v,bondlen);
1948		_c[(NumAtoms())*3]   = v.x();
1949		_c[(NumAtoms())*3+1] = v.y();
1950		_c[(NumAtoms())*3+2] = v.z();
1951	<	}
1951	>	}
1952		h = NewAtom();
1953		h->SetType("H");
1954		h->SetAtomicNum(1);
1955		AddBond(atom->GetIdx(),h->GetIdx(),1);
1956		h->SetCoordPtr(&_c);
1957	<	}
1958	<	}
1957	>	}
1958	>	}
1959
1960		DecrementMod();
1961		SetConformer(0);
#	Line 1823 | Line 1964 | bool OBMol::AddHydrogens(OBAtom *atom)
1964		//_flags &= (~(OB_PCHARGE_MOL|OB_ATOMTYPES_MOL));
1965
1966		return(true);
1967	<	}
1967	>	}
1968
1969	<	bool OBMol::CorrectForPH()
1970	<	{
1969	>	bool OBMol::CorrectForPH()
1970	>	{
1971		if (IsCorrectedForPH())
1972	<	return(true);
1972	>	return(true);
1973		phmodel.CorrectForPH(*this);
1974
1975	<	obErrorLog.ThrowError(__FUNCTION__,
1975	>	obErrorLog.ThrowError(__func__,
1976		"Ran OpenBabel::CorrectForPH", obAuditMsg);
1977
1978		return(true);
1979	<	}
1979	>	}
1980
1981	<	//! \brief set spin multiplicity for H-deficient atoms
1982	<	bool OBMol::AssignSpinMultiplicity()
1983	<	{
1981	>	//! \brief set spin multiplicity for H-deficient atoms
1982	>	bool OBMol::AssignSpinMultiplicity()
1983	>	{
1984		if (HasSpinMultiplicityAssigned())
1985	<	return(true);
1985	>	return(true);
1986
1987		SetSpinMultiplicityAssigned();
1988
1989	<	obErrorLog.ThrowError(__FUNCTION__,
1989	>	obErrorLog.ThrowError(__func__,
1990		"Ran OpenBabel::AssignSpinMultiplicity",
1991		obAuditMsg);
1992
#	Line 1857 | Line 1998 | bool OBMol::AssignSpinMultiplicity()
1998		//Any discrepancy with the expected atom valency is because it is a radical of some sort
1999		//Also adjust the ImplicitValence for radical atoms
2000		for (atom = BeginAtom(k);atom;atom = NextAtom(k))
2001	<	{
2001	>	{
2002
2003	<	if (!atom->IsHydrogen() && atom->ExplicitHydrogenCount()!=0)
2004	<	{
2003	>	if (!atom->IsHydrogen() && atom->ExplicitHydrogenCount()!=0)
2004	>	{
2005		diff=atom->GetImplicitValence() - (atom->GetHvyValence() + atom->ExplicitHydrogenCount());
2006		if (diff)
2007	<	atom->SetSpinMultiplicity(diff+1);//radicals =2; all carbenes =3
2008	<	}
2007	>	atom->SetSpinMultiplicity(diff+1);//radicals =2; all carbenes =3
2008	>	}
2009
2010	<	//Jan05        mult=atom->GetSpinMultiplicity();
2011	<	//        if(mult) //radical or carbene
2012	<	//            atom->DecrementImplicitValence();
2013	<	//        if(mult==1 || mult==3) //e.g.singlet or triplet carbene
2014	<	//            atom->DecrementImplicitValence();
2015	<	}
2010	>	//Jan05        mult=atom->GetSpinMultiplicity();
2011	>	//        if(mult) //radical or carbene
2012	>	//            atom->DecrementImplicitValence();
2013	>	//        if(mult==1 || mult==3) //e.g.singlet or triplet carbene
2014	>	//            atom->DecrementImplicitValence();
2015	>	}
2016		//end CM
2017
2018		vector<OBNodeBase*>::iterator i;
2019		unsigned int spin = 1;
2020
2021		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
2022	<	{
2022	>	{
2023		if (atom->GetSpinMultiplicity() > 1)
2024	<	spin += atom->GetSpinMultiplicity() - 1;
2025	<	}
2024	>	spin += atom->GetSpinMultiplicity() - 1;
2025	>	}
2026		_totalSpin = spin;
2027
2028		return (true);
2029	<	}
2029	>	}
2030
2031
2032	<	// Not used anywhere internally -- likely predates OBBase code
2033	<	// static void ResetVisit(OBMol &mol,vector<int> &visit,int depth)
2034	<	// {
2035	<	//     OBBond *bond;
2036	<	//     vector<OBEdgeBase*>::iterator i;
2032	>	// Not used anywhere internally -- likely predates OBBase code
2033	>	// static void ResetVisit(OBMol &mol,vector<int> &visit,int depth)
2034	>	// {
2035	>	//     OBBond *bond;
2036	>	//     vector<OBEdgeBase*>::iterator i;
2037
2038	<	//     for (bond = mol.BeginBond(i);bond;bond = mol.NextBond(i))
2039	<	//         if (bond->IsAromatic() && visit[bond->GetIdx()] >= depth)
2040	<	//             visit[bond->GetIdx()] = 0;
2041	<	// }
2038	>	//     for (bond = mol.BeginBond(i);bond;bond = mol.NextBond(i))
2039	>	//         if (bond->IsAromatic() && visit[bond->GetIdx()] >= depth)
2040	>	//             visit[bond->GetIdx()] = 0;
2041	>	// }
2042
2043	<	static int ValenceSum(OBAtom *atom)
2044	<	{
2043	>	static int ValenceSum(OBAtom *atom)
2044	>	{
2045		int count = atom->GetImplicitValence();
2046
2047		OBBond *bond;
2048		vector<OBEdgeBase*>::iterator i;
2049		for (bond = atom->BeginBond(i);bond;bond = atom->NextBond(i))
2050	<	if (bond->IsKDouble())
2051	<	count++;
2050	>	if (bond->IsKDouble())
2051	>	count++;
2052
2053		return(count);
2054	<	}
2054	>	}
2055
2056	<	static bool KekulePropagate(OBAtom *atom,vector<int> &visit,vector<int> &ival,int depth)
2057	<	{
2056	>	static bool KekulePropagate(OBAtom *atom,vector<int> &visit,vector<int> &ival,int depth)
2057	>	{
2058		int count = 0;
2059		OBBond *bond;
2060		vector<OBEdgeBase*>::iterator i;
2061		for (bond = atom->BeginBond(i);bond;bond = atom->NextBond(i))
2062	<	if (!visit[bond->GetIdx()])
2063	<	count++;
2062	>	if (!visit[bond->GetIdx()])
2063	>	count++;
2064
2065		if (!count)
2066	<	return(ValenceSum(atom) == ival[atom->GetIdx()]);
2066	>	return(ValenceSum(atom) == ival[atom->GetIdx()]);
2067
2068		bool result = true;
2069		OBAtom *nbr;
2070
2071		if (ValenceSum(atom) >= ival[atom->GetIdx()])
2072	<	{
2072	>	{
2073		for (nbr = atom->BeginNbrAtom(i);nbr;nbr = atom->NextNbrAtom(i))
2074	<	if (nbr->IsAromatic() && !visit[(*i)->GetIdx()])
2074	>	if (nbr->IsAromatic() && !visit[(*i)->GetIdx()])
2075		{
2076	<	visit[(*i)->GetIdx()] = depth;
2077	<	((OBBond*)*i)->SetKSingle();
2078	<	result = KekulePropagate(nbr,visit,ival,depth);
2079	<	if (result)
2080	<	break;
2081	<	//            if (!result) break;
2076	>	visit[(*i)->GetIdx()] = depth;
2077	>	((OBBond*)*i)->SetKSingle();
2078	>	result = KekulePropagate(nbr,visit,ival,depth);
2079	>	if (result)
2080	>	break;
2081	>	//            if (!result) break;
2082		}
2083	<	}
2083	>	}
2084		else if (count == 1)
2085	<	for (nbr = atom->BeginNbrAtom(i);nbr;nbr = atom->NextNbrAtom(i))
2086	<	if (nbr->IsAromatic() && !visit[(*i)->GetIdx()])
2087	<	{
2088	<	visit[(*i)->GetIdx()] = depth;
2089	<	((OBBond*)*i)->SetKDouble();
2090	<	result = KekulePropagate(nbr,visit,ival,depth);
2091	<	//break;
2092	<	if (result)
2093	<	break;
2094	<	}
2085	>	for (nbr = atom->BeginNbrAtom(i);nbr;nbr = atom->NextNbrAtom(i))
2086	>	if (nbr->IsAromatic() && !visit[(*i)->GetIdx()])
2087	>	{
2088	>	visit[(*i)->GetIdx()] = depth;
2089	>	((OBBond*)*i)->SetKDouble();
2090	>	result = KekulePropagate(nbr,visit,ival,depth);
2091	>	//break;
2092	>	if (result)
2093	>	break;
2094	>	}
2095		return(result);
2096	<	}
2096	>	}
2097
2098	<	int GetCurrentValence(OBAtom *atom)
2099	<	{
2098	>	int GetCurrentValence(OBAtom *atom)
2099	>	{
2100		int count = atom->GetImplicitValence();
2101
2102		OBBond *bond;
2103		vector<OBEdgeBase*>::iterator i;
2104		for (bond = atom->BeginBond(i);bond;bond = atom->NextBond(i))
2105	<	{
2105	>	{
2106		if (bond->IsKDouble())
2107	<	count++;
2107	>	count++;
2108		else if (bond->IsKTriple())
2109	<	count += 2;
2109	>	count += 2;
2110		//      else if (bond->IsSingle()) count++;
2111		//      else if (bond->IsDouble()) count += 2;
2112		//      else if (bond->IsTriple()) count += 3;
2113	<	}
2113	>	}
2114		return(count);
2115	<	}
2115	>	}
2116
2117	<	bool ExpandKekule(OBMol &mol, vector<OBNodeBase*> &va,
2118	<	vector<OBNodeBase*>::iterator i,
2119	<	vector<int> &maxv,bool secondpass)
2120	<	{
2117	>	bool ExpandKekule(OBMol &mol, vector<OBNodeBase*> &va,
2118	>	vector<OBNodeBase*>::iterator i,
2119	>	vector<int> &maxv,bool secondpass)
2120	>	{
2121		if (i == va.end())
2122	<	{
2122	>	{
2123		//check to see that the ideal valence has been achieved for all atoms
2124		vector<OBNodeBase*>::iterator j;
2125		for (j = va.begin();j != va.end();j++)
2126	<	{
2126	>	{
2127		//let erroneously aromatic carboxylates pass
2128		if (((OBAtom*)*j)->IsOxygen() && ((OBAtom*)*j)->GetValence() == 1)
2129	<	continue;
2129	>	continue;
2130		if (GetCurrentValence((OBAtom*)*j) != maxv[(*j)->GetIdx()])
2131	<	{
2131	>	{
2132		//            cout << " ExpandKekule atom: " << ((OBAtom*)*j)->GetIdx()
2133		//                 << " valence is " << (GetCurrentValence((OBAtom*)*j))
2134		//                 << " should be " << maxv[(*j)->GetIdx()] << endl;
2135		return(false);
2136	<	}
2137	<	}
2136	>	}
2137	>	}
2138		return(true);
2139	<	}
2139	>	}
2140
2141		//jump to next atom in list if current atom doesn't have any attached
2142		//aromatic bonds
#	Line 2004 | Line 2145 | bool ExpandKekule(OBMol &mol, vector<OBNodeBase*> &va,
2145		vector<OBEdgeBase*>::iterator j;
2146		bool done = true;
2147		for (bond = atom->BeginBond(j);bond;bond = atom->NextBond(j))
2148	<	if (bond->GetBO() == 5)
2148	>	if (bond->GetBO() == 5)
2149		{
2150	<	done = false;
2151	<	break;
2150	>	done = false;
2151	>	break;
2152		}
2153		if (done)
2154	<	return(ExpandKekule(mol,va,i+1,maxv,secondpass));
2154	>	return(ExpandKekule(mol,va,i+1,maxv,secondpass));
2155
2156		//store list of attached aromatic atoms
2157		OBAtom *nbr;
2158		vector<OBEdgeBase*> vb;
2159		for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
2160	<	if ((*j)->GetBO() == 5)
2160	>	if ((*j)->GetBO() == 5)
2161		{
2162	<	vb.push_back(*j);
2163	<	((OBBond *)*j)->SetBO(1);
2164	<	((OBBond *)*j)->SetKSingle();
2162	>	vb.push_back(*j);
2163	>	((OBBond *)*j)->SetBO(1);
2164	>	((OBBond *)*j)->SetKSingle();
2165		}
2166
2167		//try setting a double bond
2168		if (GetCurrentValence(atom) < maxv[atom->GetIdx()])
2169	<	{
2169	>	{
2170		for (j = vb.begin();j != vb.end();j++)
2171	<	{
2171	>	{
2172		nbr = ((OBBond *)*j)->GetNbrAtom(atom);
2173		if (GetCurrentValence(nbr) <= maxv[nbr->GetIdx()])
2174	<	{
2174	>	{
2175		((OBBond*)*j)->SetKDouble();
2176		((OBBond*)*j)->SetBO(2);
2177		if (ExpandKekule(mol,va,i+1,maxv,secondpass))
2178	<	return(true);
2178	>	return(true);
2179		((OBBond*)*j)->SetKSingle();
2180		((OBBond*)*j)->SetBO(1);
2181	<	}
2182	<	}
2181	>	}
2182	>	}
2183
2184		if (secondpass && atom->IsNitrogen() && atom->GetFormalCharge() == 0 &&
2185	<	atom->GetImplicitValence() == 2)
2186	<	{
2185	>	atom->GetImplicitValence() == 2)
2186	>	{
2187		atom->IncrementImplicitValence();
2188		if (ExpandKekule(mol,va,i+1,maxv,secondpass))
2189	<	return(true);
2189	>	return(true);
2190		atom->DecrementImplicitValence();
2191	<	}
2192	<	}
2191	>	}
2192	>	}
2193		else  //full valence - no double bond to set
2194	<	{
2194	>	{
2195		if (ExpandKekule(mol,va,i+1,maxv,secondpass))
2196	<	return(true);
2196	>	return(true);
2197
2198		bool trycharge = false;
2199		if (secondpass && atom->GetFormalCharge() == 0)
2200	<	{
2200	>	{
2201		if (atom->IsNitrogen() && atom->GetHvyValence() == 3)
2202	<	trycharge = true;
2202	>	trycharge = true;
2203		if (atom->IsOxygen() && atom->GetHvyValence() == 2)
2204	<	trycharge = true;
2204	>	trycharge = true;
2205		if (atom->IsSulfur() && atom->GetHvyValence() == 2)
2206	<	trycharge = true;
2207	<	}
2206	>	trycharge = true;
2207	>	}
2208
2209		if (trycharge) //attempt to charge up O,N,S to make a valid kekule form
2210	<	{
2210	>	{
2211		maxv[atom->GetIdx()]++;
2212		atom->SetFormalCharge(1);
2213		for (j = vb.begin();j != vb.end();j++)
2214	<	{
2214	>	{
2215		nbr = ((OBBond*)*j)->GetNbrAtom(atom);
2216		if (GetCurrentValence(nbr) <= maxv[nbr->GetIdx()])
2217	<	{
2217	>	{
2218		((OBBond*)*j)->SetKDouble();
2219		((OBBond*)*j)->SetBO(2);
2220		if (ExpandKekule(mol,va,i+1,maxv,secondpass))
2221	<	return(true);
2221	>	return(true);
2222		((OBBond*)*j)->SetKSingle();
2223		((OBBond*)*j)->SetBO(1);
2224	<	}
2225	<	}
2224	>	}
2225	>	}
2226		maxv[atom->GetIdx()]--;
2227		atom->SetFormalCharge(0);
2228	<	}
2228	>	}
2229
2230		if (secondpass && atom->IsNitrogen() && atom->GetFormalCharge() == 0 &&
2231	<	atom->GetImplicitValence() == 2) //try protonating the nitrogen
2232	<	{
2231	>	atom->GetImplicitValence() == 2) //try protonating the nitrogen
2232	>	{
2233		atom->IncrementImplicitValence();
2234		if (ExpandKekule(mol,va,i+1,maxv,secondpass))
2235	<	return(true);
2235	>	return(true);
2236		atom->DecrementImplicitValence();
2237	<	}
2238	<	}
2237	>	}
2238	>	}
2239
2240		//failed to find a valid solution - reset attached bonds
2241		for (j = vb.begin();j != vb.end();j++)
2242	<	{
2242	>	{
2243		((OBBond*)*j)->SetKSingle();
2244		((OBBond*)*j)->SetBO(5);
2245	<	}
2245	>	}
2246
2247		return(false);
2248	<	}
2248	>	}
2249
2250	<	void CorrectBadResonanceForm(OBMol &mol)
2251	<	{
2250	>	void CorrectBadResonanceForm(OBMol &mol)
2251	>	{
2252		string s;
2253		OBBond *b1,*b2,*b3;
2254		OBAtom *a1,*a2,*a3,*a4;
2255		vector<vector<int> > mlist;
2256		vector<vector<int> >::iterator i;
2257
2258	<	obErrorLog.ThrowError(__FUNCTION__,
2258	>	obErrorLog.ThrowError(__func__,
2259		"Ran OpenBabel::CorrectBadResonanceForm", obAuditMsg);
2260
2261		OBSmartsPattern acid;
#	Line 2122 | Line 2263 | void CorrectBadResonanceForm(OBMol &mol)
2263
2264		//carboxylic acid
2265		if (acid.Match(mol))
2266	<	{
2266	>	{
2267		mlist = acid.GetUMapList();
2268		for (i = mlist.begin();i != mlist.end();i++)
2269	<	{
2269	>	{
2270		a1 = mol.GetAtom((*i)[0]);
2271		a2 = mol.GetAtom((*i)[1]);
2272		a3 = mol.GetAtom((*i)[2]);
2273		b1 = a2->GetBond(a1);
2274		b2 = a2->GetBond(a3);
2275		if (!b1 || !b2)
2276	<	continue;
2276	>	continue;
2277		b1->SetKDouble();
2278		b2->SetKSingle();
2279	<	}
2280	<	}
2279	>	}
2280	>	}
2281
2282		//phosphonic acid
2283		OBSmartsPattern phosphate;
2284		phosphate.Init("[p]([oD1])([oD1])([oD1])[#6,#8]");
2285		if (phosphate.Match(mol))
2286	<	{
2286	>	{
2287		mlist = phosphate.GetUMapList();
2288		for (i = mlist.begin();i != mlist.end();i++)
2289	<	{
2289	>	{
2290		a1 = mol.GetAtom((*i)[0]);
2291		a2 = mol.GetAtom((*i)[1]);
2292		a3 = mol.GetAtom((*i)[2]);
#	Line 2155 | Line 2296 | void CorrectBadResonanceForm(OBMol &mol)
2296		b3 = a1->GetBond(a4);
2297
2298		if (!b1 || !b2 || !b3)
2299	<	continue;
2299	>	continue;
2300		b1->SetKDouble();
2301		b2->SetKSingle();
2302		b3->SetKSingle();
2303	<	}
2304	<	}
2303	>	}
2304	>	}
2305
2306		//amidene and guanidine
2307		OBSmartsPattern amidene;
2308		amidene.Init("[nD1]c([nD1])*");
2309		if (amidene.Match(mol))
2310	<	{
2310	>	{
2311		mlist = amidene.GetUMapList();
2312		for (i = mlist.begin();i != mlist.end();i++)
2313	<	{
2313	>	{
2314		a1 = mol.GetAtom((*i)[0]);
2315		a2 = mol.GetAtom((*i)[1]);
2316		a3 = mol.GetAtom((*i)[2]);
2317		b1 = a2->GetBond(a1);
2318		b2 = a2->GetBond(a3);
2319		if (!b1 || !b2)
2320	<	continue;
2320	>	continue;
2321		b1->SetKDouble();
2322		b2->SetKSingle();
2323	<	}
2324	<	}
2325	<	}
2323	>	}
2324	>	}
2325	>	}
2326
2327	<	bool OBMol::PerceiveKekuleBonds()
2328	<	{
2327	>	bool OBMol::PerceiveKekuleBonds()
2328	>	{
2329		if (HasKekulePerceived())
2330	<	return(true);
2330	>	return(true);
2331		SetKekulePerceived();
2332
2333		OBBond *bond;
#	Line 2198 | Line 2339 | bool OBMol::PerceiveKekuleBonds()
2339		vector<bool> varo;
2340		varo.resize(NumAtoms()+1,false);
2341		for (bond = BeginBond(i);bond;bond = NextBond(i))
2342	<	switch (bond->GetBO())
2342	>	switch (bond->GetBO())
2343		{
2344		case 2:
2345	<	bond->SetKDouble();
2346	<	break;
2345	>	bond->SetKDouble();
2346	>	break;
2347		case 3:
2348	<	bond->SetKTriple();
2349	<	break;
2348	>	bond->SetKTriple();
2349	>	break;
2350		case 5:
2351
2352	<	bond->SetKSingle();
2353	<	if (bond->IsInRing())
2352	>	bond->SetKSingle();
2353	>	if (bond->IsInRing())
2354		{
2355	<	varo[bond->GetBeginAtomIdx()] = true;
2356	<	varo[bond->GetEndAtomIdx()]   = true;
2357	<	done = false;
2355	>	varo[bond->GetBeginAtomIdx()] = true;
2356	>	varo[bond->GetEndAtomIdx()]   = true;
2357	>	done = false;
2358		}
2359	<	else
2360	<	badResonanceForm = true;
2359	>	else
2360	>	badResonanceForm = true;
2361
2362	<	break;
2362	>	break;
2363
2364		default:
2365	<	bond->SetKSingle();
2366	<	break;
2365	>	bond->SetKSingle();
2366	>	break;
2367		}
2368
2369		if (badResonanceForm)
2370	<	CorrectBadResonanceForm(*this);
2370	>	CorrectBadResonanceForm(*this);
2371
2372		if (done)
2373	<	return(true);
2373	>	return(true);
2374
2375		//set the maximum valence for each aromatic atom
2376		OBAtom *atom,*nbr;
#	Line 2238 | Line 2379 | bool OBMol::PerceiveKekuleBonds()
2379		maxv.resize(NumAtoms()+1);
2380
2381		for (atom = BeginAtom(j);atom;atom = NextAtom(j))
2382	<	if (varo[atom->GetIdx()])
2382	>	if (varo[atom->GetIdx()])
2383		{
2384	<	switch (atom->GetAtomicNum())
2384	>	switch (atom->GetAtomicNum())
2385		{
2386		case 6:
2387	<	maxv[atom->GetIdx()] = 4;
2388	<	break;
2387	>	maxv[atom->GetIdx()] = 4;
2388	>	break;
2389		case 8:
2390		case 16:
2391		case 34:
2392		case 52:
2393	<	maxv[atom->GetIdx()] = 2;
2394	<	break;
2393	>	maxv[atom->GetIdx()] = 2;
2394	>	break;
2395		case 7:
2396		case 15:
2397		case 33:
2398	<	maxv[atom->GetIdx()] = 3;
2399	<	break;
2398	>	maxv[atom->GetIdx()] = 3;
2399	>	break;
2400		}
2401	<	//correct valence for formal charges
2402	<	if (atom->IsCarbon())
2403	<	maxv[atom->GetIdx()] -= abs(atom->GetFormalCharge());
2404	<	else
2405	<	maxv[atom->GetIdx()] += atom->GetFormalCharge();
2401	>	//correct valence for formal charges
2402	>	if (atom->IsCarbon())
2403	>	maxv[atom->GetIdx()] -= abs(atom->GetFormalCharge());
2404	>	else
2405	>	maxv[atom->GetIdx()] += atom->GetFormalCharge();
2406
2407	<	if (atom->IsNitrogen() || atom->IsSulfur())
2408	<	for (nbr = atom->BeginNbrAtom(i);nbr;nbr = atom->NextNbrAtom(i))
2409	<	if (nbr->IsOxygen() && (*i)->GetBO() == 2)
2410	<	maxv[atom->GetIdx()] += 2;
2407	>	if (atom->IsNitrogen() || atom->IsSulfur())
2408	>	for (nbr = atom->BeginNbrAtom(i);nbr;nbr = atom->NextNbrAtom(i))
2409	>	if (nbr->IsOxygen() && (*i)->GetBO() == 2)
2410	>	maxv[atom->GetIdx()] += 2;
2411		}
2412
2413		bool result = true;
#	Line 2274 | Line 2415 | bool OBMol::PerceiveKekuleBonds()
2415		used.resize(NumAtoms()+1);
2416		vector<OBNodeBase*> va,curr,next;
2417		for (atom = BeginAtom(j);atom;atom = NextAtom(j))
2418	<	if (varo[atom->GetIdx()] && !used[atom->GetIdx()])
2418	>	if (varo[atom->GetIdx()] && !used[atom->GetIdx()])
2419		{
2420	<	va.clear();
2421	<	va.push_back(atom);
2422	<	curr.clear();
2423	<	curr.push_back(atom);
2424	<	used[atom->GetIdx()] = true;
2420	>	va.clear();
2421	>	va.push_back(atom);
2422	>	curr.clear();
2423	>	curr.push_back(atom);
2424	>	used[atom->GetIdx()] = true;
2425
2426	<	for (;!curr.empty();)
2426	>	for (;!curr.empty();)
2427		{
2428	<	next.clear();
2429	<	for (k = curr.begin();k != curr.end();k++)
2430	<	for (nbr = ((OBAtom*)*k)->BeginNbrAtom(i);nbr;nbr = ((OBAtom*)*k)->NextNbrAtom(i))
2431	<	if (varo[nbr->GetIdx()] && !used[nbr->GetIdx()])
2432	<	{
2433	<	used[nbr->GetIdx()] = true;
2434	<	next.push_back(nbr);
2435	<	va.push_back(nbr);
2436	<	}
2437	<	curr = next;
2428	>	next.clear();
2429	>	for (k = curr.begin();k != curr.end();k++)
2430	>	for (nbr = ((OBAtom*)*k)->BeginNbrAtom(i);nbr;nbr = ((OBAtom*)*k)->NextNbrAtom(i))
2431	>	if (varo[nbr->GetIdx()] && !used[nbr->GetIdx()])
2432	>	{
2433	>	used[nbr->GetIdx()] = true;
2434	>	next.push_back(nbr);
2435	>	va.push_back(nbr);
2436	>	}
2437	>	curr = next;
2438		}
2439
2440	<	//try it first without protonating aromatic nitrogens
2441	<	if (!ExpandKekule(*this,va,va.begin(),maxv,false) &&
2442	<	!ExpandKekule(*this,va,va.begin(),maxv,true))
2440	>	//try it first without protonating aromatic nitrogens
2441	>	if (!ExpandKekule(*this,va,va.begin(),maxv,false) &&
2442	>	!ExpandKekule(*this,va,va.begin(),maxv,true))
2443		{
2444	<	result = false;
2445	<	//            cerr << " Died on atom " << atom->GetIdx() << endl;
2444	>	result = false;
2445	>	//              cerr << " Died on atom " << atom->GetIdx() << endl;
2446		}
2447		}
2448
2449		if (!result)
2450	<	{
2451	<	//        cerr << "Kekulization Error = " << GetTitle() << endl;
2450	>	{
2451	>	//        cerr << "Kekulization Error = " << GetTitle() << endl;
2452		//exit(0);
2453	<	}
2453	>	}
2454
2455		return(result);
2456	<	}
2456	>	}
2457
2458	<	bool OBMol::Kekulize()
2459	<	{
2458	>	bool OBMol::Kekulize()
2459	>	{
2460		OBBond *bond;
2461		vector<OBEdgeBase*>::iterator i;
2462		// Not quite sure why this is here -GRH 2003
2463		//  if (NumAtoms() > 255) return(false);
2464
2465	<	obErrorLog.ThrowError(__FUNCTION__,
2465	>	obErrorLog.ThrowError(__func__,
2466		"Ran OpenBabel::Kekulize", obAuditMsg);
2467
2468		for (bond = BeginBond(i);bond;bond = NextBond(i))
2469	<	if (bond->IsKSingle())
2470	<	bond->SetBO(1);
2471	<	else if (bond->IsKDouble())
2472	<	bond->SetBO(2);
2473	<	else if (bond->IsKTriple())
2474	<	bond->SetBO(3);
2469	>	if (bond->IsKSingle())
2470	>	bond->SetBO(1);
2471	>	else if (bond->IsKDouble())
2472	>	bond->SetBO(2);
2473	>	else if (bond->IsKTriple())
2474	>	bond->SetBO(3);
2475
2476		return(true);
2477	<	}
2477	>	}
2478
2479	<	bool OBMol::DeleteAtom(OBAtom *atom)
2480	<	{
2479	>	bool OBMol::DeleteAtom(OBAtom *atom)
2480	>	{
2481		if (atom->IsHydrogen())
2482	<	return(DeleteHydrogen(atom));
2482	>	return(DeleteHydrogen(atom));
2483
2484		BeginModify();
2485		//don't need to do anything with coordinates b/c
#	Line 2349 | Line 2490 | bool OBMol::DeleteAtom(OBAtom *atom)
2490		vector<OBEdgeBase*> vdb;
2491		vector<OBEdgeBase*>::iterator j;
2492		for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j))
2493	<	vdb.push_back(*j);
2493	>	vdb.push_back(*j);
2494
2495		for (j = vdb.begin();j != vdb.end();j++)
2496	<	DeleteBond((OBBond *)*j); //delete bonds
2496	>	DeleteBond((OBBond *)*j); //delete bonds
2497
2498		_vatom.erase(_vatom.begin()+(atom->GetIdx()-1));
2499		DestroyAtom(atom);
#	Line 2362 | Line 2503 | bool OBMol::DeleteAtom(OBAtom *atom)
2503		int idx;
2504		vector<OBNodeBase*>::iterator i;
2505		for (idx=1,atom = BeginAtom(i);atom;atom = NextAtom(i),idx++)
2506	<	atom->SetIdx(idx);
2506	>	atom->SetIdx(idx);
2507
2508		EndModify();
2509
2510		return(true);
2511	<	}
2511	>	}
2512
2513	<	bool OBMol::DeleteResidue(OBResidue *residue)
2514	<	{
2513	>	bool OBMol::DeleteResidue(OBResidue *residue)
2514	>	{
2515		unsigned short idx = residue->GetIdx();
2516		for ( unsigned short i = idx ; i < _residue.size() ; i++ )
2517	<	_residue[i]->SetIdx(i-1);
2517	>	_residue[i]->SetIdx(i-1);
2518
2519		_residue.erase(_residue.begin() + idx);
2520
2521		if (residue)
2522	<	{
2522	>	{
2523		delete residue;
2524		residue = NULL;
2525	<	}
2525	>	}
2526
2527		return(true);
2528	<	}
2528	>	}
2529
2530	<	bool OBMol::AddBond(int first,int second,int order,int stereo,int insertpos)
2531	<	{
2530	>	bool OBMol::AddBond(int first,int second,int order,int stereo,int insertpos)
2531	>	{
2532	>	if (first == second)
2533	>	return(false);
2534	>
2535		BeginModify();
2536	<
2536	>
2537		if ((unsigned)first <= NumAtoms() && (unsigned)second <= NumAtoms()
2538	<	&& !GetBond(first, second))
2539	<	//atoms exist and bond doesn't
2540	<	{
2541	<	OBBond *bond = CreateBond();
2542	<	if (!bond)
2543	<	{
2544	<	EndModify();
2538	>	&& !GetBond(first, second))
2539	>	//atoms exist and bond doesn't
2540	>	{
2541	>	OBBond *bond = CreateBond();
2542	>	if (!bond)
2543	>	{
2544	>	EndModify();
2545		return(false);
2546	<	}
2546	>	}
2547
2548		OBAtom *bgn,*end;
2549		bgn = GetAtom(first);
2550		end = GetAtom(second);
2551		if (!bgn || !end)
2552	<	{
2553	<	obErrorLog.ThrowError(__FUNCTION__, "Unable to add bond - invalid atom index", obDebug);
2552	>	{
2553	>	obErrorLog.ThrowError(__func__, "Unable to add bond - invalid atom index", obDebug);
2554		return(false);
2555	<	}
2555	>	}
2556		bond->Set(_nbonds,bgn,end,order,stereo);
2557		bond->SetParent(this);
2558
2559		//set aromatic flags if it has the appropriate order
2560		if (order == 5)
2561	<	{
2561	>	{
2562		bond->SetAromatic();
2563		bgn->SetAromatic();
2564		end->SetAromatic();
2565	<	}
2565	>	}
2566
2567		#define OBBondIncrement 100
2568		if (_vbond.empty() || _nbonds+1 >= (signed)_vbond.size())
2569	<	{
2569	>	{
2570		_vbond.resize(_nbonds+OBBondIncrement);
2571		vector<OBEdgeBase*>::iterator i;
2572		for (i = _vbond.begin(),i+=(_nbonds+1);i != _vbond.end();i++)
2573	<	*i = (OBEdgeBase*)NULL;
2574	<	}
2573	>	*i = (OBEdgeBase*)NULL;
2574	>	}
2575		#undef  OBBondIncrement
2576
2577		_vbond[_nbonds] = (OBEdgeBase*)bond;
2578		_nbonds++;
2579
2580		if (insertpos == -1)
2581	<	{
2581	>	{
2582		bgn->AddBond(bond);
2583		end->AddBond(bond);
2584	<	}
2584	>	}
2585		else
2586	<	{
2586	>	{
2587		if (insertpos >= static_cast<int>(bgn->GetValence())
2588	<	) bgn->AddBond(bond);
2588	>	) bgn->AddBond(bond);
2589		else //need to insert the bond for the connectivity order to be preserved
2590	<	{    //otherwise stereochemistry gets screwed up
2590	>	{    //otherwise stereochemistry gets screwed up
2591		vector<OBEdgeBase*>::iterator bi;
2592		bgn->BeginNbrAtom(bi);
2593		bi += insertpos;
2594		bgn->InsertBond(bi,bond);
2595	<	}
2595	>	}
2596		end->AddBond(bond);
2597	<	}
2598	<	}
2597	>	}
2598	>	}
2599		else //at least one atom doesn't exist yet - add to bond_q
2600	<	SetData(new OBVirtualBond(first,second,order,stereo));
2600	>	SetData(new OBVirtualBond(first,second,order,stereo));
2601
2602		EndModify();
2603		return(true);
2604	<	}
2604	>	}
2605
2606	<	bool OBMol::AddBond(OBBond &bond)
2607	<	{
2606	>	bool OBMol::AddBond(OBBond &bond)
2607	>	{
2608		return(AddBond(bond.GetBeginAtomIdx(),
2609		bond.GetEndAtomIdx(),
2610		bond.GetBO(),
2611		bond.GetFlags()));
2612	<	}
2612	>	}
2613
2614	<	bool OBMol::DeleteBond(OBBond *bond)
2615	<	{
2614	>	bool OBMol::DeleteBond(OBBond *bond)
2615	>	{
2616		BeginModify();
2617
2618		(bond->GetBeginAtom())->DeleteBond(bond);
#	Line 2480 | Line 2624 | bool OBMol::DeleteBond(OBBond *bond)
2624		vector<OBEdgeBase*>::iterator i;
2625		int j;
2626		for (bond = BeginBond(i),j=0;bond;bond = NextBond(i),j++)
2627	<	bond->SetIdx(j);
2627	>	bond->SetIdx(j);
2628
2629		_nbonds--;
2630		EndModify();
2631		return(true);
2632	<	}
2632	>	}
2633
2634	<	void OBMol::Align(OBAtom *a1,OBAtom *a2,vector3 &p1,vector3 &p2)
2635	<	{
2634	>	void OBMol::Align(OBAtom *a1,OBAtom *a2,vector3 &p1,vector3 &p2)
2635	>	{
2636		vector<int> children;
2637
2638	<	obErrorLog.ThrowError(__FUNCTION__,
2638	>	obErrorLog.ThrowError(__func__,
2639		"Ran OpenBabel::Align", obAuditMsg);
2640
2641		//find which atoms to rotate
#	Line 2514 | Line 2658 | void OBMol::Align(OBAtom *a1,OBAtom *a2,vector3 &p1,ve
2658		OBAtom *atom;
2659		vector<int>::iterator i;
2660		for (i = children.begin();i != children.end();i++)
2661	<	{
2661	>	{
2662		atom = GetAtom(*i);
2663		v = atom->GetVector();
2664		v -= a1->GetVector();
2665		v *= m;   //rotate the point
2666		v += p1;  //translate the vector
2667		atom->SetVector(v);
2668	<	}
2668	>	}
2669		//set a1 = p1
2670		a1->SetVector(p1);
2671	<	}
2671	>	}
2672
2673	<	void OBMol::ToInertialFrame()
2674	<	{
2673	>	void OBMol::ToInertialFrame()
2674	>	{
2675		double m[9];
2676		for (int i = 0;i < NumConformers();i++)
2677	<	ToInertialFrame(i,m);
2678	<	}
2677	>	ToInertialFrame(i,m);
2678	>	}
2679
2680	<	void OBMol::ToInertialFrame(int conf,double *rmat)
2681	<	{
2680	>	void OBMol::ToInertialFrame(int conf,double *rmat)
2681	>	{
2682		unsigned int i;
2683		double x,y,z;
2684		double mi;
2685		double mass = 0.0;
2686		double center[3],m[3][3];
2687
2688	<	obErrorLog.ThrowError(__FUNCTION__,
2688	>	obErrorLog.ThrowError(__func__,
2689		"Ran OpenBabel::ToInertialFrame", obAuditMsg);
2690
2691		for (i = 0;i < 3;i++)
2692	<	memset(&m[i],'\0',sizeof(double)*3);
2692	>	memset(&m[i],'\0',sizeof(double)*3);
2693		memset(center,'\0',sizeof(double)*3);
2694
2695		SetConformer(conf);
#	Line 2553 | Line 2697 | void OBMol::ToInertialFrame(int conf,double *rmat)
2697		vector<OBNodeBase*>::iterator j;
2698		//find center of mass
2699		for (atom = BeginAtom(j);atom;atom = NextAtom(j))
2700	<	{
2700	>	{
2701		mi = atom->GetAtomicMass();
2702		center[0] += mi*atom->x();
2703		center[1] += mi*atom->y();
2704		center[2] += mi*atom->z();
2705		mass += mi;
2706	<	}
2706	>	}
2707
2708		center[0] /= mass;
2709		center[1] /= mass;
#	Line 2567 | Line 2711 | void OBMol::ToInertialFrame(int conf,double *rmat)
2711
2712		//calculate inertial tensor
2713		for (atom = BeginAtom(j);atom;atom = NextAtom(j))
2714	<	{
2714	>	{
2715		x = atom->x()-center[0];
2716		y = atom->y()-center[1];
2717		z = atom->z()-center[2];
#	Line 2582 | Line 2726 | void OBMol::ToInertialFrame(int conf,double *rmat)
2726		m[2][0] -= mi*x*z;
2727		m[2][1] -= mi*y*z;
2728		m[2][2] += mi*(x*x+y*y);
2729	<	}
2729	>	}
2730
2731		/* find rotation matrix for moment of inertia */
2732		ob_make_rmat(m,rmat);
#	Line 2590 | Line 2734 | void OBMol::ToInertialFrame(int conf,double *rmat)
2734		/* rotate all coordinates */
2735		double *c = GetConformer(conf);
2736		for(i=0; i < NumAtoms();i++)
2737	<	{
2737	>	{
2738		x = c[i*3]-center[0];
2739		y = c[i*3+1]-center[1];
2740		z = c[i*3+2]-center[2];
2741		c[i*3]   = x*rmat[0] + y*rmat[1] + z*rmat[2];
2742		c[i*3+1] = x*rmat[3] + y*rmat[4] + z*rmat[5];
2743		c[i*3+2] = x*rmat[6] + y*rmat[7] + z*rmat[8];
2744	<	}
2745	<	}
2744	>	}
2745	>	}
2746
2747	<	/*NF
2748	<	istream& operator>> (istream &ifs, OBMol &mol)
2605	<	{
2606	<	bool retcode = OBFileFormat::ReadMolecule(ifs, mol);
2607	<
2608	<	if (!retcode)
2609	<	{
2610	<	if (mol.GetMod())
2611	<	mol.EndModify();
2612	<	mol.Clear();
2613	<	}
2614	<
2615	<	return(ifs);
2616	<	}
2617	<
2618	<	ostream& operator<< (ostream &ofs, OBMol &mol)
2619	<	{
2620	<	OBFileFormat::WriteMolecule(ofs, mol);
2621	<	return(ofs);
2622	<	}
2623	<	*/
2624	<
2625	<	OBMol::OBMol()
2626	<	{
2747	>	OBMol::OBMol()
2748	>	{
2749		_natoms = _nbonds = 0;
2750		_mod = 0;
2751		_energy = 0.0;
#	Line 2638 | Line 2760 | OBMol::OBMol()
2760		_vconf.clear();
2761		_autoPartialCharge = true;
2762		_autoFormalCharge = true;
2763	<	}
2763	>	}
2764
2765	<	OBMol::OBMol(const OBMol &mol) :
2766	<	OBGraphBase()
2767	<	{
2765	>	OBMol::OBMol(const OBMol &mol) :
2766	>	OBGraphBase()
2767	>	{
2768		_natoms = _nbonds = 0;
2769		_mod = 0;
2770		_totalCharge = 0;
#	Line 2657 | Line 2779 | OBMol::OBMol(const OBMol &mol) :
2779		_autoFormalCharge = true;
2780		//NF  _compressed = false;
2781		*this = mol;
2782	<	}
2782	>	}
2783
2784	<	OBMol::~OBMol()
2785	<	{
2784	>	OBMol::~OBMol()
2785	>	{
2786		OBAtom    *atom;
2787		OBBond    *bond;
2788		OBResidue *residue;
#	Line 2668 | Line 2790 | OBMol::~OBMol()
2790		vector<OBEdgeBase*>::iterator j;
2791		vector<OBResidue*>::iterator r;
2792		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
2793	<	DestroyAtom(atom);
2793	>	DestroyAtom(atom);
2794		for (bond = BeginBond(j);bond;bond = NextBond(j))
2795	<	DestroyBond(bond);
2795	>	DestroyBond(bond);
2796		for (residue = BeginResidue(r);residue;residue = NextResidue(r))
2797	<	delete residue;
2797	>	delete residue;
2798
2799		//clear out the multiconformer data
2800		vector<double*>::iterator k;
2801		for (k = _vconf.begin();k != _vconf.end();k++)
2802	<	delete [] *k;
2802	>	delete [] *k;
2803		_vconf.clear();
2804
2805		if (!_vdata.empty())
2806	<	{
2806	>	{
2807		vector<OBGenericData*>::iterator m;
2808		for (m = _vdata.begin();m != _vdata.end();m++)
2809	<	delete *m;
2809	>	delete *m;
2810		_vdata.clear();
2811	<	}
2812	<	}
2811	>	}
2812	>	}
2813
2814	<	bool OBMol::HasData(string &s)
2815	<	{
2814	>	bool OBMol::HasData(string &s)
2815	>	{
2816		if (_vdata.empty())
2817	<	return(false);
2817	>	return(false);
2818
2819		vector<OBGenericData*>::iterator i;
2820
2821		for (i = _vdata.begin();i != _vdata.end();i++)
2822	<	if ((*i)->GetAttribute() == s)
2823	<	return(true);
2822	>	if ((*i)->GetAttribute() == s)
2823	>	return(true);
2824
2825		return(false);
2826	<	}
2826	>	}
2827
2828	<	bool OBMol::HasData(const char *s)
2829	<	//returns true if the generic attribute/value pair exists
2830	<	{
2828	>	bool OBMol::HasData(const char *s)
2829	>	//returns true if the generic attribute/value pair exists
2830	>	{
2831		if (_vdata.empty())
2832	<	return(false);
2832	>	return(false);
2833
2834		vector<OBGenericData*>::iterator i;
2835
2836		for (i = _vdata.begin();i != _vdata.end();i++)
2837	<	if ((*i)->GetAttribute() == s)
2838	<	return(true);
2837	>	if ((*i)->GetAttribute() == s)
2838	>	return(true);
2839
2840		return(false);
2841	<	}
2841	>	}
2842
2843
2844	<	bool OBMol::HasData(unsigned int dt)
2845	<	//returns true if the generic attribute/value pair exists
2846	<	{
2844	>	bool OBMol::HasData(unsigned int dt)
2845	>	//returns true if the generic attribute/value pair exists
2846	>	{
2847		if (_vdata.empty())
2848	<	return(false);
2848	>	return(false);
2849
2850		vector<OBGenericData*>::iterator i;
2851
2852		for (i = _vdata.begin();i != _vdata.end();i++)
2853	<	if ((*i)->GetDataType() == dt)
2854	<	return(true);
2853	>	if ((*i)->GetDataType() == dt)
2854	>	return(true);
2855
2856		return(false);
2857	<	}
2857	>	}
2858
2859	<	//! Returns the value given an attribute name
2860	<	OBGenericData *OBMol::GetData(string &s)
2861	<	{
2859	>	//! Returns the value given an attribute name
2860	>	OBGenericData *OBMol::GetData(string &s)
2861	>	{
2862		vector<OBGenericData*>::iterator i;
2863
2864		for (i = _vdata.begin();i != _vdata.end();i++)
2865	<	if ((*i)->GetAttribute() == s)
2866	<	return(*i);
2865	>	if ((*i)->GetAttribute() == s)
2866	>	return(*i);
2867
2868		return(NULL);
2869	<	}
2869	>	}
2870
2871	<	//! Returns the value given an attribute name
2872	<	OBGenericData *OBMol::GetData(const char *s)
2873	<	{
2871	>	//! Returns the value given an attribute name
2872	>	OBGenericData *OBMol::GetData(const char *s)
2873	>	{
2874		vector<OBGenericData*>::iterator i;
2875
2876		for (i = _vdata.begin();i != _vdata.end();i++)
2877	<	if ((*i)->GetAttribute() == s)
2878	<	return(*i);
2877	>	if ((*i)->GetAttribute() == s)
2878	>	return(*i);
2879
2880		return(NULL);
2881	<	}
2881	>	}
2882
2883	<	OBGenericData *OBMol::GetData(unsigned int dt)
2884	<	{
2883	>	OBGenericData *OBMol::GetData(unsigned int dt)
2884	>	{
2885		vector<OBGenericData*>::iterator i;
2886		for (i = _vdata.begin();i != _vdata.end();i++)
2887	<	if ((*i)->GetDataType() == dt)
2888	<	return(*i);
2887	>	if ((*i)->GetDataType() == dt)
2888	>	return(*i);
2889		return(NULL);
2890	<	}
2890	>	}
2891
2892	<	void OBMol::DeleteData(unsigned int dt)
2893	<	{
2892	>	void OBMol::DeleteData(unsigned int dt)
2893	>	{
2894		vector<OBGenericData*> vdata;
2895		vector<OBGenericData*>::iterator i;
2896		for (i = _vdata.begin();i != _vdata.end();i++)
2897	<	if ((*i)->GetDataType() == dt)
2898	<	delete *i;
2899	<	else
2900	<	vdata.push_back(*i);
2897	>	if ((*i)->GetDataType() == dt)
2898	>	delete *i;
2899	>	else
2900	>	vdata.push_back(*i);
2901		_vdata = vdata;
2902	<	}
2902	>	}
2903
2904	<	void OBMol::DeleteData(vector<OBGenericData*> &vg)
2905	<	{
2904	>	void OBMol::DeleteData(vector<OBGenericData*> &vg)
2905	>	{
2906		vector<OBGenericData*> vdata;
2907		vector<OBGenericData*>::iterator i,j;
2908
2909		bool del;
2910		for (i = _vdata.begin();i != _vdata.end();i++)
2911	<	{
2911	>	{
2912		del = false;
2913		for (j = vg.begin();j != vg.end();j++)
2914	<	if (*i == *j)
2914	>	if (*i == *j)
2915		{
2916	<	del = true;
2917	<	break;
2916	>	del = true;
2917	>	break;
2918		}
2919		if (del)
2920	<	delete *i;
2920	>	delete *i;
2921		else
2922	<	vdata.push_back(*i);
2923	<	}
2922	>	vdata.push_back(*i);
2923	>	}
2924		_vdata = vdata;
2925	<	}
2925	>	}
2926
2927	<	void OBMol::DeleteData(OBGenericData *gd)
2928	<	{
2927	>	void OBMol::DeleteData(OBGenericData *gd)
2928	>	{
2929		vector<OBGenericData*>::iterator i;
2930		for (i = _vdata.begin();i != _vdata.end();i++)
2931	<	if (*i == gd)
2931	>	if (*i == gd)
2932		{
2933	<	delete *i;
2934	<	_vdata.erase(i);
2933	>	delete *i;
2934	>	_vdata.erase(i);
2935		}
2936
2937	<	}
2937	>	}
2938
2939	<	bool OBMol::HasNonZeroCoords()
2940	<	{
2939	>	bool OBMol::HasNonZeroCoords()
2940	>	{
2941		OBAtom *atom;
2942		vector<OBNodeBase*>::iterator i;
2943
2944		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
2945	<	if (atom->GetVector() != VZero)
2946	<	return(true);
2945	>	if (atom->GetVector() != VZero)
2946	>	return(true);
2947
2948		return(false);
2949	<	}
2949	>	}
2950
2951	<	bool OBMol::Has2D()
2952	<	{
2951	>	bool OBMol::Has2D()
2952	>	{
2953		bool hasX,hasY;
2954		OBAtom *atom;
2955		vector<OBNodeBase*>::iterator i;
2956
2957		hasX = hasY = false;
2958		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
2959	<	{
2959	>	{
2960		if (!hasX && !IsNearZero(atom->x()))
2961	<	hasX = true;
2961	>	hasX = true;
2962		if (!hasY && !IsNearZero(atom->y()))
2963	<	hasY = true;
2963	>	hasY = true;
2964
2965		if (hasX && hasY)
2966	<	return(true);
2967	<	}
2966	>	return(true);
2967	>	}
2968		return(false);
2969	<	}
2969	>	}
2970
2971	<	bool OBMol::Has3D()
2972	<	{
2971	>	bool OBMol::Has3D()
2972	>	{
2973		bool hasX,hasY,hasZ;
2974		OBAtom *atom;
2975		vector<OBNodeBase*>::iterator i;
2976
2977		hasX = hasY = hasZ = false;
2978		if (this->_c == NULL)
2979	<	return(false);
2979	>	return(false);
2980		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
2981	<	{
2981	>	{
2982		if (!hasX && !IsNearZero(atom->x()))
2983	<	hasX = true;
2983	>	hasX = true;
2984		if (!hasY && !IsNearZero(atom->y()))
2985	<	hasY = true;
2985	>	hasY = true;
2986		if (!hasZ && !IsNearZero(atom->z()))
2987	<	hasZ = true;
2987	>	hasZ = true;
2988
2989		if (hasX && hasY && hasZ)
2990	<	return(true);
2991	<	}
2990	>	return(true);
2991	>	}
2992		return(false);
2993	<	}
2993	>	}
2994
2995	<	bool OBMol::IsChiral()
2996	<	{
2995	>	bool OBMol::IsChiral()
2996	>	{
2997		OBAtom *atom;
2998		vector<OBNodeBase*>::iterator i;
2999
3000		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
3001	<	if ((atom->IsCarbon() || atom->IsNitrogen()) && atom->GetHvyValence() > 2 && atom->IsChiral())
3002	<	return(true);
3001	>	if ((atom->IsCarbon() || atom->IsNitrogen()) && atom->GetHvyValence() > 2 && atom->IsChiral())
3002	>	return(true);
3003
3004		return(false);
3005	<	}
3005	>	}
3006
3007	<	//! Renumber the atoms in this molecule according to the order in the supplied
3008	<	//! vector. This will return without action if the supplied vector is empty or
3009	<	//! does not have the same number of atoms as the molecule.
3010	<	void OBMol::RenumberAtoms(vector<OBNodeBase*> &v)
3011	<	{
3007	>	//! Renumber the atoms in this molecule according to the order in the supplied
3008	>	//! vector. This will return without action if the supplied vector is empty or
3009	>	//! does not have the same number of atoms as the molecule.
3010	>	void OBMol::RenumberAtoms(vector<OBNodeBase*> &v)
3011	>	{
3012		if (Empty())
3013	<	return;
3013	>	return;
3014
3015	<	obErrorLog.ThrowError(__FUNCTION__,
3015	>	obErrorLog.ThrowError(__func__,
3016		"Ran OpenBabel::RenumberAtoms", obAuditMsg);
3017
3018		OBAtom *atom;
#	Line 2916 | Line 3038 | void OBMol::RenumberAtoms(vector<OBNodeBase*> &v)
3038		double *ctmp = new double [NumAtoms()*3];
3039
3040		for (j = 0;j < NumConformers();j++)
3041	<	{
3041	>	{
3042		c = GetConformer(j);
3043		for (k=0,i = va.begin();i != va.end();i++,k++)
3044	<	memcpy((char*)&ctmp[k*3],(char*)&c[((OBAtom*)*i)->GetCIdx()],sizeof(double)*3);
3044	>	memcpy((char*)&ctmp[k*3],(char*)&c[((OBAtom*)*i)->GetCIdx()],sizeof(double)*3);
3045		memcpy((char*)c,(char*)ctmp,sizeof(double)*3*NumAtoms());
3046	<	}
3046	>	}
3047
3048		for (k=1,i = va.begin();i != va.end();i++,k++)
3049	<	(*i)->SetIdx(k);
3049	>	(*i)->SetIdx(k);
3050
3051		delete [] ctmp;
3052
3053		_vatom.clear();
3054		for (i = va.begin();i != va.end();i++)
3055	<	_vatom.push_back(*i);
3056	<	}
3055	>	_vatom.push_back(*i);
3056	>	}
3057
3058		#ifdef REMOVE_LATER
3059	<	bool CompareBonds(const OBEdgeBase *a,const OBEdgeBase *b)
3060	<	{
3059	>	bool CompareBonds(const OBEdgeBase *a,const OBEdgeBase *b)
3060	>	{
3061		if (a->GetBgn()->GetIdx() == b->GetBgn()->GetIdx())
3062	<	return(a->GetEnd()->GetIdx() < b->GetEnd()->GetIdx());
3062	>	return(a->GetEnd()->GetIdx() < b->GetEnd()->GetIdx());
3063
3064		//return((a->GetBgn())->GetIdx() < (b->GetBgn())->GetIdx());
3065	<	}
3065	>	}
3066		#endif
3067
3068	<	bool WriteTitles(ostream &ofs, OBMol &mol)
3069	<	{
3068	>	bool WriteTitles(ostream &ofs, OBMol &mol)
3069	>	{
3070		ofs << mol.GetTitle() << endl;
3071		return true;
3072	<	}
3072	>	}
3073
3074	<	/*! This method adds single bonds between all atoms
3075	<	closer than their combined atomic covalent radii,
3076	<	then "cleans up" making sure bonded atoms are not
3077	<	closer than 0.4A and the atom does not exceed its valence.
3078	<	It implements blue-obelisk:rebondFrom3DCoordinates.
3074	>	/*! This method adds single bonds between all atoms
3075	>	closer than their combined atomic covalent radii,
3076	>	then "cleans up" making sure bonded atoms are not
3077	>	closer than 0.4A and the atom does not exceed its valence.
3078	>	It implements blue-obelisk:rebondFrom3DCoordinates.
3079
3080	<	*/
3081	<	void OBMol::ConnectTheDots(void)
3082	<	{
3080	>	*/
3081	>	void OBMol::ConnectTheDots(void)
3082	>	{
3083		if (Empty())
3084	<	return;
3084	>	return;
3085		if (_dimension != 3) return; // not useful on non-3D structures
3086
3087	<	obErrorLog.ThrowError(__FUNCTION__,
3087	>	obErrorLog.ThrowError(__func__,
3088		"Ran OpenBabel::ConnectTheDots", obAuditMsg);
3089
3090		int j,k,max;
#	Line 2977 | Line 3099 | void OBMol::ConnectTheDots(void)
3099		rad.resize(_natoms);
3100
3101		for (j = 0, atom = BeginAtom(i) ; atom ; atom = NextAtom(i), j++)
3102	<	{
3102	>	{
3103		(atom->GetVector()).Get(&c[j*3]);
3104		pair<OBAtom*,double> entry(atom, atom->GetVector().z());
3105		zsortedAtoms.push_back(entry);
3106	<	}
3106	>	}
3107		sort(zsortedAtoms.begin(), zsortedAtoms.end(), SortAtomZ);
3108
3109		max = zsortedAtoms.size();
3110
3111		for ( j = 0 ; j < max ; j++ )
3112	<	{
3112	>	{
3113		atom   = zsortedAtoms[j].first;
3114		rad[j] = etab.GetCovalentRad(atom->GetAtomicNum());
3115		zsorted.push_back(atom->GetIdx()-1);
3116	<	}
3116	>	}
3117
3118		int idx1, idx2;
3119		double d2,cutoff,zd;
3120		for (j = 0 ; j < max ; j++)
3121	<	{
3121	>	{
3122		idx1 = zsorted[j];
3123		for (k = j + 1 ; k < max ; k++ )
3124	<	{
3124	>	{
3125		idx2 = zsorted[k];
3126
3127		// bonded if closer than elemental Rcov + tolerance
#	Line 3007 | Line 3129 | void OBMol::ConnectTheDots(void)
3129
3130		zd  = SQUARE(c[idx1*3+2] - c[idx2*3+2]);
3131		if (zd > 25.0 )
3132	<	break; // bigger than max cutoff
3132	>	break; // bigger than max cutoff
3133
3134		d2  = SQUARE(c[idx1*3]   - c[idx2*3]);
3135		d2 += SQUARE(c[idx1*3+1] - c[idx2*3+1]);
3136		d2 += zd;
3137
3138		if (d2 > cutoff)
3139	<	continue;
3139	>	continue;
3140		if (d2 < 0.40)
3141	<	continue;
3141	>	continue;
3142
3143		atom = GetAtom(idx1+1);
3144		nbr  = GetAtom(idx2+1);
3145
3146		if (atom->IsConnected(nbr))
3147	<	continue;
3147	>	continue;
3148		if (atom->IsHydrogen() && nbr->IsHydrogen())
3149	<	continue;
3149	>	continue;
3150
3151		AddBond(idx1+1,idx2+1,1);
3152	<	}
3153	<	}
3152	>	}
3153	>	}
3154
3155		// If between BeginModify and EndModify, coord pointers are NULL
3156		// setup molecule to handle current coordinates
3157
3158		if (_c == NULL)
3159	<	{
3159	>	{
3160		_c = c;
3161		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
3162	<	atom->SetCoordPtr(&_c);
3162	>	atom->SetCoordPtr(&_c);
3163		_vconf.push_back(c);
3164		unset = true;
3165	<	}
3165	>	}
3166
3167		// Cleanup -- delete long bonds that exceed max valence
3168		OBBond *maxbond, *bond;
3169		double maxlength;
3170		vector<OBEdgeBase*>::iterator l;
3171		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
3172	<	{
3172	>	{
3173		while (atom->BOSum() > static_cast<unsigned int>(etab.GetMaxBonds(atom->GetAtomicNum()))
3174	<	|| atom->SmallestBondAngle() < 45.0)
3175	<	{
3174	>	|| atom->SmallestBondAngle() < 45.0)
3175	>	{
3176		maxbond = atom->BeginBond(l);
3177		maxlength = maxbond->GetLength();
3178		for (bond = atom->BeginBond(l);bond;bond = atom->NextBond(l))
3179	<	{
3179	>	{
3180		if (bond->GetLength() > maxlength)
3181	<	{
3181	>	{
3182		maxbond = bond;
3183		maxlength = bond->GetLength();
3184	<	}
3185	<	}
3184	>	}
3185	>	}
3186		DeleteBond(maxbond);
3187	<	}
3188	<	}
3187	>	}
3188	>	}
3189
3190		if (unset)
3191	<	{
3191	>	{
3192		_c = NULL;
3193		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
3194	<	atom->ClearCoordPtr();
3194	>	atom->ClearCoordPtr();
3195		_vconf.resize(_vconf.size()-1);
3196	<	}
3196	>	}
3197
3198		delete [] c;
3199	<	}
3199	>	}
3200
3201	<	/*! This method uses bond angles and geometries from current
3202	<	connectivity to guess atom types and then filling empty valences
3203	<	with multiple bonds. It currently has a pass to detect some
3204	<	frequent functional groups. It still needs a pass to detect aromatic
3205	<	rings to "clean up." */
3206	<	void OBMol::PerceiveBondOrders()
3207	<	{
3201	>	/*! This method uses bond angles and geometries from current
3202	>	connectivity to guess atom types and then filling empty valences
3203	>	with multiple bonds. It currently has a pass to detect some
3204	>	frequent functional groups. It still needs a pass to detect aromatic
3205	>	rings to "clean up." */
3206	>	void OBMol::PerceiveBondOrders()
3207	>	{
3208		if (Empty())
3209	<	return;
3209	>	return;
3210		if (_dimension != 3) return; // not useful on non-3D structures
3211
3212	<	obErrorLog.ThrowError(__FUNCTION__,
3212	>	obErrorLog.ThrowError(__func__,
3213		"Ran OpenBabel::PerceiveBondOrders", obAuditMsg);
3214
3215		OBAtom *atom, *b, *c;
#	Line 3100 | Line 3222 | void OBMol::PerceiveBondOrders()
3222
3223		// Pass 1: Assign estimated hybridization based on avg. angles
3224		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
3225	<	{
3225	>	{
3226		angle = atom->AverageBondAngle();
3227
3228		if (angle > 155.0)
3229	<	atom->SetHyb(1);
3229	>	atom->SetHyb(1);
3230		else if ( angle <= 155.0 && angle > 115)
3231	<	atom->SetHyb(2);
3232	<	} // pass 1
3231	>	atom->SetHyb(2);
3232	>	} // pass 1
3233
3234		// Make sure upcoming calls to GetHyb() don't kill these temporary values
3235		SetHybridizationPerceived();
#	Line 3122 | Line 3244 | void OBMol::PerceiveBondOrders()
3244		double torsions = 0.0;
3245
3246		if (!HasSSSRPerceived())
3247	<	FindSSSR();
3247	>	FindSSSR();
3248		rlist = GetSSSR();
3249		for (ringit = rlist.begin(); ringit != rlist.end(); ringit++)
3250	<	{
3250	>	{
3251		if ((*ringit)->PathSize() == 5)
3252	<	{
3252	>	{
3253		path = (*ringit)->_path;
3254		torsions =
3255	<	( fabs(GetTorsion(path[0], path[1], path[2], path[3])) +
3256	<	fabs(GetTorsion(path[1], path[2], path[3], path[4])) +
3257	<	fabs(GetTorsion(path[2], path[3], path[4], path[0])) +
3258	<	fabs(GetTorsion(path[3], path[4], path[0], path[1])) +
3259	<	fabs(GetTorsion(path[4], path[0], path[1], path[2])) ) / 5.0;
3255	>	( fabs(GetTorsion(path[0], path[1], path[2], path[3])) +
3256	>	fabs(GetTorsion(path[1], path[2], path[3], path[4])) +
3257	>	fabs(GetTorsion(path[2], path[3], path[4], path[0])) +
3258	>	fabs(GetTorsion(path[3], path[4], path[0], path[1])) +
3259	>	fabs(GetTorsion(path[4], path[0], path[1], path[2])) ) / 5.0;
3260		if (torsions <= 7.5)
3261	<	{
3261	>	{
3262		for (unsigned int ringAtom = 0; ringAtom != path.size(); ringAtom++)
3263	<	{
3263	>	{
3264		b = GetAtom(path[ringAtom]);
3265		// if an aromatic ring atom has valence 3, it is already set
3266		// to sp2 because the average angles should be 120 anyway
3267		// so only look for valence 2
3268		if (b->GetValence() == 2)
3269	<	b->SetHyb(2);
3270	<	}
3271	<	}
3272	<	}
3269	>	b->SetHyb(2);
3270	>	}
3271	>	}
3272	>	}
3273		else if ((*ringit)->PathSize() == 6)
3274	<	{
3274	>	{
3275		path = (*ringit)->_path;
3276		torsions =
3277	<	( fabs(GetTorsion(path[0], path[1], path[2], path[3])) +
3278	<	fabs(GetTorsion(path[1], path[2], path[3], path[4])) +
3279	<	fabs(GetTorsion(path[2], path[3], path[4], path[5])) +
3280	<	fabs(GetTorsion(path[3], path[4], path[5], path[0])) +
3281	<	fabs(GetTorsion(path[4], path[5], path[0], path[1])) +
3282	<	fabs(GetTorsion(path[5], path[0], path[1], path[2])) ) / 6.0;
3277	>	( fabs(GetTorsion(path[0], path[1], path[2], path[3])) +
3278	>	fabs(GetTorsion(path[1], path[2], path[3], path[4])) +
3279	>	fabs(GetTorsion(path[2], path[3], path[4], path[5])) +
3280	>	fabs(GetTorsion(path[3], path[4], path[5], path[0])) +
3281	>	fabs(GetTorsion(path[4], path[5], path[0], path[1])) +
3282	>	fabs(GetTorsion(path[5], path[0], path[1], path[2])) ) / 6.0;
3283		if (torsions <= 12.0)
3284	<	{
3284	>	{
3285		for (unsigned int ringAtom = 0; ringAtom != path.size(); ringAtom++)
3286	<	{
3286	>	{
3287		b = GetAtom(path[ringAtom]);
3288		if (b->GetValence() == 2 || b->GetValence() == 3)
3289	<	b->SetHyb(2);
3290	<	}
3291	<	}
3292	<	}
3293	<	}
3289	>	b->SetHyb(2);
3290	>	}
3291	>	}
3292	>	}
3293	>	}
3294
3295		// Pass 3: "Antialiasing" If an atom marked as sp hybrid isn't
3296		//          bonded to another or an sp2 hybrid isn't bonded
#	Line 3176 | Line 3298 | void OBMol::PerceiveBondOrders()
3298		//          mark them to a lower hybridization for now
3299		bool openNbr;
3300		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
3301	<	{
3301	>	{
3302		if (atom->GetHyb() == 2 || atom->GetHyb() == 1)
3303	<	{
3303	>	{
3304		openNbr = false;
3305		for (b = atom->BeginNbrAtom(j); b; b = atom->NextNbrAtom(j))
3306	<	{
3306	>	{
3307		if (b->GetHyb() < 3 || b->GetValence() == 1)
3308	<	{
3308	>	{
3309		openNbr = true;
3310		break;
3311	<	}
3312	<	}
3311	>	}
3312	>	}
3313		if (!openNbr && atom->GetHyb() == 2)
3314	<	atom->SetHyb(3);
3314	>	atom->SetHyb(3);
3315		else if (!openNbr && atom->GetHyb() == 1)
3316	<	atom->SetHyb(2);
3317	<	}
3318	<	} // pass 3
3316	>	atom->SetHyb(2);
3317	>	}
3318	>	} // pass 3
3319
3320		// Pass 4: Check for known functional group patterns and assign bonds
3321		//         to the canonical form
#	Line 3214 | Line 3336 | void OBMol::PerceiveBondOrders()
3336		bool typed; // has this ring been typed?
3337		unsigned int loop, loopSize;
3338		for (ringit = rlist.begin(); ringit != rlist.end(); ringit++)
3339	<	{
3339	>	{
3340		typed = false;
3341		loopSize = (*ringit)->PathSize();
3342		if (loopSize == 5 || loopSize == 6)
3343	<	{
3343	>	{
3344		path = (*ringit)->_path;
3345		for(loop = 0; loop < loopSize; loop++)
3346	<	{
3346	>	{
3347		atom = GetAtom(path[loop]);
3348		if(atom->HasBondOfOrder(2) || atom->HasBondOfOrder(3)
3349	<	|| atom->GetHyb() != 2)
3350	<	{
3349	>	|| atom->GetHyb() != 2)
3350	>	{
3351		typed = true;
3352		break;
3353	<	}
3354	<	}
3353	>	}
3354	>	}
3355
3356		if (!typed)
3357	<	for(loop = 0; loop < loopSize; loop++)
3357	>	for(loop = 0; loop < loopSize; loop++)
3358		{
3359	<	//          cout << " set aromatic " << path[loop] << endl;
3360	<	(GetBond(path[loop], path[(loop+1) % loopSize]))->SetBO(5);
3361	<	(GetBond(path[loop], path[(loop+1) % loopSize]))->UnsetKekule();
3359	>	//            cout << " set aromatic " << path[loop] << endl;
3360	>	(GetBond(path[loop], path[(loop+1) % loopSize]))->SetBO(5);
3361	>	(GetBond(path[loop], path[(loop+1) % loopSize]))->UnsetKekule();
3362		}
3363	<	}
3364	<	}
3363	>	}
3364	>	}
3365		_flags &= (~(OB_KEKULE_MOL));
3366		Kekulize();
3367
#	Line 3251 | Line 3373 | void OBMol::PerceiveBondOrders()
3373		double maxElNeg, shortestBond, currentElNeg;
3374
3375		for (atom = BeginAtom(i) ; atom ; atom = NextAtom(i))
3376	<	{
3377	<	// if atoms have the same electronegativity, make sure those with shorter bonds
3378	<	// are handled first (helps with assignment of conjugated single/double bonds)
3379	<	shortestBond = 1.0e5f;
3380	<	for (b = atom->BeginNbrAtom(j); b; b = atom->NextNbrAtom(j))
3381	<	{
3382	<	if (b->GetAtomicNum()!=1) shortestBond =
3383	<	min(shortestBond,(atom->GetBond(b))->GetLength());
3384	<	}
3385	<	pair<OBAtom*,double> entry(atom,
3386	<	etab.GetElectroNeg(atom->GetAtomicNum())*1e6+shortestBond);
3376	>	{
3377	>	// if atoms have the same electronegativity, make sure those with shorter bonds
3378	>	// are handled first (helps with assignment of conjugated single/double bonds)
3379	>	shortestBond = 1.0e5f;
3380	>	for (b = atom->BeginNbrAtom(j); b; b = atom->NextNbrAtom(j))
3381	>	{
3382	>	if (b->GetAtomicNum()!=1) shortestBond =
3383	>	min(shortestBond,(atom->GetBond(b))->GetLength());
3384	>	}
3385	>	pair<OBAtom*,double> entry(atom,
3386	>	etab.GetElectroNeg(atom->GetAtomicNum())*1e6+shortestBond);
3387
3388	<	sortedAtoms.push_back(entry);
3389	<	}
3388	>	sortedAtoms.push_back(entry);
3389	>	}
3390		sort(sortedAtoms.begin(), sortedAtoms.end(), SortAtomZ);
3391
3392		max = sortedAtoms.size();
3393
3394		for (iter = 0 ; iter < max ; iter++ )
3395	<	{
3395	>	{
3396		atom = sortedAtoms[iter].first;
3397		if ( (atom->GetHyb() == 1 || atom->GetValence() == 1)
3398	<	&& atom->BOSum() + 2  <= static_cast<unsigned int>(etab.GetMaxBonds(atom->GetAtomicNum()))
3399	<	)
3400	<	{
3398	>	&& atom->BOSum() + 2  <= static_cast<unsigned int>(etab.GetMaxBonds(atom->GetAtomicNum()))
3399	>	)
3400	>	{
3401		// loop through the neighbors looking for a hybrid or terminal atom
3402		// (and pick the one with highest electronegativity first)
3403		// *or* pick a neighbor that's a terminal atom
3404
3405		if (atom->HasNonSingleBond() ||
3406	<	(atom->GetAtomicNum() == 7 && atom->BOSum() + 2 > 3))
3407	<	continue;
3406	>	(atom->GetAtomicNum() == 7 && atom->BOSum() + 2 > 3))
3407	>	continue;
3408
3409		maxElNeg = 0.0;
3410		shortestBond = 5000.0;
3411		c = NULL;
3412		for (b = atom->BeginNbrAtom(j); b; b = atom->NextNbrAtom(j))
3413	<	{
3413	>	{
3414		currentElNeg = etab.GetElectroNeg(b->GetAtomicNum());
3415		if ( (b->GetHyb() == 1 || b->GetValence() == 1)
3416	<	&& b->BOSum() + 2 <= static_cast<unsigned int>(etab.GetMaxBonds(b->GetAtomicNum()))
3417	<	&& (currentElNeg > maxElNeg ||
3418	<	(IsNear(currentElNeg,maxElNeg)
3419	<	&& (atom->GetBond(b))->GetLength() < shortestBond)) )
3420	<	{
3416	>	&& b->BOSum() + 2 <= static_cast<unsigned int>(etab.GetMaxBonds(b->GetAtomicNum()))
3417	>	&& (currentElNeg > maxElNeg ||
3418	>	(IsNear(currentElNeg,maxElNeg)
3419	>	&& (atom->GetBond(b))->GetLength() < shortestBond)) )
3420	>	{
3421		if (b->HasNonSingleBond() ||
3422	<	(b->GetAtomicNum() == 7 && b->BOSum() + 2 > 3))
3423	<	continue;
3422	>	(b->GetAtomicNum() == 7 && b->BOSum() + 2 > 3))
3423	>	continue;
3424
3425		shortestBond = (atom->GetBond(b))->GetLength();
3426		maxElNeg = etab.GetElectroNeg(b->GetAtomicNum());
3427		c = b; // save this atom for later use
3428	<	}
3429	<	}
3428	>	}
3429	>	}
3430		if (c)
3431	<	(atom->GetBond(c))->SetBO(3);
3432	<	}
3431	>	(atom->GetBond(c))->SetBO(3);
3432	>	}
3433		else if ( (atom->GetHyb() == 2 || atom->GetValence() == 1)
3434		&& atom->BOSum() + 1 <= static_cast<unsigned int>(etab.GetMaxBonds(atom->GetAtomicNum())) )
3435	<	{
3435	>	{
3436		// as above
3437		if (atom->HasNonSingleBond() ||
3438	<	(atom->GetAtomicNum() == 7 && atom->BOSum() + 1 > 3))
3439	<	continue;
3438	>	(atom->GetAtomicNum() == 7 && atom->BOSum() + 1 > 3))
3439	>	continue;
3440
3441		maxElNeg = 0.0;
3442		shortestBond = 5000.0;
3443		c = NULL;
3444		for (b = atom->BeginNbrAtom(j); b; b = atom->NextNbrAtom(j))
3445	<	{
3445	>	{
3446		currentElNeg = etab.GetElectroNeg(b->GetAtomicNum());
3447		if ( (b->GetHyb() == 2 || b->GetValence() == 1)
3448		&& b->BOSum() + 1 <= static_cast<unsigned int>(etab.GetMaxBonds(b->GetAtomicNum()))
#	Line 3331 | Line 3453 | void OBMol::PerceiveBondOrders()
3453		&& (((atom->GetBond(b))->GetLength() < shortestBond)
3454		&& (!atom->IsInRing() || !c || !c->IsInRing() || b->IsInRing()))
3455		|| (atom->IsInRing() && c && !c->IsInRing() && b->IsInRing()))))
3456	<	{
3456	>	{
3457		if (b->HasNonSingleBond() ||
3458	<	(b->GetAtomicNum() == 7 && b->BOSum() + 1 > 3))
3459	<	continue;
3458	>	(b->GetAtomicNum() == 7 && b->BOSum() + 1 > 3))
3459	>	continue;
3460
3461		shortestBond = (atom->GetBond(b))->GetLength();
3462		maxElNeg = etab.GetElectroNeg(b->GetAtomicNum());
3463		c = b; // save this atom for later use
3464	<	}
3465	<	}
3464	>	}
3465	>	}
3466		if (c)
3467	<	(atom->GetBond(c))->SetBO(2);
3468	<	}
3469	<	} // pass 6
3467	>	(atom->GetBond(c))->SetBO(2);
3468	>	}
3469	>	} // pass 6
3470
3471		// Now let the atom typer go to work again
3472		_flags &= (~(OB_HYBRID_MOL));
#	Line 3353 | Line 3475 | void OBMol::PerceiveBondOrders()
3475		_flags &= (~(OB_ATOMTYPES_MOL));
3476		_flags &= (~(OB_IMPVAL_MOL));
3477		//  EndModify(true); // "nuke" perceived data
3478	<	}
3478	>	}
3479
3480	<	void OBMol::Center()
3481	<	{
3480	>	void OBMol::Center()
3481	>	{
3482		int j,size;
3483		double *c,x,y,z,fsize;
3484
3485		size = NumAtoms();
3486		fsize = -1.0/(double)NumAtoms();
3487
3488	<	obErrorLog.ThrowError(__FUNCTION__,
3488	>	obErrorLog.ThrowError(__func__,
3489		"Ran OpenBabel::Center", obAuditMsg);
3490
3491		vector<double*>::iterator i;
3492		for (i = _vconf.begin();i != _vconf.end();i++)
3493	<	{
3493	>	{
3494		c = *i;
3495		x = y = z = 0.0;
3496		for (j = 0;j < size;j++)
3497	<	{
3497	>	{
3498		x += c[j*3];
3499		y += c[j*3+1];
3500		z += c[j*3+2];
3501	<	}
3501	>	}
3502		x *= fsize;
3503		y *= fsize;
3504		z *= fsize;
3505
3506		for (j = 0;j < size;j++)
3507	<	{
3507	>	{
3508		c[j*3]+=x;
3509		c[j*3+1]+=y;
3510		c[j*3+2]+=z;
3511	<	}
3512	<	}
3511	>	}
3512	>	}
3513
3514	<	}
3514	>	}
3515
3516	<	vector3 OBMol::Center(int nconf)
3517	<	{
3518	<	obErrorLog.ThrowError(__FUNCTION__,
3516	>	vector3 OBMol::Center(int nconf)
3517	>	{
3518	>	obErrorLog.ThrowError(__func__,
3519		"Ran OpenBabel::Center", obAuditMsg);
3520
3521		SetConformer(nconf);
#	Line 3403 | Line 3525 | vector3 OBMol::Center(int nconf)
3525
3526		double x=0.0,y=0.0,z=0.0;
3527		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
3528	<	{
3528	>	{
3529		x += atom->x();
3530		y += atom->y();
3531		z += atom->z();
3532	<	}
3532	>	}
3533
3534		x /= (double)NumAtoms();
3535		y /= (double)NumAtoms();
#	Line 3417 | Line 3539 | vector3 OBMol::Center(int nconf)
3539		vector3 v(x,y,z);
3540
3541		for (atom = BeginAtom(i);atom;atom = NextAtom(i))
3542	<	{
3542	>	{
3543		vtmp = atom->GetVector() - v;
3544		atom->SetVector(vtmp);
3545	<	}
3545	>	}
3546
3547		return(v);
3548	<	}
3548	>	}
3549
3550
3551	<	/*! this method adds the vector v to all atom positions in all conformers */
3552	<	void OBMol::Translate(const vector3 &v)
3553	<	{
3551	>	/*! this method adds the vector v to all atom positions in all conformers */
3552	>	void OBMol::Translate(const vector3 &v)
3553	>	{
3554		for (int i = 0;i < NumConformers();i++)
3555	<	Translate(v,i);
3556	<	}
3555	>	Translate(v,i);
3556	>	}
3557
3558	<	/*! this method adds the vector v to all atom positions in the
3559	<	conformer nconf. If nconf == OB_CURRENT_CONFORMER, then the atom
3560	<	positions in the current conformer are translated. */
3561	<	void OBMol::Translate(const vector3 &v,int nconf)
3562	<	{
3563	<	obErrorLog.ThrowError(__FUNCTION__,
3558	>	/*! this method adds the vector v to all atom positions in the
3559	>	conformer nconf. If nconf == OB_CURRENT_CONFORMER, then the atom
3560	>	positions in the current conformer are translated. */
3561	>	void OBMol::Translate(const vector3 &v,int nconf)
3562	>	{
3563	>	obErrorLog.ThrowError(__func__,
3564		"Ran OpenBabel::Translate", obAuditMsg);
3565
3566		int i,size;
#	Line 3450 | Line 3572 | void OBMol::Translate(const vector3 &v,int nconf)
3572		z = v.z();
3573		size = NumAtoms();
3574		for (i = 0;i < size;i++)
3575	<	{
3575	>	{
3576		c[i*3  ] += x;
3577		c[i*3+1] += y;
3578		c[i*3+2] += z;
3579	<	}
3580	<	}
3579	>	}
3580	>	}
3581
3582	<	void OBMol::Rotate(const double u[3][3])
3583	<	{
3582	>	void OBMol::Rotate(const double u[3][3])
3583	>	{
3584		int i,j,k;
3585		double m[9];
3586		for (k=0,i = 0;i < 3;i++)
3587	<	for (j = 0;j < 3;j++)
3588	<	m[k++] = u[i][j];
3587	>	for (j = 0;j < 3;j++)
3588	>	m[k++] = u[i][j];
3589
3590		for (i = 0;i < NumConformers();i++)
3591	<	Rotate(m,i);
3592	<	}
3591	>	Rotate(m,i);
3592	>	}
3593
3594	<	void OBMol::Rotate(const double m[9])
3595	<	{
3594	>	void OBMol::Rotate(const double m[9])
3595	>	{
3596		for (int i = 0;i < NumConformers();i++)
3597	<	Rotate(m,i);
3598	<	}
3597	>	Rotate(m,i);
3598	>	}
3599
3600	<	void OBMol::Rotate(const double m[9],int nconf)
3601	<	{
3600	>	void OBMol::Rotate(const double m[9],int nconf)
3601	>	{
3602		int i,size;
3603		double x,y,z;
3604		double *c = (nconf == OB_CURRENT_CONFORMER)? _c : GetConformer(nconf);
3605
3606	<	obErrorLog.ThrowError(__FUNCTION__,
3606	>	obErrorLog.ThrowError(__func__,
3607		"Ran OpenBabel::Rotate", obAuditMsg);
3608
3609		size = NumAtoms();
3610		for (i = 0;i < size;i++)
3611	<	{
3611	>	{
3612		x = c[i*3  ];
3613		y = c[i*3+1];
3614		z = c[i*3+2];
3615		c[i*3  ] = m[0]*x + m[1]*y + m[2]*z;
3616		c[i*3+1] = m[3]*x + m[4]*y + m[5]*z;
3617		c[i*3+2] = m[6]*x + m[7]*y + m[8]*z;
3618	<	}
3619	<	}
3618	>	}
3619	>	}
3620
3621
3622	<	void OBMol::SetConformers(vector<double*> &v)
3623	<	{
3622	>	void OBMol::SetConformers(vector<double*> &v)
3623	>	{
3624		vector<double*>::iterator i;
3625		for (i = _vconf.begin();i != _vconf.end();i++)
3626	<	delete [] *i;
3626	>	delete [] *i;
3627
3628		_vconf = v;
3629		_c = (_vconf.empty()) ? NULL : _vconf[0];
3630
3631	<	}
3631	>	}
3632
3633	<	void OBMol::CopyConformer(double *c,int idx)
3634	<	{
3635	<	//    obAssert(!_vconf.empty() && (unsigned)idx < _vconf.size());
3633	>	void OBMol::CopyConformer(double *c,int idx)
3634	>	{
3635	>	//    obAssert(!_vconf.empty() && (unsigned)idx < _vconf.size());
3636		memcpy((char*)_vconf[idx],(char*)c,sizeof(double)*3*NumAtoms());
3637	<	}
3637	>	}
3638
3639	<	// void OBMol::CopyConformer(double *c,int idx)
3640	<	// {
3641	<	//   obAssert(!_vconf.empty() && (unsigned)idx < _vconf.size());
3639	>	// void OBMol::CopyConformer(double *c,int idx)
3640	>	// {
3641	>	//   obAssert(!_vconf.empty() && (unsigned)idx < _vconf.size());
3642
3643	<	//   unsigned int i;
3644	<	//   for (i = 0;i < NumAtoms();i++)
3645	<	//     {
3646	<	//       _vconf[idx][i*3  ] = (double)c[i*3  ];
3647	<	//       _vconf[idx][i*3+1] = (double)c[i*3+1];
3648	<	//       _vconf[idx][i*3+2] = (double)c[i*3+2];
3649	<	//     }
3650	<	// }
3643	>	//   unsigned int i;
3644	>	//   for (i = 0;i < NumAtoms();i++)
3645	>	//     {
3646	>	//       _vconf[idx][i*3  ] = (double)c[i*3  ];
3647	>	//       _vconf[idx][i*3+1] = (double)c[i*3+1];
3648	>	//       _vconf[idx][i*3+2] = (double)c[i*3+2];
3649	>	//     }
3650	>	// }
3651
3652	<	void OBMol::DeleteConformer(int idx)
3653	<	{
3652	>	void OBMol::DeleteConformer(int idx)
3653	>	{
3654		if (idx < 0 || idx >= (signed)_vconf.size())
3655	<	return;
3655	>	return;
3656
3657		delete [] _vconf[idx];
3658		_vconf.erase((_vconf.begin()+idx));
3659	<	}
3659	>	}
3660
3661	<	///Converts for instance [N+]([O-])=O to N(=O)=O
3662	<	bool OBMol::ConvertDativeBonds()
3663	<	{
3664	<	obErrorLog.ThrowError(__FUNCTION__,
3665	<	"Ran OpenBabel::ConvertDativeBonds", obAuditMsg);
3661	>	///Converts for instance [N+]([O-])=O to N(=O)=O
3662	>	bool OBMol::ConvertDativeBonds()
3663	>	{
3664	>	obErrorLog.ThrowError(__func__,
3665	>	"Ran OpenBabel::ConvertDativeBonds", obAuditMsg);
3666
3667	<	//Look for + and - charges on adjacent atoms
3668	<	OBAtom* patom;
3669	<	vector<OBNodeBase*>::iterator i;
3670	<	for (patom = BeginAtom(i);patom;patom = NextAtom(i))
3671	<	{
3672	<	vector<OBEdgeBase*>::iterator itr;
3673	<	OBBond *pbond;
3674	<	for (pbond = patom->BeginBond(itr);patom->GetFormalCharge() && pbond;pbond = patom->NextBond(itr))
3675	<	{
3676	<	OBAtom* pNbratom = pbond->GetNbrAtom(patom);
3677	<	int chg1 = patom->GetFormalCharge();
3678	<	int chg2 = pNbratom->GetFormalCharge();
3679	<	if((chg1>0 && chg2<0)|| (chg1<0 && chg2>0))
3680	<	{
3681	<	//dative bond. Reduce charges and increase bond order
3682	<	if(chg1>0)
3683	<	--chg1;
3684	<	else
3685	<	++chg1;
3686	<	patom->SetFormalCharge(chg1);
3687	<	if(chg2>0)
3688	<	--chg2;
3689	<	else
3690	<	++chg2;
3691	<	pNbratom->SetFormalCharge(chg2);
3692	<	pbond->SetBO(pbond->GetBO()+1);
3693	<	}
3694	<	}
3695	<	}
3696	<	return true;
3697	<	}
3667	>	//Look for + and - charges on adjacent atoms
3668	>	OBAtom* patom;
3669	>	vector<OBNodeBase*>::iterator i;
3670	>	for (patom = BeginAtom(i);patom;patom = NextAtom(i))
3671	>	{
3672	>	vector<OBEdgeBase*>::iterator itr;
3673	>	OBBond *pbond;
3674	>	for (pbond = patom->BeginBond(itr);patom->GetFormalCharge() && pbond;pbond = patom->NextBond(itr))
3675	>	{
3676	>	OBAtom* pNbratom = pbond->GetNbrAtom(patom);
3677	>	int chg1 = patom->GetFormalCharge();
3678	>	int chg2 = pNbratom->GetFormalCharge();
3679	>	if((chg1>0 && chg2<0)|| (chg1<0 && chg2>0))
3680	>	{
3681	>	//dative bond. Reduce charges and increase bond order
3682	>	if(chg1>0)
3683	>	--chg1;
3684	>	else
3685	>	++chg1;
3686	>	patom->SetFormalCharge(chg1);
3687	>	if(chg2>0)
3688	>	--chg2;
3689	>	else
3690	>	++chg2;
3691	>	pNbratom->SetFormalCharge(chg2);
3692	>	pbond->SetBO(pbond->GetBO()+1);
3693	>	}
3694	>	}
3695	>	}
3696	>	return true;
3697	>	}
3698
3699	<	OBAtom *OBMol::BeginAtom(vector<OBNodeBase*>::iterator &i)
3700	<	{
3699	>	OBAtom *OBMol::BeginAtom(vector<OBNodeBase*>::iterator &i)
3700	>	{
3701		i = _vatom.begin();
3702		return((i == _vatom.end()) ? (OBAtom*)NULL : (OBAtom*)*i);
3703	<	}
3703	>	}
3704
3705	<	OBAtom *OBMol::NextAtom(vector<OBNodeBase*>::iterator &i)
3706	<	{
3705	>	OBAtom *OBMol::NextAtom(vector<OBNodeBase*>::iterator &i)
3706	>	{
3707		i++;
3708		return((i == _vatom.end()) ? (OBAtom*)NULL : (OBAtom*)*i);
3709	<	}
3709	>	}
3710
3711	<	OBBond *OBMol::BeginBond(vector<OBEdgeBase*>::iterator &i)
3712	<	{
3711	>	OBBond *OBMol::BeginBond(vector<OBEdgeBase*>::iterator &i)
3712	>	{
3713		i = _vbond.begin();
3714		return((i == _vbond.end()) ? (OBBond*)NULL : (OBBond*)*i);
3715	<	}
3715	>	}
3716
3717	<	OBBond *OBMol::NextBond(vector<OBEdgeBase*>::iterator &i)
3718	<	{
3717	>	OBBond *OBMol::NextBond(vector<OBEdgeBase*>::iterator &i)
3718	>	{
3719		i++;
3720		return((i == _vbond.end()) ? (OBBond*)NULL : (OBBond*)*i);
3721	<	}
3721	>	}
3722
3723		} // end namespace OpenBabel
3724

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines