ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/brains/MoleculeCreator.cpp
(Generate patch)

Comparing trunk/src/brains/MoleculeCreator.cpp (file contents):
Revision 2045 by gezelter, Sat Apr 5 20:56:01 2014 UTC vs.
Revision 2046 by gezelter, Tue Dec 2 22:11:04 2014 UTC

# Line 56 | Line 56
56   #include "primitives/GhostTorsion.hpp"
57   #include "types/AtomType.hpp"
58   #include "types/FixedBondType.hpp"
59 + #include "types/BondTypeParser.hpp"
60 + #include "types/BendTypeParser.hpp"
61 + #include "types/TorsionTypeParser.hpp"
62 + #include "types/InversionTypeParser.hpp"
63   #include "utils/simError.h"
64   #include "utils/StringUtils.hpp"
65  
# Line 96 | Line 100 | namespace OpenMD {
100      int nBonds = molStamp->getNBonds();
101  
102      for (int i = 0; i < nBonds; ++i) {
103 <      currentBondStamp = molStamp->getBondStamp(i);
103 >      currentBondStamp = molStamp->getBondStamp(i);        
104        bond = createBond(ff, mol, currentBondStamp, localIndexMan);
105        mol->addBond(bond);
106      }
# Line 304 | Line 308 | namespace OpenMD {
308    Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
309                                      BondStamp* stamp,
310                                      LocalIndexManager* localIndexMan) {
311 <    BondType* bondType;
311 >    BondTypeParser btParser;        
312 >    BondType* bondType = NULL;
313      Atom* atomA;
314      Atom* atomB;
315      
# Line 312 | Line 317 | namespace OpenMD {
317      atomB = mol->getAtomAt(stamp->getB());
318      
319      assert( atomA && atomB);
320 <    
321 <    bondType = ff->getBondType(atomA->getType(), atomB->getType());
320 >
321 >    if (stamp->hasOverride()) {
322  
323 <    if (bondType == NULL) {
324 <      sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
325 <              atomA->getType().c_str(),
326 <              atomB->getType().c_str());
327 <      
328 <      painCave.isFatal = 1;
329 <      simError();
323 >      try {
324 >        bondType = btParser.parseTypeAndPars(stamp->getOverrideType(),
325 >                                             stamp->getOverridePars() );
326 >      }
327 >      catch( OpenMDException e) {
328 >        sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
329 >                "for molecule %s\n",
330 >                e.what(), mol->getType().c_str() );
331 >        painCave.isFatal = 1;
332 >        simError();
333 >      }
334 >
335 >    } else {
336 >      bondType = ff->getBondType(atomA->getType(), atomB->getType());
337 >
338 >      if (bondType == NULL) {
339 >        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
340 >                atomA->getType().c_str(),
341 >                atomB->getType().c_str());
342 >        
343 >        painCave.isFatal = 1;
344 >        simError();
345 >      }
346      }
347 +    
348      Bond* bond = new Bond(atomA, atomB, bondType);
349  
350      //set the local index of this bond, the global index will be set later
# Line 342 | Line 364 | namespace OpenMD {
364    Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol,
365                                      BendStamp* stamp,
366                                      LocalIndexManager* localIndexMan) {
367 <    Bend* bend = NULL;
368 <    std::vector<int> bendAtoms = stamp->getMembers();
367 >    BendTypeParser btParser;
368 >    BendType* bendType = NULL;
369 >    Bend* bend = NULL;
370 >    
371 >    std::vector<int> bendAtoms = stamp->getMembers();
372      if (bendAtoms.size() == 3) {
373        Atom* atomA = mol->getAtomAt(bendAtoms[0]);
374        Atom* atomB = mol->getAtomAt(bendAtoms[1]);
375        Atom* atomC = mol->getAtomAt(bendAtoms[2]);
376        
377 <      assert( atomA && atomB && atomC);
378 <      
379 <      BendType* bendType = ff->getBendType(atomA->getType().c_str(),
355 <                                           atomB->getType().c_str(),
356 <                                           atomC->getType().c_str());
357 <      
358 <      if (bendType == NULL) {
359 <        sprintf(painCave.errMsg,
360 <                "Can not find Matching Bend Type for[%s, %s, %s]",
361 <                atomA->getType().c_str(),
362 <                atomB->getType().c_str(),
363 <                atomC->getType().c_str());
377 >      assert( atomA && atomB && atomC );
378 >
379 >      if (stamp->hasOverride()) {
380          
381 <        painCave.isFatal = 1;
382 <        simError();
381 >        try {
382 >          bendType = btParser.parseTypeAndPars(stamp->getOverrideType(),
383 >                                               stamp->getOverridePars() );
384 >        }
385 >        catch( OpenMDException e) {
386 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
387 >                  "for molecule %s\n",
388 >                  e.what(), mol->getType().c_str() );
389 >          painCave.isFatal = 1;
390 >          simError();
391 >        }
392 >      } else {
393 >        
394 >        bendType = ff->getBendType(atomA->getType().c_str(),
395 >                                             atomB->getType().c_str(),
396 >                                             atomC->getType().c_str());
397 >      
398 >        if (bendType == NULL) {
399 >          sprintf(painCave.errMsg,
400 >                  "Can not find Matching Bend Type for[%s, %s, %s]",
401 >                  atomA->getType().c_str(),
402 >                  atomB->getType().c_str(),
403 >                  atomC->getType().c_str());
404 >          
405 >          painCave.isFatal = 1;
406 >          simError();
407 >        }
408        }
409        
410        bend = new Bend(atomA, atomB, atomC, bendType);
411 +      
412      } else if ( bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
413        int ghostIndex = stamp->getGhostVectorSource();
414 <      int normalIndex = ghostIndex != bendAtoms[0] ? bendAtoms[0] : bendAtoms[1];
414 >      int normalIndex = ghostIndex != bendAtoms[0] ?
415 >        bendAtoms[0] : bendAtoms[1];
416        Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
417        DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
418        if (ghostAtom == NULL) {
# Line 377 | Line 420 | namespace OpenMD {
420          painCave.isFatal = 1;
421          simError();
422        }
380 –                
381 –      BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
423  
424 <      if (bendType == NULL) {
425 <        sprintf(painCave.errMsg,
426 <                "Can not find Matching Bend Type for[%s, %s, %s]",
427 <                normalAtom->getType().c_str(),
428 <                ghostAtom->getType().c_str(),
429 <                "GHOST");
430 <
431 <        painCave.isFatal = 1;
432 <        simError();
424 >      if (stamp->hasOverride()) {
425 >        
426 >        try {
427 >          bendType = btParser.parseTypeAndPars(stamp->getOverrideType(),
428 >                                               stamp->getOverridePars() );
429 >        }
430 >        catch( OpenMDException e) {
431 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
432 >                  "for molecule %s\n",
433 >                  e.what(), mol->getType().c_str() );
434 >          painCave.isFatal = 1;
435 >          simError();
436 >        }
437 >      } else {
438 >      
439 >        bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(),
440 >                                   "GHOST");
441 >        
442 >        if (bendType == NULL) {
443 >          sprintf(painCave.errMsg,
444 >                  "Can not find Matching Bend Type for[%s, %s, %s]",
445 >                  normalAtom->getType().c_str(),
446 >                  ghostAtom->getType().c_str(),
447 >                  "GHOST");
448 >          
449 >          painCave.isFatal = 1;
450 >          simError();
451 >        }
452        }
453        
454        bend = new GhostBend(normalAtom, ghostAtom, bendType);      
455        
456      }
457 <
457 >    
458      //set the local index of this bend, the global index will be set later
459      bend->setLocalIndex(localIndexMan->getNextBendIndex());
460 <
460 >    
461      // The rule for naming a bend is: MoleculeName_Bend_Integer
462      // The first part is the name of the molecule
463      // The second part is always fixed as "Bend"
464      // The third part is the index of the bend defined in meta-data file
465      // For example, Butane_Bend_0 is a valid Bend name in a butane molecule
466 <
466 >    
467      std::string s = OpenMD_itoa(mol->getNBends(), 10);
468      bend->setName(mol->getType() + "_Bend_" + s.c_str());    
469      return bend;
# Line 413 | Line 473 | namespace OpenMD {
473                                            TorsionStamp* stamp,
474                                            LocalIndexManager* localIndexMan) {
475  
476 +    TorsionTypeParser ttParser;
477 +    TorsionType* torsionType = NULL;
478      Torsion* torsion = NULL;
479 +
480      std::vector<int> torsionAtoms = stamp->getMembers();
481      if (torsionAtoms.size() < 3) {
482          return torsion;
# Line 426 | Line 489 | namespace OpenMD {
489      if (torsionAtoms.size() == 4) {
490        Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
491  
492 <      assert(atomA && atomB && atomC && atomD);
492 >      assert(atomA && atomB && atomC && atomD );
493 >
494 >      if (stamp->hasOverride()) {
495          
496 <      TorsionType* torsionType = ff->getTorsionType(atomA->getType(),
497 <                                                    atomB->getType(),
498 <                                                    atomC->getType(),
499 <                                                    atomD->getType());
500 <      if (torsionType == NULL) {
501 <        sprintf(painCave.errMsg,
502 <                "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
503 <                atomA->getType().c_str(),
504 <                atomB->getType().c_str(),
505 <                atomC->getType().c_str(),
506 <                atomD->getType().c_str());
496 >        try {
497 >          torsionType = ttParser.parseTypeAndPars(stamp->getOverrideType(),
498 >                                                  stamp->getOverridePars() );
499 >        }
500 >        catch( OpenMDException e) {
501 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
502 >                  "for molecule %s\n",
503 >                  e.what(), mol->getType().c_str() );
504 >          painCave.isFatal = 1;
505 >          simError();
506 >        }
507 >      } else {
508 >
509          
510 <        painCave.isFatal = 1;
511 <        simError();
510 >        torsionType = ff->getTorsionType(atomA->getType(),
511 >                                         atomB->getType(),
512 >                                         atomC->getType(),
513 >                                         atomD->getType());
514 >        if (torsionType == NULL) {
515 >          sprintf(painCave.errMsg,
516 >                  "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
517 >                  atomA->getType().c_str(),
518 >                  atomB->getType().c_str(),
519 >                  atomC->getType().c_str(),
520 >                  atomD->getType().c_str());
521 >          
522 >          painCave.isFatal = 1;
523 >          simError();
524 >        }
525        }
526        
527        torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
528 <    }
449 <    else {
528 >    } else {
529        
530        DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(stamp->getGhostVectorSource()));
531        if (dAtom == NULL) {
# Line 454 | Line 533 | namespace OpenMD {
533          painCave.isFatal = 1;
534          simError();
535        }        
536 <      
537 <      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
459 <                                                    atomC->getType(), "GHOST");
460 <      
461 <      if (torsionType == NULL) {
462 <        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
463 <                atomA->getType().c_str(),
464 <                atomB->getType().c_str(),
465 <                atomC->getType().c_str(),
466 <                "GHOST");
536 >
537 >      if (stamp->hasOverride()) {
538          
539 <        painCave.isFatal = 1;
540 <        simError();
541 <      }
539 >        try {
540 >          torsionType = ttParser.parseTypeAndPars(stamp->getOverrideType(),
541 >                                                  stamp->getOverridePars() );
542 >        }
543 >        catch( OpenMDException e) {
544 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
545 >                  "for molecule %s\n",
546 >                  e.what(), mol->getType().c_str() );
547 >          painCave.isFatal = 1;
548 >          simError();
549 >        }
550 >      } else {              
551 >        torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
552 >                                         atomC->getType(), "GHOST");
553        
554 +        if (torsionType == NULL) {
555 +          sprintf(painCave.errMsg,
556 +                  "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
557 +                  atomA->getType().c_str(),
558 +                  atomB->getType().c_str(),
559 +                  atomC->getType().c_str(),
560 +                  "GHOST");
561 +          
562 +          painCave.isFatal = 1;
563 +          simError();
564 +        }
565 +      }
566 +
567        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
568      }
569  
# Line 490 | Line 585 | namespace OpenMD {
585    Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol,
586                                                InversionStamp* stamp,
587                                                LocalIndexManager* localIndexMan) {
588 <    
588 >
589 >    InversionTypeParser itParser;
590 >    InversionType* inversionType = NULL;
591      Inversion* inversion = NULL;
592 +    
593      int center = stamp->getCenter();
594      std::vector<int> satellites = stamp->getSatellites();
595      if (satellites.size() != 3) {
# Line 504 | Line 602 | namespace OpenMD {
602      Atom* atomD = mol->getAtomAt(satellites[2]);
603        
604      assert(atomA && atomB && atomC && atomD);
507 –    
508 –    InversionType* inversionType = ff->getInversionType(atomA->getType(),
509 –                                                        atomB->getType(),
510 –                                                        atomC->getType(),
511 –                                                        atomD->getType());
605  
606 <    if (inversionType == NULL) {
514 <      sprintf(painCave.errMsg, "No Matching Inversion Type for[%s, %s, %s, %s]\n"
515 <              "\t(May not be a problem: not all inversions are parametrized)\n",
516 <              atomA->getType().c_str(),
517 <              atomB->getType().c_str(),
518 <              atomC->getType().c_str(),
519 <              atomD->getType().c_str());
606 >    if (stamp->hasOverride()) {
607        
608 <      painCave.isFatal = 0;
609 <      painCave.severity = OPENMD_INFO;
610 <      simError();
611 <      return NULL;
608 >      try {
609 >        inversionType = itParser.parseTypeAndPars(stamp->getOverrideType(),
610 >                                                  stamp->getOverridePars() );
611 >      }
612 >      catch( OpenMDException e) {
613 >        sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
614 >                "for molecule %s\n",
615 >                e.what(), mol->getType().c_str() );
616 >        painCave.isFatal = 1;
617 >        simError();
618 >      }
619      } else {
620        
621 +      inversionType = ff->getInversionType(atomA->getType(),
622 +                                           atomB->getType(),
623 +                                           atomC->getType(),
624 +                                           atomD->getType());
625 +      
626 +      if (inversionType == NULL) {
627 +        sprintf(painCave.errMsg,
628 +                "No Matching Inversion Type for[%s, %s, %s, %s]\n"
629 +                "\t(May not be a problem: not all inversions are parametrized)\n",
630 +                atomA->getType().c_str(),
631 +                atomB->getType().c_str(),
632 +                atomC->getType().c_str(),
633 +                atomD->getType().c_str());
634 +        
635 +        painCave.isFatal = 0;
636 +        painCave.severity = OPENMD_INFO;
637 +        simError();
638 +      }
639 +    }
640 +    if (inversionType != NULL) {
641 +      
642        inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType);
643 <
643 >      
644        // set the local index of this inversion, the global index will
645        // be set later
646        inversion->setLocalIndex(localIndexMan->getNextInversionIndex());
647 <
647 >      
648        // The rule for naming an inversion is: MoleculeName_Inversion_Integer
649        // The first part is the name of the molecule
650        // The second part is always fixed as "Inversion"
# Line 540 | Line 655 | namespace OpenMD {
655        std::string s = OpenMD_itoa(mol->getNInversions(), 10);
656        inversion->setName(mol->getType() + "_Inversion_" + s.c_str());
657        return inversion;
658 +    } else {
659 +      return NULL;
660      }
661    }
662    

Diff Legend

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