| 1 | 
< | 
 /* | 
| 1 | 
> | 
/* | 
| 2 | 
  | 
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. | 
| 3 | 
  | 
 * | 
| 4 | 
  | 
 * The University of Notre Dame grants you ("Licensee") a | 
| 39 | 
  | 
 * such damages. | 
| 40 | 
  | 
 */ | 
| 41 | 
  | 
  | 
| 42 | 
< | 
 /** | 
| 43 | 
< | 
  * @file ForceField.cpp | 
| 44 | 
< | 
  * @author tlin | 
| 45 | 
< | 
  * @date 11/04/2004 | 
| 46 | 
< | 
  * @time 22:51am | 
| 47 | 
< | 
  * @version 1.0 | 
| 48 | 
< | 
  */ | 
| 42 | 
> | 
/** | 
| 43 | 
> | 
 * @file ForceField.cpp | 
| 44 | 
> | 
 * @author tlin | 
| 45 | 
> | 
 * @date 11/04/2004 | 
| 46 | 
> | 
 * @time 22:51am | 
| 47 | 
> | 
 * @version 1.0 | 
| 48 | 
> | 
 */ | 
| 49 | 
  | 
   | 
| 50 | 
  | 
#include "UseTheForce/ForceField.hpp" | 
| 51 | 
  | 
#include "utils/simError.h" | 
| 52 | 
+ | 
#include "UseTheForce/DarkSide/atype_interface.h" | 
| 53 | 
+ | 
#include "UseTheForce/DarkSide/fForceOptions_interface.h" | 
| 54 | 
+ | 
#include "UseTheForce/DarkSide/switcheroo_interface.h" | 
| 55 | 
  | 
namespace oopse { | 
| 56 | 
  | 
 | 
| 57 | 
< | 
ForceField::ForceField() {  | 
| 57 | 
> | 
  ForceField::ForceField() {  | 
| 58 | 
  | 
    char* tempPath;  | 
| 59 | 
  | 
    tempPath = getenv("FORCE_PARAM_PATH"); | 
| 60 | 
  | 
 | 
| 61 | 
  | 
    if (tempPath == NULL) { | 
| 62 | 
< | 
        //convert a macro from compiler to a string in c++ | 
| 63 | 
< | 
        STR_DEFINE(ffPath_, FRC_PATH ); | 
| 62 | 
> | 
      //convert a macro from compiler to a string in c++ | 
| 63 | 
> | 
      STR_DEFINE(ffPath_, FRC_PATH ); | 
| 64 | 
  | 
    } else { | 
| 65 | 
< | 
        ffPath_ = tempPath; | 
| 65 | 
> | 
      ffPath_ = tempPath; | 
| 66 | 
  | 
    } | 
| 67 | 
< | 
} | 
| 67 | 
> | 
  } | 
| 68 | 
  | 
 | 
| 69 | 
< | 
AtomType* ForceField::getAtomType(const std::string &at) { | 
| 69 | 
> | 
 | 
| 70 | 
> | 
  ForceField::~ForceField() { | 
| 71 | 
> | 
    deleteAtypes(); | 
| 72 | 
> | 
    deleteSwitch(); | 
| 73 | 
> | 
  } | 
| 74 | 
> | 
 | 
| 75 | 
> | 
  AtomType* ForceField::getAtomType(const std::string &at) { | 
| 76 | 
  | 
    std::vector<std::string> keys; | 
| 77 | 
  | 
    keys.push_back(at); | 
| 78 | 
  | 
    return atomTypeCont_.find(keys); | 
| 79 | 
< | 
} | 
| 79 | 
> | 
  } | 
| 80 | 
  | 
 | 
| 81 | 
< | 
BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) { | 
| 81 | 
> | 
  BondType* ForceField::getBondType(const std::string &at1,  | 
| 82 | 
> | 
                                    const std::string &at2) { | 
| 83 | 
  | 
    std::vector<std::string> keys; | 
| 84 | 
  | 
    keys.push_back(at1); | 
| 85 | 
  | 
    keys.push_back(at2);     | 
| 87 | 
  | 
    //try exact match first | 
| 88 | 
  | 
    BondType* bondType = bondTypeCont_.find(keys); | 
| 89 | 
  | 
    if (bondType) { | 
| 90 | 
< | 
        return bondType; | 
| 90 | 
> | 
      return bondType; | 
| 91 | 
  | 
    } else { | 
| 92 | 
< | 
        //if no exact match found, try wild card match | 
| 93 | 
< | 
        return bondTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 92 | 
> | 
      //if no exact match found, try wild card match | 
| 93 | 
> | 
      return bondTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 94 | 
  | 
    } | 
| 95 | 
+ | 
  } | 
| 96 | 
  | 
 | 
| 97 | 
< | 
} | 
| 98 | 
< | 
 | 
| 99 | 
< | 
BendType* ForceField::getBendType(const std::string &at1, const std::string &at2, | 
| 89 | 
< | 
                        const std::string &at3) { | 
| 97 | 
> | 
  BendType* ForceField::getBendType(const std::string &at1,  | 
| 98 | 
> | 
                                    const std::string &at2, | 
| 99 | 
> | 
                                    const std::string &at3) { | 
| 100 | 
  | 
    std::vector<std::string> keys; | 
| 101 | 
  | 
    keys.push_back(at1); | 
| 102 | 
  | 
    keys.push_back(at2);     | 
| 105 | 
  | 
    //try exact match first | 
| 106 | 
  | 
    BendType* bendType = bendTypeCont_.find(keys); | 
| 107 | 
  | 
    if (bendType) { | 
| 108 | 
< | 
        return bendType; | 
| 108 | 
> | 
      return bendType; | 
| 109 | 
  | 
    } else { | 
| 110 | 
< | 
        //if no exact match found, try wild card match | 
| 111 | 
< | 
        return bendTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 110 | 
> | 
      //if no exact match found, try wild card match | 
| 111 | 
> | 
      return bendTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 112 | 
  | 
    } | 
| 113 | 
< | 
} | 
| 113 | 
> | 
  } | 
| 114 | 
  | 
 | 
| 115 | 
< | 
TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2, | 
| 116 | 
< | 
                              const std::string &at3, const std::string &at4) { | 
| 115 | 
> | 
  TorsionType* ForceField::getTorsionType(const std::string &at1,  | 
| 116 | 
> | 
                                          const std::string &at2, | 
| 117 | 
> | 
                                          const std::string &at3,  | 
| 118 | 
> | 
                                          const std::string &at4) { | 
| 119 | 
  | 
    std::vector<std::string> keys; | 
| 120 | 
  | 
    keys.push_back(at1); | 
| 121 | 
  | 
    keys.push_back(at2);     | 
| 124 | 
  | 
 | 
| 125 | 
  | 
    TorsionType* torsionType = torsionTypeCont_.find(keys); | 
| 126 | 
  | 
    if (torsionType) { | 
| 127 | 
< | 
        return torsionType; | 
| 127 | 
> | 
      return torsionType; | 
| 128 | 
  | 
    } else { | 
| 129 | 
< | 
        //if no exact match found, try wild card match | 
| 130 | 
< | 
        return torsionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 129 | 
> | 
      //if no exact match found, try wild card match | 
| 130 | 
> | 
      return torsionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 131 | 
  | 
    } | 
| 132 | 
  | 
     | 
| 133 | 
  | 
    return torsionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 134 | 
+ | 
  } | 
| 135 | 
  | 
 | 
| 136 | 
< | 
} | 
| 124 | 
< | 
 | 
| 125 | 
< | 
BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){  | 
| 136 | 
> | 
  NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) { | 
| 137 | 
  | 
    std::vector<std::string> keys; | 
| 138 | 
  | 
    keys.push_back(at1); | 
| 139 | 
  | 
    keys.push_back(at2);     | 
| 140 | 
+ | 
     | 
| 141 | 
+ | 
    //try exact match first | 
| 142 | 
+ | 
    NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys); | 
| 143 | 
+ | 
    if (nbiType) { | 
| 144 | 
+ | 
      return nbiType; | 
| 145 | 
+ | 
    } else { | 
| 146 | 
+ | 
      //if no exact match found, try wild card match | 
| 147 | 
+ | 
      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 148 | 
+ | 
    }     | 
| 149 | 
+ | 
  } | 
| 150 | 
+ | 
   | 
| 151 | 
+ | 
  BondType* ForceField::getExactBondType(const std::string &at1,  | 
| 152 | 
+ | 
                                         const std::string &at2){  | 
| 153 | 
+ | 
    std::vector<std::string> keys; | 
| 154 | 
+ | 
    keys.push_back(at1); | 
| 155 | 
+ | 
    keys.push_back(at2);     | 
| 156 | 
  | 
    return bondTypeCont_.find(keys); | 
| 157 | 
< | 
} | 
| 158 | 
< | 
 | 
| 159 | 
< | 
BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2, | 
| 160 | 
< | 
                            const std::string &at3){  | 
| 157 | 
> | 
  } | 
| 158 | 
> | 
   | 
| 159 | 
> | 
  BendType* ForceField::getExactBendType(const std::string &at1,  | 
| 160 | 
> | 
                                         const std::string &at2, | 
| 161 | 
> | 
                                         const std::string &at3){  | 
| 162 | 
  | 
    std::vector<std::string> keys; | 
| 163 | 
  | 
    keys.push_back(at1); | 
| 164 | 
  | 
    keys.push_back(at2);     | 
| 165 | 
  | 
    keys.push_back(at3);     | 
| 166 | 
  | 
    return bendTypeCont_.find(keys); | 
| 167 | 
< | 
} | 
| 168 | 
< | 
 | 
| 169 | 
< | 
TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2, | 
| 170 | 
< | 
                                  const std::string &at3, const std::string &at4){  | 
| 167 | 
> | 
  } | 
| 168 | 
> | 
   | 
| 169 | 
> | 
  TorsionType* ForceField::getExactTorsionType(const std::string &at1,  | 
| 170 | 
> | 
                                               const std::string &at2, | 
| 171 | 
> | 
                                               const std::string &at3,  | 
| 172 | 
> | 
                                               const std::string &at4){  | 
| 173 | 
  | 
    std::vector<std::string> keys; | 
| 174 | 
  | 
    keys.push_back(at1); | 
| 175 | 
  | 
    keys.push_back(at2);     | 
| 176 | 
  | 
    keys.push_back(at3);     | 
| 177 | 
  | 
    keys.push_back(at4);    | 
| 178 | 
  | 
    return torsionTypeCont_.find(keys); | 
| 179 | 
< | 
} | 
| 180 | 
< | 
bool ForceField::addAtomType(const std::string &at, AtomType* atomType) { | 
| 179 | 
> | 
  } | 
| 180 | 
> | 
 | 
| 181 | 
> | 
  NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){  | 
| 182 | 
  | 
    std::vector<std::string> keys; | 
| 183 | 
+ | 
    keys.push_back(at1); | 
| 184 | 
+ | 
    keys.push_back(at2);     | 
| 185 | 
+ | 
    return nonBondedInteractionTypeCont_.find(keys); | 
| 186 | 
+ | 
  } | 
| 187 | 
+ | 
 | 
| 188 | 
+ | 
 | 
| 189 | 
+ | 
  bool ForceField::addAtomType(const std::string &at, AtomType* atomType) { | 
| 190 | 
+ | 
    std::vector<std::string> keys; | 
| 191 | 
  | 
    keys.push_back(at); | 
| 192 | 
  | 
    return atomTypeCont_.add(keys, atomType); | 
| 193 | 
< | 
} | 
| 193 | 
> | 
  } | 
| 194 | 
  | 
 | 
| 195 | 
< | 
bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) { | 
| 195 | 
> | 
  bool ForceField::addBondType(const std::string &at1, const std::string &at2, | 
| 196 | 
> | 
                               BondType* bondType) { | 
| 197 | 
  | 
    std::vector<std::string> keys; | 
| 198 | 
  | 
    keys.push_back(at1); | 
| 199 | 
  | 
    keys.push_back(at2);     | 
| 200 | 
< | 
    return bondTypeCont_.add(keys, bondType); | 
| 201 | 
< | 
 | 
| 202 | 
< | 
} | 
| 203 | 
< | 
 | 
| 204 | 
< | 
bool ForceField::addBendType(const std::string &at1, const std::string &at2, | 
| 165 | 
< | 
                        const std::string &at3, BendType* bendType) { | 
| 200 | 
> | 
    return bondTypeCont_.add(keys, bondType);     | 
| 201 | 
> | 
  } | 
| 202 | 
> | 
   | 
| 203 | 
> | 
  bool ForceField::addBendType(const std::string &at1, const std::string &at2, | 
| 204 | 
> | 
                               const std::string &at3, BendType* bendType) { | 
| 205 | 
  | 
    std::vector<std::string> keys; | 
| 206 | 
  | 
    keys.push_back(at1); | 
| 207 | 
  | 
    keys.push_back(at2);     | 
| 208 | 
  | 
    keys.push_back(at3);     | 
| 209 | 
  | 
    return bendTypeCont_.add(keys, bendType); | 
| 210 | 
< | 
} | 
| 211 | 
< | 
 | 
| 212 | 
< | 
bool ForceField::addTorsionType(const std::string &at1, const std::string &at2, | 
| 213 | 
< | 
                              const std::string &at3, const std::string &at4, TorsionType* torsionType) { | 
| 210 | 
> | 
  } | 
| 211 | 
> | 
   | 
| 212 | 
> | 
  bool ForceField::addTorsionType(const std::string &at1,  | 
| 213 | 
> | 
                                  const std::string &at2, | 
| 214 | 
> | 
                                  const std::string &at3,  | 
| 215 | 
> | 
                                  const std::string &at4,  | 
| 216 | 
> | 
                                  TorsionType* torsionType) { | 
| 217 | 
  | 
    std::vector<std::string> keys; | 
| 218 | 
  | 
    keys.push_back(at1); | 
| 219 | 
  | 
    keys.push_back(at2);     | 
| 220 | 
  | 
    keys.push_back(at3);     | 
| 221 | 
  | 
    keys.push_back(at4);     | 
| 222 | 
  | 
    return torsionTypeCont_.add(keys, torsionType); | 
| 223 | 
< | 
} | 
| 223 | 
> | 
  } | 
