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

  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;
  nCutoffGroups = theInit.nCutoffGroups;

  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;

  //creat a global index set of atoms which belong to cutoff group. 
  //it will be use to speed up the query whether an atom belongs to cutoff group or not
  cutoffAtomSet.clear();

  for(curCutoffGroup = beginCutoffGroup(iterCutoff); curCutoffGroup != NULL;
                                                  curCutoffGroup = nextCutoffGroup(iterCutoff)){

      for(cutoffAtom = curCutoffGroup->beginAtom(iterAtom); cutoffAtom != NULL; 
                                           cutoffAtom = curCutoffGroup->beginAtom(iterAtom)){
#ifdef IS_MPI
        atomIndex = cutoffAtom->getGlobalIndex();
#else
        atomIndex = cutoffAtom->getIndex();
#endif
        cutoffAtomSet.insert(atomIndex);
     }// end for(cutoffAtom)    
  }//end for(curCutoffGroup)
  
}

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