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

Comparing trunk/src/utils/MoLocator.cpp (property svn:keywords):
Revision 483 by gezelter, Tue Apr 12 21:28:07 2005 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

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

Diff Legend

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