ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/brains/Thermo.cpp
(Generate patch)

Comparing trunk/src/brains/Thermo.cpp (file contents):
Revision 1666 by chuckv, Wed Dec 14 20:21:54 2011 UTC vs.
Revision 2046 by gezelter, Tue Dec 2 22:11:04 2014 UTC

# Line 32 | Line 32
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]  Vardeman & Gezelter, in progress (2009).
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, 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  
42 – #include <math.h>
43 – #include <iostream>
44 –
43   #ifdef IS_MPI
44   #include <mpi.h>
45   #endif //is_mpi
46 +
47 + #include <math.h>
48 + #include <iostream>
49  
50   #include "brains/Thermo.hpp"
51   #include "primitives/Molecule.hpp"
52   #include "utils/simError.h"
53   #include "utils/PhysicalConstants.hpp"
54 + #include "types/FixedChargeAdapter.hpp"
55 + #include "types/FluctuatingChargeAdapter.hpp"
56 + #include "types/MultipoleAdapter.hpp"
57 + #ifdef HAVE_QHULL
58 + #include "math/ConvexHull.hpp"
59 + #include "math/AlphaHull.hpp"
60 + #endif
61  
62 + using namespace std;
63   namespace OpenMD {
64  
65 <  RealType Thermo::getKinetic() {
66 <    SimInfo::MoleculeIterator miter;
58 <    std::vector<StuntDouble*>::iterator iiter;
59 <    Molecule* mol;
60 <    StuntDouble* integrableObject;
61 <    Vector3d vel;
62 <    Vector3d angMom;
63 <    Mat3x3d I;
64 <    int i;
65 <    int j;
66 <    int k;
67 <    RealType mass;
68 <    RealType kinetic = 0.0;
69 <    RealType kinetic_global = 0.0;
65 >  RealType Thermo::getTranslationalKinetic() {
66 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
67  
68 <    for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
69 <      for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL;
70 <     integrableObject = mol->nextIntegrableObject(iiter)) {
71 <
72 <  mass = integrableObject->getMass();
73 <  vel = integrableObject->getVel();
74 <
75 <  kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
76 <
77 <  if (integrableObject->isDirectional()) {
78 <    angMom = integrableObject->getJ();
79 <    I = integrableObject->getI();
80 <
81 <    if (integrableObject->isLinear()) {
82 <      i = integrableObject->linearAxis();
83 <      j = (i + 1) % 3;
84 <      k = (i + 2) % 3;
85 <      kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
86 <    } else {
87 <      kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1)
88 <        + angMom[2]*angMom[2]/I(2, 2);
68 >    if (!snap->hasTranslationalKineticEnergy) {
69 >      SimInfo::MoleculeIterator miter;
70 >      vector<StuntDouble*>::iterator iiter;
71 >      Molecule* mol;
72 >      StuntDouble* sd;    
73 >      Vector3d vel;
74 >      RealType mass;
75 >      RealType kinetic(0.0);
76 >      
77 >      for (mol = info_->beginMolecule(miter); mol != NULL;
78 >           mol = info_->nextMolecule(miter)) {
79 >        
80 >        for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
81 >             sd = mol->nextIntegrableObject(iiter)) {
82 >          
83 >          mass = sd->getMass();
84 >          vel = sd->getVel();
85 >          
86 >          kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
87 >          
88 >        }
89 >      }
90 >      
91 > #ifdef IS_MPI
92 >      MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE,
93 >                    MPI_SUM, MPI_COMM_WORLD);
94 > #endif
95 >      
96 >      kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert;
97 >      
98 >      
99 >      snap->setTranslationalKineticEnergy(kinetic);
100      }
101 +    return snap->getTranslationalKineticEnergy();
102    }
103  
104 +  RealType Thermo::getRotationalKinetic() {
105 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
106 +
107 +    if (!snap->hasRotationalKineticEnergy) {
108 +      SimInfo::MoleculeIterator miter;
109 +      vector<StuntDouble*>::iterator iiter;
110 +      Molecule* mol;
111 +      StuntDouble* sd;    
112 +      Vector3d angMom;
113 +      Mat3x3d I;
114 +      int i, j, k;
115 +      RealType kinetic(0.0);
116 +      
117 +      for (mol = info_->beginMolecule(miter); mol != NULL;
118 +           mol = info_->nextMolecule(miter)) {
119 +        
120 +        for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
121 +             sd = mol->nextIntegrableObject(iiter)) {
122 +          
123 +          if (sd->isDirectional()) {
124 +            angMom = sd->getJ();
125 +            I = sd->getI();
126 +            
127 +            if (sd->isLinear()) {
128 +              i = sd->linearAxis();
129 +              j = (i + 1) % 3;
130 +              k = (i + 2) % 3;
131 +              kinetic += angMom[j] * angMom[j] / I(j, j)
132 +                + angMom[k] * angMom[k] / I(k, k);
133 +            } else {                        
134 +              kinetic += angMom[0]*angMom[0]/I(0, 0)
135 +                + angMom[1]*angMom[1]/I(1, 1)
136 +                + angMom[2]*angMom[2]/I(2, 2);
137 +            }
138 +          }          
139 +        }
140        }
141 +      
142 + #ifdef IS_MPI
143 +      MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE,
144 +                    MPI_SUM, MPI_COMM_WORLD);
145 + #endif
146 +      
147 +      kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert;
148 +          
149 +      snap->setRotationalKineticEnergy(kinetic);
150      }
151 +    return snap->getRotationalKineticEnergy();
152 +  }
153  
154 < #ifdef IS_MPI
154 >      
155  
156 <    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM,
157 <                  MPI_COMM_WORLD);
102 <    kinetic = kinetic_global;
156 >  RealType Thermo::getKinetic() {
157 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
158  
159 < #endif //is_mpi
160 <
161 <    kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert;
162 <
163 <    return kinetic;
159 >    if (!snap->hasKineticEnergy) {
160 >      RealType ke = getTranslationalKinetic() + getRotationalKinetic();
161 >      snap->setKineticEnergy(ke);
162 >    }
163 >    return snap->getKineticEnergy();
164    }
165  
166    RealType Thermo::getPotential() {
112 –    RealType potential = 0.0;
113 –    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
114 –    RealType shortRangePot_local =  curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
167  
168 <    // Get total potential for entire system from MPI.
168 >    // ForceManager computes the potential and stores it in the
169 >    // Snapshot.  All we have to do is report it.
170  
171 < #ifdef IS_MPI
171 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
172 >    return snap->getPotentialEnergy();
173 >  }
174  
175 <    MPI_Allreduce(&shortRangePot_local, &potential, 1, MPI_REALTYPE, MPI_SUM,
121 <                  MPI_COMM_WORLD);
122 <    potential += curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
175 >  RealType Thermo::getTotalEnergy() {
176  
177 < #else
177 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
178  
179 <    potential = shortRangePot_local + curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
179 >    if (!snap->hasTotalEnergy) {
180 >      snap->setTotalEnergy(this->getKinetic() + this->getPotential());
181 >    }
182  
183 < #endif // is_mpi
129 <
130 <    return potential;
183 >    return snap->getTotalEnergy();
184    }
185  
186 <  RealType Thermo::getTotalE() {
134 <    RealType total;
186 >  RealType Thermo::getTemperature() {
187  
188 <    total = this->getKinetic() + this->getPotential();
137 <    return total;
138 <  }
188 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
189  
190 <  RealType Thermo::getTemperature() {
190 >    if (!snap->hasTemperature) {
191  
192 <    RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* PhysicalConstants::kb );
193 <    return temperature;
144 <  }
192 >      RealType temperature = ( 2.0 * this->getKinetic() )
193 >        / (info_->getNdf()* PhysicalConstants::kb );
194  
195 <  RealType Thermo::getVolume() {
196 <    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
197 <    return curSnapshot->getVolume();
195 >      snap->setTemperature(temperature);
196 >    }
197 >    
198 >    return snap->getTemperature();
199    }
200  
201 <  RealType Thermo::getPressure() {
201 >  RealType Thermo::getElectronicTemperature() {
202 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
203  
204 <    // Relies on the calculation of the full molecular pressure tensor
204 >    if (!snap->hasElectronicTemperature) {
205 >      
206 >      SimInfo::MoleculeIterator miter;
207 >      vector<Atom*>::iterator iiter;
208 >      Molecule* mol;
209 >      Atom* atom;    
210 >      RealType cvel;
211 >      RealType cmass;
212 >      RealType kinetic(0.0);
213 >      RealType eTemp;
214 >      
215 >      for (mol = info_->beginMolecule(miter); mol != NULL;
216 >           mol = info_->nextMolecule(miter)) {
217 >        
218 >        for (atom = mol->beginFluctuatingCharge(iiter); atom != NULL;
219 >             atom = mol->nextFluctuatingCharge(iiter)) {
220 >          
221 >          cmass = atom->getChargeMass();
222 >          cvel = atom->getFlucQVel();
223 >          
224 >          kinetic += cmass * cvel * cvel;
225 >          
226 >        }
227 >      }
228 >    
229 > #ifdef IS_MPI
230 >      MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE,
231 >                    MPI_SUM, MPI_COMM_WORLD);
232 > #endif
233  
234 +      kinetic *= 0.5;
235 +      eTemp =  (2.0 * kinetic) /
236 +        (info_->getNFluctuatingCharges() * PhysicalConstants::kb );            
237 +    
238 +      snap->setElectronicTemperature(eTemp);
239 +    }
240  
241 <    Mat3x3d tensor;
242 <    RealType pressure;
241 >    return snap->getElectronicTemperature();
242 >  }
243  
159 –    tensor = getPressureTensor();
244  
245 <    pressure = PhysicalConstants::pressureConvert * (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0;
246 <
247 <    return pressure;
245 >  RealType Thermo::getVolume() {
246 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
247 >    return snap->getVolume();
248    }
249  
250 <  RealType Thermo::getPressure(int direction) {
250 >  RealType Thermo::getPressure() {
251 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
252  
253 <    // Relies on the calculation of the full molecular pressure tensor
254 <
255 <
256 <    Mat3x3d tensor;
257 <    RealType pressure;
258 <
259 <    tensor = getPressureTensor();
260 <
261 <    pressure = PhysicalConstants::pressureConvert * tensor(direction, direction);
262 <
263 <    return pressure;
253 >    if (!snap->hasPressure) {
254 >      // Relies on the calculation of the full molecular pressure tensor
255 >      
256 >      Mat3x3d tensor;
257 >      RealType pressure;
258 >      
259 >      tensor = getPressureTensor();
260 >      
261 >      pressure = PhysicalConstants::pressureConvert *
262 >        (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0;
263 >      
264 >      snap->setPressure(pressure);
265 >    }
266 >    
267 >    return snap->getPressure();    
268    }
269  
270    Mat3x3d Thermo::getPressureTensor() {
271      // returns pressure tensor in units amu*fs^-2*Ang^-1
272      // routine derived via viral theorem description in:
273      // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
274 <    Mat3x3d pressureTensor;
186 <    Mat3x3d p_local(0.0);
187 <    Mat3x3d p_global(0.0);
274 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
275  
276 <    SimInfo::MoleculeIterator i;
190 <    std::vector<StuntDouble*>::iterator j;
191 <    Molecule* mol;
192 <    StuntDouble* integrableObject;
193 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
194 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
195 <     integrableObject = mol->nextIntegrableObject(j)) {
276 >    if (!snap->hasPressureTensor) {
277  
278 <  RealType mass = integrableObject->getMass();
279 <  Vector3d vcom = integrableObject->getVel();
280 <  p_local += mass * outProduct(vcom, vcom);
278 >      Mat3x3d pressureTensor;
279 >      Mat3x3d p_tens(0.0);
280 >      RealType mass;
281 >      Vector3d vcom;
282 >      
283 >      SimInfo::MoleculeIterator i;
284 >      vector<StuntDouble*>::iterator j;
285 >      Molecule* mol;
286 >      StuntDouble* sd;    
287 >      for (mol = info_->beginMolecule(i); mol != NULL;
288 >           mol = info_->nextMolecule(i)) {
289 >        
290 >        for (sd = mol->beginIntegrableObject(j); sd != NULL;
291 >             sd = mol->nextIntegrableObject(j)) {
292 >          
293 >          mass = sd->getMass();
294 >          vcom = sd->getVel();
295 >          p_tens += mass * outProduct(vcom, vcom);        
296 >        }
297 >      }
298 >      
299 > #ifdef IS_MPI
300 >      MPI_Allreduce(MPI_IN_PLACE, p_tens.getArrayPointer(), 9,
301 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
302 > #endif
303 >      
304 >      RealType volume = this->getVolume();
305 >      Mat3x3d stressTensor = snap->getStressTensor();
306 >      
307 >      pressureTensor =  (p_tens +
308 >                         PhysicalConstants::energyConvert * stressTensor)/volume;
309 >      
310 >      snap->setPressureTensor(pressureTensor);
311 >    }
312 >    return snap->getPressureTensor();
313 >  }
314 >
315 >
316 >
317 >
318 >  Vector3d Thermo::getSystemDipole() {
319 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
320 >
321 >    if (!snap->hasSystemDipole) {
322 >      SimInfo::MoleculeIterator miter;
323 >      vector<Atom*>::iterator aiter;
324 >      Molecule* mol;
325 >      Atom* atom;
326 >      RealType charge;
327 >      Vector3d ri(0.0);
328 >      Vector3d dipoleVector(0.0);
329 >      Vector3d nPos(0.0);
330 >      Vector3d pPos(0.0);
331 >      RealType nChg(0.0);
332 >      RealType pChg(0.0);
333 >      int nCount = 0;
334 >      int pCount = 0;
335 >      
336 >      RealType chargeToC = 1.60217733e-19;
337 >      RealType angstromToM = 1.0e-10;
338 >      RealType debyeToCm = 3.33564095198e-30;
339 >      
340 >      for (mol = info_->beginMolecule(miter); mol != NULL;
341 >           mol = info_->nextMolecule(miter)) {
342 >        
343 >        for (atom = mol->beginAtom(aiter); atom != NULL;
344 >             atom = mol->nextAtom(aiter)) {
345 >          
346 >          charge = 0.0;
347 >          
348 >          FixedChargeAdapter fca = FixedChargeAdapter(atom->getAtomType());
349 >          if ( fca.isFixedCharge() ) {
350 >            charge = fca.getCharge();
351 >          }
352 >          
353 >          FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atom->getAtomType());
354 >          if ( fqa.isFluctuatingCharge() ) {
355 >            charge += atom->getFlucQPos();
356 >          }
357 >          
358 >          charge *= chargeToC;
359 >          
360 >          ri = atom->getPos();
361 >          snap->wrapVector(ri);
362 >          ri *= angstromToM;
363 >          
364 >          if (charge < 0.0) {
365 >            nPos += ri;
366 >            nChg -= charge;
367 >            nCount++;
368 >          } else if (charge > 0.0) {
369 >            pPos += ri;
370 >            pChg += charge;
371 >            pCount++;
372 >          }
373 >          
374 >          if (atom->isDipole()) {
375 >            dipoleVector += atom->getDipole() * debyeToCm;
376 >          }
377 >        }
378        }
379 +      
380 +      
381 + #ifdef IS_MPI
382 +      MPI_Allreduce(MPI_IN_PLACE, &pChg, 1, MPI_REALTYPE,
383 +                    MPI_SUM, MPI_COMM_WORLD);
384 +      MPI_Allreduce(MPI_IN_PLACE, &nChg, 1, MPI_REALTYPE,
385 +                    MPI_SUM, MPI_COMM_WORLD);
386 +      
387 +      MPI_Allreduce(MPI_IN_PLACE, &pCount, 1, MPI_INTEGER,
388 +                    MPI_SUM, MPI_COMM_WORLD);
389 +      MPI_Allreduce(MPI_IN_PLACE, &nCount, 1, MPI_INTEGER,
390 +                    MPI_SUM, MPI_COMM_WORLD);
391 +      
392 +      MPI_Allreduce(MPI_IN_PLACE, pPos.getArrayPointer(), 3,
393 +                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
394 +      MPI_Allreduce(MPI_IN_PLACE, nPos.getArrayPointer(), 3,
395 +                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
396 +
397 +      MPI_Allreduce(MPI_IN_PLACE, dipoleVector.getArrayPointer(),
398 +                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
399 + #endif
400 +      
401 +      // first load the accumulated dipole moment (if dipoles were present)
402 +      Vector3d boxDipole = dipoleVector;
403 +      // now include the dipole moment due to charges
404 +      // use the lesser of the positive and negative charge totals
405 +      RealType chg_value = nChg <= pChg ? nChg : pChg;
406 +      
407 +      // find the average positions
408 +      if (pCount > 0 && nCount > 0 ) {
409 +        pPos /= pCount;
410 +        nPos /= nCount;
411 +      }
412 +      
413 +      // dipole is from the negative to the positive (physics notation)
414 +      boxDipole += (pPos - nPos) * chg_value;
415 +      snap->setSystemDipole(boxDipole);
416      }
417  
418 +    return snap->getSystemDipole();
419 +  }
420 +
421 +
422 +  Mat3x3d Thermo::getSystemQuadrupole() {
423 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
424 +
425 +    if (!snap->hasSystemQuadrupole) {
426 +      SimInfo::MoleculeIterator miter;
427 +      vector<Atom*>::iterator aiter;
428 +      Molecule* mol;
429 +      Atom* atom;
430 +      RealType charge;
431 +      Vector3d ri(0.0);
432 +      Vector3d dipole(0.0);
433 +      Mat3x3d qpole(0.0);
434 +      
435 +      RealType chargeToC = 1.60217733e-19;
436 +      RealType angstromToM = 1.0e-10;
437 +      RealType debyeToCm = 3.33564095198e-30;
438 +      
439 +      for (mol = info_->beginMolecule(miter); mol != NULL;
440 +           mol = info_->nextMolecule(miter)) {
441 +        
442 +        for (atom = mol->beginAtom(aiter); atom != NULL;
443 +             atom = mol->nextAtom(aiter)) {
444 +
445 +          ri = atom->getPos();
446 +          snap->wrapVector(ri);
447 +          ri *= angstromToM;
448 +          
449 +          charge = 0.0;
450 +          
451 +          FixedChargeAdapter fca = FixedChargeAdapter(atom->getAtomType());
452 +          if ( fca.isFixedCharge() ) {
453 +            charge = fca.getCharge();
454 +          }
455 +          
456 +          FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atom->getAtomType());
457 +          if ( fqa.isFluctuatingCharge() ) {
458 +            charge += atom->getFlucQPos();
459 +          }
460 +          
461 +          charge *= chargeToC;
462 +          
463 +          qpole += 0.5 * charge * outProduct(ri, ri);
464 +
465 +          MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
466 +          
467 +          if ( ma.isDipole() ) {
468 +            dipole = atom->getDipole() * debyeToCm;
469 +            qpole += 0.5 * outProduct( dipole, ri );
470 +            qpole += 0.5 * outProduct( ri, dipole );
471 +          }
472 +
473 +          if ( ma.isQuadrupole() ) {
474 +            qpole += atom->getQuadrupole() * debyeToCm * angstromToM;          
475 +          }
476 +        }
477 +      }
478 +        
479   #ifdef IS_MPI
480 <    MPI_Allreduce(p_local.getArrayPointer(), p_global.getArrayPointer(), 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
481 < #else
482 <    p_global = p_local;
483 < #endif // is_mpi
480 >      MPI_Allreduce(MPI_IN_PLACE, qpole.getArrayPointer(),
481 >                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
482 > #endif
483 >      
484 >      snap->setSystemQuadrupole(qpole);
485 >    }
486 >    
487 >    return snap->getSystemQuadrupole();
488 >  }
489  
490 <    RealType volume = this->getVolume();
491 <    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
492 <    Mat3x3d tau = curSnapshot->statData.getTau();
490 >  // Returns the Heat Flux Vector for the system
491 >  Vector3d Thermo::getHeatFlux(){
492 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
493 >    SimInfo::MoleculeIterator miter;
494 >    vector<StuntDouble*>::iterator iiter;
495 >    Molecule* mol;
496 >    StuntDouble* sd;    
497 >    RigidBody::AtomIterator ai;
498 >    Atom* atom;      
499 >    Vector3d vel;
500 >    Vector3d angMom;
501 >    Mat3x3d I;
502 >    int i;
503 >    int j;
504 >    int k;
505 >    RealType mass;
506  
507 <    pressureTensor =  (p_global + PhysicalConstants::energyConvert* tau)/volume;
507 >    Vector3d x_a;
508 >    RealType kinetic;
509 >    RealType potential;
510 >    RealType eatom;
511 >    // Convective portion of the heat flux
512 >    Vector3d heatFluxJc = V3Zero;
513  
514 <    return pressureTensor;
514 >    /* Calculate convective portion of the heat flux */
515 >    for (mol = info_->beginMolecule(miter); mol != NULL;
516 >         mol = info_->nextMolecule(miter)) {
517 >      
518 >      for (sd = mol->beginIntegrableObject(iiter);
519 >           sd != NULL;
520 >           sd = mol->nextIntegrableObject(iiter)) {
521 >        
522 >        mass = sd->getMass();
523 >        vel = sd->getVel();
524 >
525 >        kinetic = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
526 >        
527 >        if (sd->isDirectional()) {
528 >          angMom = sd->getJ();
529 >          I = sd->getI();
530 >
531 >          if (sd->isLinear()) {
532 >            i = sd->linearAxis();
533 >            j = (i + 1) % 3;
534 >            k = (i + 2) % 3;
535 >            kinetic += angMom[j] * angMom[j] / I(j, j)
536 >              + angMom[k] * angMom[k] / I(k, k);
537 >          } else {                        
538 >            kinetic += angMom[0]*angMom[0]/I(0, 0)
539 >              + angMom[1]*angMom[1]/I(1, 1)
540 >              + angMom[2]*angMom[2]/I(2, 2);
541 >          }
542 >        }
543 >
544 >        potential = 0.0;
545 >
546 >        if (sd->isRigidBody()) {
547 >          RigidBody* rb = dynamic_cast<RigidBody*>(sd);
548 >          for (atom = rb->beginAtom(ai); atom != NULL;
549 >               atom = rb->nextAtom(ai)) {
550 >            potential +=  atom->getParticlePot();
551 >          }          
552 >        } else {
553 >          potential = sd->getParticlePot();
554 >        }
555 >
556 >        potential *= PhysicalConstants::energyConvert; // amu A^2/fs^2
557 >        // The potential may not be a 1/2 factor
558 >        eatom = (kinetic + potential)/2.0;  // amu A^2/fs^2
559 >        heatFluxJc[0] += eatom*vel[0]; // amu A^3/fs^3
560 >        heatFluxJc[1] += eatom*vel[1]; // amu A^3/fs^3
561 >        heatFluxJc[2] += eatom*vel[2]; // amu A^3/fs^3
562 >      }
563 >    }
564 >
565 >    /* The J_v vector is reduced in the forceManager so everyone has
566 >     *  the global Jv. Jc is computed over the local atoms and must be
567 >     *  reduced among all processors.
568 >     */
569 > #ifdef IS_MPI
570 >    MPI_Allreduce(MPI_IN_PLACE, &heatFluxJc[0], 3, MPI_REALTYPE,
571 >                  MPI_SUM, MPI_COMM_WORLD);
572 > #endif
573 >    
574 >    // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3
575 >
576 >    Vector3d heatFluxJv = currSnapshot->getConductiveHeatFlux() *
577 >      PhysicalConstants::energyConvert;
578 >        
579 >    // Correct for the fact the flux is 1/V (Jc + Jv)
580 >    return (heatFluxJv + heatFluxJc) / this->getVolume(); // amu / fs^3
581    }
582  
583  
584 <  void Thermo::saveStat(){
585 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
221 <    Stats& stat = currSnapshot->statData;
584 >  Vector3d Thermo::getComVel(){
585 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
586  
587 <    stat[Stats::KINETIC_ENERGY] = getKinetic();
224 <    stat[Stats::POTENTIAL_ENERGY] = getPotential();
225 <    stat[Stats::TOTAL_ENERGY] = stat[Stats::KINETIC_ENERGY]  + stat[Stats::POTENTIAL_ENERGY] ;
226 <    stat[Stats::TEMPERATURE] = getTemperature();
227 <    stat[Stats::PRESSURE] = getPressure();
228 <    stat[Stats::VOLUME] = getVolume();
587 >    if (!snap->hasCOMvel) {
588  
589 <    Mat3x3d tensor =getPressureTensor();
590 <    stat[Stats::PRESSURE_TENSOR_XX] = tensor(0, 0);
591 <    stat[Stats::PRESSURE_TENSOR_XY] = tensor(0, 1);
592 <    stat[Stats::PRESSURE_TENSOR_XZ] = tensor(0, 2);
593 <    stat[Stats::PRESSURE_TENSOR_YX] = tensor(1, 0);
594 <    stat[Stats::PRESSURE_TENSOR_YY] = tensor(1, 1);
595 <    stat[Stats::PRESSURE_TENSOR_YZ] = tensor(1, 2);
596 <    stat[Stats::PRESSURE_TENSOR_ZX] = tensor(2, 0);
597 <    stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);
598 <    stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);
599 <    Vector3d GKappa_t = getThermalHelfand();
600 <    stat[Stats::THERMAL_HELFANDMOMENT_X] = GKappa_t.x();
601 <    stat[Stats::THERMAL_HELFANDMOMENT_Y] = GKappa_t.y();
602 <    stat[Stats::THERMAL_HELFANDMOMENT_Z] = GKappa_t.z();
589 >      SimInfo::MoleculeIterator i;
590 >      Molecule* mol;
591 >      
592 >      Vector3d comVel(0.0);
593 >      RealType totalMass(0.0);
594 >      
595 >      for (mol = info_->beginMolecule(i); mol != NULL;
596 >           mol = info_->nextMolecule(i)) {
597 >        RealType mass = mol->getMass();
598 >        totalMass += mass;
599 >        comVel += mass * mol->getComVel();
600 >      }  
601 >      
602 > #ifdef IS_MPI
603 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
604 >                    MPI_SUM, MPI_COMM_WORLD);
605 >      MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3,
606 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
607 > #endif
608 >      
609 >      comVel /= totalMass;
610 >      snap->setCOMvel(comVel);
611 >    }
612 >    return snap->getCOMvel();
613 >  }
614  
615 <    Globals* simParams = info_->getSimParams();
615 >  Vector3d Thermo::getCom(){
616 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
617  
618 <    if (simParams->haveTaggedAtomPair() &&
619 <        simParams->havePrintTaggedPairDistance()) {
620 <      if ( simParams->getPrintTaggedPairDistance()) {
618 >    if (!snap->hasCOM) {
619 >      
620 >      SimInfo::MoleculeIterator i;
621 >      Molecule* mol;
622 >      
623 >      Vector3d com(0.0);
624 >      RealType totalMass(0.0);
625 >      
626 >      for (mol = info_->beginMolecule(i); mol != NULL;
627 >           mol = info_->nextMolecule(i)) {
628 >        RealType mass = mol->getMass();
629 >        totalMass += mass;
630 >        com += mass * mol->getCom();
631 >      }  
632 >      
633 > #ifdef IS_MPI
634 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
635 >                    MPI_SUM, MPI_COMM_WORLD);
636 >      MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3,
637 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
638 > #endif
639 >      
640 >      com /= totalMass;
641 >      snap->setCOM(com);
642 >    }
643 >    return snap->getCOM();
644 >  }        
645  
646 <        std::pair<int, int> tap = simParams->getTaggedAtomPair();
647 <        Vector3d pos1, pos2, rab;
646 >  /**
647 >   * Returns center of mass and center of mass velocity in one
648 >   * function call.
649 >   */  
650 >  void Thermo::getComAll(Vector3d &com, Vector3d &comVel){
651 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
652  
653 +    if (!(snap->hasCOM && snap->hasCOMvel)) {
654 +
655 +      SimInfo::MoleculeIterator i;
656 +      Molecule* mol;
657 +      
658 +      RealType totalMass(0.0);
659 +      
660 +      com = 0.0;
661 +      comVel = 0.0;
662 +      
663 +      for (mol = info_->beginMolecule(i); mol != NULL;
664 +           mol = info_->nextMolecule(i)) {
665 +        RealType mass = mol->getMass();
666 +        totalMass += mass;
667 +        com += mass * mol->getCom();
668 +        comVel += mass * mol->getComVel();          
669 +      }  
670 +      
671   #ifdef IS_MPI
672 <        std::cerr << "tap = " << tap.first << "  " << tap.second << std::endl;
672 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
673 >                    MPI_SUM, MPI_COMM_WORLD);
674 >      MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3,
675 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
676 >      MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3,
677 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
678 > #endif
679 >      
680 >      com /= totalMass;
681 >      comVel /= totalMass;
682 >      snap->setCOM(com);
683 >      snap->setCOMvel(comVel);
684 >    }    
685 >    com = snap->getCOM();
686 >    comVel = snap->getCOMvel();
687 >    return;
688 >  }        
689 >  
690 >  /**
691 >   * \brief Return inertia tensor for entire system and angular momentum
692 >   *  Vector.
693 >   *
694 >   *
695 >   *
696 >   *    [  Ixx -Ixy  -Ixz ]
697 >   * I =| -Iyx  Iyy  -Iyz |
698 >   *    [ -Izx -Iyz   Izz ]
699 >   */
700 >  void Thermo::getInertiaTensor(Mat3x3d &inertiaTensor,
701 >                                Vector3d &angularMomentum){
702  
703 <  int mol1 = info_->getGlobalMolMembership(tap.first);
704 <  int mol2 = info_->getGlobalMolMembership(tap.second);
705 <        std::cerr << "mols = " << mol1 << " " << mol2 << std::endl;
703 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
704 >    
705 >    if (!(snap->hasInertiaTensor && snap->hasCOMw)) {
706 >      
707 >      RealType xx = 0.0;
708 >      RealType yy = 0.0;
709 >      RealType zz = 0.0;
710 >      RealType xy = 0.0;
711 >      RealType xz = 0.0;
712 >      RealType yz = 0.0;
713 >      Vector3d com(0.0);
714 >      Vector3d comVel(0.0);
715 >      
716 >      getComAll(com, comVel);
717 >      
718 >      SimInfo::MoleculeIterator i;
719 >      Molecule* mol;
720 >      
721 >      Vector3d thisq(0.0);
722 >      Vector3d thisv(0.0);
723 >      
724 >      RealType thisMass = 0.0;
725 >      
726 >      for (mol = info_->beginMolecule(i); mol != NULL;
727 >           mol = info_->nextMolecule(i)) {
728 >        
729 >        thisq = mol->getCom()-com;
730 >        thisv = mol->getComVel()-comVel;
731 >        thisMass = mol->getMass();
732 >        // Compute moment of intertia coefficients.
733 >        xx += thisq[0]*thisq[0]*thisMass;
734 >        yy += thisq[1]*thisq[1]*thisMass;
735 >        zz += thisq[2]*thisq[2]*thisMass;
736 >        
737 >        // compute products of intertia
738 >        xy += thisq[0]*thisq[1]*thisMass;
739 >        xz += thisq[0]*thisq[2]*thisMass;
740 >        yz += thisq[1]*thisq[2]*thisMass;
741 >        
742 >        angularMomentum += cross( thisq, thisv ) * thisMass;            
743 >      }
744 >      
745 >      inertiaTensor(0,0) = yy + zz;
746 >      inertiaTensor(0,1) = -xy;
747 >      inertiaTensor(0,2) = -xz;
748 >      inertiaTensor(1,0) = -xy;
749 >      inertiaTensor(1,1) = xx + zz;
750 >      inertiaTensor(1,2) = -yz;
751 >      inertiaTensor(2,0) = -xz;
752 >      inertiaTensor(2,1) = -yz;
753 >      inertiaTensor(2,2) = xx + yy;
754 >      
755 > #ifdef IS_MPI
756 >      MPI_Allreduce(MPI_IN_PLACE, inertiaTensor.getArrayPointer(),
757 >                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
758 >      MPI_Allreduce(MPI_IN_PLACE,
759 >                    angularMomentum.getArrayPointer(), 3,
760 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
761 > #endif
762 >      
763 >      snap->setCOMw(angularMomentum);
764 >      snap->setInertiaTensor(inertiaTensor);
765 >    }
766 >    
767 >    angularMomentum = snap->getCOMw();
768 >    inertiaTensor = snap->getInertiaTensor();
769 >    
770 >    return;
771 >  }
772  
773 +
774 +  Mat3x3d Thermo::getBoundingBox(){
775 +    
776 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
777 +    
778 +    if (!(snap->hasBoundingBox)) {
779 +      
780 +      SimInfo::MoleculeIterator i;
781 +      Molecule::RigidBodyIterator ri;
782 +      Molecule::AtomIterator ai;
783 +      Molecule* mol;
784 +      RigidBody* rb;
785 +      Atom* atom;
786 +      Vector3d pos, bMax, bMin;
787 +      int index = 0;
788 +      
789 +      for (mol = info_->beginMolecule(i); mol != NULL;
790 +           mol = info_->nextMolecule(i)) {
791 +        
792 +        //change the positions of atoms which belong to the rigidbodies
793 +        for (rb = mol->beginRigidBody(ri); rb != NULL;
794 +             rb = mol->nextRigidBody(ri)) {          
795 +          rb->updateAtoms();
796 +        }
797 +        
798 +        for(atom = mol->beginAtom(ai); atom != NULL;
799 +            atom = mol->nextAtom(ai)) {
800 +          
801 +          pos = atom->getPos();
802 +
803 +          if (index == 0) {
804 +            bMax = pos;
805 +            bMin = pos;
806 +          } else {
807 +            for (int i = 0; i < 3; i++) {
808 +              bMax[i] = max(bMax[i], pos[i]);
809 +              bMin[i] = min(bMin[i], pos[i]);
810 +            }
811 +          }
812 +          index++;
813 +        }
814 +      }
815 +      
816 + #ifdef IS_MPI
817 +      MPI_Allreduce(MPI_IN_PLACE, &bMax[0], 3, MPI_REALTYPE,
818 +                    MPI_MAX, MPI_COMM_WORLD);
819 +
820 +      MPI_Allreduce(MPI_IN_PLACE, &bMin[0], 3, MPI_REALTYPE,
821 +                    MPI_MIN, MPI_COMM_WORLD);
822 + #endif
823 +      Mat3x3d bBox = Mat3x3d(0.0);
824 +      for (int i = 0; i < 3; i++) {          
825 +        bBox(i,i) = bMax[i] - bMin[i];
826 +      }
827 +      snap->setBoundingBox(bBox);
828 +    }
829 +    
830 +    return snap->getBoundingBox();    
831 +  }
832 +  
833 +  
834 +  // Returns the angular momentum of the system
835 +  Vector3d Thermo::getAngularMomentum(){
836 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
837 +    
838 +    if (!snap->hasCOMw) {
839 +      
840 +      Vector3d com(0.0);
841 +      Vector3d comVel(0.0);
842 +      Vector3d angularMomentum(0.0);
843 +      
844 +      getComAll(com, comVel);
845 +      
846 +      SimInfo::MoleculeIterator i;
847 +      Molecule* mol;
848 +      
849 +      Vector3d thisr(0.0);
850 +      Vector3d thisp(0.0);
851 +      
852 +      RealType thisMass;
853 +      
854 +      for (mol = info_->beginMolecule(i); mol != NULL;
855 +           mol = info_->nextMolecule(i)) {
856 +        thisMass = mol->getMass();
857 +        thisr = mol->getCom() - com;
858 +        thisp = (mol->getComVel() - comVel) * thisMass;
859 +        
860 +        angularMomentum += cross( thisr, thisp );      
861 +      }  
862 +      
863 + #ifdef IS_MPI
864 +      MPI_Allreduce(MPI_IN_PLACE,
865 +                    angularMomentum.getArrayPointer(), 3,
866 +                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
867 + #endif
868 +      
869 +      snap->setCOMw(angularMomentum);
870 +    }
871 +    
872 +    return snap->getCOMw();
873 +  }
874 +  
875 +  
876 +  /**
877 +   * Returns the Volume of the system based on a ellipsoid with
878 +   * semi-axes based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3
879 +   * where R_i are related to the principle inertia moments
880 +   *  R_i = sqrt(C*I_i/N), this reduces to
881 +   *  V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)).
882 +   * See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
883 +   */
884 +  RealType Thermo::getGyrationalVolume(){
885 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
886 +    
887 +    if (!snap->hasGyrationalVolume) {
888 +      
889 +      Mat3x3d intTensor;
890 +      RealType det;
891 +      Vector3d dummyAngMom;
892 +      RealType sysconstants;
893 +      RealType geomCnst;
894 +      RealType volume;
895 +      
896 +      geomCnst = 3.0/2.0;
897 +      /* Get the inertial tensor and angular momentum for free*/
898 +      getInertiaTensor(intTensor, dummyAngMom);
899 +      
900 +      det = intTensor.determinant();
901 +      sysconstants = geomCnst / (RealType)(info_->getNGlobalIntegrableObjects());
902 +      volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(det);
903 +
904 +      snap->setGyrationalVolume(volume);
905 +    }
906 +    return snap->getGyrationalVolume();
907 +  }
908 +  
909 +  void Thermo::getGyrationalVolume(RealType &volume, RealType &detI){
910 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
911 +
912 +    if (!(snap->hasInertiaTensor && snap->hasGyrationalVolume)) {
913 +    
914 +      Mat3x3d intTensor;
915 +      Vector3d dummyAngMom;
916 +      RealType sysconstants;
917 +      RealType geomCnst;
918 +      
919 +      geomCnst = 3.0/2.0;
920 +      /* Get the inertia tensor and angular momentum for free*/
921 +      this->getInertiaTensor(intTensor, dummyAngMom);
922 +      
923 +      detI = intTensor.determinant();
924 +      sysconstants = geomCnst/(RealType)(info_->getNGlobalIntegrableObjects());
925 +      volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(detI);
926 +      snap->setGyrationalVolume(volume);
927 +    } else {
928 +      volume = snap->getGyrationalVolume();
929 +      detI = snap->getInertiaTensor().determinant();
930 +    }
931 +    return;
932 +  }
933 +  
934 +  RealType Thermo::getTaggedAtomPairDistance(){
935 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
936 +    Globals* simParams = info_->getSimParams();
937 +    
938 +    if (simParams->haveTaggedAtomPair() &&
939 +        simParams->havePrintTaggedPairDistance()) {
940 +      if ( simParams->getPrintTaggedPairDistance()) {
941 +        
942 +        pair<int, int> tap = simParams->getTaggedAtomPair();
943 +        Vector3d pos1, pos2, rab;
944 +        
945 + #ifdef IS_MPI        
946 +        int mol1 = info_->getGlobalMolMembership(tap.first);
947 +        int mol2 = info_->getGlobalMolMembership(tap.second);
948 +
949          int proc1 = info_->getMolToProc(mol1);
950          int proc2 = info_->getMolToProc(mol2);
951  
952 <        std::cerr << " procs = " << proc1 << " " <<proc2 <<std::endl;
265 <
266 <  RealType data[3];
952 >        RealType data[3];
953          if (proc1 == worldRank) {
954            StuntDouble* sd1 = info_->getIOIndexToIntegrableObject(tap.first);
269 –          std::cerr << " on proc " << proc1 << ", sd1 has global index= " << sd1->getGlobalIndex() << std::endl;
955            pos1 = sd1->getPos();
956            data[0] = pos1.x();
957            data[1] = pos1.y();
958 <          data[2] = pos1.z();
958 >          data[2] = pos1.z();          
959            MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
960          } else {
961            MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
962            pos1 = Vector3d(data);
963          }
964  
280 –
965          if (proc2 == worldRank) {
966            StuntDouble* sd2 = info_->getIOIndexToIntegrableObject(tap.second);
283 –          std::cerr << " on proc " << proc2 << ", sd2 has global index= " << sd2->getGlobalIndex() << std::endl;
967            pos2 = sd2->getPos();
968            data[0] = pos2.x();
969            data[1] = pos2.y();
970 <          data[2] = pos2.z();
970 >          data[2] = pos2.z();  
971            MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
972          } else {
973            MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
# Line 295 | Line 978 | namespace OpenMD {
978          StuntDouble* at2 = info_->getIOIndexToIntegrableObject(tap.second);
979          pos1 = at1->getPos();
980          pos2 = at2->getPos();
981 < #endif
981 > #endif        
982          rab = pos2 - pos1;
983          currSnapshot->wrapVector(rab);
984 <        stat[Stats::TAGGED_PAIR_DISTANCE] =  rab.length();
984 >        return rab.length();
985        }
986 +      return 0.0;    
987      }
988 <
305 <    /**@todo need refactorying*/
306 <    //Conserved Quantity is set by integrator and time is set by setTime
307 <
988 >    return 0.0;
989    }
990  
991 <
992 <
993 < Vector3d Thermo::getBoxDipole() {
994 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
995 <    SimInfo::MoleculeIterator miter;
996 <    std::vector<Atom*>::iterator aiter;
997 <    Molecule* mol;
998 <    Atom* atom;
999 <    RealType charge;
1000 <    RealType moment(0.0);
1001 <    Vector3d ri(0.0);
1002 <    Vector3d dipoleVector(0.0);
1003 <    Vector3d nPos(0.0);
1004 <    Vector3d pPos(0.0);
1005 <    RealType nChg(0.0);
325 <    RealType pChg(0.0);
326 <    int nCount = 0;
327 <    int pCount = 0;
328 <
329 <    RealType chargeToC = 1.60217733e-19;
330 <    RealType angstromToM = 1.0e-10;    RealType debyeToCm = 3.33564095198e-30;
331 <
332 <    for (mol = info_->beginMolecule(miter); mol != NULL;
333 <         mol = info_->nextMolecule(miter)) {
334 <
335 <      for (atom = mol->beginAtom(aiter); atom != NULL;
336 <           atom = mol->nextAtom(aiter)) {
337 <
338 <        if (atom->isCharge() ) {
339 <          charge = 0.0;
340 <          GenericData* data = atom->getAtomType()->getPropertyByName("Charge");
341 <          if (data != NULL) {
342 <
343 <            charge = (dynamic_cast<DoubleGenericData*>(data))->getData();
344 <            charge *= chargeToC;
345 <
346 <            ri = atom->getPos();
347 <            currSnapshot->wrapVector(ri);
348 <            ri *= angstromToM;
349 <
350 <            if (charge < 0.0) {
351 <              nPos += ri;
352 <              nChg -= charge;
353 <              nCount++;
354 <            } else if (charge > 0.0) {
355 <              pPos += ri;
356 <              pChg += charge;
357 <              pCount++;
358 <            }
359 <          }
360 <        }
361 <
362 <        if (atom->isDipole() ) {
363 <          Vector3d u_i = atom->getElectroFrame().getColumn(2);
364 <          GenericData* data = dynamic_cast<DirectionalAtomType*>(atom->getAtomType())->getPropertyByName("Dipole");
365 <          if (data != NULL) {
366 <            moment = (dynamic_cast<DoubleGenericData*>(data))->getData();
367 <
368 <            moment *= debyeToCm;
369 <            dipoleVector += u_i * moment;
370 <          }
371 <        }
991 >  RealType Thermo::getHullVolume(){
992 > #ifdef HAVE_QHULL    
993 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
994 >    if (!snap->hasHullVolume) {
995 >      Hull* surfaceMesh_;
996 >      
997 >      Globals* simParams = info_->getSimParams();
998 >      const std::string ht = simParams->getHULL_Method();
999 >      
1000 >      if (ht == "Convex") {
1001 >        surfaceMesh_ = new ConvexHull();
1002 >      } else if (ht == "AlphaShape") {
1003 >        surfaceMesh_ = new AlphaHull(simParams->getAlpha());
1004 >      } else {
1005 >        return 0.0;
1006        }
1007 +      
1008 +      // Build a vector of stunt doubles to determine if they are
1009 +      // surface atoms
1010 +      std::vector<StuntDouble*> localSites_;
1011 +      Molecule* mol;
1012 +      StuntDouble* sd;
1013 +      SimInfo::MoleculeIterator i;
1014 +      Molecule::IntegrableObjectIterator  j;
1015 +      
1016 +      for (mol = info_->beginMolecule(i); mol != NULL;
1017 +           mol = info_->nextMolecule(i)) {          
1018 +        for (sd = mol->beginIntegrableObject(j);
1019 +             sd != NULL;
1020 +             sd = mol->nextIntegrableObject(j)) {  
1021 +          localSites_.push_back(sd);
1022 +        }
1023 +      }  
1024 +      
1025 +      // Compute surface Mesh
1026 +      surfaceMesh_->computeHull(localSites_);
1027 +      snap->setHullVolume(surfaceMesh_->getVolume());
1028 +      
1029 +      delete surfaceMesh_;
1030      }
1031 <
1032 <
376 < #ifdef IS_MPI
377 <    RealType pChg_global, nChg_global;
378 <    int pCount_global, nCount_global;
379 <    Vector3d pPos_global, nPos_global, dipVec_global;
380 <
381 <    MPI_Allreduce(&pChg, &pChg_global, 1, MPI_REALTYPE, MPI_SUM,
382 <                  MPI_COMM_WORLD);
383 <    pChg = pChg_global;
384 <    MPI_Allreduce(&nChg, &nChg_global, 1, MPI_REALTYPE, MPI_SUM,
385 <                  MPI_COMM_WORLD);
386 <    nChg = nChg_global;
387 <    MPI_Allreduce(&pCount, &pCount_global, 1, MPI_INTEGER, MPI_SUM,
388 <                  MPI_COMM_WORLD);
389 <    pCount = pCount_global;
390 <    MPI_Allreduce(&nCount, &nCount_global, 1, MPI_INTEGER, MPI_SUM,
391 <                  MPI_COMM_WORLD);
392 <    nCount = nCount_global;
393 <    MPI_Allreduce(pPos.getArrayPointer(), pPos_global.getArrayPointer(), 3,
394 <                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
395 <    pPos = pPos_global;
396 <    MPI_Allreduce(nPos.getArrayPointer(), nPos_global.getArrayPointer(), 3,
397 <                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
398 <    nPos = nPos_global;
399 <    MPI_Allreduce(dipoleVector.getArrayPointer(),
400 <                  dipVec_global.getArrayPointer(), 3,
401 <                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
402 <    dipoleVector = dipVec_global;
403 < #endif //is_mpi
404 <
405 <    // first load the accumulated dipole moment (if dipoles were present)
406 <    Vector3d boxDipole = dipoleVector;
407 <    // now include the dipole moment due to charges
408 <    // use the lesser of the positive and negative charge totals
409 <    RealType chg_value = nChg <= pChg ? nChg : pChg;
410 <
411 <    // find the average positions
412 <    if (pCount > 0 && nCount > 0 ) {
413 <      pPos /= pCount;
414 <      nPos /= nCount;
415 <    }
416 <
417 <    // dipole is from the negative to the positive (physics notation)
418 <    boxDipole += (pPos - nPos) * chg_value;
419 <
420 <    return boxDipole;
421 <  }
422 <
423 < Vector3d Thermo::getThermalHelfand() {
424 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
425 <    SimInfo::MoleculeIterator miter;
426 <    std::vector<Atom*>::iterator aiter;
427 <    Molecule* mol;
428 <    Atom* atom;
429 <    RealType mass;
430 <    Vector3d velocity;
431 <    Vector3d x_a;
432 <    RealType kinetic;
433 <    RealType potential;
434 <    RealType eatom;
435 <    RealType AvgE_a_ = 0;
436 <    Vector3d GKappa_t = V3Zero;
437 <    Vector3d ThermalHelfandMoment;
438 <
439 <    for (mol = info_->beginMolecule(miter); mol != NULL;
440 <         mol = info_->nextMolecule(miter)) {
441 <
442 <      for (atom = mol->beginAtom(aiter); atom != NULL;
443 <           atom = mol->nextAtom(aiter)) {
444 <
445 <        mass = atom->getMass();
446 <        velocity = atom->getVel();
447 <        kinetic = mass * (velocity[0]*velocity[0] + velocity[1]*velocity[1] +
448 <                                   velocity[2]*velocity[2]) / PhysicalConstants::energyConvert;
449 <        potential =  atom->getParticlePot();
450 <        eatom += (kinetic + potential)/2.0;
451 <      }
452 <    }
453 <
454 <   int natoms = info_->getNGlobalAtoms();
455 < #ifdef IS_MPI
456 <
457 <    MPI_Allreduce(&eatom, &AvgE_a_, 1, MPI_REALTYPE, MPI_SUM,
458 <                  MPI_COMM_WORLD);
1031 >    
1032 >    return snap->getHullVolume();
1033   #else
1034 <    AvgE_a_ = eatom;
1034 >    return 0.0;
1035   #endif
1036 <    AvgE_a_ = AvgE_a_/RealType(natoms);
1036 >  }
1037  
464 –    for (mol = info_->beginMolecule(miter); mol != NULL;
465 –         mol = info_->nextMolecule(miter)) {
1038  
1039 <      for (atom = mol->beginAtom(aiter); atom != NULL;
468 <           atom = mol->nextAtom(aiter)) {
469 <
470 <        /* We think that x_a is relative to the total box and should be a wrapped coordinate */
471 <        x_a = atom->getPos();
472 <        currSnapshot->wrapVector(x_a);
473 <        potential =  atom->getParticlePot();
474 <        velocity = atom->getVel();
475 <        kinetic = mass * (velocity[0]*velocity[0] + velocity[1]*velocity[1] +
476 <                           velocity[2]*velocity[2]) / PhysicalConstants::energyConvert;
477 <        eatom += (kinetic + potential)/2.0
478 <        GKappa_t += x_a*(eatom-AvgE_a_);
479 <        }
480 <      }
481 < #ifdef IS_MPI
482 <     MPI_Allreduce(GKappa_t.getArrayPointer(), ThermalHelfandMoment.getArrayPointer(), 3,
483 <                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
484 < #else
485 <     ThermalHelfandMoment = GKappa_t;
486 < #endif
487 <     return ThermalHelfandMoment;
488 <
489 < }
490 <
491 <
492 <
493 < } //end namespace OpenMD
1039 > }

Diff Legend

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