# | Line 53 | Line 53 | |
---|---|---|
53 | #include "brains/SimInfo.hpp" | |
54 | #include "math/Vector3.hpp" | |
55 | #include "primitives/Molecule.hpp" | |
56 | + | #include "primitives/StuntDouble.hpp" |
57 | #include "UseTheForce/fCutoffPolicy.h" | |
58 | #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" | |
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" | |
# | Line 66 | Line 68 | |
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" | |
# | Line 89 | Line 92 | namespace oopse { | |
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; | |
# | Line 600 | Line 604 | namespace oopse { | |
604 | /** @deprecate */ | |
605 | int isError = 0; | |
606 | ||
607 | + | setupCutoff(); |
608 | + | |
609 | setupElectrostaticSummationMethod( isError ); | |
610 | setupSwitchingFunction(); | |
611 | setupAccumulateBoxDipole(); | |
# | Line 610 | Line 616 | namespace oopse { | |
616 | painCave.isFatal = 1; | |
617 | simError(); | |
618 | } | |
613 | – | |
614 | – | |
615 | – | setupCutoff(); |
619 | ||
620 | calcNdf(); | |
621 | calcNdfRaw(); | |
# | Line 664 | Line 667 | namespace oopse { | |
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"){ | |
# | Line 687 | Line 692 | namespace oopse { | |
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(); | |
# | Line 764 | Line 771 | namespace oopse { | |
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; | |
# | Line 783 | Line 793 | namespace oopse { | |
793 | fInfo_.SIM_uses_SF = useSF; | |
794 | fInfo_.SIM_uses_SP = useSP; | |
795 | fInfo_.SIM_uses_BoxDipole = useBoxDipole; | |
796 | < | |
787 | < | if( myMethod == "REACTION_FIELD") { |
788 | < | |
789 | < | if (simParams_->haveDielectric()) { |
790 | < | fInfo_.dielect = simParams_->getDielectric(); |
791 | < | } else { |
792 | < | sprintf(painCave.errMsg, |
793 | < | "SimSetup Error: No Dielectric constant was set.\n" |
794 | < | "\tYou are trying to use Reaction Field without" |
795 | < | "\tsetting a dielectric constant!\n"); |
796 | < | painCave.isFatal = 1; |
797 | < | simError(); |
798 | < | } |
799 | < | } |
800 | < | |
796 | > | fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; |
797 | } | |
798 | ||
799 | void SimInfo::setupFortranSim() { | |
# | Line 882 | Line 878 | namespace oopse { | |
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 | ||
# | Line 1073 | Line 1077 | namespace oopse { | |
1077 | int sm = UNDAMPED; | |
1078 | RealType alphaVal; | |
1079 | RealType dielectric; | |
1080 | < | |
1080 | > | |
1081 | errorOut = isError; | |
1078 | – | alphaVal = simParams_->getDampingAlpha(); |
1079 | – | dielectric = simParams_->getDielectric(); |
1082 | ||
1083 | if (simParams_->haveElectrostaticSummationMethod()) { | |
1084 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
# | Line 1093 | Line 1095 | namespace oopse { | |
1095 | if (myMethod == "SHIFTED_FORCE") { | |
1096 | esm = SHIFTED_FORCE; | |
1097 | } else { | |
1098 | < | if (myMethod == "REACTION_FIELD") { |
1098 | > | if (myMethod == "REACTION_FIELD") { |
1099 | esm = REACTION_FIELD; | |
1100 | + | dielectric = simParams_->getDielectric(); |
1101 | + | if (!simParams_->haveDielectric()) { |
1102 | + | // throw warning |
1103 | + | sprintf( painCave.errMsg, |
1104 | + | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" |
1105 | + | "\tA default value of %f will be used for the dielectric.\n", dielectric); |
1106 | + | painCave.isFatal = 0; |
1107 | + | simError(); |
1108 | + | } |
1109 | } else { | |
1110 | // throw error | |
1111 | sprintf( painCave.errMsg, | |
# | Line 1121 | Line 1132 | namespace oopse { | |
1132 | if (myScreen == "DAMPED") { | |
1133 | sm = DAMPED; | |
1134 | if (!simParams_->haveDampingAlpha()) { | |
1135 | < | //throw error |
1135 | > | // first set a cutoff dependent alpha value |
1136 | > | // we assume alpha depends linearly with rcut from 0 to 20.5 ang |
1137 | > | alphaVal = 0.5125 - rcut_* 0.025; |
1138 | > | // for values rcut > 20.5, alpha is zero |
1139 | > | if (alphaVal < 0) alphaVal = 0; |
1140 | > | |
1141 | > | // throw warning |
1142 | sprintf( painCave.errMsg, | |
1143 | "SimInfo warning: dampingAlpha was not specified in the input file.\n" | |
1144 | < | "\tA default value of %f (1/ang) will be used.\n", alphaVal); |
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, | |
# | Line 1435 | Line 1455 | namespace oopse { | |
1455 | ||
1456 | return angularMomentum; | |
1457 | } | |
1438 | – | |
1458 | ||
1459 | + | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { |
1460 | + | return IOIndexToIntegrableObject.at(index); |
1461 | + | } |
1462 | + | |
1463 | + | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { |
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_); |
1507 | + | sdByGlobalIndex_ = v; |
1508 | + | } |
1509 | + | |
1510 | + | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { |
1511 | + | //assert(index < nAtoms_ + nRigidBodies_); |
1512 | + | return sdByGlobalIndex_.at(index); |
1513 | + | } |
1514 | + | */ |
1515 | }//end namespace oopse | |
1516 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |