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 1908 by gezelter, Fri Jul 19 21:25:45 2013 UTC vs.
Revision 1979 by gezelter, Sat Apr 5 20:56:01 2014 UTC

# Line 97 | Line 97 | namespace OpenMD {
97  
98      for (int i = 0; i < nBonds; ++i) {
99        currentBondStamp = molStamp->getBondStamp(i);
100 <      bond = createBond(ff, mol, currentBondStamp);
100 >      bond = createBond(ff, mol, currentBondStamp, localIndexMan);
101        mol->addBond(bond);
102      }
103  
# Line 107 | Line 107 | namespace OpenMD {
107      int nBends = molStamp->getNBends();
108      for (int i = 0; i < nBends; ++i) {
109        currentBendStamp = molStamp->getBendStamp(i);
110 <      bend = createBend(ff, mol, currentBendStamp);
110 >      bend = createBend(ff, mol, currentBendStamp, localIndexMan);
111        mol->addBend(bend);
112      }
113  
# Line 117 | Line 117 | namespace OpenMD {
117      int nTorsions = molStamp->getNTorsions();
118      for (int i = 0; i < nTorsions; ++i) {
119        currentTorsionStamp = molStamp->getTorsionStamp(i);
120 <      torsion = createTorsion(ff, mol, currentTorsionStamp);
120 >      torsion = createTorsion(ff, mol, currentTorsionStamp, localIndexMan);
121        mol->addTorsion(torsion);
122      }
123  
# Line 127 | Line 127 | namespace OpenMD {
127      int nInversions = molStamp->getNInversions();
128      for (int i = 0; i < nInversions; ++i) {
129        currentInversionStamp = molStamp->getInversionStamp(i);
130 <      inversion = createInversion(ff, mol, currentInversionStamp);
130 >      inversion = createInversion(ff, mol, currentInversionStamp,
131 >                                  localIndexMan);
132        if (inversion != NULL ) {
133          mol->addInversion(inversion);
134        }
# Line 139 | Line 140 | namespace OpenMD {
140      int nCutoffGroups = molStamp->getNCutoffGroups();
141      for (int i = 0; i < nCutoffGroups; ++i) {
142        currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i);
143 <      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp, localIndexMan);
143 >      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp,
144 >                                      localIndexMan);
145        mol->addCutoffGroup(cutoffGroup);
146      }
147  
# Line 173 | Line 175 | namespace OpenMD {
175        cutoffGroup = createCutoffGroup(mol, *fai, localIndexMan);
176        mol->addCutoffGroup(cutoffGroup);
177      }
178 <    //create constraints
178 >
179 >    //create bonded constraintPairs:
180      createConstraintPair(mol);
181 +
182 +    //create non-bonded constraintPairs
183 +    for (int i = 0; i < molStamp->getNConstraints(); ++i) {
184 +      ConstraintStamp* cStamp = molStamp->getConstraintStamp(i);
185 +      Atom* atomA;
186 +      Atom* atomB;
187 +      
188 +      atomA = mol->getAtomAt(cStamp->getA());
189 +      atomB = mol->getAtomAt(cStamp->getB());
190 +      assert( atomA && atomB );
191 +      
192 +      RealType distance;
193 +      bool printConstraintForce = false;
194 +
195 +      if (!cStamp->haveConstrainedDistance()) {
196 +        sprintf(painCave.errMsg,
197 +                "Constraint Error: A non-bond constraint was specified\n"
198 +                "\twithout providing a value for the constrainedDistance.\n");
199 +        painCave.isFatal = 1;
200 +        simError();      
201 +      } else {
202 +        distance = cStamp->getConstrainedDistance();
203 +      }
204 +
205 +      if (cStamp->havePrintConstraintForce()) {
206 +        printConstraintForce = cStamp->getPrintConstraintForce();
207 +      }
208 +    
209 +      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);
214 +    }
215 +
216 +    // now create the constraint elements:
217 +
218      createConstraintElem(mol);
219      
220      // Does this molecule stamp define a total constrained charge value?
# Line 226 | Line 266 | namespace OpenMD {
266    RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp,
267                                                Molecule* mol,
268                                                RigidBodyStamp* rbStamp,
269 <                                              LocalIndexManager* localIndexMan) {
269 >                                              LocalIndexManager* localIndexMan){
270      Atom* atom;
271      int nAtoms;
272      Vector3d refCoor;
# Line 250 | Line 290 | namespace OpenMD {
290      //set the local index of this rigid body, global index will be set later
291      rb->setLocalIndex(localIndexMan->getNextRigidBodyIndex());
292  
293 <    //the rule for naming rigidbody MoleculeName_RB_Integer
294 <    //The first part is the name of the molecule
295 <    //The second part is alway fixed as "RB"
296 <    //The third part is the index of the rigidbody defined in meta-data file
297 <    //For example, Butane_RB_0 is a valid rigid body name of butane molecule
298 <    /**@todo replace itoa by lexi_cast */
293 >    // The rule for naming a rigidbody is: MoleculeName_RB_Integer
294 >    // The first part is the name of the molecule
295 >    // The second part is always fixed as "RB"
296 >    // The third part is the index of the rigidbody defined in meta-data file
297 >    // For example, Butane_RB_0 is a valid rigid body name of butane molecule
298 >
299      std::string s = OpenMD_itoa(mol->getNRigidBodies(), 10);
300      rb->setType(mol->getType() + "_RB_" + s.c_str());
261 –
301      return rb;
302    }    
303  
304    Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
305 <                                    BondStamp* stamp) {
305 >                                    BondStamp* stamp,
306 >                                    LocalIndexManager* localIndexMan) {
307      BondType* bondType;
308      Atom* atomA;
309      Atom* atomB;
# Line 283 | Line 323 | namespace OpenMD {
323        painCave.isFatal = 1;
324        simError();
325      }
326 <    return new Bond(atomA, atomB, bondType);    
326 >    Bond* bond = new Bond(atomA, atomB, bondType);
327 >
328 >    //set the local index of this bond, the global index will be set later
329 >    bond->setLocalIndex(localIndexMan->getNextBondIndex());
330 >
331 >    // The rule for naming a bond is: MoleculeName_Bond_Integer
332 >    // The first part is the name of the molecule
333 >    // The second part is always fixed as "Bond"
334 >    // The third part is the index of the bond defined in meta-data file
335 >    // For example, Butane_bond_0 is a valid Bond name in a butane molecule
336 >
337 >    std::string s = OpenMD_itoa(mol->getNBonds(), 10);
338 >    bond->setName(mol->getType() + "_Bond_" + s.c_str());
339 >    return bond;    
340    }    
341    
342    Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol,
343 <                                    BendStamp* stamp) {
343 >                                    BendStamp* stamp,
344 >                                    LocalIndexManager* localIndexMan) {
345      Bend* bend = NULL;
346      std::vector<int> bendAtoms = stamp->getMembers();
347      if (bendAtoms.size() == 3) {
# Line 302 | Line 356 | namespace OpenMD {
356                                             atomC->getType().c_str());
357        
358        if (bendType == NULL) {
359 <        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
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());
# Line 326 | Line 381 | namespace OpenMD {
381        BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
382  
383        if (bendType == NULL) {
384 <        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
384 >        sprintf(painCave.errMsg,
385 >                "Can not find Matching Bend Type for[%s, %s, %s]",
386                  normalAtom->getType().c_str(),
387                  ghostAtom->getType().c_str(),
388                  "GHOST");
# Line 338 | Line 394 | namespace OpenMD {
394        bend = new GhostBend(normalAtom, ghostAtom, bendType);      
395        
396      }
397 <    
397 >
398 >    //set the local index of this bend, the global index will be set later
399 >    bend->setLocalIndex(localIndexMan->getNextBendIndex());
400 >
401 >    // The rule for naming a bend is: MoleculeName_Bend_Integer
402 >    // The first part is the name of the molecule
403 >    // The second part is always fixed as "Bend"
404 >    // The third part is the index of the bend defined in meta-data file
405 >    // For example, Butane_Bend_0 is a valid Bend name in a butane molecule
406 >
407 >    std::string s = OpenMD_itoa(mol->getNBends(), 10);
408 >    bend->setName(mol->getType() + "_Bend_" + s.c_str());    
409      return bend;
410    }    
411  
412    Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
413 <                                          TorsionStamp* stamp) {
413 >                                          TorsionStamp* stamp,
414 >                                          LocalIndexManager* localIndexMan) {
415  
416      Torsion* torsion = NULL;
417      std::vector<int> torsionAtoms = stamp->getMembers();
# Line 365 | Line 433 | namespace OpenMD {
433                                                      atomC->getType(),
434                                                      atomD->getType());
435        if (torsionType == NULL) {
436 <        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
436 >        sprintf(painCave.errMsg,
437 >                "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
438                  atomA->getType().c_str(),
439                  atomB->getType().c_str(),
440                  atomC->getType().c_str(),
# Line 402 | Line 471 | namespace OpenMD {
471        
472        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
473      }
474 +
475 +    //set the local index of this torsion, the global index will be set later
476 +    torsion->setLocalIndex(localIndexMan->getNextTorsionIndex());
477      
478 +    // The rule for naming a torsion is: MoleculeName_Torsion_Integer
479 +    // The first part is the name of the molecule
480 +    // The second part is always fixed as "Torsion"
481 +    // The third part is the index of the torsion defined in meta-data file
482 +    // For example, Butane_Torsion_0 is a valid Torsion name in a
483 +    // butane molecule
484 +
485 +    std::string s = OpenMD_itoa(mol->getNTorsions(), 10);
486 +    torsion->setName(mol->getType() + "_Torsion_" + s.c_str());
487      return torsion;
488    }    
489  
490    Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol,
491 <                                              InversionStamp* stamp) {
491 >                                              InversionStamp* stamp,
492 >                                              LocalIndexManager* localIndexMan) {
493      
494      Inversion* inversion = NULL;
495      int center = stamp->getCenter();
# Line 443 | Line 525 | namespace OpenMD {
525      } else {
526        
527        inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType);
528 +
529 +      // set the local index of this inversion, the global index will
530 +      // be set later
531 +      inversion->setLocalIndex(localIndexMan->getNextInversionIndex());
532 +
533 +      // The rule for naming an inversion is: MoleculeName_Inversion_Integer
534 +      // The first part is the name of the molecule
535 +      // The second part is always fixed as "Inversion"
536 +      // The third part is the index of the inversion defined in meta-data file
537 +      // For example, Benzene_Inversion_0 is a valid Inversion name in a
538 +      // Benzene molecule
539 +
540 +      std::string s = OpenMD_itoa(mol->getNInversions(), 10);
541 +      inversion->setName(mol->getType() + "_Inversion_" + s.c_str());
542        return inversion;
543      }
544    }
# Line 486 | Line 582 | namespace OpenMD {
582      //add bond constraints
583      Molecule::BondIterator bi;
584      Bond* bond;
585 +    ConstraintPair* cPair;
586 +
587      for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
588          
589        BondType* bt = bond->getBondType();
590  
493 –      //class Parent1 {};
494 –      //class Child1 : public Parent {};
495 –      //class Child2 : public Parent {};
496 –      //Child1* ch1 = new Child1();
497 –      //Child2* ch2 = dynamic_cast<Child2*>(ch1);
498 –      //the dynamic_cast is succeed in above line. A compiler bug?        
499 –
591        if (typeid(FixedBondType) == typeid(*bt)) {
592          FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
593  
594          ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
595          ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
596 <        ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
597 <        mol->addConstraintPair(consPair);
596 >        cPair = new ConstraintPair(consElemA, consElemB,
597 >                                   fbt->getEquilibriumBondLength(), false);
598 >        mol->addConstraintPair(cPair);
599        }
600      }
601  
# Line 515 | Line 607 | namespace OpenMD {
607      ConstraintPair* consPair;
608      Molecule::ConstraintPairIterator cpi;
609      std::set<StuntDouble*> sdSet;
610 <    for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
610 >    for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
611 >         consPair = mol->nextConstraintPair(cpi)) {
612  
613        StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
614        if (sdSet.find(sdA) == sdSet.end()){
615          sdSet.insert(sdA);
616          mol->addConstraintElem(new ConstraintElem(sdA));
617        }
618 <
618 >      
619        StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
620        if (sdSet.find(sdB) == sdSet.end()){
621          sdSet.insert(sdB);
622          mol->addConstraintElem(new ConstraintElem(sdB));
623 <      }
531 <        
623 >      }      
624      }
533 –
625    }
535 –    
626   }

Diff Legend

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