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 "
 342                 "using either \n\t(1)\t\"md-vv\" or \"md-vvak\" integrators with \"andersen\" or "
 343                 "\"andersen-interval\" thermostat, or \n\t(2)\t\"sd\",\"sd1\" or \"bd\" integrators\n");
 344
 345     if (ir->opts.ngtc > 1)
 346         gmx_fatal(FARGS,"** OpenMM Error ** : Multiple temperature coupling groups are not supported. \n");
 347
 348     if (ir->epc != etcNO)
 349         gmx_fatal(FARGS,"** OpenMM Error ** : Pressure coupling is not supported. \n");
 350
 351     if (ir->opts.annealing[0])
 352         gmx_fatal(FARGS,"** OpenMM Error ** : Simulated annealing is not supported. \n");
 353
 354     if (ir->eConstrAlg != econtSHAKE)
 355         warning_note("** OpenMM Warning ** : Constraints in OpenMM are done by a combination "
 356                 "of SHAKE, SETTLE and CCMA. Accuracy is based on the SHAKE tolerance set "
 357                 "by the \"shake_tol\" option.");
 358
 359     if (ir->nwall != 0)
 360         gmx_fatal(FARGS,"** OpenMM Error ** : Walls are not supported.\n");
 361
 362     if (ir->ePull != epullNO)
 363         gmx_fatal(FARGS,"** OpenMM Error ** : Pulling is not supported.\n");
 364
 365     if (top->idef.il[F_DISRES].nr > 0)
 366         gmx_fatal(FARGS,"** OpenMM Error ** : Distant restraints are not supported.\n");
 367
 368     if (top->idef.il[F_ORIRES].nr > 0)
 369         gmx_fatal(FARGS,"** OpenMM Error ** : Orientation restraints are not supported.\n");
 370
 371     if (top->idef.il[F_ANGRES].nr > 0)
 372         gmx_fatal(FARGS,"** OpenMM Error ** : Angle restraints are not supported.\n");
 373
 374     if (top->idef.il[F_DIHRES].nr > 0)
 375         gmx_fatal(FARGS,"** OpenMM Error ** : Dihedral restraints are not supported.\n");
 376
 377     if (ir->efep != efepNO)
 378         gmx_fatal(FARGS,"** OpenMM Error ** : Free energy calculations are not supported.\n");
 379
 380     if (ir->opts.ngacc > 1)
 381         gmx_fatal(FARGS,"** OpenMM Error ** : Non-equilibrium MD (accelerated groups) are "
 382                 "not supported. \n");
 383
 384     if (IR_ELEC_FIELD(*ir))
 385         gmx_fatal(FARGS,"** OpenMM Error ** : Electric fields are not supported. \n");
 386
 387     if (ir->bQMMM)
 388         gmx_fatal(FARGS,"** OpenMM Error ** : QMMM calculations are not supported. \n");
 389
 390     if (ir->rcoulomb != ir->rvdw)
 391         gmx_fatal(FARGS,"** OpenMM Error ** : OpenMM uses a single cutoff for both Coulomb "
 392                 "and VdW interactions. Please set rcoulomb equal to rvdw. \n");
 393
 394 }
 395
 396 /**
 397  * Initialize OpenMM, and return a pointer to the OpenMMData object containing the System, Integrator, and Context.
 398  */
 399 void* openmm_init(FILE *fplog, const char *platformOptStr,
 400                     t_commrec *cr,t_inputrec *ir,
 401                     gmx_mtop_t *top_global, gmx_localtop_t *top,
 402                     t_mdatoms *mdatoms, t_forcerec *fr,t_state *state)
 403 {
 404     static bool hasLoadedPlugins = false;
 405     string usedPluginDir;
 406
 407     try {
 408     if (!hasLoadedPlugins)
 409     {
 410         vector<string> loadedPlugins;
 411     /*  Look for OpenMM plugins at various locations (listed in order of priority):
 412         - on the path in OPENMM_PLUGIN_DIR environment variable if this is specified
 413         - on the path in the OPENMM_PLUGIN_DIR macro that is set by the build script
 414         - at the default location assumed by OpenMM
 415     */
 416         /* env var */
 417         char *pluginDir = getenv("OPENMM_PLUGIN_DIR");
 418         trim(pluginDir);
 419         /* no env var or empty */
 420         if (pluginDir != NULL && *pluginDir != '\0')
 421         {
 422             loadedPlugins = Platform::loadPluginsFromDirectory(pluginDir);
 423             if (loadedPlugins.size() > 0)
 424             {
 425                 hasLoadedPlugins = true;
 426                 usedPluginDir = pluginDir;
 427             }
 428             else
 429             {
 430                 gmx_fatal(FARGS, "The directory provided in the OPENMM_PLUGIN_DIR environment variable "
 431                     "(%s) does not contain valid OpenMM plugins. Check your OpenMM installation!", pluginDir);
 432             }
 433         }
 434
 435         /* macro set at build time  */
 436 #ifdef OPENMM_PLUGIN_DIR
 437         if (!hasLoadedPlugins)
 438         {
 439             loadedPlugins = Platform::loadPluginsFromDirectory(OPENMM_PLUGIN_DIR);
 440             if (loadedPlugins.size() > 0)
 441             {
 442                 hasLoadedPlugins = true;
 443                 usedPluginDir = OPENMM_PLUGIN_DIR;
 444             }
 445         }
 446 #endif
 447         /* default loocation */
 448         if (!hasLoadedPlugins)
 449         {
 450             loadedPlugins = Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
 451             if (loadedPlugins.size() > 0)
 452             {
 453                 hasLoadedPlugins = true;
 454                 usedPluginDir = Platform::getDefaultPluginsDirectory();
 455             }
 456         }
 457
 458         /* if there are still no plugins loaded there won't be any */
 459         if (!hasLoadedPlugins)
 460         {
 461             gmx_fatal(FARGS, "No OpenMM plugins were found! You can provide the"
 462                 " plugin directory in the OPENMM_PLUGIN_DIR environment variable.", pluginDir);
 463         }
 464
 465         fprintf(fplog, "\nPlugins loaded from directory %s:\t", usedPluginDir.c_str());
 466         for (int i = 0; i < loadedPlugins.size(); i++)
 467         {
 468             fprintf(fplog, "%s, ", loadedPlugins[i].c_str());
 469         }
 470         fprintf(fplog, "\n");
 471     }
 472
 473     /* parse option string */
 474     GmxOpenMMPlatformOptions *opt = new GmxOpenMMPlatformOptions(platformOptStr);
 475
 476     /* check wheter Gromacs options compatibility with OpenMM */
 477     checkGmxOptions(ir, top);
 478
 479     /* check GPU support */
 480     char s[1000];
 481     is_supported_cuda_gpu(atoi(opt->getOptionValue("deviceid").c_str()), s);
 482
 483     // Create the system.
 484     const t_idef& idef = top->idef;
 485     const int numAtoms = top_global->natoms;
 486     const int numConstraints = idef.il[F_CONSTR].nr/3;
 487     const int numSettle = idef.il[F_SETTLE].nr/2;
 488     const int numBonds = idef.il[F_BONDS].nr/3;
 489     const int numAngles = idef.il[F_ANGLES].nr/4;
 490     const int numPeriodic = idef.il[F_PDIHS].nr/5;
 491     const int numRB = idef.il[F_RBDIHS].nr/5;
 492     const int num14 = idef.il[F_LJ14].nr/3;
 493     System* sys = new System();
 494     if (ir->nstcomm > 0)
 495         sys->addForce(new CMMotionRemover(ir->nstcomm));
 496
 497     // Set bonded force field terms.
 498     const int* bondAtoms = (int*) idef.il[F_BONDS].iatoms;
 499     HarmonicBondForce* bondForce = new HarmonicBondForce();
 500     sys->addForce(bondForce);
 501     int offset = 0;
 502     for (int i = 0; i < numBonds; ++i) {
 503         int type = bondAtoms[offset++];
 504         int atom1 = bondAtoms[offset++];
 505         int atom2 = bondAtoms[offset++];
 506         bondForce->addBond(atom1, atom2, idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
 507     }
 508     const int* angleAtoms = (int*) idef.il[F_ANGLES].iatoms;
 509     HarmonicAngleForce* angleForce = new HarmonicAngleForce();
 510     sys->addForce(angleForce);
 511     offset = 0;
 512     for (int i = 0; i < numAngles; ++i) {
 513         int type = angleAtoms[offset++];
 514         int atom1 = angleAtoms[offset++];
 515         int atom2 = angleAtoms[offset++];
 516         int atom3 = angleAtoms[offset++];
 517         angleForce->addAngle(atom1, atom2, atom3, idef.iparams[type].harmonic.rA*M_PI/180.0, idef.iparams[type].harmonic.krA);
 518     }
 519     const int* periodicAtoms = (int*) idef.il[F_PDIHS].iatoms;
 520     PeriodicTorsionForce* periodicForce = new PeriodicTorsionForce();
 521     sys->addForce(periodicForce);
 522     offset = 0;
 523     for (int i = 0; i < numPeriodic; ++i) {
 524         int type = periodicAtoms[offset++];
 525         int atom1 = periodicAtoms[offset++];
 526         int atom2 = periodicAtoms[offset++];
 527         int atom3 = periodicAtoms[offset++];
 528         int atom4 = periodicAtoms[offset++];
 529         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);
 530     }
 531     const int* rbAtoms = (int*) idef.il[F_RBDIHS].iatoms;
 532     RBTorsionForce* rbForce = new RBTorsionForce();
 533     sys->addForce(rbForce);
 534     offset = 0;
 535     for (int i = 0; i < numRB; ++i) {
 536         int type = rbAtoms[offset++];
 537         int atom1 = rbAtoms[offset++];
 538         int atom2 = rbAtoms[offset++];
 539         int atom3 = rbAtoms[offset++];
 540         int atom4 = rbAtoms[offset++];
 541         rbForce->addTorsion(atom1, atom2, atom3, atom4, idef.iparams[type].rbdihs.rbcA[0], idef.iparams[type].rbdihs.rbcA[1],
 542                 idef.iparams[type].rbdihs.rbcA[2], idef.iparams[type].rbdihs.rbcA[3], idef.iparams[type].rbdihs.rbcA[4], idef.iparams[type].rbdihs.rbcA[5]);
 543     }
 544
 545     // Set nonbonded parameters and masses.
 546
 547     int ntypes = fr->ntype;
 548     int* types = mdatoms->typeA;
 549     real* nbfp = fr->nbfp;
 550     real* charges = mdatoms->chargeA;
 551     real* masses = mdatoms->massT;
 552     NonbondedForce* nonbondedForce = new NonbondedForce();
 553     sys->addForce(nonbondedForce);
 554
 555     if (ir->rcoulomb == 0)
 556             nonbondedForce->setNonbondedMethod(NonbondedForce::NoCutoff);
 557
 558     else {
 559             switch (ir->coulombtype)
 560             {
 561                 case eelRF: // TODO what is the correct condition?
 562                     if (ir->ePBC == epbcXYZ) {
 563                          nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
 564                          sys->setPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0), Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));
 565                     }
 566                     else if (ir->ePBC == epbcNONE)
 567                         nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
 568                     else
 569                         gmx_fatal(FARGS,"** OpenMM Error ** : Only full periodic boundary conditions (pbc = xyz), or none (pbc = no) are supported\n");
 570                     nonbondedForce->setCutoffDistance(ir->rcoulomb);
 571                     break;
 572
 573              case eelEWALD:
 574                      if (ir->ewald_geometry == eewg3DC)
 575                          gmx_fatal(FARGS,"** OpenMM Error ** : Only Ewald 3D geometry is supported\n");
 576                      if (ir->epsilon_surface != 0)
 577                          gmx_fatal(FARGS,"** OpenMM Error ** : Dipole correction in Ewald summation is not supported\n");
 578                      nonbondedForce->setNonbondedMethod(NonbondedForce::Ewald);
 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                  case eelPME:
 584                      nonbondedForce->setNonbondedMethod(NonbondedForce::PME);
 585                      nonbondedForce->setCutoffDistance(ir->rcoulomb);
 586                      sys->setPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0), Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));
 587                      break;
 588
 589                  default:
 590                      gmx_fatal(FARGS,"Internal error: you should not see this message, it that the"
 591                     "electrosatics option check failed. Please report this error!");
 592                 }
 593     }
 594
 595     for (int i = 0; i < numAtoms; ++i) {
 596         real c6 = nbfp[types[i]*2*ntypes+types[i]*2];
 597         real c12 = nbfp[types[i]*2*ntypes+types[i]*2+1];
 598         if (c12 <= 0)
 599             nonbondedForce->addParticle(charges[i], 1.0, 0.0);
 600         else
 601             nonbondedForce->addParticle(charges[i], pow(c12/c6, (real) (1.0/6.0)), c6*c6/(4.0*c12));
 602         sys->addParticle(masses[i]);
 603     }
 604
 605     // Build a table of all exclusions.
 606
 607     vector<set<int> > exclusions(numAtoms);
 608     for (int i = 0; i < numAtoms; i++) {
 609         int start = top->excls.index[i];
 610         int end = top->excls.index[i+1];
 611         for (int j = start; j < end; j++)
 612             exclusions[i].insert(top->excls.a[j]);
 613     }
 614
 615     // Record the 1-4 interactions, and remove them from the list of exclusions.
 616
 617     const int* nb14Atoms = (int*) idef.il[F_LJ14].iatoms;
 618     offset = 0;
 619     for (int i = 0; i < num14; ++i) {
 620         int type = nb14Atoms[offset++];
 621         int atom1 = nb14Atoms[offset++];
 622         int atom2 = nb14Atoms[offset++];
 623         real c6 = idef.iparams[type].lj14.c6A;
 624         real c12 = idef.iparams[type].lj14.c12A;
 625         real sigma, epsilon;
 626         if (c12 <= 0) {
 627             epsilon = (real) 0.0;
 628             sigma = (real) 1.0;
 629         }
 630         else {
 631             epsilon = (real) ((c6*c6)/(4.0*c12));
 632             sigma = (real) pow(c12/c6, (real) (1.0/6.0));
 633         }
 634         nonbondedForce->addException(atom1, atom2, fr->fudgeQQ*charges[atom1]*charges[atom2], sigma, epsilon);
 635         exclusions[atom1].erase(atom2);
 636         exclusions[atom2].erase(atom1);
 637     }
 638
 639     // Record exclusions.
 640
 641     for (int i = 0; i < numAtoms; i++) {
 642         for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter) {
 643             if (i < *iter)
 644                 nonbondedForce->addException(i, *iter, 0.0, 1.0, 0.0);
 645         }
 646     }
 647
 648     // Add GBSA if needed.
 649     t_atoms        atoms;
 650     atoms    = gmx_mtop_global_atoms(top_global);
 651
 652     if (ir->implicit_solvent == eisGBSA) {
 653             GBSAOBCForce* gbsa = new GBSAOBCForce();
 654             sys->addForce(gbsa);
 655             gbsa->setSoluteDielectric(ir->epsilon_r);
 656             gbsa->setSolventDielectric(ir->gb_epsilon_solvent);
 657             gbsa->setCutoffDistance(nonbondedForce->getCutoffDistance());
 658             if (nonbondedForce->getNonbondedMethod() == NonbondedForce::NoCutoff)
 659                 gbsa->setNonbondedMethod(GBSAOBCForce::NoCutoff);
 660             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
 661                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
 662             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffPeriodic)
 663                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
 664             else
 665                         gmx_fatal(FARGS,"** OpenMM Error ** : Only Reaction-Field electrostatics is supported with OBC/GBSA.\n");
 666
 667             for (int i = 0; i < numAtoms; ++i)
 668                 gbsa->addParticle(charges[i],
 669                                   top_global->atomtypes.gb_radius[atoms.atom[i].type],
 670                                   top_global->atomtypes.S_hct[atoms.atom[i].type]);
 671     }
 672
 673     // Set constraints.
 674
 675     const int* constraintAtoms = (int*) idef.il[F_CONSTR].iatoms;
 676     offset = 0;
 677     for (int i = 0; i < numConstraints; ++i) {
 678         int type = constraintAtoms[offset++];
 679         int atom1 = constraintAtoms[offset++];
 680         int atom2 = constraintAtoms[offset++];
 681         sys->addConstraint(atom1, atom2, idef.iparams[type].constr.dA);
 682     }
 683     const int* settleAtoms = (int*) idef.il[F_SETTLE].iatoms;
 684     offset = 0;
 685     for (int i = 0; i < numSettle; ++i) {
 686         int type = settleAtoms[offset++];
 687         int oxygen = settleAtoms[offset++];
 688         sys->addConstraint(oxygen, oxygen+1, idef.iparams[type].settle.doh);
 689         sys->addConstraint(oxygen, oxygen+2, idef.iparams[type].settle.doh);
 690         sys->addConstraint(oxygen+1, oxygen+2, idef.iparams[type].settle.dhh);
 691     }
 692
 693     // Create an integrator for simulating the system.
 694
 695     real friction = (ir->opts.tau_t[0] == 0.0 ? 0.0 : 1.0/ir->opts.tau_t[0]);
 696     Integrator* integ;
 697     if (ir->eI == eiVV || ir->eI == eiVVAK) {
 698         integ = new VerletIntegrator(ir->delta_t);
 699         if ( ir->etc == etcANDERSEN) {
 700            real collisionFreq = ir->opts.tau_t[0] / 1000; /* tau_t (ps) / 1000 = collisionFreq (fs^-1) */
 701            AndersenThermostat* thermostat = new AndersenThermostat(ir->opts.ref_t[0], friction); /* TODO test this  */
 702            sys->addForce(thermostat);
 703         }
 704     }
 705
 706     else if (ir->eI == eiBD) {
 707         integ = new BrownianIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
 708         static_cast<BrownianIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed); /* TODO test this */
 709     }
 710     else if (EI_SD(ir->eI)) {
 711         integ = new LangevinIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
 712         static_cast<LangevinIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed); /* TODO test this */
 713     }
 714     else {
 715         gmx_fatal(FARGS, "** OpenMM Error ** : Unsupported integrator requested."
 716                 "Available integrators are: md-vv/md-vvak, sd/sd1, and bd. \n");
 717     }
 718
 719     integ->setConstraintTolerance(ir->shake_tol);
 720
 721     // Create a context and initialize it.
 722     Context* context = NULL;
 723     if (platformOptStr == NULL || platformOptStr == "")
 724     {
 725         context = new Context(*sys, *integ);
 726     }
 727     else
 728     {
 729         // Find which platform is it.
 730         for (int i = 0; i < Platform::getNumPlatforms() && context == NULL; i++)
 731         {
 732             if (isStringEqNCase(opt->getOptionValue("platform"), Platform::getPlatform(i).getName()))
 733             {
 734                 Platform& platform = Platform::getPlatform(i);
 735                 // set standard properties
 736                 platform.setPropertyDefaultValue("CudaDevice", opt->getOptionValue("deviceid"));
 737                 // TODO add extra properties
 738                 context = new Context(*sys, *integ, platform);
 739             }
 740         }
 741         if (context == NULL)
 742         {
 743             gmx_fatal(FARGS, "The requested platform \"%s\" could not be found.\n", opt->getOptionValue("platform").c_str());
 744         }
 745     }
 746
 747     Platform& platform = context->getPlatform();
 748     fprintf(fplog, "Gromacs will use the OpenMM platform: %s\n", platform.getName().c_str());
 749
 750     const vector<string>& properties = platform.getPropertyNames();
 751     if (debug)
 752     {
 753         for (int i = 0; i < properties.size(); i++)
 754         {
 755             printf(">> %s: %s\n", properties[i].c_str(), platform.getPropertyValue(*context, properties[i]).c_str());
 756             fprintf(fplog, ">> %s: %s\n", properties[i].c_str(), platform.getPropertyValue(*context, properties[i]).c_str());
 757         }
 758     }
 759
 760     int devId;
 761     if (!from_string<int>(devId, platform.getPropertyValue(*context, "CudaDevice"), std::dec))
 762     {
 763         gmx_fatal(FARGS, "Internal error: couldn't determine the device selected by OpenMM");
 764     }
 765
 766     /* check GPU compatibility */
 767     char gpuname[STRLEN];
 768     if (!is_supported_cuda_gpu(-1, gpuname))
 769     {
 770         if (!gmx_strcasecmp(opt->getOptionValue("force-device").c_str(), "yes"))
 771         {
 772             sprintf(warn_buf, "Non-supported GPU selected (#%d, %s), forced continuing.\n"
 773                     "Note, that the simulation can be slow or it migth even crash.", devId, gpuname);
 774             fprintf(fplog, "%s", warn_buf);
 775             warning_note(NULL);
 776         }
 777         else
 778         {
 779             gmx_fatal(FARGS, "The selected GPU (#%d, %s) is not supported by Gromacs! "
 780                     "Most probably you have a low-end GPU which would not perform well, or new hardware that"
 781                     "has not been tested yet with Gromacs-OpenMM. If you still want to try using this "
 782                     "device use the force=on option.", devId, gpuname);
 783         }
 784     }
 785     else
 786     {
 787         fprintf(fplog, "Gromacs will run on the GPU #%d (%s).\n", devId, gpuname);
 788     }
 789
 790     /* do the pre-simulation memtest */
 791     run_memtest(fplog, -1, "Pre", opt);
 792
 793     vector<Vec3> pos(numAtoms);
 794     vector<Vec3> vel(numAtoms);
 795     for (int i = 0; i < numAtoms; ++i) {
 796         pos[i] = Vec3(state->x[i][0], state->x[i][1], state->x[i][2]);
 797         vel[i] = Vec3(state->v[i][0], state->v[i][1], state->v[i][2]);
 798     }
 799     context->setPositions(pos);
 800     context->setVelocities(vel);
 801
 802     // Return a structure containing the system, integrator, and context.
 803     OpenMMData* data = new OpenMMData();
 804     data->system = sys;
 805     data->integrator = integ;
 806     data->context = context;
 807     data->removeCM = (ir->nstcomm > 0);
 808     data->platformOpt = opt;
 809     return data;
 810
 811     } catch (std::exception& e) {
 812         gmx_fatal(FARGS, "OpenMM exception caught while initializating: %s\n", e.what());
 813     }
 814 }
 815
 816 /**
 817  * Integrate one step.
 818  *
 819  * @param data the OpenMMData object created by openmm_init().
 820  */
 821 void openmm_take_one_step(void* data)
 822 {
 823     // static int step = 0; printf("----> taking step #%d\n", step++);
 824     try {
 825         static_cast<OpenMMData*>(data)->integrator->step(1);
 826     } catch (std::exception& e) {
 827         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s\n", e.what());
 828     }
 829 }
 830
 831 /**
 832  * Integrate n steps.
 833  *
 834  * @param data the OpenMMData object created by openmm_init().
 835  */
 836 void openmm_take_steps(void* data, int nstep)
 837 {
 838     try {
 839         static_cast<OpenMMData*>(data)->integrator->step(nstep);
 840     } catch (std::exception& e) {
 841         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s\n", e.what());
 842     }
 843 }
 844
 845 /**
 846  * Clean up all the data structures used by OpenMM.
 847  *
 848  * @param log file pointer
 849  * @param data the OpenMMData object created by openmm_init().
 850  */
 851 void openmm_cleanup(FILE* fplog, void* data)
 852 {
 853     OpenMMData* d = static_cast<OpenMMData*>(data);
 854     run_memtest(fplog, -1, "Post", d->platformOpt);
 855     delete d->system;
 856     delete d->integrator;
 857     delete d->context;
 858     delete d->platformOpt;
 859     delete d;
 860 }
 861
 862 /**
 863  * Copy the current state information (positions, velocities, and forces) into the data structures used
 864  * by Gromacs.
 865  *
 866  * @param data the OpenMMData object created by openmm_init().
 867  */
 868 void openmm_copy_state(void *data,
 869                         t_state *state, double *time,
 870                         rvec f[], gmx_enerdata_t *enerd,
 871                         bool includePos, bool includeVel, bool includeForce, bool includeEnergy)
 872 {
 873     int types = 0;
 874     if (includePos)
 875         types += State::Positions;
 876     if (includeVel)
 877         types += State::Velocities;
 878     if (includeForce)
 879         types += State::Forces;
 880     if (includeEnergy)
 881         types += State::Energy;
 882     if (types == 0)
 883         return;
 884     try {
 885         State currentState = static_cast<OpenMMData*>(data)->context->getState(types);
 886         int numAtoms =  static_cast<OpenMMData*>(data)->system->getNumParticles();
 887         if (includePos) {
 888             for (int i = 0; i < numAtoms; i++) {
 889                 Vec3 x = currentState.getPositions()[i];
 890                 state->x[i][0] = x[0];
 891                 state->x[i][1] = x[1];
 892                 state->x[i][2] = x[2];
 893             }
 894         }
 895         if (includeVel) {
 896             for (int i = 0; i < numAtoms; i++) {
 897                 Vec3 v = currentState.getVelocities()[i];
 898                 state->v[i][0] = v[0];
 899                 state->v[i][1] = v[1];
 900                 state->v[i][2] = v[2];
 901             }
 902         }
 903         if (includeForce) {
 904             for (int i = 0; i < numAtoms; i++) {
 905                 Vec3 force = currentState.getForces()[i];
 906                 f[i][0] = force[0];
 907                 f[i][1] = force[1];
 908                 f[i][2] = force[2];
 909             }
 910         }
 911         if (includeEnergy) {
 912             int numConstraints = static_cast<OpenMMData*>(data)->system->getNumConstraints();
 913             int dof = 3*numAtoms-numConstraints;
 914             if (static_cast<OpenMMData*>(data)->removeCM)
 915                 dof -= 3;
 916             enerd->term[F_EPOT] = currentState.getPotentialEnergy();
 917             enerd->term[F_EKIN] = currentState.getKineticEnergy();
 918             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
 919             enerd->term[F_TEMP] = 2.0*enerd->term[F_EKIN]/dof/BOLTZ;
 920         }
 921         *time = currentState.getTime();
 922     } catch (std::exception& e) {
 923         gmx_fatal(FARGS, "OpenMM exception caught while retrieving state information: %s\n", e.what());
 924     }
 925 }