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 809 by gezelter, Mon Dec 12 19:32:50 2005 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 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, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41    
42   /**
# Line 59 | Line 59
59   #include "utils/simError.h"
60   #include "utils/StringUtils.hpp"
61  
62 < namespace oopse {
62 > namespace OpenMD {
63    
64 <  Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp,
65 <                                            int stampId, int globalIndex, LocalIndexManager* localIndexMan) {
66 <
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      
70      //create atoms
# Line 83 | Line 84 | namespace oopse {
84  
85      for (int i = 0; i < nRigidbodies; ++i) {
86        currentRigidBodyStamp = molStamp->getRigidBodyStamp(i);
87 <      rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
87 >      rb = createRigidBody(molStamp, mol, currentRigidBodyStamp,
88 >                           localIndexMan);
89        mol->addRigidBody(rb);
90      }
91 <
91 >    
92      //create bonds
93      Bond* bond;
94      BondStamp* currentBondStamp;
# Line 118 | Line 120 | namespace oopse {
120        mol->addTorsion(torsion);
121      }
122  
123 +    //create inversions
124 +    Inversion* inversion;
125 +    InversionStamp* currentInversionStamp;
126 +    int nInversions = molStamp->getNInversions();
127 +    for (int i = 0; i < nInversions; ++i) {
128 +      currentInversionStamp = molStamp->getInversionStamp(i);
129 +      inversion = createInversion(ff, mol, currentInversionStamp);
130 +      if (inversion != NULL ) {
131 +        mol->addInversion(inversion);
132 +      }
133 +    }
134 +
135      //create cutoffGroups
136      CutoffGroup* cutoffGroup;
137      CutoffGroupStamp* currentCutoffGroupStamp;
# Line 141 | Line 155 | namespace oopse {
155      Molecule::CutoffGroupIterator ci;
156      CutoffGroup* cg;
157      
158 <    for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
159 <
158 >    for (cg = mol->beginCutoffGroup(ci); cg != NULL;
159 >         cg = mol->nextCutoffGroup(ci)) {
160 >      
161        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
162 <           //erase the atoms belong to cutoff groups from freeAtoms vector
163 <           freeAtoms.erase(std::remove(freeAtoms.begin(), freeAtoms.end(), atom), freeAtoms.end());
164 <      }
165 <
162 >        //erase the atoms belong to cutoff groups from freeAtoms vector
163 >        freeAtoms.erase(std::remove(freeAtoms.begin(), freeAtoms.end(), atom),
164 >                        freeAtoms.end());
165 >      }      
166      }      
167      
168 <    //loop over the free atoms and then create one cutoff group for every single free atom
168 >    // loop over the free atoms and then create one cutoff group for
169 >    // every single free atom
170      
171      for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
172        cutoffGroup = createCutoffGroup(mol, *fai);
# Line 162 | Line 178 | namespace oopse {
178      
179      //the construction of this molecule is finished
180      mol->complete();
181 <
181 >    
182      return mol;
183    }    
184  
185  
186 <  Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp,
187 <                                    LocalIndexManager* localIndexMan) {
186 >  Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol,
187 >                                    AtomStamp* stamp,
188 >                                    LocalIndexManager* localIndexMan) {
189      AtomType * atomType;
190      Atom* atom;
191  
192      atomType =  ff->getAtomType(stamp->getType());
193 <
193 >    
194      if (atomType == NULL) {
195        sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
196                stamp->getType().c_str());
# Line 206 | Line 223 | namespace oopse {
223  
224      return atom;
225    }
226 <
227 <  RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
226 >  
227 >  RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp,
228 >                                              Molecule* mol,
229                                                RigidBodyStamp* rbStamp,
230                                                LocalIndexManager* localIndexMan) {
231      Atom* atom;
# Line 218 | Line 236 | namespace oopse {
236      RigidBody* rb = new RigidBody();
237      nAtoms = rbStamp->getNMembers();    
238      for (int i = 0; i < nAtoms; ++i) {
239 <      //rbStamp->getMember(i) return the local index of current atom inside the molecule.
240 <      //It is not the same as local index of atom which is the index of atom at DataStorage class
239 >      //rbStamp->getMember(i) return the local index of current atom
240 >      //inside the molecule.  It is not the same as local index of
241 >      //atom which is the index of atom at DataStorage class
242        atom = mol->getAtomAt(rbStamp->getMemberAt(i));
243        atomStamp= molStamp->getAtomStamp(rbStamp->getMemberAt(i));    
244        rb->addAtom(atom, atomStamp);
245      }
246  
247 <    //after all of the atoms are added, we need to calculate the reference coordinates
247 >    //after all of the atoms are added, we need to calculate the
248 >    //reference coordinates
249      rb->calcRefCoords();
250  
251      //set the local index of this rigid body, global index will be set later
# Line 237 | Line 257 | namespace oopse {
257      //The third part is the index of the rigidbody defined in meta-data file
258      //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);
260 >    std::string s = OpenMD_itoa(mol->getNRigidBodies(), 10);
261      rb->setType(mol->getType() + "_RB_" + s.c_str());
262  
263      return rb;
264    }    
265  
266 <  Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
266 >  Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
267 >                                    BondStamp* stamp) {
268      BondType* bondType;
269      Atom* atomA;
270      Atom* atomB;
271 <
271 >    
272      atomA = mol->getAtomAt(stamp->getA());
273      atomB = mol->getAtomAt(stamp->getB());
274 <
274 >    
275      assert( atomA && atomB);
276      
277      bondType = ff->getBondType(atomA->getType(), atomB->getType());
# Line 259 | Line 280 | namespace oopse {
280        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
281                atomA->getType().c_str(),
282                atomB->getType().c_str());
283 <
283 >      
284        painCave.isFatal = 1;
285        simError();
286      }
287      return new Bond(atomA, atomB, bondType);    
288    }    
289 <
290 <  Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
289 >  
290 >  Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol,
291 >                                    BendStamp* stamp) {
292      Bend* bend = NULL;
293      std::vector<int> bendAtoms = stamp->getMembers();
294      if (bendAtoms.size() == 3) {
295        Atom* atomA = mol->getAtomAt(bendAtoms[0]);
296        Atom* atomB = mol->getAtomAt(bendAtoms[1]);
297        Atom* atomC = mol->getAtomAt(bendAtoms[2]);
298 <
298 >      
299        assert( atomA && atomB && atomC);
300 <
301 <      BendType* bendType = ff->getBendType(atomA->getType().c_str(), atomB->getType().c_str(), atomC->getType().c_str());
302 <
300 >      
301 >      BendType* bendType = ff->getBendType(atomA->getType().c_str(),
302 >                                           atomB->getType().c_str(),
303 >                                           atomC->getType().c_str());
304 >      
305        if (bendType == NULL) {
306          sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
307                  atomA->getType().c_str(),
308                  atomB->getType().c_str(),
309                  atomC->getType().c_str());
310 <
310 >        
311          painCave.isFatal = 1;
312          simError();
313        }
314 <
314 >      
315        bend = new Bend(atomA, atomB, atomC, bendType);
316      } else if ( bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
317        int ghostIndex = stamp->getGhostVectorSource();
# Line 311 | Line 335 | namespace oopse {
335          painCave.isFatal = 1;
336          simError();
337        }
338 <        
338 >      
339        bend = new GhostBend(normalAtom, ghostAtom, bendType);      
340 <
340 >      
341      }
342      
343      return bend;
344    }    
345  
346 <  Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
346 >  Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
347 >                                          TorsionStamp* stamp) {
348  
349      Torsion* torsion = NULL;
350      std::vector<int> torsionAtoms = stamp->getMembers();
# Line 336 | Line 361 | namespace oopse {
361  
362        assert(atomA && atomB && atomC && atomD);
363          
364 <      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
365 <                                                    atomC->getType(), atomD->getType());
366 <
364 >      TorsionType* torsionType = ff->getTorsionType(atomA->getType(),
365 >                                                    atomB->getType(),
366 >                                                    atomC->getType(),
367 >                                                    atomD->getType());
368        if (torsionType == NULL) {
369          sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
370                  atomA->getType().c_str(),
371                  atomB->getType().c_str(),
372                  atomC->getType().c_str(),
373                  atomD->getType().c_str());
374 <
374 >        
375          painCave.isFatal = 1;
376          simError();
377        }
378 <        
378 >      
379        torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
380      }
381      else {
382 <
382 >      
383        DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(stamp->getGhostVectorSource()));
384        if (dAtom == NULL) {
385          sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
386          painCave.isFatal = 1;
387          simError();
388        }        
389 <
389 >      
390        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
391                                                      atomC->getType(), "GHOST");
392 <
392 >      
393        if (torsionType == NULL) {
394          sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
395                  atomA->getType().c_str(),
396                  atomB->getType().c_str(),
397                  atomC->getType().c_str(),
398                  "GHOST");
399 <
399 >        
400          painCave.isFatal = 1;
401          simError();
402        }
403 <        
403 >      
404        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
405      }
406 <
406 >    
407      return torsion;
408    }    
409  
410 +  Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol,
411 +                                              InversionStamp* stamp) {
412 +    
413 +    Inversion* inversion = NULL;
414 +    int center = stamp->getCenter();
415 +    std::vector<int> satellites = stamp->getSatellites();
416 +    if (satellites.size() != 3) {
417 +        return inversion;
418 +    }
419 +
420 +    Atom* atomA = mol->getAtomAt(center);
421 +    Atom* atomB = mol->getAtomAt(satellites[0]);
422 +    Atom* atomC = mol->getAtomAt(satellites[1]);
423 +    Atom* atomD = mol->getAtomAt(satellites[2]);
424 +      
425 +    assert(atomA && atomB && atomC && atomD);
426 +    
427 +    InversionType* inversionType = ff->getInversionType(atomA->getType(),
428 +                                                        atomB->getType(),
429 +                                                        atomC->getType(),
430 +                                                        atomD->getType());
431 +
432 +    if (inversionType == NULL) {
433 +      sprintf(painCave.errMsg, "No Matching Inversion Type for[%s, %s, %s, %s]\n"
434 +              "\t(May not be a problem: not all inversions are parametrized)\n",
435 +              atomA->getType().c_str(),
436 +              atomB->getType().c_str(),
437 +              atomC->getType().c_str(),
438 +              atomD->getType().c_str());
439 +      
440 +      painCave.isFatal = 0;
441 +      painCave.severity = OPENMD_INFO;
442 +      simError();
443 +      return NULL;
444 +    } else {
445 +      
446 +      inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType);
447 +      return inversion;
448 +    }
449 +  }
450 +  
451 +
452    CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
453      int nAtoms;
454      CutoffGroup* cg;

Diff Legend

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