src/kernel/openmm_wrapper.cpp

   1 #include <cmath>
   2 #include <set>
   3 #include <iostream>
   4 #include <sstream>
   5 #include <fstream>
   6 #include <map>
   7 #include <vector>
   8 #include <cctype>
   9 #include <algorithm>
  10
  11 using namespace std;
  12
  13 #include "OpenMM.h"
  14
  15 #include "gmx_fatal.h"
  16 #include "typedefs.h"
  17 #include "mdrun.h"
  18 #include "physics.h"
  19 #include "string2.h"
  20 #include "gmx_gpu_utils.h"
  21 #include "mtop_util.h"
  22
  23 #include "openmm_wrapper.h"
  24
  25 using namespace OpenMM;
  26
  27 #define MEM_ERR_MSG(str) \
  28     "The %s-simulation GPU memory test detected errors. As memory errors would cause incorrect " \
  29     "simulation results, gromacs has aborted execution.\n Make sure that your GPU's memory is not " \
  30     "overclocked and that the device is properly cooled.\n", (str)
  31
  32 /** Convert string to type T. */
  33 template <class T>
  34 static bool from_string(T& t, const string& s, ios_base& (*f)(ios_base&))
  35 {
  36   istringstream iss(s);
  37   return !(iss >> f >> t).fail();
  38 }
  39
  40 /**
  41  * Split string around a given delimiter.
  42  */
  43 static vector<string> split(const string &s, char delim)
  44 {
  45     vector<string> elems;
  46     stringstream ss(s);
  47     string item;
  48     while(getline(ss, item, delim))
  49     {
  50         if (item.length() != 0)
  51             elems.push_back(item);
  52     }
  53     return elems;
  54 }
  55
  56 /**
  57  * Split a string "option=value" into "option" and "value" strings.
  58  */
  59 static void split_opt_val(const string &s, string &opt, string &val)
  60 {
  61     size_t eqPos = s.find('=');
  62     if (eqPos != string::npos)
  63     {
  64         opt = s.substr(0, eqPos);
  65         if (eqPos != s.length())  val = s.substr(eqPos+1);
  66     }
  67 }
  68
  69 /**
  70  * Compare two strings ignoring case.
  71  */
  72 static bool isStringEqNCase(const string s1, const string s2)
  73 {
  74     return (gmx_strncasecmp(s1.c_str(), s2.c_str(), max(s1.length(), s2.length())) == 0);
  75 }
  76
  77 /**
  78  * Conver string to upper case.
  79  */
  80 static string toUpper(const string &s)
  81 {
  82     string stmp(s);
  83     std::transform(stmp.begin(), stmp.end(), stmp.begin(), static_cast < int(*)(int) > (toupper));
  84     return stmp;
  85 }
  86
  87 #define SIZEOF_PLATFORMS    1
  88 #define SIZEOF_MEMTESTS     3
  89 #define SIZEOF_DEVICEIDS    1
  90 #define SIZEOF_FORCE_DEV    2
  91
  92 /** Possible options of the mdrun -device parameter. */
  93 static const char *devOptStrings[] = { "platform", "deviceid", "memtest", "force-device" };
  94 enum devOpt {
  95         PLATFORM     = 0,
  96         DEVICEID     = 1,
  97         MEMTEST      = 2,
  98         FORCE_DEVICE = 3
  99 };
 100
 101 /**
 102  * Parse the options string provided to mdrun.
 103  */
 104 class GmxOpenMMPlatformOptions
 105 {
 106 public:
 107   GmxOpenMMPlatformOptions(const char *optionString);
 108 //    GmxOpenMMPlatformOptions(string &optionString);
 109     ~GmxOpenMMPlatformOptions() { options.clear(); }
 110     string getOptionValue(const string &optName);
 111     void remOption(const string &optName) { options.erase(toUpper(optName)); }
 112 private:
 113     void setOption(const string &opt, const string &val);
 114     map<string, string> options;
 115     static const char * const platforms[SIZEOF_PLATFORMS];
 116     static const char * const memtests[SIZEOF_MEMTESTS];
 117     static const char * const deviceid[SIZEOF_DEVICEIDS];
 118     static const char * const force_dev[SIZEOF_FORCE_DEV];
 119 };
 120
 121 /* possible values of the parameters that will not be passed to OMM
 122    first value will always be the default used in case if the option is not specified */
 123 const char * const GmxOpenMMPlatformOptions::platforms[SIZEOF_PLATFORMS] = { "Cuda" /*,"OpenCL"*/ };
 124 const char * const GmxOpenMMPlatformOptions::memtests[SIZEOF_MEMTESTS]  = { "15", "full", "off" }; /* also valid any positive integer >15 */
 125 const char * const GmxOpenMMPlatformOptions::deviceid[SIZEOF_DEVICEIDS]  = { "0" };       /* also valid any positive integer */
 126 const char * const GmxOpenMMPlatformOptions::force_dev[SIZEOF_FORCE_DEV]   = { "no", "yes" };
 127
 128 GmxOpenMMPlatformOptions::GmxOpenMMPlatformOptions(const char *optionString)
 129 {
 130     // set default values
 131     setOption("platform",       platforms[0]);
 132     setOption("memtest",        memtests[0]);
 133     setOption("deviceid",       deviceid[0]);
 134     setOption("force-device",   force_dev[0]);
 135
 136     string opt(optionString);
 137
 138     // remove all the whitespaces
 139     opt.erase(remove_if(opt.begin(), opt.end(), ::isspace), opt.end());
 140     // tokenize around ","-s
 141     vector<string> tokens = split(opt, ',');
 142
 143     for (vector<string>::iterator it = tokens.begin(); it != tokens.end(); ++it)
 144     {
 145         string opt = "", val = "";
 146         split_opt_val(*it, opt, val);
 147
 148         if (isStringEqNCase(opt, "platform"))
 149         {
 150             setOption(opt, val);
 151             continue;
 152         }
 153
 154         if (isStringEqNCase(opt, "memtest"))
 155         {
 156             if (!isStringEqNCase(val, "full") && !isStringEqNCase(val, "off")) /* the value has to be an integer >15 */
 157             {
 158                 int secs;
 159                 if (!from_string<int>(secs, val, std::dec))
 160                 {
 161                     gmx_fatal(FARGS, "Invalid value for option memtestoption: \"%s\"!\n", val.c_str());
 162                 }
 163                 if (secs < 15)
 164                 {
 165                     gmx_fatal(FARGS, "Incorrect value for memtest option (%d). Memtest needs to run for at least 15s!\n", secs);
 166                 }
 167             }
 168             setOption(opt, val);
 169             continue;
 170         }
 171
 172         if (isStringEqNCase(opt, "deviceid"))
 173         {
 174             int id;
 175             if (!from_string<int>(id, val, std::dec) )
 176             {
 177                 gmx_fatal(FARGS, "Invalid device id: \"%s\"!\n", val.c_str());
 178             }
 179             setOption(opt, val);
 180             continue;
 181         }
 182
 183         if (isStringEqNCase(opt,"force-device"))
 184         {
 185             if (!isStringEqNCase(val, "yes") && !isStringEqNCase(val, "no"))
 186             {
 187                 gmx_fatal(FARGS, "Invalid OpenMM force option: \"%s\"!\n", val.c_str());
 188             }
 189             setOption(opt, val);
 190             continue;
 191         }
 192
 193         // if we got till here something went wrong
 194         gmx_fatal(FARGS, "Invalid OpenMM platform option: \"%s\"!\n", (*it).c_str());
 195     }
 196
 197 /*    cout << "== platform = " << getOptionValue("platform") << endl
 198          << "== deviceID = " << getOptionValue("deviceid") << endl
 199          << "== memtest  = " << getOptionValue("memtest") << endl
 200          << "== force    = " << getOptionValue("force-device") << endl;
 201 */
 202 }
 203
 204 string GmxOpenMMPlatformOptions::getOptionValue(const string &optName)
 205 {
 206     if (options.find(toUpper(optName)) != options.end())
 207     {
 208         return options[toUpper(optName)];
 209     }
 210     else
 211     {
 212         return NULL;
 213     }
 214 }
 215
 216 void GmxOpenMMPlatformOptions::setOption(const string &opt, const string &val)
 217 {
 218     options[toUpper(opt)] = val;
 219 }
 220
 221 class OpenMMData
 222 {
 223 public:
 224     System* system;
 225     Context* context;
 226     Integrator* integrator;
 227     bool removeCM;
 228     GmxOpenMMPlatformOptions *platformOpt;
 229 };
 230
 231 void run_memtest(FILE* fplog, int devId, const char* pre_post, GmxOpenMMPlatformOptions *opt)
 232 {
 233     // TODO I guess the log entries should go on stdout as well, shouldn't they? [sz]
 234     int which_test;
 235     int res;
 236     const char * test_type = opt->getOptionValue("memtest").c_str();
 237     if (!gmx_strcasecmp(test_type, "off"))
 238     {
 239             which_test = 0;
 240     }
 241     else
 242     {
 243         if (!gmx_strcasecmp(test_type, "full"))
 244         {
 245             which_test = 2;
 246         }
 247         else
 248         {
 249              from_string<int>(which_test, test_type, std::dec);
 250         }
 251     }
 252
 253     switch (which_test)
 254     {
 255         case 0: /* no memtest */
 256
 257             sprintf(warn_buf, "%s-simulation GPU memtest skipped. Note, that faulty memory can cause "
 258                     "incorrect results!\n", pre_post);
 259             fprintf(fplog, "%s", warn_buf);
 260             warning_note(NULL);
 261             break; /* case 0 */
 262
 263         case 1: /* quick memtest */
 264             fprintf(fplog,  "%s-simulation %s GPU memtest in progress...\n", pre_post, test_type);
 265             fprintf(stdout, "\n%s-simulation %s GPU memtest in progress...", pre_post, test_type);
 266             fflush(fplog); fflush(stdout);
 267             if (do_quick_memtest(-1) != 0)
 268             {
 269                 gmx_fatal(FARGS,MEM_ERR_MSG(pre_post));
 270
 271             }
 272             else
 273             {
 274                 fprintf(fplog,  "Memory test completed without errors.\n");
 275                 fprintf(stdout, "done, no errors detected\n");
 276             }
 277             break; /* case 1 */
 278
 279         case 2: /* full memtest */
 280             fprintf(fplog,  "%s-simulation %s memtest in progress...\n", pre_post, test_type);
 281             fprintf(stdout, "\n%s-simulation %s memtest in progress...", pre_post, test_type);
 282             fflush(fplog); fflush(stdout);
 283             if (do_full_memtest(-1) != 0)
 284             {
 285                 gmx_fatal(FARGS, MEM_ERR_MSG(pre_post) );
 286
 287             }
 288             else
 289             {
 290                 fprintf(fplog, "Memory test completed without errors.\n");
 291                 fprintf(stdout, "done, no errors detected\n");
 292             }
 293             break; /* case 2 */
 294
 295         default: /* timed memtest */
 296             fprintf(fplog,  "%s-simulation memtest for ~%ds in progress...\n", pre_post, which_test);
 297             fprintf(stdout, "\n%s-simulation memtest for ~%ds in progress...", pre_post, which_test);
 298             fflush(fplog); fflush(stdout);
 299             if (do_timed_memtest(-1, which_test) != 0)
 300             {
 301                 gmx_fatal(FARGS, MEM_ERR_MSG(pre_post));
 302
 303             }
 304             else
 305             {
 306                 fprintf(fplog, "Memory test completed without errors.\n");
 307                 fprintf(stdout, "done, no errors detected.\n");
 308             }
 309     }
 310     fflush(fplog); fflush(stdout);
 311 }
 312
 313 void checkGmxOptions(t_inputrec *ir, gmx_localtop_t *top)
 314 {
 315     // Abort if unsupported critical options are present
 316
 317     /* Integrator */
 318     if ( (ir->eI !=  eiVV)   &&
 319          (ir->eI !=  eiVVAK) &&
 320          (ir->eI !=  eiSD1)  &&
 321          (ir->eI !=  eiSD2)  &&
 322          (ir->eI !=  eiBD)
 323        )
 324     {
 325         gmx_fatal(FARGS, "** OpenMM Error ** : Unsupported integrator requested."
 326                 "Available integrators are: md-vv/md-vvak, sd/sd1, and bd. \n");
 327     }
 328
 329     /* Electroctstics */
 330     if ((ir->coulombtype != eelPME) && (ir->coulombtype != eelRF) &&  (ir->coulombtype != eelEWALD))
 331     {
 332         gmx_fatal(FARGS,"** OpenMM Error ** : Unsupported method for electrostatics."
 333             "Use NoCutoff (i.e. rcoulomb = rvdw = 0 ),Reaction-Field, Ewald or PME.\n");
 334     }
 335
 336     /* TODO what is this code checking? */
 337     if (ir->bPeriodicMols)
 338         gmx_fatal(FARGS,"** OpenMM Error ** : Systems with periodic molecules ares not supported\n");
 339
 340     if ( (ir->etc != etcNO) && (ir->etc != etcANDERSEN) && (ir->etc != etcANDERSENINTERVAL))
 341         gmx_fatal(FARGS,"** OpenMM Error ** : Temperature coupling can be achieved by using either \n\t(1)\t\"md-vv\" or \"md-vvak\" integrators with \"andersen\" or \"andersen-interval\" thermostat, or \n\t(2)\t\"sd\",\"sd1\" or \"bd\" integrators\n");
 342
 343         if (ir->opts.ngtc > 1)
 344                 gmx_fatal(FARGS,"** OpenMM Error ** : Multiple temperature coupling groups are not supported. \n");
 345
 346     if (ir->epc != etcNO)
 347         gmx_fatal(FARGS,"** OpenMM Error ** : Pressure coupling is not supported. \n");
 348
 349     if (ir->opts.annealing[0])
 350         gmx_fatal(FARGS,"** OpenMM Error ** : Simulated annealing is not supported. \n");
 351
 352     if (ir->eConstrAlg != econtSHAKE)
 353         warning_note("** OpenMM Warning ** : Constraints in OpenMM are done by a combination of SHAKE, SETTLE and CCMA. Accuracy is based on the SHAKE tolerance set by the \"shake_tol\" option.");
 354
 355     if (ir->nwall != 0)
 356         gmx_fatal(FARGS,"** OpenMM Error ** : Walls are not supported.\n");
 357
 358     if (ir->ePull != epullNO)
 359         gmx_fatal(FARGS,"** OpenMM Error ** : Pulling is not supported.\n");
 360
 361     if (top->idef.il[F_DISRES].nr > 0)
 362         gmx_fatal(FARGS,"** OpenMM Error ** : Distant restraints are not supported.\n");
 363
 364     if (top->idef.il[F_ORIRES].nr > 0)
 365         gmx_fatal(FARGS,"** OpenMM Error ** : Orientation restraints are not supported.\n");
 366
 367     if (top->idef.il[F_ANGRES].nr > 0)
 368         gmx_fatal(FARGS,"** OpenMM Error ** : Angle restraints are not supported.\n");
 369
 370     if (top->idef.il[F_DIHRES].nr > 0)
 371         gmx_fatal(FARGS,"** OpenMM Error ** : Dihedral restraints are not supported.\n");
 372
 373     if (ir->efep != efepNO)
 374         gmx_fatal(FARGS,"** OpenMM Error ** : Free energy calculations are not supported.\n");
 375
 376     if (ir->opts.ngacc > 1)
 377         gmx_fatal(FARGS,"** OpenMM Error ** : Non-equilibrium MD (accelerated groups) are not supported. \n");
 378
 379     if (IR_ELEC_FIELD(*ir))
 380         gmx_fatal(FARGS,"** OpenMM Error ** : Electric fields are not supported. \n");
 381
 382     if (ir->bQMMM)
 383         gmx_fatal(FARGS,"** OpenMM Error ** : QMMM calculations are not supported. \n");
 384
 385     if (ir->rcoulomb != ir->rvdw)
 386         gmx_fatal(FARGS,"** OpenMM Error ** : OpenMM uses a single cutoff for both Coulomb and VdW interactions. Please set rcoulomb equal to rvdw. \n");
 387
 388 }
 389
 390 /**
 391  * Initialize OpenMM, and return a pointer to the OpenMMData object containing the System, Integrator, and Context.
 392  */
 393 void* openmm_init(FILE *fplog, const char *platformOptStr,
 394                     t_commrec *cr,t_inputrec *ir,
 395                     gmx_mtop_t *top_global, gmx_localtop_t *top,
 396                     t_mdatoms *mdatoms, t_forcerec *fr,t_state *state)
 397 {
 398     static bool hasLoadedPlugins = false;
 399     string usedPluginDir;
 400
 401     try {
 402     if (!hasLoadedPlugins)
 403     {
 404         vector<string> loadedPlugins;
 405     /*  Look for OpenMM plugins at various locations (listed in order of priority):
 406         - on the path in OPENMM_PLUGIN_DIR environment variable if this is specified
 407         - on the path in the OPENMM_PLUGIN_DIR macro that is set by the build script
 408         - at the default location assumed by OpenMM
 409     */
 410         /* env var */
 411         char *pluginDir = getenv("OPENMM_PLUGIN_DIR");
 412         trim(pluginDir);
 413         /* no env var or empty */
 414         if (pluginDir != NULL && *pluginDir != '\0')
 415         {
 416             loadedPlugins = Platform::loadPluginsFromDirectory(pluginDir);
 417             if (loadedPlugins.size() > 0)
 418             {
 419                 hasLoadedPlugins = true;
 420                 usedPluginDir = pluginDir;
 421             }
 422             else
 423             {
 424                 gmx_fatal(FARGS, "The directory provided in the OPENMM_PLUGIN_DIR environment variable "
 425                     "(%s) does not contain valid OpenMM plugins. Check your OpenMM installation!", pluginDir);
 426             }
 427         }
 428
 429         /* macro set at build time  */
 430 #ifdef OPENMM_PLUGIN_DIR
 431         if (!hasLoadedPlugins)
 432         {
 433             loadedPlugins = Platform::loadPluginsFromDirectory(OPENMM_PLUGIN_DIR);
 434             if (loadedPlugins.size() > 0)
 435             {
 436                 hasLoadedPlugins = true;
 437                 usedPluginDir = OPENMM_PLUGIN_DIR;
 438             }
 439         }
 440 #endif
 441         /* default loocation */
 442         if (!hasLoadedPlugins)
 443         {
 444             loadedPlugins = Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
 445             if (loadedPlugins.size() > 0)
 446             {
 447                 hasLoadedPlugins = true;
 448                 usedPluginDir = Platform::getDefaultPluginsDirectory();
 449             }
 450         }
 451
 452         /* if there are still no plugins loaded there won't be any */
 453         if (!hasLoadedPlugins)
 454         {
 455             gmx_fatal(FARGS, "No OpenMM plugins were found! You can provide the"
 456                 " plugin directory in the OPENMM_PLUGIN_DIR environment variable.", pluginDir);
 457         }
 458
 459         fprintf(fplog, "\nPlugins loaded from directory %s:\t", usedPluginDir.c_str());
 460         for (int i = 0; i < loadedPlugins.size(); i++)
 461         {
 462             fprintf(fplog, "%s, ", loadedPlugins[i].c_str());
 463         }
 464         fprintf(fplog, "\n");
 465     }
 466
 467     /* parse option string */
 468     GmxOpenMMPlatformOptions *opt = new GmxOpenMMPlatformOptions(platformOptStr);
 469
 470     /* check wheter Gromacs options compatibility with OpenMM */
 471     checkGmxOptions(ir, top);
 472
 473     /* check GPU support */
 474     char s[1000];
 475     is_supported_cuda_gpu(atoi(opt->getOptionValue("deviceid").c_str()), s);
 476
 477     // Create the system.
 478     const t_idef& idef = top->idef;
 479     const int numAtoms = top_global->natoms;
 480     const int numConstraints = idef.il[F_CONSTR].nr/3;
 481     const int numSettle = idef.il[F_SETTLE].nr/2;
 482     const int numBonds = idef.il[F_BONDS].nr/3;
 483     const int numAngles = idef.il[F_ANGLES].nr/4;
 484     const int numPeriodic = idef.il[F_PDIHS].nr/5;
 485     const int numRB = idef.il[F_RBDIHS].nr/5;
 486     const int num14 = idef.il[F_LJ14].nr/3;
 487     System* sys = new System();
 488     if (ir->nstcomm > 0)
 489         sys->addForce(new CMMotionRemover(ir->nstcomm));
 490
 491     // Set bonded force field terms.
 492     const int* bondAtoms = (int*) idef.il[F_BONDS].iatoms;
 493     HarmonicBondForce* bondForce = new HarmonicBondForce();
 494     sys->addForce(bondForce);
 495     int offset = 0;
 496     for (int i = 0; i < numBonds; ++i) {
 497         int type = bondAtoms[offset++];
 498         int atom1 = bondAtoms[offset++];
 499         int atom2 = bondAtoms[offset++];
 500         bondForce->addBond(atom1, atom2, idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
 501     }
 502     const int* angleAtoms = (int*) idef.il[F_ANGLES].iatoms;
 503     HarmonicAngleForce* angleForce = new HarmonicAngleForce();
 504     sys->addForce(angleForce);
 505     offset = 0;
 506     for (int i = 0; i < numAngles; ++i) {
 507         int type = angleAtoms[offset++];
 508         int atom1 = angleAtoms[offset++];
 509         int atom2 = angleAtoms[offset++];
 510         int atom3 = angleAtoms[offset++];
 511         angleForce->addAngle(atom1, atom2, atom3, idef.iparams[type].harmonic.rA*M_PI/180.0, idef.iparams[type].harmonic.krA);
 512     }
 513     const int* periodicAtoms = (int*) idef.il[F_PDIHS].iatoms;
 514     PeriodicTorsionForce* periodicForce = new PeriodicTorsionForce();
 515     sys->addForce(periodicForce);
 516     offset = 0;
 517     for (int i = 0; i < numPeriodic; ++i) {
 518         int type = periodicAtoms[offset++];
 519         int atom1 = periodicAtoms[offset++];
 520         int atom2 = periodicAtoms[offset++];
 521         int atom3 = periodicAtoms[offset++];
 522         int atom4 = periodicAtoms[offset++];
 523         periodicForce->addTorsion(atom1, atom2, atom3, atom4, idef.iparams[type].pdihs.mult, idef.iparams[type].pdihs.phiA*M_PI/180.0, idef.iparams[type].pdihs.cpA);
 524     }
 525     const int* rbAtoms = (int*) idef.il[F_RBDIHS].iatoms;
 526     RBTorsionForce* rbForce = new RBTorsionForce();
 527     sys->addForce(rbForce);
 528     offset = 0;
 529     for (int i = 0; i < numRB; ++i) {
 530         int type = rbAtoms[offset++];
 531         int atom1 = rbAtoms[offset++];
 532         int atom2 = rbAtoms[offset++];
 533         int atom3 = rbAtoms[offset++];
 534         int atom4 = rbAtoms[offset++];
 535         rbForce->addTorsion(atom1, atom2, atom3, atom4, idef.iparams[type].rbdihs.rbcA[0], idef.iparams[type].rbdihs.rbcA[1],
 536                 idef.iparams[type].rbdihs.rbcA[2], idef.iparams[type].rbdihs.rbcA[3], idef.iparams[type].rbdihs.rbcA[4], idef.iparams[type].rbdihs.rbcA[5]);
 537     }
 538
 539     // Set nonbonded parameters and masses.
 540
 541     int ntypes = fr->ntype;
 542     int* types = mdatoms->typeA;
 543     real* nbfp = fr->nbfp;
 544     real* charges = mdatoms->chargeA;
 545     real* masses = mdatoms->massT;
 546     NonbondedForce* nonbondedForce = new NonbondedForce();
 547     sys->addForce(nonbondedForce);
 548
 549     if (ir->rcoulomb == 0)
 550             nonbondedForce->setNonbondedMethod(NonbondedForce::NoCutoff);
 551
 552     else {
 553             switch (ir->coulombtype)
 554         {
 555              case eelRF: // TODO what is the correct condition?
 556                      if (ir->ePBC == epbcXYZ) {
 557                          nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
 558                          sys->setPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0), Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));
 559                      }
 560                      else if (ir->ePBC == epbcNONE)
 561                          nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
 562                      else
 563                          gmx_fatal(FARGS,"** OpenMM Error ** : Only full periodic boundary conditions (pbc = xyz), or none (pbc = no) are supported\n");
 564                      nonbondedForce->setCutoffDistance(ir->rcoulomb);
 565                          break;
 566
 567              case eelEWALD:
 568                      if (ir->ewald_geometry == eewg3DC)
 569                          gmx_fatal(FARGS,"** OpenMM Error ** : Only Ewald 3D geometry is supported\n");
 570                      if (ir->epsilon_surface != 0)
 571                          gmx_fatal(FARGS,"** OpenMM Error ** : Dipole correction in Ewald summation is not supported\n");
 572                      nonbondedForce->setNonbondedMethod(NonbondedForce::Ewald);
 573                      nonbondedForce->setCutoffDistance(ir->rcoulomb);
 574                      sys->setPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0), Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));
 575                      break;
 576
 577                  case eelPME:
 578                      nonbondedForce->setNonbondedMethod(NonbondedForce::PME);
 579                      nonbondedForce->setCutoffDistance(ir->rcoulomb);
 580                      sys->setPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0), Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));
 581                      break;
 582
 583                  default:
 584                      gmx_fatal(FARGS,"Internal error: you should not see this message, it that the"
 585                     "electrosatics option check failed. Please report this error!");
 586                 }
 587     }
 588
 589     for (int i = 0; i < numAtoms; ++i) {
 590         real c6 = nbfp[types[i]*2*ntypes+types[i]*2];
 591         real c12 = nbfp[types[i]*2*ntypes+types[i]*2+1];
 592         if (c12 <= 0)
 593             nonbondedForce->addParticle(charges[i], 1.0, 0.0);
 594         else
 595             nonbondedForce->addParticle(charges[i], pow(c12/c6, (real) (1.0/6.0)), c6*c6/(4.0*c12));
 596         sys->addParticle(masses[i]);
 597     }
 598
 599     // Build a table of all exclusions.
 600
 601     vector<set<int> > exclusions(numAtoms);
 602     for (int i = 0; i < numAtoms; i++) {
 603         int start = top->excls.index[i];
 604         int end = top->excls.index[i+1];
 605         for (int j = start; j < end; j++)
 606             exclusions[i].insert(top->excls.a[j]);
 607     }
 608
 609     // Record the 1-4 interactions, and remove them from the list of exclusions.
 610
 611     const int* nb14Atoms = (int*) idef.il[F_LJ14].iatoms;
 612     offset = 0;
 613     for (int i = 0; i < num14; ++i) {
 614         int type = nb14Atoms[offset++];
 615         int atom1 = nb14Atoms[offset++];
 616         int atom2 = nb14Atoms[offset++];
 617         real c6 = idef.iparams[type].lj14.c6A;
 618         real c12 = idef.iparams[type].lj14.c12A;
 619         real sigma, epsilon;
 620         if (c12 <= 0) {
 621             epsilon = (real) 0.0;
 622             sigma = (real) 1.0;
 623         }
 624         else {
 625             epsilon = (real) ((c6*c6)/(4.0*c12));
 626             sigma = (real) pow(c12/c6, (real) (1.0/6.0));
 627         }
 628         nonbondedForce->addException(atom1, atom2, fr->fudgeQQ*charges[atom1]*charges[atom2], sigma, epsilon);
 629         exclusions[atom1].erase(atom2);
 630         exclusions[atom2].erase(atom1);
 631     }
 632
 633     // Record exclusions.
 634
 635     for (int i = 0; i < numAtoms; i++) {
 636         for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter) {
 637             if (i < *iter)
 638                 nonbondedForce->addException(i, *iter, 0.0, 1.0, 0.0);
 639         }
 640     }
 641
 642     // Add GBSA if needed.
 643     t_atoms        atoms;
 644     atoms    = gmx_mtop_global_atoms(top_global);
 645
 646     if (ir->implicit_solvent == eisGBSA) {
 647             GBSAOBCForce* gbsa = new GBSAOBCForce();
 648             sys->addForce(gbsa);
 649             gbsa->setSoluteDielectric(ir->epsilon_r);
 650             gbsa->setSolventDielectric(ir->gb_epsilon_solvent);
 651             gbsa->setCutoffDistance(nonbondedForce->getCutoffDistance());
 652             if (nonbondedForce->getNonbondedMethod() == NonbondedForce::NoCutoff)
 653                 gbsa->setNonbondedMethod(GBSAOBCForce::NoCutoff);
 654             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
 655                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
 656             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffPeriodic)
 657                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
 658             else
 659                         gmx_fatal(FARGS,"** OpenMM Error ** : Only Reaction-Field electrostatics is supported with OBC/GBSA.\n");
 660
 661             for (int i = 0; i < numAtoms; ++i)
 662                 gbsa->addParticle(charges[i],
 663                                   top_global->atomtypes.gb_radius[atoms.atom[i].type],
 664                                   top_global->atomtypes.S_hct[atoms.atom[i].type]);
 665     }
 666
 667     // Set constraints.
 668
 669     const int* constraintAtoms = (int*) idef.il[F_CONSTR].iatoms;
 670     offset = 0;
 671     for (int i = 0; i < numConstraints; ++i) {
 672         int type = constraintAtoms[offset++];
 673         int atom1 = constraintAtoms[offset++];
 674         int atom2 = constraintAtoms[offset++];
 675         sys->addConstraint(atom1, atom2, idef.iparams[type].constr.dA);
 676     }
 677     const int* settleAtoms = (int*) idef.il[F_SETTLE].iatoms;
 678     offset = 0;
 679     for (int i = 0; i < numSettle; ++i) {
 680         int type = settleAtoms[offset++];
 681         int oxygen = settleAtoms[offset++];
 682         sys->addConstraint(oxygen, oxygen+1, idef.iparams[type].settle.doh);
 683         sys->addConstraint(oxygen, oxygen+2, idef.iparams[type].settle.doh);
 684         sys->addConstraint(oxygen+1, oxygen+2, idef.iparams[type].settle.dhh);
 685     }
 686
 687     // Create an integrator for simulating the system.
 688
 689     real friction = (ir->opts.tau_t[0] == 0.0 ? 0.0 : 1.0/ir->opts.tau_t[0]);
 690     Integrator* integ;
 691     if (ir->eI == eiVV || ir->eI == eiVVAK) {
 692         integ = new VerletIntegrator(ir->delta_t);
 693         if ( ir->etc == etcANDERSEN) {
 694            real collisionFreq = ir->opts.tau_t[0] / 1000; /* tau_t (ps) / 1000 = collisionFreq (fs^-1) */
 695            AndersenThermostat* thermostat = new AndersenThermostat(ir->opts.ref_t[0], friction); /* TODO test this  */
 696            sys->addForce(thermostat);
 697         }
 698     }
 699
 700     else if (ir->eI == eiBD) {
 701         integ = new BrownianIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
 702         static_cast<BrownianIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed); /* TODO test this */
 703     }
 704     else if (EI_SD(ir->eI)) {
 705         integ = new LangevinIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
 706         static_cast<LangevinIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed); /* TODO test this */
 707     }
 708     else {
 709         gmx_fatal(FARGS, "** OpenMM Error ** : Unsupported integrator requested."
 710                 "Available integrators are: md-vv/md-vvak, sd/sd1, and bd. \n");
 711     }
 712
 713     integ->setConstraintTolerance(ir->shake_tol);
 714
 715     // Create a context and initialize it.
 716     Context* context = NULL;
 717     if (platformOptStr == NULL || platformOptStr == "")
 718     {
 719         context = new Context(*sys, *integ);
 720     }
 721     else
 722     {
 723         // Find which platform is it.
 724         for (int i = 0; i < Platform::getNumPlatforms() && context == NULL; i++)
 725         {
 726             if (isStringEqNCase(opt->getOptionValue("platform"), Platform::getPlatform(i).getName()))
 727             {
 728                 Platform& platform = Platform::getPlatform(i);
 729                 // set standard properties
 730                 platform.setPropertyDefaultValue("CudaDevice", opt->getOptionValue("deviceid"));
 731                 // TODO add extra properties
 732                 context = new Context(*sys, *integ, platform);
 733             }
 734         }
 735         if (context == NULL)
 736         {
 737             gmx_fatal(FARGS, "The requested platform \"%s\" could not be found.\n", opt->getOptionValue("platform").c_str());
 738         }
 739     }
 740
 741     Platform& platform = context->getPlatform();
 742     fprintf(fplog, "Gromacs will use the OpenMM platform: %s\n", platform.getName().c_str());
 743
 744     const vector<string>& properties = platform.getPropertyNames();
 745     if (debug)
 746     {
 747         for (int i = 0; i < properties.size(); i++)
 748         {
 749             printf(">> %s: %s\n", properties[i].c_str(), platform.getPropertyValue(*context, properties[i]).c_str());
 750             fprintf(fplog, ">> %s: %s\n", properties[i].c_str(), platform.getPropertyValue(*context, properties[i]).c_str());
 751         }
 752     }
 753
 754     int devId;
 755     if (!from_string<int>(devId, platform.getPropertyValue(*context, "CudaDevice"), std::dec))
 756     {
 757         gmx_fatal(FARGS, "Internal error: couldn't determine the device selected by OpenMM");
 758     }
 759
 760     /* check GPU compatibility */
 761     char gpuname[STRLEN];
 762     if (!is_supported_cuda_gpu(-1, gpuname))
 763     {
 764         if (!gmx_strcasecmp(opt->getOptionValue("force-device").c_str(), "yes"))
 765         {
 766             sprintf(warn_buf, "Non-supported GPU selected (#%d, %s), forced continuing.\n"
 767                     "Note, that the simulation can be slow or it migth even crash.", devId, gpuname);
 768             fprintf(fplog, "%s", warn_buf);
 769             warning_note(NULL);
 770         }
 771         else
 772         {
 773             gmx_fatal(FARGS, "The selected GPU (#%d, %s) is not supported by Gromacs! "
 774                     "Most probably you have a low-end GPU which would not perform well, or new hardware that"
 775                     "has not been tested yet with Gromacs-OpenMM. If you still want to try using this "
 776                     "device use the force=on option.", devId, gpuname);
 777         }
 778     }
 779     else
 780     {
 781         fprintf(fplog, "Gromacs will run on the GPU #%d (%s).\n", devId, gpuname);
 782     }
 783
 784     /* do the pre-simulation memtest */
 785     run_memtest(fplog, -1, "Pre", opt);
 786
 787     vector<Vec3> pos(numAtoms);
 788     vector<Vec3> vel(numAtoms);
 789     for (int i = 0; i < numAtoms; ++i) {
 790         pos[i] = Vec3(state->x[i][0], state->x[i][1], state->x[i][2]);
 791         vel[i] = Vec3(state->v[i][0], state->v[i][1], state->v[i][2]);
 792     }
 793     context->setPositions(pos);
 794     context->setVelocities(vel);
 795
 796     // Return a structure containing the system, integrator, and context.
 797     OpenMMData* data = new OpenMMData();
 798     data->system = sys;
 799     data->integrator = integ;
 800     data->context = context;
 801     data->removeCM = (ir->nstcomm > 0);
 802     data->platformOpt = opt;
 803     return data;
 804
 805     } catch (std::exception& e) {
 806         gmx_fatal(FARGS, "OpenMM exception caught while initializating: %s\n", e.what());
 807     }
 808 }
 809
 810 /**
 811  * Integrate one step.
 812  *
 813  * @param data the OpenMMData object created by openmm_init().
 814  */
 815 void openmm_take_one_step(void* data)
 816 {
 817     // static int step = 0; printf("----> taking step #%d\n", step++);
 818     try {
 819         static_cast<OpenMMData*>(data)->integrator->step(1);
 820     } catch (std::exception& e) {
 821         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s\n", e.what());
 822     }
 823 }
 824
 825 /**
 826  * Integrate n steps.
 827  *
 828  * @param data the OpenMMData object created by openmm_init().
 829  */
 830 void openmm_take_steps(void* data, int nstep)
 831 {
 832     try {
 833         static_cast<OpenMMData*>(data)->integrator->step(nstep);
 834     } catch (std::exception& e) {
 835         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s\n", e.what());
 836     }
 837 }
 838
 839 /**
 840  * Clean up all the data structures used by OpenMM.
 841  *
 842  * @param log file pointer
 843  * @param data the OpenMMData object created by openmm_init().
 844  */
 845 void openmm_cleanup(FILE* fplog, void* data)
 846 {
 847     OpenMMData* d = static_cast<OpenMMData*>(data);
 848     run_memtest(fplog, -1, "Post", d->platformOpt);
 849     delete d->system;
 850     delete d->integrator;
 851     delete d->context;
 852     delete d->platformOpt;
 853     delete d;
 854 }
 855
 856 /**
 857  * Copy the current state information (positions, velocities, and forces) into the data structures used
 858  * by Gromacs.
 859  *
 860  * @param data the OpenMMData object created by openmm_init().
 861  */
 862 void openmm_copy_state(void *data,
 863                         t_state *state, double *time,
 864                         rvec f[], gmx_enerdata_t *enerd,
 865                         bool includePos, bool includeVel, bool includeForce, bool includeEnergy)
 866 {
 867     int types = 0;
 868     if (includePos)
 869         types += State::Positions;
 870     if (includeVel)
 871         types += State::Velocities;
 872     if (includeForce)
 873         types += State::Forces;
 874     if (includeEnergy)
 875         types += State::Energy;
 876     if (types == 0)
 877         return;
 878     try {
 879         State currentState = static_cast<OpenMMData*>(data)->context->getState(types);
 880         int numAtoms =  static_cast<OpenMMData*>(data)->system->getNumParticles();
 881         if (includePos) {
 882             for (int i = 0; i < numAtoms; i++) {
 883                 Vec3 x = currentState.getPositions()[i];
 884                 state->x[i][0] = x[0];
 885                 state->x[i][1] = x[1];
 886                 state->x[i][2] = x[2];
 887             }
 888         }
 889         if (includeVel) {
 890             for (int i = 0; i < numAtoms; i++) {
 891                 Vec3 v = currentState.getVelocities()[i];
 892                 state->v[i][0] = v[0];
 893                 state->v[i][1] = v[1];
 894                 state->v[i][2] = v[2];
 895             }
 896         }
 897         if (includeForce) {
 898             for (int i = 0; i < numAtoms; i++) {
 899                 Vec3 force = currentState.getForces()[i];
 900                 f[i][0] = force[0];
 901                 f[i][1] = force[1];
 902                 f[i][2] = force[2];
 903             }
 904         }
 905         if (includeEnergy) {
 906             int numConstraints = static_cast<OpenMMData*>(data)->system->getNumConstraints();
 907             int dof = 3*numAtoms-numConstraints;
 908             if (static_cast<OpenMMData*>(data)->removeCM)
 909                 dof -= 3;
 910             enerd->term[F_EPOT] = currentState.getPotentialEnergy();
 911             enerd->term[F_EKIN] = currentState.getKineticEnergy();
 912             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
 913             enerd->term[F_TEMP] = 2.0*enerd->term[F_EKIN]/dof/BOLTZ;
 914         }
 915         *time = currentState.getTime();
 916     } catch (std::exception& e) {
 917         gmx_fatal(FARGS, "OpenMM exception caught while retrieving state information: %s\n", e.what());
 918     }
 919 }