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 273 by tim, Tue Jan 25 17:45:23 2005 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# 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 < /**
43 <  * @file MoleculeCreator.cpp
44 <  * @author tlin
45 <  * @date 11/04/2004
46 <  * @time 13:44am
47 <  * @version 1.0
48 <  */
42 > /**
43 > * @file MoleculeCreator.cpp
44 > * @author tlin
45 > * @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"
# Line 58 | Line 59
59   #include "utils/simError.h"
60   #include "utils/StringUtils.hpp"
61  
62 < namespace oopse {
63 <
64 < Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp,
65 <    int stampId, int globalIndex, LocalIndexManager* localIndexMan) {
66 <
67 <    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getID());
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      
70      //create atoms
71      Atom* atom;
72      AtomStamp* currentAtomStamp;
73      int nAtom = molStamp->getNAtoms();
74      for (int i = 0; i < nAtom; ++i) {
75 <        currentAtomStamp = molStamp->getAtom(i);
76 <        atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
77 <        mol->addAtom(atom);
75 >      currentAtomStamp = molStamp->getAtomStamp(i);
76 >      atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
77 >      mol->addAtom(atom);
78      }
79  
80      //create rigidbodies
# Line 81 | Line 83 | Molecule* MoleculeCreator::createMolecule(ForceField*
83      int nRigidbodies = molStamp->getNRigidBodies();
84  
85      for (int i = 0; i < nRigidbodies; ++i) {
86 <        currentRigidBodyStamp = molStamp->getRigidBody(i);
87 <        rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
88 <        mol->addRigidBody(rb);
86 >      currentRigidBodyStamp = molStamp->getRigidBodyStamp(i);
87 >      rb = createRigidBody(molStamp, mol, currentRigidBodyStamp,
88 >                           localIndexMan);
89 >      mol->addRigidBody(rb);
90      }
91 <
91 >    
92      //create bonds
93      Bond* bond;
94      BondStamp* currentBondStamp;
95      int nBonds = molStamp->getNBonds();
96  
97      for (int i = 0; i < nBonds; ++i) {
98 <        currentBondStamp = molStamp->getBond(i);
99 <        bond = createBond(ff, mol, currentBondStamp);
100 <        mol->addBond(bond);
98 >      currentBondStamp = molStamp->getBondStamp(i);
99 >      bond = createBond(ff, mol, currentBondStamp);
100 >      mol->addBond(bond);
101      }
102  
103      //create bends
# Line 102 | Line 105 | Molecule* MoleculeCreator::createMolecule(ForceField*
105      BendStamp* currentBendStamp;
106      int nBends = molStamp->getNBends();
107      for (int i = 0; i < nBends; ++i) {
108 <        currentBendStamp = molStamp->getBend(i);
109 <        bend = createBend(ff, mol, currentBendStamp);
110 <        mol->addBend(bend);
108 >      currentBendStamp = molStamp->getBendStamp(i);
109 >      bend = createBend(ff, mol, currentBendStamp);
110 >      mol->addBend(bend);
111      }
112  
113      //create torsions
# Line 112 | Line 115 | Molecule* MoleculeCreator::createMolecule(ForceField*
115      TorsionStamp* currentTorsionStamp;
116      int nTorsions = molStamp->getNTorsions();
117      for (int i = 0; i < nTorsions; ++i) {
118 <        currentTorsionStamp = molStamp->getTorsion(i);
119 <        torsion = createTorsion(ff, mol, currentTorsionStamp);
120 <        mol->addTorsion(torsion);
118 >      currentTorsionStamp = molStamp->getTorsionStamp(i);
119 >      torsion = createTorsion(ff, mol, currentTorsionStamp);
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->getCutoffGroup(i);
141 <        cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
142 <        mol->addCutoffGroup(cutoffGroup);
140 >      currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i);
141 >      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
142 >      mol->addCutoffGroup(cutoffGroup);
143      }
144  
145      //every free atom is a cutoff group    
146 <    std::set<Atom*> allAtoms;
147 <     Molecule::AtomIterator ai;
146 >    std::vector<Atom*> freeAtoms;
147 >    std::vector<Atom*>::iterator ai;
148 >    std::vector<Atom*>::iterator fai;
149  
150      //add all atoms into allAtoms set
151 <    for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
152 <        allAtoms.insert(atom);
151 >    for(atom = mol->beginAtom(fai); atom != NULL; atom = mol->nextAtom(fai)) {
152 >      freeAtoms.push_back(atom);
153      }
154  
155      Molecule::CutoffGroupIterator ci;
156      CutoffGroup* cg;
141 –    std::set<Atom*> cutoffAtoms;    
157      
158 <    //add all of the atoms belong to cutoff groups into cutoffAtoms set
159 <    for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
160 <
161 <        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
162 <            cutoffAtoms.insert(atom);
163 <        }
164 <
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),
164 >                        freeAtoms.end());
165 >      }      
166      }      
167      
168 <    //find all free atoms (which do not belong to cutoff groups)  
169 <    //performs the "difference" operation from set theory,  the output range contains a copy of every
170 <    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
155 <    //[cutoffAtoms.begin(), cutoffAtoms.end()).
156 <    std::vector<Atom*> freeAtoms;    
157 <    std::set_difference(allAtoms.begin(), allAtoms.end(), cutoffAtoms.begin(), cutoffAtoms.end(),
158 <                            std::back_inserter(freeAtoms));
159 <
160 <    if (freeAtoms.size() != allAtoms.size() - cutoffAtoms.size()) {
161 <        //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
162 <        sprintf(painCave.errMsg, "Atoms in cutoff groups are not in the atom list of the same molecule");
163 <
164 <        painCave.isFatal = 1;
165 <        simError();        
166 <    }
167 <
168 <    //loop over the free atoms and then create one cutoff group for every single free atom
169 <    std::vector<Atom*>::iterator fai;
170 <
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);
173 <        mol->addCutoffGroup(cutoffGroup);
172 >      cutoffGroup = createCutoffGroup(mol, *fai);
173 >      mol->addCutoffGroup(cutoffGroup);
174      }
175      //create constraints
176      createConstraintPair(mol);
# Line 178 | Line 178 | Molecule* MoleculeCreator::createMolecule(ForceField*
178      
179      //the construction of this molecule is finished
180      mol->complete();
181 <
181 >    
182      return mol;
183 < }    
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());
195 >      sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
196 >              stamp->getType().c_str());
197  
198 <        painCave.isFatal = 1;
199 <        simError();
198 >      painCave.isFatal = 1;
199 >      simError();
200      }
201      
202      //below code still have some kind of hard-coding smell
203      if (atomType->isDirectional()){
204      
205 <        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
205 >      DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
206          
207 <        if (dAtomType == NULL) {
208 <            sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
207 >      if (dAtomType == NULL) {
208 >        sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
209  
210 <            painCave.isFatal = 1;
211 <            simError();
212 <        }
210 >        painCave.isFatal = 1;
211 >        simError();
212 >      }
213  
214 <        DirectionalAtom* dAtom;
215 <        dAtom = new DirectionalAtom(dAtomType);
216 <        atom = dAtom;    
214 >      DirectionalAtom* dAtom;
215 >      dAtom = new DirectionalAtom(dAtomType);
216 >      atom = dAtom;    
217      }
218      else{
219 <        atom = new Atom(atomType);
219 >      atom = new Atom(atomType);
220      }
221  
222      atom->setLocalIndex(localIndexMan->getNextAtomIndex());
223  
224      return atom;
225 < }
226 <
227 < RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
228 <                                                                                    RigidBodyStamp* rbStamp,
229 <                                                                                    LocalIndexManager* localIndexMan) {
225 >  }
226 >  
227 >  RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp,
228 >                                              Molecule* mol,
229 >                                              RigidBodyStamp* rbStamp,
230 >                                              LocalIndexManager* localIndexMan) {
231      Atom* atom;
232      int nAtoms;
233      Vector3d refCoor;
# Line 234 | Line 236 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
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
241 <        atom = mol->getAtomAt(rbStamp->getMember(i));
242 <        atomStamp= molStamp->getAtom(rbStamp->getMember(i));    
243 <        rb->addAtom(atom, atomStamp);
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 253 | Line 257 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
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 <    rb->setType(mol->getType() + "_RB_" + toString(mol->getNRigidBodies()));
261 <    
260 >    std::string s = OpenMD_itoa(mol->getNRigidBodies(), 10);
261 >    rb->setType(mol->getType() + "_RB_" + s.c_str());
262 >
263      return rb;
264 < }    
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());
278  
279      if (bondType == NULL) {
280 <        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
281 <                   atomA->getType().c_str(),
282 <                   atomB->getType().c_str());
283 <
284 <        painCave.isFatal = 1;
285 <        simError();
280 >      sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
281 >              atomA->getType().c_str(),
282 >              atomB->getType().c_str());
283 >      
284 >      painCave.isFatal = 1;
285 >      simError();
286      }
287      return new Bond(atomA, atomB, bondType);    
288 < }    
288 >  }    
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 >      
299 >      assert( atomA && atomB && atomC);
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 >        
311 >        painCave.isFatal = 1;
312 >        simError();
313 >      }
314 >      
315 >      bend = new Bend(atomA, atomB, atomC, bendType);
316 >    } else if ( bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
317 >      int ghostIndex = stamp->getGhostVectorSource();
318 >      int normalIndex = ghostIndex != bendAtoms[0] ? bendAtoms[0] : bendAtoms[1];
319 >      Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
320 >      DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
321 >      if (ghostAtom == NULL) {
322 >        sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
323 >        painCave.isFatal = 1;
324 >        simError();
325 >      }
326 >                
327 >      BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
328  
329 < Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
330 <    bool isGhostBend = false;
331 <    int ghostIndex;
329 >      if (bendType == NULL) {
330 >        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
331 >                normalAtom->getType().c_str(),
332 >                ghostAtom->getType().c_str(),
333 >                "GHOST");
334  
335 +        painCave.isFatal = 1;
336 +        simError();
337 +      }
338 +      
339 +      bend = new GhostBend(normalAtom, ghostAtom, bendType);      
340 +      
341 +    }
342      
343 <    //
344 <    if (stamp->haveExtras()){
291 <        LinkedAssign* extras = stamp->getExtras();
292 <        LinkedAssign* currentExtra = extras;
343 >    return bend;
344 >  }    
345  
346 <        while (currentExtra != NULL){
347 <            if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
296 <                switch (currentExtra->getType()){
297 <                case 0:
298 <                    ghostIndex = currentExtra->getInt();
299 <                    isGhostBend = true;
300 <                    break;
346 >  Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
347 >                                          TorsionStamp* stamp) {
348  
349 <                default:
350 <                sprintf(painCave.errMsg,
351 <                "SimSetup Error: ghostVectorSource must be an int.\n");
352 <                painCave.isFatal = 1;
306 <                simError();
307 <                }
308 <            } else{
309 <                sprintf(painCave.errMsg,
310 <                "SimSetup Error: unhandled bend assignment:\n");
311 <                painCave.isFatal = 1;
312 <                simError();
313 <            }
314 <            currentExtra = currentExtra->getNext();
315 <        }
316 <        
349 >    Torsion* torsion = NULL;
350 >    std::vector<int> torsionAtoms = stamp->getMembers();
351 >    if (torsionAtoms.size() < 3) {
352 >        return torsion;
353      }
354  
355 <    if (isGhostBend) {
355 >    Atom* atomA = mol->getAtomAt(torsionAtoms[0]);
356 >    Atom* atomB = mol->getAtomAt(torsionAtoms[1]);
357 >    Atom* atomC = mol->getAtomAt(torsionAtoms[2]);
358  
359 <        int indexA = stamp->getA();
360 <        int indexB= stamp->getB();
359 >    if (torsionAtoms.size() == 4) {
360 >      Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
361  
362 <        assert(indexA != indexB);
325 <
326 <        int normalIndex;
327 <        if (indexA == ghostIndex) {
328 <            normalIndex = indexB;
329 <        } else if (indexB == ghostIndex) {
330 <            normalIndex = indexA;
331 <        }
362 >      assert(atomA && atomB && atomC && atomD);
363          
364 <        Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
365 <        DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
366 <        if (ghostAtom == NULL) {
367 <            sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
368 <            painCave.isFatal = 1;
369 <            simError();
370 <        }
371 <                
372 <        BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
373 <
343 <        if (bendType == NULL) {
344 <            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
345 <                       normalAtom->getType().c_str(),
346 <                       ghostAtom->getType().c_str(),
347 <                       "GHOST");
348 <
349 <            painCave.isFatal = 1;
350 <            simError();
351 <        }
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          
375 <        return new GhostBend(normalAtom, ghostAtom, bendType);      
376 <
377 <    } else {
378 <            
379 <        Atom* atomA = mol->getAtomAt(stamp->getA());
380 <        Atom* atomB = mol->getAtomAt(stamp->getB());
381 <        Atom* atomC = mol->getAtomAt(stamp->getC());
382 <
383 <        assert( atomA && atomB && atomC);
375 >        painCave.isFatal = 1;
376 >        simError();
377 >      }
378 >      
379 >      torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
380 >    }
381 >    else {
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 >      
390 >      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
391 >                                                    atomC->getType(), "GHOST");
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          
400 <        BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
401 <
402 <        if (bendType == NULL) {
403 <            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
404 <                       atomA->getType().c_str(),
368 <                       atomB->getType().c_str(),
369 <                       atomC->getType().c_str());
370 <
371 <            painCave.isFatal = 1;
372 <            simError();
373 <        }
374 <
375 <        return new Bend(atomA, atomB, atomC, bendType);      
400 >        painCave.isFatal = 1;
401 >        simError();
402 >      }
403 >      
404 >      torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
405      }
406 < }    
406 >    
407 >    return torsion;
408 >  }    
409  
410 < Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
411 <
412 <    Atom* atomA = mol->getAtomAt(stamp->getA());
413 <    Atom* atomB = mol->getAtomAt(stamp->getB());
414 <    Atom* atomC = mol->getAtomAt(stamp->getC());
415 <    Torsion* torsion;
416 <
417 <    if (stamp->getD() != -1) {
387 <        Atom* atomD = mol->getAtomAt(stamp->getD());
388 <
389 <        assert(atomA && atomB && atomC && atomD);
390 <        
391 <        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
392 <                                                           atomC->getType(), atomD->getType());
393 <
394 <        if (torsionType == NULL) {
395 <            sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
396 <                       atomA->getType().c_str(),
397 <                       atomB->getType().c_str(),
398 <                       atomC->getType().c_str(),
399 <                       atomD->getType().c_str());
400 <
401 <            painCave.isFatal = 1;
402 <            simError();
403 <        }
404 <        
405 <        torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
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      }
407 –    else {
419  
420 <        DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(atomC);
421 <        if (dAtom == NULL) {
422 <            sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
423 <            painCave.isFatal = 1;
424 <            simError();
425 <        }        
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 <        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
433 <                                                           atomC->getType(), "GHOST");
434 <
435 <        if (torsionType == NULL) {
436 <            sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
437 <                       atomA->getType().c_str(),
438 <                       atomB->getType().c_str(),
439 <                       atomC->getType().c_str(),
440 <                       "GHOST");
441 <
442 <            painCave.isFatal = 1;
443 <            simError();
444 <        }
445 <        
446 <        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
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 <    return torsion;
434 < }    
435 <
436 < CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
452 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
453      int nAtoms;
454      CutoffGroup* cg;
455      Atom* atom;
# Line 441 | Line 457 | CutoffGroup* MoleculeCreator::createCutoffGroup(Molecu
457      
458      nAtoms = stamp->getNMembers();
459      for (int i =0; i < nAtoms; ++i) {
460 <        atom = mol->getAtomAt(stamp->getMember(i));
461 <        assert(atom);
462 <        cg->addAtom(atom);
460 >      atom = mol->getAtomAt(stamp->getMemberAt(i));
461 >      assert(atom);
462 >      cg->addAtom(atom);
463      }
464  
465      return cg;
466 < }    
466 >  }    
467  
468 < CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
468 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
469      CutoffGroup* cg;
470      cg  = new CutoffGroup();
471      cg->addAtom(atom);
472      return cg;
473 < }
473 >  }
474  
475 < void MoleculeCreator::createConstraintPair(Molecule* mol) {
475 >  void MoleculeCreator::createConstraintPair(Molecule* mol) {
476  
477      //add bond constraints
478      Molecule::BondIterator bi;
479      Bond* bond;
480      for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
481          
482 <        BondType* bt = bond->getBondType();
482 >      BondType* bt = bond->getBondType();
483  
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?        
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  
491 <        if (typeid(FixedBondType) == typeid(*bt)) {
492 <            FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
491 >      if (typeid(FixedBondType) == typeid(*bt)) {
492 >        FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
493  
494 <            ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
495 <            ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
496 <            ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
497 <            mol->addConstraintPair(consPair);
498 <        }
494 >        ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
495 >        ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
496 >        ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
497 >        mol->addConstraintPair(consPair);
498 >      }
499      }
500  
501      //rigidbody -- rigidbody constraint is not support yet
502 < }
502 >  }
503  
504 < void MoleculeCreator::createConstraintElem(Molecule* mol) {
504 >  void MoleculeCreator::createConstraintElem(Molecule* mol) {
505  
506      ConstraintPair* consPair;
507      Molecule::ConstraintPairIterator cpi;
508      std::set<StuntDouble*> sdSet;
509      for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
510  
511 <        StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
512 <        if (sdSet.find(sdA) == sdSet.end()){
513 <            sdSet.insert(sdA);
514 <            mol->addConstraintElem(new ConstraintElem(sdA));
515 <        }
511 >      StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
512 >      if (sdSet.find(sdA) == sdSet.end()){
513 >        sdSet.insert(sdA);
514 >        mol->addConstraintElem(new ConstraintElem(sdA));
515 >      }
516  
517 <        StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
518 <        if (sdSet.find(sdB) == sdSet.end()){
519 <            sdSet.insert(sdB);
520 <            mol->addConstraintElem(new ConstraintElem(sdB));
521 <        }
517 >      StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
518 >      if (sdSet.find(sdB) == sdSet.end()){
519 >        sdSet.insert(sdB);
520 >        mol->addConstraintElem(new ConstraintElem(sdB));
521 >      }
522          
523      }
524  
525 < }
525 >  }
526      
527   }

Diff Legend

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