ViewVC Help

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

root/OpenMD/trunk/src/utils/MoLocator.cpp

Comparing trunk/src/utils/MoLocator.cpp (file contents):
Revision 490 by tim, Wed Apr 13 18:41:17 2005 UTC vs.
Revision 1442 by gezelter, Mon May 10 17:28:26 2010 UTC

#	Line 1 | Line 1
1	<	/*
1	>	/*
2		* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3		*
4		* The University of Notre Dame grants you ("Licensee") a
#	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		#include <iostream>
#	Line 47 | Line 47
47		#include "utils/simError.h"
48		#include "utils/MoLocator.hpp"
49		#include "types/AtomType.hpp"
50	–	namespace oopse {
51	–	MoLocator::MoLocator( MoleculeStamp* theStamp, ForceField* theFF){
50
51	<	myStamp = theStamp;
52	<	myFF = theFF;
53	<	nIntegrableObjects = myStamp->getNIntegrable();
54	<	calcRef();
55	<	}
51	>	namespace OpenMD {
52	>	MoLocator::MoLocator( MoleculeStamp* theStamp, ForceField* theFF){
53	>
54	>	myStamp = theStamp;
55	>	myFF = theFF;
56	>	nIntegrableObjects = myStamp->getNIntegrable();
57	>	calcRef();
58	>	}
59	>
60	>	void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
61
59	–	void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
62		Vector3d newCoor;
63		Vector3d curRefCoor;
64		RotMat3x3d rotMat = latVec2RotMat(ort);
65	<
65	>
66		if(mol->getNIntegrableObjects() != nIntegrableObjects){
67	<	sprintf( painCave.errMsg,
68	<	"MoLocator error.\n"
69	<	"  The number of integrable objects of MoleculeStamp is not the same as  that of Molecule\n");
70	<	painCave.isFatal = 1;
71	<	simError();
67	>	sprintf( painCave.errMsg,
68	>	"MoLocator error.\n"
69	>	"  The number of integrable objects of MoleculeStamp is not the same as  that of Molecule\n");
70	>	painCave.isFatal = 1;
71	>	simError();
72		}
73	<
73	>
74		Molecule::IntegrableObjectIterator ii;
75		StuntDouble* integrableObject;
76		int i;
77		for (integrableObject = mol->beginIntegrableObject(ii), i = 0; integrableObject != NULL;
78	<	integrableObject = mol->nextIntegrableObject(ii), ++i) {
79	<
80	<	newCoor = rotMat * refCoords[i];
81	<	newCoor += offset;
82	<
83	<	integrableObject->setPos( newCoor);
84	<	integrableObject->setVel(V3Zero);
85	<
86	<	if(integrableObject->isDirectional()){
87	<	integrableObject->setA(rotMat * integrableObject->getA());
88	<	integrableObject->setJ(V3Zero);
89	<	}
78	>	integrableObject = mol->nextIntegrableObject(ii), ++i) {
79	>
80	>	newCoor = rotMat * refCoords[i];
81	>	newCoor += offset;
82	>
83	>	integrableObject->setPos(newCoor);
84	>	integrableObject->setVel(V3Zero);
85	>
86	>	if(integrableObject->isDirectional()){
87	>	integrableObject->setA(rotMat * integrableObject->getA());
88	>	integrableObject->setJ(V3Zero);
89	>	}
90		}
91	<	}
90	<
91	<	void MoLocator::calcRef( void ){
92	<	AtomStamp* currAtomStamp;
93	<	int nAtoms;
94	<	int nRigidBodies;
95	<	std::vector<double> mass;
96	<	Vector3d coor;
97	<	Vector3d refMolCom;
98	<	int nAtomsInRb;
99	<	double totMassInRb;
100	<	double currAtomMass;
101	<	double molMass;
91	>	}
92
93	<	nAtoms= myStamp->getNAtoms();
94	<	nRigidBodies = myStamp->getNRigidBodies();
95	<
96	<	for(size_t i=0; i<nAtoms; i++){
97	<
98	<	currAtomStamp = myStamp->getAtom(i);
99	<
100	<	if( !currAtomStamp->havePosition() ){
101	<	sprintf( painCave.errMsg,
102	<	"MoLocator error.\n"
103	<	"  Component %s, atom %s does not have a position specified.\n"
104	<	"  This means MoLocator cannot initalize it's position.\n",
105	<	myStamp->getID(),
106	<	currAtomStamp->getType() );
107	<
108	<	painCave.isFatal = 1;
109	<	simError();
93	>	void MoLocator::calcRef( void ){
94	>	AtomStamp* currAtomStamp;
95	>	RigidBodyStamp* rbStamp;
96	>	int nAtoms;
97	>	int nRigidBodies;
98	>	std::vector<RealType> mass;
99	>	Vector3d coor;
100	>	Vector3d refMolCom;
101	>	int nAtomsInRb;
102	>	RealType totMassInRb;
103	>	RealType currAtomMass;
104	>	RealType molMass;
105	>
106	>	nAtoms= myStamp->getNAtoms();
107	>	nRigidBodies = myStamp->getNRigidBodies();
108	>
109	>	for(size_t i=0; i<nAtoms; i++){
110	>
111	>	currAtomStamp = myStamp->getAtomStamp(i);
112	>
113	>	if( !currAtomStamp->havePosition() ){
114	>	sprintf( painCave.errMsg,
115	>	"MoLocator error.\n"
116	>	"  Component %s, atom %s does not have a position specified.\n"
117	>	"  This means MoLocator cannot initalize it's position.\n",
118	>	myStamp->getName().c_str(),
119	>	currAtomStamp->getType().c_str());
120	>
121	>	painCave.isFatal = 1;
122	>	simError();
123	>	}
124	>
125	>	//if atom belongs to rigidbody, just skip it
126	>	if(myStamp->isAtomInRigidBody(i))
127	>	continue;
128	>	//get mass and the reference coordinate
129	>	else{
130	>	currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);
131	>	mass.push_back(currAtomMass);
132	>	coor.x() = currAtomStamp->getPosX();
133	>	coor.y() = currAtomStamp->getPosY();
134	>	coor.z() = currAtomStamp->getPosZ();
135	>	refCoords.push_back(coor);
136	>
137	>	}
138		}
139	<
140	<	//if atom belongs to rigidbody, just skip it
141	<	if(myStamp->isAtomInRigidBody(i))
142	<	continue;
143	<	//get mass and the reference coordinate
144	<	else{
145	<	currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);
146	<	mass.push_back(currAtomMass);
147	<	coor.x() = currAtomStamp->getPosX();
148	<	coor.y() = currAtomStamp->getPosY();
149	<	coor.z() = currAtomStamp->getPosZ();
139	>
140	>	for(int i = 0; i < nRigidBodies; i++){
141	>
142	>	rbStamp = myStamp->getRigidBodyStamp(i);
143	>	nAtomsInRb = rbStamp->getNMembers();
144	>
145	>	coor.x() = 0.0;
146	>	coor.y() = 0.0;
147	>	coor.z() = 0.0;
148	>	totMassInRb = 0.0;
149	>
150	>	for(int j = 0; j < nAtomsInRb; j++){
151	>
152	>	currAtomStamp = myStamp->getAtomStamp(rbStamp->getMemberAt(j));
153	>	currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);
154	>	totMassInRb +=  currAtomMass;
155	>
156	>	coor.x() += currAtomStamp->getPosX() * currAtomMass;
157	>	coor.y() += currAtomStamp->getPosY() * currAtomMass;
158	>	coor.z() += currAtomStamp->getPosZ() * currAtomMass;
159	>	}
160	>
161	>	mass.push_back(totMassInRb);
162	>	coor /= totMassInRb;
163		refCoords.push_back(coor);
133	–
164		}
165	<	}
166	<
167	<	for(int i = 0; i < nRigidBodies; i++){
168	<	coor.x() = 0;
169	<	coor.y() = 0;
170	<	coor.z() = 0;
171	<	totMassInRb = 0;
172	<
173	<	for(int j = 0; j < nAtomsInRb; j++){
174	<
175	<	currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);
146	<	totMassInRb +=  currAtomMass;
147	<
148	<	coor.x() += currAtomStamp->getPosX() * currAtomMass;
149	<	coor.y() += currAtomStamp->getPosY() * currAtomMass;
150	<	coor.z() += currAtomStamp->getPosZ() * currAtomMass;
165	>
166	>
167	>	//calculate the reference center of mass
168	>	molMass = 0;
169	>	refMolCom.x() = 0;
170	>	refMolCom.y() = 0;
171	>	refMolCom.z() = 0;
172	>
173	>	for(int i = 0; i < nIntegrableObjects; i++){
174	>	refMolCom += refCoords[i] * mass[i];
175	>	molMass += mass[i];
176		}
177	<
178	<	mass.push_back(totMassInRb);
179	<	coor /= totMassInRb;
180	<	refCoords.push_back(coor);
177	>
178	>	refMolCom /= molMass;
179	>
180	>	//move the reference center of mass to (0,0,0) and adjust the reference coordinate
181	>	//of the integrabel objects
182	>	for(int i = 0; i < nIntegrableObjects; i++)
183	>	refCoords[i] -= refMolCom;
184		}
157	–
158	–
159	–	//calculate the reference center of mass
160	–	molMass = 0;
161	–	refMolCom.x() = 0;
162	–	refMolCom.y() = 0;
163	–	refMolCom.z() = 0;
185
186	<	for(int i = 0; i < nIntegrableObjects; i++){
187	<	refMolCom += refCoords[i] * mass[i];
167	<	molMass += mass[i];
168	<	}
169	<
170	<	refMolCom /= molMass;
171	<
172	<	//move the reference center of mass to (0,0,0) and adjust the reference coordinate
173	<	//of the integrabel objects
174	<	for(int i = 0; i < nIntegrableObjects; i++)
175	<	refCoords[i] -= refMolCom;
176	<	}
177	<
178	<
179	<
180	<	double getAtomMass(const std::string& at, ForceField* myFF) {
181	<	double mass;
186	>	RealType getAtomMass(const std::string& at, ForceField* myFF) {
187	>	RealType mass;
188		AtomType* atomType= myFF->getAtomType(at);
189		if (atomType != NULL) {
190	<	mass =     atomType->getMass();
190	>	mass =     atomType->getMass();
191		} else {
192	<	mass = 0.0;
193	<	std::cerr << "Can not find AtomType: " << at << std::endl;
192	>	mass = 0.0;
193	>	std::cerr << "Can not find AtomType: " << at << std::endl;
194		}
195		return mass;
196	<	}
197	<
198	<	double getMolMass(MoleculeStamp *molStamp, ForceField *myFF) {
196	>	}
197	>
198	>	RealType getMolMass(MoleculeStamp *molStamp, ForceField *myFF) {
199		int nAtoms;
200	<	double totMass = 0;
200	>	RealType totMass = 0;
201		nAtoms = molStamp->getNAtoms();
202	<
202	>
203		for(size_t i = 0; i < nAtoms; i++) {
204	<	AtomStamp *currAtomStamp = molStamp->getAtom(i);
205	<	totMass += getAtomMass(currAtomStamp->getType(), myFF);
204	>	AtomStamp *currAtomStamp = molStamp->getAtomStamp(i);
205	>	totMass += getAtomMass(currAtomStamp->getType(), myFF);
206		}
207		return totMass;
208	<	}
209	<	RotMat3x3d latVec2RotMat(const Vector3d& lv){
210	<
211	<	double theta =acos(lv[2]);
212	<	double phi = atan2(lv[1], lv[0]);
213	<	double psi = 0;
214	<
208	>	}
209	>	RotMat3x3d latVec2RotMat(const Vector3d& lv){
210	>
211	>	RealType theta =acos(lv[2]);
212	>	RealType phi = atan2(lv[1], lv[0]);
213	>	RealType psi = 0;
214	>
215		return RotMat3x3d(phi, theta, psi);
216	<
216	>
217	>	}
218		}
212	–	}
219

Comparing trunk/src/utils/MoLocator.cpp (property svn:keywords):
Revision 490 by tim, Wed Apr 13 18:41:17 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