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 1710 by gezelter, Fri May 18 21:44:02 2012 UTC vs.
Revision 1764 by gezelter, Tue Jul 3 18:32:27 2012 UTC

# Line 51 | 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;
72 <    
73 <    for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
74 <      for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL;
75 <           integrableObject = mol->nextIntegrableObject(iiter)) {
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 <    }
99 <    
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,
103 <                  MPI_COMM_WORLD);
104 <    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() {
114 –    RealType potential = 0.0;
115 –    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
116 –    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,
123 <                  MPI_COMM_WORLD);
124 <    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
131 <
132 <    return potential;
181 >    return snap->getTotalEnergy();
182    }
183  
135 –  RealType Thermo::getTotalE() {
136 –    RealType total;
137 –
138 –    total = this->getKinetic() + this->getPotential();
139 –    return total;
140 –  }
141 –
184    RealType Thermo::getTemperature() {
143 –    
144 –    RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* PhysicalConstants::kb );
145 –    return temperature;
146 –  }
185  
186 <  RealType Thermo::getVolume() {
149 <    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
150 <    return curSnapshot->getVolume();
151 <  }
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  
168 –  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;
174 <    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;
188 <    Mat3x3d p_local(0.0);
189 <    Mat3x3d p_global(0.0);
272 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
273  
274 <    SimInfo::MoleculeIterator i;
192 <    std::vector<StuntDouble*>::iterator j;
193 <    Molecule* mol;
194 <    StuntDouble* integrableObject;    
195 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
196 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
197 <           integrableObject = mol->nextIntegrableObject(j)) {
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->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);      
239 <    stat[Stats::PRESSURE_TENSOR_ZX] = tensor(2, 0);      
240 <    stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);      
241 <    stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);      
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        
261 –        std::cerr << "tap = " << tap.first << "  " << tap.second << std::endl;
262 –
820          int mol1 = info_->getGlobalMolMembership(tap.first);
821          int mol2 = info_->getGlobalMolMembership(tap.second);
265 –        std::cerr << "mols = " << mol1 << " " << mol2 << std::endl;
822  
823          int proc1 = info_->getMolToProc(mol1);
824          int proc2 = info_->getMolToProc(mol2);
825  
270 –        std::cerr << " procs = " << proc1 << " " <<proc2 <<std::endl;
271 –
826          RealType data[3];
827          if (proc1 == worldRank) {
828            StuntDouble* sd1 = info_->getIOIndexToIntegrableObject(tap.first);
275 –          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 283 | Line 836 | namespace OpenMD {
836            pos1 = Vector3d(data);
837          }
838  
286 –
839          if (proc2 == worldRank) {
840            StuntDouble* sd2 = info_->getIOIndexToIntegrableObject(tap.second);
289 –          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 304 | 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 <      
311 <    /**@todo need refactorying*/
312 <    //Conserved Quantity is set by integrator and time is set by setTime
313 <    
862 >    return 0.0;
863    }
864  
865 <
866 <  Vector3d Thermo::getBoxDipole() {
318 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
319 <    SimInfo::MoleculeIterator miter;
320 <    std::vector<Atom*>::iterator aiter;
321 <    Molecule* mol;
322 <    Atom* atom;
323 <    RealType charge;
324 <    RealType moment(0.0);
325 <    Vector3d ri(0.0);
326 <    Vector3d dipoleVector(0.0);
327 <    Vector3d nPos(0.0);
328 <    Vector3d pPos(0.0);
329 <    RealType nChg(0.0);
330 <    RealType pChg(0.0);
331 <    int nCount = 0;
332 <    int pCount = 0;
333 <
334 <    RealType chargeToC = 1.60217733e-19;
335 <    RealType angstromToM = 1.0e-10;
336 <    RealType debyeToCm = 3.33564095198e-30;
865 >  RealType Thermo::getHullVolume(){
866 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
867      
868 <    for (mol = info_->beginMolecule(miter); mol != NULL;
339 <         mol = info_->nextMolecule(miter)) {
868 >    if (!snap->hasHullVolume) {
869  
870 <      for (atom = mol->beginAtom(aiter); atom != NULL;
342 <           atom = mol->nextAtom(aiter)) {
343 <        
344 <        if (atom->isCharge() ) {
345 <          charge = 0.0;
346 <          GenericData* data = atom->getAtomType()->getPropertyByName("Charge");
347 <          if (data != NULL) {
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;
358 <              nChg -= charge;
359 <              nCount++;
360 <            } else if (charge > 0.0) {
361 <              pPos += ri;
362 <              pChg += charge;
363 <              pCount++;
364 <            }                      
365 <          }
366 <        }
367 <        
368 <        MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
369 <        if (ma.isDipole() ) {
370 <          Vector3d u_i = atom->getElectroFrame().getColumn(2);
371 <          moment = ma.getDipoleMoment();
372 <          moment *= debyeToCm;
373 <          dipoleVector += u_i * moment;
374 <        }
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
380 <    RealType pChg_global, nChg_global;
381 <    int pCount_global, nCount_global;
382 <    Vector3d pPos_global, nPos_global, dipVec_global;
383 <
384 <    MPI_Allreduce(&pChg, &pChg_global, 1, MPI_REALTYPE, MPI_SUM,
385 <                  MPI_COMM_WORLD);
386 <    pChg = pChg_global;
387 <    MPI_Allreduce(&nChg, &nChg_global, 1, MPI_REALTYPE, MPI_SUM,
388 <                  MPI_COMM_WORLD);
389 <    nChg = nChg_global;
390 <    MPI_Allreduce(&pCount, &pCount_global, 1, MPI_INTEGER, MPI_SUM,
391 <                  MPI_COMM_WORLD);
392 <    pCount = pCount_global;
393 <    MPI_Allreduce(&nCount, &nCount_global, 1, MPI_INTEGER, MPI_SUM,
394 <                  MPI_COMM_WORLD);
395 <    nCount = nCount_global;
396 <    MPI_Allreduce(pPos.getArrayPointer(), pPos_global.getArrayPointer(), 3,
397 <                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
398 <    pPos = pPos_global;
399 <    MPI_Allreduce(nPos.getArrayPointer(), nPos_global.getArrayPointer(), 3,
400 <                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
401 <    nPos = nPos_global;
402 <    MPI_Allreduce(dipoleVector.getArrayPointer(),
403 <                  dipVec_global.getArrayPointer(), 3,
404 <                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
405 <    dipoleVector = dipVec_global;
406 < #endif //is_mpi
407 <
408 <    // first load the accumulated dipole moment (if dipoles were present)
409 <    Vector3d boxDipole = dipoleVector;
410 <    // now include the dipole moment due to charges
411 <    // use the lesser of the positive and negative charge totals
412 <    RealType chg_value = nChg <= pChg ? nChg : pChg;
413 <      
414 <    // find the average positions
415 <    if (pCount > 0 && nCount > 0 ) {
416 <      pPos /= pCount;
417 <      nPos /= nCount;
418 <    }
419 <
420 <    // dipole is from the negative to the positive (physics notation)
421 <    boxDipole += (pPos - nPos) * chg_value;
422 <
423 <    return boxDipole;
424 <  }
425 < } //end namespace OpenMD
904 >    return snap->getHullVolume();
905 >  }  
906 > }

Diff Legend

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