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){
  Vector3d newCoor;
  Vector3d velocity(0.0, 0.0, 0.0);
  Vector3d angMomentum(0.0, 0.0, 0.0);
  double quaternion[4];
  vector<StuntDouble*> myIntegrableObjects;
  double rotMat[3][3];
  
  quaternion[0] = 1.0;
  quaternion[1] = 0.0;
  quaternion[2] = 0.0;
  quaternion[3] = 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
    matVecMul3(rotMat, refCoords[i].vec, newCoor.vec);    
    newCoor +=  offset;

    myIntegrableObjects[i]->setPos( newCoor.vec);
    myIntegrableObjects[i]->setVel(velocity.vec);

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

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