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

Comparing trunk/src/integrators/Velocitizer.cpp (file contents):
Revision 388 by tim, Tue Mar 1 23:02:33 2005 UTC vs.
Revision 557 by chuckv, Tue May 31 22:31:54 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 MTRand(seedValue);
67 >      seedValue = simParams->getSeed();
68 >      randNumGen_ = new SeqRandNumGen(seedValue);
69      }else {
70 <        randNumGen_ = new MTRand();
70 >      randNumGen_ = new SeqRandNumGen();
71      }    
72   #else
73      if (simParams->haveSeed()) {
74 <        seedValue = simParams->getSeed();
75 <        randNumGen_ = new ParallelRandNumGen(seedValue);
74 >      seedValue = simParams->getSeed();
75 >      randNumGen_ = new ParallelRandNumGen(seedValue);
76      }else {
77 <        randNumGen_ = new ParallelRandNumGen();
77 >      randNumGen_ = new ParallelRandNumGen();
78      }    
79   #endif
80 < }
81 <
82 < Velocitizer::~Velocitizer() {
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 85 | 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 160 | 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 +      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 +      std::cerr << "Angular Momentum after is "
238 +         << angularMomentum <<  std::endl;
239  
240 +      
241 +   }
242 +  
243 +  
244 +  
245 +  
246   }
172 –
173 – }

Diff Legend

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