ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/integrators/VelocityVerletIntegrator.cpp
Revision: 3501
Committed: Fri May 8 19:47:05 2009 UTC (16 years, 2 months ago) by skuang
File size: 8292 byte(s)
Log Message:
Bug fixes and sanity checks for RNEMD (nBins should be even), check to make
sure we aren't selecting inappropriate numbers of integrable Objects

File Contents

# User Rev Content
1 gezelter 2204 /*
2 gezelter 1930 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     *
4     * The University of Notre Dame grants you ("Licensee") a
5     * non-exclusive, royalty free, license to use, modify and
6     * redistribute this software in source and binary code form, provided
7     * that the following conditions are met:
8     *
9     * 1. Acknowledgement of the program authors must be made in any
10     * publication of scientific results based in part on use of the
11     * program. An acceptable form of acknowledgement is citation of
12     * the article in which the program was described (Matthew
13     * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14     * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15     * Parallel Simulation Engine for Molecular Dynamics,"
16     * J. Comput. Chem. 26, pp. 252-271 (2005))
17     *
18     * 2. Redistributions of source code must retain the above copyright
19     * notice, this list of conditions and the following disclaimer.
20     *
21     * 3. Redistributions in binary form must reproduce the above copyright
22     * notice, this list of conditions and the following disclaimer in the
23     * documentation and/or other materials provided with the
24     * distribution.
25     *
26     * This software is provided "AS IS," without a warranty of any
27     * kind. All express or implied conditions, representations and
28     * warranties, including any implied warranty of merchantability,
29     * fitness for a particular purpose or non-infringement, are hereby
30     * excluded. The University of Notre Dame and its licensors shall not
31     * be liable for any damages suffered by licensee as a result of
32     * using, modifying or distributing the software or its
33     * derivatives. In no event will the University of Notre Dame or its
34     * licensors be liable for any lost revenue, profit or data, or for
35     * direct, indirect, special, consequential, incidental or punitive
36     * damages, however caused and regardless of the theory of liability,
37     * arising out of the use of or inability to use software, even if the
38     * University of Notre Dame has been advised of the possibility of
39     * such damages.
40     */
41    
42     /**
43     * @file VelocityVerletIntegrator.cpp
44     * @author tlin
45     * @date 11/09/2004
46     * @time 16:16am
47     * @version 1.0
48     */
49    
50     #include "integrators/VelocityVerletIntegrator.hpp"
51     #include "integrators/DLM.hpp"
52 tim 1960 #include "utils/StringUtils.hpp"
53 gezelter 1930
54     namespace oopse {
55 chrisfen 2101 VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info), rotAlgo(NULL) {
56 gezelter 1930 dt2 = 0.5 * dt;
57     rotAlgo = new DLM();
58     rattle = new Rattle(info);
59 chrisfen 2101 }
60    
61     VelocityVerletIntegrator::~VelocityVerletIntegrator() {
62 gezelter 1930 delete rotAlgo;
63     delete rattle;
64 chrisfen 2101 }
65    
66     void VelocityVerletIntegrator::initialize(){
67    
68 gezelter 1930 forceMan_->init();
69 chrisfen 2101
70 chrisfen 2879 // remove center of mass drift velocity (in case we passed in a
71     // configuration that was drifting)
72 gezelter 1930 velocitizer_->removeComDrift();
73 chrisfen 2101
74 gezelter 1930 // initialize the forces before the first step
75     calcForce(true, true);
76 chrisfen 2101
77 chrisfen 2879 // execute the constraint algorithm to make sure that the system is
78     // constrained at the very beginning
79 gezelter 1930 if (info_->getNGlobalConstraints() > 0) {
80 chrisfen 2101 rattle->constraintA();
81     calcForce(true, true);
82 chrisfen 2879 rattle->constraintB();
83     //copy the current snapshot to previous snapshot
84     info_->getSnapshotManager()->advance();
85 gezelter 1930 }
86 chrisfen 2101
87 gezelter 1930 if (needVelocityScaling) {
88 chrisfen 2101 velocitizer_->velocitize(targetScalingTemp);
89 gezelter 1930 }
90    
91     dumpWriter = createDumpWriter();
92 chrisfen 2101
93 gezelter 1930 statWriter = createStatWriter();
94 chrisfen 2879
95     dumpWriter->writeDumpAndEor();
96    
97 chrisfen 2101 if (simParams->getUseSolidThermInt()) {
98     restWriter = createRestWriter();
99 chrisfen 2879 restWriter->writeZAngFile();
100 chrisfen 2101 }
101    
102 gezelter 1930 //save statistics, before writeStat, we must save statistics
103     thermo.saveStat();
104     saveConservedQuantity();
105 skuang 3501 if (simParams->getUseRNEMD())
106     rnemd_->getStatus();
107    
108 gezelter 1930 statWriter->writeStat(currentSnapshot_->statData);
109 chrisfen 2101
110 gezelter 1930 currSample = sampleTime + currentSnapshot_->getTime();
111 gezelter 3488 currStatus = statusTime + currentSnapshot_->getTime();
112 tim 2243 currThermal = thermalTime + currentSnapshot_->getTime();
113     if (needReset) {
114     currReset = resetTime + currentSnapshot_->getTime();
115     }
116 gezelter 3488 if (simParams->getUseRNEMD()){
117     currRNEMD = RNEMD_swapTime + currentSnapshot_->getTime();
118     }
119 gezelter 1930 needPotential = false;
120     needStress = false;
121 chrisfen 2101
122 gezelter 2204 }
123 chrisfen 2101
124 gezelter 2204 void VelocityVerletIntegrator::doIntegrate() {
125 chrisfen 2101
126    
127 gezelter 2204 initialize();
128 chrisfen 2101
129 gezelter 2204 while (currentSnapshot_->getTime() < runTime) {
130 gezelter 1930
131 gezelter 2204 preStep();
132 chrisfen 2101
133 gezelter 2204 integrateStep();
134 chrisfen 2101
135 gezelter 2204 postStep();
136 chrisfen 2101
137 gezelter 2204 }
138 chrisfen 2101
139 gezelter 2204 finalize();
140 chrisfen 2101
141 gezelter 2204 }
142 gezelter 1930
143    
144 gezelter 2204 void VelocityVerletIntegrator::preStep() {
145 tim 2759 RealType difference = currentSnapshot_->getTime() + dt - currStatus;
146 chrisfen 2101
147 gezelter 2204 if (difference > 0 || fabs(difference) < oopse::epsilon) {
148     needPotential = true;
149     needStress = true;
150     }
151 chrisfen 2101 }
152 gezelter 1930
153 gezelter 2204 void VelocityVerletIntegrator::postStep() {
154 tim 2448
155 gezelter 2204 //save snapshot
156     info_->getSnapshotManager()->advance();
157 chrisfen 2101
158 gezelter 2204 //increase time
159     currentSnapshot_->increaseTime(dt);
160 chrisfen 2879
161 gezelter 2204 if (needVelocityScaling) {
162     if (currentSnapshot_->getTime() >= currThermal) {
163     velocitizer_->velocitize(targetScalingTemp);
164     currThermal += thermalTime;
165     }
166 chrisfen 2101 }
167 gezelter 3488 if (useRNEMD) {
168     if (currentSnapshot_->getTime() >= currRNEMD) {
169     rnemd_->doSwap();
170     currRNEMD += RNEMD_swapTime;
171     }
172     }
173    
174 gezelter 2204 if (currentSnapshot_->getTime() >= currSample) {
175     dumpWriter->writeDumpAndEor();
176 gezelter 3488
177 gezelter 2204 if (simParams->getUseSolidThermInt())
178 chrisfen 2879 restWriter->writeZAngFile();
179 gezelter 3488
180 gezelter 2204 currSample += sampleTime;
181     }
182 gezelter 3488
183 gezelter 2204 if (currentSnapshot_->getTime() >= currStatus) {
184     //save statistics, before writeStat, we must save statistics
185     thermo.saveStat();
186     saveConservedQuantity();
187 skuang 3501
188 skuang 3498 if (simParams->getUseRNEMD())
189     rnemd_->getStatus();
190 skuang 3501
191     statWriter->writeStat(currentSnapshot_->statData);
192 gezelter 3488
193 gezelter 2204 needPotential = false;
194     needStress = false;
195     currStatus += statusTime;
196     }
197 gezelter 3488
198     if (needReset && currentSnapshot_->getTime() >= currReset) {
199     resetIntegrator();
200     currReset += resetTime;
201     }
202 gezelter 2204 }
203 gezelter 1930
204    
205 gezelter 2204 void VelocityVerletIntegrator::finalize() {
206     dumpWriter->writeEor();
207 chrisfen 2101
208 gezelter 2204 if (simParams->getUseSolidThermInt()) {
209 chrisfen 2991 restWriter->writeZAngFile();
210 gezelter 2204 delete restWriter;
211     restWriter = NULL;
212     }
213 chrisfen 2101
214 gezelter 2204 delete dumpWriter;
215     delete statWriter;
216 chrisfen 2101
217 gezelter 2204 dumpWriter = NULL;
218     statWriter = NULL;
219 chrisfen 2101
220 gezelter 2204 }
221 gezelter 1930
222 gezelter 2204 void VelocityVerletIntegrator::integrateStep() {
223 chrisfen 2101
224 gezelter 2204 moveA();
225     calcForce(needPotential, needStress);
226     moveB();
227     }
228 gezelter 1930
229    
230 gezelter 2204 void VelocityVerletIntegrator::calcForce(bool needPotential,
231     bool needStress) {
232     forceMan_->calcForces(needPotential, needStress);
233     }
234 gezelter 1930
235 gezelter 2204 DumpWriter* VelocityVerletIntegrator::createDumpWriter() {
236     return new DumpWriter(info_);
237     }
238 gezelter 1930
239 gezelter 2204 StatWriter* VelocityVerletIntegrator::createStatWriter() {
240 tim 2380
241     std::string statFileFormatString = simParams->getStatFileFormat();
242     StatsBitSet mask = parseStatFileFormat(statFileFormatString);
243    
244 gezelter 2204 // if solidThermInt is true, add extra information to the statfile
245     if (simParams->getUseSolidThermInt()){
246     mask.set(Stats::VRAW);
247     mask.set(Stats::VHARM);
248     }
249 tim 2238
250 chrisfen 2879 if (simParams->havePrintPressureTensor() &&
251     simParams->getPrintPressureTensor()){
252 gezelter 3448 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 tim 2238 }
262    
263 chrisfen 2917 if (simParams->getAccumulateBoxDipole()) {
264     mask.set(Stats::BOX_DIPOLE_X);
265     mask.set(Stats::BOX_DIPOLE_Y);
266     mask.set(Stats::BOX_DIPOLE_Z);
267     }
268 gezelter 3488
269 gezelter 3448 if (simParams->haveTaggedAtomPair() && simParams->havePrintTaggedPairDistance()) {
270     if (simParams->getPrintTaggedPairDistance()) {
271     mask.set(Stats::TAGGED_PAIR_DISTANCE);
272     }
273     }
274 gezelter 3488
275     if (simParams->getUseRNEMD()) {
276     mask.set(Stats::RNEMD_SWAP_TOTAL);
277     }
278 gezelter 3448
279 chrisfen 2917
280 tim 2380 return new StatWriter(info_->getStatFileName(), mask);
281 chrisfen 2101 }
282 gezelter 1930
283 gezelter 2204 RestWriter* VelocityVerletIntegrator::createRestWriter(){
284     return new RestWriter(info_);
285     }
286 gezelter 1930
287 chrisfen 2101
288 gezelter 1930 } //end namespace oopse

Properties

Name Value
svn:executable *