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 1979 by gezelter, Sat Apr 5 20:56:01 2014 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 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 180 | Line 184 | namespace OpenMD {
184      createConstraintPair(mol);
185  
186      //create non-bonded constraintPairs
187 <    for (int i = 0; i < molStamp->getNConstraints(); ++i) {
187 >    for (std::size_t i = 0; i < molStamp->getNConstraints(); ++i) {
188        ConstraintStamp* cStamp = molStamp->getConstraintStamp(i);
189        Atom* atomA;
190        Atom* atomB;
# Line 189 | Line 193 | namespace OpenMD {
193        atomB = mol->getAtomAt(cStamp->getB());
194        assert( atomA && atomB );
195        
192 –      RealType distance;
196        bool printConstraintForce = false;
197  
198 +      if (cStamp->havePrintConstraintForce()) {
199 +        printConstraintForce = cStamp->getPrintConstraintForce();
200 +      }
201 +
202        if (!cStamp->haveConstrainedDistance()) {
203          sprintf(painCave.errMsg,
204                  "Constraint Error: A non-bond constraint was specified\n"
# Line 199 | Line 206 | namespace OpenMD {
206          painCave.isFatal = 1;
207          simError();      
208        } else {
209 <        distance = cStamp->getConstrainedDistance();
210 <      }
211 <
212 <      if (cStamp->havePrintConstraintForce()) {
213 <        printConstraintForce = cStamp->getPrintConstraintForce();
214 <      }
215 <    
216 <      ConstraintElem* consElemA = new ConstraintElem(atomA);
210 <      ConstraintElem* consElemB = new ConstraintElem(atomB);
211 <      ConstraintPair* cPair = new ConstraintPair(consElemA, consElemB, distance,
212 <                                                 printConstraintForce);
213 <      mol->addConstraintPair(cPair);
209 >        RealType distance = cStamp->getConstrainedDistance();
210 >        ConstraintElem* consElemA = new ConstraintElem(atomA);
211 >        ConstraintElem* consElemB = new ConstraintElem(atomB);
212 >        ConstraintPair* cPair = new ConstraintPair(consElemA, consElemB,
213 >                                                   distance,
214 >                                                   printConstraintForce);
215 >        mol->addConstraintPair(cPair);
216 >      }  
217      }
218  
219      // now create the constraint elements:
# Line 224 | Line 227 | namespace OpenMD {
227        mol->setConstrainTotalCharge( molStamp->getConstrainTotalCharge() );
228      }
229  
230 <    //the construction of this molecule is finished
230 >    // The construction of this molecule is finished:
231      mol->complete();
232      
233      return mol;
# Line 304 | Line 307 | namespace OpenMD {
307    Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
308                                      BondStamp* stamp,
309                                      LocalIndexManager* localIndexMan) {
310 <    BondType* bondType;
310 >    BondTypeParser btParser;        
311 >    BondType* bondType = NULL;
312      Atom* atomA;
313      Atom* atomB;
314      
# Line 312 | Line 316 | namespace OpenMD {
316      atomB = mol->getAtomAt(stamp->getB());
317      
318      assert( atomA && atomB);
315 –    
316 –    bondType = ff->getBondType(atomA->getType(), atomB->getType());
319  
320 <    if (bondType == NULL) {
321 <      sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
322 <              atomA->getType().c_str(),
323 <              atomB->getType().c_str());
324 <      
325 <      painCave.isFatal = 1;
326 <      simError();
320 >    if (stamp->hasOverride()) {
321 >
322 >      try {
323 >        bondType = btParser.parseTypeAndPars(stamp->getOverrideType(),
324 >                                             stamp->getOverridePars() );
325 >      }
326 >      catch( OpenMDException e) {
327 >        sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
328 >                "for molecule %s\n",
329 >                e.what(), mol->getType().c_str() );
330 >        painCave.isFatal = 1;
331 >        simError();
332 >      }
333 >
334 >    } else {
335 >      bondType = ff->getBondType(atomA->getType(), atomB->getType());
336 >
337 >      if (bondType == NULL) {
338 >        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
339 >                atomA->getType().c_str(),
340 >                atomB->getType().c_str());
341 >        
342 >        painCave.isFatal = 1;
343 >        simError();
344 >      }
345      }
346 +    
347      Bond* bond = new Bond(atomA, atomB, bondType);
348  
349      //set the local index of this bond, the global index will be set later
# Line 342 | Line 363 | namespace OpenMD {
363    Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol,
364                                      BendStamp* stamp,
365                                      LocalIndexManager* localIndexMan) {
366 <    Bend* bend = NULL;
367 <    std::vector<int> bendAtoms = stamp->getMembers();
366 >    BendTypeParser btParser;
367 >    BendType* bendType = NULL;
368 >    Bend* bend = NULL;
369 >    
370 >    std::vector<int> bendAtoms = stamp->getMembers();
371      if (bendAtoms.size() == 3) {
372        Atom* atomA = mol->getAtomAt(bendAtoms[0]);
373        Atom* atomB = mol->getAtomAt(bendAtoms[1]);
374        Atom* atomC = mol->getAtomAt(bendAtoms[2]);
375        
376 <      assert( atomA && atomB && atomC);
377 <      
378 <      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());
376 >      assert( atomA && atomB && atomC );
377 >
378 >      if (stamp->hasOverride()) {
379          
380 <        painCave.isFatal = 1;
381 <        simError();
380 >        try {
381 >          bendType = btParser.parseTypeAndPars(stamp->getOverrideType(),
382 >                                               stamp->getOverridePars() );
383 >        }
384 >        catch( OpenMDException e) {
385 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
386 >                  "for molecule %s\n",
387 >                  e.what(), mol->getType().c_str() );
388 >          painCave.isFatal = 1;
389 >          simError();
390 >        }
391 >      } else {
392 >        
393 >        bendType = ff->getBendType(atomA->getType().c_str(),
394 >                                             atomB->getType().c_str(),
395 >                                             atomC->getType().c_str());
396 >      
397 >        if (bendType == NULL) {
398 >          sprintf(painCave.errMsg,
399 >                  "Can not find Matching Bend Type for[%s, %s, %s]",
400 >                  atomA->getType().c_str(),
401 >                  atomB->getType().c_str(),
402 >                  atomC->getType().c_str());
403 >          
404 >          painCave.isFatal = 1;
405 >          simError();
406 >        }
407        }
408        
409        bend = new Bend(atomA, atomB, atomC, bendType);
410 +      
411      } else if ( bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
412        int ghostIndex = stamp->getGhostVectorSource();
413 <      int normalIndex = ghostIndex != bendAtoms[0] ? bendAtoms[0] : bendAtoms[1];
413 >      int normalIndex = ghostIndex != bendAtoms[0] ?
414 >        bendAtoms[0] : bendAtoms[1];
415        Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
416        DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
417        if (ghostAtom == NULL) {
# Line 377 | Line 419 | namespace OpenMD {
419          painCave.isFatal = 1;
420          simError();
421        }
380 –                
381 –      BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
422  
423 <      if (bendType == NULL) {
424 <        sprintf(painCave.errMsg,
425 <                "Can not find Matching Bend Type for[%s, %s, %s]",
426 <                normalAtom->getType().c_str(),
427 <                ghostAtom->getType().c_str(),
428 <                "GHOST");
429 <
430 <        painCave.isFatal = 1;
431 <        simError();
423 >      if (stamp->hasOverride()) {
424 >        
425 >        try {
426 >          bendType = btParser.parseTypeAndPars(stamp->getOverrideType(),
427 >                                               stamp->getOverridePars() );
428 >        }
429 >        catch( OpenMDException e) {
430 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
431 >                  "for molecule %s\n",
432 >                  e.what(), mol->getType().c_str() );
433 >          painCave.isFatal = 1;
434 >          simError();
435 >        }
436 >      } else {
437 >      
438 >        bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(),
439 >                                   "GHOST");
440 >        
441 >        if (bendType == NULL) {
442 >          sprintf(painCave.errMsg,
443 >                  "Can not find Matching Bend Type for[%s, %s, %s]",
444 >                  normalAtom->getType().c_str(),
445 >                  ghostAtom->getType().c_str(),
446 >                  "GHOST");
447 >          
448 >          painCave.isFatal = 1;
449 >          simError();
450 >        }
451        }
452        
453        bend = new GhostBend(normalAtom, ghostAtom, bendType);      
454        
455      }
456 <
456 >    
457      //set the local index of this bend, the global index will be set later
458      bend->setLocalIndex(localIndexMan->getNextBendIndex());
459 <
459 >    
460      // The rule for naming a bend is: MoleculeName_Bend_Integer
461      // The first part is the name of the molecule
462      // The second part is always fixed as "Bend"
463      // The third part is the index of the bend defined in meta-data file
464      // For example, Butane_Bend_0 is a valid Bend name in a butane molecule
465 <
465 >    
466      std::string s = OpenMD_itoa(mol->getNBends(), 10);
467      bend->setName(mol->getType() + "_Bend_" + s.c_str());    
468      return bend;
# Line 413 | Line 472 | namespace OpenMD {
472                                            TorsionStamp* stamp,
473                                            LocalIndexManager* localIndexMan) {
474  
475 +    TorsionTypeParser ttParser;
476 +    TorsionType* torsionType = NULL;
477      Torsion* torsion = NULL;
478 +
479      std::vector<int> torsionAtoms = stamp->getMembers();
480      if (torsionAtoms.size() < 3) {
481          return torsion;
# Line 426 | Line 488 | namespace OpenMD {
488      if (torsionAtoms.size() == 4) {
489        Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
490  
491 <      assert(atomA && atomB && atomC && atomD);
491 >      assert(atomA && atomB && atomC && atomD );
492 >
493 >      if (stamp->hasOverride()) {
494          
495 <      TorsionType* torsionType = ff->getTorsionType(atomA->getType(),
496 <                                                    atomB->getType(),
497 <                                                    atomC->getType(),
498 <                                                    atomD->getType());
499 <      if (torsionType == NULL) {
500 <        sprintf(painCave.errMsg,
501 <                "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
502 <                atomA->getType().c_str(),
503 <                atomB->getType().c_str(),
504 <                atomC->getType().c_str(),
505 <                atomD->getType().c_str());
495 >        try {
496 >          torsionType = ttParser.parseTypeAndPars(stamp->getOverrideType(),
497 >                                                  stamp->getOverridePars() );
498 >        }
499 >        catch( OpenMDException e) {
500 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
501 >                  "for molecule %s\n",
502 >                  e.what(), mol->getType().c_str() );
503 >          painCave.isFatal = 1;
504 >          simError();
505 >        }
506 >      } else {
507 >
508          
509 <        painCave.isFatal = 1;
510 <        simError();
509 >        torsionType = ff->getTorsionType(atomA->getType(),
510 >                                         atomB->getType(),
511 >                                         atomC->getType(),
512 >                                         atomD->getType());
513 >        if (torsionType == NULL) {
514 >          sprintf(painCave.errMsg,
515 >                  "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
516 >                  atomA->getType().c_str(),
517 >                  atomB->getType().c_str(),
518 >                  atomC->getType().c_str(),
519 >                  atomD->getType().c_str());
520 >          
521 >          painCave.isFatal = 1;
522 >          simError();
523 >        }
524        }
525        
526        torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
527 <    }
449 <    else {
527 >    } else {
528        
529        DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(stamp->getGhostVectorSource()));
530        if (dAtom == NULL) {
# Line 454 | Line 532 | namespace OpenMD {
532          painCave.isFatal = 1;
533          simError();
534        }        
535 <      
536 <      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");
535 >
536 >      if (stamp->hasOverride()) {
537          
538 <        painCave.isFatal = 1;
539 <        simError();
540 <      }
538 >        try {
539 >          torsionType = ttParser.parseTypeAndPars(stamp->getOverrideType(),
540 >                                                  stamp->getOverridePars() );
541 >        }
542 >        catch( OpenMDException e) {
543 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
544 >                  "for molecule %s\n",
545 >                  e.what(), mol->getType().c_str() );
546 >          painCave.isFatal = 1;
547 >          simError();
548 >        }
549 >      } else {              
550 >        torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
551 >                                         atomC->getType(), "GHOST");
552        
553 +        if (torsionType == NULL) {
554 +          sprintf(painCave.errMsg,
555 +                  "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
556 +                  atomA->getType().c_str(),
557 +                  atomB->getType().c_str(),
558 +                  atomC->getType().c_str(),
559 +                  "GHOST");
560 +          
561 +          painCave.isFatal = 1;
562 +          simError();
563 +        }
564 +      }
565 +
566        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
567      }
568  
# Line 490 | Line 584 | namespace OpenMD {
584    Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol,
585                                                InversionStamp* stamp,
586                                                LocalIndexManager* localIndexMan) {
587 <    
587 >
588 >    InversionTypeParser itParser;
589 >    InversionType* inversionType = NULL;
590      Inversion* inversion = NULL;
591 +    
592      int center = stamp->getCenter();
593      std::vector<int> satellites = stamp->getSatellites();
594      if (satellites.size() != 3) {
# Line 504 | Line 601 | namespace OpenMD {
601      Atom* atomD = mol->getAtomAt(satellites[2]);
602        
603      assert(atomA && atomB && atomC && atomD);
507 –    
508 –    InversionType* inversionType = ff->getInversionType(atomA->getType(),
509 –                                                        atomB->getType(),
510 –                                                        atomC->getType(),
511 –                                                        atomD->getType());
604  
605 <    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());
605 >    if (stamp->hasOverride()) {
606        
607 <      painCave.isFatal = 0;
608 <      painCave.severity = OPENMD_INFO;
609 <      simError();
610 <      return NULL;
607 >      try {
608 >        inversionType = itParser.parseTypeAndPars(stamp->getOverrideType(),
609 >                                                  stamp->getOverridePars() );
610 >      }
611 >      catch( OpenMDException e) {
612 >        sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
613 >                "for molecule %s\n",
614 >                e.what(), mol->getType().c_str() );
615 >        painCave.isFatal = 1;
616 >        simError();
617 >      }
618      } else {
619        
620 +      inversionType = ff->getInversionType(atomA->getType(),
621 +                                           atomB->getType(),
622 +                                           atomC->getType(),
623 +                                           atomD->getType());
624 +      
625 +      if (inversionType == NULL) {
626 +        sprintf(painCave.errMsg,
627 +                "No Matching Inversion Type for[%s, %s, %s, %s]\n"
628 +                "\t(May not be a problem: not all inversions are parametrized)\n",
629 +                atomA->getType().c_str(),
630 +                atomB->getType().c_str(),
631 +                atomC->getType().c_str(),
632 +                atomD->getType().c_str());
633 +        
634 +        painCave.isFatal = 0;
635 +        painCave.severity = OPENMD_INFO;
636 +        simError();
637 +      }
638 +    }
639 +    if (inversionType != NULL) {
640 +      
641        inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType);
642 <
642 >      
643        // set the local index of this inversion, the global index will
644        // be set later
645        inversion->setLocalIndex(localIndexMan->getNextInversionIndex());
646 <
646 >      
647        // The rule for naming an inversion is: MoleculeName_Inversion_Integer
648        // The first part is the name of the molecule
649        // The second part is always fixed as "Inversion"
# Line 540 | Line 654 | namespace OpenMD {
654        std::string s = OpenMD_itoa(mol->getNInversions(), 10);
655        inversion->setName(mol->getType() + "_Inversion_" + s.c_str());
656        return inversion;
657 +    } else {
658 +      return NULL;
659      }
660    }
661    

Diff Legend

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