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 1492 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 2204 by gezelter, Fri Apr 15 22:04:00 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 "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;
68 < }
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;
21 <  
22 <  if( myAtoms != NULL ){
23 <    for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
24 <    delete[] myAtoms;
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;
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;
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;
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 <  
105 < }
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();
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 <  }
101 <
102 <  // Rigid Body forces and torques are done after the fortran force loop
103 <
104 < }
105 <
174 >    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
175 >      mass += sd->getMass();
176 >    }
177  
178 < double Molecule::getPotential( void ){
108 <  
109 <  int i;
110 <  double myPot = 0.0;
178 >    return mass;
179  
112  for(i=0; i<myRigidBodies.size(); i++) {
113    myRigidBodies[i]->updateAtoms();
180    }
115  
116  for(i=0; i<nBonds; i++){
117    myPot += myBonds[i]->get_potential();
118  }
181  
182 <  for(i=0; i<nBends; i++){
183 <    myPot += myBends[i]->get_potential();
184 <  }
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 <  for(i=0; i<nTorsions; i++){
125 <    myPot += myTorsions[i]->get_potential();
126 <  }
195 >    com /= totalMass;
196  
197 <  return myPot;
129 < }
130 <
131 < void Molecule::printMe( void ){
132 <  
133 <  int i;
134 <
135 <  for(i=0; i<nBonds; i++){
136 <    myBonds[i]->printMe();
197 >    return com;
198    }
199  
200 <  for(i=0; i<nBends; i++){
201 <    myBends[i]->printMe();
202 <  }
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  
143  for(i=0; i<nTorsions; i++){
144    myTorsions[i]->printMe();
208    }
209  
210 < }
211 <
212 < void Molecule::moveCOM(double delta[3]){
213 <  double aPos[3];
214 <  int i, j;
215 <
216 <  for(i=0; i<myIntegrableObjects.size(); i++) {
217 <    if(myIntegrableObjects[i] != NULL ) {
218 <      
219 <      myIntegrableObjects[i]->getPos( aPos );
220 <      
158 <      for (j=0; j< 3; j++)
159 <        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 <    }
174 < }
175 <
176 < void Molecule::atoms2rigidBodies( void ) {
177 <  int i;
178 <  for (i = 0; i < myRigidBodies.size(); i++) {
179 <    myRigidBodies[i]->calcForcesAndTorques();  
225 >    return velCom;
226    }
181 }
227  
228 < void Molecule::getCOM( double COM[3] ) {
228 >  double Molecule::getPotential() {
229  
230 <  double mass, mtot;
231 <  double aPos[3];
232 <  int i, j;
230 >    Bond* bond;
231 >    Bend* bend;
232 >    Torsion* torsion;
233 >    Molecule::BondIterator bondIter;;
234 >    Molecule::BendIterator  bendIter;
235 >    Molecule::TorsionIterator  torsionIter;
236  
237 <  for (j=0; j<3; j++)
190 <    COM[j] = 0.0;
237 >    double potential = 0.0;
238  
239 <  mtot   = 0.0;
239 >    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
240 >      potential += bond->getPotential();
241 >    }
242  
243 <  for (i=0; i < myIntegrableObjects.size(); i++) {
244 <    if (myIntegrableObjects[i] != NULL) {
243 >    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
244 >      potential += bend->getPotential();
245 >    }
246  
247 <      mass = myIntegrableObjects[i]->getMass();
248 <      mtot   += mass;
199 <      
200 <      myIntegrableObjects[i]->getPos( aPos );
201 <
202 <      for( j = 0; j < 3; j++)
203 <        COM[j] += aPos[j] * mass;
204 <
247 >    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
248 >      potential += torsion->getPotential();
249      }
206  }
250  
251 <  for (j = 0; j < 3; j++)
209 <    COM[j] /= mtot;
210 < }
251 >    return potential;
252  
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
235    }
253    }
254  
255 <  for (j=0; j<3; j++)
256 <    COMvel[j] /= mtot;
257 <
258 <  return mtot;
259 <
260 < }
261 <
262 < double Molecule::getTotalMass()
263 < {
264 <
265 <  double totalMass;
266 <  
250 <  totalMass = 0;
251 <  for(int i =0; i < myIntegrableObjects.size(); i++){
252 <    totalMass += myIntegrableObjects[i]->getMass();
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 <  return totalMass;
256 < }
269 > }//end namespace oopse

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines