# | Line 53 | Line 53 | |
---|---|---|
53 | #include "math/Vector3.hpp" | |
54 | #include "primitives/Molecule.hpp" | |
55 | #include "UseTheForce/fCutoffPolicy.h" | |
56 | < | #include "UseTheForce/fCoulombicCorrection.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 82 | 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 107 | Line 111 | namespace oopse { | |
111 | } | |
112 | ||
113 | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; | |
114 | + | |
115 | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; | |
116 | ||
117 | //calculate atoms in rigid bodies | |
# | Line 123 | 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 465 | Line 473 | namespace oopse { | |
473 | /** @deprecate */ | |
474 | int isError = 0; | |
475 | ||
476 | < | setupCoulombicCorrection( isError ); |
476 | > | setupElectrostaticSummationMethod( isError ); |
477 | > | setupSwitchingFunction(); |
478 | ||
479 | if(isError){ | |
480 | sprintf( painCave.errMsg, | |
# | Line 510 | Line 519 | namespace oopse { | |
519 | int useLennardJones = 0; | |
520 | int useElectrostatic = 0; | |
521 | int useEAM = 0; | |
522 | + | int useSC = 0; |
523 | int useCharge = 0; | |
524 | int useDirectional = 0; | |
525 | int useDipole = 0; | |
# | Line 521 | Line 531 | namespace oopse { | |
531 | int useDirectionalAtom = 0; | |
532 | int useElectrostatics = 0; | |
533 | //usePBC and useRF are from simParams | |
534 | < | int usePBC = simParams_->getPBC(); |
535 | < | int useRF = simParams_->getUseRF(); |
534 | > | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); |
535 | > | int useRF; |
536 | > | int useSF; |
537 | > | std::string myMethod; |
538 | > | |
539 | > | // set the useRF logical |
540 | > | useRF = 0; |
541 | > | useSF = 0; |
542 | > | |
543 | ||
544 | + | if (simParams_->haveElectrostaticSummationMethod()) { |
545 | + | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
546 | + | toUpper(myMethod); |
547 | + | if (myMethod == "REACTION_FIELD") { |
548 | + | useRF=1; |
549 | + | } else { |
550 | + | if (myMethod == "SHIFTED_FORCE") { |
551 | + | useSF = 1; |
552 | + | } |
553 | + | } |
554 | + | } |
555 | + | |
556 | //loop over all of the atom types | |
557 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
558 | useLennardJones |= (*i)->isLennardJones(); | |
559 | useElectrostatic |= (*i)->isElectrostatic(); | |
560 | useEAM |= (*i)->isEAM(); | |
561 | + | useSC |= (*i)->isSC(); |
562 | useCharge |= (*i)->isCharge(); | |
563 | useDirectional |= (*i)->isDirectional(); | |
564 | useDipole |= (*i)->isDipole(); | |
# | Line 579 | Line 609 | namespace oopse { | |
609 | temp = useEAM; | |
610 | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
611 | ||
612 | + | temp = useSC; |
613 | + | MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
614 | + | |
615 | temp = useShape; | |
616 | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
617 | ||
# | Line 588 | Line 621 | namespace oopse { | |
621 | temp = useRF; | |
622 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
623 | ||
624 | < | temp = useUW; |
625 | < | MPI_Allreduce(&temp, &useUW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
624 | > | temp = useSF; |
625 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
626 | ||
594 | – | temp = useDW; |
595 | – | MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
596 | – | |
627 | #endif | |
628 | ||
629 | fInfo_.SIM_uses_PBC = usePBC; | |
# | Line 606 | Line 636 | namespace oopse { | |
636 | fInfo_.SIM_uses_StickyPower = useStickyPower; | |
637 | fInfo_.SIM_uses_GayBerne = useGayBerne; | |
638 | fInfo_.SIM_uses_EAM = useEAM; | |
639 | + | fInfo_.SIM_uses_SC = useSC; |
640 | fInfo_.SIM_uses_Shapes = useShape; | |
641 | fInfo_.SIM_uses_FLARB = useFLARB; | |
642 | fInfo_.SIM_uses_RF = useRF; | |
643 | + | fInfo_.SIM_uses_SF = useSF; |
644 | ||
645 | < | if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { |
646 | < | |
645 | > | if( myMethod == "REACTION_FIELD") { |
646 | > | |
647 | if (simParams_->haveDielectric()) { | |
648 | fInfo_.dielect = simParams_->getDielectric(); | |
649 | } else { | |
# | Line 621 | Line 653 | namespace oopse { | |
653 | "\tsetting a dielectric constant!\n"); | |
654 | painCave.isFatal = 1; | |
655 | simError(); | |
656 | < | } |
625 | < | |
626 | < | } else { |
627 | < | fInfo_.dielect = 0.0; |
656 | > | } |
657 | } | |
658 | ||
659 | } | |
# | Line 660 | Line 689 | namespace oopse { | |
689 | ||
690 | totalMass = cg->getMass(); | |
691 | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | |
692 | < | mfact.push_back(atom->getMass()/totalMass); |
692 | > | // Check for massless groups - set mfact to 1 if true |
693 | > | if (totalMass != 0) |
694 | > | mfact.push_back(atom->getMass()/totalMass); |
695 | > | else |
696 | > | mfact.push_back( 1.0 ); |
697 | } | |
698 | ||
699 | } | |
# | Line 796 | Line 829 | namespace oopse { | |
829 | ||
830 | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | |
831 | ||
832 | < | if (!simParams_->haveRcut()){ |
832 | > | if (!simParams_->haveCutoffRadius()){ |
833 | sprintf(painCave.errMsg, | |
834 | "SimCreator Warning: No value was set for the cutoffRadius.\n" | |
835 | "\tOOPSE will use a default value of 15.0 angstroms" | |
# | Line 805 | Line 838 | namespace oopse { | |
838 | simError(); | |
839 | rcut = 15.0; | |
840 | } else{ | |
841 | < | rcut = simParams_->getRcut(); |
841 | > | rcut = simParams_->getCutoffRadius(); |
842 | } | |
843 | ||
844 | < | if (!simParams_->haveRsw()){ |
844 | > | if (!simParams_->haveSwitchingRadius()){ |
845 | sprintf(painCave.errMsg, | |
846 | "SimCreator Warning: No value was set for switchingRadius.\n" | |
847 | "\tOOPSE will use a default value of\n" | |
848 | < | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
848 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
849 | painCave.isFatal = 0; | |
850 | simError(); | |
851 | < | rsw = 0.95 * rcut; |
851 | > | rsw = 0.85 * rcut; |
852 | } else{ | |
853 | < | rsw = simParams_->getRsw(); |
853 | > | rsw = simParams_->getSwitchingRadius(); |
854 | } | |
855 | ||
856 | } else { | |
857 | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in | |
858 | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used | |
859 | ||
860 | < | if (simParams_->haveRcut()) { |
861 | < | rcut = simParams_->getRcut(); |
860 | > | if (simParams_->haveCutoffRadius()) { |
861 | > | rcut = simParams_->getCutoffRadius(); |
862 | } else { | |
863 | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system | |
864 | rcut = calcMaxCutoffRadius(); | |
865 | } | |
866 | ||
867 | < | if (simParams_->haveRsw()) { |
868 | < | rsw = simParams_->getRsw(); |
867 | > | if (simParams_->haveSwitchingRadius()) { |
868 | > | rsw = simParams_->getSwitchingRadius(); |
869 | } else { | |
870 | rsw = rcut; | |
871 | } | |
# | Line 849 | Line 882 | namespace oopse { | |
882 | int cp = TRADITIONAL_CUTOFF_POLICY; | |
883 | if (simParams_->haveCutoffPolicy()) { | |
884 | std::string myPolicy = simParams_->getCutoffPolicy(); | |
885 | + | toUpper(myPolicy); |
886 | if (myPolicy == "MIX") { | |
887 | cp = MIX_CUTOFF_POLICY; | |
888 | } else { | |
# | Line 867 | Line 901 | namespace oopse { | |
901 | } | |
902 | } | |
903 | } | |
904 | + | |
905 | + | |
906 | + | if (simParams_->haveSkinThickness()) { |
907 | + | double skinThickness = simParams_->getSkinThickness(); |
908 | + | } |
909 | + | |
910 | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); | |
911 | + | // also send cutoff notification to electrostatics |
912 | + | setElectrostaticCutoffRadius(&rcut_, &rsw_); |
913 | } | |
914 | ||
915 | < | void SimInfo::setupCoulombicCorrection( int isError ) { |
915 | > | void SimInfo::setupElectrostaticSummationMethod( int isError ) { |
916 | ||
917 | int errorOut; | |
918 | < | int cc = NONE; |
918 | > | int esm = NONE; |
919 | > | int sm = UNDAMPED; |
920 | double alphaVal; | |
921 | + | double dielectric; |
922 | ||
923 | errorOut = isError; | |
924 | + | alphaVal = simParams_->getDampingAlpha(); |
925 | + | dielectric = simParams_->getDielectric(); |
926 | ||
927 | < | if (simParams_->haveCoulombicCorrection()) { |
928 | < | std::string myCorrection = simParams_->getCoulombicCorrection(); |
929 | < | if (myCorrection == "NONE") { |
930 | < | cc = NONE; |
927 | > | if (simParams_->haveElectrostaticSummationMethod()) { |
928 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
929 | > | toUpper(myMethod); |
930 | > | if (myMethod == "NONE") { |
931 | > | esm = NONE; |
932 | } else { | |
933 | < | if (myCorrection == "UNDAMPED_WOLF") { |
934 | < | cc = UNDAMPED_WOLF; |
933 | > | if (myMethod == "SWITCHING_FUNCTION") { |
934 | > | esm = SWITCHING_FUNCTION; |
935 | } else { | |
936 | < | if (myCorrection == "WOLF") { |
937 | < | cc = WOLF; |
938 | < | if (!simParams_->haveDampingAlpha()) { |
939 | < | //throw error |
940 | < | sprintf( painCave.errMsg, |
894 | < | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Wolf Coulombic Correction.", simParams_->getDampingAlpha()); |
895 | < | painCave.isFatal = 0; |
896 | < | simError(); |
897 | < | } |
898 | < | alphaVal = simParams_->getDampingAlpha(); |
899 | < | } else { |
900 | < | if (myCorrection == "REACTION_FIELD") { |
901 | < | cc = REACTION_FIELD; |
936 | > | if (myMethod == "SHIFTED_POTENTIAL") { |
937 | > | esm = SHIFTED_POTENTIAL; |
938 | > | } else { |
939 | > | if (myMethod == "SHIFTED_FORCE") { |
940 | > | esm = SHIFTED_FORCE; |
941 | } else { | |
942 | < | // throw error |
943 | < | sprintf( painCave.errMsg, |
944 | < | "SimInfo error: Unknown coulombicCorrection. (Input file specified %s .)\n\tcoulombicCorrection must be one of: \"none\", \"undamped_wolf\", \"wolf\", or \"reaction_field\".", myCorrection.c_str() ); |
945 | < | painCave.isFatal = 1; |
946 | < | simError(); |
947 | < | } |
948 | < | } |
942 | > | if (myMethod == "REACTION_FIELD") { |
943 | > | esm = REACTION_FIELD; |
944 | > | } else { |
945 | > | // throw error |
946 | > | sprintf( painCave.errMsg, |
947 | > | "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() ); |
948 | > | painCave.isFatal = 1; |
949 | > | simError(); |
950 | > | } |
951 | > | } |
952 | > | } |
953 | } | |
954 | } | |
955 | } | |
956 | < | initFortranFF( &fInfo_.SIM_uses_RF, &cc, &alphaVal, &errorOut ); |
956 | > | |
957 | > | if (simParams_->haveElectrostaticScreeningMethod()) { |
958 | > | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); |
959 | > | toUpper(myScreen); |
960 | > | if (myScreen == "UNDAMPED") { |
961 | > | sm = UNDAMPED; |
962 | > | } else { |
963 | > | if (myScreen == "DAMPED") { |
964 | > | sm = DAMPED; |
965 | > | if (!simParams_->haveDampingAlpha()) { |
966 | > | //throw error |
967 | > | sprintf( painCave.errMsg, |
968 | > | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal); |
969 | > | painCave.isFatal = 0; |
970 | > | simError(); |
971 | > | } |
972 | > | } else { |
973 | > | // throw error |
974 | > | sprintf( painCave.errMsg, |
975 | > | "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); |
976 | > | painCave.isFatal = 1; |
977 | > | simError(); |
978 | > | } |
979 | > | } |
980 | > | } |
981 | > | |
982 | > | // let's pass some summation method variables to fortran |
983 | > | setElectrostaticSummationMethod( &esm ); |
984 | > | setScreeningMethod( &sm ); |
985 | > | setDampingAlpha( &alphaVal ); |
986 | > | setReactionFieldDielectric( &dielectric ); |
987 | > | initFortranFF( &esm, &errorOut ); |
988 | } | |
989 | ||
990 | + | void SimInfo::setupSwitchingFunction() { |
991 | + | int ft = CUBIC; |
992 | + | |
993 | + | if (simParams_->haveSwitchingFunctionType()) { |
994 | + | std::string funcType = simParams_->getSwitchingFunctionType(); |
995 | + | toUpper(funcType); |
996 | + | if (funcType == "CUBIC") { |
997 | + | ft = CUBIC; |
998 | + | } else { |
999 | + | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
1000 | + | ft = FIFTH_ORDER_POLY; |
1001 | + | } else { |
1002 | + | // throw error |
1003 | + | sprintf( painCave.errMsg, |
1004 | + | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); |
1005 | + | painCave.isFatal = 1; |
1006 | + | simError(); |
1007 | + | } |
1008 | + | } |
1009 | + | } |
1010 | + | |
1011 | + | // send switching function notification to switcheroo |
1012 | + | setFunctionType(&ft); |
1013 | + | |
1014 | + | } |
1015 | + | |
1016 | void SimInfo::addProperty(GenericData* genData) { | |
1017 | properties_.addProperty(genData); | |
1018 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |