ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/brains/MoleculeCreator.cpp
Revision: 3432
Committed: Mon Jul 14 12:35:58 2008 UTC (16 years, 9 months ago) by gezelter
File size: 17524 byte(s)
Log Message:
Changes for implementing Amber force field:  Added Inversions and
worked on BaseAtomTypes so that they'd function with the fortran side.

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
369 if (torsionType == NULL) {
370 sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
371 atomA->getType().c_str(),
372 atomB->getType().c_str(),
373 atomC->getType().c_str(),
374 atomD->getType().c_str());
375
376 painCave.isFatal = 1;
377 simError();
378 }
379
380 torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);
381 }
382 else {
383
384 DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(stamp->getGhostVectorSource()));
385 if (dAtom == NULL) {
386 sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
387 painCave.isFatal = 1;
388 simError();
389 }
390
391 TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
392 atomC->getType(), "GHOST");
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 "GHOST");
400
401 painCave.isFatal = 1;
402 simError();
403 }
404
405 torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);
406 }
407
408 return torsion;
409 }
410
411 Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol,
412 InversionStamp* stamp) {
413
414 Inversion* inversion = NULL;
415 int center = stamp->getCenter();
416 std::vector<int> satellites = stamp->getSatellites();
417 if (satellites.size() != 3) {
418 return inversion;
419 }
420
421 Atom* atomA = mol->getAtomAt(center);
422 Atom* atomB = mol->getAtomAt(satellites[0]);
423 Atom* atomC = mol->getAtomAt(satellites[1]);
424 Atom* atomD = mol->getAtomAt(satellites[2]);
425
426 assert(atomA && atomB && atomC && atomD);
427
428 InversionType* inversionType = ff->getInversionType(atomA->getType(),
429 atomB->getType(),
430 atomC->getType(),
431 atomD->getType());
432
433 if (inversionType == NULL) {
434 sprintf(painCave.errMsg, "No Matching Inversion Type for[%s, %s, %s, %s]\n"
435 "\t(May not be a problem: not all inversions are parametrized)\n",
436 atomA->getType().c_str(),
437 atomB->getType().c_str(),
438 atomC->getType().c_str(),
439 atomD->getType().c_str());
440
441 painCave.isFatal = 0;
442 painCave.severity = OOPSE_INFO;
443 simError();
444 return NULL;
445 } else {
446
447 inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType);
448 return inversion;
449 }
450 }
451
452
453 CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
454 int nAtoms;
455 CutoffGroup* cg;
456 Atom* atom;
457 cg = new CutoffGroup();
458
459 nAtoms = stamp->getNMembers();
460 for (int i =0; i < nAtoms; ++i) {
461 atom = mol->getAtomAt(stamp->getMemberAt(i));
462 assert(atom);
463 cg->addAtom(atom);
464 }
465
466 return cg;
467 }
468
469 CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
470 CutoffGroup* cg;
471 cg = new CutoffGroup();
472 cg->addAtom(atom);
473 return cg;
474 }
475
476 void MoleculeCreator::createConstraintPair(Molecule* mol) {
477
478 //add bond constraints
479 Molecule::BondIterator bi;
480 Bond* bond;
481 for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
482
483 BondType* bt = bond->getBondType();
484
485 //class Parent1 {};
486 //class Child1 : public Parent {};
487 //class Child2 : public Parent {};
488 //Child1* ch1 = new Child1();
489 //Child2* ch2 = dynamic_cast<Child2*>(ch1);
490 //the dynamic_cast is succeed in above line. A compiler bug?
491
492 if (typeid(FixedBondType) == typeid(*bt)) {
493 FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
494
495 ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
496 ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());
497 ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
498 mol->addConstraintPair(consPair);
499 }
500 }
501
502 //rigidbody -- rigidbody constraint is not support yet
503 }
504
505 void MoleculeCreator::createConstraintElem(Molecule* mol) {
506
507 ConstraintPair* consPair;
508 Molecule::ConstraintPairIterator cpi;
509 std::set<StuntDouble*> sdSet;
510 for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
511
512 StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();
513 if (sdSet.find(sdA) == sdSet.end()){
514 sdSet.insert(sdA);
515 mol->addConstraintElem(new ConstraintElem(sdA));
516 }
517
518 StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();
519 if (sdSet.find(sdB) == sdSet.end()){
520 sdSet.insert(sdB);
521 mol->addConstraintElem(new ConstraintElem(sdB));
522 }
523
524 }
525
526 }
527
528 }

Properties

Name Value
svn:executable *