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

Diff Legend

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