| 224 | 
  | 
 | 
| 225 | 
< | 
double ForceField::getRcutFromAtomType(AtomType* at) { | 
| 225 | 
> | 
  bool ForceField::addNonBondedInteractionType(const std::string &at1,  | 
| 226 | 
> | 
                                               const std::string &at2,  | 
| 227 | 
> | 
                                               NonBondedInteractionType* nbiType) { | 
| 228 | 
> | 
    std::vector<std::string> keys; | 
| 229 | 
> | 
    keys.push_back(at1); | 
| 230 | 
> | 
    keys.push_back(at2);     | 
| 231 | 
> | 
    return nonBondedInteractionTypeCont_.add(keys, nbiType); | 
| 232 | 
> | 
  } | 
| 233 | 
> | 
   | 
| 234 | 
> | 
  RealType ForceField::getRcutFromAtomType(AtomType* at) { | 
| 235 | 
  | 
    /**@todo */ | 
| 236 | 
  | 
    GenericData* data; | 
| 237 | 
< | 
    double rcut = 0.0; | 
| 238 | 
< | 
 | 
| 237 | 
> | 
    RealType rcut = 0.0; | 
| 238 | 
> | 
     | 
| 239 | 
  | 
    if (at->isLennardJones()) { | 
| 240 | 
< | 
        data = at->getPropertyByName("LennardJones"); | 
| 241 | 
< | 
        if (data != NULL) { | 
| 242 | 
< | 
            LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data); | 
| 243 | 
< | 
 | 
| 244 | 
< | 
            if (ljData != NULL) { | 
| 245 | 
< | 
                LJParam ljParam = ljData->getData(); | 
| 246 | 
< | 
 | 
| 247 | 
< | 
                //by default use 2.5*sigma as cutoff radius | 
| 248 | 
< | 
                rcut = 2.5 * ljParam.sigma; | 
| 249 | 
< | 
                 | 
| 250 | 
< | 
            } else { | 
| 251 | 
< | 
                    sprintf( painCave.errMsg, | 
| 252 | 
< | 
                           "Can not cast GenericData to LJParam\n"); | 
| 253 | 
< | 
                    painCave.severity = OOPSE_ERROR; | 
| 254 | 
< | 
                    painCave.isFatal = 1; | 
| 255 | 
< | 
                    simError();           | 
| 256 | 
< | 
            }             | 
| 257 | 
< | 
        } else { | 
| 258 | 
< | 
            sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n"); | 
| 259 | 
< | 
            painCave.severity = OOPSE_ERROR; | 
| 260 | 
< | 
            painCave.isFatal = 1; | 
| 261 | 
< | 
            simError();           | 
| 262 | 
< | 
        } | 
| 240 | 
> | 
      data = at->getPropertyByName("LennardJones"); | 
| 241 | 
> | 
      if (data != NULL) { | 
| 242 | 
> | 
        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data); | 
| 243 | 
> | 
         | 
| 244 | 
> | 
        if (ljData != NULL) { | 
| 245 | 
> | 
          LJParam ljParam = ljData->getData(); | 
| 246 | 
> | 
           | 
| 247 | 
> | 
          //by default use 2.5*sigma as cutoff radius | 
| 248 | 
> | 
          rcut = 2.5 * ljParam.sigma; | 
| 249 | 
> | 
           | 
| 250 | 
> | 
        } else { | 
| 251 | 
> | 
          sprintf( painCave.errMsg, | 
| 252 | 
> | 
                   "Can not cast GenericData to LJParam\n"); | 
| 253 | 
> | 
          painCave.severity = OOPSE_ERROR; | 
| 254 | 
> | 
          painCave.isFatal = 1; | 
| 255 | 
> | 
          simError();           | 
| 256 | 
> | 
        }             | 
| 257 | 
> | 
      } else { | 
| 258 | 
> | 
        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n"); | 
| 259 | 
> | 
        painCave.severity = OOPSE_ERROR; | 
| 260 | 
> | 
        painCave.isFatal = 1; | 
| 261 | 
> | 
        simError();           | 
| 262 | 
> | 
      } | 
| 263 | 
  | 
    } | 
