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 507 by gezelter, Fri Apr 15 22:04:00 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 – #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 > namespace oopse {
54 >  
55 >  Velocitizer::Velocitizer(SimInfo* info) : info_(info) {
56 >    
57 >    int seedValue;
58 >    Globals * simParams = info->getSimParams();
59 >    
60 > #ifndef IS_MPI
61 >    if (simParams->haveSeed()) {
62 >      seedValue = simParams->getSeed();
63 >      randNumGen_ = new SeqRandNumGen(seedValue);
64 >    }else {
65 >      randNumGen_ = new SeqRandNumGen();
66 >    }    
67 > #else
68 >    if (simParams->haveSeed()) {
69 >      seedValue = simParams->getSeed();
70 >      randNumGen_ = new ParallelRandNumGen(seedValue);
71 >    }else {
72 >      randNumGen_ = new ParallelRandNumGen();
73 >    }    
74 > #endif
75 >  }
76 >  
77 >  Velocitizer::~Velocitizer() {
78 >    delete randNumGen_;
79 >  }
80 >  
81 >  void Velocitizer::velocitize(double temperature) {
82      Vector3d aVel;
83      Vector3d aJ;
84      Mat3x3d I;
# Line 59 | Line 91 | void Velocitizer::velocitize(double temperature) {
91      const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
92      double av2;
93      double kebar;
94 <
94 >    
95      SimInfo::MoleculeIterator i;
96      Molecule::IntegrableObjectIterator j;
97      Molecule * mol;
98      StuntDouble * integrableObject;
67 –
99      
100 < #ifndef IS_MPI
101 <    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 <
100 >    
101 >    
102      kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf());
103 <
103 >    
104      for( mol = info_->beginMolecule(i); mol != NULL;
105 <        mol = info_->nextMolecule(i) ) {
106 <        for( integrableObject = mol->beginIntegrableObject(j);
107 <            integrableObject != NULL;
108 <            integrableObject = mol->nextIntegrableObject(j) ) {
109 <
110 <            // uses equipartition theory to solve for vbar in angstrom/fs
111 <
112 <            av2 = 2.0 * kebar / integrableObject->getMass();
113 <            vbar = sqrt(av2);
114 <
115 <            // picks random velocities from a gaussian distribution
116 <            // centered on vbar
117 <
118 <            for( int k = 0; k < 3; k++ ) {
119 <                aVel[k] = vbar * randNumGen.randNorm(0.0, 1.0);
120 <            }
121 <
122 <            integrableObject->setVel(aVel);
123 <
124 <            if (integrableObject->isDirectional()) {
125 <                I = integrableObject->getI();
126 <
127 <                if (integrableObject->isLinear()) {
128 <                    l = integrableObject->linearAxis();
129 <                    m = (l + 1) % 3;
130 <                    n = (l + 2) % 3;
131 <
132 <                    aJ[l] = 0.0;
133 <                    vbar = sqrt(2.0 * kebar * I(m, m));
134 <                    aJ[m] = vbar * randNumGen.randNorm(0.0, 1.0);
135 <                    vbar = sqrt(2.0 * kebar * I(n, n));
136 <                    aJ[n] = vbar * randNumGen.randNorm(0.0, 1.0);
137 <                } else {
138 <                    for( int k = 0; k < 3; k++ ) {
139 <                        vbar = sqrt(2.0 * kebar * I(k, k));
140 <                        aJ[k] = vbar * randNumGen.randNorm(0.0, 1.0);
141 <                    }
142 <                } // else isLinear
143 <
144 <                integrableObject->setJ(aJ);
145 <            }     //isDirectional
146 <        }
105 >         mol = info_->nextMolecule(i) ) {
106 >      for( integrableObject = mol->beginIntegrableObject(j);
107 >           integrableObject != NULL;
108 >           integrableObject = mol->nextIntegrableObject(j) ) {
109 >        
110 >        // uses equipartition theory to solve for vbar in angstrom/fs
111 >        
112 >        av2 = 2.0 * kebar / integrableObject->getMass();
113 >        vbar = sqrt(av2);
114 >        
115 >        // picks random velocities from a gaussian distribution
116 >        // centered on vbar
117 >        
118 >        for( int k = 0; k < 3; k++ ) {
119 >          aVel[k] = vbar * randNumGen_->randNorm(0.0, 1.0);
120 >        }
121 >        
122 >        integrableObject->setVel(aVel);
123 >        
124 >        if (integrableObject->isDirectional()) {
125 >          I = integrableObject->getI();
126 >          
127 >          if (integrableObject->isLinear()) {
128 >            l = integrableObject->linearAxis();
129 >            m = (l + 1) % 3;
130 >            n = (l + 2) % 3;
131 >            
132 >            aJ[l] = 0.0;
133 >            vbar = sqrt(2.0 * kebar * I(m, m));
134 >            aJ[m] = vbar * randNumGen_->randNorm(0.0, 1.0);
135 >            vbar = sqrt(2.0 * kebar * I(n, n));
136 >            aJ[n] = vbar * randNumGen_->randNorm(0.0, 1.0);
137 >          } else {
138 >            for( int k = 0; k < 3; k++ ) {
139 >              vbar = sqrt(2.0 * kebar * I(k, k));
140 >              aJ[k] = vbar *randNumGen_->randNorm(0.0, 1.0);
141 >            }
142 >          } // else isLinear
143 >          
144 >          integrableObject->setJ(aJ);
145 >        }     //isDirectional
146 >      }
147      }             //end for (mol = beginMolecule(i); ...)
148 <
149 <
150 <
148 >    
149 >    
150 >    
151      removeComDrift();
152 <
153 < }
154 <
155 <
156 <
157 < void Velocitizer::removeComDrift() {
152 >    
153 >  }
154 >  
155 >  
156 >  
157 >  void Velocitizer::removeComDrift() {
158      // Get the Center of Mass drift velocity.
159      Vector3d vdrift = info_->getComVel();
160      
# Line 141 | Line 166 | void Velocitizer::removeComDrift() {
166      //  Corrects for the center of mass drift.
167      // sums all the momentum and divides by total mass.
168      for( mol = info_->beginMolecule(i); mol != NULL;
169 <        mol = info_->nextMolecule(i) ) {
170 <        for( integrableObject = mol->beginIntegrableObject(j);
171 <            integrableObject != NULL;
172 <            integrableObject = mol->nextIntegrableObject(j) ) {
173 <            integrableObject->setVel(integrableObject->getVel() - vdrift);
174 <        }
169 >         mol = info_->nextMolecule(i) ) {
170 >      for( integrableObject = mol->beginIntegrableObject(j);
171 >           integrableObject != NULL;
172 >           integrableObject = mol->nextIntegrableObject(j) ) {
173 >        integrableObject->setVel(integrableObject->getVel() - vdrift);
174 >      }
175      }
176 <
176 >    
177 >  }
178 >  
179   }
153 –
154 – }

Diff Legend

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