ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/brains/MoleculeCreator.cpp

Comparing trunk/src/brains/MoleculeCreator.cpp (file contents):
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 1442 by gezelter, Mon May 10 17:28:26 2010 UTC

#	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		/**
#	Line 48 | Line 48
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);
75	>	currentAtomStamp = molStamp->getAtomStamp(i);
76		atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
77		mol->addAtom(atom);
78		}
#	Line 81 | Line 83 | namespace oopse {
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);
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);
98	>	currentBondStamp = molStamp->getBondStamp(i);
99		bond = createBond(ff, mol, currentBondStamp);
100		mol->addBond(bond);
101		}
#	Line 102 | Line 105 | namespace oopse {
105		BendStamp* currentBendStamp;
106		int nBends = molStamp->getNBends();
107		for (int i = 0; i < nBends; ++i) {
108	<	currentBendStamp = molStamp->getBend(i);
108	>	currentBendStamp = molStamp->getBendStamp(i);
109		bend = createBend(ff, mol, currentBendStamp);
110		mol->addBend(bend);
111		}
#	Line 112 | Line 115 | namespace oopse {
115		TorsionStamp* currentTorsionStamp;
116		int nTorsions = molStamp->getNTorsions();
117		for (int i = 0; i < nTorsions; ++i) {
118	<	currentTorsionStamp = molStamp->getTorsion(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->getCutoffGroup(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::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	<
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	<	cutoffAtoms.insert(atom);
163	<	}
164	<
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);
#	Line 178 | Line 178 | namespace oopse {
178
179		//the construction of this molecule is finished
180		mol->complete();
181	<
181	>
182		return mol;
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());
196	>	stamp->getType().c_str());
197
198		painCave.isFatal = 1;
199		simError();
#	Line 222 | Line 223 | namespace oopse {
223
224		return atom;
225		}
226	<
227	<	RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
226	>
227	>	RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp,
228	>	Molecule* mol,
229		RigidBodyStamp* rbStamp,
230		LocalIndexManager* localIndexMan) {
231		Atom* atom;
#	Line 234 | Line 236 | namespace oopse {
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));
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 | namespace oopse {
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);
260	>	std::string s = OpenMD_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, 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());
#	Line 275 | Line 280 | namespace oopse {
280		sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
281		atomA->getType().c_str(),
282		atomB->getType().c_str());
283	<
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, BendStamp* stamp) {
291	<	bool isGhostBend = false;
292	<	int ghostIndex;
293	<
294	<
295	<	//
296	<	if (stamp->haveExtras()){
297	<	LinkedAssign* extras = stamp->getExtras();
298	<	LinkedAssign* currentExtra = extras;
299	<
300	<	while (currentExtra != NULL){
301	<	if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
302	<	switch (currentExtra->getType()){
303	<	case 0:
304	<	ghostIndex = currentExtra->getInt();
305	<	isGhostBend = true;
306	<	break;
307	<
308	<	default:
309	<	sprintf(painCave.errMsg,
305	<	"SimSetup Error: ghostVectorSource must be an int.\n");
306	<	painCave.isFatal = 1;
307	<	simError();
308	<	}
309	<	} else{
310	<	sprintf(painCave.errMsg,
311	<	"SimSetup Error: unhandled bend assignment:\n");
312	<	painCave.isFatal = 1;
313	<	simError();
314	<	}
315	<	currentExtra = currentExtra->getNext();
316	<	}
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	<	}
312	<
320	<	if (isGhostBend) {
321	<
322	<	int indexA = stamp->getA();
323	<	int indexB= stamp->getB();
324	<
325	<	assert(indexA != indexB);
326	<
327	<	int normalIndex;
328	<	if (indexA == ghostIndex) {
329	<	normalIndex = indexB;
330	<	} else if (indexB == ghostIndex) {
331	<	normalIndex = indexA;
311	>	painCave.isFatal = 1;
312	>	simError();
313		}
314	<
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) {
#	Line 350 | Line 335 | namespace oopse {
335		painCave.isFatal = 1;
336		simError();
337		}
338	<
339	<	return new GhostBend(normalAtom, ghostAtom, bendType);
338	>
339	>	bend = new GhostBend(normalAtom, ghostAtom, bendType);
340	>
341	>	}
342	>
343	>	return bend;
344	>	}
345
346	<	} else {
347	<
358	<	Atom* atomA = mol->getAtomAt(stamp->getA());
359	<	Atom* atomB = mol->getAtomAt(stamp->getB());
360	<	Atom* atomC = mol->getAtomAt(stamp->getC());
346	>	Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
347	>	TorsionStamp* stamp) {
348
349	<	assert( atomA && atomB && atomC);
350	<
351	<	BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
352	<
366	<	if (bendType == NULL) {
367	<	sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
368	<	atomA->getType().c_str(),
369	<	atomB->getType().c_str(),
370	<	atomC->getType().c_str());
371	<
372	<	painCave.isFatal = 1;
373	<	simError();
374	<	}
375	<
376	<	return new Bend(atomA, atomB, atomC, bendType);
349	>	Torsion* torsion = NULL;
350	>	std::vector<int> torsionAtoms = stamp->getMembers();
351	>	if (torsionAtoms.size() < 3) {
352	>	return torsion;
353		}
378	–	}
354
355	<	Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
355	>	Atom* atomA = mol->getAtomAt(torsionAtoms[0]);
356	>	Atom* atomB = mol->getAtomAt(torsionAtoms[1]);
357	>	Atom* atomC = mol->getAtomAt(torsionAtoms[2]);
358
359	<	Atom* atomA = mol->getAtomAt(stamp->getA());
360	<	Atom* atomB = mol->getAtomAt(stamp->getB());
384	<	Atom* atomC = mol->getAtomAt(stamp->getC());
385	<	Torsion* torsion;
359	>	if (torsionAtoms.size() == 4) {
360	>	Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
361
387	–	if (stamp->getD() != -1) {
388	–	Atom* atomD = mol->getAtomAt(stamp->getD());
389	–
362		assert(atomA && atomB && atomC && atomD);
363
364	<	TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
365	<	atomC->getType(), atomD->getType());
366	<
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	<
374	>
375		painCave.isFatal = 1;
376		simError();
377		}
378	<
378	>
379		torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);
380		}
381		else {
382	<
383	<	DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(atomC);
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	<
389	>
390		TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
391		atomC->getType(), "GHOST");
392	<
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	<
399	>
400		painCave.isFatal = 1;
401		simError();
402		}
403	<
403	>
404		torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);
405		}
406	<
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 = 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		CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
453		int nAtoms;
454		CutoffGroup* cg;
#	Line 442 | Line 457 | namespace oopse {
457
458		nAtoms = stamp->getNMembers();
459		for (int i =0; i < nAtoms; ++i) {
460	<	atom = mol->getAtomAt(stamp->getMember(i));
460	>	atom = mol->getAtomAt(stamp->getMemberAt(i));
461		assert(atom);
462		cg->addAtom(atom);
463		}

Comparing trunk/src/brains/MoleculeCreator.cpp (property svn:keywords):
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 1442 by gezelter, Mon May 10 17:28:26 2010 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

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