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 1290 by cli2, Wed Sep 10 19:51:45 2008 UTC vs.
Revision 1979 by gezelter, Sat Apr 5 20:56:01 2014 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 + * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 + * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42    
43   /**
44   * @file MoleculeCreator.cpp
45   * @author tlin
46   * @date 11/04/2004
46 * @time 13:44am
47   * @version 1.0
48   */
49  
# Line 54 | Line 54
54   #include "brains/MoleculeCreator.hpp"
55   #include "primitives/GhostBend.hpp"
56   #include "primitives/GhostTorsion.hpp"
57 < #include "types/DirectionalAtomType.hpp"
57 > #include "types/AtomType.hpp"
58   #include "types/FixedBondType.hpp"
59   #include "utils/simError.h"
60   #include "utils/StringUtils.hpp"
61  
62 < namespace oopse {
62 > namespace OpenMD {
63    
64    Molecule* MoleculeCreator::createMolecule(ForceField* ff,
65                                              MoleculeStamp *molStamp,
66                                              int stampId, int globalIndex,
67                                              LocalIndexManager* localIndexMan) {
68 <    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName());
69 <    
68 >    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName(),
69 >                                 molStamp->getRegion() );
70 >
71      //create atoms
72      Atom* atom;
73      AtomStamp* currentAtomStamp;
# Line 96 | Line 97 | namespace oopse {
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 106 | Line 107 | namespace oopse {
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 116 | Line 117 | namespace oopse {
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 126 | Line 127 | namespace oopse {
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 138 | Line 140 | namespace oopse {
140      int nCutoffGroups = molStamp->getNCutoffGroups();
141      for (int i = 0; i < nCutoffGroups; ++i) {
142        currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i);
143 <      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
143 >      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp,
144 >                                      localIndexMan);
145        mol->addCutoffGroup(cutoffGroup);
146      }
147  
# Line 169 | Line 172 | namespace oopse {
172      // every single free atom
173      
174      for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
175 <      cutoffGroup = createCutoffGroup(mol, *fai);
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?
221 +    // If so, let the created molecule know about it.
222 +
223 +    if (molStamp->haveConstrainTotalCharge() ) {
224 +      mol->setConstrainTotalCharge( molStamp->getConstrainTotalCharge() );
225 +    }
226 +
227      //the construction of this molecule is finished
228      mol->complete();
229      
# Line 198 | Line 246 | namespace oopse {
246        painCave.isFatal = 1;
247        simError();
248      }
249 <    
249 >
250      //below code still have some kind of hard-coding smell
251      if (atomType->isDirectional()){
204    
205      DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
206        
207      if (dAtomType == NULL) {
208        sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
252  
210        painCave.isFatal = 1;
211        simError();
212      }
213
253        DirectionalAtom* dAtom;
254 <      dAtom = new DirectionalAtom(dAtomType);
254 >      dAtom = new DirectionalAtom(atomType);
255        atom = dAtom;    
256      }
257      else{
# Line 227 | Line 266 | namespace oopse {
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 251 | Line 290 | namespace oopse {
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
259 <    /**@todo replace itoa by lexi_cast */
260 <    std::string s = OOPSE_itoa(mol->getNRigidBodies(), 10);
261 <    rb->setType(mol->getType() + "_RB_" + s.c_str());
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());
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 284 | Line 323 | namespace oopse {
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 303 | Line 356 | namespace oopse {
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 327 | Line 381 | namespace oopse {
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 339 | Line 394 | namespace oopse {
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 366 | Line 433 | namespace oopse {
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 403 | Line 471 | namespace oopse {
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 438 | Line 519 | namespace oopse {
519                atomD->getType().c_str());
520        
521        painCave.isFatal = 0;
522 <      painCave.severity = OOPSE_INFO;
522 >      painCave.severity = OPENMD_INFO;
523        simError();
524        return NULL;
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    }
545    
546  
547 <  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
547 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol,
548 >                                                  CutoffGroupStamp* stamp,
549 >                                                  LocalIndexManager* localIndexMan) {
550      int nAtoms;
551      CutoffGroup* cg;
552      Atom* atom;
# Line 461 | Line 558 | namespace oopse {
558        assert(atom);
559        cg->addAtom(atom);
560      }
561 <
561 >    
562 >    //set the local index of this cutoffGroup, global index will be set later
563 >    cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
564 >    
565      return cg;
566    }    
567 <
568 <  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
567 >  
568 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom,
569 >                                                  LocalIndexManager* localIndexMan) {
570      CutoffGroup* cg;
571      cg  = new CutoffGroup();
572      cg->addAtom(atom);
573 +
574 +    //set the local index of this cutoffGroup, global index will be set later
575 +    cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
576 +
577      return cg;
578    }
579  
# Line 477 | Line 582 | namespace oopse {
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  
484      //class Parent1 {};
485      //class Child1 : public Parent {};
486      //class Child2 : public Parent {};
487      //Child1* ch1 = new Child1();
488      //Child2* ch2 = dynamic_cast<Child2*>(ch1);
489      //the dynamic_cast is succeed in above line. A compiler bug?        
490
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 506 | Line 607 | namespace oopse {
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 <      }
522 <        
623 >      }      
624      }
524
625    }
526    
626   }

Comparing trunk/src/brains/MoleculeCreator.cpp (property svn:keywords):
Revision 1290 by cli2, Wed Sep 10 19:51:45 2008 UTC vs.
Revision 1979 by gezelter, Sat Apr 5 20:56:01 2014 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines