ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Molecule.cpp
Revision: 446
Committed: Thu Apr 3 20:19:50 2003 UTC (21 years, 3 months ago) by gezelter
File size: 2887 byte(s)
Log Message:
Starting work on NPT

File Contents

# 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     void Molecule::moveCOM(double* delta){
124     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     double* Molecule::getCOM() {
142    
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     return COM;
168     }