| 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" | 
| 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" | 
| 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 | 
  | 
 | 
| 97 | 
  | 
      MoleculeStamp* molStamp; | 
| 98 | 
  | 
      int nMolWithSameStamp; | 
| 602 | 
  | 
    //setup fortran force field | 
| 603 | 
  | 
    /** @deprecate */     | 
| 604 | 
  | 
    int isError = 0; | 
| 605 | 
+ | 
     | 
| 606 | 
+ | 
    setupCutoff(); | 
| 607 | 
  | 
     | 
| 608 | 
  | 
    setupElectrostaticSummationMethod( isError ); | 
| 609 | 
  | 
    setupSwitchingFunction(); | 
| 610 | 
+ | 
    setupAccumulateBoxDipole(); | 
| 611 | 
  | 
 | 
| 612 | 
  | 
    if(isError){ | 
| 613 | 
  | 
      sprintf( painCave.errMsg, | 
| 615 | 
  | 
      painCave.isFatal = 1; | 
| 616 | 
  | 
      simError(); | 
| 617 | 
  | 
    } | 
| 612 | 
– | 
   | 
| 613 | 
– | 
     | 
| 614 | 
– | 
    setupCutoff(); | 
| 618 | 
  | 
 | 
| 619 | 
  | 
    calcNdf(); | 
| 620 | 
  | 
    calcNdfRaw(); | 
| 664 | 
  | 
    int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | 
| 665 | 
  | 
    int useRF; | 
| 666 | 
  | 
    int useSF; | 
| 667 | 
+ | 
    int useSP; | 
| 668 | 
+ | 
    int useBoxDipole; | 
| 669 | 
  | 
    std::string myMethod; | 
| 670 | 
  | 
 | 
| 671 | 
  | 
    // set the useRF logical | 
| 672 | 
  | 
    useRF = 0; | 
| 673 | 
  | 
    useSF = 0; | 
| 674 | 
+ | 
    useSP = 0; | 
| 675 | 
  | 
 | 
| 676 | 
  | 
 | 
| 677 | 
  | 
    if (simParams_->haveElectrostaticSummationMethod()) { | 
| 678 | 
  | 
      std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 679 | 
  | 
      toUpper(myMethod); | 
| 680 | 
< | 
      if (myMethod == "REACTION_FIELD") { | 
| 681 | 
< | 
        useRF=1; | 
| 682 | 
< | 
      } else { | 
| 683 | 
< | 
        if (myMethod == "SHIFTED_FORCE") { | 
| 684 | 
< | 
          useSF = 1; | 
| 685 | 
< | 
        } | 
| 680 | 
> | 
      if (myMethod == "REACTION_FIELD"){ | 
| 681 | 
> | 
        useRF = 1; | 
| 682 | 
> | 
      } else if (myMethod == "SHIFTED_FORCE"){ | 
| 683 | 
> | 
        useSF = 1; | 
| 684 | 
> | 
      } else if (myMethod == "SHIFTED_POTENTIAL"){ | 
| 685 | 
> | 
        useSP = 1; | 
| 686 | 
  | 
      } | 
| 687 | 
  | 
    } | 
| 688 | 
+ | 
     | 
| 689 | 
+ | 
    if (simParams_->haveAccumulateBoxDipole())  | 
| 690 | 
+ | 
      if (simParams_->getAccumulateBoxDipole()) | 
| 691 | 
+ | 
        useBoxDipole = 1; | 
| 692 | 
  | 
 | 
| 693 | 
  | 
    //loop over all of the atom types | 
| 694 | 
  | 
    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 759 | 
  | 
    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);     | 
| 760 | 
  | 
 | 
| 761 | 
  | 
    temp = useSF; | 
| 762 | 
< | 
    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);     | 
| 762 | 
> | 
    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    | 
| 763 | 
  | 
 | 
| 764 | 
+ | 
    temp = useSP; | 
| 765 | 
+ | 
    MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 766 | 
+ | 
 | 
| 767 | 
+ | 
    temp = useBoxDipole; | 
