OOPSE/libmdtools/Molecule.cpp

#include <cstdlib>


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

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


void Molecule::initialize( molInit &theInit ){

  nAtoms = theInit.nAtoms;
  nMembers = nAtoms;
  nBonds = theInit.nBonds;
  nBends = theInit.nBends;
  nTorsions = theInit.nTorsions;
  nExcludes = theInit.nExcludes;
  nOriented = theInit.nOriented;

  myAtoms = theInit.myAtoms;
  myBonds = theInit.myBonds;
  myBends = theInit.myBends;
  myTorsions = theInit.myTorsions;
  myExcludes = theInit.myExcludes;
    
}

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

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


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 x, y, z;
  int i;

  for(i=0; i<nAtoms; i++) {
    if(myAtoms[i] != NULL ) {

      x = myAtoms[i]->getX() + delta[0];
      y = myAtoms[i]->getY() + delta[1];
      z = myAtoms[i]->getZ() + delta[2];

      myAtoms[i]->setX(x);
      myAtoms[i]->setY(y);
      myAtoms[i]->setZ(z);
    }
  }
}

void Molecule::getCOM( double COM[3] ) {

  double mass, mtot;
  int i;

  COM[0] = 0.0;
  COM[1] = 0.0;
  COM[2] = 0.0;
  mtot   = 0.0;

  for (i=0; i < nAtoms; i++) {
    if (myAtoms[i] != NULL) {

      mass = myAtoms[i]->getMass();
      mtot   += mass;
      COM[0] += myAtoms[i]->getX() * mass;
      COM[1] += myAtoms[i]->getY() * mass;
      COM[2] += myAtoms[i]->getZ() * mass;

    }
  }

  COM[0] /= mtot;
  COM[1] /= mtot;
  COM[2] /= mtot;
 
}

double Molecule::getCOMvel( double COMvel[3] ) {

  double mass, mtot;
  int i;

  COMvel[0] = 0.0;
  COMvel[1] = 0.0;
  COMvel[2] = 0.0;
  mtot   = 0.0;

  for (i=0; i < nAtoms; i++) {
    if (myAtoms[i] != NULL) {

      mass = myAtoms[i]->getMass();
      mtot   += mass;
      COMvel[0] += myAtoms[i]->get_vx() * mass;
      COMvel[1] += myAtoms[i]->get_vy() * mass;
      COMvel[2] += myAtoms[i]->get_vz() * mass;

    }
  }

  COMvel[0] /= mtot;
  COMvel[1] /= mtot;
  COMvel[2] /= mtot;
 
  return mtot;

}

