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 770 by tim, Fri Dec 2 15:38:03 2005 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# Line 6 | Line 6
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
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
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.
# Line 37 | Line 28
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   /**
# Line 50 | Line 50
50   #include "brains/ForceManager.hpp"
51   #include "primitives/Molecule.hpp"
52   #include "UseTheForce/doForces_interface.h"
53 < #define __C
53 > #define __OPENMD_C
54   #include "UseTheForce/DarkSide/fInteractionMap.h"
55   #include "utils/simError.h"
56 + #include "primitives/Bond.hpp"
57   #include "primitives/Bend.hpp"
58 < #include "primitives/Bend.hpp"
59 < namespace oopse {
58 > #include "primitives/Torsion.hpp"
59 > #include "primitives/Inversion.hpp"
60 > namespace OpenMD {
61  
60 – /*
61 –  struct BendOrderStruct {
62 –    Bend* bend;
63 –    BendDataSet dataSet;
64 –  };
65 –  struct TorsionOrderStruct {
66 –    Torsion* torsion;
67 –    TorsionDataSet dataSet;
68 –  };
69 –
70 –  bool  BendSortFunctor(const BendOrderStruct& b1, const BendOrderStruct& b2) {
71 –    return b1.dataSet.deltaV < b2.dataSet.deltaV;
72 –  }
73 –
74 –  bool  TorsionSortFunctor(const TorsionOrderStruct& t1, const TorsionOrderStruct& t2) {
75 –    return t1.dataSet.deltaV < t2.dataSet.deltaV;
76 –  }
77 –  */
62    void ForceManager::calcForces(bool needPotential, bool needStress) {
63 <
63 >    
64      if (!info_->isFortranInitialized()) {
65        info_->update();
66      }
67 <
67 >    
68      preCalculation();
69      
70      calcShortRangeInteraction();
71  
72      calcLongRangeInteraction(needPotential, needStress);
73  
74 <    postCalculation();
91 <
92 < /*
93 <    std::vector<BendOrderStruct> bendOrderStruct;
94 <    for(std::map<Bend*, BendDataSet>::iterator i = bendDataSets.begin(); i != bendDataSets.end(); ++i) {
95 <        BendOrderStruct tmp;
96 <        tmp.bend= const_cast<Bend*>(i->first);
97 <        tmp.dataSet = i->second;
98 <        bendOrderStruct.push_back(tmp);
99 <    }
100 <
101 <    std::vector<TorsionOrderStruct> torsionOrderStruct;
102 <    for(std::map<Torsion*, TorsionDataSet>::iterator j = torsionDataSets.begin(); j != torsionDataSets.end(); ++j) {
103 <        TorsionOrderStruct tmp;
104 <        tmp.torsion = const_cast<Torsion*>(j->first);
105 <        tmp.dataSet = j->second;
106 <        torsionOrderStruct.push_back(tmp);
107 <    }
74 >    postCalculation(needStress);
75      
109 –    std::sort(bendOrderStruct.begin(), bendOrderStruct.end(), std::ptr_fun(BendSortFunctor));
110 –    std::sort(torsionOrderStruct.begin(), torsionOrderStruct.end(), std::ptr_fun(TorsionSortFunctor));
111 –    std::cout << "bend" << std::endl;
112 –    for (std::vector<BendOrderStruct>::iterator k = bendOrderStruct.begin(); k != bendOrderStruct.end(); ++k) {
113 –        Bend* bend = k->bend;
114 –        std::cout << "atom1=" <<bend->getAtomA()->getGlobalIndex() << ",atom2 = "<< bend->getAtomB()->getGlobalIndex() << ",atom3="<<bend->getAtomC()->getGlobalIndex() << " ";
115 –        std::cout << "deltaV=" << k->dataSet.deltaV << ",p_theta=" << k->dataSet.prev.angle <<",p_pot=" << k->dataSet.prev.potential<< ",c_theta=" << k->dataSet.curr.angle << ", c_pot = " << k->dataSet.curr.potential <<std::endl;
116 –    }
117 –    std::cout << "torsio" << std::endl;
118 –    for (std::vector<TorsionOrderStruct>::iterator l = torsionOrderStruct.begin(); l != torsionOrderStruct.end(); ++l) {
119 –        Torsion* torsion = l->torsion;
120 –        std::cout << "atom1=" <<torsion->getAtomA()->getGlobalIndex() << ",atom2 = "<< torsion->getAtomB()->getGlobalIndex() << ",atom3="<<torsion->getAtomC()->getGlobalIndex() << ",atom4="<<torsion->getAtomD()->getGlobalIndex()<< " ";
121 –        std::cout << "deltaV=" << l->dataSet.deltaV << ",p_theta=" << l->dataSet.prev.angle <<",p_pot=" << l->dataSet.prev.potential<< ",c_theta=" << l->dataSet.curr.angle << ", c_pot = " << l->dataSet.curr.potential <<std::endl;
122 –    }
123 –   */
76    }
77 <
77 >  
78    void ForceManager::preCalculation() {
79      SimInfo::MoleculeIterator mi;
80      Molecule* mol;
# Line 133 | Line 85 | namespace oopse {
85      
86      // forces are zeroed here, before any are accumulated.
87      // NOTE: do not rezero the forces in Fortran.
88 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
88 >    
89 >    for (mol = info_->beginMolecule(mi); mol != NULL;
90 >         mol = info_->nextMolecule(mi)) {
91        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
92          atom->zeroForcesAndTorques();
93        }
94 <        
94 >          
95        //change the positions of atoms which belong to the rigidbodies
96 <      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
96 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
97 >           rb = mol->nextRigidBody(rbIter)) {
98          rb->zeroForcesAndTorques();
99        }        
100 +          
101      }
102      
103 +    // Zero out the stress tensor
104 +    tau *= 0.0;
105 +    
106    }
107 <
107 >  
108    void ForceManager::calcShortRangeInteraction() {
109      Molecule* mol;
110      RigidBody* rb;
111      Bond* bond;
112      Bend* bend;
113      Torsion* torsion;
114 +    Inversion* inversion;
115      SimInfo::MoleculeIterator mi;
116      Molecule::RigidBodyIterator rbIter;
117      Molecule::BondIterator bondIter;;
118      Molecule::BendIterator  bendIter;
119      Molecule::TorsionIterator  torsionIter;
120 <    double bondPotential = 0.0;
121 <    double bendPotential = 0.0;
122 <    double torsionPotential = 0.0;
120 >    Molecule::InversionIterator  inversionIter;
121 >    RealType bondPotential = 0.0;
122 >    RealType bendPotential = 0.0;
123 >    RealType torsionPotential = 0.0;
124 >    RealType inversionPotential = 0.0;
125  
126      //calculate short range interactions    
127 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
127 >    for (mol = info_->beginMolecule(mi); mol != NULL;
128 >         mol = info_->nextMolecule(mi)) {
129  
130        //change the positions of atoms which belong to the rigidbodies
131 <      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
132 <          rb->updateAtoms();
131 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
132 >           rb = mol->nextRigidBody(rbIter)) {
133 >        rb->updateAtoms();
134        }
135  
136 <      for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
136 >      for (bond = mol->beginBond(bondIter); bond != NULL;
137 >           bond = mol->nextBond(bondIter)) {
138          bond->calcForce();
139          bondPotential += bond->getPotential();
140        }
141  
142 <
143 <      for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
144 <
145 <          double angle;
146 <            bend->calcForce(angle);
147 <          double currBendPot = bend->getPotential();          
148 <            bendPotential += bend->getPotential();
149 <          std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
150 <          if (i == bendDataSets.end()) {
151 <            BendDataSet dataSet;
152 <            dataSet.prev.angle = dataSet.curr.angle = angle;
153 <            dataSet.prev.potential = dataSet.curr.potential = currBendPot;
154 <            dataSet.deltaV = 0.0;
155 <            bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
156 <          }else {
157 <            i->second.prev.angle = i->second.curr.angle;
158 <            i->second.prev.potential = i->second.curr.potential;
159 <            i->second.curr.angle = angle;
160 <            i->second.curr.potential = currBendPot;
161 <            i->second.deltaV =  fabs(i->second.curr.potential -  i->second.prev.potential);
162 <          }
163 <      }
199 <
200 <      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
201 <        double angle;
202 <          torsion->calcForce(angle);
203 <        double currTorsionPot = torsion->getPotential();
204 <          torsionPotential += torsion->getPotential();
205 <          std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
206 <          if (i == torsionDataSets.end()) {
207 <            TorsionDataSet dataSet;
208 <            dataSet.prev.angle = dataSet.curr.angle = angle;
209 <            dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
210 <            dataSet.deltaV = 0.0;
211 <            torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
212 <          }else {
213 <            i->second.prev.angle = i->second.curr.angle;
214 <            i->second.prev.potential = i->second.curr.potential;
215 <            i->second.curr.angle = angle;
216 <            i->second.curr.potential = currTorsionPot;
217 <            i->second.deltaV =  fabs(i->second.curr.potential -  i->second.prev.potential);
218 <          }      
142 >      for (bend = mol->beginBend(bendIter); bend != NULL;
143 >           bend = mol->nextBend(bendIter)) {
144 >        
145 >        RealType angle;
146 >        bend->calcForce(angle);
147 >        RealType currBendPot = bend->getPotential();          
148 >        bendPotential += bend->getPotential();
149 >        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
150 >        if (i == bendDataSets.end()) {
151 >          BendDataSet dataSet;
152 >          dataSet.prev.angle = dataSet.curr.angle = angle;
153 >          dataSet.prev.potential = dataSet.curr.potential = currBendPot;
154 >          dataSet.deltaV = 0.0;
155 >          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
156 >        }else {
157 >          i->second.prev.angle = i->second.curr.angle;
158 >          i->second.prev.potential = i->second.curr.potential;
159 >          i->second.curr.angle = angle;
160 >          i->second.curr.potential = currBendPot;
161 >          i->second.deltaV =  fabs(i->second.curr.potential -  
162 >                                   i->second.prev.potential);
163 >        }
164        }
165 +      
166 +      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
167 +           torsion = mol->nextTorsion(torsionIter)) {
168 +        RealType angle;
169 +        torsion->calcForce(angle);
170 +        RealType currTorsionPot = torsion->getPotential();
171 +        torsionPotential += torsion->getPotential();
172 +        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
173 +        if (i == torsionDataSets.end()) {
174 +          TorsionDataSet dataSet;
175 +          dataSet.prev.angle = dataSet.curr.angle = angle;
176 +          dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
177 +          dataSet.deltaV = 0.0;
178 +          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
179 +        }else {
180 +          i->second.prev.angle = i->second.curr.angle;
181 +          i->second.prev.potential = i->second.curr.potential;
182 +          i->second.curr.angle = angle;
183 +          i->second.curr.potential = currTorsionPot;
184 +          i->second.deltaV =  fabs(i->second.curr.potential -  
185 +                                   i->second.prev.potential);
186 +        }      
187 +      }      
188  
189 +      for (inversion = mol->beginInversion(inversionIter);
190 +           inversion != NULL;
191 +           inversion = mol->nextInversion(inversionIter)) {
192 +        RealType angle;
193 +        inversion->calcForce(angle);
194 +        RealType currInversionPot = inversion->getPotential();
195 +        inversionPotential += inversion->getPotential();
196 +        std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
197 +        if (i == inversionDataSets.end()) {
198 +          InversionDataSet dataSet;
199 +          dataSet.prev.angle = dataSet.curr.angle = angle;
200 +          dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
201 +          dataSet.deltaV = 0.0;
202 +          inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
203 +        }else {
204 +          i->second.prev.angle = i->second.curr.angle;
205 +          i->second.prev.potential = i->second.curr.potential;
206 +          i->second.curr.angle = angle;
207 +          i->second.curr.potential = currInversionPot;
208 +          i->second.deltaV =  fabs(i->second.curr.potential -  
209 +                                   i->second.prev.potential);
210 +        }      
211 +      }      
212      }
213      
214 <    double  shortRangePotential = bondPotential + bendPotential + torsionPotential;    
214 >    RealType  shortRangePotential = bondPotential + bendPotential +
215 >      torsionPotential +  inversionPotential;    
216      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
217      curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
218      curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
219      curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
220      curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
221 +    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
222      
223    }
224 <
225 <  void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
224 >  
225 >  void ForceManager::calcLongRangeInteraction(bool needPotential,
226 >                                              bool needStress) {
227      Snapshot* curSnapshot;
228      DataStorage* config;
229 <    double* frc;
230 <    double* pos;
231 <    double* trq;
232 <    double* A;
233 <    double* electroFrame;
234 <    double* rc;
229 >    RealType* frc;
230 >    RealType* pos;
231 >    RealType* trq;
232 >    RealType* A;
233 >    RealType* electroFrame;
234 >    RealType* rc;
235 >    RealType* particlePot;
236      
237      //get current snapshot from SimInfo
238      curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
239 <
239 >    
240      //get array pointers
241      config = &(curSnapshot->atomData);
242      frc = config->getArrayPointer(DataStorage::dslForce);
# Line 249 | Line 244 | namespace oopse {
244      trq = config->getArrayPointer(DataStorage::dslTorque);
245      A   = config->getArrayPointer(DataStorage::dslAmat);
246      electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
247 +    particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
248  
249      //calculate the center of mass of cutoff group
250      SimInfo::MoleculeIterator mi;
# Line 257 | Line 253 | namespace oopse {
253      CutoffGroup* cg;
254      Vector3d com;
255      std::vector<Vector3d> rcGroup;
256 <
256 >    
257      if(info_->getNCutoffGroups() > 0){
258 <
259 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
260 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
258 >      
259 >      for (mol = info_->beginMolecule(mi); mol != NULL;
260 >           mol = info_->nextMolecule(mi)) {
261 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
262 >            cg = mol->nextCutoffGroup(ci)) {
263            cg->getCOM(com);
264            rcGroup.push_back(com);
265          }
# Line 269 | Line 267 | namespace oopse {
267        
268        rc = rcGroup[0].getArrayPointer();
269      } else {
270 <      // center of mass of the group is the same as position of the atom  if cutoff group does not exist
270 >      // center of mass of the group is the same as position of the atom  
271 >      // if cutoff group does not exist
272        rc = pos;
273      }
275 –  
276 –    //initialize data before passing to fortran
277 –    double longRangePotential[LR_POT_TYPES];
278 –    double lrPot = 0.0;
274      
275 <    Mat3x3d tau;
275 >    //initialize data before passing to fortran
276 >    RealType longRangePotential[LR_POT_TYPES];
277 >    RealType lrPot = 0.0;
278 >    Vector3d totalDipole;
279      short int passedCalcPot = needPotential;
280      short int passedCalcStress = needStress;
281      int isError = 0;
# Line 285 | Line 283 | namespace oopse {
283      for (int i=0; i<LR_POT_TYPES;i++){
284        longRangePotential[i]=0.0; //Initialize array
285      }
286 <
287 <    doForceLoop( pos,
288 <                 rc,
289 <                 A,
290 <                 electroFrame,
291 <                 frc,
292 <                 trq,
293 <                 tau.getArrayPointer(),
294 <                 longRangePotential,
295 <                 &passedCalcPot,
296 <                 &passedCalcStress,
297 <                 &isError );
298 <
286 >    
287 >    doForceLoop(pos,
288 >                rc,
289 >                A,
290 >                electroFrame,
291 >                frc,
292 >                trq,
293 >                tau.getArrayPointer(),
294 >                longRangePotential,
295 >                particlePot,
296 >                &passedCalcPot,
297 >                &passedCalcStress,
298 >                &isError );
299 >    
300      if( isError ){
301        sprintf( painCave.errMsg,
302                 "Error returned from the fortran force calculation.\n" );
# Line 307 | Line 306 | namespace oopse {
306      for (int i=0; i<LR_POT_TYPES;i++){
307        lrPot += longRangePotential[i]; //Quick hack
308      }
309 <
309 >    
310 >    // grab the simulation box dipole moment if specified
311 >    if (info_->getCalcBoxDipole()){
312 >      getAccumulatedBoxDipole(totalDipole.getArrayPointer());
313 >      
314 >      curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0);
315 >      curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1);
316 >      curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2);
317 >    }
318 >    
319      //store the tau and long range potential    
320      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
321      curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
322      curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
315 –
316 –    curSnapshot->statData.setTau(tau);
323    }
324  
325 <
326 <  void ForceManager::postCalculation() {
325 >  
326 >  void ForceManager::postCalculation(bool needStress) {
327      SimInfo::MoleculeIterator mi;
328      Molecule* mol;
329      Molecule::RigidBodyIterator rbIter;
330      RigidBody* rb;
331 +    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
332      
333      // collect the atomic forces onto rigid bodies
334 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
335 <      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
336 <        rb->calcForcesAndTorques();
334 >    
335 >    for (mol = info_->beginMolecule(mi); mol != NULL;
336 >         mol = info_->nextMolecule(mi)) {
337 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
338 >           rb = mol->nextRigidBody(rbIter)) {
339 >        if (needStress) {          
340 >          Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
341 >          tau += rbTau;
342 >        } else{
343 >          rb->calcForcesAndTorques();
344 >        }
345        }
346      }
347  
348 +    if (needStress) {
349 + #ifdef IS_MPI
350 +      Mat3x3d tmpTau(tau);
351 +      MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
352 +                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
353 + #endif
354 +      curSnapshot->statData.setTau(tau);
355 +    }
356    }
357  
358 < } //end namespace oopse
358 > } //end namespace OpenMD

Diff Legend

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