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

Comparing branches/development/src/brains/ForceManager.cpp (file contents):
Revision 1715 by gezelter, Tue May 22 21:55:31 2012 UTC vs.
Revision 1749 by gezelter, Thu Jun 7 02:47:21 2012 UTC

# Line 390 | Line 390 | namespace OpenMD {
390        info_->prepareTopology();      
391  
392        doParticlePot_ = info_->getSimParams()->getOutputParticlePotential();
393 +      doHeatFlux_ = info_->getSimParams()->getPrintHeatFlux();
394 +      if (doHeatFlux_) doParticlePot_ = true;
395    
396      }
397  
# Line 471 | Line 473 | namespace OpenMD {
473      }
474      
475      // Zero out the stress tensor
476 <    tau *= 0.0;
477 <    
476 >    stressTensor *= 0.0;
477 >    // Zero out the heatFlux
478 >    fDecomp_->setHeatFlux( Vector3d(0.0) );    
479    }
480    
481    void ForceManager::shortRangeInteractions() {
# Line 595 | Line 598 | namespace OpenMD {
598    
599    void ForceManager::longRangeInteractions() {
600  
601 +
602      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
603      DataStorage* config = &(curSnapshot->atomData);
604      DataStorage* cgConfig = &(curSnapshot->cgData);
# Line 618 | Line 622 | namespace OpenMD {
622        // center of mass of the group is the same as position of the atom  
623        // if cutoff group does not exist
624        cgConfig->position = config->position;
625 +      cgConfig->velocity = config->velocity;
626      }
627  
628      fDecomp_->zeroWorkArrays();
629      fDecomp_->distributeData();
630      
631      int cg1, cg2, atom1, atom2, topoDist;
632 <    Vector3d d_grp, dag, d;
632 >    Vector3d d_grp, dag, d, gvel2, vel2;
633      RealType rgrpsq, rgrp, r2, r;
634      RealType electroMult, vdwMult;
635      RealType vij;
# Line 639 | Line 644 | namespace OpenMD {
644      RealType mf;
645      RealType lrPot;
646      RealType vpair;
647 +    RealType dVdFQ1(0.0);
648 +    RealType dVdFQ2(0.0);
649      potVec longRangePotential(0.0);
650      potVec workPot(0.0);
651      vector<int>::iterator ia, jb;
# Line 650 | Line 657 | namespace OpenMD {
657      idat.pot = &workPot;
658      sdat.pot = fDecomp_->getEmbeddingPotential();
659      idat.vpair = &vpair;
660 +    idat.dVdFQ1 = &dVdFQ1;
661 +    idat.dVdFQ2 = &dVdFQ2;
662      idat.f1 = &f1;
663      idat.sw = &sw;
664      idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
# Line 663 | Line 672 | namespace OpenMD {
672      } else {
673        loopStart = PAIR_LOOP;
674      }
666 –  
675      for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
676      
677        if (iLoop == loopStart) {
# Line 695 | Line 703 | namespace OpenMD {
703            
704            in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
705                                                       rgrp);
698 –          
706            atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
707            atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
708  
709 +          if (doHeatFlux_)
710 +            gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
711 +
712            for (ia = atomListRow.begin();
713                 ia != atomListRow.end(); ++ia) {            
714              atom1 = (*ia);
715 <            
715 >
716              for (jb = atomListColumn.begin();
717                   jb != atomListColumn.end(); ++jb) {              
718                atom2 = (*jb);
# Line 711 | Line 721 | namespace OpenMD {
721                  vpair = 0.0;
722                  workPot = 0.0;
723                  f1 = V3Zero;
724 +                dVdFQ1 = 0.0;
725 +                dVdFQ2 = 0.0;
726  
727                  fDecomp_->fillInteractionData(idat, atom1, atom2);
728 <                
728 >
729                  topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
730                  vdwMult = vdwScale_[topoDist];
731                  electroMult = electrostaticScale_[topoDist];
# Line 721 | Line 733 | namespace OpenMD {
733                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
734                    idat.d = &d_grp;
735                    idat.r2 = &rgrpsq;
736 +                  if (doHeatFlux_)
737 +                    vel2 = gvel2;
738                  } else {
739                    d = fDecomp_->getInteratomicVector(atom1, atom2);
740                    curSnapshot->wrapVector( d );
741                    r2 = d.lengthSquare();
742                    idat.d = &d;
743                    idat.r2 = &r2;
744 +                  if (doHeatFlux_)
745 +                    vel2 = fDecomp_->getAtomVelocityColumn(atom2);
746                  }
747                
748                  r = sqrt( *(idat.r2) );
# Line 739 | Line 755 | namespace OpenMD {
755                    fDecomp_->unpackInteractionData(idat, atom1, atom2);
756                    vij += vpair;
757                    fij += f1;
758 <                  tau -= outProduct( *(idat.d), f1);
758 >                  stressTensor -= outProduct( *(idat.d), f1);
759 >                  if (doHeatFlux_)
760 >                    fDecomp_->addToHeatFlux(*(idat.d) * dot(f1, vel2));
761                  }
762                }
763              }
# Line 752 | Line 770 | namespace OpenMD {
770                fij += fg;
771  
772                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
773 <                tau -= outProduct( *(idat.d), fg);
773 >                stressTensor -= outProduct( *(idat.d), fg);
774 >                if (doHeatFlux_)
775 >                  fDecomp_->addToHeatFlux(*(idat.d) * dot(fg, vel2));
776 >                
777                }
778            
779                for (ia = atomListRow.begin();
# Line 768 | Line 789 | namespace OpenMD {
789                      // find the distance between the atom
790                      // and the center of the cutoff group:
791                      dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
792 <                    tau -= outProduct(dag, fg);
792 >                    stressTensor -= outProduct(dag, fg);
793 >                    if (doHeatFlux_)
794 >                      fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
795                    }
796                  }
797                }
# Line 786 | Line 809 | namespace OpenMD {
809                      // find the distance between the atom
810                      // and the center of the cutoff group:
811                      dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
812 <                    tau -= outProduct(dag, fg);
812 >                    stressTensor -= outProduct(dag, fg);
813 >                    if (doHeatFlux_)
814 >                      fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
815                    }
816                  }
817                }
818              }
819              //if (!info_->usesAtomicVirial()) {
820 <            //  tau -= outProduct(d_grp, fij);
820 >            //  stressTensor -= outProduct(d_grp, fij);
821 >            //  if (doHeatFlux_)
822 >            //     fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
823              //}
824            }
825          }
# Line 830 | Line 857 | namespace OpenMD {
857  
858      lrPot = longRangePotential.sum();
859  
860 <    //store the tau and long range potential    
860 >    //store the stressTensor and long range potential    
861      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
862      curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
863      curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
# Line 851 | Line 878 | namespace OpenMD {
878        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
879             rb = mol->nextRigidBody(rbIter)) {
880          Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
881 <        tau += rbTau;
881 >        stressTensor += rbTau;
882        }
883      }
884      
885   #ifdef IS_MPI
886 <    Mat3x3d tmpTau(tau);
887 <    MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
888 <                  9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
886 >
887 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9,
888 >                              MPI::REALTYPE, MPI::SUM);
889   #endif
890 <    curSnapshot->setTau(tau);
890 >    curSnapshot->setStressTensor(stressTensor);
891 >    
892    }
893  
894   } //end namespace OpenMD

Diff Legend

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