utils/sysbuilder/MoLocator.cpp

#include <iostream>

#include <cstdlib>
#include <cmath>

#include "simError.h"
#include "MoLocator.hpp"
#include "MatVec3.h"

MoLocator::MoLocator( MoleculeStamp* theStamp, ForceFields* theFF){

  myStamp = theStamp;
  myFF = theFF;
  nIntegrableObjects = myStamp->getNIntegrable();
  calcRefCoords();
}

void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
  double newCoor[3];
  double curRefCoor[3];
  double zeroVector[3];
  vector<StuntDouble*> myIntegrableObjects;
  double rotMat[3][3];
  
  zeroVector[0] = 0.0;
  zeroVector[1] = 0.0;
  zeroVector[2] = 0.0;
  
  latVec2RotMat(ort, rotMat);
  
  myIntegrableObjects = mol->getIntegrableObjects();

  if(myIntegrableObjects.size() != nIntegrableObjects){
      sprintf( painCave.errMsg,
               "MoLocator error.\n"
               "  The number of integrable objects of MoleculeStamp is not the same as  that of Molecule\n");
          painCave.isFatal = 1;
      simError();

  }
    
  for(int i=0; i<nIntegrableObjects; i++) {

    //calculate the reference coordinate for integrable objects after rotation
    curRefCoor[0] = refCoords[i][0];
    curRefCoor[1] = refCoords[i][1];
    curRefCoor[2] = refCoords[i][2];
    
    matVecMul3(rotMat, curRefCoor, newCoor);    

    newCoor[0] +=  offset[0];
    newCoor[1] +=  offset[1];
    newCoor[2] +=  offset[2];

    myIntegrableObjects[i]->setPos( newCoor);
    myIntegrableObjects[i]->setVel(zeroVector);

    if(myIntegrableObjects[i]->isDirectional()){
      myIntegrableObjects[i]->setA(rotMat);
      myIntegrableObjects[i]->setJ(zeroVector);  
    }
  }

}
 
void MoLocator::calcRefCoords( void ){
  AtomStamp* currAtomStamp;
  int nAtoms; 
  int nRigidBodies;
  vector<double> mass;
  Vector3d coor;
  Vector3d refMolCom;  
  int nAtomsInRb;
  double totMassInRb;
  double currAtomMass;
  double molMass;
  
  nAtoms= myStamp->getNAtoms();
  nRigidBodies = myStamp->getNRigidBodies();

  for(size_t i=0; i<nAtoms; i++){

    currAtomStamp = myStamp->getAtom(i);

    if( !currAtomStamp->havePosition() ){
      sprintf( painCave.errMsg,
                  "MoLocator error.\n"
                  "  Component %s, atom %s does not have a position specified.\n"
                  "  This means MoLocator cannot initalize it's position.\n",
                  myStamp->getID(),
                  currAtomStamp->getType() );

      painCave.isFatal = 1;
      simError();
    }

    //if atom belongs to rigidbody, just skip it
    if(myStamp->isAtomInRigidBody(i))
      continue;
    //get mass and the reference coordinate 
    else{
      currAtomMass = myFF->getAtomTypeMass(currAtomStamp->getType());
      mass.push_back(currAtomMass);
      coor.x = currAtomStamp->getPosX();
      coor.y = currAtomStamp->getPosY();
      coor.z = currAtomStamp->getPosZ();
      refCoords.push_back(coor);

    }
  }

  for(int i = 0; i < nRigidBodies; i++){
    coor.x = 0;
    coor.y = 0;
    coor.z = 0;
    totMassInRb = 0;

    for(int j = 0; j < nAtomsInRb; j++){

      currAtomMass = myFF->getAtomTypeMass(currAtomStamp->getType());
      totMassInRb +=  currAtomMass;
      
      coor.x += currAtomStamp->getPosX() * currAtomMass;
      coor.y += currAtomStamp->getPosY() * currAtomMass;
      coor.z += currAtomStamp->getPosZ() * currAtomMass;
    }

    mass.push_back(totMassInRb);
    coor /= totMassInRb;
    refCoords.push_back(coor);
  }


  //calculate the reference center of mass
  molMass = 0;
  refMolCom.x = 0;
  refMolCom.y = 0;
  refMolCom.z = 0;
  
  for(int i = 0; i < nIntegrableObjects; i++){
    refMolCom += refCoords[i] * mass[i];
   molMass += mass[i];
  }
  
  refMolCom /= molMass;

  //move the reference center of mass to (0,0,0) and adjust the reference coordinate 
  //of the integrabel objects
  for(int i = 0; i < nIntegrableObjects; i++)
    refCoords[i] -= refMolCom;
}


