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 1215 by xsun, Wed Jan 23 21:22:37 2008 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 + #include "primitives/Bond.hpp"
57 + #include "primitives/Bend.hpp"
58   namespace oopse {
59  
60 < void ForceManager::calcForces(bool needPotential, bool needStress) {
61 <
60 >  void ForceManager::calcForces(bool needPotential, bool needStress) {
61 >    
62      if (!info_->isFortranInitialized()) {
63 <        info_->update();
63 >      info_->update();
64      }
65 <
65 >    
66      preCalculation();
67      
68      calcShortRangeInteraction();
69  
70      calcLongRangeInteraction(needPotential, needStress);
71  
72 <    postCalculation();
73 <        
74 < }
75 <
76 < void ForceManager::preCalculation() {
72 >    postCalculation(needStress);
73 >    
74 >  }
75 >  
76 >  void ForceManager::preCalculation() {
77      SimInfo::MoleculeIterator mi;
78      Molecule* mol;
79      Molecule::AtomIterator ai;
# Line 79 | Line 83 | void ForceManager::preCalculation() {
83      
84      // forces are zeroed here, before any are accumulated.
85      // NOTE: do not rezero the forces in Fortran.
86 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
87 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
88 <            atom->zeroForcesAndTorques();
89 <        }
90 <        
91 <        //change the positions of atoms which belong to the rigidbodies
92 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
93 <            rb->zeroForcesAndTorques();
94 <        }        
86 >
87 >    for (mol = info_->beginMolecule(mi); mol != NULL;
88 >         mol = info_->nextMolecule(mi)) {
89 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
90 >        atom->zeroForcesAndTorques();
91 >      }
92 >      
93 >      //change the positions of atoms which belong to the rigidbodies
94 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
95 >           rb = mol->nextRigidBody(rbIter)) {
96 >        rb->zeroForcesAndTorques();
97 >      }        
98      }
99      
100 < }
101 <
102 < void ForceManager::calcShortRangeInteraction() {
100 >    // Zero out the stress tensor
101 >    tau *= 0.0;
102 >    
103 >  }
104 >  
105 >  void ForceManager::calcShortRangeInteraction() {
106      Molecule* mol;
107      RigidBody* rb;
108      Bond* bond;
# Line 103 | Line 113 | void ForceManager::calcShortRangeInteraction() {
113      Molecule::BondIterator bondIter;;
114      Molecule::BendIterator  bendIter;
115      Molecule::TorsionIterator  torsionIter;
116 +    RealType bondPotential = 0.0;
117 +    RealType bendPotential = 0.0;
118 +    RealType torsionPotential = 0.0;
119  
120      //calculate short range interactions    
121 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
121 >    for (mol = info_->beginMolecule(mi); mol != NULL;
122 >         mol = info_->nextMolecule(mi)) {
123  
124 <        //change the positions of atoms which belong to the rigidbodies
125 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
126 <            rb->updateAtoms();
127 <        }
124 >      //change the positions of atoms which belong to the rigidbodies
125 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
126 >           rb = mol->nextRigidBody(rbIter)) {
127 >        rb->updateAtoms();
128 >      }
129  
130 <        for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
131 <            bond->calcForce();
132 <        }
130 >      for (bond = mol->beginBond(bondIter); bond != NULL;
131 >           bond = mol->nextBond(bondIter)) {
132 >        bond->calcForce();
133 >        bondPotential += bond->getPotential();
134 >      }
135  
136 <        for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
137 <            bend->calcForce();
136 >      for (bend = mol->beginBend(bendIter); bend != NULL;
137 >           bend = mol->nextBend(bendIter)) {
138 >        
139 >        RealType angle;
140 >        bend->calcForce(angle);
141 >        RealType currBendPot = bend->getPotential();          
142 >        bendPotential += bend->getPotential();
143 >        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
144 >        if (i == bendDataSets.end()) {
145 >          BendDataSet dataSet;
146 >          dataSet.prev.angle = dataSet.curr.angle = angle;
147 >          dataSet.prev.potential = dataSet.curr.potential = currBendPot;
148 >          dataSet.deltaV = 0.0;
149 >          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
150 >        }else {
151 >          i->second.prev.angle = i->second.curr.angle;
152 >          i->second.prev.potential = i->second.curr.potential;
153 >          i->second.curr.angle = angle;
154 >          i->second.curr.potential = currBendPot;
155 >          i->second.deltaV =  fabs(i->second.curr.potential -  
156 >                                   i->second.prev.potential);
157          }
158 <
159 <        for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
160 <            torsion->calcForce();
161 <        }
162 <
158 >      }
159 >      
160 >      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
161 >           torsion = mol->nextTorsion(torsionIter)) {
162 >        RealType angle;
163 >        torsion->calcForce(angle);
164 >        RealType currTorsionPot = torsion->getPotential();
165 >        torsionPotential += torsion->getPotential();
166 >        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
167 >        if (i == torsionDataSets.end()) {
168 >          TorsionDataSet dataSet;
169 >          dataSet.prev.angle = dataSet.curr.angle = angle;
170 >          dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
171 >          dataSet.deltaV = 0.0;
172 >          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
173 >        }else {
174 >          i->second.prev.angle = i->second.curr.angle;
175 >          i->second.prev.potential = i->second.curr.potential;
176 >          i->second.curr.angle = angle;
177 >          i->second.curr.potential = currTorsionPot;
178 >          i->second.deltaV =  fabs(i->second.curr.potential -  
179 >                                   i->second.prev.potential);
180 >        }      
181 >      }      
182      }
183      
184 <    double  shortRangePotential = 0.0;
185 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
131 <        shortRangePotential += mol->getPotential();
132 <    }
133 <
184 >    RealType  shortRangePotential = bondPotential + bendPotential +
185 >      torsionPotential;    
186      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
187      curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
188 < }
189 <
190 < void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
188 >    curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
189 >    curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
190 >    curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
191 >    
192 >  }
193 >  
194 >  void ForceManager::calcLongRangeInteraction(bool needPotential,
195 >                                              bool needStress) {
196      Snapshot* curSnapshot;
197      DataStorage* config;
198 <    double* frc;
199 <    double* pos;
200 <    double* trq;
201 <    double* A;
202 <    double* electroFrame;
203 <    double* rc;
198 >    RealType* frc;
199 >    RealType* pos;
200 >    RealType* trq;
201 >    RealType* A;
202 >    RealType* electroFrame;
203 >    RealType* rc;
204      
205      //get current snapshot from SimInfo
206      curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
207 <
207 >    
208      //get array pointers
209      config = &(curSnapshot->atomData);
210      frc = config->getArrayPointer(DataStorage::dslForce);
# Line 165 | Line 222 | void ForceManager::calcLongRangeInteraction(bool needP
222      std::vector<Vector3d> rcGroup;
223      
224      if(info_->getNCutoffGroups() > 0){
225 <
226 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
227 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
228 <            cg->getCOM(com);
229 <            rcGroup.push_back(com);
225 >      
226 >      for (mol = info_->beginMolecule(mi); mol != NULL;
227 >           mol = info_->nextMolecule(mi)) {
228 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
229 >            cg = mol->nextCutoffGroup(ci)) {
230 >          cg->getCOM(com);
231 >          rcGroup.push_back(com);
232          }
233 <    }// end for (mol)
233 >      }// end for (mol)
234        
235 <        rc = rcGroup[0].getArrayPointer();
235 >      rc = rcGroup[0].getArrayPointer();
236      } else {
237 <        // center of mass of the group is the same as position of the atom  if cutoff group does not exist
238 <        rc = pos;
237 >      // center of mass of the group is the same as position of the atom  
238 >      // if cutoff group does not exist
239 >      rc = pos;
240      }
241 <  
241 >    
242      //initialize data before passing to fortran
243 <    double longRangePotential = 0.0;
244 <    Mat3x3d tau;
243 >    RealType longRangePotential[LR_POT_TYPES];
244 >    RealType lrPot = 0.0;
245 >    Vector3d totalDipole;
246      short int passedCalcPot = needPotential;
247      short int passedCalcStress = needStress;
248      int isError = 0;
249  
250 <    doForceLoop( pos,
251 <            rc,
252 <            A,
253 <            electroFrame,
254 <            frc,
255 <            trq,
256 <            tau.getArrayPointer(),
257 <            &longRangePotential,
258 <            &passedCalcPot,
259 <            &passedCalcStress,
260 <            &isError );
261 <
250 >    for (int i=0; i<LR_POT_TYPES;i++){
251 >      longRangePotential[i]=0.0; //Initialize array
252 >    }
253 >    
254 >    doForceLoop(pos,
255 >                rc,
256 >                A,
257 >                electroFrame,
258 >                frc,
259 >                trq,
260 >                tau.getArrayPointer(),
261 >                longRangePotential,
262 >                &passedCalcPot,
263 >                &passedCalcStress,
264 >                &isError );
265 >    
266      if( isError ){
267 <        sprintf( painCave.errMsg,
268 <             "Error returned from the fortran force calculation.\n" );
269 <        painCave.isFatal = 1;
270 <        simError();
267 >      sprintf( painCave.errMsg,
268 >               "Error returned from the fortran force calculation.\n" );
269 >      painCave.isFatal = 1;
270 >      simError();
271      }
272 <
272 >    for (int i=0; i<LR_POT_TYPES;i++){
273 >      lrPot += longRangePotential[i]; //Quick hack
274 >    }
275 >    
276 >    // grab the simulation box dipole moment if specified
277 >    if (info_->getCalcBoxDipole()){
278 >      getAccumulatedBoxDipole(totalDipole.getArrayPointer());
279 >      
280 >      curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0);
281 >      curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1);
282 >      curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2);
283 >    }
284 >    
285      //store the tau and long range potential    
286 <    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential;
287 <    curSnapshot->statData.setTau(tau);
288 < }
286 >    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
287 >    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
288 >    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
289 >  }
290  
291 <
292 < void ForceManager::postCalculation() {
291 >  
292 >  void ForceManager::postCalculation(bool needStress) {
293      SimInfo::MoleculeIterator mi;
294      Molecule* mol;
295      Molecule::RigidBodyIterator rbIter;
296      RigidBody* rb;
297 +    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
298      
299      // collect the atomic forces onto rigid bodies
300 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
301 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
302 <            rb->calcForcesAndTorques();
300 >    
301 >    for (mol = info_->beginMolecule(mi); mol != NULL;
302 >         mol = info_->nextMolecule(mi)) {
303 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
304 >           rb = mol->nextRigidBody(rbIter)) {
305 >        if (needStress) {          
306 >          Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
307 >          tau += rbTau;
308 >        } else{
309 >          rb->calcForcesAndTorques();
310          }
311 +      }
312      }
313  
314 < }
314 >    if (needStress) {
315 > #ifdef IS_MPI
316 >      Mat3x3d tmpTau(tau);
317 >      MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
318 >                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
319 > #endif
320 >      curSnapshot->statData.setTau(tau);
321 >    }
322 >  }
323  
324   } //end namespace oopse

Diff Legend

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