# | 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) { |
95 | > | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false), |
96 | > | useAtomicVirial_(true) { |
97 | ||
98 | MoleculeStamp* molStamp; | |
99 | int nMolWithSameStamp; | |
# | Line 599 | Line 603 | namespace oopse { | |
603 | //setup fortran force field | |
604 | /** @deprecate */ | |
605 | int isError = 0; | |
606 | + | |
607 | + | setupCutoff(); |
608 | ||
609 | setupElectrostaticSummationMethod( isError ); | |
610 | setupSwitchingFunction(); | |
611 | + | setupAccumulateBoxDipole(); |
612 | ||
613 | if(isError){ | |
614 | sprintf( painCave.errMsg, | |
# | Line 609 | Line 616 | namespace oopse { | |
616 | painCave.isFatal = 1; | |
617 | simError(); | |
618 | } | |
612 | – | |
613 | – | |
614 | – | setupCutoff(); |
619 | ||
620 | calcNdf(); | |
621 | calcNdfRaw(); | |
# | Line 661 | Line 665 | namespace oopse { | |
665 | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | |
666 | int useRF; | |
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; |
684 | < | } else { |
685 | < | if (myMethod == "SHIFTED_FORCE") { |
686 | < | useSF = 1; |
687 | < | } |
682 | > | if (myMethod == "REACTION_FIELD"){ |
683 | > | useRF = 1; |
684 | > | } else if (myMethod == "SHIFTED_FORCE"){ |
685 | > | useSF = 1; |
686 | > | } else if (myMethod == "SHIFTED_POTENTIAL"){ |
687 | > | useSP = 1; |
688 | } | |
689 | } | |
690 | + | |
691 | + | if (simParams_->haveAccumulateBoxDipole()) |
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 749 | Line 763 | namespace oopse { | |
763 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
764 | ||
765 | temp = useSF; | |
766 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
766 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
767 | ||
768 | + | temp = useSP; |
769 | + | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
770 | + | |
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 768 | Line 791 | namespace oopse { | |
791 | fInfo_.SIM_uses_FLARB = useFLARB; | |
792 | fInfo_.SIM_uses_RF = useRF; | |
793 | fInfo_.SIM_uses_SF = useSF; | |
794 | < | |
795 | < | if( myMethod == "REACTION_FIELD") { |
796 | < | |
774 | < | if (simParams_->haveDielectric()) { |
775 | < | fInfo_.dielect = simParams_->getDielectric(); |
776 | < | } else { |
777 | < | sprintf(painCave.errMsg, |
778 | < | "SimSetup Error: No Dielectric constant was set.\n" |
779 | < | "\tYou are trying to use Reaction Field without" |
780 | < | "\tsetting a dielectric constant!\n"); |
781 | < | painCave.isFatal = 1; |
782 | < | simError(); |
783 | < | } |
784 | < | } |
785 | < | |
794 | > | fInfo_.SIM_uses_SP = useSP; |
795 | > | fInfo_.SIM_uses_BoxDipole = useBoxDipole; |
796 | > | fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; |
797 | } | |
798 | ||
799 | void SimInfo::setupFortranSim() { | |
# | Line 867 | 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 1058 | Line 1077 | namespace oopse { | |
1077 | int sm = UNDAMPED; | |
1078 | RealType alphaVal; | |
1079 | RealType dielectric; | |
1080 | < | |
1080 | > | |
1081 | errorOut = isError; | |
1063 | – | alphaVal = simParams_->getDampingAlpha(); |
1064 | – | dielectric = simParams_->getDielectric(); |
1082 | ||
1083 | if (simParams_->haveElectrostaticSummationMethod()) { | |
1084 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
# | Line 1078 | 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 1106 | 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 1161 | Line 1196 | namespace oopse { | |
1196 | ||
1197 | } | |
1198 | ||
1199 | + | void SimInfo::setupAccumulateBoxDipole() { |
1200 | + | |
1201 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
1202 | + | if ( simParams_->haveAccumulateBoxDipole() ) |
1203 | + | if ( simParams_->getAccumulateBoxDipole() ) { |
1204 | + | setAccumulateBoxDipole(); |
1205 | + | calcBoxDipole_ = true; |
1206 | + | } |
1207 | + | |
1208 | + | } |
1209 | + | |
1210 | void SimInfo::addProperty(GenericData* genData) { | |
1211 | properties_.addProperty(genData); | |
1212 | } | |
# | Line 1410 | Line 1456 | namespace oopse { | |
1456 | return angularMomentum; | |
1457 | } | |
1458 | ||
1459 | < | |
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 |