# | 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 44 | Line 44 | |
44 | * @file ForceManager.cpp | |
45 | * @author tlin | |
46 | * @date 11/09/2004 | |
47 | – | * @time 10:39am |
47 | * @version 1.0 | |
48 | */ | |
49 | ||
# | Line 68 | Line 67 | namespace OpenMD { | |
67 | using namespace std; | |
68 | namespace OpenMD { | |
69 | ||
70 | < | ForceManager::ForceManager(SimInfo * info) : info_(info) { |
70 | > | ForceManager::ForceManager(SimInfo * info) : info_(info), switcher_(NULL), |
71 | > | initialized_(false) { |
72 | forceField_ = info_->getForceField(); | |
73 | interactionMan_ = new InteractionManager(); | |
74 | fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_); | |
75 | + | thermo = new Thermo(info_); |
76 | } | |
77 | ||
78 | + | ForceManager::~ForceManager() { |
79 | + | perturbations_.clear(); |
80 | + | |
81 | + | delete switcher_; |
82 | + | delete interactionMan_; |
83 | + | delete fDecomp_; |
84 | + | delete thermo; |
85 | + | } |
86 | + | |
87 | /** | |
88 | * setupCutoffs | |
89 | * | |
# | Line 88 | Line 98 | namespace OpenMD { | |
98 | * simulation for suggested cutoff values (e.g. 2.5 * sigma). | |
99 | * Use the maximum suggested value that was found. | |
100 | * | |
101 | < | * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, |
102 | < | * or SHIFTED_POTENTIAL) |
101 | > | * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, TAYLOR_SHIFTED, |
102 | > | * SHIFTED_POTENTIAL, or EWALD_FULL) |
103 | * If cutoffMethod was explicitly set, use that choice. | |
104 | * If cutoffMethod was not explicitly set, use SHIFTED_FORCE | |
105 | * | |
# | Line 118 | Line 128 | namespace OpenMD { | |
128 | else | |
129 | mdFileVersion = 0; | |
130 | ||
131 | + | // We need the list of simulated atom types to figure out cutoffs |
132 | + | // as well as long range corrections. |
133 | + | |
134 | + | set<AtomType*>::iterator i; |
135 | + | set<AtomType*> atomTypes_; |
136 | + | atomTypes_ = info_->getSimulatedAtomTypes(); |
137 | + | |
138 | if (simParams_->haveCutoffRadius()) { | |
139 | rCut_ = simParams_->getCutoffRadius(); | |
140 | } else { | |
# | Line 132 | Line 149 | namespace OpenMD { | |
149 | rCut_ = 12.0; | |
150 | } else { | |
151 | RealType thisCut; | |
152 | < | set<AtomType*>::iterator i; |
136 | < | set<AtomType*> atomTypes; |
137 | < | atomTypes = info_->getSimulatedAtomTypes(); |
138 | < | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
152 | > | for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) { |
153 | thisCut = interactionMan_->getSuggestedCutoffRadius((*i)); | |
154 | rCut_ = max(thisCut, rCut_); | |
155 | } | |
# | Line 157 | Line 171 | namespace OpenMD { | |
171 | stringToCutoffMethod["SWITCHED"] = SWITCHED; | |
172 | stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL; | |
173 | stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; | |
174 | + | stringToCutoffMethod["TAYLOR_SHIFTED"] = TAYLOR_SHIFTED; |
175 | + | stringToCutoffMethod["EWALD_FULL"] = EWALD_FULL; |
176 | ||
177 | if (simParams_->haveCutoffMethod()) { | |
178 | string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); | |
# | Line 166 | Line 182 | namespace OpenMD { | |
182 | sprintf(painCave.errMsg, | |
183 | "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n" | |
184 | "\tShould be one of: " | |
185 | < | "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", |
185 | > | "HARD, SWITCHED, SHIFTED_POTENTIAL, TAYLOR_SHIFTED,\n" |
186 | > | "\tSHIFTED_FORCE, or EWALD_FULL\n", |
187 | cutMeth.c_str()); | |
188 | painCave.isFatal = 1; | |
189 | painCave.severity = OPENMD_ERROR; | |
# | Line 210 | Line 227 | namespace OpenMD { | |
227 | cutoffMethod_ = SHIFTED_POTENTIAL; | |
228 | } else if (myMethod == "SHIFTED_FORCE") { | |
229 | cutoffMethod_ = SHIFTED_FORCE; | |
230 | + | } else if (myMethod == "TAYLOR_SHIFTED") { |
231 | + | cutoffMethod_ = TAYLOR_SHIFTED; |
232 | + | } else if (myMethod == "EWALD_FULL") { |
233 | + | cutoffMethod_ = EWALD_FULL; |
234 | } | |
235 | ||
236 | if (simParams_->haveSwitchingRadius()) | |
237 | rSwitch_ = simParams_->getSwitchingRadius(); | |
238 | ||
239 | < | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { |
239 | > | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE" || |
240 | > | myMethod == "TAYLOR_SHIFTED" || myMethod == "EWALD_FULL") { |
241 | if (simParams_->haveSwitchingRadius()){ | |
242 | sprintf(painCave.errMsg, | |
243 | "ForceManager::setupCutoffs : DEPRECATED ERROR MESSAGE\n" | |
# | Line 370 | Line 392 | namespace OpenMD { | |
392 | } | |
393 | switcher_->setSwitchType(sft_); | |
394 | switcher_->setSwitch(rSwitch_, rCut_); | |
373 | – | interactionMan_->setSwitchingRadius(rSwitch_); |
395 | } | |
396 | ||
397 | ||
# | Line 394 | Line 415 | namespace OpenMD { | |
415 | doParticlePot_ = info_->getSimParams()->getOutputParticlePotential(); | |
416 | doHeatFlux_ = info_->getSimParams()->getPrintHeatFlux(); | |
417 | if (doHeatFlux_) doParticlePot_ = true; | |
418 | + | |
419 | + | doElectricField_ = info_->getSimParams()->getOutputElectricField(); |
420 | ||
421 | } | |
422 | ||
# | Line 429 | Line 452 | namespace OpenMD { | |
452 | perturbations_.push_back(eField); | |
453 | } | |
454 | ||
455 | + | usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions(); |
456 | + | |
457 | fDecomp_->distributeInitialData(); | |
458 | < | |
458 | > | |
459 | initialized_ = true; | |
460 | < | |
460 | > | |
461 | } | |
462 | < | |
462 | > | |
463 | void ForceManager::calcForces() { | |
464 | ||
465 | if (!initialized_) initialize(); | |
466 | < | |
466 | > | |
467 | preCalculation(); | |
468 | shortRangeInteractions(); | |
469 | longRangeInteractions(); | |
# | Line 637 | Line 662 | namespace OpenMD { | |
662 | ||
663 | void ForceManager::longRangeInteractions() { | |
664 | ||
640 | – | |
665 | Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); | |
666 | DataStorage* config = &(curSnapshot->atomData); | |
667 | DataStorage* cgConfig = &(curSnapshot->cgData); | |
668 | ||
669 | + | |
670 | //calculate the center of mass of cutoff group | |
671 | ||
672 | SimInfo::MoleculeIterator mi; | |
# | Line 674 | Line 699 | namespace OpenMD { | |
699 | RealType vij; | |
700 | Vector3d fij, fg, f1; | |
701 | tuple3<RealType, RealType, RealType> cuts; | |
702 | < | RealType rCutSq; |
702 | > | RealType rCut, rCutSq, rListSq; |
703 | bool in_switching_region; | |
704 | RealType sw, dswdr, swderiv; | |
705 | < | vector<int> atomListColumn, atomListRow, atomListLocal; |
705 | > | vector<int> atomListColumn, atomListRow; |
706 | InteractionData idat; | |
707 | SelfData sdat; | |
708 | RealType mf; | |
# | Line 685 | Line 710 | namespace OpenMD { | |
710 | RealType dVdFQ1(0.0); | |
711 | RealType dVdFQ2(0.0); | |
712 | potVec longRangePotential(0.0); | |
713 | + | potVec reciprocalPotential(0.0); |
714 | potVec workPot(0.0); | |
715 | potVec exPot(0.0); | |
716 | + | Vector3d eField1(0.0); |
717 | + | Vector3d eField2(0.0); |
718 | vector<int>::iterator ia, jb; | |
719 | ||
720 | int loopStart, loopEnd; | |
721 | < | |
721 | > | |
722 | > | idat.rcut = &rCut; |
723 | idat.vdwMult = &vdwMult; | |
724 | idat.electroMult = &electroMult; | |
725 | idat.pot = &workPot; | |
# | Line 700 | Line 729 | namespace OpenMD { | |
729 | idat.vpair = &vpair; | |
730 | idat.dVdFQ1 = &dVdFQ1; | |
731 | idat.dVdFQ2 = &dVdFQ2; | |
732 | + | idat.eField1 = &eField1; |
733 | + | idat.eField2 = &eField2; |
734 | idat.f1 = &f1; | |
735 | idat.sw = &sw; | |
736 | idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false; | |
737 | < | idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false; |
737 | > | idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE || cutoffMethod_ == TAYLOR_SHIFTED) ? true : false; |
738 | idat.doParticlePot = doParticlePot_; | |
739 | + | idat.doElectricField = doElectricField_; |
740 | sdat.doParticlePot = doParticlePot_; | |
741 | ||
742 | loopEnd = PAIR_LOOP; | |
# | Line 717 | Line 749 | namespace OpenMD { | |
749 | ||
750 | if (iLoop == loopStart) { | |
751 | bool update_nlist = fDecomp_->checkNeighborList(); | |
752 | < | if (update_nlist) |
753 | < | neighborList = fDecomp_->buildNeighborList(); |
754 | < | } |
752 | > | if (update_nlist) { |
753 | > | if (!usePeriodicBoundaryConditions_) |
754 | > | Mat3x3d bbox = thermo->getBoundingBox(); |
755 | > | fDecomp_->buildNeighborList(neighborList_); |
756 | > | } |
757 | > | } |
758 | ||
759 | < | for (vector<pair<int, int> >::iterator it = neighborList.begin(); |
760 | < | it != neighborList.end(); ++it) { |
759 | > | for (vector<pair<int, int> >::iterator it = neighborList_.begin(); |
760 | > | it != neighborList_.end(); ++it) { |
761 | ||
762 | cg1 = (*it).first; | |
763 | cg2 = (*it).second; | |
764 | ||
765 | < | cuts = fDecomp_->getGroupCutoffs(cg1, cg2); |
765 | > | fDecomp_->getGroupCutoffs(cg1, cg2, rCut, rCutSq, rListSq); |
766 | ||
767 | d_grp = fDecomp_->getIntergroupVector(cg1, cg2); | |
768 | ||
769 | < | curSnapshot->wrapVector(d_grp); |
769 | > | // already wrapped in the getIntergroupVector call: |
770 | > | // curSnapshot->wrapVector(d_grp); |
771 | rgrpsq = d_grp.lengthSquare(); | |
736 | – | rCutSq = cuts.second; |
772 | ||
773 | if (rgrpsq < rCutSq) { | |
739 | – | idat.rcut = &cuts.first; |
774 | if (iLoop == PAIR_LOOP) { | |
775 | vij = 0.0; | |
776 | < | fij = V3Zero; |
776 | > | fij.zero(); |
777 | > | eField1.zero(); |
778 | > | eField2.zero(); |
779 | } | |
780 | ||
781 | in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, | |
# | Line 764 | Line 800 | namespace OpenMD { | |
800 | vpair = 0.0; | |
801 | workPot = 0.0; | |
802 | exPot = 0.0; | |
803 | < | f1 = V3Zero; |
803 | > | f1.zero(); |
804 | dVdFQ1 = 0.0; | |
805 | dVdFQ2 = 0.0; | |
806 | ||
# | Line 791 | Line 827 | namespace OpenMD { | |
827 | ||
828 | r = sqrt( *(idat.r2) ); | |
829 | idat.rij = &r; | |
830 | < | |
830 | > | |
831 | if (iLoop == PREPAIR_LOOP) { | |
832 | interactionMan_->doPrePair(idat); | |
833 | } else { | |
# | Line 814 | Line 850 | namespace OpenMD { | |
850 | fij += fg; | |
851 | ||
852 | if (atomListRow.size() == 1 && atomListColumn.size() == 1) { | |
853 | < | stressTensor -= outProduct( *(idat.d), fg); |
854 | < | if (doHeatFlux_) |
855 | < | fDecomp_->addToHeatFlux(*(idat.d) * dot(fg, vel2)); |
856 | < | |
853 | > | if (!fDecomp_->skipAtomPair(atomListRow[0], |
854 | > | atomListColumn[0], |
855 | > | cg1, cg2)) { |
856 | > | stressTensor -= outProduct( *(idat.d), fg); |
857 | > | if (doHeatFlux_) |
858 | > | fDecomp_->addToHeatFlux(*(idat.d) * dot(fg, vel2)); |
859 | > | } |
860 | } | |
861 | ||
862 | for (ia = atomListRow.begin(); | |
# | Line 884 | Line 923 | namespace OpenMD { | |
923 | } | |
924 | } | |
925 | } | |
926 | < | |
926 | > | |
927 | // collects pairwise information | |
928 | fDecomp_->collectData(); | |
929 | + | if (cutoffMethod_ == EWALD_FULL) { |
930 | + | interactionMan_->doReciprocalSpaceSum(reciprocalPotential); |
931 | + | } |
932 | ||
933 | if (info_->requiresSelfCorrection()) { | |
934 | for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { | |
# | Line 899 | Line 941 | namespace OpenMD { | |
941 | fDecomp_->collectSelfData(); | |
942 | ||
943 | longRangePotential = *(fDecomp_->getEmbeddingPotential()) + | |
944 | < | *(fDecomp_->getPairwisePotential()); |
944 | > | *(fDecomp_->getPairwisePotential()) + reciprocalPotential; |
945 | ||
946 | curSnapshot->setLongRangePotential(longRangePotential); | |
947 | ||
# | Line 908 | Line 950 | namespace OpenMD { | |
950 | ||
951 | } | |
952 | ||
911 | – | |
953 | void ForceManager::postCalculation() { | |
954 | ||
955 | vector<Perturbation*>::iterator pi; | |
# | Line 939 | Line 980 | namespace OpenMD { | |
980 | #endif | |
981 | curSnapshot->setStressTensor(stressTensor); | |
982 | ||
983 | + | if (info_->getSimParams()->getUseLongRangeCorrections()) { |
984 | + | /* |
985 | + | RealType vol = curSnapshot->getVolume(); |
986 | + | RealType Elrc(0.0); |
987 | + | RealType Wlrc(0.0); |
988 | + | |
989 | + | set<AtomType*>::iterator i; |
990 | + | set<AtomType*>::iterator j; |
991 | + | |
992 | + | RealType n_i, n_j; |
993 | + | RealType rho_i, rho_j; |
994 | + | pair<RealType, RealType> LRI; |
995 | + | |
996 | + | for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) { |
997 | + | n_i = RealType(info_->getGlobalCountOfType(*i)); |
998 | + | rho_i = n_i / vol; |
999 | + | for (j = atomTypes_.begin(); j != atomTypes_.end(); ++j) { |
1000 | + | n_j = RealType(info_->getGlobalCountOfType(*j)); |
1001 | + | rho_j = n_j / vol; |
1002 | + | |
1003 | + | LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) ); |
1004 | + | |
1005 | + | Elrc += n_i * rho_j * LRI.first; |
1006 | + | Wlrc -= rho_i * rho_j * LRI.second; |
1007 | + | } |
1008 | + | } |
1009 | + | Elrc *= 2.0 * NumericConstant::PI; |
1010 | + | Wlrc *= 2.0 * NumericConstant::PI; |
1011 | + | |
1012 | + | RealType lrp = curSnapshot->getLongRangePotential(); |
1013 | + | curSnapshot->setLongRangePotential(lrp + Elrc); |
1014 | + | stressTensor += Wlrc * SquareMatrix3<RealType>::identity(); |
1015 | + | curSnapshot->setStressTensor(stressTensor); |
1016 | + | */ |
1017 | + | |
1018 | + | } |
1019 | } | |
1020 | < | } //end namespace OpenMD |
1020 | > | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |