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 2071 by gezelter, Sat Mar 7 21:41:51 2015 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	>	std::vector<RealType> mass;
100	>	Vector3d coor;
101	>	Vector3d refMolCom;
102	>	RealType totMassInRb;
103	>	RealType currAtomMass;
104	>	RealType molMass;
105	>
106	>	std::size_t nAtoms= myStamp->getNAtoms();
107	>	std::size_t nRigidBodies = myStamp->getNRigidBodies();
108	>
109	>	for(std::size_t i = 0; i < nAtoms; i++){
110	>
111	>	currAtomStamp = myStamp->getAtomStamp(i);
112	>
113	>	if( !currAtomStamp->havePosition() ){
114	>	sprintf( painCave.errMsg,
115	>	"MoLocator::calcRef error.\n"
116	>	"\tComponent %s, atom %s does not have a position specified.\n"
117	>	"\tThis 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		}
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{
139
140	<	mass.push_back(currAtomMass);
141	<	coor.x() = currAtomStamp->getPosX();
142	<	coor.y() = currAtomStamp->getPosY();
143	<	coor.z() = currAtomStamp->getPosZ();
140	>	for(std::size_t i = 0; i < nRigidBodies; i++){
141	>
142	>	rbStamp = myStamp->getRigidBodyStamp(i);
143	>	std::size_t 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(std::size_t 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(std::size_t 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
181	>	//reference coordinate of the integrabel objects
182	>	for(std::size_t 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 MoLocator::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) {
199	<	int nAtoms;
200	<	double totMass = 0;
196	>	}
197	>
198	>	RealType MoLocator::getMolMass(MoleculeStamp *molStamp, ForceField *myFF) {
199	>	unsigned int nAtoms;
200	>	RealType totMass = 0;
201		nAtoms = molStamp->getNAtoms();
202	<
203	<	for(size_t i = 0; i < nAtoms; i++) {
204	<	AtomStamp *currAtomStamp = molStamp->getAtom(i);
205	<	totMass += getAtomMass(currAtomStamp->getType(), myFF);
202	>
203	>	for(std::size_t i = 0; i < nAtoms; i++) {
204	>	AtomStamp *currAtomStamp = molStamp->getAtomStamp(i);
205	>	totMass += getAtomMass(currAtomStamp->getType(), myFF);
206		}
207		return totMass;
208	<	}
203	<	RotMat3x3d latVec2RotMat(const Vector3d& lv){
208	>	}
209
210	<	double theta =acos(lv[2]);
211	<	double phi = atan2(lv[1], lv[0]);
212	<	double psi = 0;
213	<
214	<	return RotMat3x3d(phi, theta, psi);
215	<
210	>	RotMat3x3d MoLocator::latVec2RotMat(const Vector3d& lv){
211	>
212	>	RealType theta =acos(lv[2]);
213	>	RealType phi = atan2(lv[1], lv[0]);
214	>	RealType psi = 0;
215	>
216	>	return RotMat3x3d(phi, theta, psi);
217	>	}
218		}
212	–	}
219

Comparing trunk/src/utils/MoLocator.cpp (property svn:keywords):
Revision 483 by gezelter, Tue Apr 12 21:28:07 2005 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

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

Diff Legend

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