void Molecule::atomicRollCall(int* molMembership) {
  int i, which;

  for (i=0; i < nAtoms; i++) {
    if (myAtoms[i] != NULL) {
      which = myAtoms[i]->getIndex();
      molMembership[which] = globalIndex;
    }
  }
}
Revision:	489
Committed:	Thu Apr 10 20:08:56 2003 UTC (21 years, 3 months ago) by mmeineke
File size:	3625 byte(s)
Log Message:	added a globalIndex counter to Molecule
#	User	Rev	Content
1	mmeineke	427	#include <cstdlib>
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
19
20			Molecule::~Molecule( void ){
21			int i;
22
23			if( myAtoms != NULL ){
24			for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
25			delete[] myAtoms;
26			}
27
28			if( myBonds != NULL ){
29			for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i];
30			delete[] myBonds;
31			}
32
33			if( myBends != NULL ){
34			for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i];
35			delete[] myBends;
36			}
37
38			if( myTorsions != NULL ){
39			for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i];
40			delete[] myTorsions;
41			}
42
43			if( myExcludes != NULL ){
44			for(i=0; i<nExcludes; i++) if(myExcludes[i] != NULL ) delete myExcludes[i];
45			delete[] myExcludes;
46			}
47			}
48
49
50			void Molecule::initialize( molInit &theInit ){
51
52			nAtoms = theInit.nAtoms;
53			nMembers = nAtoms;
54			nBonds = theInit.nBonds;
55			nBends = theInit.nBends;
56			nTorsions = theInit.nTorsions;
57			nExcludes = theInit.nExcludes;
58			nOriented = theInit.nOriented;
59
60			myAtoms = theInit.myAtoms;
61			myBonds = theInit.myBonds;
62			myBends = theInit.myBends;
63			myTorsions = theInit.myTorsions;
64	mmeineke	428	myExcludes = theInit.myExcludes;
65	mmeineke	426
66			}
67
68			void Molecule::calcForces( void ){
69
70			int i;
71
72			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			}
84
85
86	mmeineke	428	double Molecule::getPotential( void ){
87	mmeineke	426
88			int i;
89			double myPot = 0.0;
90
91			for(i=0; i<nBonds; i++){
92			myPot += myBonds[i]->get_potential();
93			}
94
95			for(i=0; i<nBends; i++){
96			myPot += myBends[i]->get_potential();
97			}
98
99			for(i=0; i<nTorsions; i++){
100			myPot += myTorsions[i]->get_potential();
101			}
102
103			return myPot;
104			}
105	mmeineke	435
106			void Molecule::printMe( void ){
107
108			int i;
109
110			for(i=0; i<nBonds; i++){
111			myBonds[i]->printMe();
112			}
113
114			for(i=0; i<nBends; i++){
115			myBends[i]->printMe();
116			}
117
118			for(i=0; i<nTorsions; i++){
119			myTorsions[i]->printMe();
120			}
121			}
122	gezelter	446
123	mmeineke	449	void Molecule::moveCOM(double delta[3]){
124	gezelter	446	double x, y, z;
125			int i;
126
127			for(i=0; i<nAtoms; i++) {
128			if(myAtoms[i] != NULL ) {
129
130			x = myAtoms[i]->getX() + delta[0];
131			y = myAtoms[i]->getY() + delta[1];
132			z = myAtoms[i]->getZ() + delta[2];
133
134			myAtoms[i]->setX(x);
135			myAtoms[i]->setY(y);
136			myAtoms[i]->setZ(z);
137			}
138			}
139			}
140
141	mmeineke	449	void Molecule::getCOM( double COM[3] ) {
142	gezelter	446
143			double mass, mtot;
144			int i;
145
146			COM[0] = 0.0;
147			COM[1] = 0.0;
148			COM[2] = 0.0;
149			mtot = 0.0;
150
151			for (i=0; i < nAtoms; i++) {
152			if (myAtoms[i] != NULL) {
153
154			mass = myAtoms[i]->getMass();
155			mtot += mass;
156			COM[0] += myAtoms[i]->getX() * mass;
157			COM[1] += myAtoms[i]->getY() * mass;
158			COM[2] += myAtoms[i]->getZ() * mass;
159
160			}
161			}
162
163			COM[0] /= mtot;
164			COM[1] /= mtot;
165			COM[2] /= mtot;
166
167			}
168	gezelter	468
169	gezelter	475	double Molecule::getCOMvel( double COMvel[3] ) {
170	gezelter	468
171	gezelter	475	double mass, mtot;
172	gezelter	468	int i;
173
174			COMvel[0] = 0.0;
175			COMvel[1] = 0.0;
176			COMvel[2] = 0.0;
177			mtot = 0.0;
178
179			for (i=0; i < nAtoms; i++) {
180			if (myAtoms[i] != NULL) {
181	mmeineke	489
182	gezelter	468	mass = myAtoms[i]->getMass();
183			mtot += mass;
184			COMvel[0] += myAtoms[i]->get_vx() * mass;
185			COMvel[1] += myAtoms[i]->get_vy() * mass;
186			COMvel[2] += myAtoms[i]->get_vz() * mass;
187
188			}
189			}
190
191			COMvel[0] /= mtot;
192			COMvel[1] /= mtot;
193			COMvel[2] /= mtot;
194
195	gezelter	475	return mtot;
196
197	gezelter	468	}
198	gezelter	483
199			void Molecule::atomicRollCall(int* molMembership) {
200			int i, which;
201
202			for (i=0; i < nAtoms; i++) {
203			if (myAtoms[i] != NULL) {
204			which = myAtoms[i]->getIndex();
205	mmeineke	489	molMembership[which] = globalIndex;
206	gezelter	483	}
207			}
208			}