ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/integrators/VelocityVerletIntegrator.cpp
(Generate patch)

Comparing branches/development/src/integrators/VelocityVerletIntegrator.cpp (file contents):
Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1878 by gezelter, Thu Jun 13 14:26:09 2013 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   /**
44   * @file VelocityVerletIntegrator.cpp
45   * @author tlin
46   * @date 11/09/2004
46 – * @time 16:16am
47   * @version 1.0
48   */
49  
# Line 53 | Line 53 | namespace OpenMD {
53   #include "utils/ProgressBar.hpp"
54  
55   namespace OpenMD {
56 <  VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info), rotAlgo(NULL) {
56 >  VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info) {
57      dt2 = 0.5 * dt;
58 –    rotAlgo = new DLM();
59 –    rattle = new Rattle(info);
58    }
59    
60    VelocityVerletIntegrator::~VelocityVerletIntegrator() {
63 –    delete rotAlgo;
64 –    delete rattle;
61    }
62    
63    void VelocityVerletIntegrator::initialize(){
64      
65 <    forceMan_->init();
65 >    forceMan_->initialize();
66      
67      // remove center of mass drift velocity (in case we passed in a
68      // configuration that was drifting)
69      velocitizer_->removeComDrift();
70 +
71 +    // find the initial fluctuating charges.
72 +    flucQ_->initialize();
73      
74      // initialize the forces before the first step
75      calcForce();
# Line 78 | Line 77 | namespace OpenMD {
77      // execute the constraint algorithm to make sure that the system is
78      // constrained at the very beginning  
79      if (info_->getNGlobalConstraints() > 0) {
80 <      rattle->constraintA();
80 >      rattle_->constraintA();
81        calcForce();
82 <      rattle->constraintB();      
82 >      rattle_->constraintB();      
83        //copy the current snapshot to previous snapshot
84        info_->getSnapshotManager()->advance();
85      }
# Line 98 | Line 97 | namespace OpenMD {
97      progressBar = new ProgressBar();
98  
99      //save statistics, before writeStat,  we must save statistics
101 –    thermo.saveStat();
100      saveConservedQuantity();
101 <    if (simParams->getUseRNEMD())
101 >    stats->collectStats();
102 >
103 >    if (simParams->getRNEMDParameters()->getUseRNEMD())
104        rnemd_->getStarted();
105  
106 <    statWriter->writeStat(currentSnapshot_->statData);
106 >    statWriter->writeStat();
107      
108 <    currSample = sampleTime + currentSnapshot_->getTime();
109 <    currStatus =  statusTime + currentSnapshot_->getTime();
110 <    currThermal = thermalTime + currentSnapshot_->getTime();
108 >    currSample = sampleTime + snap->getTime();
109 >    currStatus =  statusTime + snap->getTime();
110 >    currThermal = thermalTime + snap->getTime();
111      if (needReset) {
112 <      currReset = resetTime + currentSnapshot_->getTime();
112 >      currReset = resetTime + snap->getTime();
113      }
114 <    if (simParams->getUseRNEMD()){
115 <      currRNEMD = RNEMD_exchangeTime + currentSnapshot_->getTime();
114 >    if (simParams->getRNEMDParameters()->getUseRNEMD()){
115 >      currRNEMD = RNEMD_exchangeTime + snap->getTime();
116      }
117      needPotential = false;
118      needStress = false;      
# Line 121 | Line 121 | namespace OpenMD {
121  
122    void VelocityVerletIntegrator::doIntegrate() {
123    
124 –  
124      initialize();
125    
126 <    while (currentSnapshot_->getTime() < runTime) {
127 <    
128 <      preStep();
129 <    
130 <      integrateStep();
132 <    
133 <      postStep();
134 <      
135 <    }
136 <  
126 >    while (snap->getTime() < runTime) {
127 >      preStep();    
128 >      integrateStep();    
129 >      postStep();      
130 >    }  
131      finalize();
138 –  
132    }
133  
134  
135    void VelocityVerletIntegrator::preStep() {
136 <    RealType difference = currentSnapshot_->getTime() + dt - currStatus;
136 >    RealType difference = snap->getTime() + dt - currStatus;
137    
138      if (difference > 0 || fabs(difference) < OpenMD::epsilon) {
139        needPotential = true;
# Line 150 | Line 143 | namespace OpenMD {
143  
144    void VelocityVerletIntegrator::postStep() {
145  
153 –    //save snapshot
154 –    info_->getSnapshotManager()->advance();
155 –  
156 –    //increase time
157 –    currentSnapshot_->increaseTime(dt);        
146    
147      if (needVelocityScaling) {
148 <      if (currentSnapshot_->getTime() >= currThermal) {
148 >      if (snap->getTime() >= currThermal) {
149          velocitizer_->velocitize(targetScalingTemp);
150          currThermal += thermalTime;
151        }
152      }
153      if (useRNEMD) {
154 <      if (currentSnapshot_->getTime() >= currRNEMD) {
154 >      if (snap->getTime() >= currRNEMD) {
155          rnemd_->doRNEMD();
156          currRNEMD += RNEMD_exchangeTime;
157        }
158        rnemd_->collectData();
159      }
160      
161 <    if (currentSnapshot_->getTime() >= currSample) {
161 >    if (snap->getTime() >= currSample) {
162        dumpWriter->writeDumpAndEor();
163        
164        currSample += sampleTime;
165      }
166      
167 <    if (currentSnapshot_->getTime() >= currStatus) {
167 >    if (snap->getTime() >= currStatus) {
168        //save statistics, before writeStat,  we must save statistics
181 –      thermo.saveStat();
169        saveConservedQuantity();
170 +      stats->collectStats();
171  
172 <      if (simParams->getUseRNEMD()) {
173 <        rnemd_->getStatus();
172 >      if (simParams->getRNEMDParameters()->getUseRNEMD()) {
173 >        rnemd_->writeOutputFile();
174        }
175  
176 <      statWriter->writeStat(currentSnapshot_->statData);
176 >      statWriter->writeStat();
177  
178 <      progressBar->setStatus(currentSnapshot_->getTime(), runTime);
178 >      progressBar->setStatus(snap->getTime(), runTime);
179        progressBar->update();
180  
181        needPotential = false;
# Line 195 | Line 183 | namespace OpenMD {
183        currStatus += statusTime;
184      }
185      
186 <    if (needReset && currentSnapshot_->getTime() >= currReset) {    
186 >    if (needReset && snap->getTime() >= currReset) {    
187        resetIntegrator();
188        currReset += resetTime;
189      }        
190 +    //save snapshot
191 +    info_->getSnapshotManager()->advance();
192 +  
193 +    //increase time
194 +    snap->increaseTime(dt);        
195 +
196    }
197  
198  
199    void VelocityVerletIntegrator::finalize() {
200      dumpWriter->writeEor();
201 <  
201 >    if (simParams->getRNEMDParameters()->getUseRNEMD()) {
202 >      rnemd_->writeOutputFile();
203 >    }
204 >
205      delete dumpWriter;
206      delete statWriter;
207    
208      dumpWriter = NULL;
209      statWriter = NULL;
213 –  
210    }
211  
212    void VelocityVerletIntegrator::integrateStep() {
# Line 223 | Line 219 | namespace OpenMD {
219  
220    void VelocityVerletIntegrator::calcForce() {
221      forceMan_->calcForces();
222 +    flucQ_->applyConstraints();
223    }
224  
225    DumpWriter* VelocityVerletIntegrator::createDumpWriter() {
# Line 230 | Line 227 | namespace OpenMD {
227    }
228  
229    StatWriter* VelocityVerletIntegrator::createStatWriter() {
233 –
234 –    std::string statFileFormatString = simParams->getStatFileFormat();
235 –    StatsBitSet mask = parseStatFileFormat(statFileFormatString);
236 –  
237 –    // if we're doing a thermodynamic integration, we'll want the raw
238 –    // potential as well as the full potential:
239 –
240 –
241 –    if (simParams->getUseThermodynamicIntegration())
242 –      mask.set(Stats::VRAW);
243 –
244 –    // if we've got restraints turned on, we'll also want a report of the
245 –    // total harmonic restraints
246 –    if (simParams->getUseRestraints()){
247 –      mask.set(Stats::VHARM);
248 –    }
249 –
250 –    if (simParams->havePrintPressureTensor() &&
251 –        simParams->getPrintPressureTensor()){
252 –      mask.set(Stats::PRESSURE_TENSOR_XX);
253 –      mask.set(Stats::PRESSURE_TENSOR_XY);
254 –      mask.set(Stats::PRESSURE_TENSOR_XZ);
255 –      mask.set(Stats::PRESSURE_TENSOR_YX);
256 –      mask.set(Stats::PRESSURE_TENSOR_YY);
257 –      mask.set(Stats::PRESSURE_TENSOR_YZ);
258 –      mask.set(Stats::PRESSURE_TENSOR_ZX);
259 –      mask.set(Stats::PRESSURE_TENSOR_ZY);
260 –      mask.set(Stats::PRESSURE_TENSOR_ZZ);
261 –    }
230      
231 <    if (simParams->getAccumulateBoxDipole()) {
232 <      mask.set(Stats::BOX_DIPOLE_X);
265 <      mask.set(Stats::BOX_DIPOLE_Y);
266 <      mask.set(Stats::BOX_DIPOLE_Z);
267 <    }
268 <  
269 <    if (simParams->haveTaggedAtomPair() && simParams->havePrintTaggedPairDistance()) {
270 <      if (simParams->getPrintTaggedPairDistance()) {
271 <        mask.set(Stats::TAGGED_PAIR_DISTANCE);
272 <      }
273 <    }
274 <
275 <    if (simParams->getUseRNEMD()) {
276 <      mask.set(Stats::RNEMD_EXCHANGE_TOTAL);
277 <    }
231 >    stats = new Stats(info_);
232 >    statWriter = new StatWriter(info_->getStatFileName(), stats);
233      
234 <
280 <    return new StatWriter(info_->getStatFileName(), mask);
234 >    return statWriter;
235    }
236  
283 –
237   } //end namespace OpenMD

Diff Legend

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