OOPSE/libmdtools/Molecule.cpp

#include <stdlib.h>


#include "Molecule.hpp"
#include "simError.h"


Molecule::Molecule( void ){

  myAtoms = NULL;
  myBonds = NULL;
  myBends = NULL;
  myTorsions = NULL;
}

Molecule::~Molecule( void ){
  int i;
  CutoffGroup* cg;
  vector<CutoffGroup*>::iterator iter;
  
  if( myAtoms != NULL ){
    for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
    delete[] myAtoms;
  }

  if( myBonds != NULL ){
    for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i];
    delete[] myBonds;
  }

  if( myBends != NULL ){
    for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i];
    delete[] myBends;
  }

  if( myTorsions != NULL ){
    for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i];
    delete[] myTorsions;
  }

  for(cg = beginCutoffGroup(iter);  cg != NULL; cg = nextCutoffGroup(iter))
    delete cg;
  myCutoffGroups.clear();
  
}


void Molecule::initialize( molInit &theInit ){

  CutoffGroup* curCutoffGroup;
  vector<CutoffGroup*>::iterator iterCutoff;
  Atom* cutoffAtom;
  vector<Atom*>::iterator iterAtom;
  int atomIndex;
  GenericData* gdata;
  ConsRbData* rbData;
  RigidBody* oldRb;
  
  nAtoms = theInit.nAtoms;
  nMembers = nAtoms;
  nBonds = theInit.nBonds;
  nBends = theInit.nBends;
  nTorsions = theInit.nTorsions;
  nRigidBodies = theInit.nRigidBodies;
  nOriented = theInit.nOriented;

  myAtoms = theInit.myAtoms;
  myBonds = theInit.myBonds;
  myBends = theInit.myBends;
  myTorsions = theInit.myTorsions;
  myRigidBodies = theInit.myRigidBodies;

  myIntegrableObjects = theInit.myIntegrableObjects;

  for (int i = 0; i < myRigidBodies.size(); i++){
      myRigidBodies[i]->calcRefCoords();
    //just a quick hack
    
    gdata = myRigidBodies[i]->getProperty("OldState");
    if(gdata != NULL){
      rbData = dynamic_cast<ConsRbData*>(gdata);
      if(rbData ==NULL)
        cerr << "dynamic_cast to ConsRbData Error in Molecule::initialize()" << endl;
      else{
        oldRb = rbData->getData();
        oldRb->calcRefCoords();
      }
    }//end if(gata != NULL)
    
  }//end for(int i = 0; i < myRigidBodies.size(); i++)

  myCutoffGroups = theInit.myCutoffGroups;
  nCutoffGroups = myCutoffGroups.size();

  myConstraintPairs = theInit.myConstraintPairs;
  
}

void Molecule::calcForces( void ){
  
  int i;
  double com[3];

  for(i=0; i<myRigidBodies.size(); i++) {
    myRigidBodies[i]->updateAtoms();
  }

  //calculate the center of mass of the molecule
  //getCOM(com);  
  //for(int i = 0; i < nAtoms; i ++)
  //  myAtoms[i]->setRc(com);  
  

  for(i=0; i<nBonds; i++){
    myBonds[i]->calc_forces();
  }

  for(i=0; i<nBends; i++){
    myBends[i]->calc_forces();
  }

  for(i=0; i<nTorsions; i++){
    myTorsions[i]->calc_forces();
  }

  // Rigid Body forces and torques are done after the fortran force loop

}


double Molecule::getPotential( void ){
  
  int i;
  double myPot = 0.0;

  for(i=0; i<myRigidBodies.size(); i++) {
    myRigidBodies[i]->updateAtoms();
  }
   
  for(i=0; i<nBonds; i++){
    myPot += myBonds[i]->get_potential();
  }

  for(i=0; i<nBends; i++){
    myPot += myBends[i]->get_potential();
  }

  for(i=0; i<nTorsions; i++){
    myPot += myTorsions[i]->get_potential();
  }

  return myPot;
}

