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 1725 by tim, Wed Nov 10 22:01:06 2004 UTC vs.
Revision 1901 by tim, Tue Jan 4 22:18:36 2005 UTC

# Line 32 | Line 32
32   */
33  
34   #include "integrators/VelocityVerletIntegrator.hpp"
35 + #include "integrators/DLM.hpp"
36  
37   namespace oopse {
38 < VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) { }
38 > VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info), rotAlgo(NULL) {
39 >    dt2 = 0.5 * dt;
40 >    rotAlgo = new DLM();
41 >    rattle = new Rattle(info);
42 > }
43  
44 < VelocityVerletIntegrator::~VelocityVerletIntegrator() { }
44 > VelocityVerletIntegrator::~VelocityVerletIntegrator() {
45 >    delete rotAlgo;
46 >    delete rattle;
47 > }
48  
49 < void VelocityVerletIntegrator::integrate() {
42 <    double runTime = info->run_time;
49 > void VelocityVerletIntegrator::initialize(){
50  
44 –    double sampleTime = info->sampleTime;
45 –    double statusTime = info->statusTime;
46 –    double thermalTime = info->thermalTime;
47 –    double resetTime = info->resetTime;
48 –
49 –    double difference;
50 –    double currSample;
51 –    double currThermal;
52 –    double currStatus;
53 –    double currReset;
54 –
55 –    int calcPot,
56 –    calcStress;
57 –
58 –    tStats = new Thermo(info);
59 –    statOut = new StatWriter(info);
60 –    dumpOut = new DumpWriter(info);
61 –
62 –    atoms = info->atoms;
63 –
64 –    fullStep_ = info->dt;
65 –    halfStep_ = 0.5 * fullStep_;
66 –
67 –    readyCheck();
68 –
51      // remove center of mass drift velocity (in case we passed in a configuration
52      // that was drifting
53 <    tStats->removeCOMdrift();
53 >    velocitizer_->removeComDrift();
54  
55      // initialize the forces before the first step
56 +    calcForce(true, true);
57  
75 –    calcForce(1, 1);
76 –
58      //execute constraint algorithm to make sure at the very beginning the system is constrained  
59 <    if (nConstrained) {
60 <        constrainA();
61 <        calcForce(1, 1);
62 <        constrainB();
63 <    }
59 >    //if (nConstrained) {
60 >    //    constrainA();
61 >    //    calcForce(true, true);
62 >    //    constrainB();
63 >    //}
64  
65 <    if (info->setTemp) {
66 <        thermalize();
65 >    if (needVelocityScaling) {
66 >        velocitizer_->velocitize(targetScalingTemp);
67      }
68 +    
69 +    dumpWriter = createDumpWriter();
70 +    statWriter = createStatWriter();
71  
72 <    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();
72 >    dumpWriter->writeDump();
73  
74 <    dumpOut->writeDump(info->getTime());
75 <    statOut->writeStat(info->getTime());
74 >    //save statistics, before writeStat,  we must save statistics
75 >    thermo.saveStat();
76 >    saveConservedQuantity();
77 >    statWriter->writeStat(currentSnapshot_->statData);
78  
79 < #ifdef IS_MPI
79 >    currSample = sampleTime + currentSnapshot_->getTime();
80 >    currStatus =  statusTime + currentSnapshot_->getTime();;
81 >    currThermal = thermalTime +  + currentSnapshot_->getTime();
82 >    needPotential = false;
83 >    needStress = false;      
84 >  
85 > }
86 >        
87 > void VelocityVerletIntegrator::doIntegrate() {
88  
100 –    strcpy(checkPointMsg, "The integrator is ready to go.");
101 –    MPIcheckPoint();
89  
90 < #endif // is_mpi
90 >    initialize();
91  
92 <    while (info->getTime() < runTime && !stopIntegrator()) {
93 <        difference = info->getTime() + fullStep_ - currStatus;
92 >    while (currentSnapshot_->getTime() < runTime) {
93 >        
94 >        preStep();
95  
96 <        if (difference > 0 || fabs(difference) < 1e - 4) {
97 <            calcPot = 1;
98 <            calcStress = 1;
111 <        }
96 >        integrateStep();
97 >        
98 >        postStep();
99  
100 < #ifdef PROFILE
100 >    }
101  
102 <        startProfile(pro1);
102 >    finalize();
103 >    
104 > }
105  
117 – #endif
106  
107 <        integrateStep(calcPot, calcStress);
107 > void VelocityVerletIntegrator::preStep() {
108 >        double difference = currentSnapshot_->getTime() + dt - currStatus;
109  
110 < #ifdef PROFILE
110 >        if (difference > 0 || fabs(difference) < oopse::epsilon) {
111 >            needPotential = true;
112 >            needStress = true;  
113 >        }
114  
115 <        endProfile(pro1);
115 > }
116  
117 <        startProfile(pro2);
117 > void VelocityVerletIntegrator::postStep() {
118  
119 < #endif // profile
119 >        //save snapshot
120 >        info_->getSnapshotManager()->advance();
121  
122 <        info->incrTime(fullStep_);
123 <
124 <        if (info->setTemp) {
125 <            if (info->getTime() >= currThermal) {
126 <                thermalize();
122 >        //increase time
123 >        currentSnapshot_->increaseTime(dt);        
124 >        
125 >        if (needVelocityScaling) {
126 >            if (currentSnapshot_->getTime() >= currThermal) {
127 >                velocitizer_->velocitize(targetScalingTemp);
128                  currThermal += thermalTime;
129              }
130          }
131  
132 <        if (info->getTime() >= currSample) {
133 <            dumpOut->writeDump(info->getTime());
132 >        if (currentSnapshot_->getTime() >= currSample) {
133 >            dumpWriter->writeDump();
134              currSample += sampleTime;
135          }
136  
137 <        if (info->getTime() >= currStatus) {
138 <            statOut->writeStat(info->getTime());
139 <            calcPot = 0;
140 <            calcStress = 0;
137 >        if (currentSnapshot_->getTime() >= currStatus) {
138 >            //save statistics, before writeStat,  we must save statistics
139 >            thermo.saveStat();
140 >            saveConservedQuantity();
141 >            statWriter->writeStat(currentSnapshot_->statData);
142 >
143 >            needPotential = false;
144 >            needStress = false;
145              currStatus += statusTime;
146          }
147  
148 <        if (info->resetIntegrator) {
149 <            if (info->getTime() >= currReset) {
152 <                this->resetIntegrator();
153 <                currReset += resetTime;
154 <            }
155 <        }
148 >        
149 > }
150  
157 – #ifdef PROFILE
151  
152 <        endProfile(pro2);
152 > void VelocityVerletIntegrator::finalize() {
153  
154 < #endif //profile
154 >    dumpWriter->writeDump();
155  
156 < #ifdef IS_MPI
156 >    delete dumpWriter;
157 >    delete statWriter;
158  
159 <        strcpy(checkPointMsg, "successfully took a time step.");
160 <        MPIcheckPoint();
159 >    dumpWriter = NULL;
160 >    statWriter = NULL;
161 >    
162 > }
163  
164 < #endif // is_mpi
164 > void VelocityVerletIntegrator::integrateStep() {
165  
166 <    }
167 <
168 <    dumpOut->writeFinal(info->getTime());
173 <
174 <    delete dumpOut;
175 <    delete statOut;
166 >    moveA();
167 >    calcForce(needPotential, needStress);
168 >    moveB();
169   }
170  
178 – void VelocityVerletIntegrator::integrateStep() { }
171  
172 <
173 < void VelocityVerletIntegrator::thermalize() {
174 <
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 <
172 > void VelocityVerletIntegrator::calcForce(bool needPotential,
173 >                                         bool needStress) {
174 >    forceMan_->calcForces(needPotential, needStress);
175   }
176  
177 < void VelocityVerletIntegrator::velocitize(double temperature) {
178 <    Vector3d aVel;
196 <    Vector3d aJ;
197 <    Mat3x3d I;
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.
204 <    double av2;
205 <    double kebar;
206 <
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 <    kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf());
214 <
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 <
221 <            // uses equipartition theory to solve for vbar in angstrom/fs
222 <
223 <            av2 = 2.0 * kebar / integrableObject->getMass();
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++ ) {
230 <                aVel[k] = vbar * gaussStream.getGaussian();
231 <            }
232 <
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 <
243 <                    aJ[l] = 0.0;
244 <                    vbar = sqrt(2.0 * kebar * I(m, m));
245 <                    aJ[m] = vbar * gaussStream.getGaussian();
246 <                    vbar = sqrt(2.0 * kebar * I(n, n));
247 <                    aJ[n] = vbar * gaussStream.getGaussian();
248 <                } else {
249 <                    for( int k = 0; k < 3; k++ ) {
250 <                        vbar = sqrt(2.0 * kebar * I(k, k));
251 <                        aJ[k] = vbar * gaussStream.getGaussian();
252 <                    }
253 <                } // else isLinear
254 <
255 <                integrableObject->setJ(aJ);
256 <            }     //isDirectional
257 <        }
258 <    }             //end for (mol = beginMolecule(i); ...)
259 <
260 <    // Get the Center of Mass drift velocity.
261 <    vdrift = info_->getComVel();
262 <
263 <    removeComDrift(vdrift);
264 <
177 > DumpWriter* VelocityVerletIntegrator::createDumpWriter() {
178 >    return new DumpWriter(info_, info_->getDumpFileName());
179   }
180  
181 < void VelocityVerletIntegrator::calcForce(bool needPotential,
182 <                                         bool needStress) { }
269 <
270 < void VelocityVerletIntegrator::removeComDrift(const Vector3d& vdrift) {
271 <
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 <
181 > StatWriter* VelocityVerletIntegrator::createStatWriter() {
182 >    return new StatWriter(info_->getStatFileName());
183   }
184  
185 +
186   } //end namespace oopse

Diff Legend

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