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 |
9 |
> |
* 1. Redistributions of source code must retain the above copyright |
10 |
|
* notice, this list of conditions and the following disclaimer. |
11 |
|
* |
12 |
< |
* 3. Redistributions in binary form must reproduce the above copyright |
12 |
> |
* 2. Redistributions in binary form must reproduce the above copyright |
13 |
|
* notice, this list of conditions and the following disclaimer in the |
14 |
|
* documentation and/or other materials provided with the |
15 |
|
* distribution. |
28 |
|
* arising out of the use of or inability to use software, even if the |
29 |
|
* University of Notre Dame has been advised of the possibility of |
30 |
|
* such damages. |
31 |
+ |
* |
32 |
+ |
* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your |
33 |
+ |
* research, please cite the appropriate papers when you publish your |
34 |
+ |
* work. Good starting points are: |
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] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 |
+ |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 |
|
*/ |
42 |
|
|
43 |
|
/** |
54 |
|
#include <vector> |
55 |
|
|
56 |
|
#include "brains/DataStorage.hpp" |
57 |
+ |
#include "nonbonded/NonBondedInteraction.hpp" |
58 |
|
#include "brains/Stats.hpp" |
57 |
– |
#include "UseTheForce/DarkSide/simulation_interface.h" |
59 |
|
|
60 |
< |
namespace oopse{ |
60 |
> |
namespace OpenMD{ |
61 |
|
|
62 |
+ |
struct FrameData { |
63 |
+ |
int id; /**< identification number of the snapshot */ |
64 |
+ |
RealType currentTime; /**< current time */ |
65 |
+ |
Mat3x3d hmat; /**< axes of the periodic box in matrix form */ |
66 |
+ |
Mat3x3d invHmat; /**< the inverse of the Hmat matrix */ |
67 |
+ |
bool orthoRhombic; /**< is this an orthorhombic periodic box? */ |
68 |
+ |
RealType volume; /**< total volume of this frame */ |
69 |
+ |
RealType pressure; /**< pressure of this frame */ |
70 |
+ |
RealType totalEnergy; /**< total energy of this frame */ |
71 |
+ |
RealType kineticEnergy; /**< kinetic energy of this frame */ |
72 |
+ |
RealType potentialEnergy; /**< potential energy of this frame */ |
73 |
+ |
RealType shortRangePotential; /**< short-range contributions to the potential*/ |
74 |
+ |
RealType longRangePotential; /**< long-range contributions to the potential */ |
75 |
+ |
RealType bondPotential; /**< bonded contribution to the potential */ |
76 |
+ |
RealType bendPotential; /**< angle-bending contribution to the potential */ |
77 |
+ |
RealType torsionPotential; /**< dihedral (torsion angle) contribution to the potential */ |
78 |
+ |
RealType inversionPotential; /**< inversion (planarity) contribution to the potential */ |
79 |
+ |
potVec lrPotentials; /**< breakdown of long-range potentials by family */ |
80 |
+ |
RealType temperature; /**< temperature of this frame */ |
81 |
+ |
RealType chi; /**< thermostat velocity */ |
82 |
+ |
RealType integralOfChiDt; /**< the actual thermostat */ |
83 |
+ |
RealType electronicTemperature; /**< temperature of the electronic degrees of freedom */ |
84 |
+ |
RealType chiQ; /**< fluctuating charge thermostat velocity */ |
85 |
+ |
RealType integralOfChiQDt; /**< the actual fluctuating charge thermostat */ |
86 |
+ |
Mat3x3d eta; /**< barostat matrix */ |
87 |
+ |
Vector3d COM; /**< location of center of mass */ |
88 |
+ |
Vector3d COMvel; /**< system center of mass velocity */ |
89 |
+ |
Vector3d COMw; /**< system center of mass angular velocity */ |
90 |
+ |
Mat3x3d stressTensor; /**< stress tensor */ |
91 |
+ |
Mat3x3d pressureTensor; /**< pressure tensor */ |
92 |
+ |
Vector3d systemDipole; /**< total system dipole moment */ |
93 |
+ |
Vector3d conductiveHeatFlux; /**< heat flux vector (conductive only) */ |
94 |
+ |
}; |
95 |
+ |
|
96 |
+ |
|
97 |
|
/** |
98 |
|
* @class Snapshot Snapshot.hpp "brains/Snapshot.hpp" |
99 |
|
* @brief Snapshot class is a repository class for storing dynamic data during |
100 |
|
* Simulation |
101 |
< |
* Every snapshot class will contain one DataStorage for atoms and one DataStorage |
101 |
> |
* Every snapshot class will contain one DataStorage for atoms and one DataStorage |
102 |
|
* for rigid bodies. |
103 |
|
*/ |
104 |
|
class Snapshot { |
105 |
|
public: |
106 |
|
|
107 |
< |
Snapshot(int nAtoms, int nRigidbodies) : atomData(nAtoms), rigidbodyData(nRigidbodies), |
108 |
< |
currentTime_(0), orthoRhombic_(0), chi_(0.0), integralOfChiDt_(0.0), eta_(0.0), id_(-1) { |
109 |
< |
|
107 |
> |
Snapshot(int nAtoms, int nRigidbodies, |
108 |
> |
int nCutoffGroups) : atomData(nAtoms), |
109 |
> |
rigidbodyData(nRigidbodies), |
110 |
> |
cgData(nCutoffGroups, DataStorage::dslPosition), |
111 |
> |
orthoTolerance_(1e-6), hasCOM_(false), hasVolume_(false){ |
112 |
> |
|
113 |
> |
frameData.id = -1; |
114 |
> |
frameData.currentTime = 0; |
115 |
> |
frameData.hmat = Mat3x3d(0.0); |
116 |
> |
frameData.invHmat = Mat3x3d(0.0); |
117 |
> |
frameData.orthoRhombic = false; |
118 |
> |
frameData.volume = 0.0; |
119 |
> |
frameData.pressure = 0.0; |
120 |
> |
frameData.totalEnergy = 0.0; |
121 |
> |
frameData.kineticEnergy = 0.0; |
122 |
> |
frameData.potentialEnergy = 0.0; |
123 |
> |
frameData.temperature = 0.0; |
124 |
> |
frameData.chi = 0.0; |
125 |
> |
frameData.integralOfChiDt = 0.0; |
126 |
> |
frameData.electronicTemperature = 0.0; |
127 |
> |
frameData.chiQ = 0.0; |
128 |
> |
frameData.integralOfChiQDt = 0.0; |
129 |
> |
frameData.eta = Mat3x3d(0.0); |
130 |
> |
frameData.COM = V3Zero; |
131 |
> |
frameData.COMvel = V3Zero; |
132 |
> |
frameData.COMw = V3Zero; |
133 |
> |
frameData.stressTensor = Mat3x3d(0.0); |
134 |
> |
frameData.pressureTensor = Mat3x3d(0.0); |
135 |
> |
frameData.systemDipole = Vector3d(0.0); |
136 |
> |
frameData.conductiveHeatFlux = Vector3d(0.0, 0.0, 0.0); |
137 |
|
} |
138 |
|
|
139 |
< |
Snapshot(int nAtoms, int nRigidbodies, int storageLayout) |
140 |
< |
: atomData(nAtoms, storageLayout), rigidbodyData(nRigidbodies, storageLayout), |
141 |
< |
currentTime_(0), orthoRhombic_(0), chi_(0.0), integralOfChiDt_(0.0), eta_(0.0), id_(-1) { |
142 |
< |
|
143 |
< |
} |
144 |
< |
|
139 |
> |
Snapshot(int nAtoms, int nRigidbodies, int nCutoffGroups, |
140 |
> |
int storageLayout) : atomData(nAtoms, storageLayout), |
141 |
> |
rigidbodyData(nRigidbodies, storageLayout), |
142 |
> |
cgData(nCutoffGroups, DataStorage::dslPosition), |
143 |
> |
orthoTolerance_(1e-6), |
144 |
> |
hasCOM_(false), |
145 |
> |
hasVolume_(false) { |
146 |
> |
frameData.id = -1; |
147 |
> |
frameData.currentTime = 0; |
148 |
> |
frameData.hmat = Mat3x3d(0.0); |
149 |
> |
frameData.invHmat = Mat3x3d(0.0); |
150 |
> |
frameData.orthoRhombic = false; |
151 |
> |
frameData.volume = 0.0; |
152 |
> |
frameData.pressure = 0.0; |
153 |
> |
frameData.totalEnergy = 0.0; |
154 |
> |
frameData.kineticEnergy = 0.0; |
155 |
> |
frameData.potentialEnergy = 0.0; |
156 |
> |
frameData.temperature = 0.0; |
157 |
> |
frameData.chi = 0.0; |
158 |
> |
frameData.integralOfChiDt = 0.0; |
159 |
> |
frameData.electronicTemperature = 0.0; |
160 |
> |
frameData.chiQ = 0.0; |
161 |
> |
frameData.integralOfChiQDt = 0.0; |
162 |
> |
frameData.eta = Mat3x3d(0.0); |
163 |
> |
frameData.COM = V3Zero; |
164 |
> |
frameData.COMvel = V3Zero; |
165 |
> |
frameData.COMw = V3Zero; |
166 |
> |
frameData.stressTensor = Mat3x3d(0.0); |
167 |
> |
frameData.pressureTensor = Mat3x3d(0.0); |
168 |
> |
frameData.systemDipole = V3Zero; |
169 |
> |
frameData.conductiveHeatFlux = Vector3d(0.0, 0.0, 0.0); |
170 |
> |
} |
171 |
> |
|
172 |
|
/** Returns the id of this Snapshot */ |
173 |
|
int getID() { |
174 |
< |
return id_; |
174 |
> |
return frameData.id; |
175 |
|
} |
176 |
|
|
177 |
|
/** Sets the id of this Snapshot */ |
178 |
|
void setID(int id) { |
179 |
< |
id_ = id; |
179 |
> |
frameData.id = id; |
180 |
|
} |
181 |
|
|
182 |
|
int getSize() { |
193 |
|
return rigidbodyData.getSize(); |
194 |
|
} |
195 |
|
|
196 |
+ |
/** Returns the number of rigid bodies */ |
197 |
+ |
int getNumberOfCutoffGroups() { |
198 |
+ |
return cgData.getSize(); |
199 |
+ |
} |
200 |
+ |
|
201 |
|
/** Returns the H-Matrix */ |
202 |
|
Mat3x3d getHmat() { |
203 |
< |
return hmat_; |
203 |
> |
return frameData.hmat; |
204 |
|
} |
205 |
|
|
206 |
|
/** Sets the H-Matrix */ |
207 |
|
void setHmat(const Mat3x3d& m); |
208 |
|
|
209 |
|
RealType getVolume() { |
210 |
< |
return hmat_.determinant(); |
210 |
> |
if (hasVolume_){ |
211 |
> |
return frameData.volume; |
212 |
> |
}else{ |
213 |
> |
return frameData.hmat.determinant(); |
214 |
> |
} |
215 |
|
} |
216 |
|
|
217 |
+ |
void setVolume(RealType volume){ |
218 |
+ |
hasVolume_=true; |
219 |
+ |
frameData.volume = volume; |
220 |
+ |
} |
221 |
+ |
|
222 |
|
/** Returns the inverse H-Matrix */ |
223 |
|
Mat3x3d getInvHmat() { |
224 |
< |
return invHmat_; |
224 |
> |
return frameData.invHmat; |
225 |
|
} |
226 |
|
|
227 |
|
/** Wrapping the vector according to periodic boundary condition*/ |
228 |
|
void wrapVector(Vector3d& v); |
229 |
+ |
/** Scaling a vector to multiples of the periodic box */ |
230 |
+ |
Vector3d scaleVector(Vector3d &v); |
231 |
|
|
232 |
+ |
|
233 |
+ |
Vector3d getCOM(); |
234 |
+ |
Vector3d getCOMvel(); |
235 |
+ |
Vector3d getCOMw(); |
236 |
|
|
237 |
|
RealType getTime() { |
238 |
< |
return currentTime_; |
238 |
> |
return frameData.currentTime; |
239 |
|
} |
240 |
|
|
241 |
|
void increaseTime(RealType dt) { |
243 |
|
} |
244 |
|
|
245 |
|
void setTime(RealType time) { |
246 |
< |
currentTime_ =time; |
246 |
> |
frameData.currentTime =time; |
247 |
|
//time at statData is redundant |
248 |
< |
statData[Stats::TIME] = currentTime_; |
248 |
> |
statData[Stats::TIME] = frameData.currentTime; |
249 |
|
} |
250 |
|
|
251 |
|
RealType getChi() { |
252 |
< |
return chi_; |
252 |
> |
return frameData.chi; |
253 |
|
} |
254 |
|
|
255 |
|
void setChi(RealType chi) { |
256 |
< |
chi_ = chi; |
256 |
> |
frameData.chi = chi; |
257 |
|
} |
258 |
|
|
259 |
|
RealType getIntegralOfChiDt() { |
260 |
< |
return integralOfChiDt_; |
260 |
> |
return frameData.integralOfChiDt; |
261 |
|
} |
262 |
|
|
263 |
|
void setIntegralOfChiDt(RealType integralOfChiDt) { |
264 |
< |
integralOfChiDt_ = integralOfChiDt; |
264 |
> |
frameData.integralOfChiDt = integralOfChiDt; |
265 |
|
} |
266 |
|
|
267 |
+ |
RealType getChiElectronic() { |
268 |
+ |
return frameData.chiQ; |
269 |
+ |
} |
270 |
+ |
|
271 |
+ |
void setChiElectronic(RealType chiQ) { |
272 |
+ |
frameData.chiQ = chiQ; |
273 |
+ |
} |
274 |
+ |
|
275 |
+ |
RealType getIntegralOfChiElectronicDt() { |
276 |
+ |
return frameData.integralOfChiQDt; |
277 |
+ |
} |
278 |
+ |
|
279 |
+ |
void setIntegralOfChiElectronicDt(RealType integralOfChiQDt) { |
280 |
+ |
frameData.integralOfChiQDt = integralOfChiQDt; |
281 |
+ |
} |
282 |
+ |
|
283 |
+ |
|
284 |
+ |
void setOrthoTolerance(RealType orthoTolerance) { |
285 |
+ |
orthoTolerance_ = orthoTolerance; |
286 |
+ |
} |
287 |
+ |
|
288 |
|
Mat3x3d getEta() { |
289 |
< |
return eta_; |
289 |
> |
return frameData.eta; |
290 |
|
} |
291 |
|
|
292 |
|
void setEta(const Mat3x3d& eta) { |
293 |
< |
eta_ = eta; |
293 |
> |
frameData.eta = eta; |
294 |
|
} |
295 |
< |
|
295 |
> |
|
296 |
> |
Mat3x3d getStressTensor() { |
297 |
> |
return frameData.stressTensor; |
298 |
> |
} |
299 |
> |
|
300 |
> |
void setStressTensor(const Mat3x3d& stressTensor) { |
301 |
> |
frameData.stressTensor = stressTensor; |
302 |
> |
} |
303 |
> |
|
304 |
> |
Vector3d getConductiveHeatFlux() { |
305 |
> |
return frameData.conductiveHeatFlux; |
306 |
> |
} |
307 |
> |
|
308 |
> |
void setConductiveHeatFlux(const Vector3d& heatFlux) { |
309 |
> |
frameData.conductiveHeatFlux = heatFlux; |
310 |
> |
} |
311 |
> |
|
312 |
> |
bool hasCOM() { |
313 |
> |
return hasCOM_; |
314 |
> |
} |
315 |
> |
|
316 |
> |
void setCOMprops(const Vector3d& COM, const Vector3d& COMvel, const Vector3d& COMw) { |
317 |
> |
frameData.COM = COM; |
318 |
> |
frameData.COMvel = COMvel; |
319 |
> |
frameData.COMw = COMw; |
320 |
> |
hasCOM_ = true; |
321 |
> |
} |
322 |
> |
|
323 |
|
DataStorage atomData; |
324 |
|
DataStorage rigidbodyData; |
325 |
+ |
DataStorage cgData; |
326 |
+ |
FrameData frameData; |
327 |
|
Stats statData; |
168 |
– |
|
169 |
– |
private: |
170 |
– |
RealType currentTime_; |
328 |
|
|
329 |
< |
Mat3x3d hmat_; |
330 |
< |
Mat3x3d invHmat_; |
331 |
< |
int orthoRhombic_; |
332 |
< |
|
176 |
< |
RealType chi_; |
177 |
< |
RealType integralOfChiDt_; |
178 |
< |
Mat3x3d eta_; |
179 |
< |
|
180 |
< |
int id_; /**< identification number of the snapshot */ |
329 |
> |
private: |
330 |
> |
RealType orthoTolerance_; |
331 |
> |
bool hasCOM_; |
332 |
> |
bool hasVolume_; |
333 |
|
}; |
334 |
|
|
335 |
|
typedef DataStorage (Snapshot::*DataStoragePointer); |