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 770 by tim, Fri Dec 2 15:38:03 2005 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 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  
50   #include <cassert>
51 + #include <typeinfo>
52   #include <set>
53  
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, 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 <    
69 >
70      //create atoms
71      Atom* atom;
72      AtomStamp* currentAtomStamp;
# Line 82 | 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 117 | 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;
138      int nCutoffGroups = molStamp->getNCutoffGroups();
139      for (int i = 0; i < nCutoffGroups; ++i) {
140        currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i);
141 <      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
141 >      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp, localIndexMan);
142        mol->addCutoffGroup(cutoffGroup);
143      }
144  
# Line 140 | 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);
172 >      cutoffGroup = createCutoffGroup(mol, *fai, localIndexMan);
173        mol->addCutoffGroup(cutoffGroup);
174      }
175      //create constraints
176      createConstraintPair(mol);
177      createConstraintElem(mol);
178      
179 +    // Does this molecule stamp define a total constrained charge value?
180 +    // If so, let the created molecule know about it.
181 +
182 +    if (molStamp->haveConstrainTotalCharge() ) {
183 +      mol->setConstrainTotalCharge( molStamp->getConstrainTotalCharge() );
184 +    }
185 +
186      //the construction of this molecule is finished
187      mol->complete();
188 <
188 >    
189      return mol;
190    }    
191  
192  
193 <  Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp,
194 <                                    LocalIndexManager* localIndexMan) {
193 >  Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol,
194 >                                    AtomStamp* stamp,
195 >                                    LocalIndexManager* localIndexMan) {
196      AtomType * atomType;
197      Atom* atom;
198  
199      atomType =  ff->getAtomType(stamp->getType());
200 <
200 >    
201      if (atomType == NULL) {
202        sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
203                stamp->getType().c_str());
# Line 180 | Line 205 | namespace oopse {
205        painCave.isFatal = 1;
206        simError();
207      }
208 <    
208 >
209      //below code still have some kind of hard-coding smell
210      if (atomType->isDirectional()){
186    
187      DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
188        
189      if (dAtomType == NULL) {
190        sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
211  
192        painCave.isFatal = 1;
193        simError();
194      }
195
212        DirectionalAtom* dAtom;
213 <      dAtom = new DirectionalAtom(dAtomType);
213 >      dAtom = new DirectionalAtom(atomType);
214        atom = dAtom;    
215      }
216      else{
# Line 205 | Line 221 | namespace oopse {
221  
222      return atom;
223    }
224 <
225 <  RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
224 >  
225 >  RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp,
226 >                                              Molecule* mol,
227                                                RigidBodyStamp* rbStamp,
228                                                LocalIndexManager* localIndexMan) {
229      Atom* atom;
# Line 217 | Line 234 | namespace oopse {
234      RigidBody* rb = new RigidBody();
235      nAtoms = rbStamp->getNMembers();    
236      for (int i = 0; i < nAtoms; ++i) {
237 <      //rbStamp->getMember(i) return the local index of current atom inside the molecule.
238 <      //It is not the same as local index of atom which is the index of atom at DataStorage class
237 >      //rbStamp->getMember(i) return the local index of current atom
238 >      //inside the molecule.  It is not the same as local index of
239 >      //atom which is the index of atom at DataStorage class
240        atom = mol->getAtomAt(rbStamp->getMemberAt(i));
241        atomStamp= molStamp->getAtomStamp(rbStamp->getMemberAt(i));    
242        rb->addAtom(atom, atomStamp);
243      }
244  
245 <    //after all of the atoms are added, we need to calculate the reference coordinates
245 >    //after all of the atoms are added, we need to calculate the
246 >    //reference coordinates
247      rb->calcRefCoords();
248  
249      //set the local index of this rigid body, global index will be set later
# Line 236 | Line 255 | namespace oopse {
255      //The third part is the index of the rigidbody defined in meta-data file
256      //For example, Butane_RB_0 is a valid rigid body name of butane molecule
257      /**@todo replace itoa by lexi_cast */
258 <    std::string s = OOPSE_itoa(mol->getNRigidBodies(), 10);
258 >    std::string s = OpenMD_itoa(mol->getNRigidBodies(), 10);
259      rb->setType(mol->getType() + "_RB_" + s.c_str());
260  
261      return rb;
262    }    
263  
264 <  Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
264 >  Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
265 >                                    BondStamp* stamp) {
266      BondType* bondType;
267      Atom* atomA;
268      Atom* atomB;
269 <
269 >    
270      atomA = mol->getAtomAt(stamp->getA());
271      atomB = mol->getAtomAt(stamp->getB());
272 <
272 >    
273      assert( atomA && atomB);
274      
275      bondType = ff->getBondType(atomA->getType(), atomB->getType());
# Line 258 | Line 278 | namespace oopse {
278        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
279                atomA->getType().c_str(),
280                atomB->getType().c_str());
281 <
281 >      
282        painCave.isFatal = 1;
283        simError();
284      }
285      return new Bond(atomA, atomB, bondType);    
286    }    
287 <
288 <  Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
287 >  
288 >  Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol,
289 >                                    BendStamp* stamp) {
290      Bend* bend = NULL;
291      std::vector<int> bendAtoms = stamp->getMembers();
292      if (bendAtoms.size() == 3) {
293        Atom* atomA = mol->getAtomAt(bendAtoms[0]);
294        Atom* atomB = mol->getAtomAt(bendAtoms[1]);
295        Atom* atomC = mol->getAtomAt(bendAtoms[2]);
296 <
296 >      
297        assert( atomA && atomB && atomC);
298 <
299 <      BendType* bendType = ff->getBendType(atomA->getType().c_str(), atomB->getType().c_str(), atomC->getType().c_str());
300 <
298 >      
299 >      BendType* bendType = ff->getBendType(atomA->getType().c_str(),
300 >                                           atomB->getType().c_str(),
301 >                                           atomC->getType().c_str());
302 >      
303        if (bendType == NULL) {
304          sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
305                  atomA->getType().c_str(),
306                  atomB->getType().c_str(),
307                  atomC->getType().c_str());
308 <
308 >        
309          painCave.isFatal = 1;
310          simError();
311        }
312 <
312 >      
313        bend = new Bend(atomA, atomB, atomC, bendType);
314      } else if ( bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
315        int ghostIndex = stamp->getGhostVectorSource();
# Line 310 | Line 333 | namespace oopse {
333          painCave.isFatal = 1;
334          simError();
335        }
336 <        
336 >      
337        bend = new GhostBend(normalAtom, ghostAtom, bendType);      
338 <
338 >      
339      }
340      
341      return bend;
342    }    
343  
344 <  Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
344 >  Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
345 >                                          TorsionStamp* stamp) {
346  
347      Torsion* torsion = NULL;
348      std::vector<int> torsionAtoms = stamp->getMembers();
# Line 335 | Line 359 | namespace oopse {
359  
360        assert(atomA && atomB && atomC && atomD);
361          
362 <      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
363 <                                                    atomC->getType(), atomD->getType());
364 <
362 >      TorsionType* torsionType = ff->getTorsionType(atomA->getType(),
363 >                                                    atomB->getType(),
364 >                                                    atomC->getType(),
365 >                                                    atomD->getType());
366        if (torsionType == NULL) {
367          sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
368                  atomA->getType().c_str(),
369                  atomB->getType().c_str(),
370                  atomC->getType().c_str(),
371                  atomD->getType().c_str());
372 <
372 >        
373          painCave.isFatal = 1;
374          simError();
375        }
376 <        
376 >      
377        torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
378      }
379      else {
380 <
380 >      
381        DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(stamp->getGhostVectorSource()));
382        if (dAtom == NULL) {
383          sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
384          painCave.isFatal = 1;
385          simError();
386        }        
387 <
387 >      
388        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
389                                                      atomC->getType(), "GHOST");
390 <
390 >      
391        if (torsionType == NULL) {
392          sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
393                  atomA->getType().c_str(),
394                  atomB->getType().c_str(),
395                  atomC->getType().c_str(),
396                  "GHOST");
397 <
397 >        
398          painCave.isFatal = 1;
399          simError();
400        }
401 <        
401 >      
402        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
403      }
404 <
404 >    
405      return torsion;
406    }    
407  
408 <  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
408 >  Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol,
409 >                                              InversionStamp* stamp) {
410 >    
411 >    Inversion* inversion = NULL;
412 >    int center = stamp->getCenter();
413 >    std::vector<int> satellites = stamp->getSatellites();
414 >    if (satellites.size() != 3) {
415 >        return inversion;
416 >    }
417 >
418 >    Atom* atomA = mol->getAtomAt(center);
419 >    Atom* atomB = mol->getAtomAt(satellites[0]);
420 >    Atom* atomC = mol->getAtomAt(satellites[1]);
421 >    Atom* atomD = mol->getAtomAt(satellites[2]);
422 >      
423 >    assert(atomA && atomB && atomC && atomD);
424 >    
425 >    InversionType* inversionType = ff->getInversionType(atomA->getType(),
426 >                                                        atomB->getType(),
427 >                                                        atomC->getType(),
428 >                                                        atomD->getType());
429 >
430 >    if (inversionType == NULL) {
431 >      sprintf(painCave.errMsg, "No Matching Inversion Type for[%s, %s, %s, %s]\n"
432 >              "\t(May not be a problem: not all inversions are parametrized)\n",
433 >              atomA->getType().c_str(),
434 >              atomB->getType().c_str(),
435 >              atomC->getType().c_str(),
436 >              atomD->getType().c_str());
437 >      
438 >      painCave.isFatal = 0;
439 >      painCave.severity = OPENMD_INFO;
440 >      simError();
441 >      return NULL;
442 >    } else {
443 >      
444 >      inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType);
445 >      return inversion;
446 >    }
447 >  }
448 >  
449 >
450 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol,
451 >                                                  CutoffGroupStamp* stamp,
452 >                                                  LocalIndexManager* localIndexMan) {
453      int nAtoms;
454      CutoffGroup* cg;
455      Atom* atom;
# Line 392 | Line 461 | namespace oopse {
461        assert(atom);
462        cg->addAtom(atom);
463      }
464 <
464 >    
465 >    //set the local index of this cutoffGroup, global index will be set later
466 >    cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
467 >    
468      return cg;
469    }    
470 <
471 <  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
470 >  
471 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom,
472 >                                                  LocalIndexManager* localIndexMan) {
473      CutoffGroup* cg;
474      cg  = new CutoffGroup();
475      cg->addAtom(atom);
476 +
477 +    //set the local index of this cutoffGroup, global index will be set later
478 +    cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
479 +
480      return cg;
481    }
482  

Comparing trunk/src/brains/MoleculeCreator.cpp (property svn:keywords):
Revision 770 by tim, Fri Dec 2 15:38:03 2005 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines