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 2252 by chuckv, Mon May 30 14:01:52 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 +    
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  
209 < }
209 >      
210 >      inertiaTensor.inverse();
211 >      
212 >      
213 >      omega = inertiaTensor*angularMomentum;
214 >      
215 >      SimInfo::MoleculeIterator i;
216 >      Molecule::IntegrableObjectIterator j;
217 >      Molecule * mol;
218 >      StuntDouble * integrableObject;
219 >      Vector3d tempComPos;
220 >      
221 >      //  Corrects for the center of mass angular drift.
222 >      // sums all the angular momentum and divides by total mass.
223 >      for( mol = info_->beginMolecule(i); mol != NULL;
224 >           mol = info_->nextMolecule(i) ) {
225 >         for( integrableObject = mol->beginIntegrableObject(j);
226 >              integrableObject != NULL;
227 >              integrableObject = mol->nextIntegrableObject(j) ) {
228 >            tempComPos = integrableObject->getPos()-com;
229 >            integrableObject->setVel((integrableObject->getVel() - vdrift)-cross(omega,tempComPos));
230 >         }
231 >      }
232 >      
233 >      angularMomentum = info_->getAngularMomentum();
234 >      std::cerr << "Angular Momentum after is "
235 >         << angularMomentum <<  std::endl;
236  
237 +      
238 +   }
239 +  
240 +  
241 +  
242 +  
243   }

Diff Legend

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