# | Line 52 | Line 52 | |
---|---|---|
52 | #include "brains/SimInfo.hpp" | |
53 | #include "math/Vector3.hpp" | |
54 | #include "primitives/Molecule.hpp" | |
55 | + | #include "UseTheForce/fCutoffPolicy.h" |
56 | + | #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" |
57 | + | #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" |
58 | + | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" |
59 | #include "UseTheForce/doForces_interface.h" | |
60 | + | #include "UseTheForce/DarkSide/electrostatic_interface.h" |
61 | #include "UseTheForce/notifyCutoffs_interface.h" | |
62 | + | #include "UseTheForce/DarkSide/switcheroo_interface.h" |
63 | #include "utils/MemoryUtils.hpp" | |
64 | #include "utils/simError.h" | |
65 | #include "selection/SelectionManager.hpp" | |
# | Line 80 | Line 86 | namespace oopse { | |
86 | MoleculeStamp* molStamp; | |
87 | int nMolWithSameStamp; | |
88 | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | |
89 | < | int nGroups = 0; //total cutoff groups defined in meta-data file |
89 | > | int nGroups = 0; //total cutoff groups defined in meta-data file |
90 | CutoffGroupStamp* cgStamp; | |
91 | RigidBodyStamp* rbStamp; | |
92 | int nRigidAtoms = 0; | |
# | Line 105 | Line 111 | namespace oopse { | |
111 | } | |
112 | ||
113 | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; | |
114 | + | |
115 | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; | |
116 | ||
117 | //calculate atoms in rigid bodies | |
# | Line 121 | Line 128 | namespace oopse { | |
128 | ||
129 | } | |
130 | ||
131 | < | //every free atom (atom does not belong to cutoff groups) is a cutoff group |
132 | < | //therefore the total number of cutoff groups in the system is equal to |
133 | < | //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data |
134 | < | //file plus the number of cutoff groups defined in meta-data file |
131 | > | //every free atom (atom does not belong to cutoff groups) is a cutoff |
132 | > | //group therefore the total number of cutoff groups in the system is |
133 | > | //equal to the total number of atoms minus number of atoms belong to |
134 | > | //cutoff group defined in meta-data file plus the number of cutoff |
135 | > | //groups defined in meta-data file |
136 | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | |
137 | ||
138 | < | //every free atom (atom does not belong to rigid bodies) is an integrable object |
139 | < | //therefore the total number of integrable objects in the system is equal to |
140 | < | //the total number of atoms minus number of atoms belong to rigid body defined in meta-data |
141 | < | //file plus the number of rigid bodies defined in meta-data file |
142 | < | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_; |
143 | < | |
138 | > | //every free atom (atom does not belong to rigid bodies) is an |
139 | > | //integrable object therefore the total number of integrable objects |
140 | > | //in the system is equal to the total number of atoms minus number of |
141 | > | //atoms belong to rigid body defined in meta-data file plus the number |
142 | > | //of rigid bodies defined in meta-data file |
143 | > | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms |
144 | > | + nGlobalRigidBodies_; |
145 | > | |
146 | nGlobalMols_ = molStampIds_.size(); | |
147 | ||
148 | #ifdef IS_MPI | |
# | Line 462 | Line 472 | namespace oopse { | |
472 | //setup fortran force field | |
473 | /** @deprecate */ | |
474 | int isError = 0; | |
475 | < | initFortranFF( &fInfo_.SIM_uses_RF , &isError ); |
475 | > | |
476 | > | setupElectrostaticSummationMethod( isError ); |
477 | > | setupSwitchingFunction(); |
478 | > | |
479 | if(isError){ | |
480 | sprintf( painCave.errMsg, | |
481 | "ForceField error: There was an error initializing the forceField in fortran.\n" ); | |
# | Line 517 | Line 530 | namespace oopse { | |
530 | int useDirectionalAtom = 0; | |
531 | int useElectrostatics = 0; | |
532 | //usePBC and useRF are from simParams | |
533 | < | int usePBC = simParams_->getPBC(); |
534 | < | int useRF = simParams_->getUseRF(); |
533 | > | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); |
534 | > | int useRF; |
535 | > | int useSF; |
536 | > | std::string myMethod; |
537 | > | |
538 | > | // set the useRF logical |
539 | > | useRF = 0; |
540 | > | useSF = 0; |
541 | > | |
542 | > | |
543 | > | if (simParams_->haveElectrostaticSummationMethod()) { |
544 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
545 | > | toUpper(myMethod); |
546 | > | if (myMethod == "REACTION_FIELD") { |
547 | > | useRF=1; |
548 | > | } else { |
549 | > | if (myMethod == "SHIFTED_FORCE") { |
550 | > | useSF = 1; |
551 | > | } |
552 | > | } |
553 | > | } |
554 | ||
555 | //loop over all of the atom types | |
556 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
# | Line 583 | Line 615 | namespace oopse { | |
615 | ||
616 | temp = useRF; | |
617 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
618 | < | |
618 | > | |
619 | > | temp = useSF; |
620 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
621 | > | |
622 | #endif | |
623 | ||
624 | fInfo_.SIM_uses_PBC = usePBC; | |
# | Line 599 | Line 634 | namespace oopse { | |
634 | fInfo_.SIM_uses_Shapes = useShape; | |
635 | fInfo_.SIM_uses_FLARB = useFLARB; | |
636 | fInfo_.SIM_uses_RF = useRF; | |
637 | + | fInfo_.SIM_uses_SF = useSF; |
638 | ||
639 | < | if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { |
640 | < | |
639 | > | if( myMethod == "REACTION_FIELD") { |
640 | > | |
641 | if (simParams_->haveDielectric()) { | |
642 | fInfo_.dielect = simParams_->getDielectric(); | |
643 | } else { | |
# | Line 611 | Line 647 | namespace oopse { | |
647 | "\tsetting a dielectric constant!\n"); | |
648 | painCave.isFatal = 1; | |
649 | simError(); | |
650 | < | } |
615 | < | |
616 | < | } else { |
617 | < | fInfo_.dielect = 0.0; |
650 | > | } |
651 | } | |
652 | ||
653 | } | |
# | Line 650 | Line 683 | namespace oopse { | |
683 | ||
684 | totalMass = cg->getMass(); | |
685 | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | |
686 | < | mfact.push_back(atom->getMass()/totalMass); |
686 | > | // Check for massless groups - set mfact to 1 if true |
687 | > | if (totalMass != 0) |
688 | > | mfact.push_back(atom->getMass()/totalMass); |
689 | > | else |
690 | > | mfact.push_back( 1.0 ); |
691 | } | |
692 | ||
693 | } | |
# | Line 786 | Line 823 | namespace oopse { | |
823 | ||
824 | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | |
825 | ||
826 | < | if (!simParams_->haveRcut()){ |
826 | > | if (!simParams_->haveCutoffRadius()){ |
827 | sprintf(painCave.errMsg, | |
828 | "SimCreator Warning: No value was set for the cutoffRadius.\n" | |
829 | "\tOOPSE will use a default value of 15.0 angstroms" | |
# | Line 795 | Line 832 | namespace oopse { | |
832 | simError(); | |
833 | rcut = 15.0; | |
834 | } else{ | |
835 | < | rcut = simParams_->getRcut(); |
835 | > | rcut = simParams_->getCutoffRadius(); |
836 | } | |
837 | ||
838 | < | if (!simParams_->haveRsw()){ |
838 | > | if (!simParams_->haveSwitchingRadius()){ |
839 | sprintf(painCave.errMsg, | |
840 | "SimCreator Warning: No value was set for switchingRadius.\n" | |
841 | "\tOOPSE will use a default value of\n" | |
842 | < | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
842 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
843 | painCave.isFatal = 0; | |
844 | simError(); | |
845 | < | rsw = 0.95 * rcut; |
845 | > | rsw = 0.85 * rcut; |
846 | } else{ | |
847 | < | rsw = simParams_->getRsw(); |
847 | > | rsw = simParams_->getSwitchingRadius(); |
848 | } | |
849 | ||
850 | } else { | |
851 | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in | |
852 | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used | |
853 | ||
854 | < | if (simParams_->haveRcut()) { |
855 | < | rcut = simParams_->getRcut(); |
854 | > | if (simParams_->haveCutoffRadius()) { |
855 | > | rcut = simParams_->getCutoffRadius(); |
856 | } else { | |
857 | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system | |
858 | rcut = calcMaxCutoffRadius(); | |
859 | } | |
860 | ||
861 | < | if (simParams_->haveRsw()) { |
862 | < | rsw = simParams_->getRsw(); |
861 | > | if (simParams_->haveSwitchingRadius()) { |
862 | > | rsw = simParams_->getSwitchingRadius(); |
863 | } else { | |
864 | rsw = rcut; | |
865 | } | |
# | Line 830 | Line 867 | namespace oopse { | |
867 | } | |
868 | } | |
869 | ||
870 | < | void SimInfo::setupCutoff() { |
870 | > | void SimInfo::setupCutoff() { |
871 | getCutoff(rcut_, rsw_); | |
872 | double rnblist = rcut_ + 1; // skin of neighbor list | |
873 | ||
874 | //Pass these cutoff radius etc. to fortran. This function should be called once and only once | |
875 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist); |
875 | > | |
876 | > | int cp = TRADITIONAL_CUTOFF_POLICY; |
877 | > | if (simParams_->haveCutoffPolicy()) { |
878 | > | std::string myPolicy = simParams_->getCutoffPolicy(); |
879 | > | toUpper(myPolicy); |
880 | > | if (myPolicy == "MIX") { |
881 | > | cp = MIX_CUTOFF_POLICY; |
882 | > | } else { |
883 | > | if (myPolicy == "MAX") { |
884 | > | cp = MAX_CUTOFF_POLICY; |
885 | > | } else { |
886 | > | if (myPolicy == "TRADITIONAL") { |
887 | > | cp = TRADITIONAL_CUTOFF_POLICY; |
888 | > | } else { |
889 | > | // throw error |
890 | > | sprintf( painCave.errMsg, |
891 | > | "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); |
892 | > | painCave.isFatal = 1; |
893 | > | simError(); |
894 | > | } |
895 | > | } |
896 | > | } |
897 | > | } |
898 | > | |
899 | > | |
900 | > | if (simParams_->haveSkinThickness()) { |
901 | > | double skinThickness = simParams_->getSkinThickness(); |
902 | > | } |
903 | > | |
904 | > | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); |
905 | > | // also send cutoff notification to electrostatics |
906 | > | setElectrostaticCutoffRadius(&rcut_, &rsw_); |
907 | > | } |
908 | > | |
909 | > | void SimInfo::setupElectrostaticSummationMethod( int isError ) { |
910 | > | |
911 | > | int errorOut; |
912 | > | int esm = NONE; |
913 | > | int sm = UNDAMPED; |
914 | > | double alphaVal; |
915 | > | double dielectric; |
916 | > | |
917 | > | errorOut = isError; |
918 | > | alphaVal = simParams_->getDampingAlpha(); |
919 | > | dielectric = simParams_->getDielectric(); |
920 | > | |
921 | > | if (simParams_->haveElectrostaticSummationMethod()) { |
922 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
923 | > | toUpper(myMethod); |
924 | > | if (myMethod == "NONE") { |
925 | > | esm = NONE; |
926 | > | } else { |
927 | > | if (myMethod == "SWITCHING_FUNCTION") { |
928 | > | esm = SWITCHING_FUNCTION; |
929 | > | } else { |
930 | > | if (myMethod == "SHIFTED_POTENTIAL") { |
931 | > | esm = SHIFTED_POTENTIAL; |
932 | > | } else { |
933 | > | if (myMethod == "SHIFTED_FORCE") { |
934 | > | esm = SHIFTED_FORCE; |
935 | > | } else { |
936 | > | if (myMethod == "REACTION_FIELD") { |
937 | > | esm = REACTION_FIELD; |
938 | > | } else { |
939 | > | // throw error |
940 | > | sprintf( painCave.errMsg, |
941 | > | "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() ); |
942 | > | painCave.isFatal = 1; |
943 | > | simError(); |
944 | > | } |
945 | > | } |
946 | > | } |
947 | > | } |
948 | > | } |
949 | > | } |
950 | > | |
951 | > | if (simParams_->haveElectrostaticScreeningMethod()) { |
952 | > | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); |
953 | > | toUpper(myScreen); |
954 | > | if (myScreen == "UNDAMPED") { |
955 | > | sm = UNDAMPED; |
956 | > | } else { |
957 | > | if (myScreen == "DAMPED") { |
958 | > | sm = DAMPED; |
959 | > | if (!simParams_->haveDampingAlpha()) { |
960 | > | //throw error |
961 | > | sprintf( painCave.errMsg, |
962 | > | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal); |
963 | > | painCave.isFatal = 0; |
964 | > | simError(); |
965 | > | } |
966 | > | } else { |
967 | > | // throw error |
968 | > | sprintf( painCave.errMsg, |
969 | > | "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); |
970 | > | painCave.isFatal = 1; |
971 | > | simError(); |
972 | > | } |
973 | > | } |
974 | > | } |
975 | > | |
976 | > | // let's pass some summation method variables to fortran |
977 | > | setElectrostaticSummationMethod( &esm ); |
978 | > | setScreeningMethod( &sm ); |
979 | > | setDampingAlpha( &alphaVal ); |
980 | > | setReactionFieldDielectric( &dielectric ); |
981 | > | initFortranFF( &esm, &errorOut ); |
982 | > | } |
983 | > | |
984 | > | void SimInfo::setupSwitchingFunction() { |
985 | > | int ft = CUBIC; |
986 | > | |
987 | > | if (simParams_->haveSwitchingFunctionType()) { |
988 | > | std::string funcType = simParams_->getSwitchingFunctionType(); |
989 | > | toUpper(funcType); |
990 | > | if (funcType == "CUBIC") { |
991 | > | ft = CUBIC; |
992 | > | } else { |
993 | > | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
994 | > | ft = FIFTH_ORDER_POLY; |
995 | > | } else { |
996 | > | // throw error |
997 | > | sprintf( painCave.errMsg, |
998 | > | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); |
999 | > | painCave.isFatal = 1; |
1000 | > | simError(); |
1001 | > | } |
1002 | > | } |
1003 | > | } |
1004 | > | |
1005 | > | // send switching function notification to switcheroo |
1006 | > | setFunctionType(&ft); |
1007 | > | |
1008 | } | |
1009 | ||
1010 | void SimInfo::addProperty(GenericData* genData) { |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |