ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/brains/ForceManager.cpp
(Generate patch)

Comparing trunk/src/brains/ForceManager.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 664 by chuckv, Wed Oct 12 21:57:16 2005 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 39 | Line 39
39   * such damages.
40   */
41  
42 < /**
43 <  * @file ForceManager.cpp
44 <  * @author tlin
45 <  * @date 11/09/2004
46 <  * @time 10:39am
47 <  * @version 1.0
48 <  */
42 > /**
43 > * @file ForceManager.cpp
44 > * @author tlin
45 > * @date 11/09/2004
46 > * @time 10:39am
47 > * @version 1.0
48 > */
49  
50   #include "brains/ForceManager.hpp"
51   #include "primitives/Molecule.hpp"
52   #include "UseTheForce/doForces_interface.h"
53 + #define __C
54 + #include "UseTheForce/DarkSide/fInteractionMap.h"
55   #include "utils/simError.h"
56   namespace oopse {
57  
58 < void ForceManager::calcForces(bool needPotential, bool needStress) {
58 >  void ForceManager::calcForces(bool needPotential, bool needStress) {
59  
60      if (!info_->isFortranInitialized()) {
61 <        info_->update();
61 >      info_->update();
62      }
63  
64      preCalculation();
# Line 67 | Line 69 | void ForceManager::calcForces(bool needPotential, bool
69  
70      postCalculation();
71          
72 < }
72 >  }
73  
74 < void ForceManager::preCalculation() {
74 >  void ForceManager::preCalculation() {
75      SimInfo::MoleculeIterator mi;
76      Molecule* mol;
77      Molecule::AtomIterator ai;
# Line 80 | Line 82 | void ForceManager::preCalculation() {
82      // forces are zeroed here, before any are accumulated.
83      // NOTE: do not rezero the forces in Fortran.
84      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
85 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
86 <            atom->zeroForcesAndTorques();
87 <        }
85 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
86 >        atom->zeroForcesAndTorques();
87 >      }
88          
89 <        //change the positions of atoms which belong to the rigidbodies
90 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
91 <            rb->zeroForcesAndTorques();
92 <        }        
89 >      //change the positions of atoms which belong to the rigidbodies
90 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
91 >        rb->zeroForcesAndTorques();
92 >      }        
93      }
94      
95 < }
95 >  }
96  
97 < void ForceManager::calcShortRangeInteraction() {
97 >  void ForceManager::calcShortRangeInteraction() {
98      Molecule* mol;
99      RigidBody* rb;
100      Bond* bond;
# Line 107 | Line 109 | void ForceManager::calcShortRangeInteraction() {
109      //calculate short range interactions    
110      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
111  
112 <        //change the positions of atoms which belong to the rigidbodies
113 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
114 <            rb->updateAtoms();
115 <        }
112 >      //change the positions of atoms which belong to the rigidbodies
113 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
114 >        rb->updateAtoms();
115 >      }
116  
117 <        for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
118 <            bond->calcForce();
119 <        }
117 >      for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
118 >        bond->calcForce();
119 >      }
120  
121 <        for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
122 <            bend->calcForce();
123 <        }
121 >      for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
122 >        bend->calcForce();
123 >      }
124  
125 <        for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
126 <            torsion->calcForce();
127 <        }
125 >      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
126 >        torsion->calcForce();
127 >      }
128  
129      }
130      
131      double  shortRangePotential = 0.0;
132      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
133 <        shortRangePotential += mol->getPotential();
133 >      shortRangePotential += mol->getPotential();
134      }
135  
136      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
137      curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
138 < }
138 >  }
139  
140 < void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
140 >  void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
141      Snapshot* curSnapshot;
142      DataStorage* config;
143      double* frc;
# Line 163 | Line 165 | void ForceManager::calcLongRangeInteraction(bool needP
165      CutoffGroup* cg;
166      Vector3d com;
167      std::vector<Vector3d> rcGroup;
166 –    
167 –    if(info_->getNCutoffGroups() > 0){
168  
169 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
169 >    if(info_->getNCutoffGroups() > 0){
170 >
171 >      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
172          for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
173 <            cg->getCOM(com);
174 <            rcGroup.push_back(com);
173 >          cg->getCOM(com);
174 >          rcGroup.push_back(com);
175          }
176 <    }// end for (mol)
176 >      }// end for (mol)
177        
178 <        rc = rcGroup[0].getArrayPointer();
178 >      rc = rcGroup[0].getArrayPointer();
179      } else {
180 <        // center of mass of the group is the same as position of the atom  if cutoff group does not exist
181 <        rc = pos;
180 >      // center of mass of the group is the same as position of the atom  if cutoff group does not exist
181 >      rc = pos;
182      }
183    
184      //initialize data before passing to fortran
185 <    double longRangePotential = 0.0;
185 >    double longRangePotential[LR_POT_TYPES];
186 >    double lrPot = 0.0;
187 >    
188      Mat3x3d tau;
189      short int passedCalcPot = needPotential;
190      short int passedCalcStress = needStress;
191      int isError = 0;
192  
193 +    for (int i=0; i<LR_POT_TYPES;i++){
194 +      longRangePotential[i]=0.0; //Initialize array
195 +    }
196 +
197 +
198 +
199      doForceLoop( pos,
200 <            rc,
201 <            A,
202 <            electroFrame,
203 <            frc,
204 <            trq,
205 <            tau.getArrayPointer(),
206 <            &longRangePotential,
207 <            &passedCalcPot,
208 <            &passedCalcStress,
209 <            &isError );
200 >                 rc,
201 >                 A,
202 >                 electroFrame,
203 >                 frc,
204 >                 trq,
205 >                 tau.getArrayPointer(),
206 >                 longRangePotential,
207 >                 &passedCalcPot,
208 >                 &passedCalcStress,
209 >                 &isError );
210  
211      if( isError ){
212 <        sprintf( painCave.errMsg,
213 <             "Error returned from the fortran force calculation.\n" );
214 <        painCave.isFatal = 1;
215 <        simError();
212 >      sprintf( painCave.errMsg,
213 >               "Error returned from the fortran force calculation.\n" );
214 >      painCave.isFatal = 1;
215 >      simError();
216      }
217 +    for (int i=0; i<LR_POT_TYPES;i++){
218 +      lrPot += longRangePotential[i]; //Quick hack
219 +    }
220  
221      //store the tau and long range potential    
222 <    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential;
222 >    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
223 >    //  curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential;
224      curSnapshot->statData.setTau(tau);
225 < }
225 >  }
226  
227  
228 < void ForceManager::postCalculation() {
228 >  void ForceManager::postCalculation() {
229      SimInfo::MoleculeIterator mi;
230      Molecule* mol;
231      Molecule::RigidBodyIterator rbIter;
# Line 219 | Line 233 | void ForceManager::postCalculation() {
233      
234      // collect the atomic forces onto rigid bodies
235      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
236 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
237 <            rb->calcForcesAndTorques();
238 <        }
236 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
237 >        rb->calcForcesAndTorques();
238 >      }
239      }
240  
241 < }
241 >  }
242  
243   } //end namespace oopse

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines