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 1737 by tim, Wed Nov 10 22:50:03 2004 UTC vs.
Revision 1738 by tim, Sat Nov 13 05:08:12 2004 UTC

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

Diff Legend

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