ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/primitives/Molecule.cpp
(Generate patch)

Comparing trunk/src/primitives/Molecule.cpp (file contents):
Revision 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 1211 by gezelter, Wed Jan 23 16:38:22 2008 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 <
57 <
58 < Molecule::Molecule( void ){
59 <
11 <  myAtoms = NULL;
12 <  myBonds = NULL;
13 <  myBends = NULL;
14 <  myTorsions = NULL;
15 < }
16 <
17 < Molecule::~Molecule( void ){
18 <  int i;
19 <  CutoffGroup* cg;
20 <  vector<CutoffGroup*>::iterator iter;
56 > namespace oopse {
57 >  Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
58 >    : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
59 >  }
60    
61 <  if( myAtoms != NULL ){
62 <    for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
63 <    delete[] myAtoms;
61 >  Molecule::~Molecule() {
62 >    
63 >    MemoryUtils::deletePointers(atoms_);
64 >    MemoryUtils::deletePointers(bonds_);
65 >    MemoryUtils::deletePointers(bends_);
66 >    MemoryUtils::deletePointers(torsions_);
67 >    MemoryUtils::deletePointers(rigidBodies_);
68 >    MemoryUtils::deletePointers(cutoffGroups_);
69 >    MemoryUtils::deletePointers(constraintPairs_);
70 >    MemoryUtils::deletePointers(constraintElems_);
71 >    // integrableObjects_ don't own the objects
72 >    integrableObjects_.clear();
73 >    
74    }
75 <
76 <  if( myBonds != NULL ){
77 <    for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i];
78 <    delete[] myBonds;
75 >  
76 >  void Molecule::addAtom(Atom* atom) {
77 >    if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
78 >      atoms_.push_back(atom);
79 >    }
80    }
81 <
82 <  if( myBends != NULL ){
83 <    for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i];
84 <    delete[] myBends;
81 >  
82 >  void Molecule::addBond(Bond* bond) {
83 >    if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
84 >      bonds_.push_back(bond);
85 >    }
86    }
87 <
88 <  if( myTorsions != NULL ){
89 <    for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i];
90 <    delete[] myTorsions;
87 >  
88 >  void Molecule::addBend(Bend* bend) {
89 >    if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
90 >      bends_.push_back(bend);
91 >    }
92    }
41
42  for(cg = beginCutoffGroup(iter);  cg != NULL; cg = nextCutoffGroup(iter))
43    delete cg;
44  myCutoffGroups.clear();
93    
94 < }
95 <
96 <
97 < void Molecule::initialize( molInit &theInit ){
98 <
99 <  CutoffGroup* curCutoffGroup;
52 <  vector<CutoffGroup*>::iterator iterCutoff;
53 <  Atom* cutoffAtom;
54 <  vector<Atom*>::iterator iterAtom;
55 <  int atomIndex;
94 >  void Molecule::addTorsion(Torsion* torsion) {
95 >    if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
96 >        torsions_.end()) {
97 >      torsions_.push_back(torsion);
98 >    }
99 >  }
100    
101 <  nAtoms = theInit.nAtoms;
102 <  nMembers = nAtoms;
103 <  nBonds = theInit.nBonds;
104 <  nBends = theInit.nBends;
105 <  nTorsions = theInit.nTorsions;
106 <  nRigidBodies = theInit.nRigidBodies;
63 <  nOriented = theInit.nOriented;
64 <
65 <  myAtoms = theInit.myAtoms;
66 <  myBonds = theInit.myBonds;
67 <  myBends = theInit.myBends;
68 <  myTorsions = theInit.myTorsions;
69 <  myRigidBodies = theInit.myRigidBodies;
70 <
71 <  myIntegrableObjects = theInit.myIntegrableObjects;
72 <
73 <  for (int i = 0; i < myRigidBodies.size(); i++)
74 <      myRigidBodies[i]->calcRefCoords();
75 <
76 <  myCutoffGroups = theInit.myCutoffGroups;
77 <  nCutoffGroups = myCutoffGroups.size();
78 <
79 < }
80 <
81 < void Molecule::calcForces( void ){
101 >  void Molecule::addRigidBody(RigidBody *rb) {
102 >    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
103 >        rigidBodies_.end()) {
104 >      rigidBodies_.push_back(rb);
105 >    }
106 >  }
107    
108 <  int i;
109 <  double com[3];
110 <
111 <  for(i=0; i<myRigidBodies.size(); i++) {
112 <    myRigidBodies[i]->updateAtoms();
108 >  void Molecule::addCutoffGroup(CutoffGroup* cp) {
109 >    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
110 >        cutoffGroups_.end()) {
111 >      cutoffGroups_.push_back(cp);
112 >    }    
113    }
114 <
115 <  for(i=0; i<nBonds; i++){
116 <    myBonds[i]->calc_forces();
114 >  
115 >  void Molecule::addConstraintPair(ConstraintPair* cp) {
116 >    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
117 >        constraintPairs_.end()) {
118 >      constraintPairs_.push_back(cp);
119 >    }    
120    }
93
94  for(i=0; i<nBends; i++){
95    myBends[i]->calc_forces();
96  }
97
98  for(i=0; i<nTorsions; i++){
99    myTorsions[i]->calc_forces();
100  }
101
102  // Rigid Body forces and torques are done after the fortran force loop
103
104 }
105
106
107 double Molecule::getPotential( void ){
121    
122 <  int i;
123 <  double myPot = 0.0;
124 <
125 <  for(i=0; i<myRigidBodies.size(); i++) {
126 <    myRigidBodies[i]->updateAtoms();
122 >  void Molecule::addConstraintElem(ConstraintElem* cp) {
123 >    if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
124 >        constraintElems_.end()) {
125 >      constraintElems_.push_back(cp);
126 >    }
127    }
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;
129 }
130
131 void Molecule::printMe( void ){
128    
129 <  int i;
129 >  void Molecule::complete() {
130 >    
131 >    std::set<Atom*> rigidAtoms;
132 >    RigidBody* rb;
133 >    std::vector<RigidBody*>::iterator rbIter;
134  
135 <  for(i=0; i<nBonds; i++){
136 <    myBonds[i]->printMe();
137 <  }
135 >    
136 >    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
137 >      rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
138 >    }
139 >    
140 >    Atom* atom;
141 >    AtomIterator ai;
142 >    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
143 >      
144 >      if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
145  
146 <  for(i=0; i<nBends; i++){
147 <    myBends[i]->printMe();
141 <  }
146 >        // If an atom does not belong to a rigid body, it is an
147 >        // integrable object
148  
149 <  for(i=0; i<nTorsions; i++){
150 <    myTorsions[i]->printMe();
149 >        integrableObjects_.push_back(*ai);
150 >      }
151 >    }
152 >    
153 >    //find all free atoms (which do not belong to rigid bodies)  
154 >    // performs the "difference" operation from set theory, the output
155 >    // range contains a copy of every element that is contained in
156 >    // [allAtoms.begin(), allAtoms.end()) and not contained in
157 >    // [rigidAtoms.begin(), rigidAtoms.end()).
158 >    //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
159 >    //                        std::back_inserter(integrableObjects_));
160 >
161 >    //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
162 >    //    //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
163 >    //    sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
164 >    //
165 >    //    painCave.isFatal = 1;
166 >    //    simError();        
167 >    //}
168 >    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
169 >      integrableObjects_.push_back(rb);
170 >    }
171 >    //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
172    }
173  
174 < }
175 <
176 < void Molecule::moveCOM(double delta[3]){
177 <  double aPos[3];
178 <  int i, j;
179 <
180 <  for(i=0; i<myIntegrableObjects.size(); i++) {
181 <    if(myIntegrableObjects[i] != NULL ) {
155 <      
156 <      myIntegrableObjects[i]->getPos( aPos );
157 <      
158 <      for (j=0; j< 3; j++)
159 <        aPos[j] += delta[j];
160 <
161 <      myIntegrableObjects[i]->setPos( aPos );
174 >  RealType Molecule::getMass() {
175 >    StuntDouble* sd;
176 >    std::vector<StuntDouble*>::iterator i;
177 >    RealType mass = 0.0;
178 >    
179 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
180 >           nextIntegrableObject(i)){
181 >      mass += sd->getMass();
182      }
183 +    
184 +    return mass;    
185    }
186  
187 <  for(i=0; i<myRigidBodies.size(); i++) {
188 <
189 <      myRigidBodies[i]->getPos( aPos );
190 <
191 <      for (j=0; j< 3; j++)
192 <        aPos[j] += delta[j];
193 <      
194 <      myRigidBodies[i]->setPos( aPos );
187 >  Vector3d Molecule::getCom() {
188 >    StuntDouble* sd;
189 >    std::vector<StuntDouble*>::iterator i;
190 >    Vector3d com;
191 >    RealType totalMass = 0;
192 >    RealType mass;
193 >    
194 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
195 >           nextIntegrableObject(i)){
196 >      mass = sd->getMass();
197 >      totalMass += mass;
198 >      com += sd->getPos() * mass;    
199      }
200 < }
200 >    
201 >    com /= totalMass;
202  
203 < void Molecule::atoms2rigidBodies( void ) {
177 <  int i;
178 <  for (i = 0; i < myRigidBodies.size(); i++) {
179 <    myRigidBodies[i]->calcForcesAndTorques();  
203 >    return com;
204    }
181 }
205  
206 < void Molecule::getCOM( double COM[3] ) {
206 >  void Molecule::moveCom(const Vector3d& delta) {
207 >    StuntDouble* sd;
208 >    std::vector<StuntDouble*>::iterator i;
209 >    
210 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
211 >           nextIntegrableObject(i)){
212 >      sd->setPos(sd->getPos() + delta);
213 >    }    
214 >  }
215  
216 <  double mass, mtot;
217 <  double aPos[3];
218 <  int i, j;
219 <
220 <  for (j=0; j<3; j++)
221 <    COM[j] = 0.0;
222 <
223 <  mtot   = 0.0;
224 <
225 <  for (i=0; i < myIntegrableObjects.size(); i++) {
226 <    if (myIntegrableObjects[i] != NULL) {
227 <
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 <
216 >  Vector3d Molecule::getComVel() {
217 >    StuntDouble* sd;
218 >    std::vector<StuntDouble*>::iterator i;
219 >    Vector3d velCom;
220 >    RealType totalMass = 0;
221 >    RealType mass;
222 >    
223 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
224 >           nextIntegrableObject(i)){
225 >      mass = sd->getMass();
226 >      totalMass += mass;
227 >      velCom += sd->getVel() * mass;    
228      }
229 +    
230 +    velCom /= totalMass;
231 +    
232 +    return velCom;
233    }
234  
235 <  for (j = 0; j < 3; j++)
209 <    COM[j] /= mtot;
210 < }
235 >  RealType Molecule::getPotential() {
236  
237 < double Molecule::getCOMvel( double COMvel[3] ) {
237 >    Bond* bond;
238 >    Bend* bend;
239 >    Torsion* torsion;
240 >    Molecule::BondIterator bondIter;;
241 >    Molecule::BendIterator  bendIter;
242 >    Molecule::TorsionIterator  torsionIter;
243  
244 <  double mass, mtot;
215 <  double aVel[3];
216 <  int i, j;
244 >    RealType potential = 0.0;
245  
246 +    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
247 +      potential += bond->getPotential();
248 +    }
249  
250 <  for (j=0; j<3; j++)
251 <    COMvel[j] = 0.0;
250 >    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
251 >      potential += bend->getPotential();
252 >    }
253  
254 <  mtot   = 0.0;
255 <
256 <  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 <
254 >    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
255 >           nextTorsion(torsionIter)) {
256 >      potential += torsion->getPotential();
257      }
258 +    
259 +    return potential;
260 +    
261    }
237
238  for (j=0; j<3; j++)
239    COMvel[j] /= mtot;
240
241  return mtot;
242
243 }
244
245 double Molecule::getTotalMass()
246 {
247
248  double totalMass;
262    
263 <  totalMass = 0;
264 <  for(int i =0; i < myIntegrableObjects.size(); i++){
265 <    totalMass += myIntegrableObjects[i]->getMass();
263 >  std::ostream& operator <<(std::ostream& o, Molecule& mol) {
264 >    o << std::endl;
265 >    o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
266 >    o << mol.getNAtoms() << " atoms" << std::endl;
267 >    o << mol.getNBonds() << " bonds" << std::endl;
268 >    o << mol.getNBends() << " bends" << std::endl;
269 >    o << mol.getNTorsions() << " torsions" << std::endl;
270 >    o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
271 >    o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
272 >    o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
273 >    o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
274 >    return o;
275    }
276 <
277 <  return totalMass;
256 < }
276 >  
277 > }//end namespace oopse

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines