ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/applications/atom2md/oopseformat.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/applications/atom2md/oopseformat.cpp (file contents):
Revision 3429 by gezelter, Thu Jul 3 12:45:30 2008 UTC vs.
Revision 3446 by cli2, Wed Sep 10 19:51:45 2008 UTC

# Line 19 | Line 19 | GNU General Public License for more details.
19   #include <openbabel/obiter.h>
20   #include <openbabel/mol.h>
21   #include <openbabel/chains.h>
22 + #include <openbabel/data.h>
23   #include <fstream>
24  
25   #include "utils/StringUtils.hpp"
# Line 60 | Line 61 | namespace OpenBabel
61      OBMol* createMolFromFragment(OBMol& mol, vector<int>& fragment);
62      void WriteMDFile(vector<OBMol*> mols, vector<int> numMols, ostream& os,
63                       OBMol& mol, vector<int>& indices);
64 +    void CalcBoundingBox(OBMol &mol,
65 +                         double &min_x, double &max_x,
66 +                         double &min_y, double &max_y,
67 +                         double &min_z, double &max_z);
68 +      
69    };
70    
71    //Make an instance of the format class
# Line 187 | Line 193 | namespace OpenBabel
193      unsigned int i;
194      const int BUFFLEN = 1024;
195      char buffer[BUFFLEN];
196 <    string str, str1;
196 >    string str, str1, str2, str3;
197      OBAtom *a, *b, *c, *d;    
198      bool molIsWater = false;
199      OBResidue *r;
200 <    int resKey;
200 >    int resKey, myserial;
201      char type_name[10];
202      char *element_name;
203      int res_num;
204 <    OBChainsParser* chainParser = new OBChainsParser();
205 <
204 >    OBChainsParser* chainParser = new OBChainsParser();  
205 >    double min_x, max_x, min_y, max_y, min_z, max_z; /* Edges of bounding box */
206      
207      os << "<OOPSE version=4>" << endl;
208      os << "  <MetaData>" << endl << endl;
# Line 247 | Line 253 | namespace OpenBabel
253              
254              // If we know what residue we've got, the specific atom name can
255              // be used to help specify partial charges.
256 +
257 +            //resdat.SetResName(resName);
258              
259 +            // atom type from residue:
260              str = r->GetAtomID(&*atom);
261 +          
262 +            // arginine has separate indices for chemically-identical
263 +            // nitrogen atoms:
264 +            if (resName.compare("ARG") == 0) {
265 +              if (str.compare("NH1") == 0 || str.compare("NH2") == 0) {
266 +                str = "NH";
267 +              }
268 +            }
269 +            if (resName.compare("VAL") == 0) {
270 +              if (str.compare("CG1") == 0 || str.compare("CG2") == 0) {
271 +                str = "CG";
272 +              }
273 +            }
274 +            if (resName.compare("LEU") == 0) {
275 +              if (str.compare("CD1") == 0 || str.compare("CD2") == 0) {
276 +                str = "CD";
277 +              }
278 +            }
279 +            if (resName.compare("ASP") == 0) {
280 +              if (str.compare("OD1") == 0 || str.compare("OD2") == 0) {
281 +                str = "OD";
282 +              }
283 +            }
284 +            if (resName.compare("GLU") == 0) {
285 +              if (str.compare("OE1") == 0 || str.compare("OE2") == 0) {
286 +                str = "OE";
287 +              }
288 +            }
289 +            if (resName.compare("TYR") == 0) {
290 +              if (str.compare("CD1") == 0 || str.compare("CD2") == 0) {
291 +                str = "CD";
292 +              }
293 +              if (str.compare("CE1") == 0 || str.compare("CE2") == 0) {
294 +                str = "CE";
295 +              }
296 +            }
297 +            
298 +
299 +            if ((&*atom)->IsHydrogen()) {
300 +               FOR_NBORS_OF_ATOM(nbr, *atom) {
301 +                 str2 = r->GetAtomID(&*nbr);
302 +                 size_t startpos = str2.find_first_not_of(" ");
303 +                 size_t endpos = str2.find_last_not_of(" ");
304 +                 if ((endpos - startpos) < 1) {
305 +                   // if the bonded atom type has only one character (i.e. N)
306 +                   // then the hydrogen will be labeled "HN" to show what
307 +                   // kind of proton it is:
308 +                   str3 = str2;
309 +                 } else {
310 +                   if (str2.compare("OH") == 0) {
311 +                      str3 = "O";
312 +                   } else {
313 +                     // When the bonded atom type is more specific, we drop
314 +                     // the first character:  i.e. H bonded to OG1 is HG1 type:
315 +                     str3 = str2.substr(startpos+1, endpos-startpos);
316 +                   }
317 +                 }
318 +                str = "H" + str3;
319 +               }
320 +               // same problem with arginine NH atoms, but now for connected hydrogens
321 +               if (resName.compare("ARG") == 0) {
322 +                 if (str.compare("HH1") == 0 || str.compare("HH2") == 0) {
323 +                   str = "HH";
324 +                 }
325 +               }
326 +               if (resName.compare("VAL") == 0) {
327 +                 if (str.compare("HG1") == 0 || str.compare("HG2") == 0) {
328 +                   str = "HG";
329 +                 }
330 +               }
331 +               if (resName.compare("LEU") == 0) {
332 +                 if (str.compare("HD1") == 0 || str.compare("HD2") == 0) {
333 +                   str = "HD";
334 +                 }
335 +               }
336 +               if (resName.compare("TYR") == 0) {
337 +                 if (str.compare("HD1") == 0 || str.compare("HD2") == 0) {
338 +                   str = "HD";
339 +                 }
340 +                 if (str.compare("HE1") == 0 || str.compare("HE2") == 0) {
341 +                   str = "HE";
342 +                 }
343 +               }
344 +
345 +            }
346 +
347 +            // atom type from residue table:
348 +            //resdat.LookupType(str, str2, hyb);
349              size_t startpos = str.find_first_not_of(" ");
350              size_t endpos = str.find_last_not_of(" ");
351              str = str.substr( startpos, endpos-startpos+1 );
# Line 275 | Line 372 | namespace OpenBabel
372            else
373              os << "members(" << b2 <<  ", " << b1 << "); ";
374            
375 <          os << "}\n";
375 >          os << "}" << endl;
376          }
377          
378 <        os << "\n";
379 <
283 <        std::vector<int> possibleInversion;
284 <        FOR_ATOMS_OF_MOL(atom, *pmol) {
285 <          possibleInversion.clear();
286 <          
287 <          FOR_NBORS_OF_ATOM(nbor, &*atom) {
288 <            possibleInversion.push_back(atomMap[&(*nbor)]);          
289 <          }
290 <          
291 <          if (possibleInversion.size() == 3) {
292 <
293 <            os << "  inversion { ";
294 <            os << "center(" << atomMap[&(*atom)] << "); ";
295 <            os << "}\n";
296 <          }
297 <          
298 <        }
378 >        os << endl;
379 >      
380          os << "}" << endl;
381          os << endl;
382        }
383      }
384      
385      os << endl;
386 <    
306 <    
386 >        
387      for(i=0; i < mols.size(); ++i) {
388        OBMol* pmol = mols[i];      
389        os << "component{" << endl;
# Line 324 | Line 404 | namespace OpenBabel
404      sprintf(buffer, "        Time: %.10g", 0.0);
405      
406      os << buffer << endl;
327 –    
407  
408 <    // should really compute a bounding box here:
409 <    sprintf(buffer, "        Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}", 100.0, 0.0, 0.0, 0.0, 100.0, 0.0, 0.0, 0.0, 100.0);
408 >    CalcBoundingBox(mol, min_x, max_x, min_y, max_y, min_z, max_z);
409 >
410 >    // still to do: should compute a bounding box here
411 >    sprintf(buffer, "        Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}",
412 >            max_x - min_x, 0.0, 0.0, 0.0, max_y - min_y, 0.0, 0.0, 0.0, max_z - min_z);
413      
414      os << buffer << endl;
415      os << "    </FrameData>" << endl;
# Line 339 | Line 421 | namespace OpenBabel
421  
422      for(vector<int>::iterator i = indices.begin();i != indices.end(); ++i) {    
423        atom = mol.GetAtom(*i);
424 <      sprintf(buffer, "%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e", *i - 1, "pv", atom->GetX(), atom->GetY(), atom->GetZ(), 0.0, 0.0, 0.0);
424 >      sprintf(buffer, "%10d %7s %18.10g %18.10g %18.10g %13e %13e %13e", *i - 1,
425 >              "pv", atom->GetX(), atom->GetY(), atom->GetZ(), 0.0, 0.0, 0.0);
426        os << buffer << endl;
427      }
428      os << "    </StuntDoubles>" << endl;
429      os << "  </Snapshot>" << endl;
430      os << "</OOPSE>" << endl;
431    }
432 <  
432 >
433 >  void OOPSEFormat::CalcBoundingBox(OBMol &mol,
434 >                                    double &min_x, double &max_x,
435 >                                    double &min_y, double &max_y,
436 >                                    double &min_z, double &max_z
437 >                                    )
438 >  {
439 >    /* ---- Init bounding-box variables ---- */
440 >    min_x = (double) 0.0;
441 >    max_x = (double) 0.0;
442 >    min_y = (double) 0.0;
443 >    max_y = (double) 0.0;
444 >    min_z = (double) 0.0;
445 >    max_z = (double) 0.0;
446 >    
447 >    /* ---- Check all atoms ---- */
448 >    for(unsigned int i = 1; i <= mol.NumAtoms(); ++i)
449 >      {
450 >        
451 >        /* ---- Get a pointer to ith atom ---- */
452 >        OBAtom *atom = mol.GetAtom(i);
453 >        
454 >        /* ---- Check for minimal/maximal x-position ---- */
455 >        if (atom -> GetX() < min_x)
456 >          min_x = atom -> GetX();
457 >        if (atom -> GetX() > max_x)
458 >          max_x = atom -> GetX();
459 >        
460 >        /* ---- Check for minimal/maximal y-position ---- */
461 >        if (atom -> GetY() < min_y)
462 >          min_y = atom -> GetY();
463 >        if (atom -> GetY() > max_y)
464 >          max_y = atom -> GetY();
465 >        
466 >        /* ---- Check for minimal/maximal z-position ---- */
467 >        if (atom -> GetZ() < min_z)
468 >          min_z = atom -> GetZ();
469 >        if (atom -> GetZ() > max_z)
470 >          max_z = atom -> GetZ();
471 >        
472 >      }    
473 >  }
474   } //namespace OpenBabel
475  

Diff Legend

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