1 // -*- C++ -*-
2
3 /*
4 * Gnome Chemistry Utils
5 * libs/gcu/chain.cc
6 *
7 * Copyright (C) 2001-2011 Jean Bréfort <jean.brefort@normalesup.org>
8 *
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License as
11 * published by the Free Software Foundation; either version 3 of the
12 * License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, write to the Free Software
21 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
22 * USA
23 */
24
25 #include "config.h"
26 #include "atom.h"
27 #include "bond.h"
28 #include "chain.h"
29 #include "cycle.h"
30 #include "molecule.h"
31 #include "document.h"
32 #include <glib/gi18n-lib.h>
33
34 using namespace std;
35
36 namespace gcu {
37
Chain(Bond * pBond,Atom * pAtom,TypeId Type)38 Chain::Chain (Bond* pBond, Atom* pAtom, TypeId Type): Object (Type)
39 {
40 Atom *pAtom0;
41 if (pAtom)
42 pAtom0 = (Atom*) pBond->GetAtom (pAtom);
43 else {
44 pAtom0 = (Atom*) pBond->GetAtom (1);
45 pAtom = (Atom*) pBond->GetAtom (0);
46 }
47 m_Bonds[pAtom].fwd = pBond;
48 m_Bonds[pAtom0].rev = pBond;
49 }
50
Chain(Molecule * molecule,Atom * pAtom,TypeId Type)51 Chain::Chain (Molecule* molecule, Atom* pAtom, TypeId Type): Object (Type)
52 {
53 m_Molecule = molecule;
54 if (pAtom) {
55 FindCycles (pAtom);
56 }
57 }
58
59 /*
60 * Add a bond in an existing molecule and update cycles
61 * Implementation might have to be changed
62 */
Chain(Molecule * molecule,Bond * pBond,TypeId Type)63 Chain::Chain (Molecule* molecule, Bond* pBond, TypeId Type): Object (Type)
64 {
65 m_Molecule = molecule;
66 if (pBond) {
67 Atom *pAtom0, *pAtom;
68 pAtom0 = (Atom*) pBond->GetAtom (0);
69 m_Bonds[pAtom0].fwd = pBond;
70 pAtom = (Atom*) pBond->GetAtom (1);
71 m_Bonds[pAtom].rev = pBond;
72 map<gcu::Atom*, gcu::Bond*>::iterator i;
73 Bond* pBond0 = (Bond*) pAtom->GetFirstBond (i);
74 while (pBond0) {
75 if ((pBond0 != pBond) && FindCycle (pAtom, pBond0))
76 break;
77 pBond0 = (Bond*) pAtom->GetNextBond (i);
78 }
79 }
80 }
81
~Chain()82 Chain::~Chain ()
83 {
84 m_Bonds.clear ();
85 }
86
FindCycle(Atom * pAtom,Bond * pBond)87 bool Chain::FindCycle (Atom* pAtom, Bond* pBond)
88 {
89 Atom* pAtom1 = (Atom*) pBond->GetAtom (pAtom);
90 if (m_Bonds[pAtom1].fwd != NULL) {
91 if (m_Bonds[pAtom1].rev != NULL)
92 return false;
93 Cycle* pCycle = new Cycle (m_Molecule);
94 pCycle->m_Bonds[pAtom1].rev = pBond;
95 pCycle->m_Bonds[pAtom1].fwd = m_Bonds[pAtom1].fwd;
96 pCycle->m_Bonds[pAtom].fwd = pBond;
97 pCycle->m_Bonds[pAtom].rev= m_Bonds[pAtom].rev;
98 m_Bonds[pAtom].rev->AddCycle (pCycle);
99 pBond->AddCycle (pCycle);
100 while (pBond = pCycle->m_Bonds[pAtom1].fwd, pAtom1 = (Atom*) pBond->GetAtom (pAtom1), pAtom != pAtom1) {
101 pCycle->m_Bonds[pAtom1].rev = pBond;
102 pCycle->m_Bonds[pAtom1].fwd = m_Bonds[pAtom1].fwd;
103 pBond->AddCycle (pCycle);
104 }
105 pCycle->Simplify (); //to reduce size of fused cycles
106 m_Molecule->m_Cycles.push_back (pCycle);
107 return true;
108 }
109 m_Bonds[pAtom].fwd = pBond;
110 m_Bonds[pAtom1].rev = pBond;
111 map<gcu::Atom*, gcu::Bond*>::iterator i;
112 Bond* pBond1 = (Bond*) pAtom1->GetFirstBond (i);
113 while (pBond1) {
114 if ((pBond1 != pBond) && FindCycle (pAtom1, pBond1))
115 return true;
116 pBond1 = (Bond*) pAtom1->GetNextBond (i);
117 }
118 m_Bonds[pAtom].fwd = NULL;
119 m_Bonds.erase (pAtom1);
120 return false;
121 }
122
FindCycles(Atom * pAtom)123 void Chain::FindCycles (Atom* pAtom)
124 {
125 map<gcu::Atom*, gcu::Bond*>::iterator i;
126 Bond* pBond = (Bond*) pAtom->GetFirstBond (i);
127 Atom* pAtom0;
128 Molecule *mol;
129 while (pBond != NULL) {
130 m_Bonds[pAtom].fwd = pBond;
131 pAtom0 = (Atom*) pBond->GetAtom (pAtom);
132 if ((mol = static_cast < Molecule * > (pBond->GetMolecule ())) != m_Molecule) {
133 if (mol)
134 mol->ClearCycles ();
135 m_Molecule->AddChild (pBond);
136 }
137 if ((pAtom0)->GetMolecule () != m_Molecule) {
138 if (pAtom0->GetMolecule () != m_Molecule)
139 m_Molecule->AddChild (pAtom0);
140 FindCycles (pAtom0);
141 } else {
142 if (m_Bonds[pAtom0].fwd != NULL) {
143 Bond* pBond0 = m_Bonds[pAtom0].fwd;
144 if (pAtom != pBond0->GetAtom (pAtom0)) {
145 //Cycle found
146 Cycle* pCycle = new Cycle (m_Molecule);
147 pCycle->m_Bonds[pAtom0].rev = pBond;
148 pCycle->m_Bonds[pAtom0].fwd = pBond0;
149 pBond0->AddCycle (pCycle);
150 while (pAtom != pAtom0) {
151 pAtom0 = (Atom*) pBond0->GetAtom (pAtom0);
152 pCycle->m_Bonds[pAtom0].rev = pBond0;
153 pBond0 = m_Bonds[pAtom0].fwd;
154 pCycle->m_Bonds[pAtom0].fwd = pBond0;
155 pBond0->AddCycle (pCycle);
156 }
157 pCycle->Simplify (); //to reduce size of fused cycles
158 m_Molecule->m_Cycles.push_back (pCycle);
159 }
160 }
161 }
162 pBond = (Bond*) pAtom->GetNextBond (i);
163 }
164 m_Bonds.erase (pAtom);
165 }
166
Reverse()167 void Chain::Reverse ()
168 {
169 map<Atom*, ChainElt>::iterator i, end = m_Bonds.end ();
170 Bond* pBond;
171 for (i = m_Bonds.begin (); i != end; i++) {
172 pBond = (*i).second.fwd;
173 (*i).second.fwd = (*i).second.rev;
174 (*i).second.rev = pBond;
175 }
176 }
177
Erase(Atom * pAtom1,Atom * pAtom2)178 void Chain::Erase (Atom* pAtom1, Atom* pAtom2)
179 {
180 //This function is not safe
181 Atom *pAtom = (Atom*) m_Bonds[pAtom1].fwd->GetAtom (pAtom1), *pAtom0;
182 m_Bonds[pAtom1].fwd = NULL;
183 while (pAtom != pAtom2) {
184 pAtom = (Atom*) m_Bonds[pAtom].fwd->GetAtom (pAtom0 = pAtom);
185 m_Bonds.erase (pAtom0);
186 }
187 m_Bonds[pAtom2].rev = NULL;
188 }
189
Insert(Atom * pAtom1,Atom * pAtom2,Chain & chain)190 void Chain::Insert (Atom* pAtom1, Atom* pAtom2, Chain& chain)
191 {
192 //This function is not safe
193 m_Bonds[pAtom1].fwd = chain.m_Bonds[pAtom1].fwd;
194 Atom *pAtom = (Atom*) m_Bonds[pAtom1].fwd->GetAtom (pAtom1);
195 while (pAtom != pAtom2) {
196 m_Bonds[pAtom] = chain.m_Bonds[pAtom];
197 pAtom = (Atom*) m_Bonds[pAtom].fwd->GetAtom (pAtom);
198 }
199 m_Bonds[pAtom2].rev = chain.m_Bonds[pAtom2].rev;
200 }
201
Extract(Atom * pAtom1,Atom * pAtom2,Chain & chain)202 void Chain::Extract (Atom* pAtom1, Atom* pAtom2, Chain& chain)
203 {
204 chain.m_Bonds.clear ();
205 if (m_Bonds[pAtom1].fwd == NULL) {
206 if (m_Bonds[pAtom1].rev == NULL)
207 m_Bonds.erase (pAtom1); //pAtom1 is not in the chain
208 return;
209 }
210 chain.m_Bonds[pAtom1].fwd = m_Bonds[pAtom1].fwd;
211 chain.m_Bonds[pAtom1].rev = NULL;
212 Atom *pAtom = (Atom*) chain.m_Bonds[pAtom1].fwd->GetAtom (pAtom1);
213 while (pAtom != pAtom2) {
214 chain.m_Bonds[pAtom] = m_Bonds[pAtom];
215 if (m_Bonds[pAtom].fwd == NULL)
216 return; //Chain never reach pAtom2
217 pAtom = (Atom*)m_Bonds[pAtom].fwd->GetAtom (pAtom);
218 }
219 chain.m_Bonds[pAtom2].rev = m_Bonds[pAtom2].rev;
220 chain.m_Bonds[pAtom2].fwd = NULL;
221 }
222
GetUnsaturations()223 unsigned Chain::GetUnsaturations ()
224 {
225 unsigned n = 0;
226 std::map<Atom*, ChainElt>::iterator i, end = m_Bonds.end ();
227 for (i = m_Bonds.begin (); i != end; i++)
228 if ((*i).second.fwd && ((*i).second.fwd->GetOrder () > 1))
229 n += 1;
230 return n;
231 }
232
GetHeteroatoms()233 unsigned Chain::GetHeteroatoms ()
234 {
235 unsigned n = 0;
236 std::map<Atom*, ChainElt>::iterator i, end = m_Bonds.end ();
237 for (i = m_Bonds.begin (); i != end; i++)
238 if ((*i).first->GetZ () != 6)
239 n += 1;
240 return n;
241 }
242
AddBond(Atom * start,Atom * end)243 void Chain::AddBond (Atom* start, Atom* end)
244 {
245 Bond* pBond = (Bond*) start->GetBond (end);
246 m_Bonds[start].fwd = pBond;
247 m_Bonds[end].rev = pBond;
248 }
249
Contains(Atom * pAtom)250 bool Chain::Contains (Atom* pAtom)
251 {
252 if ((m_Bonds[pAtom].fwd == NULL) && (m_Bonds[pAtom].rev == NULL)) {
253 m_Bonds.erase (pAtom);
254 return false;
255 }
256 return true;
257 }
258
Contains(Bond * pBond)259 bool Chain::Contains (Bond* pBond)
260 {
261 Atom *pAtom = (Atom*) pBond->GetAtom (0);
262 if ((m_Bonds[pAtom].fwd == NULL) && (m_Bonds[pAtom].rev == NULL)) {
263 m_Bonds.erase (pAtom);
264 return false;
265 }
266 if ((m_Bonds[pAtom].fwd == pBond) && (m_Bonds[pAtom].rev == pBond))
267 return true;
268 return false;
269 }
270
GetLength()271 unsigned Chain::GetLength ()
272 {
273 unsigned n = 0;
274 std::map<Atom*, ChainElt>::iterator i, end = m_Bonds.end ();
275 for (i = m_Bonds.begin(); i != end; i++)
276 if ((*i).second.fwd)
277 n++;
278 return n;
279 }
280
GetMeanBondLength()281 double Chain::GetMeanBondLength ()
282 {
283 unsigned n = 0;
284 double l = 0;
285 std::map<Atom*, ChainElt>::iterator i, end = m_Bonds.end ();
286 for (i = m_Bonds.begin (); i != end; i++)
287 if ((*i).second.fwd) {
288 n++;
289 l += (*i).second.fwd->Get2DLength ();
290 }
291 return l / n;
292 }
293
GetFirstAtom()294 Atom* Chain::GetFirstAtom ()
295 {
296 std::map < Atom*, ChainElt >::iterator i = m_Bonds.begin ();
297 if (GetType () == CycleType)
298 return (*i).first;
299 Atom *res = (*i).first, *cur = res;
300 while (cur) {
301 cur = m_Bonds[res].rev->GetAtom (res);
302 if (cur == NULL)
303 break;
304 res = cur;
305 }
306 return res;
307 }
308
GetNextAtom(Atom * pAtom)309 Atom* Chain::GetNextAtom (Atom* pAtom)
310 {
311 return (Atom*) m_Bonds[pAtom].fwd->GetAtom (pAtom);
312 }
313
Name()314 std::string Chain::Name ()
315 {
316 return _("Chain");
317 }
318
BuildLength(unsigned * cycle_size,unsigned * cycle_pos)319 unsigned Chain::BuildLength (unsigned *cycle_size, unsigned *cycle_pos)
320 {
321 // searching the longest chain starting from the first bond, and stopping at any cycle
322 // cycle_size is the largest encountered cycle
323 unsigned length = 0;
324 unsigned nb;
325 unsigned max_cycle_size = 0;
326 unsigned min_cycle_pos = 0;
327 // searching from there
328 std::map < Atom *, ChainElt >::iterator i, end = m_Bonds.end ();
329 std::map < Atom *, Bond * >::iterator b;
330 Bond *bond, *last_bond = NULL;
331 Atom *atom = NULL;
332 for (i = m_Bonds.begin(); i != end; i++) {
333 if ((*i).second.fwd)
334 length++;
335 else {
336 atom = (*i).first;
337 last_bond = (*i).second.rev;
338 }
339 }
340 while (atom && (nb = atom->GetBondsNumber ()) != 1)
341 switch (atom->GetBondsNumber ()) {
342 case 2:
343 bond = atom->GetFirstBond (b);
344 if (bond == last_bond)
345 bond = atom->GetNextBond (b);
346 m_Bonds[atom].fwd = bond;
347 atom = bond->GetAtom (atom);
348 m_Bonds[atom].rev = bond;
349 last_bond = bond;
350 length++;
351 break;
352 default: {
353 // we reached a ramification
354 // exploring each chain
355 unsigned rcycle_size = 0, rcycle_pos = 0, rlength = 0;
356 Chain *longchain = NULL;
357 for (bond = atom->GetFirstBond (b); bond; bond = atom->GetNextBond (b)) {
358 if (bond == last_bond)
359 continue;
360 if (bond->IsCyclic ()) {
361 if (min_cycle_pos == 0)
362 min_cycle_pos = length;
363 if (min_cycle_pos == length) {
364 std::list < Cycle * >::iterator c;
365 Cycle * cycle = bond->GetFirstCycle (c, NULL);
366 while (cycle) {
367 if (cycle->GetLength () > max_cycle_size)
368 max_cycle_size = cycle->GetLength ();
369 cycle = bond->GetNextCycle (c, NULL);
370 }
371 }
372 continue;
373 }
374 // if we are there, we have normal chain that we explore
375 unsigned cycle_pos = 0, cycle_size = 0, length;
376 Chain *chain = new Chain (bond, atom);
377 length = chain->BuildLength (&cycle_size, &cycle_pos);
378 if (length > rlength) {
379 if (longchain)
380 delete longchain;
381 longchain = chain;
382 if ((cycle_pos = !0) && (cycle_pos < rcycle_pos || rcycle_pos == 0)) {
383 rcycle_pos = cycle_pos;
384 rcycle_size = cycle_size;
385 }
386 if (cycle_pos == rcycle_pos && cycle_size > rcycle_size)
387 rcycle_size = cycle_size;
388 }
389 }
390 if (longchain) {
391 Append (*longchain);
392 if (min_cycle_pos == 0) {
393 min_cycle_pos = rcycle_pos;
394 max_cycle_size = rcycle_size;
395 }
396 length += rlength;
397 delete longchain;
398 }
399 atom = NULL;
400 break;
401 }
402 }
403 //and now ending
404 if (cycle_size)
405 *cycle_size = max_cycle_size;
406 if (cycle_pos)
407 *cycle_pos = min_cycle_pos;
408 return length;
409 }
410
Append(Chain & chain)411 void Chain::Append (Chain& chain)
412 {
413 std::map < Atom *, ChainElt >::iterator i;
414 for (i = m_Bonds.begin (); (*i).second.fwd; i++);
415 Atom *atom = (*i).first;
416 if (chain.m_Bonds.find (atom) == chain.m_Bonds.end ())
417 return;
418 m_Bonds[atom].fwd = chain.m_Bonds[atom].fwd;
419 atom = m_Bonds[atom].fwd->GetAtom (atom);
420 while (chain.m_Bonds[atom].fwd) {
421 m_Bonds[atom] = chain.m_Bonds[atom];
422 atom = m_Bonds[atom].fwd->GetAtom (atom);
423 }
424
425 }
426
427 } // namespace gcp
428