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 1390 by gezelter, Wed Nov 25 20:02:06 2009 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 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
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 > /*
55 > #include <cstdio>
56 > #include <iostream>
57 > */
58 > /*End remove me*/
59 >
60 > namespace OpenMD {
61 >  
62 >  Velocitizer::Velocitizer(SimInfo* info) : info_(info) {
63 >    
64 >    int seedValue;
65 >    Globals * simParams = info->getSimParams();
66 >    
67 > #ifndef IS_MPI
68 >    if (simParams->haveSeed()) {
69 >      seedValue = simParams->getSeed();
70 >      randNumGen_ = new SeqRandNumGen(seedValue);
71 >    }else {
72 >      randNumGen_ = new SeqRandNumGen();
73 >    }    
74 > #else
75 >    if (simParams->haveSeed()) {
76 >      seedValue = simParams->getSeed();
77 >      randNumGen_ = new ParallelRandNumGen(seedValue);
78 >    }else {
79 >      randNumGen_ = new ParallelRandNumGen();
80 >    }    
81 > #endif
82 >  }
83 >  
84 >  Velocitizer::~Velocitizer() {
85 >    delete randNumGen_;
86 >  }
87 >  
88 >  void Velocitizer::velocitize(RealType temperature) {
89      Vector3d aVel;
90      Vector3d aJ;
91      Mat3x3d I;
# Line 55 | Line 93 | void Velocitizer::velocitize(double temperature) {
93      int m;
94      int n;
95      Vector3d vdrift;
96 <    double vbar;
96 >    RealType vbar;
97      /**@todo refactory kb */
98 <    const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
99 <    double av2;
100 <    double kebar;
101 <
98 >    const RealType kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
99 >    RealType av2;
100 >    RealType kebar;
101 >    
102 >    Globals * simParams = info_->getSimParams();
103 >    
104      SimInfo::MoleculeIterator i;
105      Molecule::IntegrableObjectIterator j;
106      Molecule * mol;
107      StuntDouble * integrableObject;
68 –    gaussianSPRNG gaussStream(info_->getSeed());
108  
109      kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf());
71 –
110      for( mol = info_->beginMolecule(i); mol != NULL;
111 <        mol = info_->nextMolecule(i) ) {
112 <        for( integrableObject = mol->beginIntegrableObject(j);
113 <            integrableObject != NULL;
114 <            integrableObject = mol->nextIntegrableObject(j) ) {
115 <
116 <            // uses equipartition theory to solve for vbar in angstrom/fs
117 <
118 <            av2 = 2.0 * kebar / integrableObject->getMass();
119 <            vbar = sqrt(av2);
120 <
121 <            // picks random velocities from a gaussian distribution
122 <            // centered on vbar
123 <
124 <            for( int k = 0; k < 3; k++ ) {
125 <                aVel[k] = vbar * gaussStream.getGaussian();
126 <            }
127 <
128 <            integrableObject->setVel(aVel);
129 <
130 <            if (integrableObject->isDirectional()) {
131 <                I = integrableObject->getI();
132 <
133 <                if (integrableObject->isLinear()) {
134 <                    l = integrableObject->linearAxis();
135 <                    m = (l + 1) % 3;
136 <                    n = (l + 2) % 3;
137 <
138 <                    aJ[l] = 0.0;
139 <                    vbar = sqrt(2.0 * kebar * I(m, m));
140 <                    aJ[m] = vbar * gaussStream.getGaussian();
141 <                    vbar = sqrt(2.0 * kebar * I(n, n));
142 <                    aJ[n] = vbar * gaussStream.getGaussian();
143 <                } else {
144 <                    for( int k = 0; k < 3; k++ ) {
145 <                        vbar = sqrt(2.0 * kebar * I(k, k));
146 <                        aJ[k] = vbar * gaussStream.getGaussian();
147 <                    }
148 <                } // else isLinear
149 <
150 <                integrableObject->setJ(aJ);
151 <            }     //isDirectional
114 <        }
111 >         mol = info_->nextMolecule(i) ) {
112 >      for( integrableObject = mol->beginIntegrableObject(j);
113 >           integrableObject != NULL;
114 >           integrableObject = mol->nextIntegrableObject(j) ) {
115 >        
116 >        // uses equipartition theory to solve for vbar in angstrom/fs
117 >        
118 >        av2 = 2.0 * kebar / integrableObject->getMass();
119 >        vbar = sqrt(av2);
120 >        
121 >        // picks random velocities from a gaussian distribution
122 >        // centered on vbar
123 >        
124 >        for( int k = 0; k < 3; k++ ) {
125 >          aVel[k] = vbar * randNumGen_->randNorm(0.0, 1.0);
126 >        }
127 >        integrableObject->setVel(aVel);
128 >        
129 >        if (integrableObject->isDirectional()) {
130 >          I = integrableObject->getI();
131 >          
132 >          if (integrableObject->isLinear()) {
133 >            l = integrableObject->linearAxis();
134 >            m = (l + 1) % 3;
135 >            n = (l + 2) % 3;
136 >            
137 >            aJ[l] = 0.0;
138 >            vbar = sqrt(2.0 * kebar * I(m, m));
139 >            aJ[m] = vbar * randNumGen_->randNorm(0.0, 1.0);
140 >            vbar = sqrt(2.0 * kebar * I(n, n));
141 >            aJ[n] = vbar * randNumGen_->randNorm(0.0, 1.0);
142 >          } else {
143 >            for( int k = 0; k < 3; k++ ) {
144 >              vbar = sqrt(2.0 * kebar * I(k, k));
145 >              aJ[k] = vbar *randNumGen_->randNorm(0.0, 1.0);
146 >            }
147 >          } // else isLinear
148 >          
149 >          integrableObject->setJ(aJ);
150 >        }     //isDirectional
151 >      }
152      }             //end for (mol = beginMolecule(i); ...)
153 <
154 <
155 <
153 >    
154 >    
155 >    
156      removeComDrift();
157 <
158 < }
159 <
160 <
161 <
162 < void Velocitizer::removeComDrift() {
157 >    // Remove angular drift if we are not using periodic boundary conditions.
158 >    if(!simParams->getUsePeriodicBoundaryConditions()) removeAngularDrift();
159 >    
160 >  }
161 >  
162 >  
163 >  
164 >  void Velocitizer::removeComDrift() {
165      // Get the Center of Mass drift velocity.
166      Vector3d vdrift = info_->getComVel();
167      
# Line 134 | Line 173 | void Velocitizer::removeComDrift() {
173      //  Corrects for the center of mass drift.
174      // sums all the momentum and divides by total mass.
175      for( mol = info_->beginMolecule(i); mol != NULL;
176 <        mol = info_->nextMolecule(i) ) {
177 <        for( integrableObject = mol->beginIntegrableObject(j);
178 <            integrableObject != NULL;
179 <            integrableObject = mol->nextIntegrableObject(j) ) {
180 <            integrableObject->setVel(integrableObject->getVel() - vdrift);
181 <        }
176 >         mol = info_->nextMolecule(i) ) {
177 >      for( integrableObject = mol->beginIntegrableObject(j);
178 >           integrableObject != NULL;
179 >           integrableObject = mol->nextIntegrableObject(j) ) {
180 >        integrableObject->setVel(integrableObject->getVel() - vdrift);
181 >      }
182      }
183 <
183 >    
184 >  }
185 >  
186 >  
187 >  void Velocitizer::removeAngularDrift() {
188 >    // Get the Center of Mass drift velocity.
189 >      
190 >    Vector3d vdrift;
191 >    Vector3d com;
192 >      
193 >    info_->getComAll(com,vdrift);
194 >        
195 >    Mat3x3d inertiaTensor;
196 >    Vector3d angularMomentum;
197 >    Vector3d omega;
198 >      
199 >      
200 >      
201 >    info_->getInertiaTensor(inertiaTensor,angularMomentum);
202 >    // We now need the inverse of the inertia tensor.
203 >    /*
204 >    std::cerr << "Angular Momentum before is "
205 >              << angularMomentum <<  std::endl;
206 >    std::cerr << "Inertia Tensor before is "
207 >              << inertiaTensor <<  std::endl;
208 >    */
209 >    inertiaTensor =inertiaTensor.inverse();
210 >    /*
211 >    std::cerr << "Inertia Tensor after inverse is "
212 >              << inertiaTensor <<  std::endl;
213 >    */
214 >    omega = inertiaTensor*angularMomentum;
215 >      
216 >    SimInfo::MoleculeIterator i;
217 >    Molecule::IntegrableObjectIterator j;
218 >    Molecule * mol;
219 >    StuntDouble * integrableObject;
220 >    Vector3d tempComPos;
221 >      
222 >    //  Corrects for the center of mass angular drift.
223 >    // sums all the angular momentum and divides by total mass.
224 >    for( mol = info_->beginMolecule(i); mol != NULL;
225 >         mol = info_->nextMolecule(i) ) {
226 >      for( integrableObject = mol->beginIntegrableObject(j);
227 >           integrableObject != NULL;
228 >           integrableObject = mol->nextIntegrableObject(j) ) {
229 >        tempComPos = integrableObject->getPos()-com;
230 >        integrableObject->setVel((integrableObject->getVel() - vdrift)-cross(omega,tempComPos));
231 >      }
232 >    }
233 >      
234 >    angularMomentum = info_->getAngularMomentum();
235 >    /*
236 >    std::cerr << "Angular Momentum after is "
237 >              << angularMomentum <<  std::endl;
238 >    */
239 >  }
240 >  
241 >  
242 >  
243 >  
244   }
146 –
147 – }

Diff Legend

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