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

Comparing trunk/src/brains/Thermo.cpp (file contents):
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC vs.
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 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   */
42 –
43 – #include <math.h>
44 – #include <iostream>
42  
43   #ifdef IS_MPI
44   #include <mpi.h>
45   #endif //is_mpi
46 +
47 + #include <math.h>
48 + #include <iostream>
49  
50   #include "brains/Thermo.hpp"
51   #include "primitives/Molecule.hpp"
# Line 89 | Line 89 | namespace OpenMD {
89        }
90        
91   #ifdef IS_MPI
92 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE,
93 <                                MPI::SUM);
92 >      MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE,
93 >                    MPI_SUM, MPI_COMM_WORLD);
94   #endif
95        
96        kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert;
# Line 140 | Line 140 | namespace OpenMD {
140        }
141        
142   #ifdef IS_MPI
143 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE,
144 <                                MPI::SUM);
143 >      MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE,
144 >                    MPI_SUM, MPI_COMM_WORLD);
145   #endif
146        
147        kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert;
# Line 227 | Line 227 | namespace OpenMD {
227        }
228      
229   #ifdef IS_MPI
230 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE,
231 <                                MPI::SUM);
230 >      MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE,
231 >                    MPI_SUM, MPI_COMM_WORLD);
232   #endif
233  
234        kinetic *= 0.5;
235        eTemp =  (2.0 * kinetic) /
236 <        (info_->getNFluctuatingCharges() * PhysicalConstants::kb );
236 >        (info_->getNFluctuatingCharges() * PhysicalConstants::kb );            
237      
238        snap->setElectronicTemperature(eTemp);
239      }
# Line 297 | Line 297 | namespace OpenMD {
297        }
298        
299   #ifdef IS_MPI
300 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, p_tens.getArrayPointer(), 9,
301 <                                MPI::REALTYPE, MPI::SUM);
300 >      MPI_Allreduce(MPI_IN_PLACE, p_tens.getArrayPointer(), 9,
301 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
302   #endif
303        
304        RealType volume = this->getVolume();
# Line 324 | Line 324 | namespace OpenMD {
324        Molecule* mol;
325        Atom* atom;
326        RealType charge;
327 –      RealType moment(0.0);
327        Vector3d ri(0.0);
328        Vector3d dipoleVector(0.0);
329        Vector3d nPos(0.0);
# Line 372 | Line 371 | namespace OpenMD {
371              pCount++;
372            }
373            
374 <          MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
375 <          if (ma.isDipole() ) {
377 <            Vector3d u_i = atom->getElectroFrame().getColumn(2);
378 <            moment = ma.getDipoleMoment();
379 <            moment *= debyeToCm;
380 <            dipoleVector += u_i * moment;
374 >          if (atom->isDipole()) {
375 >            dipoleVector += atom->getDipole() * debyeToCm;
376            }
377          }
378        }
379        
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);
382 >      MPI_Allreduce(MPI_IN_PLACE, &pChg, 1, MPI_REALTYPE,
383 >                    MPI_SUM, MPI_COMM_WORLD);
384 >      MPI_Allreduce(MPI_IN_PLACE, &nChg, 1, MPI_REALTYPE,
385 >                    MPI_SUM, MPI_COMM_WORLD);
386 >      
387 >      MPI_Allreduce(MPI_IN_PLACE, &pCount, 1, MPI_INTEGER,
388 >                    MPI_SUM, MPI_COMM_WORLD);
389 >      MPI_Allreduce(MPI_IN_PLACE, &nCount, 1, MPI_INTEGER,
390 >                    MPI_SUM, MPI_COMM_WORLD);
391 >      
392 >      MPI_Allreduce(MPI_IN_PLACE, pPos.getArrayPointer(), 3,
393 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
394 >      MPI_Allreduce(MPI_IN_PLACE, nPos.getArrayPointer(), 3,
395 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
396  
397 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, pPos.getArrayPointer(), 3,
398 <                                MPI::REALTYPE, MPI::SUM);
399 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, nPos.getArrayPointer(), 3,
400 <                                MPI::REALTYPE, MPI::SUM);
401 <
402 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, dipoleVector.getArrayPointer(),
403 <                                3, MPI::REALTYPE, MPI::SUM);
397 >      MPI_Allreduce(MPI_IN_PLACE, dipoleVector.getArrayPointer(),
398 >                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
399   #endif
400        
401        // first load the accumulated dipole moment (if dipoles were present)
# Line 444 | Line 439 | namespace OpenMD {
439      RealType kinetic;
440      RealType potential;
441      RealType eatom;
447 –    RealType AvgE_a_ = 0;
442      // Convective portion of the heat flux
443      Vector3d heatFluxJc = V3Zero;
444  
# Line 504 | Line 498 | namespace OpenMD {
498       *  reduced among all processors.
499       */
500   #ifdef IS_MPI
501 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &heatFluxJc[0], 3, MPI::REALTYPE,
502 <                              MPI::SUM);
501 >    MPI_Allreduce(MPI_IN_PLACE, &heatFluxJc[0], 3, MPI_REALTYPE,
502 >                  MPI_SUM, MPI_COMM_WORLD);
503   #endif
504      
505      // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3
# Line 537 | Line 531 | namespace OpenMD {
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);
534 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
535 >                    MPI_SUM, MPI_COMM_WORLD);
536 >      MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3,
537 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
538   #endif
539        
540        comVel /= totalMass;
# Line 568 | Line 562 | namespace OpenMD {
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);
565 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
566 >                    MPI_SUM, MPI_COMM_WORLD);
567 >      MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3,
568 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
569   #endif
570        
571        com /= totalMass;
# Line 606 | Line 600 | namespace OpenMD {
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);
603 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
604 >                    MPI_SUM, MPI_COMM_WORLD);
605 >      MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3,
606 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
607 >      MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3,
608 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
609   #endif
610        
611        com /= totalMass;
# Line 625 | 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 690 | Line 684 | namespace OpenMD {
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);
687 >      MPI_Allreduce(MPI_IN_PLACE, inertiaTensor.getArrayPointer(),
688 >                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
689 >      MPI_Allreduce(MPI_IN_PLACE,
690 >                    angularMomentum.getArrayPointer(), 3,
691 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
692   #endif
693        
694        snap->setCOMw(angularMomentum);
# Line 706 | Line 700 | namespace OpenMD {
700      
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_Allreduce(MPI_IN_PLACE, &bMax[0], 3, MPI_REALTYPE,
749 +                    MPI_MAX, MPI_COMM_WORLD);
750  
751 +      MPI_Allreduce(MPI_IN_PLACE, &bMin[0], 3, MPI_REALTYPE,
752 +                    MPI_MIN, MPI_COMM_WORLD);
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 737 | Line 792 | namespace OpenMD {
792        }  
793        
794   #ifdef IS_MPI
795 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE,
796 <                                angularMomentum.getArrayPointer(), 3,
797 <                                MPI::REALTYPE, MPI::SUM);
795 >      MPI_Allreduce(MPI_IN_PLACE,
796 >                    angularMomentum.getArrayPointer(), 3,
797 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
798   #endif
799        
800        snap->setCOMw(angularMomentum);
# Line 843 | Line 898 | namespace OpenMD {
898            pos2 = sd2->getPos();
899            data[0] = pos2.x();
900            data[1] = pos2.y();
901 <          data[2] = pos2.z();          
901 >          data[2] = pos2.z();  
902            MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
903          } else {
904            MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
# Line 865 | Line 920 | namespace OpenMD {
920    }
921  
922    RealType Thermo::getHullVolume(){
868 –    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
869 –
923   #ifdef HAVE_QHULL    
924 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
925      if (!snap->hasHullVolume) {
926        Hull* surfaceMesh_;
927 <
927 >      
928        Globals* simParams = info_->getSimParams();
929        const std::string ht = simParams->getHULL_Method();
930        
# Line 902 | Line 956 | namespace OpenMD {
956        // Compute surface Mesh
957        surfaceMesh_->computeHull(localSites_);
958        snap->setHullVolume(surfaceMesh_->getVolume());
959 +      
960 +      delete surfaceMesh_;
961      }
962 +    
963      return snap->getHullVolume();
964   #else
965      return 0.0;
966   #endif
967    }
968 +
969 +
970   }

Diff Legend

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