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 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1079 by gezelter, Thu Oct 19 15:57:07 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/randomSPRNG.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 < void Velocitizer::velocitize(double temperature) {
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 55 | 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;
106 <    gaussianSPRNG gaussStream(info_->getSeed());
69 <
106 >        
107      kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf());
71 –
108      for( mol = info_->beginMolecule(i); mol != NULL;
109 <        mol = info_->nextMolecule(i) ) {
110 <        for( integrableObject = mol->beginIntegrableObject(j);
111 <            integrableObject != NULL;
112 <            integrableObject = mol->nextIntegrableObject(j) ) {
113 <
114 <            // uses equipartition theory to solve for vbar in angstrom/fs
115 <
116 <            av2 = 2.0 * kebar / integrableObject->getMass();
117 <            vbar = sqrt(av2);
118 <
119 <            // picks random velocities from a gaussian distribution
120 <            // centered on vbar
121 <
122 <            for( int k = 0; k < 3; k++ ) {
123 <                aVel[k] = vbar * gaussStream.getGaussian();
124 <            }
125 <
126 <            integrableObject->setVel(aVel);
127 <
128 <            if (integrableObject->isDirectional()) {
129 <                I = integrableObject->getI();
130 <
131 <                if (integrableObject->isLinear()) {
132 <                    l = integrableObject->linearAxis();
133 <                    m = (l + 1) % 3;
134 <                    n = (l + 2) % 3;
135 <
136 <                    aJ[l] = 0.0;
137 <                    vbar = sqrt(2.0 * kebar * I(m, m));
138 <                    aJ[m] = vbar * gaussStream.getGaussian();
139 <                    vbar = sqrt(2.0 * kebar * I(n, n));
140 <                    aJ[n] = vbar * gaussStream.getGaussian();
141 <                } else {
142 <                    for( int k = 0; k < 3; k++ ) {
143 <                        vbar = sqrt(2.0 * kebar * I(k, k));
144 <                        aJ[k] = vbar * gaussStream.getGaussian();
145 <                    }
146 <                } // else isLinear
147 <
148 <                integrableObject->setJ(aJ);
149 <            }     //isDirectional
114 <        }
109 >         mol = info_->nextMolecule(i) ) {
110 >      for( integrableObject = mol->beginIntegrableObject(j);
111 >           integrableObject != NULL;
112 >           integrableObject = mol->nextIntegrableObject(j) ) {
113 >        
114 >        // uses equipartition theory to solve for vbar in angstrom/fs
115 >        
116 >        av2 = 2.0 * kebar / integrableObject->getMass();
117 >        vbar = sqrt(av2);
118 >        
119 >        // picks random velocities from a gaussian distribution
120 >        // centered on vbar
121 >        
122 >        for( int k = 0; k < 3; k++ ) {
123 >          aVel[k] = vbar * randNumGen_->randNorm(0.0, 1.0);
124 >        }
125 >        integrableObject->setVel(aVel);
126 >        
127 >        if (integrableObject->isDirectional()) {
128 >          I = integrableObject->getI();
129 >          
130 >          if (integrableObject->isLinear()) {
131 >            l = integrableObject->linearAxis();
132 >            m = (l + 1) % 3;
133 >            n = (l + 2) % 3;
134 >            
135 >            aJ[l] = 0.0;
136 >            vbar = sqrt(2.0 * kebar * I(m, m));
137 >            aJ[m] = vbar * randNumGen_->randNorm(0.0, 1.0);
138 >            vbar = sqrt(2.0 * kebar * I(n, n));
139 >            aJ[n] = vbar * randNumGen_->randNorm(0.0, 1.0);
140 >          } else {
141 >            for( int k = 0; k < 3; k++ ) {
142 >              vbar = sqrt(2.0 * kebar * I(k, k));
143 >              aJ[k] = vbar *randNumGen_->randNorm(0.0, 1.0);
144 >            }
145 >          } // else isLinear
146 >          
147 >          integrableObject->setJ(aJ);
148 >        }     //isDirectional
149 >      }
150      }             //end for (mol = beginMolecule(i); ...)
151 <
152 <
153 <
151 >    
152 >    
153 >    
154      removeComDrift();
155 <
156 < }
157 <
158 <
159 <
160 < void Velocitizer::removeComDrift() {
155 >    // Remove angular drift if we are not using periodic boundary conditions.
156 >    if(!simParams->getUsePeriodicBoundaryConditions()) removeAngularDrift();
157 >    
158 >  }
159 >  
160 >  
161 >  
162 >  void Velocitizer::removeComDrift() {
163      // Get the Center of Mass drift velocity.
164      Vector3d vdrift = info_->getComVel();
165      
# Line 134 | Line 171 | void Velocitizer::removeComDrift() {
171      //  Corrects for the center of mass drift.
172      // sums all the momentum and divides by total mass.
173      for( mol = info_->beginMolecule(i); mol != NULL;
174 <        mol = info_->nextMolecule(i) ) {
175 <        for( integrableObject = mol->beginIntegrableObject(j);
176 <            integrableObject != NULL;
177 <            integrableObject = mol->nextIntegrableObject(j) ) {
178 <            integrableObject->setVel(integrableObject->getVel() - vdrift);
179 <        }
174 >         mol = info_->nextMolecule(i) ) {
175 >      for( integrableObject = mol->beginIntegrableObject(j);
176 >           integrableObject != NULL;
177 >           integrableObject = mol->nextIntegrableObject(j) ) {
178 >        integrableObject->setVel(integrableObject->getVel() - vdrift);
179 >      }
180      }
181 <
181 >    
182 >  }
183 >  
184 >  
185 >  void Velocitizer::removeAngularDrift() {
186 >    // Get the Center of Mass drift velocity.
187 >      
188 >    Vector3d vdrift;
189 >    Vector3d com;
190 >      
191 >    info_->getComAll(com,vdrift);
192 >        
193 >    Mat3x3d inertiaTensor;
194 >    Vector3d angularMomentum;
195 >    Vector3d omega;
196 >      
197 >      
198 >      
199 >    info_->getInertiaTensor(inertiaTensor,angularMomentum);
200 >    // We now need the inverse of the inertia tensor.
201 >    /*
202 >    std::cerr << "Angular Momentum before is "
203 >              << angularMomentum <<  std::endl;
204 >    std::cerr << "Inertia Tensor before is "
205 >              << inertiaTensor <<  std::endl;
206 >    */  
207 >    inertiaTensor =inertiaTensor.inverse();
208 >    /*
209 >    std::cerr << "Inertia Tensor after inverse is "
210 >              << inertiaTensor <<  std::endl;
211 >    */
212 >    omega = inertiaTensor*angularMomentum;
213 >      
214 >    SimInfo::MoleculeIterator i;
215 >    Molecule::IntegrableObjectIterator j;
216 >    Molecule * mol;
217 >    StuntDouble * integrableObject;
218 >    Vector3d tempComPos;
219 >      
220 >    //  Corrects for the center of mass angular drift.
221 >    // sums all the angular momentum and divides by total mass.
222 >    for( mol = info_->beginMolecule(i); mol != NULL;
223 >         mol = info_->nextMolecule(i) ) {
224 >      for( integrableObject = mol->beginIntegrableObject(j);
225 >           integrableObject != NULL;
226 >           integrableObject = mol->nextIntegrableObject(j) ) {
227 >        tempComPos = integrableObject->getPos()-com;
228 >        integrableObject->setVel((integrableObject->getVel() - vdrift)-cross(omega,tempComPos));
229 >      }
230 >    }
231 >      
232 >    angularMomentum = info_->getAngularMomentum();
233 >    /*
234 >    std::cerr << "Angular Momentum after is "
235 >              << angularMomentum <<  std::endl;
236 >    */
237 >      
238 >  }
239 >  
240 >  
241 >  
242 >  
243   }
146 –
147 – }

Diff Legend

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