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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines