ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/brains/MoleculeCreator.cpp
Revision: 1734
Committed: Fri Nov 12 07:05:43 2004 UTC (19 years, 8 months ago) by tim
File size: 11282 byte(s)
Log Message:
MoleculeCreator forget to create CutoffGroups for free atoms

File Contents

# Content
1 /*
2 * Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3 *
4 * Contact: oopse@oopse.org
5 *
6 * This program is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public License
8 * as published by the Free Software Foundation; either version 2.1
9 * of the License, or (at your option) any later version.
10 * All we ask is that proper credit is given for our work, which includes
11 * - but is not limited to - adding the above copyright notice to the beginning
12 * of your source code files, and to any copyright notice that you may distribute
13 * with programs based on this work.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU Lesser General Public License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public License
21 * along with this program; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23 *
24 */
25
26 /**
27 * @file MoleculeCreator.cpp
28 * @author tlin
29 * @date 11/04/2004
30 * @time 13:44am
31 * @version 1.0
32 */
33
34 #include "brains/MoleculeCreator.hpp"
35 #include <cassert>
36
37 namespace oopse {
38
39 Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp, int stampId, int globalIndex) {
40 Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getID());
41
42 //create atoms
43 Atom* atom;
44 AtomStamp* currentAtomStamp;
45 int nAtom = molStamp->getNAtoms();
46 for (int i = 0; i < nAtom; ++i) {
47 currentAtomStamp = molStamp->getAtom(i);
48 atom = createAtom(ff, currentAtomStamp);
49 mol->addAtom(atom);
50 }
51
52 //create rigidbodies
53 RigidBody* rb;
54 RigidBodyStamp * currentRigidBodyStamp;
55 int nRigidbodies = molStamp->getNRigidBodies();
56
57 for (int i = 0; i < nRigidbodies; ++i) {
58 currentRigidBodyStamp = molStamp->getRigidBody(i);
59 rb = createRigidBody(ff, mol, currentRigidBodyStamp);
60 mol->addRigidBody(rb);
61 }
62
63 //create bonds
64 Bond* bond;
65 BondStamp* currentBondStamp;
66 int nBonds = molStamp->getNBends();
67
68 for (int i = 0; i < nBonds; ++i) {
69 currentBondStamp = molStamp->getBend(i);
70 bond = createBond(ff, mol, currentBondStamp);
71 mol->addBond(bond);
72 }
73
74 //create bends
75 Bend* bend;
76 BendStamp* currentBendStamp;
77 int nBends = molStamp->getNBends();
78 for (int i = 0; i < nBends; ++i) {
79 currentBendStamp = molStamp->getBend(i);
80 bend = createBend(ff, mol, currentBendStamp);
81 mol->addBend(bend);
82 }
83
84 //create torsions
85 Torsion* torsion;
86 TorsionStamp* currentTorsionStamp;
87 int nTorsions = molStamp->getNTorsions();
88 for (int i = 0; i < nTorsions; ++i) {
89 currentTorsionStamp = molStamp->getTorsion(i);
90 torsion = createTorsion(ff, mol, currentTorsionStamp);
91 mol->addTorsion(torsion);
92 }
93
94 //create cutoffGroups
95 CutoffGroup* cutoffGroup;
96 CutoffGroupStamp* currentCutoffGroupStamp;
97 int nCutoffGroups = molStamp->getNCutoffGroups();
98 for (int i = 0; i < nCutoffGroups; ++i) {
99 currentCutoffGroupStamp = molStamp->getCutoffGroup(i);
100 cutoffGroup = createCutoffGroup(ff, mol, currentCutoffGroupStamp);
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
152 mol->complete();
153
154 return mol;
155 }
156
157
158 Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp,
159 LocalIndexManager* localIndexMan) {
160 AtomType * atomType;
161 Atom* atom;
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;
176 double phi;
177 double theta;
178 double psi;
179
180 dAtom = new DirectionalAtom(atomType);
181
182 // Directional Atoms have standard unit vectors which are oriented
183 // in space using the three Euler angles. We assume the standard
184 // unit vector was originally along the z axis below.
185
186 phi = stamp->getEulerPhi() * M_PI / 180.0;
187 theta = stamp->getEulerTheta() * M_PI / 180.0;
188 psi = stamp->getEulerPsi()* M_PI / 180.0;
189
190 dAtom->setUnitFrameFromEuler(phi, theta, psi);
191 atom = dAtom;
192 }
193 else{
194 atom = new Atom(atomType);
195 }
196
197 atom->setLocalIndex(localIndexMan->getNextAtomIndex());
198 }
199
200 RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
201 RigidBodyStamp* rbStamp,
202 LocalIndexManager* localIndexMan) {
203 Atom* atom;
204 int nAtoms;
205 Vector3d refCoor;
206 AtomStamp* atomStamp;
207
208 nAtoms = molStamp->getNMembers();
209
210 RigidBody* rb = new RigidBody();
211
212 for (int i = 0; i < nAtoms; ++i) {
213 //rbStamp->getMember(i) return the local index of current atom inside the molecule.
214 //It is not the same as local index of atom which is the index of atom at DataStorage class
215 atom = mol->getAtomAt(rbStamp->getMember(i));
216 atomStamp= molStamp->molStamp->getAtom(rbStamp->getMember(i));
217 rb->addAtom(atom, atomStamp);
218 }
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
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 }
235
236 Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
237 BondType* bondType;
238 Atom* atomA;
239 Atom* atomB;
240
241 atomA = mol->getAtomAt(stamp->getA());
242 atomB = mol->getAtomAt(stamp->getB());
243
244 assert( atomA && atomB);
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* bendType;
261 Atom* atomA;
262 Atom* atomB;
263 Atom* atomC;
264
265 //need to consider the ghost bend
266 atomA = mol->getAtomAt(stamp->getA());
267 atomB = mol->getAtomAt(stamp->getB());
268 atomC = mol->getAtomAt(stamp->getC());
269
270 assert( atomA && atomB && atomC);
271
272 bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
273
274 if (bendType == NULL) {
275 sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
276 atomA->getType().c_str(),
277 atomB->getType().c_str(),
278 atomC->getType().c_str());
279
280 painCave.isFatal = 1;
281 simError();
282 }
283
284 return new Bond(atomA, atomB, bendType);
285
286 }
287
288 Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
289 TorsionType* torsionType;
290 Atom* atomA;
291 Atom* atomB;
292 Atom* atomC;
293 Atom* atomD;
294
295 atomA = mol->getAtomAt(stamp->getA());
296 atomB = mol->getAtomAt(stamp->getB());
297 atomC = mol->getAtomAt(stamp->getC());
298 atomD = mol->getAtomAt(stamp->getD());
299
300 assert(atomA && atomB && atomC && atomD);
301
302 torsionType = ff->getTosionType(atomA->getType(), atomB->getType(),
303 atomC->getType(), atomD->getType());
304
305 if (torsionType == NULL) {
306 sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
307 atomA->getType().c_str(),
308 atomB->getType().c_str(),
309 atomC->getType().c_str(),
310 atomD->getType().c_str());
311
312 painCave.isFatal = 1;
313 simError();
314 }
315
316 return new Torsion(atomA, atomB, torsionType);
317 }
318
319 CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
320 int nAtoms;
321 CutoffGroup* cg;
322 Atom* atom;
323 cg = new CutoffGroup();
324
325 nAtoms = stamp->getNMembers();
326 for (int i =0; i < nAtoms; ++i) {
327 atom = mol->getAtomAt(stamp->getMember(i));
328 assert(atom);
329 cg->addAtom(atom);
330 }
331
332 return cg;
333 }
334
335 CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
336 CutoffGroup* cg;
337 cg = new CutoffGroup();
338 cg->addAtom();
339 return cg;
340 }
341 //Constraint* MoleculeCreator::createConstraint() {
342
343 //}
344
345 }

Properties

Name Value
svn:executable *