1 |
tim |
1722 |
/* |
2 |
|
|
* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project |
3 |
|
|
* |
4 |
|
|
* Contact: oopse@oopse.org |
5 |
|
|
* |
6 |
|
|
* This program is free software; you can redistribute it and/or |
7 |
|
|
* modify it under the terms of the GNU Lesser General Public License |
8 |
|
|
* as published by the Free Software Foundation; either version 2.1 |
9 |
|
|
* of the License, or (at your option) any later version. |
10 |
|
|
* All we ask is that proper credit is given for our work, which includes |
11 |
|
|
* - but is not limited to - adding the above copyright notice to the beginning |
12 |
|
|
* of your source code files, and to any copyright notice that you may distribute |
13 |
|
|
* with programs based on this work. |
14 |
|
|
* |
15 |
|
|
* This program is distributed in the hope that it will be useful, |
16 |
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of |
17 |
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
18 |
|
|
* GNU Lesser General Public License for more details. |
19 |
|
|
* |
20 |
|
|
* You should have received a copy of the GNU Lesser General Public License |
21 |
|
|
* along with this program; if not, write to the Free Software |
22 |
|
|
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
23 |
|
|
* |
24 |
|
|
*/ |
25 |
|
|
|
26 |
tim |
1725 |
/** |
27 |
|
|
* @file VelocityVerletIntegrator.cpp |
28 |
|
|
* @author tlin |
29 |
|
|
* @date 11/09/2004 |
30 |
|
|
* @time 16:16am |
31 |
|
|
* @version 1.0 |
32 |
|
|
*/ |
33 |
|
|
|
34 |
tim |
1722 |
#include "integrators/VelocityVerletIntegrator.hpp" |
35 |
|
|
|
36 |
|
|
namespace oopse { |
37 |
tim |
1756 |
VelocityVerletIntegrator::VelocityVerletIntegrator(SimInfo *info) : Integrator(info) { } |
38 |
tim |
1722 |
|
39 |
tim |
1725 |
VelocityVerletIntegrator::~VelocityVerletIntegrator() { } |
40 |
tim |
1722 |
|
41 |
tim |
1725 |
void VelocityVerletIntegrator::integrate() { |
42 |
tim |
1756 |
double runTime = info_->run_time; |
43 |
tim |
1722 |
|
44 |
tim |
1756 |
int calcPot; |
45 |
|
|
int calcStress; |
46 |
tim |
1722 |
|
47 |
|
|
|
48 |
tim |
1756 |
fullStep_ = info_->dt; |
49 |
tim |
1725 |
halfStep_ = 0.5 * fullStep_; |
50 |
tim |
1722 |
|
51 |
tim |
1725 |
readyCheck(); |
52 |
tim |
1722 |
|
53 |
tim |
1725 |
// remove center of mass drift velocity (in case we passed in a configuration |
54 |
|
|
// that was drifting |
55 |
|
|
tStats->removeCOMdrift(); |
56 |
tim |
1722 |
|
57 |
tim |
1725 |
// initialize the forces before the first step |
58 |
tim |
1722 |
|
59 |
tim |
1725 |
calcForce(1, 1); |
60 |
tim |
1722 |
|
61 |
tim |
1725 |
//execute constraint algorithm to make sure at the very beginning the system is constrained |
62 |
|
|
if (nConstrained) { |
63 |
|
|
constrainA(); |
64 |
|
|
calcForce(1, 1); |
65 |
|
|
constrainB(); |
66 |
|
|
} |
67 |
|
|
|
68 |
tim |
1756 |
if (info_->setTemp) { |
69 |
tim |
1725 |
thermalize(); |
70 |
|
|
} |
71 |
|
|
|
72 |
|
|
calcPot = 0; |
73 |
|
|
calcStress = 0; |
74 |
|
|
|
75 |
tim |
1756 |
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()); |
81 |
tim |
1725 |
|
82 |
|
|
#ifdef IS_MPI |
83 |
|
|
|
84 |
|
|
strcpy(checkPointMsg, "The integrator is ready to go."); |
85 |
|
|
MPIcheckPoint(); |
86 |
|
|
|
87 |
|
|
#endif // is_mpi |
88 |
|
|
|
89 |
tim |
1756 |
while (info_->getTime() < runTime) { |
90 |
|
|
difference = info_->getTime() + fullStep_ - currStatus; |
91 |
tim |
1725 |
|
92 |
|
|
if (difference > 0 || fabs(difference) < 1e - 4) { |
93 |
|
|
calcPot = 1; |
94 |
|
|
calcStress = 1; |
95 |
|
|
} |
96 |
|
|
|
97 |
tim |
1756 |
//notify before integrateStep |
98 |
|
|
notify() |
99 |
tim |
1725 |
integrateStep(calcPot, calcStress); |
100 |
tim |
1756 |
|
101 |
|
|
info_->incrTime(fullStep_); |
102 |
|
|
|
103 |
|
|
//notify after integratreStep |
104 |
|
|
notify(); |
105 |
|
|
|
106 |
|
|
if (info_->setTemp) { |
107 |
|
|
if (info_->getTime() >= currThermal) { |
108 |
tim |
1725 |
thermalize(); |
109 |
|
|
currThermal += thermalTime; |
110 |
|
|
} |
111 |
|
|
} |
112 |
|
|
|
113 |
tim |
1756 |
if (info_->getTime() >= currSample) { |
114 |
|
|
dumpOut->writeDump(info_->getTime()); |
115 |
tim |
1725 |
currSample += sampleTime; |
116 |
|
|
} |
117 |
|
|
|
118 |
tim |
1756 |
if (info_->getTime() >= currStatus) { |
119 |
|
|
statOut->writeStat(info_->getTime()); |
120 |
tim |
1725 |
calcPot = 0; |
121 |
|
|
calcStress = 0; |
122 |
|
|
currStatus += statusTime; |
123 |
|
|
} |
124 |
|
|
|
125 |
tim |
1756 |
if (info_->resetIntegrator) { |
126 |
|
|
if (info_->getTime() >= currReset) { |
127 |
tim |
1725 |
this->resetIntegrator(); |
128 |
|
|
currReset += resetTime; |
129 |
|
|
} |
130 |
|
|
} |
131 |
|
|
|
132 |
|
|
#ifdef IS_MPI |
133 |
|
|
|
134 |
|
|
strcpy(checkPointMsg, "successfully took a time step."); |
135 |
|
|
MPIcheckPoint(); |
136 |
|
|
|
137 |
|
|
#endif // is_mpi |
138 |
|
|
|
139 |
|
|
} |
140 |
|
|
|
141 |
tim |
1756 |
dumpOut->writeFinal(info_->getTime()); |
142 |
tim |
1725 |
|
143 |
|
|
delete dumpOut; |
144 |
|
|
delete statOut; |
145 |
tim |
1722 |
} |
146 |
|
|
|
147 |
tim |
1725 |
void VelocityVerletIntegrator::integrateStep() { } |
148 |
|
|
|
149 |
|
|
|
150 |
tim |
1722 |
void VelocityVerletIntegrator::thermalize() { |
151 |
|
|
|
152 |
tim |
1756 |
if (!info__->have_target_temp) { |
153 |
tim |
1725 |
sprintf(painCave.errMsg, |
154 |
|
|
"You can't resample the velocities without a targetTemp!\n"); |
155 |
|
|
painCave.isFatal = 1; |
156 |
|
|
painCave.severity = OOPSE_ERROR; |
157 |
|
|
simError(); |
158 |
|
|
return; |
159 |
|
|
} |
160 |
|
|
|
161 |
tim |
1722 |
} |
162 |
|
|
|
163 |
tim |
1725 |
void VelocityVerletIntegrator::calcForce(bool needPotential, |
164 |
|
|
bool needStress) { } |
165 |
tim |
1722 |
|
166 |
|
|
|
167 |
tim |
1725 |
} //end namespace oopse |