| 1 | tim | 741 | /********************************************************************* | 
| 2 |  |  | chiral.cpp - Detect chiral atoms and molecules. | 
| 3 |  |  |  | 
| 4 |  |  | Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc. | 
| 5 |  |  | Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison | 
| 6 |  |  |  | 
| 7 |  |  | This file is part of the Open Babel project. | 
| 8 |  |  | For more information, see <http://openbabel.sourceforge.net/> | 
| 9 |  |  |  | 
| 10 |  |  | This program is free software; you can redistribute it and/or modify | 
| 11 |  |  | it under the terms of the GNU General Public License as published by | 
| 12 |  |  | the Free Software Foundation version 2 of the License. | 
| 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 |  |  |  | 
| 20 |  |  | #include <list> | 
| 21 |  |  |  | 
| 22 |  |  | #include "mol.hpp" | 
| 23 |  |  | #include "obutil.hpp" | 
| 24 |  |  | #include "matrix.hpp" | 
| 25 |  |  | #include "chiral.hpp" | 
| 26 |  |  | #include "matrix3x3.hpp" | 
| 27 |  |  |  | 
| 28 |  |  | using namespace std; | 
| 29 |  |  | namespace OpenBabel | 
| 30 |  |  | { | 
| 31 |  |  |  | 
| 32 |  |  | //! Sets atom->IsChiral() to true for chiral atoms | 
| 33 |  |  | void OBMol::FindChiralCenters() | 
| 34 |  |  | { | 
| 35 |  |  | if (HasChiralityPerceived()) | 
| 36 |  |  | return; | 
| 37 |  |  | SetChiralityPerceived(); | 
| 38 |  |  |  | 
| 39 |  |  | obErrorLog.ThrowError(__FUNCTION__, | 
| 40 |  |  | "Ran OpenBabel::FindChiralCenters", obAuditMsg); | 
| 41 |  |  |  | 
| 42 |  |  | //do quick test to see if there are any possible chiral centers | 
| 43 |  |  | bool mayHaveChiralCenter=false; | 
| 44 |  |  | OBAtom *atom,*nbr; | 
| 45 |  |  | vector<OBNodeBase*>::iterator i; | 
| 46 |  |  | for (atom = BeginAtom(i);atom;atom = NextAtom(i)) | 
| 47 |  |  | if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3) | 
| 48 |  |  | { | 
| 49 |  |  | mayHaveChiralCenter=true; | 
| 50 |  |  | break; | 
| 51 |  |  | } | 
| 52 |  |  |  | 
| 53 |  |  | if (!mayHaveChiralCenter) | 
| 54 |  |  | return; | 
| 55 |  |  |  | 
| 56 |  |  | OBBond *bond; | 
| 57 |  |  | vector<OBEdgeBase*>::iterator j; | 
| 58 |  |  | for (bond = BeginBond(j);bond;bond = NextBond(j)) | 
| 59 |  |  | if (bond->IsWedge() || bond->IsHash()) | 
| 60 |  |  | (bond->GetBeginAtom())->SetChiral(); | 
| 61 |  |  |  | 
| 62 |  |  | #define INTMETHOD | 
| 63 |  |  |  | 
| 64 |  |  | #ifdef INTMETHOD | 
| 65 |  |  |  | 
| 66 |  |  | vector<unsigned int> vgid; | 
| 67 |  |  | GetGIDVector(vgid); | 
| 68 |  |  | vector<unsigned int> tlist; | 
| 69 |  |  | vector<unsigned int>::iterator k; | 
| 70 |  |  | #else //use Golender floating point method | 
| 71 |  |  |  | 
| 72 |  |  | vector<double> gp; | 
| 73 |  |  | GraphPotentials(*this,gp); | 
| 74 |  |  | vector<double> tlist; | 
| 75 |  |  | vector<double>::iterator k; | 
| 76 |  |  | #endif | 
| 77 |  |  |  | 
| 78 |  |  | bool ischiral; | 
| 79 |  |  | for (atom = BeginAtom(i);atom;atom = NextAtom(i)) | 
| 80 |  |  | if (atom->GetHyb() == 3 && atom->GetHvyValence() >= 3 && !atom->IsChiral()) | 
| 81 |  |  | { | 
| 82 |  |  | tlist.clear(); | 
| 83 |  |  | ischiral = true; | 
| 84 |  |  |  | 
| 85 |  |  | for (nbr = atom->BeginNbrAtom(j);nbr;nbr = atom->NextNbrAtom(j)) | 
| 86 |  |  | { | 
| 87 |  |  | for (k = tlist.begin();k != tlist.end();k++) | 
| 88 |  |  | #ifdef INTMETHOD | 
| 89 |  |  |  | 
| 90 |  |  | if (vgid[nbr->GetIdx()-1] == *k) | 
| 91 |  |  | #else | 
| 92 |  |  |  | 
| 93 |  |  | if (fabs(gp[nbr->GetIdx()]-*k) < 0.001) | 
| 94 |  |  | #endif | 
| 95 |  |  |  | 
| 96 |  |  | ischiral = false; | 
| 97 |  |  |  | 
| 98 |  |  | #ifdef INTMETHOD | 
| 99 |  |  |  | 
| 100 |  |  | if (ischiral) | 
| 101 |  |  | tlist.push_back(vgid[nbr->GetIdx()-1]); | 
| 102 |  |  | #else | 
| 103 |  |  |  | 
| 104 |  |  | if (ischiral) | 
| 105 |  |  | tlist.push_back(gp[nbr->GetIdx()]); | 
| 106 |  |  | #endif | 
| 107 |  |  |  | 
| 108 |  |  | else | 
| 109 |  |  | break; | 
| 110 |  |  | } | 
| 111 |  |  |  | 
| 112 |  |  | if (ischiral) | 
| 113 |  |  | atom->SetChiral(); | 
| 114 |  |  | } | 
| 115 |  |  | } | 
| 116 |  |  |  | 
| 117 |  |  | // Seems to make a vector chirality become filled with array of +/- 1 for chiral atoms. | 
| 118 |  |  | void GetChirality(OBMol &mol, std::vector<int> &chirality) | 
| 119 |  |  | { | 
| 120 |  |  | chirality.resize(mol.NumAtoms()+1); | 
| 121 |  |  | fill(chirality.begin(),chirality.end(),0); | 
| 122 |  |  |  | 
| 123 |  |  | OBAtom *atom; | 
| 124 |  |  | vector<OBNodeBase*>::iterator i; | 
| 125 |  |  | for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i)) | 
| 126 |  |  | if (atom->IsChiral()) | 
| 127 |  |  | { | 
| 128 |  |  | if (!atom->HasChiralVolume()) | 
| 129 |  |  | { | 
| 130 |  |  | double sv = CalcSignedVolume(mol,atom); | 
| 131 |  |  | if (sv < 0.0) | 
| 132 |  |  | { | 
| 133 |  |  | chirality[atom->GetIdx()-1] = -1; | 
| 134 |  |  | atom->SetNegativeStereo(); | 
| 135 |  |  | } | 
| 136 |  |  | else if (sv > 0.0) | 
| 137 |  |  | { | 
| 138 |  |  | chirality[atom->GetIdx()-1] = 1; | 
| 139 |  |  | atom->SetPositiveStereo(); | 
| 140 |  |  | } | 
| 141 |  |  | } | 
| 142 |  |  | else // already calculated signed volume (e.g., imported from somewhere) | 
| 143 |  |  | { | 
| 144 |  |  | if (atom ->IsPositiveStereo()) | 
| 145 |  |  | chirality[atom->GetIdx()-1] = 1; | 
| 146 |  |  | else | 
| 147 |  |  | chirality[atom->GetIdx()-1] = -1; | 
| 148 |  |  | } | 
| 149 |  |  | } | 
| 150 |  |  | } | 
| 151 |  |  |  | 
| 152 |  |  | int GetParity4Ref(vector<unsigned int> pref) | 
| 153 |  |  | { | 
| 154 |  |  | if(pref.size()!=4)return(-1); // should be given a vector of size 4. | 
| 155 |  |  | int parity=0; | 
| 156 |  |  | for (int i=0;i<3;i++) // do the bubble sort this many times | 
| 157 |  |  | { | 
| 158 |  |  | for(int j=0;j<3;j++) // iterate across the array 4th element has no | 
| 159 |  |  | {                    // right hand neighbour so no need to sort | 
| 160 |  |  | if (pref[j+1] < pref[j]) // compare the two neighbors | 
| 161 |  |  | { | 
| 162 |  |  | unsigned int tmp = pref[j];        // swap a[j] and a[j+1] | 
| 163 |  |  | pref[j] = pref[j+1]; | 
| 164 |  |  | pref[j+1] = tmp; | 
| 165 |  |  | parity++; // parity odd will invert stereochem | 
| 166 |  |  | } | 
| 167 |  |  | } | 
| 168 |  |  | } // End Bubble Sort | 
| 169 |  |  | return(parity%2); | 
| 170 |  |  | } | 
| 171 |  |  |  | 
| 172 |  |  | bool CorrectChirality(OBMol &mol, OBAtom *atm, atomreftype i, atomreftype o) | 
| 173 |  |  | { | 
| 174 |  |  | if (!atm->HasChiralitySpecified()) // if no chirality defined can't do any more for 0D | 
| 175 |  |  | return(false); | 
| 176 |  |  |  | 
| 177 |  |  | int parityI=0,parityO=0; | 
| 178 |  |  | OBChiralData* cd=(OBChiralData*)atm->GetData(OBGenericDataType::ChiralData); | 
| 179 |  |  | if ((cd->GetAtom4Refs(input)).size()!=4)return(false); // must have 4 refs | 
| 180 |  |  | parityI=GetParity4Ref(cd->GetAtom4Refs(i)); // Gets Atom4Refs used to define the chirality | 
| 181 |  |  | parityO=GetParity4Ref(cd->GetAtom4Refs(o));//GetsOutput parity. | 
| 182 |  |  | /* switch (CHTYPE) | 
| 183 |  |  | { | 
| 184 |  |  | case SMILES: // SMILES always uses 1234 atom refs | 
| 185 |  |  | parityO=0; // if parityO==parityI then clockwise in input = clockwise in output | 
| 186 |  |  | break; | 
| 187 |  |  | case MOLV3000: // MOLV3000 uses 1234 unless an H then 123H | 
| 188 |  |  | if (atm->GetHvyValence()==3) | 
| 189 |  |  | { | 
| 190 |  |  | OBAtom *nbr; | 
| 191 |  |  | int Hid=1000;// max Atom ID +1 should be used here | 
| 192 |  |  | vector<int> nbr_atms; | 
| 193 |  |  | vector<OBEdgeBase*>::iterator i; | 
| 194 |  |  | for (nbr = atm->BeginNbrAtom(i);nbr;nbr = atm->NextNbrAtom(i)) | 
| 195 |  |  | { | 
| 196 |  |  | if (nbr->IsHydrogen()){Hid=nbr->GetIdx();continue;} | 
| 197 |  |  | nbr_atms.push_back(nbr->GetIdx()); | 
| 198 |  |  | } | 
| 199 |  |  | sort(nbr_atms.begin(),nbr_atms.end()); | 
| 200 |  |  | nbr_atms.push_back(Hid); | 
| 201 |  |  | int tmp[4]; | 
| 202 |  |  | for(int i=0;i<4;i++){tmp[i]=nbr_atms[i];} | 
| 203 |  |  | parityO=GetParity4Ref(tmp); | 
| 204 |  |  | } | 
| 205 |  |  | else if (atm->GetHvyValence()==4) | 
| 206 |  |  | parityO=0; | 
| 207 |  |  | break; | 
| 208 |  |  | default: | 
| 209 |  |  | parityO=0; | 
| 210 |  |  | }*/ | 
| 211 |  |  | if (parityO==parityI) | 
| 212 |  |  | {//cout << "Parity is the same"<<endl; | 
| 213 |  |  | return(true); | 
| 214 |  |  | } | 
| 215 |  |  | else if(parityO!=parityI) // Need to invert the Chirality which has been set | 
| 216 |  |  | { //cout << "Parity is Opposite"<<endl; | 
| 217 |  |  | if (atm->IsClockwise()) | 
| 218 |  |  | {atm->UnsetStereo();atm->SetAntiClockwiseStereo();} | 
| 219 |  |  | else if (atm->IsAntiClockwise()) | 
| 220 |  |  | {atm->UnsetStereo();atm->SetClockwiseStereo();} | 
| 221 |  |  | else | 
| 222 |  |  | return(false); | 
| 223 |  |  | return(true); | 
| 224 |  |  | } | 
| 225 |  |  | return false; | 
| 226 |  |  | } | 
| 227 |  |  |  | 
| 228 |  |  | //! Calculate the signed volume for an atom.  If the atom has a valence of 3 | 
| 229 |  |  | //! the coordinates of an attached hydrogen are calculated | 
| 230 |  |  | //! Puts attached Hydrogen last at the moment, like mol V3000 format. | 
| 231 |  |  | double CalcSignedVolume(OBMol &mol,OBAtom *atm) | 
| 232 |  |  | { | 
| 233 |  |  | vector3 tmp_crd; | 
| 234 |  |  | vector<unsigned int> nbr_atms; | 
| 235 |  |  | vector<vector3> nbr_crds; | 
| 236 |  |  | bool use_central_atom = false,is2D=false; | 
| 237 |  |  | double hbrad = etab.CorrectedBondRad(1,0); | 
| 238 |  |  |  | 
| 239 |  |  | if (!mol.Has3D()) //give peudo Z coords if mol is 2D | 
| 240 |  |  | { | 
| 241 |  |  | vector3 v,vz(0.0,0.0,1.0); | 
| 242 |  |  | is2D = true; | 
| 243 |  |  | OBAtom *nbr; | 
| 244 |  |  | OBBond *bond; | 
| 245 |  |  | vector<OBEdgeBase*>::iterator i; | 
| 246 |  |  | for (bond = atm->BeginBond(i);bond;bond = atm->NextBond(i)) | 
| 247 |  |  | { | 
| 248 |  |  | nbr = bond->GetEndAtom(); | 
| 249 |  |  | if (nbr != atm) | 
| 250 |  |  | { | 
| 251 |  |  | v = nbr->GetVector(); | 
| 252 |  |  | if (bond->IsWedge()) | 
| 253 |  |  | v += vz; | 
| 254 |  |  | else | 
| 255 |  |  | if (bond->IsHash()) | 
| 256 |  |  | v -= vz; | 
| 257 |  |  |  | 
| 258 |  |  | nbr->SetVector(v); | 
| 259 |  |  | } | 
| 260 |  |  | else | 
| 261 |  |  | { | 
| 262 |  |  | nbr = bond->GetBeginAtom(); | 
| 263 |  |  | v = nbr->GetVector(); | 
| 264 |  |  | if (bond->IsWedge()) | 
| 265 |  |  | v -= vz; | 
| 266 |  |  | else | 
| 267 |  |  | if (bond->IsHash()) | 
| 268 |  |  | v += vz; | 
| 269 |  |  |  | 
| 270 |  |  | nbr->SetVector(v); | 
| 271 |  |  | } | 
| 272 |  |  | } | 
| 273 |  |  | } | 
| 274 |  |  |  | 
| 275 |  |  | if (atm->GetHvyValence() < 3) | 
| 276 |  |  | { | 
| 277 |  |  | #ifdef HAVE_SSTREAM | 
| 278 |  |  | stringstream errorMsg; | 
| 279 |  |  | #else | 
| 280 |  |  | strstream errorMsg; | 
| 281 |  |  | #endif | 
| 282 |  |  | errorMsg << "Cannot calculate a signed volume for an atom with a heavy atom valence of " << atm->GetHvyValence() << endl; | 
| 283 |  |  | obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obInfo); | 
| 284 |  |  | return(0.0); | 
| 285 |  |  | } | 
| 286 |  |  |  | 
| 287 |  |  | // Create a vector with the coordinates of the neighbor atoms | 
| 288 |  |  | // Also make a vector with Atom IDs | 
| 289 |  |  | OBAtom *nbr; | 
| 290 |  |  | vector<OBEdgeBase*>::iterator bint; | 
| 291 |  |  | for (nbr = atm->BeginNbrAtom(bint);nbr;nbr = atm->NextNbrAtom(bint)) | 
| 292 |  |  | { | 
| 293 |  |  | nbr_atms.push_back(nbr->GetIdx()); | 
| 294 |  |  | } | 
| 295 |  |  | // sort the neighbor atoms to insure a consistent ordering | 
| 296 |  |  | sort(nbr_atms.begin(),nbr_atms.end()); | 
| 297 |  |  | for (unsigned int i = 0; i < nbr_atms.size(); i++) | 
| 298 |  |  | { | 
| 299 |  |  | OBAtom *tmp_atm = mol.GetAtom(nbr_atms[i]); | 
| 300 |  |  | nbr_crds.push_back(tmp_atm->GetVector()); | 
| 301 |  |  | } | 
| 302 |  |  | /* | 
| 303 |  |  | // If we have three heavy atoms we need to calculate the position of the fourth | 
| 304 |  |  | if (atm->GetHvyValence() == 3) | 
| 305 |  |  | { | 
| 306 |  |  | double bondlen = hbrad+etab.CorrectedBondRad(atm->GetAtomicNum(),atm->GetHyb()); | 
| 307 |  |  | atm->GetNewBondVector(tmp_crd,bondlen); | 
| 308 |  |  | nbr_crds.push_back(tmp_crd); | 
| 309 |  |  | } | 
| 310 |  |  | */ | 
| 311 |  |  | for(int j=0;j < nbr_crds.size();j++) // Checks for a neighbour having 0 co-ords (added hydrogen etc) | 
| 312 |  |  | { | 
| 313 |  |  | if (nbr_crds[j]==0 && use_central_atom==false)use_central_atom=true; | 
| 314 |  |  | else if (nbr_crds[j]==0) | 
| 315 |  |  | { | 
| 316 |  |  | obErrorLog.ThrowError(__FUNCTION__, "More than 2 neighbours have 0 co-ords when attempting 3D chiral calculation", obInfo); | 
| 317 |  |  | } | 
| 318 |  |  | } | 
| 319 |  |  |  | 
| 320 |  |  | // If we have three heavy atoms we can use the chiral center atom itself for the fourth | 
| 321 |  |  | // will always give same sign (for tetrahedron), magnitude will be smaller. | 
| 322 |  |  | if(nbr_atms.size()==3 || use_central_atom==true) | 
| 323 |  |  | { | 
| 324 |  |  | nbr_crds.push_back(atm->GetVector()); | 
| 325 |  |  | nbr_atms.push_back(mol.NumAtoms()+1); // meed to add largest number on end to work | 
| 326 |  |  | } | 
| 327 |  |  | OBChiralData* cd=(OBChiralData*)atm->GetData(OBGenericDataType::ChiralData); //Set the output atom4refs to the ones used | 
| 328 |  |  | if(cd==NULL) | 
| 329 |  |  | { | 
| 330 |  |  | cd = new OBChiralData; | 
| 331 |  |  | atm->SetData(cd); | 
| 332 |  |  | } | 
| 333 |  |  | cd->SetAtom4Refs(nbr_atms,calcvolume); | 
| 334 |  |  |  | 
| 335 |  |  | //re-zero psuedo-coords | 
| 336 |  |  | if (is2D) | 
| 337 |  |  | { | 
| 338 |  |  | vector3 v; | 
| 339 |  |  | OBAtom *atom; | 
| 340 |  |  | vector<OBNodeBase*>::iterator k; | 
| 341 |  |  | for (atom = mol.BeginAtom(k);atom;atom = mol.NextAtom(k)) | 
| 342 |  |  | { | 
| 343 |  |  | v = atom->GetVector(); | 
| 344 |  |  | v.SetZ(0.0); | 
| 345 |  |  | atom->SetVector(v); | 
| 346 |  |  | } | 
| 347 |  |  | } | 
| 348 |  |  |  | 
| 349 |  |  | return(signed_volume(nbr_crds[0],nbr_crds[1],nbr_crds[2],nbr_crds[3])); | 
| 350 |  |  | } | 
| 351 |  |  |  | 
| 352 |  |  | //! Calculate a signed volume given a set of 4 coordinates | 
| 353 |  |  | double signed_volume(const vector3 &a, const vector3 &b, const vector3 &c, const vector3 &d) | 
| 354 |  |  | { | 
| 355 |  |  | vector3 A,B,C; | 
| 356 |  |  | A = b-a; | 
| 357 |  |  | B = c-a; | 
| 358 |  |  | C = d-a; | 
| 359 |  |  | matrix3x3 m(A,B,C); | 
| 360 |  |  | return(m.determinant()); | 
| 361 |  |  | } | 
| 362 |  |  |  | 
| 363 |  |  | //! \brief Calculate the Graph Potentials of a molecule | 
| 364 |  |  | //! | 
| 365 |  |  | //! based on | 
| 366 |  |  | //! V.E. and Rozenblit, A.B. Golender | 
| 367 |  |  | //! <em>Logical and Combinatorial Algorithms for Drug Design</em>. \n | 
| 368 |  |  | //! For an example see: | 
| 369 |  |  | //! Walters, W. P., Yalkowsky, S. H., \em JCICS, 1996, 36(5), 1015-1017. | 
| 370 |  |  | //! <a href="http://dx.doi.org/10.1021/ci950278o">DOI: 10.1021/ci950278o</a> | 
| 371 |  |  | void GraphPotentials(OBMol &mol, std::vector<double> &pot) | 
| 372 |  |  | { | 
| 373 |  |  | double det; | 
| 374 |  |  |  | 
| 375 |  |  | vector<vector<double> > g,c,h; | 
| 376 |  |  | construct_g_matrix(mol,g); | 
| 377 |  |  | invert_matrix(g,det); | 
| 378 |  |  | construct_c_matrix(mol,c); | 
| 379 |  |  | mult_matrix(h,g,c); | 
| 380 |  |  | pot.resize(mol.NumAtoms()+1); | 
| 381 |  |  |  | 
| 382 |  |  | for (unsigned int i = 0; i < mol.NumAtoms();i++) | 
| 383 |  |  | pot[i+1] = h[i][0]; | 
| 384 |  |  | } | 
| 385 |  |  |  | 
| 386 |  |  |  | 
| 387 |  |  | //! Construct the matrix G, which puts each atoms valence+1 | 
| 388 |  |  | //! on the diagonal and and -1 on the off diagonal if two | 
| 389 |  |  | //! atoms are connected. | 
| 390 |  |  | void construct_g_matrix(OBMol &mol, std::vector<std::vector<double> > &m) | 
| 391 |  |  | { | 
| 392 |  |  | unsigned int i,j; | 
| 393 |  |  |  | 
| 394 |  |  | OBAtom *atm1,*atm2; | 
| 395 |  |  | vector<OBNodeBase*>::iterator aint,bint; | 
| 396 |  |  |  | 
| 397 |  |  | m.resize(mol.NumAtoms()); | 
| 398 |  |  | for (i = 0; i < m.size(); i++) | 
| 399 |  |  | m[i].resize(mol.NumAtoms()); | 
| 400 |  |  |  | 
| 401 |  |  | for (atm1 = mol.BeginAtom(aint),i=0;atm1;atm1 = mol.NextAtom(aint),i++) | 
| 402 |  |  | for (atm2 = mol.BeginAtom(bint),j=0;atm2;atm2 = mol.NextAtom(bint),j++) | 
| 403 |  |  | { | 
| 404 |  |  | if (i == j) | 
| 405 |  |  | { | 
| 406 |  |  | m[i][j] = atm1->GetValence() + 1; | 
| 407 |  |  | m[i][j] += (double)atm1->GetAtomicNum()/10.0; | 
| 408 |  |  | m[i][j] += (double)atm1->GetHyb()/100.0; | 
| 409 |  |  | } | 
| 410 |  |  | else | 
| 411 |  |  | { | 
| 412 |  |  | if (atm1->IsConnected(atm2)) | 
| 413 |  |  | m[i][j] = -1; | 
| 414 |  |  | else | 
| 415 |  |  | m[i][j] = 0; | 
| 416 |  |  | } | 
| 417 |  |  | } | 
| 418 |  |  | } | 
| 419 |  |  |  | 
| 420 |  |  | //! Construct the matrix C, which is simply a column vector | 
| 421 |  |  | //! consisting of the valence for each atom | 
| 422 |  |  | void construct_c_matrix(OBMol &mol,std::vector<std::vector<double > > &m) | 
| 423 |  |  | { | 
| 424 |  |  | unsigned int i; | 
| 425 |  |  | OBAtom *atm1; | 
| 426 |  |  | vector<OBNodeBase*>::iterator aint; | 
| 427 |  |  |  | 
| 428 |  |  | m.resize(mol.NumAtoms()); | 
| 429 |  |  | for (i = 0; i < m.size(); i++) | 
| 430 |  |  | m[i].resize(1); | 
| 431 |  |  | for (atm1 = mol.BeginAtom(aint),i=0;atm1;atm1 = mol.NextAtom(aint),i++) | 
| 432 |  |  | { | 
| 433 |  |  | m[i][0] = atm1->GetValence(); | 
| 434 |  |  | } | 
| 435 |  |  | } | 
| 436 |  |  |  | 
| 437 |  |  | } // namespace OpenBabel | 
| 438 |  |  |  | 
| 439 |  |  | //! \file chiral.cpp | 
| 440 |  |  | //! \brief Detect chiral atoms and molecules. |