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 1754 by jmichalk, Wed Jun 13 14:45:59 2012 UTC vs.
Revision 1780 by jmarr, Mon Aug 20 18:28:22 2012 UTC

# Line 58 | Line 58
58   #include "primitives/Torsion.hpp"
59   #include "primitives/Inversion.hpp"
60   #include "nonbonded/NonBondedInteraction.hpp"
61 + #include "perturbations/ElectricField.hpp"
62   #include "parallel/ForceMatrixDecomposition.hpp"
63  
64   #include <cstdio>
# Line 111 | Line 112 | namespace OpenMD {
112      ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
113      int mdFileVersion;
114      rCut_ = 0.0; //Needs a value for a later max() call;  
115 <
115 >    
116      if (simParams_->haveMDfileVersion())
117        mdFileVersion = simParams_->getMDfileVersion();
118      else
# Line 423 | Line 424 | namespace OpenMD {
424      electrostaticScale_[2] = fopts.getelectrostatic13scale();
425      electrostaticScale_[3] = fopts.getelectrostatic14scale();    
426      
427 +    if (info_->getSimParams()->haveElectricField()) {
428 +      ElectricField* eField = new ElectricField(info_);
429 +      perturbations_.push_back(eField);
430 +    }
431 +
432      fDecomp_->distributeInitialData();
433  
434      initialized_ = true;
# Line 449 | Line 455 | namespace OpenMD {
455      Molecule::CutoffGroupIterator ci;
456      CutoffGroup* cg;
457      
458 <    // forces are zeroed here, before any are accumulated.
458 >    // forces and potentials are zeroed here, before any are
459 >    // accumulated.
460      
461 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
462 +
463 +    snap->setBondPotential(0.0);
464 +    snap->setBendPotential(0.0);
465 +    snap->setTorsionPotential(0.0);
466 +    snap->setInversionPotential(0.0);
467 +
468 +    potVec zeroPot(0.0);
469 +    snap->setLongRangePotential(zeroPot);
470 +    snap->setExcludedPotentials(zeroPot);
471 +
472 +    snap->setRestraintPotential(0.0);
473 +    snap->setRawPotential(0.0);
474 +
475      for (mol = info_->beginMolecule(mi); mol != NULL;
476           mol = info_->nextMolecule(mi)) {
477        for(atom = mol->beginAtom(ai); atom != NULL;
# Line 586 | Line 607 | namespace OpenMD {
607          }      
608        }      
609      }
610 <    
611 <    RealType  shortRangePotential = bondPotential + bendPotential +
612 <      torsionPotential +  inversionPotential;    
610 >
611 > #ifdef IS_MPI
612 >    // Collect from all nodes.  This should eventually be moved into a
613 >    // SystemDecomposition, but this is a better place than in
614 >    // Thermo to do the collection.
615 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE,
616 >                              MPI::SUM);
617 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE,
618 >                              MPI::SUM);
619 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1,
620 >                              MPI::REALTYPE, MPI::SUM);
621 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1,
622 >                              MPI::REALTYPE, MPI::SUM);
623 > #endif
624 >
625      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
626 <    curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
627 <    curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
628 <    curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
629 <    curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
630 <    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;    
626 >
627 >    curSnapshot->setBondPotential(bondPotential);
628 >    curSnapshot->setBendPotential(bendPotential);
629 >    curSnapshot->setTorsionPotential(torsionPotential);
630 >    curSnapshot->setInversionPotential(inversionPotential);
631 >    
632 >    // RealType shortRangePotential = bondPotential + bendPotential +
633 >    //   torsionPotential +  inversionPotential;    
634 >
635 >    // curSnapshot->setShortRangePotential(shortRangePotential);
636    }
637    
638    void ForceManager::longRangeInteractions() {
# Line 643 | Line 681 | namespace OpenMD {
681      InteractionData idat;
682      SelfData sdat;
683      RealType mf;
646 –    RealType lrPot;
684      RealType vpair;
685      RealType dVdFQ1(0.0);
686      RealType dVdFQ2(0.0);
687      potVec longRangePotential(0.0);
688      potVec workPot(0.0);
689 +    potVec exPot(0.0);
690      vector<int>::iterator ia, jb;
691  
692      int loopStart, loopEnd;
# Line 656 | Line 694 | namespace OpenMD {
694      idat.vdwMult = &vdwMult;
695      idat.electroMult = &electroMult;
696      idat.pot = &workPot;
697 +    idat.excludedPot = &exPot;
698      sdat.pot = fDecomp_->getEmbeddingPotential();
699 +    sdat.excludedPot = fDecomp_->getExcludedSelfPotential();
700      idat.vpair = &vpair;
701      idat.dVdFQ1 = &dVdFQ1;
702      idat.dVdFQ2 = &dVdFQ2;
# Line 704 | Line 744 | namespace OpenMD {
744            
745            in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
746                                                       rgrp);
747 +
748            atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
749            atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
750  
# Line 718 | Line 759 | namespace OpenMD {
759                   jb != atomListColumn.end(); ++jb) {              
760                atom2 = (*jb);
761  
762 <              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
762 >              if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
763 >
764                  vpair = 0.0;
765                  workPot = 0.0;
766 +                exPot = 0.0;
767                  f1 = V3Zero;
768                  dVdFQ1 = 0.0;
769                  dVdFQ2 = 0.0;
# Line 831 | Line 874 | namespace OpenMD {
874  
875            fDecomp_->collectIntermediateData();
876  
877 <          for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
877 >          for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
878              fDecomp_->fillSelfData(sdat, atom1);
879              interactionMan_->doPreForce(sdat);
880            }
# Line 842 | Line 885 | namespace OpenMD {
885        }
886      }
887      
888 +    // collects pairwise information
889      fDecomp_->collectData();
890          
891      if (info_->requiresSelfCorrection()) {
892 <
849 <      for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
892 >      for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
893          fDecomp_->fillSelfData(sdat, atom1);
894          interactionMan_->doSelfCorrection(sdat);
895        }
853 –
896      }
897  
898 +    // collects single-atom information
899 +    fDecomp_->collectSelfData();
900 +
901      longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
902        *(fDecomp_->getPairwisePotential());
903  
904 <    lrPot = longRangePotential.sum();
904 >    curSnapshot->setLongRangePotential(longRangePotential);
905 >    
906 >    curSnapshot->setExcludedPotentials(*(fDecomp_->getExcludedSelfPotential()) +
907 >                                         *(fDecomp_->getExcludedPotential()));
908  
861 –    //store the stressTensor and long range potential    
862 –    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
863 –    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
864 –    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
909    }
910  
911    
912    void ForceManager::postCalculation() {
913 +
914 +    vector<Perturbation*>::iterator pi;
915 +    for (pi = perturbations_.begin(); pi != perturbations_.end(); ++pi) {
916 +      (*pi)->applyPerturbation();
917 +    }
918 +
919      SimInfo::MoleculeIterator mi;
920      Molecule* mol;
921      Molecule::RigidBodyIterator rbIter;
922      RigidBody* rb;
923      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
924 <    
924 >  
925      // collect the atomic forces onto rigid bodies
926      
927      for (mol = info_->beginMolecule(mi); mol != NULL;
# Line 884 | Line 934 | namespace OpenMD {
934      }
935      
936   #ifdef IS_MPI
887 –
937      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9,
938                                MPI::REALTYPE, MPI::SUM);
939   #endif
940      curSnapshot->setStressTensor(stressTensor);
941      
942    }
894 –
943   } //end namespace OpenMD

Diff Legend

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