OOPSE/libmdtools/Molecule.cpp

#include <stdlib.h>


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


Molecule::Molecule( void ){

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

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 ){

  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();

}

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

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

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

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