ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/Thermo.cpp
Revision: 1710
Committed: Fri May 18 21:44:02 2012 UTC (12 years, 11 months ago) by gezelter
File size: 14298 byte(s)
Log Message:
Added an adapter layer between the AtomType and the rest of the code to 
handle the bolt-on capabilities of new types. 

Fixed a long-standing bug in how storageLayout was being set to the maximum
possible value.

Started to add infrastructure for Polarizable and fluc-Q calculations.

File Contents

# Content
1 /*
2 * 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. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
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.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
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 #include <math.h>
44 #include <iostream>
45
46 #ifdef IS_MPI
47 #include <mpi.h>
48 #endif //is_mpi
49
50 #include "brains/Thermo.hpp"
51 #include "primitives/Molecule.hpp"
52 #include "utils/simError.h"
53 #include "utils/PhysicalConstants.hpp"
54 #include "types/MultipoleAdapter.hpp"
55
56 namespace OpenMD {
57
58 RealType Thermo::getKinetic() {
59 SimInfo::MoleculeIterator miter;
60 std::vector<StuntDouble*>::iterator iiter;
61 Molecule* mol;
62 StuntDouble* integrableObject;
63 Vector3d vel;
64 Vector3d angMom;
65 Mat3x3d I;
66 int i;
67 int j;
68 int k;
69 RealType mass;
70 RealType kinetic = 0.0;
71 RealType kinetic_global = 0.0;
72
73 for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
74 for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL;
75 integrableObject = mol->nextIntegrableObject(iiter)) {
76
77 mass = integrableObject->getMass();
78 vel = integrableObject->getVel();
79
80 kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
81
82 if (integrableObject->isDirectional()) {
83 angMom = integrableObject->getJ();
84 I = integrableObject->getI();
85
86 if (integrableObject->isLinear()) {
87 i = integrableObject->linearAxis();
88 j = (i + 1) % 3;
89 k = (i + 2) % 3;
90 kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
91 } else {
92 kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1)
93 + angMom[2]*angMom[2]/I(2, 2);
94 }
95 }
96
97 }
98 }
99
100 #ifdef IS_MPI
101
102 MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM,
103 MPI_COMM_WORLD);
104 kinetic = kinetic_global;
105
106 #endif //is_mpi
107
108 kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert;
109
110 return kinetic;
111 }
112
113 RealType Thermo::getPotential() {
114 RealType potential = 0.0;
115 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
116 RealType shortRangePot_local = curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
117
118 // Get total potential for entire system from MPI.
119
120 #ifdef IS_MPI
121
122 MPI_Allreduce(&shortRangePot_local, &potential, 1, MPI_REALTYPE, MPI_SUM,
123 MPI_COMM_WORLD);
124 potential += curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
125
126 #else
127
128 potential = shortRangePot_local + curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
129
130 #endif // is_mpi
131
132 return potential;
133 }
134
135 RealType Thermo::getTotalE() {
136 RealType total;
137
138 total = this->getKinetic() + this->getPotential();
139 return total;
140 }
141
142 RealType Thermo::getTemperature() {
143
144 RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* PhysicalConstants::kb );
145 return temperature;
146 }
147
148 RealType Thermo::getVolume() {
149 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
150 return curSnapshot->getVolume();
151 }
152
153 RealType Thermo::getPressure() {
154
155 // Relies on the calculation of the full molecular pressure tensor
156
157
158 Mat3x3d tensor;
159 RealType pressure;
160
161 tensor = getPressureTensor();
162
163 pressure = PhysicalConstants::pressureConvert * (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0;
164
165 return pressure;
166 }
167
168 RealType Thermo::getPressure(int direction) {
169
170 // Relies on the calculation of the full molecular pressure tensor
171
172
173 Mat3x3d tensor;
174 RealType pressure;
175
176 tensor = getPressureTensor();
177
178 pressure = PhysicalConstants::pressureConvert * tensor(direction, direction);
179
180 return pressure;
181 }
182
183 Mat3x3d Thermo::getPressureTensor() {
184 // returns pressure tensor in units amu*fs^-2*Ang^-1
185 // routine derived via viral theorem description in:
186 // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
187 Mat3x3d pressureTensor;
188 Mat3x3d p_local(0.0);
189 Mat3x3d p_global(0.0);
190
191 SimInfo::MoleculeIterator i;
192 std::vector<StuntDouble*>::iterator j;
193 Molecule* mol;
194 StuntDouble* integrableObject;
195 for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
196 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
197 integrableObject = mol->nextIntegrableObject(j)) {
198
199 RealType mass = integrableObject->getMass();
200 Vector3d vcom = integrableObject->getVel();
201 p_local += mass * outProduct(vcom, vcom);
202 }
203 }
204
205 #ifdef IS_MPI
206 MPI_Allreduce(p_local.getArrayPointer(), p_global.getArrayPointer(), 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
207 #else
208 p_global = p_local;
209 #endif // is_mpi
210
211 RealType volume = this->getVolume();
212 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
213 Mat3x3d tau = curSnapshot->getTau();
214
215 pressureTensor = (p_global + PhysicalConstants::energyConvert* tau)/volume;
216
217 return pressureTensor;
218 }
219
220
221 void Thermo::saveStat(){
222 Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
223 Stats& stat = currSnapshot->statData;
224
225 stat[Stats::KINETIC_ENERGY] = getKinetic();
226 stat[Stats::POTENTIAL_ENERGY] = getPotential();
227 stat[Stats::TOTAL_ENERGY] = stat[Stats::KINETIC_ENERGY] + stat[Stats::POTENTIAL_ENERGY] ;
228 stat[Stats::TEMPERATURE] = getTemperature();
229 stat[Stats::PRESSURE] = getPressure();
230 stat[Stats::VOLUME] = getVolume();
231
232 Mat3x3d tensor =getPressureTensor();
233 stat[Stats::PRESSURE_TENSOR_XX] = tensor(0, 0);
234 stat[Stats::PRESSURE_TENSOR_XY] = tensor(0, 1);
235 stat[Stats::PRESSURE_TENSOR_XZ] = tensor(0, 2);
236 stat[Stats::PRESSURE_TENSOR_YX] = tensor(1, 0);
237 stat[Stats::PRESSURE_TENSOR_YY] = tensor(1, 1);
238 stat[Stats::PRESSURE_TENSOR_YZ] = tensor(1, 2);
239 stat[Stats::PRESSURE_TENSOR_ZX] = tensor(2, 0);
240 stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);
241 stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);
242
243 // grab the simulation box dipole moment if specified
244 if (info_->getCalcBoxDipole()){
245 Vector3d totalDipole = getBoxDipole();
246 stat[Stats::BOX_DIPOLE_X] = totalDipole(0);
247 stat[Stats::BOX_DIPOLE_Y] = totalDipole(1);
248 stat[Stats::BOX_DIPOLE_Z] = totalDipole(2);
249 }
250
251 Globals* simParams = info_->getSimParams();
252
253 if (simParams->haveTaggedAtomPair() &&
254 simParams->havePrintTaggedPairDistance()) {
255 if ( simParams->getPrintTaggedPairDistance()) {
256
257 std::pair<int, int> tap = simParams->getTaggedAtomPair();
258 Vector3d pos1, pos2, rab;
259
260 #ifdef IS_MPI
261 std::cerr << "tap = " << tap.first << " " << tap.second << std::endl;
262
263 int mol1 = info_->getGlobalMolMembership(tap.first);
264 int mol2 = info_->getGlobalMolMembership(tap.second);
265 std::cerr << "mols = " << mol1 << " " << mol2 << std::endl;
266
267 int proc1 = info_->getMolToProc(mol1);
268 int proc2 = info_->getMolToProc(mol2);
269
270 std::cerr << " procs = " << proc1 << " " <<proc2 <<std::endl;
271
272 RealType data[3];
273 if (proc1 == worldRank) {
274 StuntDouble* sd1 = info_->getIOIndexToIntegrableObject(tap.first);
275 std::cerr << " on proc " << proc1 << ", sd1 has global index= " << sd1->getGlobalIndex() << std::endl;
276 pos1 = sd1->getPos();
277 data[0] = pos1.x();
278 data[1] = pos1.y();
279 data[2] = pos1.z();
280 MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
281 } else {
282 MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
283 pos1 = Vector3d(data);
284 }
285
286
287 if (proc2 == worldRank) {
288 StuntDouble* sd2 = info_->getIOIndexToIntegrableObject(tap.second);
289 std::cerr << " on proc " << proc2 << ", sd2 has global index= " << sd2->getGlobalIndex() << std::endl;
290 pos2 = sd2->getPos();
291 data[0] = pos2.x();
292 data[1] = pos2.y();
293 data[2] = pos2.z();
294 MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
295 } else {
296 MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
297 pos2 = Vector3d(data);
298 }
299 #else
300 StuntDouble* at1 = info_->getIOIndexToIntegrableObject(tap.first);
301 StuntDouble* at2 = info_->getIOIndexToIntegrableObject(tap.second);
302 pos1 = at1->getPos();
303 pos2 = at2->getPos();
304 #endif
305 rab = pos2 - pos1;
306 currSnapshot->wrapVector(rab);
307 stat[Stats::TAGGED_PAIR_DISTANCE] = rab.length();
308 }
309 }
310
311 /**@todo need refactorying*/
312 //Conserved Quantity is set by integrator and time is set by setTime
313
314 }
315
316
317 Vector3d Thermo::getBoxDipole() {
318 Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
319 SimInfo::MoleculeIterator miter;
320 std::vector<Atom*>::iterator aiter;
321 Molecule* mol;
322 Atom* atom;
323 RealType charge;
324 RealType moment(0.0);
325 Vector3d ri(0.0);
326 Vector3d dipoleVector(0.0);
327 Vector3d nPos(0.0);
328 Vector3d pPos(0.0);
329 RealType nChg(0.0);
330 RealType pChg(0.0);
331 int nCount = 0;
332 int pCount = 0;
333
334 RealType chargeToC = 1.60217733e-19;
335 RealType angstromToM = 1.0e-10;
336 RealType debyeToCm = 3.33564095198e-30;
337
338 for (mol = info_->beginMolecule(miter); mol != NULL;
339 mol = info_->nextMolecule(miter)) {
340
341 for (atom = mol->beginAtom(aiter); atom != NULL;
342 atom = mol->nextAtom(aiter)) {
343
344 if (atom->isCharge() ) {
345 charge = 0.0;
346 GenericData* data = atom->getAtomType()->getPropertyByName("Charge");
347 if (data != NULL) {
348
349 charge = (dynamic_cast<DoubleGenericData*>(data))->getData();
350 charge *= chargeToC;
351
352 ri = atom->getPos();
353 currSnapshot->wrapVector(ri);
354 ri *= angstromToM;
355
356 if (charge < 0.0) {
357 nPos += ri;
358 nChg -= charge;
359 nCount++;
360 } else if (charge > 0.0) {
361 pPos += ri;
362 pChg += charge;
363 pCount++;
364 }
365 }
366 }
367
368 MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
369 if (ma.isDipole() ) {
370 Vector3d u_i = atom->getElectroFrame().getColumn(2);
371 moment = ma.getDipoleMoment();
372 moment *= debyeToCm;
373 dipoleVector += u_i * moment;
374 }
375 }
376 }
377
378
379 #ifdef IS_MPI
380 RealType pChg_global, nChg_global;
381 int pCount_global, nCount_global;
382 Vector3d pPos_global, nPos_global, dipVec_global;
383
384 MPI_Allreduce(&pChg, &pChg_global, 1, MPI_REALTYPE, MPI_SUM,
385 MPI_COMM_WORLD);
386 pChg = pChg_global;
387 MPI_Allreduce(&nChg, &nChg_global, 1, MPI_REALTYPE, MPI_SUM,
388 MPI_COMM_WORLD);
389 nChg = nChg_global;
390 MPI_Allreduce(&pCount, &pCount_global, 1, MPI_INTEGER, MPI_SUM,
391 MPI_COMM_WORLD);
392 pCount = pCount_global;
393 MPI_Allreduce(&nCount, &nCount_global, 1, MPI_INTEGER, MPI_SUM,
394 MPI_COMM_WORLD);
395 nCount = nCount_global;
396 MPI_Allreduce(pPos.getArrayPointer(), pPos_global.getArrayPointer(), 3,
397 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
398 pPos = pPos_global;
399 MPI_Allreduce(nPos.getArrayPointer(), nPos_global.getArrayPointer(), 3,
400 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
401 nPos = nPos_global;
402 MPI_Allreduce(dipoleVector.getArrayPointer(),
403 dipVec_global.getArrayPointer(), 3,
404 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
405 dipoleVector = dipVec_global;
406 #endif //is_mpi
407
408 // first load the accumulated dipole moment (if dipoles were present)
409 Vector3d boxDipole = dipoleVector;
410 // now include the dipole moment due to charges
411 // use the lesser of the positive and negative charge totals
412 RealType chg_value = nChg <= pChg ? nChg : pChg;
413
414 // find the average positions
415 if (pCount > 0 && nCount > 0 ) {
416 pPos /= pCount;
417 nPos /= nCount;
418 }
419
420 // dipole is from the negative to the positive (physics notation)
421 boxDipole += (pPos - nPos) * chg_value;
422
423 return boxDipole;
424 }
425 } //end namespace OpenMD

Properties

Name Value
svn:keywords Author Id Revision Date