ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/primitives/Molecule.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/primitives/Molecule.cpp (file contents):
Revision 1490 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
Revision 2082 by tim, Mon Mar 7 22:39:33 2005 UTC

# Line 1 | Line 1
1 < #include <stdlib.h>
1 > /*
2 > * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 > *
4 > * The University of Notre Dame grants you ("Licensee") a
5 > * non-exclusive, royalty free, license to use, modify and
6 > * redistribute this software in source and binary code form, provided
7 > * that the following conditions are met:
8 > *
9 > * 1. Acknowledgement of the program authors must be made in any
10 > *    publication of scientific results based in part on use of the
11 > *    program.  An acceptable form of acknowledgement is citation of
12 > *    the article in which the program was described (Matthew
13 > *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 > *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 > *    Parallel Simulation Engine for Molecular Dynamics,"
16 > *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 > *
18 > * 2. Redistributions of source code must retain the above copyright
19 > *    notice, this list of conditions and the following disclaimer.
20 > *
21 > * 3. Redistributions in binary form must reproduce the above copyright
22 > *    notice, this list of conditions and the following disclaimer in the
23 > *    documentation and/or other materials provided with the
24 > *    distribution.
25 > *
26 > * This software is provided "AS IS," without a warranty of any
27 > * kind. All express or implied conditions, representations and
28 > * warranties, including any implied warranty of merchantability,
29 > * fitness for a particular purpose or non-infringement, are hereby
30 > * excluded.  The University of Notre Dame and its licensors shall not
31 > * be liable for any damages suffered by licensee as a result of
32 > * using, modifying or distributing the software or its
33 > * derivatives. In no event will the University of Notre Dame or its
34 > * licensors be liable for any lost revenue, profit or data, or for
35 > * direct, indirect, special, consequential, incidental or punitive
36 > * damages, however caused and regardless of the theory of liability,
37 > * arising out of the use of or inability to use software, even if the
38 > * University of Notre Dame has been advised of the possibility of
39 > * such damages.
40 > */
41 >
42 > /**
43 > * @file Molecule.cpp
44 > * @author    tlin
45 > * @date  10/28/2004
46 > * @version 1.0
47 > */
48  
49 + #include <algorithm>
50 + #include <set>
51  
52 < #include "Molecule.hpp"
53 < #include "simError.h"
52 > #include "primitives/Molecule.hpp"
53 > #include "utils/MemoryUtils.hpp"
54 > #include "utils/simError.h"
55  
56 + namespace oopse {
57 + Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
58 +    : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
59  
60 + }
61  
62 < Molecule::Molecule( void ){
62 > Molecule::~Molecule() {
63  
64 <  myAtoms = NULL;
65 <  myBonds = NULL;
66 <  myBends = NULL;
67 <  myTorsions = NULL;
64 >    MemoryUtils::deletePointers(atoms_);
65 >    MemoryUtils::deletePointers(bonds_);
66 >    MemoryUtils::deletePointers(bends_);
67 >    MemoryUtils::deletePointers(torsions_);
68 >    MemoryUtils::deletePointers(rigidBodies_);
69 >    MemoryUtils::deletePointers(cutoffGroups_);
70 >    MemoryUtils::deletePointers(constraintPairs_);
71 >    MemoryUtils::deletePointers(constraintElems_);
72 >    //integrableObjects_ don't own the objects
73 >    integrableObjects_.clear();
74 >    
75   }
76  
77 < Molecule::~Molecule( void ){
78 <  int i;
79 <  CutoffGroup* cg;
80 <  vector<CutoffGroup*>::iterator iter;
81 <  
22 <  if( myAtoms != NULL ){
23 <    for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
24 <    delete[] myAtoms;
25 <  }
77 > void Molecule::addAtom(Atom* atom) {
78 >    if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
79 >        atoms_.push_back(atom);
80 >    }
81 > }
82  
83 <  if( myBonds != NULL ){
84 <    for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i];
85 <    delete[] myBonds;
86 <  }
83 > void Molecule::addBond(Bond* bond) {
84 >    if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
85 >        bonds_.push_back(bond);
86 >    }
87 > }
88  
89 <  if( myBends != NULL ){
90 <    for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i];
91 <    delete[] myBends;
92 <  }
89 > void Molecule::addBend(Bend* bend) {
90 >    if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
91 >        bends_.push_back(bend);
92 >    }
93 > }
94  
95 <  if( myTorsions != NULL ){
96 <    for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i];
97 <    delete[] myTorsions;
98 <  }
95 > void Molecule::addTorsion(Torsion* torsion) {
96 >    if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
97 >        torsions_.push_back(torsion);
98 >    }
99 > }
100  
101 <  for(cg = beginCutoffGroup(iter);  cg != NULL; cg = nextCutoffGroup(iter))
102 <    delete cg;
103 <  myCutoffGroups.clear();
104 <  
101 > void Molecule::addRigidBody(RigidBody *rb) {
102 >    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
103 >        rigidBodies_.push_back(rb);
104 >    }
105   }
106  
107 + void Molecule::addCutoffGroup(CutoffGroup* cp) {
108 +    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
109 +        cutoffGroups_.push_back(cp);
110 +    }
111  
112 < void Molecule::initialize( molInit &theInit ){
112 > }
113  
114 <  CutoffGroup* curCutoffGroup;
115 <  vector<CutoffGroup*>::iterator iterCutoff;
116 <  Atom* cutoffAtom;
117 <  vector<Atom*>::iterator iterAtom;
55 <  int atomIndex;
56 <  
57 <  nAtoms = theInit.nAtoms;
58 <  nMembers = nAtoms;
59 <  nBonds = theInit.nBonds;
60 <  nBends = theInit.nBends;
61 <  nTorsions = theInit.nTorsions;
62 <  nRigidBodies = theInit.nRigidBodies;
63 <  nOriented = theInit.nOriented;
114 > void Molecule::addConstraintPair(ConstraintPair* cp) {
115 >    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) {
116 >        constraintPairs_.push_back(cp);
117 >    }
118  
119 <  myAtoms = theInit.myAtoms;
66 <  myBonds = theInit.myBonds;
67 <  myBends = theInit.myBends;
68 <  myTorsions = theInit.myTorsions;
69 <  myRigidBodies = theInit.myRigidBodies;
119 > }
120  
121 <  myIntegrableObjects = theInit.myIntegrableObjects;
121 > void Molecule::addConstraintElem(ConstraintElem* cp) {
122 >    if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == constraintElems_.end()) {
123 >        constraintElems_.push_back(cp);
124 >    }
125  
126 <  for (int i = 0; i < myRigidBodies.size(); i++)
74 <      myRigidBodies[i]->calcRefCoords();
126 > }
127  
128 <  myCutoffGroups = theInit.myCutoffGroups;
129 <  nCutoffGroups = myCutoffGroups.size();
128 > void Molecule::complete() {
129 >    
130 >    std::set<Atom*> rigidAtoms;
131 >    RigidBody* rb;
132 >    std::vector<RigidBody*>::iterator rbIter;
133  
134 < }
134 >    
135 >    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
136 >        rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
137 >    }
138  
139 < void Molecule::calcForces( void ){
140 <  
141 <  int i;
142 <  double com[3];
139 >    Atom* atom;
140 >    AtomIterator ai;
141 >    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
142 >  
143 >        if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
144 >            //if an atom does not belong to a rigid body, it is an integrable object
145 >            integrableObjects_.push_back(*ai);
146 >        }
147 >    }
148  
149 <  for(i=0; i<myRigidBodies.size(); i++) {
150 <    myRigidBodies[i]->updateAtoms();
151 <  }
149 >    //find all free atoms (which do not belong to rigid bodies)  
150 >    //performs the "difference" operation from set theory,  the output range contains a copy of every
151 >    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
152 >    //[rigidAtoms.begin(), rigidAtoms.end()).
153 >    //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
154 >    //                        std::back_inserter(integrableObjects_));
155  
156 <  for(i=0; i<nBonds; i++){
157 <    myBonds[i]->calc_forces();
158 <  }
156 >    //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
157 >    //    //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
158 >    //    sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
159 >    //
160 >    //    painCave.isFatal = 1;
161 >    //    simError();        
162 >    //}
163 >    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
164 >      integrableObjects_.push_back(rb);
165 >    }
166 >    //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
167 > }
168  
169 <  for(i=0; i<nBends; i++){
170 <    myBends[i]->calc_forces();
171 <  }
169 > double Molecule::getMass() {
170 >    StuntDouble* sd;
171 >    std::vector<StuntDouble*>::iterator i;
172 >    double mass = 0.0;
173  
174 <  for(i=0; i<nTorsions; i++){
175 <    myTorsions[i]->calc_forces();
176 <  }
174 >    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
175 >        mass += sd->getMass();
176 >    }
177  
178 <  // Rigid Body forces and torques are done after the fortran force loop
178 >    return mass;
179  
180   }
181  
182 + Vector3d Molecule::getCom() {
183 +    StuntDouble* sd;
184 +    std::vector<StuntDouble*>::iterator i;
185 +    Vector3d com;
186 +    double totalMass = 0;
187 +    double mass;
188 +    
189 +    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
190 +        mass = sd->getMass();
191 +        totalMass += mass;
192 +        com += sd->getPos() * mass;    
193 +    }
194  
195 < double Molecule::getPotential( void ){
108 <  
109 <  int i;
110 <  double myPot = 0.0;
195 >    com /= totalMass;
196  
197 <  for(i=0; i<myRigidBodies.size(); i++) {
113 <    myRigidBodies[i]->updateAtoms();
114 <  }
115 <  
116 <  for(i=0; i<nBonds; i++){
117 <    myPot += myBonds[i]->get_potential();
118 <  }
119 <
120 <  for(i=0; i<nBends; i++){
121 <    myPot += myBends[i]->get_potential();
122 <  }
123 <
124 <  for(i=0; i<nTorsions; i++){
125 <    myPot += myTorsions[i]->get_potential();
126 <  }
127 <
128 <  return myPot;
197 >    return com;
198   }
199  
200 < void Molecule::printMe( void ){
201 <  
202 <  int i;
200 > void Molecule::moveCom(const Vector3d& delta) {
201 >    StuntDouble* sd;
202 >    std::vector<StuntDouble*>::iterator i;
203 >    
204 >    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
205 >        sd->setPos(sd->getPos() + delta);
206 >    }
207  
135  for(i=0; i<nBonds; i++){
136    myBonds[i]->printMe();
137  }
138
139  for(i=0; i<nBends; i++){
140    myBends[i]->printMe();
141  }
142
143  for(i=0; i<nTorsions; i++){
144    myTorsions[i]->printMe();
145  }
146
208   }
209  
210 < void Molecule::moveCOM(double delta[3]){
211 <  double aPos[3];
212 <  int i, j;
213 <
214 <  for(i=0; i<myIntegrableObjects.size(); i++) {
215 <    if(myIntegrableObjects[i] != NULL ) {
216 <      
217 <      myIntegrableObjects[i]->getPos( aPos );
218 <      
219 <      for (j=0; j< 3; j++)
220 <        aPos[j] += delta[j];
160 <
161 <      myIntegrableObjects[i]->setPos( aPos );
210 > Vector3d Molecule::getComVel() {
211 >    StuntDouble* sd;
212 >    std::vector<StuntDouble*>::iterator i;
213 >    Vector3d velCom;
214 >    double totalMass = 0;
215 >    double mass;
216 >    
217 >    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
218 >        mass = sd->getMass();
219 >        totalMass += mass;
220 >        velCom += sd->getVel() * mass;    
221      }
163  }
222  
223 <  for(i=0; i<myRigidBodies.size(); i++) {
223 >    velCom /= totalMass;
224  
225 <      myRigidBodies[i]->getPos( aPos );
168 <
169 <      for (j=0; j< 3; j++)
170 <        aPos[j] += delta[j];
171 <      
172 <      myRigidBodies[i]->setPos( aPos );
173 <    }
225 >    return velCom;
226   }
227  
228 < void Molecule::atoms2rigidBodies( void ) {
177 <  int i;
178 <  for (i = 0; i < myRigidBodies.size(); i++) {
179 <    myRigidBodies[i]->calcForcesAndTorques();  
180 <  }
181 < }
228 > double Molecule::getPotential() {
229  
230 < void Molecule::getCOM( double COM[3] ) {
230 >    Bond* bond;
231 >    Bend* bend;
232 >    Torsion* torsion;
233 >    Molecule::BondIterator bondIter;;
234 >    Molecule::BendIterator  bendIter;
235 >    Molecule::TorsionIterator  torsionIter;
236  
237 <  double mass, mtot;
186 <  double aPos[3];
187 <  int i, j;
237 >    double potential = 0.0;
238  
239 <  for (j=0; j<3; j++)
240 <    COM[j] = 0.0;
239 >    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
240 >        potential += bond->getPotential();
241 >    }
242  
243 <  mtot   = 0.0;
244 <
194 <  for (i=0; i < myIntegrableObjects.size(); i++) {
195 <    if (myIntegrableObjects[i] != NULL) {
196 <
197 <      mass = myIntegrableObjects[i]->getMass();
198 <      mtot   += mass;
199 <      
200 <      myIntegrableObjects[i]->getPos( aPos );
201 <
202 <      for( j = 0; j < 3; j++)
203 <        COM[j] += aPos[j] * mass;
204 <
243 >    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
244 >        potential += bend->getPotential();
245      }
206  }
246  
247 <  for (j = 0; j < 3; j++)
248 <    COM[j] /= mtot;
210 < }
211 <
212 < double Molecule::getCOMvel( double COMvel[3] ) {
213 <
214 <  double mass, mtot;
215 <  double aVel[3];
216 <  int i, j;
217 <
218 <
219 <  for (j=0; j<3; j++)
220 <    COMvel[j] = 0.0;
221 <
222 <  mtot   = 0.0;
223 <
224 <  for (i=0; i < myIntegrableObjects.size(); i++) {
225 <    if (myIntegrableObjects[i] != NULL) {
226 <
227 <      mass = myIntegrableObjects[i]->getMass();
228 <      mtot   += mass;
229 <
230 <      myIntegrableObjects[i]->getVel(aVel);
231 <
232 <      for (j=0; j<3; j++)
233 <        COMvel[j] += aVel[j]*mass;
234 <
247 >    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
248 >        potential += torsion->getPotential();
249      }
236  }
250  
251 <  for (j=0; j<3; j++)
239 <    COMvel[j] /= mtot;
240 <
241 <  return mtot;
251 >    return potential;
252  
253   }
254  
255 < double Molecule::getTotalMass()
256 < {
257 <
258 <  double totalMass;
259 <  
260 <  totalMass = 0;
261 <  for(int i =0; i < myIntegrableObjects.size(); i++){
262 <    totalMass += myIntegrableObjects[i]->getMass();
263 <  }
264 <
265 <  return totalMass;
255 > std::ostream& operator <<(std::ostream& o, Molecule& mol) {
256 >    o << std::endl;
257 >    o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
258 >    o << mol.getNAtoms() << " atoms" << std::endl;
259 >    o << mol.getNBonds() << " bonds" << std::endl;
260 >    o << mol.getNBends() << " bends" << std::endl;
261 >    o << mol.getNTorsions() << " torsions" << std::endl;
262 >    o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
263 >    o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
264 >    o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
265 >    o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
266 >    return o;
267   }
268 +
269 + }//end namespace oopse

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines