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> |
370 |
|
} |
371 |
|
switcher_->setSwitchType(sft_); |
372 |
|
switcher_->setSwitch(rSwitch_, rCut_); |
372 |
– |
interactionMan_->setSwitchingRadius(rSwitch_); |
373 |
|
} |
374 |
|
|
375 |
|
|
393 |
|
doParticlePot_ = info_->getSimParams()->getOutputParticlePotential(); |
394 |
|
doHeatFlux_ = info_->getSimParams()->getPrintHeatFlux(); |
395 |
|
if (doHeatFlux_) doParticlePot_ = true; |
396 |
+ |
|
397 |
+ |
doElectricField_ = info_->getSimParams()->getOutputElectricField(); |
398 |
|
|
399 |
|
} |
400 |
|
|
425 |
|
electrostaticScale_[2] = fopts.getelectrostatic13scale(); |
426 |
|
electrostaticScale_[3] = fopts.getelectrostatic14scale(); |
427 |
|
|
428 |
+ |
if (info_->getSimParams()->haveElectricField()) { |
429 |
+ |
ElectricField* eField = new ElectricField(info_); |
430 |
+ |
perturbations_.push_back(eField); |
431 |
+ |
} |
432 |
+ |
|
433 |
|
fDecomp_->distributeInitialData(); |
434 |
|
|
435 |
|
initialized_ = true; |
682 |
|
InteractionData idat; |
683 |
|
SelfData sdat; |
684 |
|
RealType mf; |
678 |
– |
RealType lrPot; |
685 |
|
RealType vpair; |
686 |
|
RealType dVdFQ1(0.0); |
687 |
|
RealType dVdFQ2(0.0); |
688 |
|
potVec longRangePotential(0.0); |
689 |
|
potVec workPot(0.0); |
690 |
|
potVec exPot(0.0); |
691 |
+ |
Vector3d eField1(0.0); |
692 |
+ |
Vector3d eField2(0.0); |
693 |
|
vector<int>::iterator ia, jb; |
694 |
|
|
695 |
|
int loopStart, loopEnd; |
703 |
|
idat.vpair = &vpair; |
704 |
|
idat.dVdFQ1 = &dVdFQ1; |
705 |
|
idat.dVdFQ2 = &dVdFQ2; |
706 |
+ |
idat.eField1 = &eField1; |
707 |
+ |
idat.eField2 = &eField2; |
708 |
|
idat.f1 = &f1; |
709 |
|
idat.sw = &sw; |
710 |
|
idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false; |
711 |
|
idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false; |
712 |
|
idat.doParticlePot = doParticlePot_; |
713 |
+ |
idat.doElectricField = doElectricField_; |
714 |
|
sdat.doParticlePot = doParticlePot_; |
715 |
|
|
716 |
|
loopEnd = PAIR_LOOP; |
880 |
|
|
881 |
|
fDecomp_->collectIntermediateData(); |
882 |
|
|
883 |
< |
for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { |
883 |
> |
for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { |
884 |
|
fDecomp_->fillSelfData(sdat, atom1); |
885 |
|
interactionMan_->doPreForce(sdat); |
886 |
|
} |
895 |
|
fDecomp_->collectData(); |
896 |
|
|
897 |
|
if (info_->requiresSelfCorrection()) { |
898 |
< |
for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { |
898 |
> |
for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { |
899 |
|
fDecomp_->fillSelfData(sdat, atom1); |
900 |
|
interactionMan_->doSelfCorrection(sdat); |
901 |
|
} |
909 |
|
|
910 |
|
curSnapshot->setLongRangePotential(longRangePotential); |
911 |
|
|
912 |
< |
// lrPot = longRangePotential.sum(); |
912 |
> |
// collects single-atom information |
913 |
> |
fDecomp_->collectSelfData(); |
914 |
|
|
915 |
< |
// //store the long range potential |
916 |
< |
// curSnapshot->setLongRangePotential(lrPot); |
915 |
> |
longRangePotential = *(fDecomp_->getEmbeddingPotential()) + |
916 |
> |
*(fDecomp_->getPairwisePotential()); |
917 |
> |
|
918 |
> |
curSnapshot->setLongRangePotential(longRangePotential); |
919 |
|
|
920 |
|
curSnapshot->setExcludedPotentials(*(fDecomp_->getExcludedSelfPotential()) + |
921 |
|
*(fDecomp_->getExcludedPotential())); |
924 |
|
|
925 |
|
|
926 |
|
void ForceManager::postCalculation() { |
927 |
+ |
|
928 |
+ |
vector<Perturbation*>::iterator pi; |
929 |
+ |
for (pi = perturbations_.begin(); pi != perturbations_.end(); ++pi) { |
930 |
+ |
(*pi)->applyPerturbation(); |
931 |
+ |
} |
932 |
+ |
|
933 |
|
SimInfo::MoleculeIterator mi; |
934 |
|
Molecule* mol; |
935 |
|
Molecule::RigidBodyIterator rbIter; |
954 |
|
curSnapshot->setStressTensor(stressTensor); |
955 |
|
|
956 |
|
} |
937 |
– |
|
957 |
|
} //end namespace OpenMD |