| 57 |  | #include "primitives/Torsion.hpp" | 
| 58 |  | #include "primitives/Inversion.hpp" | 
| 59 |  | #include "nonbonded/NonBondedInteraction.hpp" | 
| 60 | < | #include "perturbations/ElectricField.hpp" | 
| 60 | > | #include "perturbations/UniformField.hpp" | 
| 61 |  | #include "parallel/ForceMatrixDecomposition.hpp" | 
| 62 |  |  | 
| 63 |  | #include <cstdio> | 
| 405 |  | interactionMan_->setSimInfo(info_); | 
| 406 |  | interactionMan_->initialize(); | 
| 407 |  |  | 
| 408 | < | // We want to delay the cutoffs until after the interaction | 
| 409 | < | // manager has set up the atom-atom interactions so that we can | 
| 410 | < | // query them for suggested cutoff values | 
| 408 | > | //! We want to delay the cutoffs until after the interaction | 
| 409 | > | //! manager has set up the atom-atom interactions so that we can | 
| 410 | > | //! query them for suggested cutoff values | 
| 411 |  | setupCutoffs(); | 
| 412 |  |  | 
| 413 |  | info_->prepareTopology(); | 
| 417 |  | if (doHeatFlux_) doParticlePot_ = true; | 
| 418 |  |  | 
| 419 |  | doElectricField_ = info_->getSimParams()->getOutputElectricField(); | 
| 420 | + | doSitePotential_ = info_->getSimParams()->getOutputSitePotential(); | 
| 421 |  |  | 
| 422 |  | } | 
| 423 |  |  | 
| 424 |  | ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); | 
| 425 |  |  | 
| 426 | < | // Force fields can set options on how to scale van der Waals and | 
| 427 | < | // electrostatic interactions for atoms connected via bonds, bends | 
| 428 | < | // and torsions in this case the topological distance between | 
| 429 | < | // atoms is: | 
| 430 | < | // 0 = topologically unconnected | 
| 431 | < | // 1 = bonded together | 
| 432 | < | // 2 = connected via a bend | 
| 433 | < | // 3 = connected via a torsion | 
| 426 | > | //! Force fields can set options on how to scale van der Waals and | 
| 427 | > | //! electrostatic interactions for atoms connected via bonds, bends | 
| 428 | > | //! and torsions in this case the topological distance between | 
| 429 | > | //! atoms is: | 
| 430 | > | //! 0 = topologically unconnected | 
| 431 | > | //! 1 = bonded together | 
| 432 | > | //! 2 = connected via a bend | 
| 433 | > | //! 3 = connected via a torsion | 
| 434 |  |  | 
| 435 |  | vdwScale_.reserve(4); | 
| 436 |  | fill(vdwScale_.begin(), vdwScale_.end(), 0.0); | 
| 448 |  | electrostaticScale_[2] = fopts.getelectrostatic13scale(); | 
| 449 |  | electrostaticScale_[3] = fopts.getelectrostatic14scale(); | 
| 450 |  |  | 
| 451 | < | if (info_->getSimParams()->haveElectricField()) { | 
| 452 | < | ElectricField* eField = new ElectricField(info_); | 
| 451 | > | if (info_->getSimParams()->haveUniformField()) { | 
| 452 | > | UniformField* eField = new UniformField(info_); | 
| 453 |  | perturbations_.push_back(eField); | 
| 454 |  | } | 
| 455 |  |  | 
| 638 |  | // Collect from all nodes.  This should eventually be moved into a | 
| 639 |  | // SystemDecomposition, but this is a better place than in | 
| 640 |  | // Thermo to do the collection. | 
| 641 | < | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE, | 
| 642 | < | MPI::SUM); | 
| 643 | < | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE, | 
| 644 | < | MPI::SUM); | 
| 645 | < | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1, | 
| 646 | < | MPI::REALTYPE, MPI::SUM); | 
| 647 | < | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1, | 
| 648 | < | MPI::REALTYPE, MPI::SUM); | 
| 641 | > |  | 
| 642 | > | MPI_Allreduce(MPI_IN_PLACE, &bondPotential, 1, MPI_REALTYPE, | 
| 643 | > | MPI_SUM, MPI_COMM_WORLD); | 
| 644 | > | MPI_Allreduce(MPI_IN_PLACE, &bendPotential, 1, MPI_REALTYPE, | 
| 645 | > | MPI_SUM, MPI_COMM_WORLD); | 
| 646 | > | MPI_Allreduce(MPI_IN_PLACE, &torsionPotential, 1, | 
| 647 | > | MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 648 | > | MPI_Allreduce(MPI_IN_PLACE, &inversionPotential, 1, | 
| 649 | > | MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 650 |  | #endif | 
| 651 |  |  | 
| 652 |  | Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); | 
| 668 |  | DataStorage* config = &(curSnapshot->atomData); | 
| 669 |  | DataStorage* cgConfig = &(curSnapshot->cgData); | 
| 670 |  |  | 
| 671 | + |  | 
| 672 |  | //calculate the center of mass of cutoff group | 
| 673 |  |  | 
| 674 |  | SimInfo::MoleculeIterator mi; | 
| 712 |  | RealType dVdFQ1(0.0); | 
| 713 |  | RealType dVdFQ2(0.0); | 
| 714 |  | potVec longRangePotential(0.0); | 
| 715 | + | RealType reciprocalPotential(0.0); | 
| 716 |  | potVec workPot(0.0); | 
| 717 |  | potVec exPot(0.0); | 
| 718 |  | Vector3d eField1(0.0); | 
| 719 |  | Vector3d eField2(0.0); | 
| 720 | + | RealType sPot1(0.0); | 
| 721 | + | RealType sPot2(0.0); | 
| 722 | + |  | 
| 723 |  | vector<int>::iterator ia, jb; | 
| 724 |  |  | 
| 725 |  | int loopStart, loopEnd; | 
| 735 |  | idat.dVdFQ1 = &dVdFQ1; | 
| 736 |  | idat.dVdFQ2 = &dVdFQ2; | 
| 737 |  | idat.eField1 = &eField1; | 
| 738 | < | idat.eField2 = &eField2; | 
| 738 | > | idat.eField2 = &eField2; | 
| 739 | > | idat.sPot1 = &sPot1; | 
| 740 | > | idat.sPot2 = &sPot2; | 
| 741 |  | idat.f1 = &f1; | 
| 742 |  | idat.sw = &sw; | 
| 743 |  | idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false; | 
| 744 |  | idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE || cutoffMethod_ == TAYLOR_SHIFTED) ? true : false; | 
| 745 |  | idat.doParticlePot = doParticlePot_; | 
| 746 |  | idat.doElectricField = doElectricField_; | 
| 747 | + | idat.doSitePotential = doSitePotential_; | 
| 748 |  | sdat.doParticlePot = doParticlePot_; | 
| 749 |  |  | 
| 750 |  | loopEnd = PAIR_LOOP; | 
| 779 |  | rgrpsq = d_grp.lengthSquare(); | 
| 780 |  |  | 
| 781 |  | if (rgrpsq < rCutSq) { | 
| 772 | – |  | 
| 782 |  | if (iLoop == PAIR_LOOP) { | 
| 783 |  | vij = 0.0; | 
| 784 |  | fij.zero(); | 
| 785 |  | eField1.zero(); | 
| 786 |  | eField2.zero(); | 
| 787 | + | sPot1 = 0.0; | 
| 788 | + | sPot2 = 0.0; | 
| 789 |  | } | 
| 790 |  |  | 
| 791 |  | in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, | 
| 837 |  |  | 
| 838 |  | r = sqrt( *(idat.r2) ); | 
| 839 |  | idat.rij = &r; | 
| 840 | < |  | 
| 840 | > |  | 
| 841 |  | if (iLoop == PREPAIR_LOOP) { | 
| 842 |  | interactionMan_->doPrePair(idat); | 
| 843 |  | } else { | 
| 937 |  | // collects pairwise information | 
| 938 |  | fDecomp_->collectData(); | 
| 939 |  | if (cutoffMethod_ == EWALD_FULL) { | 
| 940 | < | interactionMan_->doReciprocalSpaceSum(); | 
| 940 | > | interactionMan_->doReciprocalSpaceSum(reciprocalPotential); | 
| 941 | > |  | 
| 942 | > | curSnapshot->setReciprocalPotential(reciprocalPotential); | 
| 943 |  | } | 
| 944 |  |  | 
| 945 |  | if (info_->requiresSelfCorrection()) { | 
| 962 |  |  | 
| 963 |  | } | 
| 964 |  |  | 
| 952 | – |  | 
| 965 |  | void ForceManager::postCalculation() { | 
| 966 |  |  | 
| 967 |  | vector<Perturbation*>::iterator pi; | 
| 987 |  | } | 
| 988 |  |  | 
| 989 |  | #ifdef IS_MPI | 
| 990 | < | MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9, | 
| 991 | < | MPI::REALTYPE, MPI::SUM); | 
| 990 | > | MPI_Allreduce(MPI_IN_PLACE, stressTensor.getArrayPointer(), 9, | 
| 991 | > | MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 992 |  | #endif | 
| 993 |  | curSnapshot->setStressTensor(stressTensor); | 
| 994 |  |  |