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

Comparing:
trunk/src/brains/Thermo.cpp (file contents), Revision 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
branches/development/src/brains/Thermo.cpp (file contents), Revision 1710 by gezelter, Fri May 18 21:44:02 2012 UTC

# Line 36 | Line 36
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include <math.h>
# Line 50 | Line 51
51   #include "primitives/Molecule.hpp"
52   #include "utils/simError.h"
53   #include "utils/PhysicalConstants.hpp"
54 + #include "types/MultipoleAdapter.hpp"
55  
56   namespace OpenMD {
57  
# Line 208 | Line 210 | namespace OpenMD {
210  
211      RealType volume = this->getVolume();
212      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
213 <    Mat3x3d tau = curSnapshot->statData.getTau();
213 >    Mat3x3d tau = curSnapshot->getTau();
214  
215      pressureTensor =  (p_global + PhysicalConstants::energyConvert* tau)/volume;
216      
# Line 238 | Line 240 | namespace OpenMD {
240      stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);      
241      stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);      
242  
243 +    // grab the simulation box dipole moment if specified
244 +    if (info_->getCalcBoxDipole()){
245 +      Vector3d totalDipole = getBoxDipole();
246 +      stat[Stats::BOX_DIPOLE_X] = totalDipole(0);
247 +      stat[Stats::BOX_DIPOLE_Y] = totalDipole(1);
248 +      stat[Stats::BOX_DIPOLE_Z] = totalDipole(2);
249 +    }
250  
251      Globals* simParams = info_->getSimParams();
252  
# Line 303 | Line 312 | namespace OpenMD {
312      //Conserved Quantity is set by integrator and time is set by setTime
313      
314    }
315 +
316 +
317 +  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;
337 +    
338 +    for (mol = info_->beginMolecule(miter); mol != NULL;
339 +         mol = info_->nextMolecule(miter)) {
340 +
341 +      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) {
348 +
349 +            charge = (dynamic_cast<DoubleGenericData*>(data))->getData();
350 +            charge *= chargeToC;
351 +
352 +            ri = atom->getPos();
353 +            currSnapshot->wrapVector(ri);
354 +            ri *= angstromToM;
355 +
356 +            if (charge < 0.0) {
357 +              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 +        }
375 +      }
376 +    }
377 +    
378 +                      
379 + #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

Diff Legend

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