void Molecule::printMe( void ){
  
  int i;

  for(i=0; i<nBonds; i++){
    myBonds[i]->printMe();
  }

  for(i=0; i<nBends; i++){
    myBends[i]->printMe();
  }

  for(i=0; i<nTorsions; i++){
    myTorsions[i]->printMe();
  }

}

void Molecule::moveCOM(double delta[3]){
  double aPos[3];
  int i, j;

  for(i=0; i<myIntegrableObjects.size(); i++) {
    if(myIntegrableObjects[i] != NULL ) {
      
      myIntegrableObjects[i]->getPos( aPos );
      
      for (j=0; j< 3; j++) 
        aPos[j] += delta[j];

      myIntegrableObjects[i]->setPos( aPos );
    }
  }

  for(i=0; i<myRigidBodies.size(); i++) {

      myRigidBodies[i]->getPos( aPos );

      for (j=0; j< 3; j++) 
        aPos[j] += delta[j];
      
      myRigidBodies[i]->setPos( aPos );
    }
}

void Molecule::atoms2rigidBodies( void ) {
  int i;
  for (i = 0; i < myRigidBodies.size(); i++) {
    myRigidBodies[i]->calcForcesAndTorques();   
  }
}

void Molecule::getCOM( double COM[3] ) {

  double mass, mtot;
  double aPos[3];
  int i, j;

  for (j=0; j<3; j++) 
    COM[j] = 0.0;

  mtot   = 0.0;

  for (i=0; i < myIntegrableObjects.size(); i++) {
    if (myIntegrableObjects[i] != NULL) {

      mass = myIntegrableObjects[i]->getMass();
      mtot   += mass;
      
      myIntegrableObjects[i]->getPos( aPos );

      for( j = 0; j < 3; j++) 
        COM[j] += aPos[j] * mass;

    }
  }

  for (j = 0; j < 3; j++) 
    COM[j] /= mtot; 
}

double Molecule::getCOMvel( double COMvel[3] ) {

  double mass, mtot;
  double aVel[3];
  int i, j;


  for (j=0; j<3; j++) 
    COMvel[j] = 0.0;

  mtot   = 0.0;

  for (i=0; i < myIntegrableObjects.size(); i++) {
    if (myIntegrableObjects[i] != NULL) {

      mass = myIntegrableObjects[i]->getMass();
      mtot   += mass;

      myIntegrableObjects[i]->getVel(aVel);

      for (j=0; j<3; j++) 
        COMvel[j] += aVel[j]*mass;

    }
  }

  for (j=0; j<3; j++) 
    COMvel[j] /= mtot;
 
  return mtot;

}

