ViewVC Help

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

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

Comparing trunk/src/primitives/DirectionalAtom.cpp (file contents):
Revision 663 by gezelter, Wed Oct 12 21:00:59 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
43		#include "primitives/DirectionalAtom.hpp"
44	+	#include "types/DirectionalAdapter.hpp"
45	+	#include "types/MultipoleAdapter.hpp"
46		#include "utils/simError.h"
47	<	namespace oopse {
47	>	namespace OpenMD {
48	>
49	>	DirectionalAtom::DirectionalAtom(AtomType* dAtomType)
50	>	: Atom(dAtomType) {
51	>	objType_= otDAtom;
52
53	<	DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType)
54	<	: Atom(dAtomType){
48	<	objType_= otDAtom;
49	<	if (dAtomType->isMultipole()) {
50	<	electroBodyFrame_ = dAtomType->getElectroBodyFrame();
51	<	}
53	>	DirectionalAdapter da = DirectionalAdapter(dAtomType);
54	>	I_ = da.getI();
55
56	<	//check if one of the diagonal inertia tensor of this directional atom  is zero
57	<	int nLinearAxis = 0;
58	<	Mat3x3d inertiaTensor = getI();
59	<	for (int i = 0; i < 3; i++) {
57	<	if (fabs(inertiaTensor(i, i)) < oopse::epsilon) {
58	<	linear_ = true;
59	<	linearAxis_ = i;
60	<	++ nLinearAxis;
61	<	}
62	<	}
56	>	MultipoleAdapter ma = MultipoleAdapter(dAtomType);
57	>	if (ma.isMultipole()) {
58	>	electroBodyFrame_ = ma.getElectroBodyFrame();
59	>	}
60
61	<	if (nLinearAxis > 1) {
62	<	sprintf( painCave.errMsg,
63	<	"Directional Atom warning.\n"
64	<	"\tOOPSE found more than one axis in this directional atom with a vanishing \n"
65	<	"\tmoment of inertia.");
66	<	painCave.isFatal = 0;
67	<	simError();
61	>	// Check if one of the diagonal inertia tensor of this directional
62	>	// atom is zero:
63	>	int nLinearAxis = 0;
64	>	Mat3x3d inertiaTensor = getI();
65	>	for (int i = 0; i < 3; i++) {
66	>	if (fabs(inertiaTensor(i, i)) < OpenMD::epsilon) {
67	>	linear_ = true;
68	>	linearAxis_ = i;
69	>	++ nLinearAxis;
70		}
72	–
71		}
72
73	+	if (nLinearAxis > 1) {
74	+	sprintf( painCave.errMsg,
75	+	"Directional Atom warning.\n"
76	+	"\tOpenMD found more than one axis in this directional atom with a vanishing \n"
77	+	"\tmoment of inertia.");
78	+	painCave.isFatal = 0;
79	+	simError();
80	+	}
81	+	}
82	+
83		Mat3x3d DirectionalAtom::getI() {
84	<	return static_cast<DirectionalAtomType*>(getAtomType())->getI();
84	>	return I_;
85		}
86	<
86	>
87		void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
88		((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
89		if (atomType_->isMultipole()) {
90		((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
91		}
92		}
93	<
94	<
93	>
94	>
95		void DirectionalAtom::setA(const RotMat3x3d& a) {
96		((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
97	<
97	>
98		if (atomType_->isMultipole()) {
99		((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
100		}
101		}
102	<
102	>
103		void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
104		((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
105	<
105	>
106		if (atomType_->isMultipole()) {
107		((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
108		}
109		}
110	<
110	>
111		void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
112		setA(m *getA());
113		}
114	<
115	<	std::vector<double> DirectionalAtom::getGrad() {
116	<	std::vector<double> grad(6, 0.0);
114	>
115	>	std::vector<RealType> DirectionalAtom::getGrad() {
116	>	std::vector<RealType> grad(6, 0.0);
117		Vector3d force;
118		Vector3d torque;
119		Vector3d myEuler;
120	<	double phi, theta, psi;
121	<	double cphi, sphi, ctheta, stheta;
120	>	RealType phi, theta;
121	>	// RealType psi;
122	>	RealType cphi, sphi, ctheta, stheta;
123		Vector3d ephi;
124		Vector3d etheta;
125		Vector3d epsi;
126	<
126	>
127		force = getFrc();
128		torque =getTrq();
129		myEuler = getA().toEulerAngles();
130	<
130	>
131		phi = myEuler[0];
132		theta = myEuler[1];
133	<	psi = myEuler[2];
134	<
133	>	// psi = myEuler[2];
134	>
135		cphi = cos(phi);
136		sphi = sin(phi);
137		ctheta = cos(theta);
138		stheta = sin(theta);
139	<
139	>
140		// get unit vectors along the phi, theta and psi rotation axes
141	<
141	>
142		ephi[0] = 0.0;
143		ephi[1] = 0.0;
144		ephi[2] = 1.0;
145	<
145	>
146	>	//etheta[0] = -sphi;
147	>	//etheta[1] =  cphi;
148	>	//etheta[2] =  0.0;
149	>
150		etheta[0] = cphi;
151		etheta[1] = sphi;
152		etheta[2] = 0.0;
153	<
153	>
154		epsi[0] = stheta * cphi;
155		epsi[1] = stheta * sphi;
156		epsi[2] = ctheta;
157	<
157	>
158		//gradient is equal to -force
159		for (int j = 0 ; j<3; j++)
160		grad[j] = -force[j];
161	<
162	<	for (int j = 0; j < 3; j++ ) {
150	<
161	>
162	>	for (int j = 0; j < 3; j++ ) {
163		grad[3] -= torque[j]*ephi[j];
164		grad[4] -= torque[j]*etheta[j];
165	<	grad[5] -= torque[j]*epsi[j];
154	<
165	>	grad[5] -= torque[j]*epsi[j];
166		}
167
168		return grad;
169		}
170	<
170	>
171		void DirectionalAtom::accept(BaseVisitor* v) {
172		v->visit(this);
173	<	}
163	<
173	>	}
174		}
175

Comparing trunk/src/primitives/DirectionalAtom.cpp (property svn:keywords):
Revision 663 by gezelter, Wed Oct 12 21:00:59 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