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

Comparing branches/development/src/brains/Thermo.cpp (file contents):
Revision 1503 by gezelter, Sat Oct 2 19:54:41 2010 UTC vs.
Revision 1764 by gezelter, Tue Jul 3 18:32:27 2012 UTC

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

Diff Legend

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