double Molecule::getTotalMass()
{

  double totalMass;
  
  totalMass = 0;
  for(int i =0; i < myIntegrableObjects.size(); i++){
    totalMass += myIntegrableObjects[i]->getMass();
  }

  return totalMass;
}
Revision:	1452
Committed:	Mon Aug 23 15:11:36 2004 UTC (19 years, 10 months ago) by tim
File size:	5385 byte(s)
Log Message:	* empty log message *
#	User	Rev	Content
1	gezelter	829	#include <stdlib.h>
2	mmeineke	426
3
4			#include "Molecule.hpp"
5	chuckv	438	#include "simError.h"
6	mmeineke	426
7
8
9			Molecule::Molecule( void ){
10
11			myAtoms = NULL;
12			myBonds = NULL;
13			myBends = NULL;
14			myTorsions = NULL;
15			}
16
17			Molecule::~Molecule( void ){
18			int i;
19	tim	1167	CutoffGroup* cg;
20			vector<CutoffGroup*>::iterator iter;
21	mmeineke	426
22			if( myAtoms != NULL ){
23			for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
24			delete[] myAtoms;
25			}
26
27			if( myBonds != NULL ){
28			for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i];
29			delete[] myBonds;
30			}
31
32			if( myBends != NULL ){
33			for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i];
34			delete[] myBends;
35			}
36
37			if( myTorsions != NULL ){
38			for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i];
39			delete[] myTorsions;
40			}
41
42	tim	1167	for(cg = beginCutoffGroup(iter); cg != NULL; cg = nextCutoffGroup(iter))
43			delete cg;
44			myCutoffGroups.clear();
45
46	mmeineke	426	}
47
48
49			void Molecule::initialize( molInit &theInit ){
50	tim	1140
51	tim	1157	CutoffGroup* curCutoffGroup;
52			vector<CutoffGroup*>::iterator iterCutoff;
53			Atom* cutoffAtom;
54			vector<Atom*>::iterator iterAtom;
55			int atomIndex;
56	tim	1452	GenericData* gdata;
57			ConsRbData* rbData;
58			RigidBody* oldRb;
59	tim	1157
60	mmeineke	426	nAtoms = theInit.nAtoms;
61			nMembers = nAtoms;
62			nBonds = theInit.nBonds;
63			nBends = theInit.nBends;
64			nTorsions = theInit.nTorsions;
65	gezelter	1097	nRigidBodies = theInit.nRigidBodies;
66	mmeineke	426	nOriented = theInit.nOriented;
67
68			myAtoms = theInit.myAtoms;
69			myBonds = theInit.myBonds;
70			myBends = theInit.myBends;
71			myTorsions = theInit.myTorsions;
72	gezelter	1097	myRigidBodies = theInit.myRigidBodies;
73	gezelter	1104
74			myIntegrableObjects = theInit.myIntegrableObjects;
75	tim	1113
76	tim	1452	for (int i = 0; i < myRigidBodies.size(); i++){
77	tim	1113	myRigidBodies[i]->calcRefCoords();
78	tim	1452	//just a quick hack
79
80			gdata = myRigidBodies[i]->getProperty("OldState");
81			if(gdata != NULL){
82			rbData = dynamic_cast<ConsRbData*>(gdata);
83			if(rbData ==NULL)
84			cerr << "dynamic_cast to ConsRbData Error in Molecule::initialize()" << endl;
85			else{
86			oldRb = rbData->getData();
87			oldRb->calcRefCoords();
88			}
89			}//end if(gata != NULL)
90
91			}//end for(int i = 0; i < myRigidBodies.size(); i++)
92	tim	1136
93	tim	1157	myCutoffGroups = theInit.myCutoffGroups;
94	tim	1167	nCutoffGroups = myCutoffGroups.size();
95	tim	1234
96			myConstraintPairs = theInit.myConstraintPairs;
97	tim	1157
98	mmeineke	426	}
99
100			void Molecule::calcForces( void ){
101
102			int i;
103	tim	1136	double com[3];
104	mmeineke	426
105	gezelter	1104	for(i=0; i<myRigidBodies.size(); i++) {
106	gezelter	1097	myRigidBodies[i]->updateAtoms();
107			}
108
109	tim	1136	//calculate the center of mass of the molecule
110	tim	1157	//getCOM(com);
111			//for(int i = 0; i < nAtoms; i ++)
112			// myAtoms[i]->setRc(com);
113	tim	1136
114
115	mmeineke	426	for(i=0; i<nBonds; i++){
116			myBonds[i]->calc_forces();
117			}
118
119			for(i=0; i<nBends; i++){
120			myBends[i]->calc_forces();
121			}
122
123			for(i=0; i<nTorsions; i++){
124			myTorsions[i]->calc_forces();
125			}
126	gezelter	1097
127			// Rigid Body forces and torques are done after the fortran force loop
128
129	mmeineke	426	}
130
131
132	mmeineke	428	double Molecule::getPotential( void ){
133	mmeineke	426
134			int i;
135			double myPot = 0.0;
136	gezelter	1109
137			for(i=0; i<myRigidBodies.size(); i++) {
138			myRigidBodies[i]->updateAtoms();
139			}
140	mmeineke	426
141			for(i=0; i<nBonds; i++){
142			myPot += myBonds[i]->get_potential();
143			}
144
145			for(i=0; i<nBends; i++){
146			myPot += myBends[i]->get_potential();
147			}
148
149			for(i=0; i<nTorsions; i++){
150			myPot += myTorsions[i]->get_potential();
151			}
152
153			return myPot;
154			}
155	mmeineke	435
156			void Molecule::printMe( void ){
157
158			int i;
159
160			for(i=0; i<nBonds; i++){
161			myBonds[i]->printMe();
162			}
163
164			for(i=0; i<nBends; i++){
165			myBends[i]->printMe();
166			}
167
168			for(i=0; i<nTorsions; i++){
169			myTorsions[i]->printMe();
170			}
171	gezelter	1097
172	mmeineke	435	}
173	gezelter	446
174	mmeineke	449	void Molecule::moveCOM(double delta[3]){
175	gezelter	610	double aPos[3];
176			int i, j;
177	gezelter	446
178	tim	1113	for(i=0; i<myIntegrableObjects.size(); i++) {
179			if(myIntegrableObjects[i] != NULL ) {
180	gezelter	610
181	tim	1113	myIntegrableObjects[i]->getPos( aPos );
182	gezelter	610
183			for (j=0; j< 3; j++)
184			aPos[j] += delta[j];
185	gezelter	446
186	tim	1113	myIntegrableObjects[i]->setPos( aPos );
187	gezelter	446	}
188			}
189	gezelter	1097
190	gezelter	1104	for(i=0; i<myRigidBodies.size(); i++) {
191	gezelter	1097
192			myRigidBodies[i]->getPos( aPos );
193
194			for (j=0; j< 3; j++)
195			aPos[j] += delta[j];
196
197			myRigidBodies[i]->setPos( aPos );
198			}
199	gezelter	446	}
200
201	gezelter	1097	void Molecule::atoms2rigidBodies( void ) {
202			int i;
203	gezelter	1104	for (i = 0; i < myRigidBodies.size(); i++) {
204			myRigidBodies[i]->calcForcesAndTorques();
205	gezelter	1097	}
206			}
207
208	mmeineke	449	void Molecule::getCOM( double COM[3] ) {
209	gezelter	446
210			double mass, mtot;
211	gezelter	610	double aPos[3];
212			int i, j;
213	gezelter	446
214	gezelter	610	for (j=0; j<3; j++)
215			COM[j] = 0.0;
216
217	gezelter	446	mtot = 0.0;
218
219	tim	1113	for (i=0; i < myIntegrableObjects.size(); i++) {
220			if (myIntegrableObjects[i] != NULL) {
221	gezelter	446
222	tim	1113	mass = myIntegrableObjects[i]->getMass();
223	gezelter	446	mtot += mass;
224	gezelter	610
225	tim	1113	myIntegrableObjects[i]->getPos( aPos );
226	gezelter	446
227	gezelter	610	for( j = 0; j < 3; j++)
228			COM[j] += aPos[j] * mass;
229
230	gezelter	446	}
231			}
232
233	gezelter	610	for (j = 0; j < 3; j++)
234			COM[j] /= mtot;
235	gezelter	446	}
236	gezelter	468
237	gezelter	475	double Molecule::getCOMvel( double COMvel[3] ) {
238	gezelter	468
239	gezelter	475	double mass, mtot;
240	gezelter	607	double aVel[3];
241			int i, j;
242	gezelter	468
243	gezelter	607
244			for (j=0; j<3; j++)
245			COMvel[j] = 0.0;
246
247	gezelter	468	mtot = 0.0;
248
249	tim	1113	for (i=0; i < myIntegrableObjects.size(); i++) {
250			if (myIntegrableObjects[i] != NULL) {
251	mmeineke	489
252	tim	1113	mass = myIntegrableObjects[i]->getMass();
253	gezelter	468	mtot += mass;
254
255	tim	1113	myIntegrableObjects[i]->getVel(aVel);
256	gezelter	607
257			for (j=0; j<3; j++)
258			COMvel[j] += aVel[j]*mass;
259
260	gezelter	468	}
261			}
262
263	gezelter	607	for (j=0; j<3; j++)
264			COMvel[j] /= mtot;
265	gezelter	468
266	gezelter	475	return mtot;
267
268	gezelter	468	}
269	tim	658
270			double Molecule::getTotalMass()
271			{
272	tim	1113
273	tim	658	double totalMass;
274
275			totalMass = 0;
276	tim	1113	for(int i =0; i < myIntegrableObjects.size(); i++){
277			totalMass += myIntegrableObjects[i]->getMass();
278	tim	658	}
279
280			return totalMass;
281			}