ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/integrators/VelocityVerletIntegrator.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/integrators/VelocityVerletIntegrator.cpp (file contents):
Revision 1722 by tim, Tue Nov 9 23:11:39 2004 UTC vs.
Revision 1725 by tim, Wed Nov 10 22:01:06 2004 UTC

# Line 22 | Line 22
22   * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
23   *
24   */
25 –
26 – /**
27 –  * @file VelocityVerletIntegrator.cpp
28 –  * @author tlin
29 –  * @date 11/09/2004
30 –  * @time 16:16am
31 –  * @version 1.0
32 –  */
25  
26 + /**
27 + * @file VelocityVerletIntegrator.cpp
28 + * @author tlin
29 + * @date 11/09/2004
30 + * @time 16:16am
31 + * @version 1.0
32 + */
33 +
34   #include "integrators/VelocityVerletIntegrator.hpp"
35  
36   namespace oopse {
37 + VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) { }
38  
39 + VelocityVerletIntegrator::~VelocityVerletIntegrator() { }
40  
41 < VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo* info){
41 > void VelocityVerletIntegrator::integrate() {
42 >    double runTime = info->run_time;
43  
44 < }
44 >    double sampleTime = info->sampleTime;
45 >    double statusTime = info->statusTime;
46 >    double thermalTime = info->thermalTime;
47 >    double resetTime = info->resetTime;
48  
49 < VelocityVerletIntegrator::~VelocityVerletIntegrator() {
49 >    double difference;
50 >    double currSample;
51 >    double currThermal;
52 >    double currStatus;
53 >    double currReset;
54  
55 < }
55 >    int calcPot,
56 >    calcStress;
57  
58 < void VelocityVerletIntegrator::integrate() {
58 >    tStats = new Thermo(info);
59 >    statOut = new StatWriter(info);
60 >    dumpOut = new DumpWriter(info);
61  
62 < }
62 >    atoms = info->atoms;
63  
64 < void VelocityVerletIntegrator::integrateStep() {
64 >    fullStep_ = info->dt;
65 >    halfStep_ = 0.5 * fullStep_;
66  
67 < }
67 >    readyCheck();
68  
69 < void VelocityVerletIntegrator::moveA() {
69 >    // remove center of mass drift velocity (in case we passed in a configuration
70 >    // that was drifting
71 >    tStats->removeCOMdrift();
72  
73 < }
73 >    // initialize the forces before the first step
74  
75 < void VelocityVerletIntegrator::moveB() {
75 >    calcForce(1, 1);
76  
77 +    //execute constraint algorithm to make sure at the very beginning the system is constrained  
78 +    if (nConstrained) {
79 +        constrainA();
80 +        calcForce(1, 1);
81 +        constrainB();
82 +    }
83 +
84 +    if (info->setTemp) {
85 +        thermalize();
86 +    }
87 +
88 +    calcPot = 0;
89 +    calcStress = 0;
90 +    currSample = sampleTime + info->getTime();
91 +    currThermal = thermalTime + info->getTime();
92 +    currStatus = statusTime + info->getTime();
93 +    currReset = resetTime + info->getTime();
94 +
95 +    dumpOut->writeDump(info->getTime());
96 +    statOut->writeStat(info->getTime());
97 +
98 + #ifdef IS_MPI
99 +
100 +    strcpy(checkPointMsg, "The integrator is ready to go.");
101 +    MPIcheckPoint();
102 +
103 + #endif // is_mpi
104 +
105 +    while (info->getTime() < runTime && !stopIntegrator()) {
106 +        difference = info->getTime() + fullStep_ - currStatus;
107 +
108 +        if (difference > 0 || fabs(difference) < 1e - 4) {
109 +            calcPot = 1;
110 +            calcStress = 1;
111 +        }
112 +
113 + #ifdef PROFILE
114 +
115 +        startProfile(pro1);
116 +
117 + #endif
118 +
119 +        integrateStep(calcPot, calcStress);
120 +
121 + #ifdef PROFILE
122 +
123 +        endProfile(pro1);
124 +
125 +        startProfile(pro2);
126 +
127 + #endif // profile
128 +
129 +        info->incrTime(fullStep_);
130 +
131 +        if (info->setTemp) {
132 +            if (info->getTime() >= currThermal) {
133 +                thermalize();
134 +                currThermal += thermalTime;
135 +            }
136 +        }
137 +
138 +        if (info->getTime() >= currSample) {
139 +            dumpOut->writeDump(info->getTime());
140 +            currSample += sampleTime;
141 +        }
142 +
143 +        if (info->getTime() >= currStatus) {
144 +            statOut->writeStat(info->getTime());
145 +            calcPot = 0;
146 +            calcStress = 0;
147 +            currStatus += statusTime;
148 +        }
149 +
150 +        if (info->resetIntegrator) {
151 +            if (info->getTime() >= currReset) {
152 +                this->resetIntegrator();
153 +                currReset += resetTime;
154 +            }
155 +        }
156 +
157 + #ifdef PROFILE
158 +
159 +        endProfile(pro2);
160 +
161 + #endif //profile
162 +
163 + #ifdef IS_MPI
164 +
165 +        strcpy(checkPointMsg, "successfully took a time step.");
166 +        MPIcheckPoint();
167 +
168 + #endif // is_mpi
169 +
170 +    }
171 +
172 +    dumpOut->writeFinal(info->getTime());
173 +
174 +    delete dumpOut;
175 +    delete statOut;
176   }
177  
178 + void VelocityVerletIntegrator::integrateStep() { }
179 +
180 +
181   void VelocityVerletIntegrator::thermalize() {
182  
183 +    if (!info_->have_target_temp) {
184 +        sprintf(painCave.errMsg,
185 +                "You can't resample the velocities without a targetTemp!\n");
186 +        painCave.isFatal = 1;
187 +        painCave.severity = OOPSE_ERROR;
188 +        simError();
189 +        return;
190 +    }
191 +
192   }
193  
194 < void VelocityVerletIntegrator::velocitize() {
68 <  
194 > void VelocityVerletIntegrator::velocitize(double temperature) {
195      Vector3d aVel;
196      Vector3d aJ;
197      Mat3x3d I;
198 <    int l, m, n; // velocity randomizer loop counts
198 >    int l;
199 >    int m;
200 >    int n; // velocity randomizer loop counts
201      Vector3d vdrift;
202      double vbar;
203 <    const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
203 >    const double kb = 8.31451e - 7; // kb in amu, angstroms, fs, etc.
204      double av2;
205      double kebar;
78 –    double temperature;
206  
207 <    std::vector<Molecule*>::iterator i;
208 <    std::vector<StuntDouble*>::iterator j;
209 <    Molecule* mol;
210 <    StuntDouble* integrableObject;
207 >    std::vector<Molecule *>::iterator i;
208 >    std::vector<StuntDouble *>::iterator j;
209 >    Molecule * mol;
210 >    StuntDouble * integrableObject;
211      gaussianSPRNG gaussStream(info_->getSeed());
212  
213 <    if (!info->have_target_temp) {
87 <        sprintf( painCave.errMsg,
88 <             "You can't resample the velocities without a targetTemp!\n"
89 <             );
90 <        painCave.isFatal = 1;
91 <        painCave.severity = OOPSE_ERROR;
92 <        simError();
93 <        return;
94 <    }
213 >    kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf());
214  
215 <    temperature   = info_->target_temp;
215 >    for( mol = info_->beginMolecule(i); mol != NULL;
216 >        mol = info_->nextMolecule(i) ) {
217 >        for( integrableObject = mol->beginIntegrableObject(j);
218 >            integrableObject != NULL;
219 >            integrableObject = mol->nextIntegrableObject(j) ) {
220  
98 –    kebar = kb * temperature * info_->getNdfRaw() / ( 2.0 * info_->getNdf());
99 –
100 –    
101 –    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
102 –        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
103 –               integrableObject = mol->nextIntegrableObject(j)) {
104 –    
221              // uses equipartition theory to solve for vbar in angstrom/fs
222  
223              av2 = 2.0 * kebar / integrableObject->getMass();
224 <            vbar = sqrt( av2 );
224 >            vbar = sqrt(av2);
225  
226              // picks random velocities from a gaussian distribution
227              // centered on vbar
228  
229 <            for (int k=0; k<3; k++) {
229 >            for( int k = 0; k < 3; k++ ) {
230                  aVel[k] = vbar * gaussStream.getGaussian();
231              }
116 –            
117 –            integrableObject->setVel( aVel );
232  
233 <            if(integrableObject->isDirectional()){
233 >            integrableObject->setVel(aVel);
234  
235 +            if (integrableObject->isDirectional()) {
236                  I = integrableObject->getI();
237  
238                  if (integrableObject->isLinear()) {
239 +                    l = integrableObject->linearAxis();
240 +                    m = (l + 1) % 3;
241 +                    n = (l + 2) % 3;
242  
125 –                    l= integrableObject->linearAxis();
126 –                    m = (l+1)%3;
127 –                    n = (l+2)%3;
128 –
243                      aJ[l] = 0.0;
244 <                    vbar = sqrt( 2.0 * kebar * I(m, m) );
244 >                    vbar = sqrt(2.0 * kebar * I(m, m));
245                      aJ[m] = vbar * gaussStream.getGaussian();
246 <                    vbar = sqrt( 2.0 * kebar * I(n, n) );
246 >                    vbar = sqrt(2.0 * kebar * I(n, n));
247                      aJ[n] = vbar * gaussStream.getGaussian();
134 –
248                  } else {
249 <
250 <                    for (int k = 0 ; k < 3; k++) {
138 <                        vbar = sqrt( 2.0 * kebar * I(k, k) );
249 >                    for( int k = 0; k < 3; k++ ) {
250 >                        vbar = sqrt(2.0 * kebar * I(k, k));
251                          aJ[k] = vbar * gaussStream.getGaussian();
252 <                    }  
141 <
252 >                    }
253                  } // else isLinear
254  
255 <                integrableObject->setJ( aJ );
255 >                integrableObject->setJ(aJ);
256 >            }     //isDirectional
257 >        }
258 >    }             //end for (mol = beginMolecule(i); ...)
259  
146 –            }//isDirectional
147 –              
148 –       }
149 –    }//end for (mol = beginMolecule(i); ...)
150 –
260      // Get the Center of Mass drift velocity.
261      vdrift = info_->getComVel();
153 –  
154 –    //  Corrects for the center of mass drift.
155 –    // sums all the momentum and divides by total mass.
156 –    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
157 –        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
158 –               integrableObject = mol->nextIntegrableObject(j)) {
262  
263 <            aVel = integrableObject->getVel();
161 <            aVel -= vdrift;
162 <            integrableObject->setVel( aVel );            
163 <        }
164 <    }  
165 <    
166 < }
263 >    removeComDrift(vdrift);
264  
168 – void VelocityVerletIntegrator::calcForce(int needPotential, int needStress){
169 –
265   }
266  
267 < void VelocityVerletIntegrator::velocitize() {
267 > void VelocityVerletIntegrator::calcForce(bool needPotential,
268 >                                         bool needStress) { }
269  
270 < }
270 > void VelocityVerletIntegrator::removeComDrift(const Vector3d& vdrift) {
271  
272 < void removeComDrift(){
272 >    std::vector<Molecule *>::iterator i;
273 >    std::vector<StuntDouble *>::iterator j;
274 >    Molecule * mol;
275 >    StuntDouble * integrableObject;
276 >    
277 >    //  Corrects for the center of mass drift.
278 >    // sums all the momentum and divides by total mass.
279 >    for( mol = info_->beginMolecule(i); mol != NULL;
280 >        mol = info_->nextMolecule(i) ) {
281 >        for( integrableObject = mol->beginIntegrableObject(j);
282 >            integrableObject != NULL;
283 >            integrableObject = mol->nextIntegrableObject(j) ) {
284 >            integrableObject->setVel(integrableObject->getVel() - vdrift);
285 >        }
286 >    }
287  
288   }
289 +
290 + } //end namespace oopse

Diff Legend

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