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 1764 by gezelter, Tue Jul 3 18:32:27 2012 UTC vs.
Revision 1867 by gezelter, Mon Apr 29 17:53:48 2013 UTC

# Line 35 | Line 35
35   *                                                                      
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).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# Line 54 | Line 54
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 {
# Line 322 | Line 324 | namespace OpenMD {
324        Molecule* mol;
325        Atom* atom;
326        RealType charge;
325 –      RealType moment(0.0);
327        Vector3d ri(0.0);
328        Vector3d dipoleVector(0.0);
329        Vector3d nPos(0.0);
# Line 370 | Line 371 | namespace OpenMD {
371              pCount++;
372            }
373            
374 <          MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
375 <          if (ma.isDipole() ) {
375 <            Vector3d u_i = atom->getElectroFrame().getColumn(2);
376 <            moment = ma.getDipoleMoment();
377 <            moment *= debyeToCm;
378 <            dipoleVector += u_i * moment;
374 >          if (atom->isDipole()) {
375 >            dipoleVector += atom->getDipole() * debyeToCm;
376            }
377          }
378        }
# Line 442 | Line 439 | namespace OpenMD {
439      RealType kinetic;
440      RealType potential;
441      RealType eatom;
445 –    RealType AvgE_a_ = 0;
442      // Convective portion of the heat flux
443      Vector3d heatFluxJc = V3Zero;
444  
# Line 623 | Line 619 | namespace OpenMD {
619    }        
620    
621    /**
622 <   * Return intertia tensor for entire system and angular momentum
623 <   * Vector.
622 >   * \brief Return inertia tensor for entire system and angular momentum
623 >   *  Vector.
624     *
625     *
626     *
# Line 705 | Line 701 | namespace OpenMD {
701      return;
702    }
703  
704 +
705 +  Mat3x3d Thermo::getBoundingBox(){
706 +    
707 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
708 +    
709 +    if (!(snap->hasBoundingBox)) {
710 +      
711 +      SimInfo::MoleculeIterator i;
712 +      Molecule::RigidBodyIterator ri;
713 +      Molecule::AtomIterator ai;
714 +      Molecule* mol;
715 +      RigidBody* rb;
716 +      Atom* atom;
717 +      Vector3d pos, bMax, bMin;
718 +      int index = 0;
719 +      
720 +      for (mol = info_->beginMolecule(i); mol != NULL;
721 +           mol = info_->nextMolecule(i)) {
722 +        
723 +        //change the positions of atoms which belong to the rigidbodies
724 +        for (rb = mol->beginRigidBody(ri); rb != NULL;
725 +             rb = mol->nextRigidBody(ri)) {          
726 +          rb->updateAtoms();
727 +        }
728 +        
729 +        for(atom = mol->beginAtom(ai); atom != NULL;
730 +            atom = mol->nextAtom(ai)) {
731 +          
732 +          pos = atom->getPos();
733 +
734 +          if (index == 0) {
735 +            bMax = pos;
736 +            bMin = pos;
737 +          } else {
738 +            for (int i = 0; i < 3; i++) {
739 +              bMax[i] = max(bMax[i], pos[i]);
740 +              bMin[i] = min(bMin[i], pos[i]);
741 +            }
742 +          }
743 +          index++;
744 +        }
745 +      }
746 +      
747 + #ifdef IS_MPI
748 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bMax[0], 3, MPI::REALTYPE,
749 +                                MPI::MAX);
750 +
751 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bMin[0], 3, MPI::REALTYPE,
752 +                                MPI::MIN);
753 + #endif
754 +      Mat3x3d bBox = Mat3x3d(0.0);
755 +      for (int i = 0; i < 3; i++) {          
756 +        bBox(i,i) = bMax[i] - bMin[i];
757 +      }
758 +      snap->setBoundingBox(bBox);
759 +    }
760 +    
761 +    return snap->getBoundingBox();    
762 +  }
763 +  
764 +  
765    // Returns the angular momentum of the system
766    Vector3d Thermo::getAngularMomentum(){
767      Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 830 | Line 887 | namespace OpenMD {
887            data[0] = pos1.x();
888            data[1] = pos1.y();
889            data[2] = pos1.z();          
890 <          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
890 >          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc1);
891          } else {
892 <          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
892 >          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc1);
893            pos1 = Vector3d(data);
894          }
895  
# Line 841 | Line 898 | namespace OpenMD {
898            pos2 = sd2->getPos();
899            data[0] = pos2.x();
900            data[1] = pos2.y();
901 <          data[2] = pos2.z();          
902 <          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
901 >          data[2] = pos2.z();  
902 >          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc2);
903          } else {
904 <          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
904 >          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc2);
905            pos2 = Vector3d(data);
906          }
907   #else
# Line 864 | Line 921 | namespace OpenMD {
921  
922    RealType Thermo::getHullVolume(){
923      Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
867 –    
868 –    if (!snap->hasHullVolume) {
924  
925 + #ifdef HAVE_QHULL    
926 +    if (!snap->hasHullVolume) {
927        Hull* surfaceMesh_;
928  
929        Globals* simParams = info_->getSimParams();
# Line 900 | Line 957 | namespace OpenMD {
957        // Compute surface Mesh
958        surfaceMesh_->computeHull(localSites_);
959        snap->setHullVolume(surfaceMesh_->getVolume());
960 +
961 +      delete surfaceMesh_;
962      }
963 +
964      return snap->getHullVolume();
965 <  }  
965 > #else
966 >    return 0.0;
967 > #endif
968 >  }
969 >
970 >
971   }

Diff Legend

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