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

}


void Molecule::initialize( molInit &theInit ){
  double totMass;
  
  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();


  //the mass ratio will never change during the simulation. Thus, we could
  //just calculate it at the begining of the simulation
  totMass = getTotalMass();
  for(int i = 0; i < nAtoms; i ++)
    myAtoms[i]->setMassRatio(myAtoms[i]->getMass()/totMass);  
}

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:	1136
Committed:	Tue Apr 27 16:26:44 2004 UTC (20 years, 2 months ago) by tim
File size:	4686 byte(s)
Log Message:	add center of mass of the molecule and massRation into atom class
#	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
18			Molecule::~Molecule( void ){
19			int i;
20
21			if( myAtoms != NULL ){
22			for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
23			delete[] myAtoms;
24			}
25
26			if( myBonds != NULL ){
27			for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i];
28			delete[] myBonds;
29			}
30
31			if( myBends != NULL ){
32			for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i];
33			delete[] myBends;
34			}
35
36			if( myTorsions != NULL ){
37			for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i];
38			delete[] myTorsions;
39			}
40
41			}
42
43
44			void Molecule::initialize( molInit &theInit ){
45	tim	1136	double totMass;
46
47	mmeineke	426	nAtoms = theInit.nAtoms;
48			nMembers = nAtoms;
49			nBonds = theInit.nBonds;
50			nBends = theInit.nBends;
51			nTorsions = theInit.nTorsions;
52	gezelter	1097	nRigidBodies = theInit.nRigidBodies;
53	mmeineke	426	nOriented = theInit.nOriented;
54
55			myAtoms = theInit.myAtoms;
56			myBonds = theInit.myBonds;
57			myBends = theInit.myBends;
58			myTorsions = theInit.myTorsions;
59	gezelter	1097	myRigidBodies = theInit.myRigidBodies;
60	gezelter	1104
61			myIntegrableObjects = theInit.myIntegrableObjects;
62	tim	1113
63			for (int i = 0; i < myRigidBodies.size(); i++)
64			myRigidBodies[i]->calcRefCoords();
65	tim	1136
66
67			//the mass ratio will never change during the simulation. Thus, we could
68			//just calculate it at the begining of the simulation
69			totMass = getTotalMass();
70			for(int i = 0; i < nAtoms; i ++)
71			myAtoms[i]->setMassRatio(myAtoms[i]->getMass()/totMass);
72	mmeineke	426	}
73
74			void Molecule::calcForces( void ){
75
76			int i;
77	tim	1136	double com[3];
78	mmeineke	426
79	gezelter	1104	for(i=0; i<myRigidBodies.size(); i++) {
80	gezelter	1097	myRigidBodies[i]->updateAtoms();
81			}
82
83	tim	1136	//calculate the center of mass of the molecule
84			getCOM(com);
85			for(int i = 0; i < nAtoms; i ++)
86			myAtoms[i]->setRc(com);
87
88
89	mmeineke	426	for(i=0; i<nBonds; i++){
90			myBonds[i]->calc_forces();
91			}
92
93			for(i=0; i<nBends; i++){
94			myBends[i]->calc_forces();
95			}
96
97			for(i=0; i<nTorsions; i++){
98			myTorsions[i]->calc_forces();
99			}
100	gezelter	1097
101			// Rigid Body forces and torques are done after the fortran force loop
102
103	mmeineke	426	}
104
105
106	mmeineke	428	double Molecule::getPotential( void ){
107	mmeineke	426
108			int i;
109			double myPot = 0.0;
110	gezelter	1109
111			for(i=0; i<myRigidBodies.size(); i++) {
112			myRigidBodies[i]->updateAtoms();
113			}
114	mmeineke	426
115			for(i=0; i<nBonds; i++){
116			myPot += myBonds[i]->get_potential();
117			}
118
119			for(i=0; i<nBends; i++){
120			myPot += myBends[i]->get_potential();
121			}
122
123			for(i=0; i<nTorsions; i++){
124			myPot += myTorsions[i]->get_potential();
125			}
126
127			return myPot;
128			}
129	mmeineke	435
130			void Molecule::printMe( void ){
131
132			int i;
133
134			for(i=0; i<nBonds; i++){
135			myBonds[i]->printMe();
136			}
137
138			for(i=0; i<nBends; i++){
139			myBends[i]->printMe();
140			}
141
142			for(i=0; i<nTorsions; i++){
143			myTorsions[i]->printMe();
144			}
145	gezelter	1097
146	mmeineke	435	}
147	gezelter	446
148	mmeineke	449	void Molecule::moveCOM(double delta[3]){
149	gezelter	610	double aPos[3];
150			int i, j;
151	gezelter	446
152	tim	1113	for(i=0; i<myIntegrableObjects.size(); i++) {
153			if(myIntegrableObjects[i] != NULL ) {
154	gezelter	610
155	tim	1113	myIntegrableObjects[i]->getPos( aPos );
156	gezelter	610
157			for (j=0; j< 3; j++)
158			aPos[j] += delta[j];
159	gezelter	446
160	tim	1113	myIntegrableObjects[i]->setPos( aPos );
161	gezelter	446	}
162			}
163	gezelter	1097
164	gezelter	1104	for(i=0; i<myRigidBodies.size(); i++) {
165	gezelter	1097
166			myRigidBodies[i]->getPos( aPos );
167
168			for (j=0; j< 3; j++)
169			aPos[j] += delta[j];
170
171			myRigidBodies[i]->setPos( aPos );
172			}
173	gezelter	446	}
174
175	gezelter	1097	void Molecule::atoms2rigidBodies( void ) {
176			int i;
177	gezelter	1104	for (i = 0; i < myRigidBodies.size(); i++) {
178			myRigidBodies[i]->calcForcesAndTorques();
179	gezelter	1097	}
180			}
181
182	mmeineke	449	void Molecule::getCOM( double COM[3] ) {
183	gezelter	446
184			double mass, mtot;
185	gezelter	610	double aPos[3];
186			int i, j;
187	gezelter	446
188	gezelter	610	for (j=0; j<3; j++)
189			COM[j] = 0.0;
190
191	gezelter	446	mtot = 0.0;
192
193	tim	1113	for (i=0; i < myIntegrableObjects.size(); i++) {
194			if (myIntegrableObjects[i] != NULL) {
195	gezelter	446
196	tim	1113	mass = myIntegrableObjects[i]->getMass();
197	gezelter	446	mtot += mass;
198	gezelter	610
199	tim	1113	myIntegrableObjects[i]->getPos( aPos );
200	gezelter	446
201	gezelter	610	for( j = 0; j < 3; j++)
202			COM[j] += aPos[j] * mass;
203
204	gezelter	446	}
205			}
206
207	gezelter	610	for (j = 0; j < 3; j++)
208			COM[j] /= mtot;
209	gezelter	446	}
210	gezelter	468
211	gezelter	475	double Molecule::getCOMvel( double COMvel[3] ) {
212	gezelter	468
213	gezelter	475	double mass, mtot;
214	gezelter	607	double aVel[3];
215			int i, j;
216	gezelter	468
217	gezelter	607
218			for (j=0; j<3; j++)
219			COMvel[j] = 0.0;
220
221	gezelter	468	mtot = 0.0;
222
223	tim	1113	for (i=0; i < myIntegrableObjects.size(); i++) {
224			if (myIntegrableObjects[i] != NULL) {
225	mmeineke	489
226	tim	1113	mass = myIntegrableObjects[i]->getMass();
227	gezelter	468	mtot += mass;
228
229	tim	1113	myIntegrableObjects[i]->getVel(aVel);
230	gezelter	607
231			for (j=0; j<3; j++)
232			COMvel[j] += aVel[j]*mass;
233
234	gezelter	468	}
235			}
236
237	gezelter	607	for (j=0; j<3; j++)
238			COMvel[j] /= mtot;
239	gezelter	468
240	gezelter	475	return mtot;
241
242	gezelter	468	}
243	tim	658
244			double Molecule::getTotalMass()
245			{
246	tim	1113
247	tim	658	double totalMass;
248
249			totalMass = 0;
250	tim	1113	for(int i =0; i < myIntegrableObjects.size(); i++){
251			totalMass += myIntegrableObjects[i]->getMass();
252	tim	658	}
253
254			return totalMass;
255			}