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

Comparing branches/new_design/OOPSE-2.0/src/integrators/VelocityVerletIntegrator.cpp (file contents):
Revision 1725 by tim, Wed Nov 10 22:01:06 2004 UTC vs.
Revision 1820 by tim, Thu Dec 2 00:09:35 2004 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 >    rotAlgo = new DLM();
40 >    dt2 = 0.5 * dt;
41 > }
42  
43 < VelocityVerletIntegrator::~VelocityVerletIntegrator() { }
43 > VelocityVerletIntegrator::~VelocityVerletIntegrator() {
44 >    delete rotAlgo;
45 > }
46  
47 < void VelocityVerletIntegrator::integrate() {
42 <    double runTime = info->run_time;
47 > void VelocityVerletIntegrator::initialize(){
48  
49 <    double sampleTime = info->sampleTime;
50 <    double statusTime = info->statusTime;
51 <    double thermalTime = info->thermalTime;
52 <    double resetTime = info->resetTime;
49 >    currSample = 0.0;
50 >    currStatus = 0.0;
51 >    currThermal = 0.0;
52 >    needPotential = false;
53 >    needStress = false;      
54  
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 –
55      // remove center of mass drift velocity (in case we passed in a configuration
56      // that was drifting
57 <    tStats->removeCOMdrift();
57 >    velocitizer_->removeComDrift();
58  
59      // initialize the forces before the first step
60 +    calcForce(true, true);
61  
75 –    calcForce(1, 1);
76 –
62      //execute constraint algorithm to make sure at the very beginning the system is constrained  
63 <    if (nConstrained) {
64 <        constrainA();
65 <        calcForce(1, 1);
66 <        constrainB();
67 <    }
63 >    //if (nConstrained) {
64 >    //    constrainA();
65 >    //    calcForce(true, true);
66 >    //    constrainB();
67 >    //}
68  
69 <    if (info->setTemp) {
70 <        thermalize();
69 >    if (needVelocityScaling) {
70 >        velocitizer_->velocitize(targetScalingTemp);
71      }
72 +    
73 +    dumpWriter = createDumpWriter();
74 +    statWriter = createStatWriter();
75  
76 <    calcPot = 0;
77 <    calcStress = 0;
78 <    currSample = sampleTime + info->getTime();
79 <    currThermal = thermalTime + info->getTime();
80 <    currStatus = statusTime + info->getTime();
81 <    currReset = resetTime + info->getTime();
76 >    dumpWriter->writeDump();
77 >    statWriter->writeStat(currentSnapshot_->statData);
78 >  
79 > }
80 >        
81 > void VelocityVerletIntegrator::doIntegrate() {
82  
95 –    dumpOut->writeDump(info->getTime());
96 –    statOut->writeStat(info->getTime());
83  
84 < #ifdef IS_MPI
84 >    initialize();
85  
86 <    strcpy(checkPointMsg, "The integrator is ready to go.");
87 <    MPIcheckPoint();
86 >    while (currentSnapshot_->getTime() < runTime) {
87 >        
88 >        preStep();
89  
90 < #endif // is_mpi
90 >        integrateStep();
91 >        
92 >        postStep();
93  
94 <    while (info->getTime() < runTime && !stopIntegrator()) {
106 <        difference = info->getTime() + fullStep_ - currStatus;
94 >    }
95  
96 <        if (difference > 0 || fabs(difference) < 1e - 4) {
97 <            calcPot = 1;
98 <            calcStress = 1;
111 <        }
96 >    finalize();
97 >    
98 > }
99  
113 – #ifdef PROFILE
100  
101 <        startProfile(pro1);
101 > void VelocityVerletIntegrator::preStep() {
102 >        double difference = currentSnapshot_->getTime() + dt - currStatus;
103  
104 < #endif
104 >        if (difference > 0 || fabs(difference) < oopse::epsilon) {
105 >            needPotential = true;
106 >            needStress = true;  
107 >        }
108  
109 <        integrateStep(calcPot, calcStress);
109 > }
110  
111 < #ifdef PROFILE
111 > void VelocityVerletIntegrator::postStep() {
112  
113 <        endProfile(pro1);
113 >        //save snapshot
114 >        info_->getSnapshotManager()->advance();
115  
116 <        startProfile(pro2);
116 >        //increase time
117 >        currentSnapshot_->increaseTime(dt);        
118  
119 < #endif // profile
120 <
121 <        info->incrTime(fullStep_);
122 <
123 <        if (info->setTemp) {
124 <            if (info->getTime() >= currThermal) {
133 <                thermalize();
119 >        //save statistics
120 >        thermo.saveStat();
121 >        
122 >        if (needVelocityScaling) {
123 >            if (currentSnapshot_->getTime() >= currThermal) {
124 >                velocitizer_->velocitize(targetScalingTemp);
125                  currThermal += thermalTime;
126              }
127          }
128  
129 <        if (info->getTime() >= currSample) {
130 <            dumpOut->writeDump(info->getTime());
129 >        if (currentSnapshot_->getTime() >= currSample) {
130 >            dumpWriter->writeDump();
131              currSample += sampleTime;
132          }
133  
134 <        if (info->getTime() >= currStatus) {
135 <            statOut->writeStat(info->getTime());
136 <            calcPot = 0;
137 <            calcStress = 0;
134 >        if (currentSnapshot_->getTime() >= currStatus) {
135 >            statWriter->writeStat(currentSnapshot_->statData);
136 >            needPotential = false;
137 >            needStress = false;
138              currStatus += statusTime;
139          }
140  
141 <        if (info->resetIntegrator) {
142 <            if (info->getTime() >= currReset) {
152 <                this->resetIntegrator();
153 <                currReset += resetTime;
154 <            }
155 <        }
141 >        
142 > }
143  
157 – #ifdef PROFILE
144  
145 <        endProfile(pro2);
145 > void VelocityVerletIntegrator::finalize() {
146  
147 < #endif //profile
147 >    dumpWriter->writeDump();
148  
149 < #ifdef IS_MPI
149 >    delete dumpWriter;
150 >    delete statWriter;
151  
152 <        strcpy(checkPointMsg, "successfully took a time step.");
153 <        MPIcheckPoint();
152 >    dumpWriter = NULL;
153 >    statWriter = NULL;
154 >    
155 > }
156  
157 < #endif // is_mpi
157 > void VelocityVerletIntegrator::integrateStep() {
158  
159 <    }
160 <
161 <    dumpOut->writeFinal(info->getTime());
173 <
174 <    delete dumpOut;
175 <    delete statOut;
159 >    moveA();
160 >    calcForce(needPotential, needStress);
161 >    moveB();
162   }
163  
178 – void VelocityVerletIntegrator::integrateStep() { }
164  
165 <
166 < void VelocityVerletIntegrator::thermalize() {
167 <
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 <
165 > void VelocityVerletIntegrator::calcForce(bool needPotential,
166 >                                         bool needStress) {
167 >    forceMan_->calcForces(needPotential, needStress);
168   }
169  
170 < void VelocityVerletIntegrator::velocitize(double temperature) {
171 <    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 <
170 > DumpWriter* VelocityVerletIntegrator::createDumpWriter() {
171 >    return new DumpWriter(info_, info_->getDumpFileName());
172   }
173  
174 < void VelocityVerletIntegrator::calcForce(bool needPotential,
175 <                                         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 <
174 > StatWriter* VelocityVerletIntegrator::createStatWriter() {
175 >    return new StatWriter(info_->getStatFileName());
176   }
177  
178 +
179   } //end namespace oopse

Diff Legend

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