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 1953 by gezelter, Thu Dec 5 18:19:26 2013 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 175 | Line 179 | namespace OpenMD {
179        cutoffGroup = createCutoffGroup(mol, *fai, localIndexMan);
180        mol->addCutoffGroup(cutoffGroup);
181      }
182 <    //create constraints
182 >
183 >    //create bonded constraintPairs:
184      createConstraintPair(mol);
185 +
186 +    //create non-bonded constraintPairs
187 +    for (std::size_t i = 0; i < molStamp->getNConstraints(); ++i) {
188 +      ConstraintStamp* cStamp = molStamp->getConstraintStamp(i);
189 +      Atom* atomA;
190 +      Atom* atomB;
191 +      
192 +      atomA = mol->getAtomAt(cStamp->getA());
193 +      atomB = mol->getAtomAt(cStamp->getB());
194 +      assert( atomA && atomB );
195 +      
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"
205 +                "\twithout providing a value for the constrainedDistance.\n");
206 +        painCave.isFatal = 1;
207 +        simError();      
208 +      } else {
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:
220 +
221      createConstraintElem(mol);
222      
223      // Does this molecule stamp define a total constrained charge value?
# Line 186 | 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 266 | 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 274 | Line 316 | namespace OpenMD {
316      atomB = mol->getAtomAt(stamp->getB());
317      
318      assert( atomA && atomB);
277 –    
278 –    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 304 | 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(),
317 <                                           atomB->getType().c_str(),
318 <                                           atomC->getType().c_str());
319 <      
320 <      if (bendType == NULL) {
321 <        sprintf(painCave.errMsg,
322 <                "Can not find Matching Bend Type for[%s, %s, %s]",
323 <                atomA->getType().c_str(),
324 <                atomB->getType().c_str(),
325 <                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 339 | Line 419 | namespace OpenMD {
419          painCave.isFatal = 1;
420          simError();
421        }
342 –                
343 –      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 375 | 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 388 | 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 <    }
411 <    else {
527 >    } else {
528        
529        DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(stamp->getGhostVectorSource()));
530        if (dAtom == NULL) {
# Line 416 | Line 532 | namespace OpenMD {
532          painCave.isFatal = 1;
533          simError();
534        }        
535 <      
536 <      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
421 <                                                    atomC->getType(), "GHOST");
422 <      
423 <      if (torsionType == NULL) {
424 <        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
425 <                atomA->getType().c_str(),
426 <                atomB->getType().c_str(),
427 <                atomC->getType().c_str(),
428 <                "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 452 | 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 466 | Line 601 | namespace OpenMD {
601      Atom* atomD = mol->getAtomAt(satellites[2]);
602        
603      assert(atomA && atomB && atomC && atomD);
469 –    
470 –    InversionType* inversionType = ff->getInversionType(atomA->getType(),
471 –                                                        atomB->getType(),
472 –                                                        atomC->getType(),
473 –                                                        atomD->getType());
604  
605 <    if (inversionType == NULL) {
476 <      sprintf(painCave.errMsg, "No Matching Inversion Type for[%s, %s, %s, %s]\n"
477 <              "\t(May not be a problem: not all inversions are parametrized)\n",
478 <              atomA->getType().c_str(),
479 <              atomB->getType().c_str(),
480 <              atomC->getType().c_str(),
481 <              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 502 | 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    
# Line 544 | Line 698 | namespace OpenMD {
698      //add bond constraints
699      Molecule::BondIterator bi;
700      Bond* bond;
701 +    ConstraintPair* cPair;
702 +
703      for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
704          
705        BondType* bt = bond->getBondType();
# Line 553 | Line 709 | namespace OpenMD {
709  
710          ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
711          ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
712 <        ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
713 <        mol->addConstraintPair(consPair);
712 >        cPair = new ConstraintPair(consElemA, consElemB,
713 >                                   fbt->getEquilibriumBondLength(), false);
714 >        mol->addConstraintPair(cPair);
715        }
716      }
717  
# Line 566 | Line 723 | namespace OpenMD {
723      ConstraintPair* consPair;
724      Molecule::ConstraintPairIterator cpi;
725      std::set<StuntDouble*> sdSet;
726 <    for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
726 >    for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
727 >         consPair = mol->nextConstraintPair(cpi)) {
728  
729        StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
730        if (sdSet.find(sdA) == sdSet.end()){
731          sdSet.insert(sdA);
732          mol->addConstraintElem(new ConstraintElem(sdA));
733        }
734 <
734 >      
735        StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
736        if (sdSet.find(sdB) == sdSet.end()){
737          sdSet.insert(sdB);
738          mol->addConstraintElem(new ConstraintElem(sdB));
739 <      }
582 <        
739 >      }      
740      }
584 –
741    }
586 –    
742   }

Diff Legend

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