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 1760 by gezelter, Thu Jun 21 19:26:46 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;
165  
166 <    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
167 <    potential = curSnapshot->getShortRangePotential() + curSnapshot->getLongRangePotential();
168 <    return potential;
166 >    // ForceManager computes the potential and stores it in the
167 >    // Snapshot.  All we have to do is report it.
168 >
169 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
170 >    return snap->getPotentialEnergy();
171    }
172  
173 <  RealType Thermo::getTotalE() {
122 <    RealType total;
173 >  RealType Thermo::getTotalEnergy() {
174  
175 <    total = this->getKinetic() + this->getPotential();
176 <    return total;
175 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
176 >
177 >    if (!snap->hasTotalEnergy) {
178 >      snap->setTotalEnergy(this->getKinetic() + this->getPotential());
179 >    }
180 >
181 >    return snap->getTotalEnergy();
182    }
183  
184    RealType Thermo::getTemperature() {
185 +
186 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
187 +
188 +    if (!snap->hasTemperature) {
189 +
190 +      RealType temperature = ( 2.0 * this->getKinetic() )
191 +        / (info_->getNdf()* PhysicalConstants::kb );
192 +
193 +      snap->setTemperature(temperature);
194 +    }
195      
196 <    RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* PhysicalConstants::kb );
131 <    return temperature;
196 >    return snap->getTemperature();
197    }
198  
199    RealType Thermo::getElectronicTemperature() {
200 <    SimInfo::MoleculeIterator miter;
201 <    std::vector<Atom*>::iterator iiter;
202 <    Molecule* mol;
203 <    Atom* atom;    
204 <    RealType cvel;
205 <    RealType cmass;
206 <    RealType kinetic = 0.0;
207 <    RealType kinetic_global = 0.0;
208 <    
209 <    for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
210 <      for (atom = mol->beginFluctuatingCharge(iiter); atom != NULL;
211 <           atom = mol->nextFluctuatingCharge(iiter)) {
212 <        cmass = atom->getChargeMass();
213 <        cvel = atom->getFlucQVel();
200 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
201 >
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 <        kinetic += cmass * cvel * cvel;
217 <        
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        }
153 –    }
226      
227   #ifdef IS_MPI
228 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE,
229 +                                MPI::SUM);
230 + #endif
231  
232 <    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM,
233 <                  MPI_COMM_WORLD);
234 <    kinetic = kinetic_global;
232 >      kinetic *= 0.5;
233 >      eTemp =  (2.0 * kinetic) /
234 >        (info_->getNFluctuatingCharges() * PhysicalConstants::kb );
235 >    
236 >      snap->setElectronicTemperature(eTemp);
237 >    }
238  
239 < #endif //is_mpi
162 <
163 <    kinetic = kinetic * 0.5;
164 <    return ( 2.0 * kinetic) / (info_->getNFluctuatingCharges()* PhysicalConstants::kb );    
239 >    return snap->getElectronicTemperature();
240    }
241  
242  
168 –
169 –
243    RealType Thermo::getVolume() {
244 <    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
245 <    return curSnapshot->getVolume();
244 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
245 >    return snap->getVolume();
246    }
247  
248    RealType Thermo::getPressure() {
249 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
250  
251 <    // Relies on the calculation of the full molecular pressure tensor
252 <
253 <
254 <    Mat3x3d tensor;
255 <    RealType pressure;
256 <
257 <    tensor = getPressureTensor();
258 <
259 <    pressure = PhysicalConstants::pressureConvert * (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0;
260 <
261 <    return pressure;
262 <  }
263 <
264 <  RealType Thermo::getPressure(int direction) {
265 <
192 <    // Relies on the calculation of the full molecular pressure tensor
193 <
194 <          
195 <    Mat3x3d tensor;
196 <    RealType pressure;
197 <
198 <    tensor = getPressureTensor();
199 <
200 <    pressure = PhysicalConstants::pressureConvert * tensor(direction, direction);
201 <
202 <    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;
210 <    Mat3x3d p_local(0.0);
211 <    Mat3x3d p_global(0.0);
272 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
273  
274 <    SimInfo::MoleculeIterator i;
214 <    std::vector<StuntDouble*>::iterator j;
215 <    Molecule* mol;
216 <    StuntDouble* integrableObject;    
217 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
218 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
219 <           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 <    }
226 <    
296 >      
297   #ifdef IS_MPI
298 <    MPI_Allreduce(p_local.getArrayPointer(), p_global.getArrayPointer(), 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
299 < #else
300 <    p_global = p_local;
301 < #endif // is_mpi
302 <
303 <    RealType volume = this->getVolume();
304 <    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
305 <    Mat3x3d stressTensor = curSnapshot->getStressTensor();
306 <
307 <    pressureTensor =  (p_global +
308 <                       PhysicalConstants::energyConvert * stressTensor)/volume;
309 <    
310 <    return pressureTensor;
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 >    return snap->getPressureTensor();
311    }
312  
313  
244 –  void Thermo::saveStat(){
245 –    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
246 –    Stats& stat = currSnapshot->statData;
247 –    
248 –    stat[Stats::KINETIC_ENERGY] = getKinetic();
249 –    stat[Stats::POTENTIAL_ENERGY] = getPotential();
250 –    stat[Stats::TOTAL_ENERGY] = stat[Stats::KINETIC_ENERGY]  + stat[Stats::POTENTIAL_ENERGY] ;
251 –    stat[Stats::TEMPERATURE] = getTemperature();
252 –    stat[Stats::PRESSURE] = getPressure();
253 –    stat[Stats::VOLUME] = getVolume();      
314  
255 –    Mat3x3d tensor =getPressureTensor();
256 –    stat[Stats::PRESSURE_TENSOR_XX] = tensor(0, 0);      
257 –    stat[Stats::PRESSURE_TENSOR_XY] = tensor(0, 1);      
258 –    stat[Stats::PRESSURE_TENSOR_XZ] = tensor(0, 2);      
259 –    stat[Stats::PRESSURE_TENSOR_YX] = tensor(1, 0);      
260 –    stat[Stats::PRESSURE_TENSOR_YY] = tensor(1, 1);      
261 –    stat[Stats::PRESSURE_TENSOR_YZ] = tensor(1, 2);      
262 –    stat[Stats::PRESSURE_TENSOR_ZX] = tensor(2, 0);      
263 –    stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);      
264 –    stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);      
315  
316 <    // grab the simulation box dipole moment if specified
317 <    if (info_->getCalcBoxDipole()){
268 <      Vector3d totalDipole = getBoxDipole();
269 <      stat[Stats::BOX_DIPOLE_X] = totalDipole(0);
270 <      stat[Stats::BOX_DIPOLE_Y] = totalDipole(1);
271 <      stat[Stats::BOX_DIPOLE_Z] = totalDipole(2);
272 <    }
316 >  Vector3d Thermo::getSystemDipole() {
317 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
318  
319 <    Globals* simParams = info_->getSimParams();
320 <    // grab the heat flux if desired
321 <    if (simParams->havePrintHeatFlux()) {
322 <      if (simParams->getPrintHeatFlux()){
323 <        Vector3d heatFlux = getHeatFlux();
324 <        stat[Stats::HEATFLUX_X] = heatFlux(0);
325 <        stat[Stats::HEATFLUX_Y] = heatFlux(1);
326 <        stat[Stats::HEATFLUX_Z] = heatFlux(2);
327 <      }
328 <    }
329 <
330 <    if (simParams->haveTaggedAtomPair() &&
331 <        simParams->havePrintTaggedPairDistance()) {
332 <      if ( simParams->getPrintTaggedPairDistance()) {
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 <        std::pair<int, int> tap = simParams->getTaggedAtomPair();
343 <        Vector3d pos1, pos2, rab;
344 <
345 < #ifdef IS_MPI        
346 <        std::cerr << "tap = " << tap.first << "  " << tap.second << std::endl;
347 <
348 <        int mol1 = info_->getGlobalMolMembership(tap.first);
349 <        int mol2 = info_->getGlobalMolMembership(tap.second);
350 <        std::cerr << "mols = " << mol1 << " " << mol2 << std::endl;
351 <
352 <        int proc1 = info_->getMolToProc(mol1);
353 <        int proc2 = info_->getMolToProc(mol2);
354 <
355 <        std::cerr << " procs = " << proc1 << " " <<proc2 <<std::endl;
356 <
357 <        RealType data[3];
358 <        if (proc1 == worldRank) {
359 <          StuntDouble* sd1 = info_->getIOIndexToIntegrableObject(tap.first);
360 <          std::cerr << " on proc " << proc1 << ", sd1 has global index= " << sd1->getGlobalIndex() << std::endl;
361 <          pos1 = sd1->getPos();
362 <          data[0] = pos1.x();
363 <          data[1] = pos1.y();
364 <          data[2] = pos1.z();          
365 <          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
366 <        } else {
367 <          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
368 <          pos1 = Vector3d(data);
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          }
317 –
318 –
319 –        if (proc2 == worldRank) {
320 –          StuntDouble* sd2 = info_->getIOIndexToIntegrableObject(tap.second);
321 –          std::cerr << " on proc " << proc2 << ", sd2 has global index= " << sd2->getGlobalIndex() << std::endl;
322 –          pos2 = sd2->getPos();
323 –          data[0] = pos2.x();
324 –          data[1] = pos2.y();
325 –          data[2] = pos2.z();          
326 –          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
327 –        } else {
328 –          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
329 –          pos2 = Vector3d(data);
330 –        }
331 – #else
332 –        StuntDouble* at1 = info_->getIOIndexToIntegrableObject(tap.first);
333 –        StuntDouble* at2 = info_->getIOIndexToIntegrableObject(tap.second);
334 –        pos1 = at1->getPos();
335 –        pos2 = at2->getPos();
336 – #endif        
337 –        rab = pos2 - pos1;
338 –        currSnapshot->wrapVector(rab);
339 –        stat[Stats::TAGGED_PAIR_DISTANCE] =  rab.length();
381        }
341 –    }
382        
383 <    /**@todo need refactorying*/
384 <    //Conserved Quantity is set by integrator and time is set by setTime
385 <    
386 <  }
383 >      
384 > #ifdef 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 +      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 <  Vector3d Thermo::getBoxDipole() {
396 <    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
397 <    SimInfo::MoleculeIterator miter;
398 <    std::vector<Atom*>::iterator aiter;
353 <    Molecule* mol;
354 <    Atom* atom;
355 <    RealType charge;
356 <    RealType moment(0.0);
357 <    Vector3d ri(0.0);
358 <    Vector3d dipoleVector(0.0);
359 <    Vector3d nPos(0.0);
360 <    Vector3d pPos(0.0);
361 <    RealType nChg(0.0);
362 <    RealType pChg(0.0);
363 <    int nCount = 0;
364 <    int pCount = 0;
365 <
366 <    RealType chargeToC = 1.60217733e-19;
367 <    RealType angstromToM = 1.0e-10;
368 <    RealType debyeToCm = 3.33564095198e-30;
369 <    
370 <    for (mol = info_->beginMolecule(miter); mol != NULL;
371 <         mol = info_->nextMolecule(miter)) {
372 <
373 <      for (atom = mol->beginAtom(aiter); atom != NULL;
374 <           atom = mol->nextAtom(aiter)) {
375 <        
376 <        if (atom->isCharge() ) {
377 <          charge = 0.0;
378 <          GenericData* data = atom->getAtomType()->getPropertyByName("Charge");
379 <          if (data != NULL) {
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 <            charge = (dynamic_cast<DoubleGenericData*>(data))->getData();
401 <            charge *= chargeToC;
402 <
403 <            ri = atom->getPos();
404 <            currSnapshot->wrapVector(ri);
405 <            ri *= angstromToM;
406 <
407 <            if (charge < 0.0) {
408 <              nPos += ri;
409 <              nChg -= charge;
410 <              nCount++;
411 <            } else if (charge > 0.0) {
412 <              pPos += ri;
413 <              pChg += charge;
395 <              pCount++;
396 <            }                      
397 <          }
398 <        }
399 <        
400 <        MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
401 <        if (ma.isDipole() ) {
402 <          Vector3d u_i = atom->getElectroFrame().getColumn(2);
403 <          moment = ma.getDipoleMoment();
404 <          moment *= debyeToCm;
405 <          dipoleVector += u_i * moment;
406 <        }
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      }
409 –    
410 –                      
411 – #ifdef IS_MPI
412 –    RealType pChg_global, nChg_global;
413 –    int pCount_global, nCount_global;
414 –    Vector3d pPos_global, nPos_global, dipVec_global;
420  
421 <    MPI_Allreduce(&pChg, &pChg_global, 1, MPI_REALTYPE, MPI_SUM,
417 <                  MPI_COMM_WORLD);
418 <    pChg = pChg_global;
419 <    MPI_Allreduce(&nChg, &nChg_global, 1, MPI_REALTYPE, MPI_SUM,
420 <                  MPI_COMM_WORLD);
421 <    nChg = nChg_global;
422 <    MPI_Allreduce(&pCount, &pCount_global, 1, MPI_INTEGER, MPI_SUM,
423 <                  MPI_COMM_WORLD);
424 <    pCount = pCount_global;
425 <    MPI_Allreduce(&nCount, &nCount_global, 1, MPI_INTEGER, MPI_SUM,
426 <                  MPI_COMM_WORLD);
427 <    nCount = nCount_global;
428 <    MPI_Allreduce(pPos.getArrayPointer(), pPos_global.getArrayPointer(), 3,
429 <                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
430 <    pPos = pPos_global;
431 <    MPI_Allreduce(nPos.getArrayPointer(), nPos_global.getArrayPointer(), 3,
432 <                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
433 <    nPos = nPos_global;
434 <    MPI_Allreduce(dipoleVector.getArrayPointer(),
435 <                  dipVec_global.getArrayPointer(), 3,
436 <                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
437 <    dipoleVector = dipVec_global;
438 < #endif //is_mpi
439 <
440 <    // first load the accumulated dipole moment (if dipoles were present)
441 <    Vector3d boxDipole = dipoleVector;
442 <    // now include the dipole moment due to charges
443 <    // use the lesser of the positive and negative charge totals
444 <    RealType chg_value = nChg <= pChg ? nChg : pChg;
445 <      
446 <    // find the average positions
447 <    if (pCount > 0 && nCount > 0 ) {
448 <      pPos /= pCount;
449 <      nPos /= nCount;
450 <    }
451 <
452 <    // dipole is from the negative to the positive (physics notation)
453 <    boxDipole += (pPos - nPos) * chg_value;
454 <
455 <    return boxDipole;
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      SimInfo::MoleculeIterator miter;
428 <    std::vector<StuntDouble*>::iterator iiter;
428 >    vector<StuntDouble*>::iterator iiter;
429      Molecule* mol;
430 <    StuntDouble* integrableObject;    
430 >    StuntDouble* sd;    
431      RigidBody::AtomIterator ai;
432      Atom* atom;      
433      Vector3d vel;
# Line 484 | Line 450 | namespace OpenMD {
450      for (mol = info_->beginMolecule(miter); mol != NULL;
451           mol = info_->nextMolecule(miter)) {
452        
453 <      for (integrableObject = mol->beginIntegrableObject(iiter);
454 <           integrableObject != NULL;
455 <           integrableObject = mol->nextIntegrableObject(iiter)) {
453 >      for (sd = mol->beginIntegrableObject(iiter);
454 >           sd != NULL;
455 >           sd = mol->nextIntegrableObject(iiter)) {
456          
457 <        mass = integrableObject->getMass();
458 <        vel = integrableObject->getVel();
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 (integrableObject->isDirectional()) {
463 <          angMom = integrableObject->getJ();
464 <          I = integrableObject->getI();
462 >        if (sd->isDirectional()) {
463 >          angMom = sd->getJ();
464 >          I = sd->getI();
465  
466 <          if (integrableObject->isLinear()) {
467 <            i = integrableObject->linearAxis();
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) + angMom[k] * angMom[k] / I(k, k);
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) + angMom[1]*angMom[1]/I(1, 1)
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 (integrableObject->isRigidBody()) {
482 <          RigidBody* rb = dynamic_cast<RigidBody*>(integrableObject);
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 = integrableObject->getParticlePot();
521 <          cerr << "ppot = "  << potential << "\n";
488 >          potential = sd->getParticlePot();
489          }
490  
491          potential *= PhysicalConstants::energyConvert; // amu A^2/fs^2
# Line 530 | Line 497 | namespace OpenMD {
497        }
498      }
499  
500 <    std::cerr << "Heat flux heatFluxJc is: " << heatFluxJc << std::endl;
501 <
502 <    /* The J_v vector is reduced in fortan so everyone has the global
536 <     *  Jv. Jc is computed over the local atoms and must be reduced
537 <     *  among all processors.
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,
# Line 545 | Line 510 | namespace OpenMD {
510  
511      Vector3d heatFluxJv = currSnapshot->getConductiveHeatFlux() *
512        PhysicalConstants::energyConvert;
513 <    
549 <    std::cerr << "Heat flux Jc is: " << heatFluxJc << std::endl;
550 <    std::cerr << "Heat flux Jv is: " << heatFluxJv << std::endl;
551 <    
513 >        
514      // Correct for the fact the flux is 1/V (Jc + Jv)
515      return (heatFluxJv + heatFluxJc) / this->getVolume(); // amu / fs^3
516    }
517 < } //end namespace OpenMD
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 >        pair<int, int> tap = simParams->getTaggedAtomPair();
817 >        Vector3d pos1, pos2, rab;
818 >        
819 > #ifdef IS_MPI        
820 >        int mol1 = info_->getGlobalMolMembership(tap.first);
821 >        int mol2 = info_->getGlobalMolMembership(tap.second);
822 >
823 >        int proc1 = info_->getMolToProc(mol1);
824 >        int proc2 = info_->getMolToProc(mol2);
825 >
826 >        RealType data[3];
827 >        if (proc1 == worldRank) {
828 >          StuntDouble* sd1 = info_->getIOIndexToIntegrableObject(tap.first);
829 >          pos1 = sd1->getPos();
830 >          data[0] = pos1.x();
831 >          data[1] = pos1.y();
832 >          data[2] = pos1.z();          
833 >          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
834 >        } else {
835 >          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
836 >          pos1 = Vector3d(data);
837 >        }
838 >
839 >        if (proc2 == worldRank) {
840 >          StuntDouble* sd2 = info_->getIOIndexToIntegrableObject(tap.second);
841 >          pos2 = sd2->getPos();
842 >          data[0] = pos2.x();
843 >          data[1] = pos2.y();
844 >          data[2] = pos2.z();          
845 >          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
846 >        } else {
847 >          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
848 >          pos2 = Vector3d(data);
849 >        }
850 > #else
851 >        StuntDouble* at1 = info_->getIOIndexToIntegrableObject(tap.first);
852 >        StuntDouble* at2 = info_->getIOIndexToIntegrableObject(tap.second);
853 >        pos1 = at1->getPos();
854 >        pos2 = at2->getPos();
855 > #endif        
856 >        rab = pos2 - pos1;
857 >        currSnapshot->wrapVector(rab);
858 >        return rab.length();
859 >      }
860 >      return 0.0;    
861 >    }
862 >    return 0.0;
863 >  }
864 >
865 >  RealType Thermo::getHullVolume(){
866 >    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
867 >    
868 >    if (!snap->hasHullVolume) {
869 >
870 >      Hull* surfaceMesh_;
871 >
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 >    return snap->getHullVolume();
905 >  }  
906 > }

Diff Legend

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