OOPSE/libmdtools/Molecule.cpp

#include <stdlib.h>


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


Molecule::Molecule( void ){

  myAtoms = NULL;
  myBonds = NULL;
  myBends = NULL;
  myTorsions = NULL;
  myRigidBodies = 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;
  }

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


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

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

  for(i=0; i<nRigidBodies; 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<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<nAtoms; i++) {
    if(myAtoms[i] != NULL ) {
      
      myAtoms[i]->getPos( aPos );
      
      for (j=0; j< 3; j++) 
        aPos[j] += delta[j];

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

  for(i=0; i<nRigidBodies; i++) {

    if (myRigidBodies[i] != NULL) {
      
      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 < nRigidBodies; i++) {
    if (myRigidBodies[i] != NULL) {
      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 < nAtoms; i++) {
    if (myAtoms[i] != NULL) {

      mass = myAtoms[i]->getMass();
      mtot   += mass;
      
      myAtoms[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 < nAtoms; i++) {
    if (myAtoms[i] != NULL) {

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

      myAtoms[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()
{
  int natoms;
  Atom** atoms;
  double totalMass;
  
  natoms = getNAtoms();
  atoms = getMyAtoms();
  totalMass = 0;
  for(int i =0; i < natoms; i++){
    totalMass += atoms[i]->getMass();
  }

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