ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/brains/MoleculeCreator.cpp
Revision: 3446
Committed: Wed Sep 10 19:51:45 2008 UTC (16 years, 7 months ago) by cli2
File size: 17523 byte(s)
Log Message:
Inversion fixes and amber mostly working

File Contents

# Content
1 /*
2 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
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
19 * notice, this list of conditions and the following disclaimer.
20 *
21 * 3. Redistributions in binary form must reproduce the above copyright
22 * notice, this list of conditions and the following disclaimer in the
23 * documentation and/or other materials provided with the
24 * distribution.
25 *
26 * This software is provided "AS IS," without a warranty of any
27 * kind. All express or implied conditions, representations and
28 * warranties, including any implied warranty of merchantability,
29 * fitness for a particular purpose or non-infringement, are hereby
30 * excluded. The University of Notre Dame and its licensors shall not
31 * be liable for any damages suffered by licensee as a result of
32 * using, modifying or distributing the software or its
33 * derivatives. In no event will the University of Notre Dame or its
34 * licensors be liable for any lost revenue, profit or data, or for
35 * direct, indirect, special, consequential, incidental or punitive
36 * damages, however caused and regardless of the theory of liability,
37 * arising out of the use of or inability to use software, even if the
38 * University of Notre Dame has been advised of the possibility of
39 * such damages.
40 */
41
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"
55 #include "primitives/GhostBend.hpp"
56 #include "primitives/GhostTorsion.hpp"
57 #include "types/DirectionalAtomType.hpp"
58 #include "types/FixedBondType.hpp"
59 #include "utils/simError.h"
60 #include "utils/StringUtils.hpp"
61
62 namespace oopse {
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->getAtomStamp(i);
76 atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
77 mol->addAtom(atom);
78 }
79
80 //create rigidbodies
81 RigidBody* rb;
82 RigidBodyStamp * currentRigidBodyStamp;
83 int nRigidbodies = molStamp->getNRigidBodies();
84
85 for (int i = 0; i < nRigidbodies; ++i) {
86 currentRigidBodyStamp = molStamp->getRigidBodyStamp(i);
87 rb = createRigidBody(molStamp, mol, currentRigidBodyStamp,
88 localIndexMan);
89 mol->addRigidBody(rb);
90 }
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->getBondStamp(i);
99 bond = createBond(ff, mol, currentBondStamp);
100 mol->addBond(bond);
101 }
102
103 //create bends
104 Bend* bend;
105 BendStamp* currentBendStamp;
106 int nBends = molStamp->getNBends();
107 for (int i = 0; i < nBends; ++i) {
108 currentBendStamp = molStamp->getBendStamp(i);
109 bend = createBend(ff, mol, currentBendStamp);
110 mol->addBend(bend);
111 }
112
113 //create torsions
114 Torsion* torsion;
115 TorsionStamp* currentTorsionStamp;
116 int nTorsions = molStamp->getNTorsions();
117 for (int i = 0; i < nTorsions; ++i) {
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->getCutoffGroupStamp(i);
141 cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
142 mol->addCutoffGroup(cutoffGroup);
143 }
144
145 //every free atom is a cutoff group
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(fai); atom != NULL; atom = mol->nextAtom(fai)) {
152 freeAtoms.push_back(atom);
153 }
154
155 Molecule::CutoffGroupIterator ci;
156 CutoffGroup* cg;
157
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 // 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);
174 }
175 //create constraints
176 createConstraintPair(mol);
177 createConstraintElem(mol);
178
179 //the construction of this molecule is finished
180 mol->complete();
181
182 return mol;
183 }
184
185
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
194 if (atomType == NULL) {
195 sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
196 stamp->getType().c_str());
197
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);
206
207 if (dAtomType == NULL) {
208 sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
209
210 painCave.isFatal = 1;
211 simError();
212 }
213
214 DirectionalAtom* dAtom;
215 dAtom = new DirectionalAtom(dAtomType);
216 atom = dAtom;
217 }
218 else{
219 atom = new Atom(atomType);
220 }
221
222 atom->setLocalIndex(localIndexMan->getNextAtomIndex());
223
224 return atom;
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;
234 AtomStamp* atomStamp;
235
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
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
248 //reference coordinates
249 rb->calcRefCoords();
250
251 //set the local index of this rigid body, global index will be set later
252 rb->setLocalIndex(localIndexMan->getNextRigidBodyIndex());
253
254 //the rule for naming rigidbody MoleculeName_RB_Integer
255 //The first part is the name of the molecule
256 //The second part is alway fixed as "RB"
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);
261 rb->setType(mol->getType() + "_RB_" + s.c_str());
262
263 return rb;
264 }
265
266 Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
267 BondStamp* stamp) {
268 BondType* bondType;
269 Atom* atomA;
270 Atom* atomB;
271
272 atomA = mol->getAtomAt(stamp->getA());
273 atomB = mol->getAtomAt(stamp->getB());
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();
286 }
287 return new Bond(atomA, atomB, bondType);
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 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 return bend;
344 }
345
346 Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
347 TorsionStamp* stamp) {
348
349 Torsion* torsion = NULL;
350 std::vector<int> torsionAtoms = stamp->getMembers();
351 if (torsionAtoms.size() < 3) {
352 return torsion;
353 }
354
355 Atom* atomA = mol->getAtomAt(torsionAtoms[0]);
356 Atom* atomB = mol->getAtomAt(torsionAtoms[1]);
357 Atom* atomC = mol->getAtomAt(torsionAtoms[2]);
358
359 if (torsionAtoms.size() == 4) {
360 Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
361
362 assert(atomA && atomB && atomC && atomD);
363
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 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 painCave.isFatal = 1;
401 simError();
402 }
403
404 torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);
405 }
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 = OOPSE_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;
455 Atom* atom;
456 cg = new CutoffGroup();
457
458 nAtoms = stamp->getNMembers();
459 for (int i =0; i < nAtoms; ++i) {
460 atom = mol->getAtomAt(stamp->getMemberAt(i));
461 assert(atom);
462 cg->addAtom(atom);
463 }
464
465 return cg;
466 }
467
468 CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
469 CutoffGroup* cg;
470 cg = new CutoffGroup();
471 cg->addAtom(atom);
472 return cg;
473 }
474
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();
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?
490
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 }
499 }
500
501 //rigidbody -- rigidbody constraint is not support yet
502 }
503
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 }
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 }
522
523 }
524
525 }
526
527 }

Properties

Name Value
svn:executable *