| 768 | 
+ | 
    MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  | 
| 769 | 
+ | 
 | 
| 770 | 
  | 
#endif | 
| 771 | 
  | 
 | 
| 772 | 
  | 
    fInfo_.SIM_uses_PBC = usePBC;     | 
| 784 | 
  | 
    fInfo_.SIM_uses_FLARB = useFLARB; | 
| 785 | 
  | 
    fInfo_.SIM_uses_RF = useRF; | 
| 786 | 
  | 
    fInfo_.SIM_uses_SF = useSF; | 
| 787 | 
< | 
 | 
| 788 | 
< | 
    if( myMethod == "REACTION_FIELD") { | 
| 773 | 
< | 
       | 
| 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 | 
< | 
 | 
| 787 | 
> | 
    fInfo_.SIM_uses_SP = useSP; | 
| 788 | 
> | 
    fInfo_.SIM_uses_BoxDipole = useBoxDipole; | 
| 789 | 
  | 
  } | 
| 790 | 
  | 
 | 
| 791 | 
  | 
  void SimInfo::setupFortranSim() { | 
| 870 | 
  | 
             "succesfully sent the simulation information to fortran.\n"); | 
| 871 | 
  | 
    MPIcheckPoint(); | 
| 872 | 
  | 
#endif // is_mpi | 
| 873 | 
+ | 
 | 
| 874 | 
+ | 
    // Setup number of neighbors in neighbor list if present | 
| 875 | 
+ | 
    if (simParams_->haveNeighborListNeighbors()) { | 
| 876 | 
+ | 
      setNeighbors(simParams_->getNeighborListNeighbors()); | 
| 877 | 
+ | 
    } | 
| 878 | 
+ | 
    | 
| 879 | 
+ | 
 | 
| 880 | 
  | 
  } | 
| 881 | 
  | 
 | 
| 882 | 
  | 
 | 
| 1068 | 
  | 
    int sm = UNDAMPED; | 
| 1069 | 
  | 
    RealType alphaVal; | 
| 1070 | 
  | 
    RealType dielectric; | 
| 1071 | 
< | 
 | 
| 1071 | 
> | 
     | 
| 1072 | 
  | 
    errorOut = isError; | 
| 1063 | 
– | 
    alphaVal = simParams_->getDampingAlpha(); | 
| 1064 | 
– | 
    dielectric = simParams_->getDielectric(); | 
| 1073 | 
  | 
 | 
