ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/integrators/Velocitizer.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/integrators/Velocitizer.cpp (file contents):
Revision 2068 by tim, Tue Mar 1 19:11:47 2005 UTC vs.
Revision 2320 by chuckv, Thu Sep 22 00:04:56 2005 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 38 | Line 38
38   * University of Notre Dame has been advised of the possibility of
39   * such damages.
40   */
41 <  
41 >
42   #include "integrators/Velocitizer.hpp"
43   #include "math/SquareMatrix3.hpp"
44   #include "primitives/Molecule.hpp"
45   #include "primitives/StuntDouble.hpp"
46  
47 < namespace oopse {
47 > #ifndef IS_MPI
48 > #include "math/SeqRandNumGen.hpp"
49 > #else
50 > #include "math/ParallelRandNumGen.hpp"
51 > #endif
52  
53 < Velocitizer::Velocitizer(SimInfo* info) {
53 > /* Remove me after testing*/
54 > #include <cstdio>
55 > #include <iostream>
56 > /*End remove me*/
57  
58 + namespace oopse {
59 +  
60 +  Velocitizer::Velocitizer(SimInfo* info) : info_(info) {
61 +    
62      int seedValue;
63      Globals * simParams = info->getSimParams();
64 <
64 >    
65 > #ifndef IS_MPI
66      if (simParams->haveSeed()) {
67 <        seedValue = simParams->getSeed();
68 <        randNumGen_ = new OOPSERandNumGen(seedValue);
67 >      seedValue = simParams->getSeed();
68 >      randNumGen_ = new SeqRandNumGen(seedValue);
69      }else {
70 <        randNumGen_ = new OOPSERandNumGen();
70 >      randNumGen_ = new SeqRandNumGen();
71      }    
72 <
73 < }
74 <
75 < Velocitizer::~Velocitizer() {
72 > #else
73 >    if (simParams->haveSeed()) {
74 >      seedValue = simParams->getSeed();
75 >      randNumGen_ = new ParallelRandNumGen(seedValue);
76 >    }else {
77 >      randNumGen_ = new ParallelRandNumGen();
78 >    }    
79 > #endif
80 >  }
81 >  
82 >  Velocitizer::~Velocitizer() {
83      delete randNumGen_;
84 < }
85 <
86 < void Velocitizer::velocitize(double temperature) {
84 >  }
85 >  
86 >  void Velocitizer::velocitize(double temperature) {
87      Vector3d aVel;
88      Vector3d aJ;
89      Mat3x3d I;
# Line 77 | Line 96 | void Velocitizer::velocitize(double temperature) {
96      const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
97      double av2;
98      double kebar;
99 <
99 >    
100 >    Globals * simParams = info_->getSimParams();
101 >    
102      SimInfo::MoleculeIterator i;
103      Molecule::IntegrableObjectIterator j;
104      Molecule * mol;
105      StuntDouble * integrableObject;
106 <
107 <
108 <
106 >    
107 >    
108 >    
109      kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf());
110 <
110 >    
111      for( mol = info_->beginMolecule(i); mol != NULL;
112 <        mol = info_->nextMolecule(i) ) {
113 <        for( integrableObject = mol->beginIntegrableObject(j);
114 <            integrableObject != NULL;
115 <            integrableObject = mol->nextIntegrableObject(j) ) {
116 <
117 <            // uses equipartition theory to solve for vbar in angstrom/fs
118 <
119 <            av2 = 2.0 * kebar / integrableObject->getMass();
120 <            vbar = sqrt(av2);
121 <
122 <            // picks random velocities from a gaussian distribution
123 <            // centered on vbar
124 <
125 <            for( int k = 0; k < 3; k++ ) {
126 <                aVel[k] = vbar * randNumGen_->randNorm(0.0, 1.0);
127 <            }
128 <
129 <            integrableObject->setVel(aVel);
130 <
131 <            if (integrableObject->isDirectional()) {
132 <                I = integrableObject->getI();
133 <
134 <                if (integrableObject->isLinear()) {
135 <                    l = integrableObject->linearAxis();
136 <                    m = (l + 1) % 3;
137 <                    n = (l + 2) % 3;
138 <
139 <                    aJ[l] = 0.0;
140 <                    vbar = sqrt(2.0 * kebar * I(m, m));
141 <                    aJ[m] = vbar * randNumGen_->randNorm(0.0, 1.0);
142 <                    vbar = sqrt(2.0 * kebar * I(n, n));
143 <                    aJ[n] = vbar * randNumGen_->randNorm(0.0, 1.0);
144 <                } else {
145 <                    for( int k = 0; k < 3; k++ ) {
146 <                        vbar = sqrt(2.0 * kebar * I(k, k));
147 <                        aJ[k] = vbar *randNumGen_->randNorm(0.0, 1.0);
148 <                    }
149 <                } // else isLinear
150 <
151 <                integrableObject->setJ(aJ);
152 <            }     //isDirectional
153 <        }
112 >         mol = info_->nextMolecule(i) ) {
113 >      for( integrableObject = mol->beginIntegrableObject(j);
114 >           integrableObject != NULL;
115 >           integrableObject = mol->nextIntegrableObject(j) ) {
116 >        
117 >        // uses equipartition theory to solve for vbar in angstrom/fs
118 >        
119 >        av2 = 2.0 * kebar / integrableObject->getMass();
120 >        vbar = sqrt(av2);
121 >        
122 >        // picks random velocities from a gaussian distribution
123 >        // centered on vbar
124 >        
125 >        for( int k = 0; k < 3; k++ ) {
126 >          aVel[k] = vbar * randNumGen_->randNorm(0.0, 1.0);
127 >        }
128 >        
129 >        integrableObject->setVel(aVel);
130 >        
131 >        if (integrableObject->isDirectional()) {
132 >          I = integrableObject->getI();
133 >          
134 >          if (integrableObject->isLinear()) {
135 >            l = integrableObject->linearAxis();
136 >            m = (l + 1) % 3;
137 >            n = (l + 2) % 3;
138 >            
139 >            aJ[l] = 0.0;
140 >            vbar = sqrt(2.0 * kebar * I(m, m));
141 >            aJ[m] = vbar * randNumGen_->randNorm(0.0, 1.0);
142 >            vbar = sqrt(2.0 * kebar * I(n, n));
143 >            aJ[n] = vbar * randNumGen_->randNorm(0.0, 1.0);
144 >          } else {
145 >            for( int k = 0; k < 3; k++ ) {
146 >              vbar = sqrt(2.0 * kebar * I(k, k));
147 >              aJ[k] = vbar *randNumGen_->randNorm(0.0, 1.0);
148 >            }
149 >          } // else isLinear
150 >          
151 >          integrableObject->setJ(aJ);
152 >        }     //isDirectional
153 >      }
154      }             //end for (mol = beginMolecule(i); ...)
155 <
156 <
157 <
155 >    
156 >    
157 >    
158      removeComDrift();
159 <
160 < }
161 <
162 <
163 <
164 < void Velocitizer::removeComDrift() {
159 >    // Remove angular drift if we are not using periodic boundary conditions.
160 >    if(!simParams->getPBC()) removeAngularDrift();
161 >    
162 >  }
163 >  
164 >  
165 >  
166 >  void Velocitizer::removeComDrift() {
167      // Get the Center of Mass drift velocity.
168      Vector3d vdrift = info_->getComVel();
169      
# Line 152 | Line 175 | void Velocitizer::removeComDrift() {
175      //  Corrects for the center of mass drift.
176      // sums all the momentum and divides by total mass.
177      for( mol = info_->beginMolecule(i); mol != NULL;
178 <        mol = info_->nextMolecule(i) ) {
179 <        for( integrableObject = mol->beginIntegrableObject(j);
180 <            integrableObject != NULL;
181 <            integrableObject = mol->nextIntegrableObject(j) ) {
182 <            integrableObject->setVel(integrableObject->getVel() - vdrift);
183 <        }
178 >         mol = info_->nextMolecule(i) ) {
179 >      for( integrableObject = mol->beginIntegrableObject(j);
180 >           integrableObject != NULL;
181 >           integrableObject = mol->nextIntegrableObject(j) ) {
182 >        integrableObject->setVel(integrableObject->getVel() - vdrift);
183 >      }
184      }
185 <
185 >    
186 >  }
187 >  
188 >  
189 >   void Velocitizer::removeAngularDrift() {
190 >      // Get the Center of Mass drift velocity.
191 >      
192 >      Vector3d vdrift;
193 >      Vector3d com;
194 >      
195 >      info_->getComAll(com,vdrift);
196 >        
197 >      Mat3x3d inertiaTensor;
198 >      Vector3d angularMomentum;
199 >      Vector3d omega;
200 >      
201 >      
202 >      
203 >      info_->getInertiaTensor(inertiaTensor,angularMomentum);
204 >      // We now need the inverse of the inertia tensor.
205 >      /*      
206 >      std::cerr << "Angular Momentum before is "
207 >                << angularMomentum <<  std::endl;
208 >      std::cerr << "Inertia Tensor before is "
209 >                << inertiaTensor <<  std::endl;
210 >      */
211 >      
212 >      inertiaTensor =inertiaTensor.inverse();
213 >       std::cerr << "Inertia Tensor after inverse is "
214 >                << inertiaTensor <<  std::endl;
215 >      
216 >      omega = inertiaTensor*angularMomentum;
217 >      
218 >      SimInfo::MoleculeIterator i;
219 >      Molecule::IntegrableObjectIterator j;
220 >      Molecule * mol;
221 >      StuntDouble * integrableObject;
222 >      Vector3d tempComPos;
223 >      
224 >      //  Corrects for the center of mass angular drift.
225 >      // sums all the angular momentum and divides by total mass.
226 >      for( mol = info_->beginMolecule(i); mol != NULL;
227 >           mol = info_->nextMolecule(i) ) {
228 >         for( integrableObject = mol->beginIntegrableObject(j);
229 >              integrableObject != NULL;
230 >              integrableObject = mol->nextIntegrableObject(j) ) {
231 >            tempComPos = integrableObject->getPos()-com;
232 >            integrableObject->setVel((integrableObject->getVel() - vdrift)-cross(omega,tempComPos));
233 >         }
234 >      }
235 >      
236 >      angularMomentum = info_->getAngularMomentum();
237 >      /*
238 >      std::cerr << "Angular Momentum after is "
239 >         << angularMomentum <<  std::endl;
240 >      */
241 >      
242 >   }
243 >  
244 >  
245 >  
246 >  
247   }
164 –
165 – }

Diff Legend

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