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 381 by tim, Tue Mar 1 14:45:45 2005 UTC vs.
Revision 963 by tim, Wed May 17 21:51:42 2006 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 – #include "math/MersenneTwister.hpp"
47 – namespace oopse {
46  
47 < void Velocitizer::velocitize(double temperature) {
47 > #ifndef IS_MPI
48 > #include "math/SeqRandNumGen.hpp"
49 > #else
50 > #include "math/ParallelRandNumGen.hpp"
51 > #endif
52 >
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 >    
65 > #ifndef IS_MPI
66 >    if (simParams->haveSeed()) {
67 >      seedValue = simParams->getSeed();
68 >      randNumGen_ = new SeqRandNumGen(seedValue);
69 >    }else {
70 >      randNumGen_ = new SeqRandNumGen();
71 >    }    
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(RealType temperature) {
87      Vector3d aVel;
88      Vector3d aJ;
89      Mat3x3d I;
# Line 54 | Line 91 | void Velocitizer::velocitize(double temperature) {
91      int m;
92      int n;
93      Vector3d vdrift;
94 <    double vbar;
94 >    RealType vbar;
95      /**@todo refactory kb */
96 <    const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
97 <    double av2;
98 <    double kebar;
99 <
96 >    const RealType kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
97 >    RealType av2;
98 >    RealType kebar;
99 >    
100 >    Globals * simParams = info_->getSimParams();
101 >    
102      SimInfo::MoleculeIterator i;
103      Molecule::IntegrableObjectIterator j;
104      Molecule * mol;
105      StuntDouble * integrableObject;
67 –
106      
107 < #ifndef IS_MPI
108 <    MTRand randNumGen(info_->getSeed());
71 < #else
72 <    int nProcessors;
73 <    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
74 <    MTRand randNumGen(info_->getSeed(), nProcessors, worldRank);
75 < #endif
76 <
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->getUsePeriodicBoundaryConditions()) 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 141 | 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 >      /*
214 >       std::cerr << "Inertia Tensor after inverse is "
215 >                << inertiaTensor <<  std::endl;
216 >      */
217 >      omega = inertiaTensor*angularMomentum;
218 >      
219 >      SimInfo::MoleculeIterator i;
220 >      Molecule::IntegrableObjectIterator j;
221 >      Molecule * mol;
222 >      StuntDouble * integrableObject;
223 >      Vector3d tempComPos;
224 >      
225 >      //  Corrects for the center of mass angular drift.
226 >      // sums all the angular momentum and divides by total mass.
227 >      for( mol = info_->beginMolecule(i); mol != NULL;
228 >           mol = info_->nextMolecule(i) ) {
229 >         for( integrableObject = mol->beginIntegrableObject(j);
230 >              integrableObject != NULL;
231 >              integrableObject = mol->nextIntegrableObject(j) ) {
232 >            tempComPos = integrableObject->getPos()-com;
233 >            integrableObject->setVel((integrableObject->getVel() - vdrift)-cross(omega,tempComPos));
234 >         }
235 >      }
236 >      
237 >      angularMomentum = info_->getAngularMomentum();
238 >      /*
239 >      std::cerr << "Angular Momentum after is "
240 >         << angularMomentum <<  std::endl;
241 >      */
242 >      
243 >   }
244 >  
245 >  
246 >  
247 >  
248   }
153 –
154 – }

Diff Legend

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