| 1074 | 
  | 
    if (simParams_->haveElectrostaticSummationMethod()) { | 
| 1075 | 
  | 
      std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 1086 | 
  | 
            if (myMethod == "SHIFTED_FORCE") {             | 
| 1087 | 
  | 
              esm = SHIFTED_FORCE; | 
| 1088 | 
  | 
            } else { | 
| 1089 | 
< | 
              if (myMethod == "REACTION_FIELD") {              | 
| 1089 | 
> | 
              if (myMethod == "REACTION_FIELD") { | 
| 1090 | 
  | 
                esm = REACTION_FIELD; | 
| 1091 | 
+ | 
                dielectric = simParams_->getDielectric(); | 
| 1092 | 
+ | 
                if (!simParams_->haveDielectric()) { | 
| 1093 | 
+ | 
                  // throw warning | 
| 1094 | 
+ | 
                  sprintf( painCave.errMsg, | 
| 1095 | 
+ | 
                           "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" | 
| 1096 | 
+ | 
                           "\tA default value of %f will be used for the dielectric.\n", dielectric); | 
| 1097 | 
+ | 
                  painCave.isFatal = 0; | 
| 1098 | 
+ | 
                  simError(); | 
| 1099 | 
+ | 
                } | 
| 1100 | 
  | 
              } else { | 
| 1101 | 
  | 
                // throw error         | 
| 1102 | 
  | 
                sprintf( painCave.errMsg, | 
| 1123 | 
  | 
        if (myScreen == "DAMPED") { | 
| 1124 | 
  | 
          sm = DAMPED; | 
| 1125 | 
  | 
          if (!simParams_->haveDampingAlpha()) { | 
| 1126 | 
< | 
            //throw error | 
| 1126 | 
> | 
            // first set a cutoff dependent alpha value | 
| 1127 | 
> | 
            // we assume alpha depends linearly with rcut from 0 to 20.5 ang | 
| 1128 | 
> | 
            alphaVal = 0.5125 - rcut_* 0.025; | 
| 1129 | 
> | 
            // for values rcut > 20.5, alpha is zero | 
| 1130 | 
> | 
            if (alphaVal < 0) alphaVal = 0; | 
| 1131 | 
> | 
 | 
| 1132 | 
> | 
            // throw warning | 
| 1133 | 
  | 
            sprintf( painCave.errMsg, | 
| 1134 | 
  | 
                     "SimInfo warning: dampingAlpha was not specified in the input file.\n" | 
| 1135 | 
< | 
                     "\tA default value of %f (1/ang) will be used.\n", alphaVal); | 
| 1135 | 
> | 
                     "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); | 
| 1136 | 
  | 
            painCave.isFatal = 0; | 
| 1137 | 
  | 
            simError(); | 
| 1138 | 
+ | 
          } else { | 
| 1139 | 
+ | 
            alphaVal = simParams_->getDampingAlpha(); | 
| 1140 | 
  | 
          } | 
| 1141 | 
+ | 
           | 
| 1142 | 
  | 
        } else { | 
| 1143 | 
  | 
          // throw error         | 
| 1144 | 
  | 
          sprintf( painCave.errMsg, | 
| 1187 | 
  | 
 | 
| 1188 | 
  | 
  } | 
| 1189 | 
  | 
 | 
| 1190 | 
+ | 
  void SimInfo::setupAccumulateBoxDipole() {     | 
| 1191 | 
+ | 
 | 
| 1192 | 
+ | 
    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true | 
| 1193 | 
+ | 
    if ( simParams_->haveAccumulateBoxDipole() )  | 
| 1194 | 
+ | 
      if ( simParams_->getAccumulateBoxDipole() ) { | 
| 1195 | 
+ | 
        setAccumulateBoxDipole(); | 
| 1196 | 
+ | 
        calcBoxDipole_ = true; | 
| 1197 | 
+ | 
      } | 
| 1198 | 
+ | 
 | 
| 1199 | 
+ | 
  } | 
| 1200 | 
+ | 
 | 
| 1201 | 
  | 
  void SimInfo::addProperty(GenericData* genData) { | 
| 1202 | 
  | 
    properties_.addProperty(genData);   | 
| 1203 | 
  | 
  } | 
| 1447 | 
  | 
      return angularMomentum; | 
| 1448 | 
  | 
   } | 
| 1449 | 
  | 
    | 
| 1450 | 
< | 
    | 
| 1450 | 
> | 
  StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { | 
| 1451 | 
> | 
    return IOIndexToIntegrableObject.at(index); | 
| 1452 | 
> | 
  } | 
| 1453 | 
> | 
   | 
| 1454 | 
> | 
  void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { | 
| 1455 | 
> | 
    IOIndexToIntegrableObject= v; | 
| 1456 | 
> | 
  } | 
| 1457 | 
> | 
 | 
| 1458 | 
> | 
/* | 
| 1459 | 
> | 
   void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { | 
| 1460 | 
> | 
      assert( v.size() == nAtoms_ + nRigidBodies_); | 
| 1461 | 
> | 
      sdByGlobalIndex_ = v; | 
| 1462 | 
> | 
    } | 
| 1463 | 
> | 
 | 
| 1464 | 
> | 
    StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { | 
| 1465 | 
> | 
      //assert(index < nAtoms_ + nRigidBodies_); | 
| 1466 | 
> | 
      return sdByGlobalIndex_.at(index); | 
| 1467 | 
> | 
    }    | 
| 1468 | 
> | 
*/    | 
| 1469 | 
  | 
}//end namespace oopse | 
| 1470 | 
  | 
 |