| 55 |  | #include "UseTheForce/fCutoffPolicy.h" | 
| 56 |  | #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" | 
| 57 |  | #include "UseTheForce/doForces_interface.h" | 
| 58 | + | #include "UseTheForce/DarkSide/electrostatic_interface.h" | 
| 59 |  | #include "UseTheForce/notifyCutoffs_interface.h" | 
| 60 |  | #include "utils/MemoryUtils.hpp" | 
| 61 |  | #include "utils/simError.h" | 
| 83 |  | MoleculeStamp* molStamp; | 
| 84 |  | int nMolWithSameStamp; | 
| 85 |  | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | 
| 86 | < | int nGroups = 0;          //total cutoff groups defined in meta-data file | 
| 86 | > | int nGroups = 0;      //total cutoff groups defined in meta-data file | 
| 87 |  | CutoffGroupStamp* cgStamp; | 
| 88 |  | RigidBodyStamp* rbStamp; | 
| 89 |  | int nRigidAtoms = 0; | 
| 108 |  | } | 
| 109 |  |  | 
| 110 |  | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; | 
| 111 | + |  | 
| 112 |  | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; | 
| 113 |  |  | 
| 114 |  | //calculate atoms in rigid bodies | 
| 125 |  |  | 
| 126 |  | } | 
| 127 |  |  | 
| 128 | < | //every free atom (atom does not belong to cutoff groups) is a cutoff group | 
| 129 | < | //therefore the total number of cutoff groups in the system is equal to | 
| 130 | < | //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data | 
| 131 | < | //file plus the number of cutoff groups defined in meta-data file | 
| 128 | > | //every free atom (atom does not belong to cutoff groups) is a cutoff | 
| 129 | > | //group therefore the total number of cutoff groups in the system is | 
| 130 | > | //equal to the total number of atoms minus number of atoms belong to | 
| 131 | > | //cutoff group defined in meta-data file plus the number of cutoff | 
| 132 | > | //groups defined in meta-data file | 
| 133 |  | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | 
| 134 |  |  | 
| 135 | < | //every free atom (atom does not belong to rigid bodies) is an integrable object | 
| 136 | < | //therefore the total number of  integrable objects in the system is equal to | 
| 137 | < | //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data | 
| 138 | < | //file plus the number of  rigid bodies defined in meta-data file | 
| 139 | < | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_; | 
| 140 | < |  | 
| 135 | > | //every free atom (atom does not belong to rigid bodies) is an | 
| 136 | > | //integrable object therefore the total number of integrable objects | 
| 137 | > | //in the system is equal to the total number of atoms minus number of | 
| 138 | > | //atoms belong to rigid body defined in meta-data file plus the number | 
| 139 | > | //of rigid bodies defined in meta-data file | 
| 140 | > | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms | 
| 141 | > | + nGlobalRigidBodies_; | 
| 142 | > |  | 
| 143 |  | nGlobalMols_ = molStampIds_.size(); | 
| 144 |  |  | 
| 145 |  | #ifdef IS_MPI | 
| 526 |  | int useDirectionalAtom = 0; | 
| 527 |  | int useElectrostatics = 0; | 
| 528 |  | //usePBC and useRF are from simParams | 
| 529 | < | int usePBC = simParams_->getPBC(); | 
| 529 | > | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | 
| 530 | > | int useRF; | 
| 531 | > | std::string myMethod; | 
| 532 | > |  | 
| 533 | > | // set the useRF logical | 
| 534 | > | useRF = 0; | 
| 535 | > | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 536 | > | myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 537 | > | if (myMethod == "REACTION_FIELD") | 
| 538 | > | useRF = 1; | 
| 539 | > | } | 
| 540 |  |  | 
| 541 |  | //loop over all of the atom types | 
| 542 |  | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 599 |  | temp = useFLARB; | 
| 600 |  | MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 601 |  |  | 
| 602 | + | temp = useRF; | 
| 603 | + | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 604 | + |  | 
| 605 |  | #endif | 
| 606 |  |  | 
| 607 |  | fInfo_.SIM_uses_PBC = usePBC; | 
| 616 |  | fInfo_.SIM_uses_EAM = useEAM; | 
| 617 |  | fInfo_.SIM_uses_Shapes = useShape; | 
| 618 |  | fInfo_.SIM_uses_FLARB = useFLARB; | 
| 619 | + | fInfo_.SIM_uses_RF = useRF; | 
| 620 |  |  | 
| 621 | < | if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { | 
| 621 | > | if( fInfo_.SIM_uses_Dipoles && myMethod == "REACTION_FIELD") { | 
| 622 |  |  | 
| 623 |  | if (simParams_->haveDielectric()) { | 
| 624 |  | fInfo_.dielect = simParams_->getDielectric(); | 
| 668 |  |  | 
| 669 |  | totalMass = cg->getMass(); | 
| 670 |  | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | 
| 671 | < | mfact.push_back(atom->getMass()/totalMass); | 
| 671 | > | // Check for massless groups - set mfact to 1 if true | 
| 672 | > | if (totalMass != 0) | 
| 673 | > | mfact.push_back(atom->getMass()/totalMass); | 
| 674 | > | else | 
| 675 | > | mfact.push_back( 1.0 ); | 
| 676 |  | } | 
| 677 |  |  | 
| 678 |  | } | 
| 808 |  |  | 
| 809 |  | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 810 |  |  | 
| 811 | < | if (!simParams_->haveRcut()){ | 
| 811 | > | if (!simParams_->haveCutoffRadius()){ | 
| 812 |  | sprintf(painCave.errMsg, | 
| 813 |  | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 814 |  | "\tOOPSE will use a default value of 15.0 angstroms" | 
| 817 |  | simError(); | 
| 818 |  | rcut = 15.0; | 
| 819 |  | } else{ | 
| 820 | < | rcut = simParams_->getRcut(); | 
| 820 | > | rcut = simParams_->getCutoffRadius(); | 
| 821 |  | } | 
| 822 |  |  | 
| 823 | < | if (!simParams_->haveRsw()){ | 
| 823 | > | if (!simParams_->haveSwitchingRadius()){ | 
| 824 |  | sprintf(painCave.errMsg, | 
| 825 |  | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 826 |  | "\tOOPSE will use a default value of\n" | 
| 829 |  | simError(); | 
| 830 |  | rsw = 0.95 * rcut; | 
| 831 |  | } else{ | 
| 832 | < | rsw = simParams_->getRsw(); | 
| 832 | > | rsw = simParams_->getSwitchingRadius(); | 
| 833 |  | } | 
| 834 |  |  | 
| 835 |  | } else { | 
| 836 |  | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in | 
| 837 |  | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used | 
| 838 |  |  | 
| 839 | < | if (simParams_->haveRcut()) { | 
| 840 | < | rcut = simParams_->getRcut(); | 
| 839 | > | if (simParams_->haveCutoffRadius()) { | 
| 840 | > | rcut = simParams_->getCutoffRadius(); | 
| 841 |  | } else { | 
| 842 |  | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system | 
| 843 |  | rcut = calcMaxCutoffRadius(); | 
| 844 |  | } | 
| 845 |  |  | 
| 846 | < | if (simParams_->haveRsw()) { | 
| 847 | < | rsw  = simParams_->getRsw(); | 
| 846 | > | if (simParams_->haveSwitchingRadius()) { | 
| 847 | > | rsw  = simParams_->getSwitchingRadius(); | 
| 848 |  | } else { | 
| 849 |  | rsw = rcut; | 
| 850 |  | } | 
| 861 |  | int cp =  TRADITIONAL_CUTOFF_POLICY; | 
| 862 |  | if (simParams_->haveCutoffPolicy()) { | 
| 863 |  | std::string myPolicy = simParams_->getCutoffPolicy(); | 
| 864 | + | toUpper(myPolicy); | 
| 865 |  | if (myPolicy == "MIX") { | 
| 866 |  | cp = MIX_CUTOFF_POLICY; | 
| 867 |  | } else { | 
| 880 |  | } | 
| 881 |  | } | 
| 882 |  | } | 
| 883 | + |  | 
| 884 | + |  | 
| 885 | + | if (simParams_->haveSkinThickness()) { | 
| 886 | + | double skinThickness = simParams_->getSkinThickness(); | 
| 887 | + | } | 
| 888 | + |  | 
| 889 |  | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); | 
| 890 | + | // also send cutoff notification to electrostatics | 
| 891 | + | setElectrostaticCutoffRadius(&rcut_, &rsw_); | 
| 892 |  | } | 
| 893 |  |  | 
| 894 |  | void SimInfo::setupElectrostaticSummationMethod( int isError ) { | 
| 896 |  | int errorOut; | 
| 897 |  | int esm =  NONE; | 
| 898 |  | double alphaVal; | 
| 899 | + | double dielectric; | 
| 900 |  |  | 
| 901 |  | errorOut = isError; | 
| 902 | + | alphaVal = simParams_->getDampingAlpha(); | 
| 903 | + | dielectric = simParams_->getDielectric(); | 
| 904 |  |  | 
| 905 |  | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 906 |  | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 907 | + | toUpper(myMethod); | 
| 908 |  | if (myMethod == "NONE") { | 
| 909 |  | esm = NONE; | 
| 910 |  | } else { | 
| 916 |  | if (!simParams_->haveDampingAlpha()) { | 
| 917 |  | //throw error | 
| 918 |  | sprintf( painCave.errMsg, | 
| 919 | < | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", simParams_->getDampingAlpha()); | 
| 919 | > | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", alphaVal); | 
| 920 |  | painCave.isFatal = 0; | 
| 921 |  | simError(); | 
| 922 |  | } | 
| 887 | – | alphaVal = simParams_->getDampingAlpha(); | 
| 923 |  | } else { | 
| 924 |  | if (myMethod == "REACTION_FIELD") { | 
| 925 |  | esm = REACTION_FIELD; | 
| 934 |  | } | 
| 935 |  | } | 
| 936 |  | } | 
| 937 | < | initFortranFF( &fInfo_.SIM_uses_RF, &esm, &alphaVal, &errorOut ); | 
| 937 | > | // let's pass some summation method variables to fortran | 
| 938 | > | setElectrostaticSummationMethod( &esm ); | 
| 939 | > | setDampedWolfAlpha( &alphaVal ); | 
| 940 | > | setReactionFieldDielectric( &dielectric ); | 
| 941 | > | initFortranFF( &esm, &errorOut ); | 
| 942 |  | } | 
| 943 |  |  | 
| 944 |  | void SimInfo::addProperty(GenericData* genData) { |