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

Revision:	1432
Committed:	Thu Jul 29 03:31:50 2004 UTC (19 years, 11 months ago) by tim
File size:	3600 byte(s)
Log Message:	simpleBuilder is built but Makefile is broken
#	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	tim	1432	double molMass;
65	gezelter	1334
66	tim	1427	mass.resize(nIntegrableObjects);
67
68			nAtoms= myStamp->getNAtoms();
69			nRigidBodies = myStamp->getNRigidBodies();
70	gezelter	1334
71	tim	1427	for(size_t i=0; i<nAtoms; i++){
72
73			currAtomStamp = myStamp->getAtom(i);
74
75			if( !currAtomStamp->havePosition() ){
76	gezelter	1334	sprintf( painCave.errMsg,
77	tim	1427	"MoLocator error.\n"
78			" Component %s, atom %s does not have a position specified.\n"
79			" This means MoLocator cannot initalize it's position.\n",
80			myStamp->getID(),
81			currAtomStamp->getType() );
82
83	gezelter	1334	painCave.isFatal = 1;
84			simError();
85			}
86
87	tim	1427	//if atom belongs to rigidbody, just skip it
88			if(myStamp->isAtomInRigidBody(i))
89			continue;
90			//get mass and the reference coordinate
91			else{
92			currAtomMass = myFF->getAtomTypeMass(currAtomStamp->getType());
93			mass.push_back(currAtomMass);
94			coor.x = currAtomStamp->getPosX();
95			coor.y = currAtomStamp->getPosY();
96			coor.z = currAtomStamp->getPosZ();
97			refCoords.push_back(coor);
98
99			}
100	gezelter	1334	}
101
102	tim	1427	for(int i = 0; i < nRigidBodies; i++){
103			coor.x = 0;
104			coor.y = 0;
105			coor.z = 0;
106			totMassInRb = 0;
107	gezelter	1334
108	tim	1427	for(int j = 0; j < nAtomsInRb; j++){
109	gezelter	1334
110	tim	1427	currAtomMass = myFF->getAtomTypeMass(currAtomStamp->getType());
111			totMassInRb += currAtomMass;
112
113			coor.x += currAtomStamp->getPosX() * currAtomMass;
114			coor.y += currAtomStamp->getPosY() * currAtomMass;
115			coor.z += currAtomStamp->getPosZ() * currAtomMass;
116			}
117	gezelter	1334
118	tim	1427	mass.push_back(totMassInRb);
119			coor /= totMassInRb;
120			refCoords.push_back(coor);
121	gezelter	1334	}
122
123	tim	1432
124	tim	1427	//calculate the reference center of mass
125	tim	1432	molMass = 0;
126			refMolCom.x = 0;
127			refMolCom.y = 0;
128			refMolCom.z = 0;
129
130	tim	1427	for(int i = 0; i < nIntegrableObjects; i++){
131			refMolCom += refCoords[i] * mass[i];
132	tim	1432	molMass += mass[i];
133	gezelter	1334	}
134	tim	1429
135	tim	1432	refMolCom /= molMass;
136	gezelter	1334
137	tim	1427	//move the reference center of mass to (0,0,0) and adjust the reference coordinate
138			//of the integrabel objects
139			for(int i = 0; i < nIntegrableObjects; i++)
140			refCoords[i] -= refMolCom;
141	gezelter	1334	}
142