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

Comparing:
trunk/src/brains/ForceManager.cpp (property svn:keywords), Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
branches/development/src/brains/ForceManager.cpp (property svn:keywords), Revision 1503 by gezelter, Sat Oct 2 19:54:41 2010 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines