src/applications/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:	2
Committed:	Fri Sep 24 04:16:43 2004 UTC (20 years, 7 months ago) by gezelter
File size:	4501 byte(s)
Log Message:	Import of OOPSE v. 2.0
#	User	Rev	Content
1	gezelter	2	#include <iostream>
2
3			#include <cstdlib>
4			#include <cmath>
5
6			#include "simError.h"
7			#include "MoLocator.hpp"
8			#include "MatVec3.h"
9
10			MoLocator::MoLocator( MoleculeStamp* theStamp, ForceFields* theFF){
11
12			myStamp = theStamp;
13			myFF = theFF;
14			nIntegrableObjects = myStamp->getNIntegrable();
15			calcRefCoords();
16			}
17
18			void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
19			double newCoor[3];
20			double curRefCoor[3];
21			double zeroVector[3];
22			vector<StuntDouble*> myIntegrableObjects;
23			double rotMat[3][3];
24
25			zeroVector[0] = 0.0;
26			zeroVector[1] = 0.0;
27			zeroVector[2] = 0.0;
28
29			latVec2RotMat(ort, rotMat);
30
31			myIntegrableObjects = mol->getIntegrableObjects();
32
33			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
40			}
41
42			for(int i=0; i<nIntegrableObjects; i++) {
43
44			//calculate the reference coordinate for integrable objects after rotation
45			curRefCoor[0] = refCoords[i][0];
46			curRefCoor[1] = refCoords[i][1];
47			curRefCoor[2] = refCoords[i][2];
48
49			matVecMul3(rotMat, curRefCoor, newCoor);
50
51			newCoor[0] += offset[0];
52			newCoor[1] += offset[1];
53			newCoor[2] += offset[2];
54
55			myIntegrableObjects[i]->setPos( newCoor);
56			myIntegrableObjects[i]->setVel(zeroVector);
57
58			if(myIntegrableObjects[i]->isDirectional()){
59			myIntegrableObjects[i]->setA(rotMat);
60			myIntegrableObjects[i]->setJ(zeroVector);
61			}
62			}
63
64			}
65
66			void MoLocator::calcRefCoords( void ){
67			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			double molMass;
77
78			nAtoms= myStamp->getNAtoms();
79			nRigidBodies = myStamp->getNRigidBodies();
80
81			for(size_t i=0; i<nAtoms; i++){
82
83			currAtomStamp = myStamp->getAtom(i);
84
85			if( !currAtomStamp->havePosition() ){
86			sprintf( painCave.errMsg,
87			"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			painCave.isFatal = 1;
94			simError();
95			}
96
97			//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			}
111
112			for(int i = 0; i < nRigidBodies; i++){
113			coor.x = 0;
114			coor.y = 0;
115			coor.z = 0;
116			totMassInRb = 0;
117
118			for(int j = 0; j < nAtomsInRb; j++){
119
120			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
128			mass.push_back(totMassInRb);
129			coor /= totMassInRb;
130			refCoords.push_back(coor);
131			}
132
133
134			//calculate the reference center of mass
135			molMass = 0;
136			refMolCom.x = 0;
137			refMolCom.y = 0;
138			refMolCom.z = 0;
139
140			for(int i = 0; i < nIntegrableObjects; i++){
141			refMolCom += refCoords[i] * mass[i];
142			molMass += mass[i];
143			}
144
145			refMolCom /= molMass;
146
147			//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			}
152
153
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