void latVec2RotMat(const Vector3d& lv, double rotMat[3][3]){

  double theta, phi, psi;
  
  theta =acos(lv.z);
  phi = atan2(lv.y, lv.x);
  psi = 0;

  rotMat[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
  rotMat[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
  rotMat[0][2] = sin(theta) * sin(psi);
  
  rotMat[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
  rotMat[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
  rotMat[1][2] = sin(theta) * cos(psi);
  
  rotMat[2][0] = sin(phi) * sin(theta);
  rotMat[2][1] = -cos(phi) * sin(theta);
  rotMat[2][2] = cos(theta);
}

Revision:	1443
Committed:	Fri Jul 30 14:44:10 2004 UTC (19 years, 11 months ago) by tim
File size:	4501 byte(s)
Log Message:	remove anonymous struct from Vector3d.hpp
#	User	Rev	Content
1	gezelter	1334	#include <iostream>
2
3			#include <cstdlib>
4			#include <cmath>
5
6			#include "simError.h"
7			#include "MoLocator.hpp"
8	tim	1435	#include "MatVec3.h"
9	gezelter	1334
10	tim	1427	MoLocator::MoLocator( MoleculeStamp* theStamp, ForceFields* theFF){
11	gezelter	1334
12			myStamp = theStamp;
13	tim	1427	myFF = theFF;
14			nIntegrableObjects = myStamp->getNIntegrable();
15	gezelter	1334	calcRefCoords();
16			}
17
18	tim	1427	void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
19	tim	1443	double newCoor[3];
20			double curRefCoor[3];
21			double zeroVector[3];
22	tim	1427	vector<StuntDouble*> myIntegrableObjects;
23	tim	1435	double rotMat[3][3];
24
25	tim	1443	zeroVector[0] = 0.0;
26			zeroVector[1] = 0.0;
27			zeroVector[2] = 0.0;
28
29	tim	1435	latVec2RotMat(ort, rotMat);
30
31	tim	1427	myIntegrableObjects = mol->getIntegrableObjects();
32	gezelter	1334
33	tim	1427	if(myIntegrableObjects.size() != nIntegrableObjects){
34			sprintf( painCave.errMsg,
35			"MoLocator error.\n"
36			" The number of integrable objects of MoleculeStamp is not the same as that of Molecule\n");
37			painCave.isFatal = 1;
38			simError();
39	gezelter	1334
40	tim	1427	}
41	gezelter	1334
42	tim	1427	for(int i=0; i<nIntegrableObjects; i++) {
43	gezelter	1334
44	tim	1435	//calculate the reference coordinate for integrable objects after rotation
45	tim	1443	curRefCoor[0] = refCoords[i][0];
46			curRefCoor[1] = refCoords[i][1];
47			curRefCoor[2] = refCoords[i][2];
48
49			matVecMul3(rotMat, curRefCoor, newCoor);
50	gezelter	1334
51	tim	1443	newCoor[0] += offset[0];
52			newCoor[1] += offset[1];
53			newCoor[2] += offset[2];
54	tim	1435
55	tim	1443	myIntegrableObjects[i]->setPos( newCoor);
56			myIntegrableObjects[i]->setVel(zeroVector);
57
58	tim	1435	if(myIntegrableObjects[i]->isDirectional()){
59			myIntegrableObjects[i]->setA(rotMat);
60	tim	1443	myIntegrableObjects[i]->setJ(zeroVector);
61	tim	1435	}
62	tim	1427	}
63	gezelter	1334
64			}
65	tim	1427
66	gezelter	1334	void MoLocator::calcRefCoords( void ){
67	tim	1427	AtomStamp* currAtomStamp;
68			int nAtoms;
69			int nRigidBodies;
70			vector<double> mass;
71			Vector3d coor;
72			Vector3d refMolCom;
73			int nAtomsInRb;
74			double totMassInRb;
75			double currAtomMass;
76	tim	1432	double molMass;
77	gezelter	1334
78	tim	1427	nAtoms= myStamp->getNAtoms();
79			nRigidBodies = myStamp->getNRigidBodies();
80	gezelter	1334
81	tim	1427	for(size_t i=0; i<nAtoms; i++){
82
83			currAtomStamp = myStamp->getAtom(i);
84
85			if( !currAtomStamp->havePosition() ){
86	gezelter	1334	sprintf( painCave.errMsg,
87	tim	1427	"MoLocator error.\n"
88			" Component %s, atom %s does not have a position specified.\n"
89			" This means MoLocator cannot initalize it's position.\n",
90			myStamp->getID(),
91			currAtomStamp->getType() );
92
93	gezelter	1334	painCave.isFatal = 1;
94			simError();
95			}
96
97	tim	1427	//if atom belongs to rigidbody, just skip it
98			if(myStamp->isAtomInRigidBody(i))
99			continue;
100			//get mass and the reference coordinate
101			else{
102			currAtomMass = myFF->getAtomTypeMass(currAtomStamp->getType());
103			mass.push_back(currAtomMass);
104			coor.x = currAtomStamp->getPosX();
105			coor.y = currAtomStamp->getPosY();
106			coor.z = currAtomStamp->getPosZ();
107			refCoords.push_back(coor);
108
109			}
110	gezelter	1334	}
111
112	tim	1427	for(int i = 0; i < nRigidBodies; i++){
113			coor.x = 0;
114			coor.y = 0;
115			coor.z = 0;
116			totMassInRb = 0;
117	gezelter	1334
118	tim	1427	for(int j = 0; j < nAtomsInRb; j++){
119	gezelter	1334
120	tim	1427	currAtomMass = myFF->getAtomTypeMass(currAtomStamp->getType());
121			totMassInRb += currAtomMass;
122
123			coor.x += currAtomStamp->getPosX() * currAtomMass;
124			coor.y += currAtomStamp->getPosY() * currAtomMass;
125			coor.z += currAtomStamp->getPosZ() * currAtomMass;
126			}
127	gezelter	1334
128	tim	1427	mass.push_back(totMassInRb);
129			coor /= totMassInRb;
130			refCoords.push_back(coor);
131	gezelter	1334	}
132
133	tim	1432
134	tim	1427	//calculate the reference center of mass
135	tim	1432	molMass = 0;
136			refMolCom.x = 0;
137			refMolCom.y = 0;
138			refMolCom.z = 0;
139
140	tim	1427	for(int i = 0; i < nIntegrableObjects; i++){
141			refMolCom += refCoords[i] * mass[i];
142	tim	1432	molMass += mass[i];
143	gezelter	1334	}
144	tim	1429
145	tim	1432	refMolCom /= molMass;
146	gezelter	1334
147	tim	1427	//move the reference center of mass to (0,0,0) and adjust the reference coordinate
148			//of the integrabel objects
149			for(int i = 0; i < nIntegrableObjects; i++)
150			refCoords[i] -= refMolCom;
151	gezelter	1334	}
152
153	tim	1435
154			void latVec2RotMat(const Vector3d& lv, double rotMat[3][3]){
155
156			double theta, phi, psi;
157
158			theta =acos(lv.z);
159			phi = atan2(lv.y, lv.x);
160			psi = 0;
161
162			rotMat[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
163			rotMat[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
164			rotMat[0][2] = sin(theta) * sin(psi);
165
166			rotMat[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
167			rotMat[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
168			rotMat[1][2] = sin(theta) * cos(psi);
169
170			rotMat[2][0] = sin(phi) * sin(theta);
171			rotMat[2][1] = -cos(phi) * sin(theta);
172			rotMat[2][2] = cos(theta);
173			}
174