applications/simpleBuilder/MoLocator.cpp

#include <iostream>

#include <cstdlib>
#include <cmath>

#include "utils/simError.h"
#include "applications/simpleBuilder/MoLocator.hpp"
#include "types/AtomType.hpp"
namespace oopse {
MoLocator::MoLocator( MoleculeStamp* theStamp, ForceField* theFF){

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

void MoLocator::placeMol( const Vector3d& offset, const Vector3d& ort, Molecule* mol){
    Vector3d newCoor;
    Vector3d curRefCoor;  
    RotMat3x3d rotMat = latVec2RotMat(ort);
  
    if(mol->getNIntegrableObjects() != 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();
    }

    Molecule::IntegrableObjectIterator ii;
    StuntDouble* integrableObject;
    int i;
    for (integrableObject = mol->beginIntegrableObject(ii), i = 0; integrableObject != NULL; 
        integrableObject = mol->nextIntegrableObject(ii), ++i) { 

        newCoor = rotMat * refCoords[i];
        newCoor += offset;

        integrableObject->setPos( newCoor);
        integrableObject->setVel(V3Zero);

        if(integrableObject->isDirectional()){
          integrableObject->setA(rotMat * integrableObject->getA());
          integrableObject->setJ(V3Zero);  
        }        
    }
}
 
void MoLocator::calcRef( void ){
  AtomStamp* currAtomStamp;
  int nAtoms; 
  int nRigidBodies;
   std::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{
    
      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 = getAtomMass(currAtomStamp->getType(), myFF);
      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;
}


double getAtomMass(const std::string& at, ForceField* myFF) {
    double mass;
    AtomType* atomType= myFF->getAtomType(at);
    if (atomType != NULL) {
        mass =     atomType->getMass();
    } else {
        mass = 0.0;
        std::cerr << "Can not find AtomType: " << at << std::endl;
    }
    return mass;
}

double getMolMass(MoleculeStamp *molStamp, ForceField *myFF) {
    int nAtoms;
    double totMass = 0;
    nAtoms = molStamp->getNAtoms();

    for(size_t i = 0; i < nAtoms; i++) {
        AtomStamp *currAtomStamp = molStamp->getAtom(i);
        totMass += getAtomMass(currAtomStamp->getType(), myFF);         
    }
    return totMass;
}
RotMat3x3d latVec2RotMat(const Vector3d& lv){

    double theta =acos(lv[2]);
    double phi = atan2(lv[1], lv[0]);
    double psi = 0;

    return RotMat3x3d(phi, theta, psi);

}
}

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