src/primitives/Molecule.cpp

#include <stdlib.h>


#include "primitives/Molecule.hpp"
#include "utils/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;
  
  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();

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

}

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

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

  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:	1683
Committed:	Thu Oct 28 22:34:02 2004 UTC (19 years, 8 months ago)
File size:	4704 byte(s)
Log Message:	This commit was manufactured by cvs2svn to create branch 'new_design'.
#	User	Rev	Content
1	gezelter	1490	#include <stdlib.h>
2
3
4	tim	1492	#include "primitives/Molecule.hpp"
5			#include "utils/simError.h"
6	gezelter	1490
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			CutoffGroup* cg;
20			vector<CutoffGroup*>::iterator iter;
21
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			for(cg = beginCutoffGroup(iter); cg != NULL; cg = nextCutoffGroup(iter))
43			delete cg;
44			myCutoffGroups.clear();
45
46			}
47
48
49			void Molecule::initialize( molInit &theInit ){
50
51			CutoffGroup* curCutoffGroup;
52			vector<CutoffGroup*>::iterator iterCutoff;
53			Atom* cutoffAtom;
54			vector<Atom*>::iterator iterAtom;
55			int atomIndex;
56
57			nAtoms = theInit.nAtoms;
58			nMembers = nAtoms;
59			nBonds = theInit.nBonds;
60			nBends = theInit.nBends;
61			nTorsions = theInit.nTorsions;
62			nRigidBodies = theInit.nRigidBodies;
63			nOriented = theInit.nOriented;
64
65			myAtoms = theInit.myAtoms;
66			myBonds = theInit.myBonds;
67			myBends = theInit.myBends;
68			myTorsions = theInit.myTorsions;
69			myRigidBodies = theInit.myRigidBodies;
70
71			myIntegrableObjects = theInit.myIntegrableObjects;
72
73			for (int i = 0; i < myRigidBodies.size(); i++)
74			myRigidBodies[i]->calcRefCoords();
75
76			myCutoffGroups = theInit.myCutoffGroups;
77			nCutoffGroups = myCutoffGroups.size();
78
79			}
80
81			void Molecule::calcForces( void ){
82
83			int i;
84			double com[3];
85
86			for(i=0; i<myRigidBodies.size(); i++) {
87			myRigidBodies[i]->updateAtoms();
88			}
89
90			for(i=0; i<nBonds; i++){
91			myBonds[i]->calc_forces();
92			}
93
94			for(i=0; i<nBends; i++){
95			myBends[i]->calc_forces();
96			}
97
98			for(i=0; i<nTorsions; i++){
99			myTorsions[i]->calc_forces();
100			}
101
102			// Rigid Body forces and torques are done after the fortran force loop
103
104			}
105
106
107			double Molecule::getPotential( void ){
108
109			int i;
110			double myPot = 0.0;
111
112			for(i=0; i<myRigidBodies.size(); i++) {
113			myRigidBodies[i]->updateAtoms();
114			}
115
116			for(i=0; i<nBonds; i++){
117			myPot += myBonds[i]->get_potential();
118			}
119
120			for(i=0; i<nBends; i++){
121			myPot += myBends[i]->get_potential();
122			}
123
124			for(i=0; i<nTorsions; i++){
125			myPot += myTorsions[i]->get_potential();
126			}
127
128			return myPot;
129			}
130
131			void Molecule::printMe( void ){
132
133			int i;
134
135			for(i=0; i<nBonds; i++){
136			myBonds[i]->printMe();
137			}
138
139			for(i=0; i<nBends; i++){
140			myBends[i]->printMe();
141			}
142
143			for(i=0; i<nTorsions; i++){
144			myTorsions[i]->printMe();
145			}
146
147			}
148
149			void Molecule::moveCOM(double delta[3]){
150			double aPos[3];
151			int i, j;
152
153			for(i=0; i<myIntegrableObjects.size(); i++) {
154			if(myIntegrableObjects[i] != NULL ) {
155
156			myIntegrableObjects[i]->getPos( aPos );
157
158			for (j=0; j< 3; j++)
159			aPos[j] += delta[j];
160
161			myIntegrableObjects[i]->setPos( aPos );
162			}
163			}
164
165			for(i=0; i<myRigidBodies.size(); i++) {
166
167			myRigidBodies[i]->getPos( aPos );
168
169			for (j=0; j< 3; j++)
170			aPos[j] += delta[j];
171
172			myRigidBodies[i]->setPos( aPos );
173			}
174			}
175
176			void Molecule::atoms2rigidBodies( void ) {
177			int i;
178			for (i = 0; i < myRigidBodies.size(); i++) {
179			myRigidBodies[i]->calcForcesAndTorques();
180			}
181			}
182
183			void Molecule::getCOM( double COM[3] ) {
184
185			double mass, mtot;
186			double aPos[3];
187			int i, j;
188
189			for (j=0; j<3; j++)
190			COM[j] = 0.0;
191
192			mtot = 0.0;
193
194			for (i=0; i < myIntegrableObjects.size(); i++) {
195			if (myIntegrableObjects[i] != NULL) {
196
197			mass = myIntegrableObjects[i]->getMass();
198			mtot += mass;
199
200			myIntegrableObjects[i]->getPos( aPos );
201
202			for( j = 0; j < 3; j++)
203			COM[j] += aPos[j] * mass;
204
205			}
206			}
207
208			for (j = 0; j < 3; j++)
209			COM[j] /= mtot;
210			}
211
212			double Molecule::getCOMvel( double COMvel[3] ) {
213
214			double mass, mtot;
215			double aVel[3];
216			int i, j;
217
218
219			for (j=0; j<3; j++)
220			COMvel[j] = 0.0;
221
222			mtot = 0.0;
223
224			for (i=0; i < myIntegrableObjects.size(); i++) {
225			if (myIntegrableObjects[i] != NULL) {
226
227			mass = myIntegrableObjects[i]->getMass();
228			mtot += mass;
229
230			myIntegrableObjects[i]->getVel(aVel);
231
232			for (j=0; j<3; j++)
233			COMvel[j] += aVel[j]*mass;
234
235			}
236			}
237
238			for (j=0; j<3; j++)
239			COMvel[j] /= mtot;
240
241			return mtot;
242
243			}
244
245			double Molecule::getTotalMass()
246			{
247
248			double totalMass;
249
250			totalMass = 0;
251			for(int i =0; i < myIntegrableObjects.size(); i++){
252			totalMass += myIntegrableObjects[i]->getMass();
253			}
254
255			return totalMass;
256			}