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 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 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. Redistributions of source code must retain the above copyright
10 > *    notice, this list of conditions and the following disclaimer.
11 > *
12 > * 2. Redistributions in binary form must reproduce the above copyright
13 > *    notice, this list of conditions and the following disclaimer in the
14 > *    documentation and/or other materials provided with the
15 > *    distribution.
16 > *
17 > * This software is provided "AS IS," without a warranty of any
18 > * kind. All express or implied conditions, representations and
19 > * warranties, including any implied warranty of merchantability,
20 > * fitness for a particular purpose or non-infringement, are hereby
21 > * excluded.  The University of Notre Dame and its licensors shall not
22 > * be liable for any damages suffered by licensee as a result of
23 > * using, modifying or distributing the software or its
24 > * derivatives. In no event will the University of Notre Dame or its
25 > * licensors be liable for any lost revenue, profit or data, or for
26 > * direct, indirect, special, consequential, incidental or punitive
27 > * damages, however caused and regardless of the theory of liability,
28 > * arising out of the use of or inability to use software, even if the
29 > * University of Notre Dame has been advised of the possibility of
30 > * such damages.
31 > *
32 > * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 > * research, please cite the appropriate papers when you publish your
34 > * work.  Good starting points are:
35 > *                                                                      
36 > * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 > * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 > * [4]  Vardeman & Gezelter, in progress (2009).                        
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 <
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 OpenMD {
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(inversions_);
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 <  if( myBonds != NULL ){
78 <    for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i];
79 <    delete[] myBonds;
76 >  
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( myBends != NULL ){
84 <    for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i];
85 <    delete[] myBends;
82 >  
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( myTorsions != NULL ){
90 <    for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i];
91 <    delete[] myTorsions;
88 >  
89 >  void Molecule::addBend(Bend* bend) {
90 >    if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
91 >      bends_.push_back(bend);
92 >    }
93    }
41
42  for(cg = beginCutoffGroup(iter);  cg != NULL; cg = nextCutoffGroup(iter))
43    delete cg;
44  myCutoffGroups.clear();
94    
95 < }
96 <
97 <
98 < void Molecule::initialize( molInit &theInit ){
99 <
51 <  CutoffGroup* curCutoffGroup;
52 <  vector<CutoffGroup*>::iterator iterCutoff;
53 <  Atom* cutoffAtom;
54 <  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;
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 ){
82 <  
83 <  int i;
84 <  double com[3];
85 <
86 <  for(i=0; i<myRigidBodies.size(); i++) {
87 <    myRigidBodies[i]->updateAtoms();
95 >  void Molecule::addTorsion(Torsion* torsion) {
96 >    if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
97 >        torsions_.end()) {
98 >      torsions_.push_back(torsion);
99 >    }
100    }
101  
102 <  for(i=0; i<nBonds; i++){
103 <    myBonds[i]->calc_forces();
102 >  void Molecule::addInversion(Inversion* inversion) {
103 >    if (std::find(inversions_.begin(), inversions_.end(), inversion) ==
104 >        inversions_.end()) {
105 >      inversions_.push_back(inversion);
106 >    }
107    }
108 <
109 <  for(i=0; i<nBends; i++){
110 <    myBends[i]->calc_forces();
108 >  
109 >  void Molecule::addRigidBody(RigidBody *rb) {
110 >    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
111 >        rigidBodies_.end()) {
112 >      rigidBodies_.push_back(rb);
113 >    }
114    }
115 <
116 <  for(i=0; i<nTorsions; i++){
117 <    myTorsions[i]->calc_forces();
115 >  
116 >  void Molecule::addCutoffGroup(CutoffGroup* cp) {
117 >    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
118 >        cutoffGroups_.end()) {
119 >      cutoffGroups_.push_back(cp);
120 >    }    
121    }
101
102  // Rigid Body forces and torques are done after the fortran force loop
103
104 }
105
106
107 double Molecule::getPotential( void ){
122    
123 <  int i;
124 <  double myPot = 0.0;
125 <
126 <  for(i=0; i<myRigidBodies.size(); i++) {
127 <    myRigidBodies[i]->updateAtoms();
123 >  void Molecule::addConstraintPair(ConstraintPair* cp) {
124 >    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
125 >        constraintPairs_.end()) {
126 >      constraintPairs_.push_back(cp);
127 >    }    
128    }
129 <  
130 <  for(i=0; i<nBonds; i++){
131 <    myPot += myBonds[i]->get_potential();
129 >  
130 >  void Molecule::addConstraintElem(ConstraintElem* cp) {
131 >    if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
132 >        constraintElems_.end()) {
133 >      constraintElems_.push_back(cp);
134 >    }
135    }
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 ){
136    
137 <  int i;
137 >  void Molecule::complete() {
138 >    
139 >    std::set<Atom*> rigidAtoms;
140 >    RigidBody* rb;
141 >    std::vector<RigidBody*>::iterator rbIter;
142  
143 <  for(i=0; i<nBonds; i++){
144 <    myBonds[i]->printMe();
145 <  }
143 >    
144 >    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
145 >      rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
146 >    }
147 >    
148 >    Atom* atom;
149 >    AtomIterator ai;
150 >    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
151 >      
152 >      if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
153  
154 <  for(i=0; i<nBends; i++){
155 <    myBends[i]->printMe();
141 <  }
154 >        // If an atom does not belong to a rigid body, it is an
155 >        // integrable object
156  
157 <  for(i=0; i<nTorsions; i++){
158 <    myTorsions[i]->printMe();
157 >        integrableObjects_.push_back(*ai);
158 >      }
159 >    }
160 >    
161 >    //find all free atoms (which do not belong to rigid bodies)  
162 >    // performs the "difference" operation from set theory, the output
163 >    // range contains a copy of every element that is contained in
164 >    // [allAtoms.begin(), allAtoms.end()) and not contained in
165 >    // [rigidAtoms.begin(), rigidAtoms.end()).
166 >    //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
167 >    //                        std::back_inserter(integrableObjects_));
168 >
169 >    //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
170 >    //    //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
171 >    //    sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
172 >    //
173 >    //    painCave.isFatal = 1;
174 >    //    simError();        
175 >    //}
176 >    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
177 >      integrableObjects_.push_back(rb);
178 >    }
179 >    //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
180    }
181  
182 < }
183 <
184 < void Molecule::moveCOM(double delta[3]){
185 <  double aPos[3];
186 <  int i, j;
187 <
188 <  for(i=0; i<myIntegrableObjects.size(); i++) {
189 <    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 );
182 >  RealType Molecule::getMass() {
183 >    StuntDouble* sd;
184 >    std::vector<StuntDouble*>::iterator i;
185 >    RealType mass = 0.0;
186 >    
187 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
188 >           nextIntegrableObject(i)){
189 >      mass += sd->getMass();
190      }
191 +    
192 +    return mass;    
193    }
194  
195 <  for(i=0; i<myRigidBodies.size(); i++) {
196 <
197 <      myRigidBodies[i]->getPos( aPos );
198 <
199 <      for (j=0; j< 3; j++)
200 <        aPos[j] += delta[j];
201 <      
202 <      myRigidBodies[i]->setPos( aPos );
195 >  Vector3d Molecule::getCom() {
196 >    StuntDouble* sd;
197 >    std::vector<StuntDouble*>::iterator i;
198 >    Vector3d com;
199 >    RealType totalMass = 0;
200 >    RealType mass;
201 >    
202 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
203 >           nextIntegrableObject(i)){
204 >      mass = sd->getMass();
205 >      totalMass += mass;
206 >      com += sd->getPos() * mass;    
207      }
208 < }
208 >    
209 >    com /= totalMass;
210  
211 < void Molecule::atoms2rigidBodies( void ) {
177 <  int i;
178 <  for (i = 0; i < myRigidBodies.size(); i++) {
179 <    myRigidBodies[i]->calcForcesAndTorques();  
211 >    return com;
212    }
181 }
213  
214 < void Molecule::getCOM( double COM[3] ) {
214 >  void Molecule::moveCom(const Vector3d& delta) {
215 >    StuntDouble* sd;
216 >    std::vector<StuntDouble*>::iterator i;
217 >    
218 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
219 >           nextIntegrableObject(i)){
220 >      sd->setPos(sd->getPos() + delta);
221 >    }    
222 >  }
223  
224 <  double mass, mtot;
225 <  double aPos[3];
226 <  int i, j;
224 >  Vector3d Molecule::getComVel() {
225 >    StuntDouble* sd;
226 >    std::vector<StuntDouble*>::iterator i;
227 >    Vector3d velCom;
228 >    RealType totalMass = 0;
229 >    RealType mass;
230 >    
231 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
232 >           nextIntegrableObject(i)){
233 >      mass = sd->getMass();
234 >      totalMass += mass;
235 >      velCom += sd->getVel() * mass;    
236 >    }
237 >    
238 >    velCom /= totalMass;
239 >    
240 >    return velCom;
241 >  }
242  
243 <  for (j=0; j<3; j++)
190 <    COM[j] = 0.0;
243 >  RealType Molecule::getPotential() {
244  
245 <  mtot   = 0.0;
245 >    Bond* bond;
246 >    Bend* bend;
247 >    Torsion* torsion;
248 >    Inversion* inversion;
249 >    Molecule::BondIterator bondIter;;
250 >    Molecule::BendIterator  bendIter;
251 >    Molecule::TorsionIterator  torsionIter;
252 >    Molecule::InversionIterator  inversionIter;
253  
254 <  for (i=0; i < myIntegrableObjects.size(); i++) {
195 <    if (myIntegrableObjects[i] != NULL) {
254 >    RealType potential = 0.0;
255  
256 <      mass = myIntegrableObjects[i]->getMass();
257 <      mtot   += mass;
258 <      
200 <      myIntegrableObjects[i]->getPos( aPos );
256 >    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
257 >      potential += bond->getPotential();
258 >    }
259  
260 <      for( j = 0; j < 3; j++)
261 <        COM[j] += aPos[j] * mass;
260 >    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
261 >      potential += bend->getPotential();
262 >    }
263  
264 +    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
265 +           nextTorsion(torsionIter)) {
266 +      potential += torsion->getPotential();
267      }
268 +    
269 +    for (inversion = beginInversion(inversionIter); torsion != NULL;
270 +         inversion =  nextInversion(inversionIter)) {
271 +      potential += inversion->getPotential();
272 +    }
273 +    
274 +    return potential;
275 +    
276    }
277 +  
278 +  void Molecule::addProperty(GenericData* genData) {
279 +    properties_.addProperty(genData);  
280 +  }
281  
282 <  for (j = 0; j < 3; j++)
283 <    COM[j] /= mtot;
284 < }
282 >  void Molecule::removeProperty(const std::string& propName) {
283 >    properties_.removeProperty(propName);  
284 >  }
285  
286 < double Molecule::getCOMvel( double COMvel[3] ) {
286 >  void Molecule::clearProperties() {
287 >    properties_.clearProperties();
288 >  }
289  
290 <  double mass, mtot;
291 <  double aVel[3];
292 <  int i, j;
290 >  std::vector<std::string> Molecule::getPropertyNames() {
291 >    return properties_.getPropertyNames();  
292 >  }
293 >      
294 >  std::vector<GenericData*> Molecule::getProperties() {
295 >    return properties_.getProperties();
296 >  }
297  
298 <
299 <  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 <    }
298 >  GenericData* Molecule::getPropertyByName(const std::string& propName) {
299 >    return properties_.getPropertyByName(propName);
300    }
301  
238  for (j=0; j<3; j++)
239    COMvel[j] /= mtot;
240
241  return mtot;
302  
243 }
303  
245 double Molecule::getTotalMass()
246 {
304  
305 <  double totalMass;
306 <  
307 <  totalMass = 0;
308 <  for(int i =0; i < myIntegrableObjects.size(); i++){
309 <    totalMass += myIntegrableObjects[i]->getMass();
305 >  std::ostream& operator <<(std::ostream& o, Molecule& mol) {
306 >    o << std::endl;
307 >    o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
308 >    o << mol.getNAtoms() << " atoms" << std::endl;
309 >    o << mol.getNBonds() << " bonds" << std::endl;
310 >    o << mol.getNBends() << " bends" << std::endl;
311 >    o << mol.getNTorsions() << " torsions" << std::endl;
312 >    o << mol.getNInversions() << " inversions" << std::endl;
313 >    o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
314 >    o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
315 >    o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
316 >    o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
317 >    return o;
318    }
319 <
320 <  return totalMass;
256 < }
319 >  
320 > }//end namespace OpenMD

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines