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 1756 by tim, Thu Nov 18 23:26:27 2004 UTC vs.
Revision 1847 by tim, Sat Dec 4 05:24:07 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) : Integrator(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 <    int calcPot;
50 <    int calcStress;
49 >    currSample = 0.0;
50 >    currStatus = 0.0;
51 >    currThermal = 0.0;
52 >    needPotential = false;
53 >    needStress = false;      
54  
47
48    fullStep_ = info_->dt;
49    halfStep_ = 0.5 * fullStep_;
50
51    readyCheck();
52
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  
59    calcForce(1, 1);
60
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;
73 <    calcStress = 0;
76 >    dumpWriter->writeDump();
77  
78 <    tStats = new Thermo(info_);
79 <    statOut = new StatWriter(info_);
80 <    dumpOut = new DumpWriter(info_);
78 <    
79 <    dumpOut->writeDump(info_->getTime());
80 <    statOut->writeStat(info_->getTime());
78 >    //save statistics
79 >    thermo.saveStat();
80 >    saveConservedQuantity();
81  
82 < #ifdef IS_MPI
82 >    //before writeStat,  we must save statistics
83 >    statWriter->writeStat(currentSnapshot_->statData);
84 >  
85 > }
86 >        
87 > void VelocityVerletIntegrator::doIntegrate() {
88  
84    strcpy(checkPointMsg, "The integrator is ready to go.");
85    MPIcheckPoint();
89  
90 < #endif // is_mpi
90 >    initialize();
91  
92 <    while (info_->getTime() < runTime) {
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;
96 >        integrateStep();
97 >        
98 >        postStep();
99 >
100 >    }
101 >
102 >    finalize();
103 >    
104 > }
105 >
106 >
107 > void VelocityVerletIntegrator::preStep() {
108 >        double difference = currentSnapshot_->getTime() + dt - currStatus;
109 >
110 >        if (difference > 0 || fabs(difference) < oopse::epsilon) {
111 >            needPotential = true;
112 >            needStress = true;  
113          }
114  
115 <        //notify before integrateStep
116 <        notify()
117 <        integrateStep(calcPot, calcStress);
115 > }
116 >
117 > void VelocityVerletIntegrator::postStep() {
118 >
119 >        //save snapshot
120 >        info_->getSnapshotManager()->advance();
121 >
122 >        //increase time
123 >        currentSnapshot_->increaseTime(dt);        
124 >
125 >        //save statistics
126 >        thermo.saveStat();
127 >        saveConservedQuantity();
128          
129 <        info_->incrTime(fullStep_);
130 <        
131 <        //notify after integratreStep        
104 <        notify();
105 <        
106 <        if (info_->setTemp) {
107 <            if (info_->getTime() >= currThermal) {
108 <                thermalize();
129 >        if (needVelocityScaling) {
130 >            if (currentSnapshot_->getTime() >= currThermal) {
131 >                velocitizer_->velocitize(targetScalingTemp);
132                  currThermal += thermalTime;
133              }
134          }
135  
136 <        if (info_->getTime() >= currSample) {
137 <            dumpOut->writeDump(info_->getTime());
136 >        if (currentSnapshot_->getTime() >= currSample) {
137 >            dumpWriter->writeDump();
138              currSample += sampleTime;
139          }
140  
141 <        if (info_->getTime() >= currStatus) {
142 <            statOut->writeStat(info_->getTime());
143 <            calcPot = 0;
144 <            calcStress = 0;
141 >        if (currentSnapshot_->getTime() >= currStatus) {
142 >            statWriter->writeStat(currentSnapshot_->statData);
143 >            needPotential = false;
144 >            needStress = false;
145              currStatus += statusTime;
146          }
147  
148 <        if (info_->resetIntegrator) {
149 <            if (info_->getTime() >= currReset) {
127 <                this->resetIntegrator();
128 <                currReset += resetTime;
129 <            }
130 <        }
148 >        
149 > }
150  
132 #ifdef IS_MPI
151  
152 <        strcpy(checkPointMsg, "successfully took a time step.");
135 <        MPIcheckPoint();
152 > void VelocityVerletIntegrator::finalize() {
153  
154 < #endif // is_mpi
154 >    dumpWriter->writeDump();
155  
156 <    }
156 >    delete dumpWriter;
157 >    delete statWriter;
158  
159 <    dumpOut->writeFinal(info_->getTime());
160 <
161 <    delete dumpOut;
144 <    delete statOut;
159 >    dumpWriter = NULL;
160 >    statWriter = NULL;
161 >    
162   }
163  
164 < void VelocityVerletIntegrator::integrateStep() { }
164 > void VelocityVerletIntegrator::integrateStep() {
165  
166 +    moveA();
167 +    calcForce(needPotential, needStress);
168 +    moveB();
169 + }
170  
150 void VelocityVerletIntegrator::thermalize() {
171  
172 <    if (!info__->have_target_temp) {
173 <        sprintf(painCave.errMsg,
174 <                "You can't resample the velocities without a targetTemp!\n");
175 <        painCave.isFatal = 1;
156 <        painCave.severity = OOPSE_ERROR;
157 <        simError();
158 <        return;
159 <    }
172 > void VelocityVerletIntegrator::calcForce(bool needPotential,
173 >                                         bool needStress) {
174 >    forceMan_->calcForces(needPotential, needStress);
175 > }
176  
177 + DumpWriter* VelocityVerletIntegrator::createDumpWriter() {
178 +    return new DumpWriter(info_, info_->getDumpFileName());
179   }
180  
181 < void VelocityVerletIntegrator::calcForce(bool needPotential,
182 <                                         bool needStress) { }
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