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

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

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

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