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

Diff Legend

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