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 1711 by gezelter, Sat May 19 02:58:35 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 <      cerr << "dPP = " << doParticlePot_ << "\n";
393 >      doHeatFlux_ = info_->getSimParams()->getPrintHeatFlux();
394 >      if (doHeatFlux_) doParticlePot_ = true;
395    
396      }
397  
# Line 472 | 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 506 | Line 508 | namespace OpenMD {
508  
509        for (bond = mol->beginBond(bondIter); bond != NULL;
510             bond = mol->nextBond(bondIter)) {
511 <        bond->calcForce();
511 >        bond->calcForce(doParticlePot_);
512          bondPotential += bond->getPotential();
513        }
514  
# Line 514 | Line 516 | namespace OpenMD {
516             bend = mol->nextBend(bendIter)) {
517          
518          RealType angle;
519 <        bend->calcForce(angle);
519 >        bend->calcForce(angle, doParticlePot_);
520          RealType currBendPot = bend->getPotential();          
521          
522          bendPotential += bend->getPotential();
# Line 539 | Line 541 | namespace OpenMD {
541        for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
542             torsion = mol->nextTorsion(torsionIter)) {
543          RealType angle;
544 <        torsion->calcForce(angle);
544 >        torsion->calcForce(angle, doParticlePot_);
545          RealType currTorsionPot = torsion->getPotential();
546          torsionPotential += torsion->getPotential();
547          map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
# Line 563 | Line 565 | namespace OpenMD {
565             inversion != NULL;
566             inversion = mol->nextInversion(inversionIter)) {
567          RealType angle;
568 <        inversion->calcForce(angle);
568 >        inversion->calcForce(angle, doParticlePot_);
569          RealType currInversionPot = inversion->getPotential();
570          inversionPotential += inversion->getPotential();
571          map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
# Line 596 | 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 619 | 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 640 | 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;
652  
653      int loopStart, loopEnd;
654  
# 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 <          for (vector<int>::iterator ia = atomListRow.begin();
709 >          if (doHeatFlux_)
710 >            gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
711 >
712 >          for (ia = atomListRow.begin();
713                 ia != atomListRow.end(); ++ia) {            
714              atom1 = (*ia);
715 <            
716 <            for (vector<int>::iterator jb = atomListColumn.begin();
715 >
716 >            for (jb = atomListColumn.begin();
717                   jb != atomListColumn.end(); ++jb) {              
718                atom2 = (*jb);
719  
# 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 (vector<int>::iterator ia = atomListRow.begin();
779 >              for (ia = atomListRow.begin();
780                     ia != atomListRow.end(); ++ia) {            
781                  atom1 = (*ia);                
782                  mf = fDecomp_->getMassFactorRow(atom1);
# 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                }
798 <              for (vector<int>::iterator jb = atomListColumn.begin();
798 >              for (jb = atomListColumn.begin();
799                     jb != atomListColumn.end(); ++jb) {              
800                  atom2 = (*jb);
801                  mf = fDecomp_->getMassFactorColumn(atom2);
# 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