ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/brains/MoleculeCreator.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/brains/MoleculeCreator.cpp (file contents):
Revision 1719 by tim, Fri Nov 5 23:38:27 2004 UTC vs.
Revision 1774 by tim, Tue Nov 23 23:12:23 2004 UTC

# Line 37 | Line 37 | Molecule* MoleculeCreator::createMolecule(ForceField*
37   namespace oopse {
38  
39   Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp, int stampId, int globalIndex) {
40 <    Molecule* mol = new Molecule(stampId, globalIndex);
40 >    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getID());
41      
42      //create atoms
43      Atom* atom;
# Line 101 | Line 101 | Molecule* MoleculeCreator::createMolecule(ForceField*
101          mol->addCutoffGroup(cutoffGroup);
102      }
103  
104 +    //every free atom is a cutoff group    
105 +    std::set<Atom*> allAtoms;
106 +    typename  Molecule::AtomIterator ai;
107 +
108 +    //add all atoms into allAtoms set
109 +    for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
110 +        allAtoms.insert(atom);
111 +    }
112 +
113 +    typename Molecule::CutoffGroupIterator ci;
114 +    CutoffGroup* cg;
115 +    std::set<Atom*> cutoffAtoms;    
116 +    
117 +    //add all of the atoms belong to cutoff groups into cutoffAtoms set
118 +    for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
119 +
120 +        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
121 +            cutoffAtoms.insert(atom);
122 +        }
123 +
124 +    }      
125 +    
126 +    //find all free atoms (which do not belong to cutoff groups)  
127 +    //performs the "difference" operation from set theory,  the output range contains a copy of every
128 +    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
129 +    //[cutoffAtoms.begin(), cutoffAtoms.end()).
130 +    std::vector<Atom*> freeAtoms;    
131 +    std::set_difference(allAtoms.begin(), allAtoms.end(), cutoffAtoms.begin(), cutoffAtoms.end(),
132 +                            std::back_inserter(freeAtoms.end()));
133 +
134 +    if (freeAtoms.size() != allAtoms.size() - cutoffAtoms.size()) {
135 +        //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
136 +        sprintf(painCave.errMsg, "Atoms in cutoff groups are not in the atom list of the same molecule");
137 +
138 +        painCave.isFatal = 1;
139 +        simError();        
140 +    }
141 +
142 +    //loop over the free atoms and then create one cutoff group for every single free atom
143 +    std::vector<Atom*>::iterator fai;
144 +
145 +    for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
146 +        createCutoffGroup(mol, *fai);
147 +        mol->addCutoffGroup(cutoffGroup);
148 +    }
149      //create constraints
150  
151      //the construction of this molecule is finished
# Line 117 | Line 162 | Atom* MoleculeCreator::createAtom(ForceField* ff, Mole
162  
163      atomType =  ff->getAtomType(stamp->getType());
164  
165 +    if (bondType == NULL) {
166 +        sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
167 +                   stamp->getType());
168 +
169 +        painCave.isFatal = 1;
170 +        simError();
171 +    }
172 +    
173      //below code still have some kind of hard-coding smell
174      if (stamp->haveOrientation()){
175          DirectionalAtom* dAtom;
# Line 152 | Line 205 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
205      Vector3d refCoor;
206      AtomStamp* atomStamp;
207      
208 <    nAtoms = stamp->getNMembers();
208 >    nAtoms = molStamp->getNMembers();
209  
210      RigidBody* rb = new RigidBody();
211      
# Line 163 | Line 216 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
216          atomStamp= molStamp->molStamp->getAtom(rbStamp->getMember(i));    
217          rb->addAtom(atom, atomStamp);
218      }
219 <    
219 >
220 >    //after all of the atoms are added, we need to calculate the reference coordinates
221      rb->calcRefCoords();
222 +
223 +    //set the local index of this rigid body, global index will be set later
224      rb->setLocalIndex(localIndexMan->getNextRigidBodyIndex());
225 <    rb->setType(mol->getType() + "_" + itoa(mol.getNRigidBodies()));
225 >
226 >    //the rule for naming rigidbody MoleculeName_RB_Integer
227 >    //The first part is the name of the molecule
228 >    //The second part is alway fixed as "RB"
229 >    //The third part is the index of the rigidbody defined in meta-data file
230 >    //For example, Butane_RB_0 is a valid rigid body name of butane molecule
231 >    rb->setType(mol->getType() + "_RB_" + itoa(mol.getNRigidBodies()));
232      
233      return rb;
234   }    
# Line 183 | Line 245 | Bond* MoleculeCreator::createBond(ForceField* ff, Mole
245      
246      bondType = ff->getBondType(atomA, atomB);
247  
248 +    if (bondType == NULL) {
249 +        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
250 +                   atomA->getType().c_str(),
251 +                   atomB->getType().c_str());
252 +
253 +        painCave.isFatal = 1;
254 +        simError();
255 +    }
256      return new Bond(atomA, atomB, bondType);    
257   }    
258  
259   Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
260 <    BendType* benType;
261 <    Atom* atomA;
192 <    Atom* atomB;
193 <    Atom* atomC;
260 >    bool isGhostBend = false;
261 >    int ghostIndex;
262  
195 –    //need to consider the ghost bend
196 –    atomA = mol->getAtomAt(stamp->getA());
197 –    atomB = mol->getAtomAt(stamp->getB());
198 –    atomC = mol->getAtomAt(stamp->getC());
199 –
200 –    assert( atomA && atomB && atomC);
263      
264 <    benType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
265 <
266 <    return new Bond(atomA, atomB, benType);      
264 >    //
265 >    if (stamp->haveExtras()){
266 >        LinkedAssign* extras = currentBend->getExtras();
267 >        LinkedAssign* currentExtra = extras;
268  
269 +        while (currentExtra != NULL){
270 +            if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
271 +                switch (currentExtra->getType()){
272 +                case 0:
273 +                    ghostIndex = currentExtra->getInt();
274 +                    isGhostBend = true;
275 +                    break;
276 +
277 +                default:
278 +                sprintf(painCave.errMsg,
279 +                "SimSetup Error: ghostVectorSource must be an int.\n");
280 +                painCave.isFatal = 1;
281 +                simError();
282 +                }
283 +            } else{
284 +                sprintf(painCave.errMsg,
285 +                "SimSetup Error: unhandled bend assignment:\n");
286 +                painCave.isFatal = 1;
287 +                simError();
288 +            }
289 +            currentExtra = currentExtra->getNext();
290 +        }
291 +        
292 +    }
293 +
294 +    if (isGhostBend) {
295 +
296 +        int indexA = stamp->getA();
297 +        int indexB= stamp->getB();
298 +
299 +        assert(indexA != indexB);
300 +
301 +        int normalIndex;
302 +        if (indexA == ghostIndex) {
303 +            normalIndex = indexB;
304 +        } else if (indexB == ghostIndex) {
305 +            normalIndex = indexA;
306 +        }
307 +        
308 +        Atom* normalAtom = mol->getAtomAt(normalIndex) ;
309 +        Atom* ghostAtom = mol->getAtomAt(ghostIndex);
310 +        
311 +        BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
312 +
313 +        if (bendType == NULL) {
314 +            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
315 +                       normalAtom->getType().c_str(),
316 +                       ghostAtom->getType().c_str(),
317 +                       "GHOST");
318 +
319 +            painCave.isFatal = 1;
320 +            simError();
321 +        }
322 +
323 +        return new GhostBend(normalAtom, ghostAtom, bendType);      
324 +
325 +    } else {
326 +            
327 +        Atom* atomA = mol->getAtomAt(stamp->getA());
328 +        Atom* atomB = mol->getAtomAt(stamp->getB());
329 +        Atom* atomC = mol->getAtomAt(stamp->getC());
330 +
331 +        assert( atomA && atomB && atomC);
332 +        
333 +        BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
334 +
335 +        if (bendType == NULL) {
336 +            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
337 +                       atomA->getType().c_str(),
338 +                       atomB->getType().c_str(),
339 +                       atomC->getType().c_str());
340 +
341 +            painCave.isFatal = 1;
342 +            simError();
343 +        }
344 +
345 +        return new Bend(atomA, atomB, atomC, bendType);      
346 +    }
347   }    
348  
349   Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
# Line 222 | Line 363 | Torsion* MoleculeCreator::createTorsion(ForceField* ff
363      torsionType = ff->getTosionType(atomA->getType(), atomB->getType(),
364                                                         atomC->getType(), atomD->getType());
365  
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 +
373 +        painCave.isFatal = 1;
374 +        simError();
375 +    }
376 +    
377      return new Torsion(atomA, atomB, torsionType);      
378   }    
379  
# Line 241 | Line 393 | CutoffGroup* MoleculeCreator::createCutoffGroup(Molecu
393      return cg;
394   }    
395  
396 + CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
397 +    CutoffGroup* cg;
398 +    cg  = new CutoffGroup();
399 +    cg->addAtom();
400 +    return cg;
401 + }
402   //Constraint* MoleculeCreator::createConstraint() {
403  
404   //}

Diff Legend

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