utils/sysbuilder/MoLocator.cpp

#include <iostream>

#include <cstdlib>
#include <cmath>

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

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;

  quaternion[0] = 1.0;
  quaternion[1] = 0.0;
  quaternion[2] = 0.0;
  quaternion[3] = 0.0;

  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++) {

      newCoor = refCoords[i] + offset;
      myIntegrableObjects[i]->setPos( newCoor.vec);
      myIntegrableObjects[i]->setVel(velocity.vec);

      if(myIntegrableObjects[i]->isDirectional()){
        myIntegrableObjects[i]->setQ(quaternion);
        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 totMass;
  
  mass.resize(nIntegrableObjects);
  
  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
  for(int i = 0; i < nIntegrableObjects; i++){
    refMolCom += refCoords[i] * mass[i];
    totMass += mass[i];
  }
  
  refMolCom /= totMass;

  //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;
}

Revision:	1429
Committed:	Wed Jul 28 20:29:49 2004 UTC (21 years, 3 months ago) by tim
File size:	3531 byte(s)
Log Message:	simpleBuilder in progress
#	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
9	tim	1427	MoLocator::MoLocator( MoleculeStamp* theStamp, ForceFields* theFF){
10	gezelter	1334
11			myStamp = theStamp;
12	tim	1427	myFF = theFF;
13			nIntegrableObjects = myStamp->getNIntegrable();
14	gezelter	1334	calcRefCoords();
15			}
16
17	tim	1427	void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
18			Vector3d newCoor;
19			Vector3d velocity(0.0, 0.0, 0.0);
20			Vector3d angMomentum(0.0, 0.0, 0.0);
21			double quaternion[4];
22			vector<StuntDouble*> myIntegrableObjects;
23	gezelter	1334
24	tim	1427	quaternion[0] = 1.0;
25			quaternion[1] = 0.0;
26			quaternion[2] = 0.0;
27			quaternion[3] = 0.0;
28	gezelter	1334
29	tim	1427	myIntegrableObjects = mol->getIntegrableObjects();
30	gezelter	1334
31	tim	1427	if(myIntegrableObjects.size() != nIntegrableObjects){
32			sprintf( painCave.errMsg,
33			"MoLocator error.\n"
34			" The number of integrable objects of MoleculeStamp is not the same as that of Molecule\n");
35			painCave.isFatal = 1;
36			simError();
37	gezelter	1334
38	tim	1427	}
39	gezelter	1334
40	tim	1427	for(int i=0; i<nIntegrableObjects; i++) {
41	gezelter	1334
42	tim	1427	newCoor = refCoords[i] + offset;
43			myIntegrableObjects[i]->setPos( newCoor.vec);
44			myIntegrableObjects[i]->setVel(velocity.vec);
45	gezelter	1334
46	tim	1427	if(myIntegrableObjects[i]->isDirectional()){
47			myIntegrableObjects[i]->setQ(quaternion);
48			myIntegrableObjects[i]->setJ(angMomentum.vec);
49			}
50			}
51	gezelter	1334
52			}
53	tim	1427
54	gezelter	1334	void MoLocator::calcRefCoords( void ){
55	tim	1427	AtomStamp* currAtomStamp;
56			int nAtoms;
57			int nRigidBodies;
58			vector<double> mass;
59			Vector3d coor;
60			Vector3d refMolCom;
61			int nAtomsInRb;
62			double totMassInRb;
63			double currAtomMass;
64	gezelter	1334
65	tim	1427	double totMass;
66	gezelter	1334
67	tim	1427	mass.resize(nIntegrableObjects);
68
69			nAtoms= myStamp->getNAtoms();
70			nRigidBodies = myStamp->getNRigidBodies();
71	gezelter	1334
72	tim	1427	//
73			for(size_t i=0; i<nAtoms; i++){
74
75			currAtomStamp = myStamp->getAtom(i);
76
77			if( !currAtomStamp->havePosition() ){
78	gezelter	1334	sprintf( painCave.errMsg,
79	tim	1427	"MoLocator error.\n"
80			" Component %s, atom %s does not have a position specified.\n"
81			" This means MoLocator cannot initalize it's position.\n",
82			myStamp->getID(),
83			currAtomStamp->getType() );
84
85	gezelter	1334	painCave.isFatal = 1;
86			simError();
87			}
88
89	tim	1427	//if atom belongs to rigidbody, just skip it
90			if(myStamp->isAtomInRigidBody(i))
91			continue;
92			//get mass and the reference coordinate
93			else{
94			currAtomMass = myFF->getAtomTypeMass(currAtomStamp->getType());
95			mass.push_back(currAtomMass);
96			coor.x = currAtomStamp->getPosX();
97			coor.y = currAtomStamp->getPosY();
98			coor.z = currAtomStamp->getPosZ();
99			refCoords.push_back(coor);
100
101			}
102	gezelter	1334	}
103
104	tim	1427	for(int i = 0; i < nRigidBodies; i++){
105			coor.x = 0;
106			coor.y = 0;
107			coor.z = 0;
108			totMassInRb = 0;
109	gezelter	1334
110	tim	1427	for(int j = 0; j < nAtomsInRb; j++){
111	gezelter	1334
112	tim	1427	currAtomMass = myFF->getAtomTypeMass(currAtomStamp->getType());
113			totMassInRb += currAtomMass;
114
115			coor.x += currAtomStamp->getPosX() * currAtomMass;
116			coor.y += currAtomStamp->getPosY() * currAtomMass;
117			coor.z += currAtomStamp->getPosZ() * currAtomMass;
118			}
119	gezelter	1334
120	tim	1427	mass.push_back(totMassInRb);
121			coor /= totMassInRb;
122			refCoords.push_back(coor);
123	gezelter	1334	}
124
125	tim	1427	//calculate the reference center of mass
126			for(int i = 0; i < nIntegrableObjects; i++){
127			refMolCom += refCoords[i] * mass[i];
128			totMass += mass[i];
129	gezelter	1334	}
130	tim	1429
131			refMolCom /= totMass;
132	gezelter	1334
133	tim	1427	//move the reference center of mass to (0,0,0) and adjust the reference coordinate
134			//of the integrabel objects
135			for(int i = 0; i < nIntegrableObjects; i++)
136			refCoords[i] -= refMolCom;
137	gezelter	1334	}
138