| 59 |  | #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" | 
| 60 |  | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" | 
| 61 |  | #include "UseTheForce/doForces_interface.h" | 
| 62 | + | #include "UseTheForce/DarkSide/neighborLists_interface.h" | 
| 63 |  | #include "UseTheForce/DarkSide/electrostatic_interface.h" | 
| 64 |  | #include "UseTheForce/DarkSide/switcheroo_interface.h" | 
| 65 |  | #include "utils/MemoryUtils.hpp" | 
| 68 |  | #include "io/ForceFieldOptions.hpp" | 
| 69 |  | #include "UseTheForce/ForceField.hpp" | 
| 70 |  |  | 
| 71 | + |  | 
| 72 |  | #ifdef IS_MPI | 
| 73 |  | #include "UseTheForce/mpiComponentPlan.h" | 
| 74 |  | #include "UseTheForce/DarkSide/simParallel_interface.h" | 
| 92 |  | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | 
| 93 |  | nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0), | 
| 94 |  | nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0), | 
| 95 | < | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) { | 
| 95 | > | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false), | 
| 96 | > | useAtomicVirial_(true) { | 
| 97 |  |  | 
| 98 |  | MoleculeStamp* molStamp; | 
| 99 |  | int nMolWithSameStamp; | 
| 667 |  | int useSF; | 
| 668 |  | int useSP; | 
| 669 |  | int useBoxDipole; | 
| 670 | + |  | 
| 671 |  | std::string myMethod; | 
| 672 |  |  | 
| 673 |  | // set the useRF logical | 
| 674 |  | useRF = 0; | 
| 675 |  | useSF = 0; | 
| 676 | + | useSP = 0; | 
| 677 |  |  | 
| 678 |  |  | 
| 679 |  | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 680 |  | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 681 |  | toUpper(myMethod); | 
| 682 |  | if (myMethod == "REACTION_FIELD"){ | 
| 683 | < | useRF=1; | 
| 683 | > | useRF = 1; | 
| 684 |  | } else if (myMethod == "SHIFTED_FORCE"){ | 
| 685 |  | useSF = 1; | 
| 686 |  | } else if (myMethod == "SHIFTED_POTENTIAL"){ | 
| 692 |  | if (simParams_->getAccumulateBoxDipole()) | 
| 693 |  | useBoxDipole = 1; | 
| 694 |  |  | 
| 695 | + | useAtomicVirial_ = simParams_->getUseAtomicVirial(); | 
| 696 | + |  | 
| 697 |  | //loop over all of the atom types | 
| 698 |  | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 699 |  | useLennardJones |= (*i)->isLennardJones(); | 
| 771 |  | temp = useBoxDipole; | 
| 772 |  | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 773 |  |  | 
| 774 | + | temp = useAtomicVirial_; | 
| 775 | + | MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 776 | + |  | 
| 777 |  | #endif | 
| 778 |  |  | 
| 779 |  | fInfo_.SIM_uses_PBC = usePBC; | 
| 793 |  | fInfo_.SIM_uses_SF = useSF; | 
| 794 |  | fInfo_.SIM_uses_SP = useSP; | 
| 795 |  | fInfo_.SIM_uses_BoxDipole = useBoxDipole; | 
| 796 | + | fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; | 
| 797 |  | } | 
| 798 |  |  | 
| 799 |  | void SimInfo::setupFortranSim() { | 
| 878 |  | "succesfully sent the simulation information to fortran.\n"); | 
| 879 |  | MPIcheckPoint(); | 
| 880 |  | #endif // is_mpi | 
| 881 | + |  | 
| 882 | + | // Setup number of neighbors in neighbor list if present | 
| 883 | + | if (simParams_->haveNeighborListNeighbors()) { | 
| 884 | + | int nlistNeighbors = simParams_->getNeighborListNeighbors(); | 
| 885 | + | setNeighbors(&nlistNeighbors); | 
| 886 | + | } | 
| 887 | + |  | 
| 888 | + |  | 
| 889 |  | } | 
| 890 |  |  | 
| 891 |  |  | 
| 1144 |  | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); | 
| 1145 |  | painCave.isFatal = 0; | 
| 1146 |  | simError(); | 
| 1147 | + | } else { | 
| 1148 | + | alphaVal = simParams_->getDampingAlpha(); | 
| 1149 |  | } | 
| 1150 | + |  | 
| 1151 |  | } else { | 
| 1152 |  | // throw error | 
| 1153 |  | sprintf( painCave.errMsg, | 
| 1464 |  | IOIndexToIntegrableObject= v; | 
| 1465 |  | } | 
| 1466 |  |  | 
| 1467 | + | /* Returns the Volume of the simulation based on a ellipsoid with semi-axes | 
| 1468 | + | based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 | 
| 1469 | + | where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to | 
| 1470 | + | V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. | 
| 1471 | + | */ | 
| 1472 | + | void SimInfo::getGyrationalVolume(RealType &volume){ | 
| 1473 | + | Mat3x3d intTensor; | 
| 1474 | + | RealType det; | 
| 1475 | + | Vector3d dummyAngMom; | 
| 1476 | + | RealType sysconstants; | 
| 1477 | + | RealType geomCnst; | 
| 1478 | + |  | 
| 1479 | + | geomCnst = 3.0/2.0; | 
| 1480 | + | /* Get the inertial tensor and angular momentum for free*/ | 
| 1481 | + | getInertiaTensor(intTensor,dummyAngMom); | 
| 1482 | + |  | 
| 1483 | + | det = intTensor.determinant(); | 
| 1484 | + | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; | 
| 1485 | + | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); | 
| 1486 | + | return; | 
| 1487 | + | } | 
| 1488 | + |  | 
| 1489 | + | void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ | 
| 1490 | + | Mat3x3d intTensor; | 
| 1491 | + | Vector3d dummyAngMom; | 
| 1492 | + | RealType sysconstants; | 
| 1493 | + | RealType geomCnst; | 
| 1494 | + |  | 
| 1495 | + | geomCnst = 3.0/2.0; | 
| 1496 | + | /* Get the inertial tensor and angular momentum for free*/ | 
| 1497 | + | getInertiaTensor(intTensor,dummyAngMom); | 
| 1498 | + |  | 
| 1499 | + | detI = intTensor.determinant(); | 
| 1500 | + | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; | 
| 1501 | + | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); | 
| 1502 | + | return; | 
| 1503 | + | } | 
| 1504 |  | /* | 
| 1505 |  | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { | 
| 1506 |  | assert( v.size() == nAtoms_ + nRigidBodies_); |