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 681 by tim, Mon Oct 17 23:13:44 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;
131 >
132 >    double bondPotential = 0.0;
133 >    double bendPotential = 0.0;
134 >    double torsionPotential = 0.0;
135 >
136      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
131 –        shortRangePotential += mol->getPotential();
132 –    }
137  
138 +      for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
139 +          bondPotential += bond->getPotential();
140 +      }
141 +
142 +      for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
143 +          bendPotential += bend->getPotential();
144 +      }
145 +
146 +      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
147 +          torsionPotential += torsion->getPotential();
148 +      }
149 +
150 +    }    
151 +
152 +    double  shortRangePotential = bondPotential + bendPotential + torsionPotential;    
153      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
154      curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
155 < }
155 >    curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
156 >    curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
157 >    curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
158 >    
159 >  }
160  
161 < void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
161 >  void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
162      Snapshot* curSnapshot;
163      DataStorage* config;
164      double* frc;
# Line 163 | Line 186 | void ForceManager::calcLongRangeInteraction(bool needP
186      CutoffGroup* cg;
187      Vector3d com;
188      std::vector<Vector3d> rcGroup;
166 –    
167 –    if(info_->getNCutoffGroups() > 0){
189  
190 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
190 >    if(info_->getNCutoffGroups() > 0){
191 >
192 >      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
193          for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
194 <            cg->getCOM(com);
195 <            rcGroup.push_back(com);
194 >          cg->getCOM(com);
195 >          rcGroup.push_back(com);
196          }
197 <    }// end for (mol)
197 >      }// end for (mol)
198        
199 <        rc = rcGroup[0].getArrayPointer();
199 >      rc = rcGroup[0].getArrayPointer();
200      } else {
201 <        // center of mass of the group is the same as position of the atom  if cutoff group does not exist
202 <        rc = pos;
201 >      // center of mass of the group is the same as position of the atom  if cutoff group does not exist
202 >      rc = pos;
203      }
204    
205      //initialize data before passing to fortran
206 <    double longRangePotential = 0.0;
206 >    double longRangePotential[LR_POT_TYPES];
207 >    double lrPot = 0.0;
208 >    
209      Mat3x3d tau;
210      short int passedCalcPot = needPotential;
211      short int passedCalcStress = needStress;
212      int isError = 0;
213  
214 +    for (int i=0; i<LR_POT_TYPES;i++){
215 +      longRangePotential[i]=0.0; //Initialize array
216 +    }
217 +
218 +
219 +
220      doForceLoop( pos,
221 <            rc,
222 <            A,
223 <            electroFrame,
224 <            frc,
225 <            trq,
226 <            tau.getArrayPointer(),
227 <            &longRangePotential,
228 <            &passedCalcPot,
229 <            &passedCalcStress,
230 <            &isError );
221 >                 rc,
222 >                 A,
223 >                 electroFrame,
224 >                 frc,
225 >                 trq,
226 >                 tau.getArrayPointer(),
227 >                 longRangePotential,
228 >                 &passedCalcPot,
229 >                 &passedCalcStress,
230 >                 &isError );
231  
232      if( isError ){
233 <        sprintf( painCave.errMsg,
234 <             "Error returned from the fortran force calculation.\n" );
235 <        painCave.isFatal = 1;
236 <        simError();
233 >      sprintf( painCave.errMsg,
234 >               "Error returned from the fortran force calculation.\n" );
235 >      painCave.isFatal = 1;
236 >      simError();
237      }
238 +    for (int i=0; i<LR_POT_TYPES;i++){
239 +      lrPot += longRangePotential[i]; //Quick hack
240 +    }
241  
242      //store the tau and long range potential    
243 <    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential;
243 >    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
244 >    //  curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential;
245 >    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
246 >    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
247 >
248      curSnapshot->statData.setTau(tau);
249 < }
249 >  }
250  
251  
252 < void ForceManager::postCalculation() {
252 >  void ForceManager::postCalculation() {
253      SimInfo::MoleculeIterator mi;
254      Molecule* mol;
255      Molecule::RigidBodyIterator rbIter;
# Line 219 | Line 257 | void ForceManager::postCalculation() {
257      
258      // collect the atomic forces onto rigid bodies
259      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
260 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
261 <            rb->calcForcesAndTorques();
262 <        }
260 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
261 >        rb->calcForcesAndTorques();
262 >      }
263      }
264  
265 < }
265 >  }
266  
267   } //end namespace oopse

Diff Legend

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