| 213 | 
– | 
 | 
| 264 | 
  | 
    return rcut;     | 
| 265 | 
< | 
} | 
| 265 | 
> | 
  } | 
| 266 | 
> | 
   | 
| 267 | 
  | 
 | 
| 268 | 
< | 
 | 
| 218 | 
< | 
ifstrstream* ForceField::openForceFieldFile(const std::string& filename) { | 
| 268 | 
> | 
  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) { | 
| 269 | 
  | 
    std::string forceFieldFilename(filename); | 
| 270 | 
  | 
    ifstrstream* ffStream = new ifstrstream(); | 
| 271 | 
  | 
     | 
| 273 | 
  | 
    ffStream->open(forceFieldFilename.c_str()); | 
| 274 | 
  | 
    if(!ffStream->is_open()){ | 
| 275 | 
  | 
 | 
| 276 | 
< | 
        forceFieldFilename = ffPath_ + "/" + forceFieldFilename; | 
| 277 | 
< | 
        ffStream->open( forceFieldFilename.c_str() ); | 
| 276 | 
> | 
      forceFieldFilename = ffPath_ + "/" + forceFieldFilename; | 
| 277 | 
> | 
      ffStream->open( forceFieldFilename.c_str() ); | 
| 278 | 
  | 
 | 
| 279 | 
< | 
        //if current directory does not contain the force field file, | 
| 280 | 
< | 
        //try to open it in the path         | 
| 281 | 
< | 
        if(!ffStream->is_open()){ | 
| 279 | 
> | 
      //if current directory does not contain the force field file, | 
| 280 | 
> | 
      //try to open it in the path         | 
| 281 | 
> | 
      if(!ffStream->is_open()){ | 
| 282 | 
  | 
 | 
| 283 | 
< | 
            sprintf( painCave.errMsg, | 
| 284 | 
< | 
               "Error opening the force field parameter file:\n" | 
| 285 | 
< | 
               "\t%s\n" | 
| 286 | 
< | 
               "\tHave you tried setting the FORCE_PARAM_PATH environment " | 
| 287 | 
< | 
               "variable?\n", | 
| 288 | 
< | 
               forceFieldFilename.c_str() ); | 
| 289 | 
< | 
            painCave.severity = OOPSE_ERROR; | 
| 290 | 
< | 
            painCave.isFatal = 1; | 
| 291 | 
< | 
            simError(); | 
| 292 | 
< | 
        } | 
| 283 | 
> | 
        sprintf( painCave.errMsg, | 
| 284 | 
> | 
                 "Error opening the force field parameter file:\n" | 
| 285 | 
> | 
                 "\t%s\n" | 
| 286 | 
> | 
                 "\tHave you tried setting the FORCE_PARAM_PATH environment " | 
| 287 | 
> | 
                 "variable?\n", | 
| 288 | 
> | 
                 forceFieldFilename.c_str() ); | 
| 289 | 
> | 
        painCave.severity = OOPSE_ERROR; | 
| 290 | 
> | 
        painCave.isFatal = 1; | 
| 291 | 
> | 
        simError(); | 
| 292 | 
> | 
      } | 
| 293 | 
  | 
    }   | 
| 244 | 
– | 
 | 
| 294 | 
  | 
    return ffStream; | 
| 295 | 
+ | 
  } | 
| 296 | 
  | 
 | 
| 297 | 
< | 
} | 
| 298 | 
< | 
 | 
| 297 | 
> | 
  void ForceField::setFortranForceOptions(){ | 
| 298 | 
> | 
    ForceOptions theseFortranOptions; | 
| 299 | 
> | 
    forceFieldOptions_.makeFortranOptions(theseFortranOptions); | 
| 300 | 
> | 
    setfForceOptions(&theseFortranOptions); | 
| 301 | 
> | 
  } | 
| 302 | 
  | 
} //end namespace oopse |