ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/brains/ForceManager.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-4/src/brains/ForceManager.cpp (file contents):
Revision 1722 by tim, Tue Nov 9 23:11:39 2004 UTC vs.
Revision 1804 by tim, Tue Nov 30 19:58:25 2004 UTC

# Line 32 | Line 32
32    */
33  
34   #include "brains/ForceManager.hpp"
35 <
35 > #include "primitives/Molecule.hpp"
36 > #include "UseTheForce/doForces_interface.h"
37 > #include "utils/simError.h"
38   namespace oopse {
39  
40   void ForceManager::calcForces(bool needPotential, bool needStress) {
41  
42 <    Snapshot* curSnapshot;
43 <    DataStorage* config;
44 <    double* frc;
43 <    double* pos;
44 <    double* trq;
45 <    double* A;
46 <    double* u_l;
47 <    double* rc;
42 >    if (!info_->isFortranInitialized()) {
43 >        info_->update();
44 >    }
45  
46 <    std::vector<Molecule*>::iterator mi;
46 >    preCalculation();
47 >    
48 >    calcShortRangeInteraction();
49 >
50 >    calcLongRangeInteraction(needPotential, needStress);
51 >
52 >    postCalculation();
53 >        
54 > }
55 >
56 > void ForceManager::preCalculation() {
57 >    SimInfo::MoleculeIterator mi;
58      Molecule* mol;
59 <    std::vector<Atom*>::iterator ai;
59 >    Molecule::AtomIterator ai;
60      Atom* atom;
53 –    std::vector<Atom*>::iterator ci;
54 –    CutoffGroup* cg;
55 –    Vector3d com;
56 –    std::vector<Vector3d> rcGroup;
61  
58 –    double longRangePotential;
59 –    double tau[9];
60 –    short int passedCalcPot;
61 –    short int passedCalcStress;
62 –    int isError;    
63 –
62      // forces are zeroed here, before any are accumulated.
63      // NOTE: do not rezero the forces in Fortran.
64      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
65          for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
66              atom->zeroForces();
67          }
68 <    }// end for (mol)
68 >    }
69 >    
70 > }
71  
72 < #ifdef PROFILE
73 <    startProfile(pro7);
74 < #endif
72 > void ForceManager::calcShortRangeInteraction() {
73 >    Molecule* mol;
74 >    RigidBody* rb;
75 >    Bond* bond;
76 >    Bend* bend;
77 >    Torsion* torsion;
78 >    SimInfo::MoleculeIterator mi;
79 >    Molecule::RigidBodyIterator rbIter;
80 >    Molecule::BondIterator bondIter;;
81 >    Molecule::BendIterator  bendIter;
82 >    Molecule::TorsionIterator  torsionIter;
83 >
84      //calculate short range interactions    
85      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
77 –        mol->calcForces();
78 –    }// end for (mol)
86  
87 < #ifdef PROFILE
88 <    endProfile( pro7 );
89 < #endif
87 >        //change the positions of atoms which belong to the rigidbodies
88 >        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
89 >            rb->updateAtoms();
90 >        }
91  
92 +        for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
93 +            bond->calcForce();
94 +        }
95 +
96 +        for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
97 +            bend->calcForce();
98 +        }
99 +
100 +        for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
101 +            torsion->calcForce();
102 +        }
103 +
104 +    }
105 +    
106 + }
107 +
108 + void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
109 +    Snapshot* curSnapshot;
110 +    DataStorage* config;
111 +    double* frc;
112 +    double* pos;
113 +    double* trq;
114 +    double* A;
115 +    double* unitFrame;
116 +    double* rc;
117 +    
118      //get current snapshot from SimInfo
119 <    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot()
119 >    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
120  
121      //get array pointers
122      config = &(curSnapshot->atomData);
# Line 90 | Line 124 | void ForceManager::calcForces(bool needPotential, bool
124      pos = config->getArrayPointer(DataStorage::dslPosition);
125      trq = config->getArrayPointer(DataStorage::dslTorque);
126      A   = config->getArrayPointer(DataStorage::dslAmat);
127 <    u_l = config->getArrayPointer(DataStorage::dslUnitVector);
127 >    unitFrame = config->getArrayPointer(DataStorage::dslUnitFrame);
128  
129      //calculate the center of mass of cutoff group
130 +    SimInfo::MoleculeIterator mi;
131 +    Molecule* mol;
132 +    Molecule::CutoffGroupIterator ci;
133 +    CutoffGroup* cg;
134 +    Vector3d com;
135 +    std::vector<Vector3d> rcGroup;
136 +    
137      if(info_->getNCutoffGroups() > 0){
138  
139      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
140          for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
141 <            cg->getCOM(com)
141 >            cg->getCOM(com);
142              rcGroup.push_back(com);
143          }
144      }// end for (mol)
145        
146 <        rc = rcGroup[0]->getArrayPointer();
147 <    }
107 <    else{
146 >        rc = rcGroup[0].getArrayPointer();
147 >    } else {
148          // center of mass of the group is the same as position of the atom  if cutoff group does not exist
149          rc = pos;
150      }
151    
112 –
152      //initialize data before passing to fortran
153 <    longRangePotential = 0.0;
154 <    for (int i=0; i<9; i++) {
155 <        tau[i] = 0.0;
156 <    }    
157 <    passedCalcPot = (short int)needPotential;
119 <    passedCalcStress = (short int)needStress;
120 <    isError = 0;
153 >    double longRangePotential = 0.0;
154 >    Mat3x3d tau;
155 >    short int passedCalcPot = needPotential;
156 >    short int passedCalcStress = needStress;
157 >    int isError = 0;
158  
122 – #ifdef PROFILE
123 –    startProfile(pro8);
124 – #endif
125 –
159      doForceLoop( pos,
160              rc,
161              A,
162 <            u_l,
162 >            unitFrame,
163              frc,
164              trq,
165 <            tau,
165 >            tau.getArrayPointer(),
166              &longRangePotential,
167              &passedCalcPot,
168              &passedCalcStress,
169              &isError );
170  
138 – #ifdef PROFILE
139 –    endProfile(pro8);
140 – #endif
141 –
171      if( isError ){
172          sprintf( painCave.errMsg,
173               "Error returned from the fortran force calculation.\n" );
# Line 147 | Line 176 | void ForceManager::calcForces(bool needPotential, bool
176      }
177  
178      //store the tau and long range potential    
179 <    curSnapshot->statData[Stats::LONGRANGEPOTENTIAL] = longRangePotential;
179 >    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential;
180      curSnapshot->statData[Stats::TAUXX] = tau[0];
181      curSnapshot->statData[Stats::TAUXY] = tau[1];
182      curSnapshot->statData[Stats::TAUXZ] = tau[2];
# Line 157 | Line 186 | void ForceManager::calcForces(bool needPotential, bool
186      curSnapshot->statData[Stats::TAUZX] = tau[6];
187      curSnapshot->statData[Stats::TAUZY] = tau[7];
188      curSnapshot->statData[Stats::TAUZZ] = tau[8];
189 + }
190  
191 +
192 + void ForceManager::postCalculation() {
193 +    SimInfo::MoleculeIterator mi;
194 +    Molecule* mol;
195 +    Molecule::RigidBodyIterator rbIter;
196 +    RigidBody* rb;
197 +    
198      // collect the atomic forces onto rigid bodies
199      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
200 <        mol->atoms2rigidBodies();
200 >        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
201 >            rb->calcForcesAndTorques();
202 >        }
203      }
204  
166 – #ifdef IS_MPI
167 –    sprintf( checkPointMsg,
168 –       "returned from the force calculation.\n" );
169 –    MPIcheckPoint();
170 – #endif // is_mpi
171 –
205   }
206  
174 –
207   } //end namespace oopse

Diff Legend

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