ViewVC Help

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

root/OpenMD/trunk/src/primitives/RigidBody.cpp

Comparing trunk/src/primitives/RigidBody.cpp (file contents):
Revision 642 by tim, Mon Oct 3 15:54:23 2005 UTC vs.
Revision 1797 by gezelter, Mon Sep 10 20:58:00 2012 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]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	+	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42		#include <algorithm>
43		#include <math.h>
44		#include "primitives/RigidBody.hpp"
45		#include "utils/simError.h"
46		#include "utils/NumericConstant.hpp"
47	<	namespace oopse {
48	<
49	<	RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData), inertiaTensor_(0.0){
50	<
47	>	namespace OpenMD {
48	>
49	>	RigidBody::RigidBody() : StuntDouble(otRigidBody, &Snapshot::rigidbodyData),
50	>	inertiaTensor_(0.0){
51		}
52	<
52	>
53		void RigidBody::setPrevA(const RotMat3x3d& a) {
54		((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
55	<	//((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
56	<
56	<	for (int i =0 ; i < atoms_.size(); ++i){
55	>
56	>	for (unsigned int i = 0 ; i < atoms_.size(); ++i){
57		if (atoms_[i]->isDirectional()) {
58	<	atoms_[i]->setPrevA(a * refOrients_[i]);
58	>	atoms_[i]->setPrevA(refOrients_[i].transpose() * a);
59		}
60		}
61	<
61	>
62		}
63	<
64	<
63	>
64	>
65		void RigidBody::setA(const RotMat3x3d& a) {
66		((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
67	–	//((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
67
68	<	for (int i =0 ; i < atoms_.size(); ++i){
68	>	for (unsigned int i = 0 ; i < atoms_.size(); ++i){
69		if (atoms_[i]->isDirectional()) {
70	<	atoms_[i]->setA(a * refOrients_[i]);
70	>	atoms_[i]->setA(refOrients_[i].transpose() * a);
71		}
72		}
73		}
74	<
74	>
75		void RigidBody::setA(const RotMat3x3d& a, int snapshotNo) {
76		((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
77	+
78		//((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * sU_;
79	<
80	<	for (int i =0 ; i < atoms_.size(); ++i){
79	>
80	>	for (unsigned int i = 0 ; i < atoms_.size(); ++i){
81		if (atoms_[i]->isDirectional()) {
82	<	atoms_[i]->setA(a * refOrients_[i], snapshotNo);
82	>	atoms_[i]->setA(refOrients_[i].transpose() * a, snapshotNo);
83		}
84		}
85	<
85	>
86		}
87	<
87	>
88		Mat3x3d RigidBody::getI() {
89		return inertiaTensor_;
90		}
91	<
92	<	std::vector<double> RigidBody::getGrad() {
93	<	std::vector<double> grad(6, 0.0);
91	>
92	>	std::vector<RealType> RigidBody::getGrad() {
93	>	std::vector<RealType> grad(6, 0.0);
94		Vector3d force;
95		Vector3d torque;
96		Vector3d myEuler;
97	<	double phi, theta, psi;
98	<	double cphi, sphi, ctheta, stheta;
97	>	RealType phi, theta;
98	>	// RealType psi;
99	>	RealType cphi, sphi, ctheta, stheta;
100		Vector3d ephi;
101		Vector3d etheta;
102		Vector3d epsi;
103	<
103	>
104		force = getFrc();
105		torque =getTrq();
106		myEuler = getA().toEulerAngles();
107	<
107	>
108		phi = myEuler[0];
109		theta = myEuler[1];
110	<	psi = myEuler[2];
111	<
110	>	// psi = myEuler[2];
111	>
112		cphi = cos(phi);
113		sphi = sin(phi);
114		ctheta = cos(theta);
115		stheta = sin(theta);
116	<
116	>
117		// get unit vectors along the phi, theta and psi rotation axes
118	<
118	>
119		ephi[0] = 0.0;
120		ephi[1] = 0.0;
121		ephi[2] = 1.0;
122	<
122	>
123	>	//etheta[0] = -sphi;
124	>	//etheta[1] =  cphi;
125	>	//etheta[2] =  0.0;
126	>
127		etheta[0] = cphi;
128		etheta[1] = sphi;
129	<	etheta[2] = 0.0;
130	<
129	>	etheta[2] =  0.0;
130	>
131		epsi[0] = stheta * cphi;
132		epsi[1] = stheta * sphi;
133		epsi[2] = ctheta;
134	<
134	>
135		//gradient is equal to -force
136		for (int j = 0 ; j<3; j++)
137		grad[j] = -force[j];
138	<
138	>
139		for (int j = 0; j < 3; j++ ) {
140	<
140	>
141		grad[3] += torque[j]*ephi[j];
142		grad[4] += torque[j]*etheta[j];
143		grad[5] += torque[j]*epsi[j];
144	<
144	>
145		}
146
147		return grad;
148		}
149	<
149	>
150		void RigidBody::accept(BaseVisitor* v) {
151		v->visit(this);
152		}
153
154		/**@todo need modification */
155		void  RigidBody::calcRefCoords() {
156	<	double mtmp;
156	>	RealType mtmp;
157		Vector3d refCOM(0.0);
158		mass_ = 0.0;
159		for (std::size_t i = 0; i < atoms_.size(); ++i) {
#	Line 157 | Line 162 | namespace oopse {
162		refCOM += refCoords_[i]*mtmp;
163		}
164		refCOM /= mass_;
165	<
165	>
166		// Next, move the origin of the reference coordinate system to the COM:
167		for (std::size_t i = 0; i < atoms_.size(); ++i) {
168		refCoords_[i] -= refCOM;
#	Line 169 | Line 174 | namespace oopse {
174		Mat3x3d IAtom(0.0);
175		mtmp = atoms_[i]->getMass();
176		IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp;
177	<	double r2 = refCoords_[i].lengthSquare();
177	>	RealType r2 = refCoords_[i].lengthSquare();
178		IAtom(0, 0) += mtmp * r2;
179		IAtom(1, 1) += mtmp * r2;
180		IAtom(2, 2) += mtmp * r2;
181	<
181	>	Itmp += IAtom;
182	>
183		//project the inertial moment of directional atoms into this rigid body
184		if (atoms_[i]->isDirectional()) {
185	<	IAtom += atoms_[i]->getI();
186	<	Itmp += refOrients_[i].transpose() * IAtom * refOrients_[i];
181	<	} else {
182	<	Itmp += IAtom;
183	<	}
185	>	Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i];
186	>	}
187		}
188
189	+	//    std::cout << Itmp << std::endl;
190	+
191		//diagonalize
192		Vector3d evals;
193		Mat3x3d::diagonalize(Itmp, evals, sU_);
#	Line 194 | Line 199 | namespace oopse {
199
200		int nLinearAxis = 0;
201		for (int i = 0; i < 3; i++) {
202	<	if (fabs(evals[i]) < oopse::epsilon) {
202	>	if (fabs(evals[i]) < OpenMD::epsilon) {
203		linear_ = true;
204		linearAxis_ = i;
205		++ nLinearAxis;
#	Line 204 | Line 209 | namespace oopse {
209		if (nLinearAxis > 1) {
210		sprintf( painCave.errMsg,
211		"RigidBody error.\n"
212	<	"\tOOPSE found more than one axis in this rigid body with a vanishing \n"
212	>	"\tOpenMD found more than one axis in this rigid body with a vanishing \n"
213		"\tmoment of inertia.  This can happen in one of three ways:\n"
214		"\t 1) Only one atom was specified, or \n"
215		"\t 2) All atoms were specified at the same location, or\n"
#	Line 223 | Line 228 | namespace oopse {
228		Vector3d apos;
229		Vector3d rpos;
230		Vector3d frc(0.0);
231	<	Vector3d trq(0.0);
231	>	Vector3d trq(0.0);
232		Vector3d pos = this->getPos();
233	<	for (int i = 0; i < atoms_.size(); i++) {
233	>	for (unsigned int i = 0; i < atoms_.size(); i++) {
234
235		afrc = atoms_[i]->getFrc();
236		apos = atoms_[i]->getPos();
#	Line 243 | Line 248 | namespace oopse {
248		if (atoms_[i]->isDirectional()) {
249		atrq = atoms_[i]->getTrq();
250		trq += atrq;
251	<	}
251	>	}
252	>	}
253	>	addFrc(frc);
254	>	addTrq(trq);
255	>	}
256	>
257	>	Mat3x3d RigidBody::calcForcesAndTorquesAndVirial() {
258	>	Vector3d afrc;
259	>	Vector3d atrq;
260	>	Vector3d apos;
261	>	Vector3d rpos;
262	>	Vector3d dfrc;
263	>	Vector3d frc(0.0);
264	>	Vector3d trq(0.0);
265	>	Vector3d pos = this->getPos();
266	>	Mat3x3d tau_(0.0);
267	>
268	>	for (unsigned int i = 0; i < atoms_.size(); i++) {
269	>
270	>	afrc = atoms_[i]->getFrc();
271	>	apos = atoms_[i]->getPos();
272	>	rpos = apos - pos;
273
274	+	frc += afrc;
275	+
276	+	trq[0] += rpos[1]*afrc[2] - rpos[2]*afrc[1];
277	+	trq[1] += rpos[2]*afrc[0] - rpos[0]*afrc[2];
278	+	trq[2] += rpos[0]*afrc[1] - rpos[1]*afrc[0];
279	+
280	+	// If the atom has a torque associated with it, then we also need to
281	+	// migrate the torques onto the center of mass:
282	+
283	+	if (atoms_[i]->isDirectional()) {
284	+	atrq = atoms_[i]->getTrq();
285	+	trq += atrq;
286	+	}
287	+
288	+	tau_(0,0) -= rpos[0]*afrc[0];
289	+	tau_(0,1) -= rpos[0]*afrc[1];
290	+	tau_(0,2) -= rpos[0]*afrc[2];
291	+	tau_(1,0) -= rpos[1]*afrc[0];
292	+	tau_(1,1) -= rpos[1]*afrc[1];
293	+	tau_(1,2) -= rpos[1]*afrc[2];
294	+	tau_(2,0) -= rpos[2]*afrc[0];
295	+	tau_(2,1) -= rpos[2]*afrc[1];
296	+	tau_(2,2) -= rpos[2]*afrc[2];
297	+
298		}
299	<
300	<	setFrc(frc);
301	<	setTrq(trq);
252	<
299	>	addFrc(frc);
300	>	addTrq(trq);
301	>	return tau_;
302		}
303
304		void  RigidBody::updateAtoms() {
#	Line 271 | Line 320 | namespace oopse {
320		if (atoms_[i]->isDirectional()) {
321
322		dAtom = (DirectionalAtom *) atoms_[i];
323	<	dAtom->setA(refOrients_[i] * a);
323	>	dAtom->setA(refOrients_[i].transpose() * a);
324		}
325
326		}
#	Line 298 | Line 347 | namespace oopse {
347		if (atoms_[i]->isDirectional()) {
348
349		dAtom = (DirectionalAtom *) atoms_[i];
350	<	dAtom->setA(refOrients_[i] * a, frame);
350	>	dAtom->setA(refOrients_[i].transpose() * a, frame);
351		}
352
353		}
#	Line 328 | Line 377 | namespace oopse {
377
378
379		Vector3d velRot;
380	<	for (int i =0 ; i < refCoords_.size(); ++i) {
380	>	for (unsigned int i = 0 ; i < refCoords_.size(); ++i) {
381		atoms_[i]->setVel(rbVel + mat * refCoords_[i]);
382		}
383
#	Line 357 | Line 406 | namespace oopse {
406
407
408		Vector3d velRot;
409	<	for (int i =0 ; i < refCoords_.size(); ++i) {
409	>	for (unsigned int i = 0 ; i < refCoords_.size(); ++i) {
410		atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame);
411		}
412
#	Line 483 | Line 532 | namespace oopse {
532		"RigidBody error.\n"
533		"\tAtom %s does not have a position specified.\n"
534		"\tThis means RigidBody cannot set up reference coordinates.\n",
535	<	ats->getType() );
535	>	ats->getType().c_str() );
536		painCave.isFatal = 1;
537		simError();
538		}
#	Line 503 | Line 552 | namespace oopse {
552		"RigidBody error.\n"
553		"\tAtom %s does not have an orientation specified.\n"
554		"\tThis means RigidBody cannot set up reference orientations.\n",
555	<	ats->getType() );
555	>	ats->getType().c_str() );
556		painCave.isFatal = 1;
557		simError();
558		}

Comparing trunk/src/primitives/RigidBody.cpp (property svn:keywords):
Revision 642 by tim, Mon Oct 3 15:54:23 2005 UTC vs.
Revision 1797 by gezelter, Mon Sep 10 20:58:00 2012 UTC

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

Diff Legend

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