src/kernel/openmm_wrapper.cpp

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2010, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35
  36 /*
  37  * Note, that parts of this source code originate from the Simtk release
  38  * of OpenMM accelerated Gromacs, for more details see:
  39  * https://simtk.org/project/xml/downloads.xml?group_id=161#package_id600
  40  */
  41
  42 #ifdef HAVE_CONFIG_H
  43 #include <config.h>
  44 #endif
  45
  46 #include <types/simple.h>
  47 #include <cmath>
  48 #include <set>
  49 #include <iostream>
  50 #include <sstream>
  51 #include <fstream>
  52 #include <map>
  53 #include <vector>
  54 #include <cctype>
  55 #include <algorithm>
  56
  57 using namespace std;
  58
  59 #include "OpenMM.h"
  60
  61 #include "gmx_fatal.h"
  62 #include "typedefs.h"
  63 #include "mdrun.h"
  64 #include "physics.h"
  65 #include "string2.h"
  66 #include "gpu_utils.h"
  67 #include "mtop_util.h"
  68
  69 #include "openmm_wrapper.h"
  70
  71 using namespace OpenMM;
  72
  73 /*! \cond */
  74 #define MEM_ERR_MSG(str) \
  75     "The %s-simulation GPU memory test detected errors. As memory errors would cause incorrect " \
  76     "simulation results, gromacs has aborted execution.\n Make sure that your GPU's memory is not " \
  77     "overclocked and that the device is properly cooled.\n", (str)
  78 /*! \endcond */
  79
  80 #define COMBRULE_CHK_TOL            1e-6
  81 #define COMBRULE_SIGMA(sig1, sig2)  (((sig1) + (sig2))/2)
  82 #define COMBRULE_EPS(eps1, eps2)    (sqrt((eps1) * (eps2)))
  83
  84 /*!
  85  * \brief Convert string to integer type.
  86  * \param[in]  s    String to convert from.
  87  * \param[in]  f    Basefield format flag that takes any of the following I/O
  88  *                  manipulators: dec, hex, oct.
  89  * \param[out] t    Destination variable to convert to.
  90  */
  91 template <class T>
  92 static gmx_bool from_string(T& t, const string& s, ios_base& (*f)(ios_base&))
  93 {
  94     istringstream iss(s);
  95     return !(iss >> f >> t).fail();
  96 }
  97
  98 /*!
  99  * \brief Split string around a given delimiter.
 100  * \param[in] s      String to split.
 101  * \param[in] delim  Delimiter character.
 102  * \returns          Vector of strings found in \p s.
 103  */
 104 static vector<string> split(const string &s, char delim)
 105 {
 106     vector<string> elems;
 107     stringstream ss(s);
 108     string item;
 109     while (getline(ss, item, delim))
 110     {
 111         if (item.length() != 0)
 112             elems.push_back(item);
 113     }
 114     return elems;
 115 }
 116
 117 /*!
 118  * \brief Split a string of the form "option=value" into "option" and "value" strings.
 119  * This string corresponds to one option and the associated value from the option list
 120  * in the mdrun -device argument.
 121  *
 122  * \param[in]  s    A string containing an "option=value" pair that needs to be split up.
 123  * \param[out] opt  The name of the option.
 124  * \param[out] val  Value of the option.
 125  */
 126 static void splitOptionValue(const string &s, string &opt, string &val)
 127 {
 128     size_t eqPos = s.find('=');
 129     if (eqPos != string::npos)
 130     {
 131         opt = s.substr(0, eqPos);
 132         if (eqPos != s.length())  val = s.substr(eqPos+1);
 133     }
 134 }
 135
 136 /*!
 137  * \brief Compare two strings ignoring case.
 138  * This function is in fact a wrapper around the gromacs function gmx_strncasecmp().
 139  * \param[in] s1 String.
 140  * \param[in] s2 String.
 141  * \returns      Similarly to the C function strncasecmp(), the return value is an
 142                  integer less than, equal to, or greater than 0 if \p s1 less than,
 143                  identical to, or greater than \p s2.
 144  */
 145 static gmx_bool isStringEqNCase(const string& s1, const string& s2)
 146 {
 147     return (gmx_strncasecmp(s1.c_str(), s2.c_str(), max(s1.length(), s2.length())) == 0);
 148 }
 149
 150 /*!
 151  * \brief Convert string to upper case.
 152  *
 153  * \param[in]  s    String to convert to uppercase.
 154  * \returns         The given string converted to uppercase.
 155  */
 156 static string toUpper(const string &s)
 157 {
 158     string stmp(s);
 159     std::transform(stmp.begin(), stmp.end(), stmp.begin(), static_cast < int(*)(int) > (toupper));
 160     return stmp;
 161 }
 162
 163 /*!
 164   \name Sizes of constant device option arrays GmxOpenMMPlatformOptions#platforms,
 165   GmxOpenMMPlatformOptions#memtests, GmxOpenMMPlatformOptions#deviceid,
 166   GmxOpenMMPlatformOptions#force_dev.  */
 167 /* {@ */
 168 #define SIZEOF_PLATFORMS    2  // 2
 169 #define SIZEOF_MEMTESTS     3
 170 #define SIZEOF_DEVICEIDS    1
 171 #define SIZEOF_FORCE_DEV    2
 172
 173 #define SIZEOF_CHECK_COMBRULE 2
 174 /* @} */
 175
 176 /*! Possible platform options in the mdrun -device option. */
 177 static const char *devOptStrings[] = { "platform", "deviceid", "memtest", "force-device", "check-combrule" };
 178
 179 /*! Enumerated platform options in the mdrun -device option. */
 180 enum devOpt
 181 {
 182     PLATFORM     = 0,
 183     DEVICEID     = 1,
 184     MEMTEST      = 2,
 185     FORCE_DEVICE = 3
 186 };
 187
 188 /*!
 189  * \brief Class to extract and manage the platform options in the mdrun -device option.
 190  *
 191  */
 192 class GmxOpenMMPlatformOptions
 193 {
 194 public:
 195     GmxOpenMMPlatformOptions(const char *opt);
 196     ~GmxOpenMMPlatformOptions() { options.clear(); }
 197     string getOptionValue(const string &opt);
 198     void remOption(const string &opt);
 199     void print();
 200 private:
 201     void setOption(const string &opt, const string &val);
 202
 203     map<string, string> options; /*!< Data structure to store the option (name, value) pairs. */
 204
 205     static const char * const platforms[SIZEOF_PLATFORMS];  /*!< Available OpenMM platforms; size #SIZEOF_PLATFORMS */
 206     static const char * const memtests[SIZEOF_MEMTESTS];    /*!< Available types of memory tests, also valid
 207                                                                  any positive integer >=15; size #SIZEOF_MEMTESTS */
 208     static const char * const deviceid[SIZEOF_DEVICEIDS];   /*!< Possible values for deviceid option;
 209                                                                  also valid any positive integer; size #SIZEOF_DEVICEIDS */
 210     static const char * const force_dev[SIZEOF_FORCE_DEV];  /*!< Possible values for for force-device option;
 211                                                                  size #SIZEOF_FORCE_DEV */
 212     static const char * const check_combrule[SIZEOF_CHECK_COMBRULE]; /* XXX temporary debug feature to
 213                                                                       turn off combination rule check */
 214 };
 215
 216 const char * const GmxOpenMMPlatformOptions::platforms[SIZEOF_PLATFORMS]
 217                     = {"CUDA", "Reference"};
 218                     //= { "Reference", "CUDA" /*,"OpenCL"*/ };
 219 const char * const GmxOpenMMPlatformOptions::memtests[SIZEOF_MEMTESTS]
 220                     = { "15", "full", "off" };
 221 const char * const GmxOpenMMPlatformOptions::deviceid[SIZEOF_DEVICEIDS]
 222                     = { "0" };
 223 const char * const GmxOpenMMPlatformOptions::force_dev[SIZEOF_FORCE_DEV]
 224                     = { "no", "yes" };
 225 const char * const GmxOpenMMPlatformOptions::check_combrule[SIZEOF_CHECK_COMBRULE]
 226                     = { "yes", "no" };
 227
 228 /*!
 229  * \brief Contructor.
 230  * Takes the option list, parses it, checks the options and their values for validity.
 231  * When certain options are not provided by the user, as default value the first item
 232  * of the respective constant array is taken (GmxOpenMMPlatformOptions#platforms,
 233  * GmxOpenMMPlatformOptions#memtests, GmxOpenMMPlatformOptions#deviceid,
 234  * GmxOpenMMPlatformOptions#force_dev).
 235  * \param[in] optionString  Option list part of the mdrun -device parameter.
 236  */
 237 GmxOpenMMPlatformOptions::GmxOpenMMPlatformOptions(const char *optionString)
 238 {
 239     // set default values
 240     setOption("platform",       platforms[0]);
 241     setOption("memtest",        memtests[0]);
 242     setOption("deviceid",       deviceid[0]);
 243     setOption("force-device",   force_dev[0]);
 244     setOption("check-combrule", check_combrule[0]);
 245
 246     string opt(optionString);
 247
 248     // remove all whitespaces
 249     opt.erase(remove_if(opt.begin(), opt.end(), ::isspace), opt.end());
 250     // tokenize around ","-s
 251     vector<string> tokens = split(opt, ',');
 252
 253     for (vector<string>::iterator it = tokens.begin(); it != tokens.end(); ++it)
 254     {
 255         string opt = "", val = "";
 256         splitOptionValue(*it, opt, val);
 257
 258         if (isStringEqNCase(opt, "platform"))
 259         {
 260             /* no check, this will fail if platform does not exist when we try to set it */
 261             setOption(opt, val);
 262             continue;
 263         }
 264
 265         if (isStringEqNCase(opt, "memtest"))
 266         {
 267             /* the value has to be an integer >15(s) or "full" OR "off" */
 268             if (!isStringEqNCase(val, "full") && !isStringEqNCase(val, "off"))
 269             {
 270                 int secs;
 271                 if (!from_string<int>(secs, val, std::dec))
 272                 {
 273                     gmx_fatal(FARGS, "Invalid value for option memtest option: \"%s\"!", val.c_str());
 274                 }
 275                 if (secs < 15)
 276                 {
 277                     gmx_fatal(FARGS, "Incorrect value for memtest option (%d). "
 278                             "Memtest needs to run for at least 15s!", secs);
 279                 }
 280             }
 281             setOption(opt, val);
 282             continue;
 283         }
 284
 285         if (isStringEqNCase(opt, "deviceid"))
 286         {
 287             int id;
 288             if (!from_string<int>(id, val, std::dec) )
 289             {
 290                 gmx_fatal(FARGS, "Invalid device id: \"%s\"!", val.c_str());
 291             }
 292             setOption(opt, val);
 293             continue;
 294         }
 295
 296         if (isStringEqNCase(opt, "force-device"))
 297         {
 298             /* */
 299             if (!isStringEqNCase(val, "yes") && !isStringEqNCase(val, "no"))
 300             {
 301                 gmx_fatal(FARGS, "Invalid OpenMM force option: \"%s\"!", val.c_str());
 302             }
 303             setOption(opt, val);
 304             continue;
 305         }
 306
 307         if (isStringEqNCase(opt, "check-combrule"))
 308         {
 309             /* */
 310             if (!isStringEqNCase(val, "yes") && !isStringEqNCase(val, "no"))
 311             {
 312                 gmx_fatal(FARGS, "Invalid OpenMM force option: \"%s\"!", val.c_str());
 313             }
 314             setOption(opt, val);
 315             continue;
 316         }
 317
 318
 319         // if we got till here something went wrong
 320         gmx_fatal(FARGS, "Invalid OpenMM platform option: \"%s\"!", (*it).c_str());
 321     }
 322 }
 323
 324
 325 /*!
 326  * \brief Getter function.
 327  * \param[in] opt   Name of the option.
 328  * \returns         Returns the value associated to an option.
 329  */
 330 string GmxOpenMMPlatformOptions::getOptionValue(const string &opt)
 331 {
 332         map<string, string> :: const_iterator it = options.find(toUpper(opt));
 333         if (it != options.end())
 334     {
 335                 return it->second;
 336     }
 337     else
 338     {
 339         return NULL;
 340     }
 341 }
 342
 343 /*!
 344  * \brief Setter function - private, only used from contructor.
 345  * \param[in] opt   Name of the option.
 346  * \param[in] val   Value for the option.
 347  */
 348 void GmxOpenMMPlatformOptions::setOption(const string &opt, const string &val)
 349 {
 350     options[toUpper(opt)] = val;
 351 }
 352
 353 /*!
 354  * \brief Removes an option with its value from the map structure. If the option
 355  * does not exist, returns without any action.
 356  * \param[in] opt   Name of the option.
 357  */
 358 void GmxOpenMMPlatformOptions::remOption(const string &opt)
 359 {
 360     options.erase(toUpper(opt));
 361 }
 362
 363 /*!
 364  * \brief Print option-value pairs to a file (debugging function).
 365  */
 366 void GmxOpenMMPlatformOptions::print()
 367 {
 368     cout << ">> Platform options: " << endl
 369          << ">> platform     = " << getOptionValue("platform") << endl
 370          << ">> deviceID     = " << getOptionValue("deviceid") << endl
 371          << ">> memtest      = " << getOptionValue("memtest") << endl
 372          << ">> force-device = " << getOptionValue("force-device") << endl;
 373 }
 374
 375 /*!
 376  * \brief Container for OpenMM related data structures that represent the bridge
 377  *        between the Gromacs data-structures and the OpenMM library and is but it's
 378  *        only passed through the API functions as void to disable direct access.
 379  */
 380 class OpenMMData
 381 {
 382 public:
 383     System* system;      /*! The system to simulate. */
 384     Context* context;   /*! The OpenMM context in which the simulation is carried out. */
 385     Integrator* integrator; /*! The integrator used in the simulation. */
 386     gmx_bool removeCM;          /*! If \true remove venter of motion, false otherwise. */
 387     GmxOpenMMPlatformOptions *platformOpt; /*! Platform options. */
 388 };
 389
 390 /*!
 391  *  \brief Runs memtest on the GPU that has alreaby been initialized by OpenMM.
 392  *  \param[in] fplog    Pointer to gromacs log file.
 393  *  \param[in] devId    Device id of the GPU to run the test on.
 394                         Note: as OpenMM previously creates the context,for now this is always -1.
 395  *  \param[in] pre_post Contains either "Pre" or "Post" just to be able to differentiate in
 396  *                      stdout messages/log between memtest carried out before and after simulation.
 397  *  \param[in] opt      Pointer to platform options object.
 398  */
 399 static void runMemtest(FILE* fplog, int devId, const char* pre_post, GmxOpenMMPlatformOptions *opt)
 400 {
 401     char        strout_buf[STRLEN];
 402     int         which_test;
 403     int         res = 0;
 404     string      s = opt->getOptionValue("memtest");
 405     const char  *test_type = s.c_str();
 406
 407     if (!gmx_strcasecmp(test_type, "off"))
 408     {
 409         which_test = 0;
 410     }
 411     else
 412     {
 413         if (!gmx_strcasecmp(test_type, "full"))
 414         {
 415             which_test = 2;
 416         }
 417         else
 418         {
 419             from_string<int>(which_test, test_type, std::dec);
 420         }
 421     }
 422
 423     if (which_test < 0)
 424     {
 425         gmx_fatal(FARGS, "Amount of seconds for memetest is negative (%d). ", which_test);
 426     }
 427
 428     switch (which_test)
 429     {
 430         case 0: /* no memtest */
 431             sprintf(strout_buf, "%s-simulation GPU memtest skipped. Note, that faulty memory can cause "
 432                 "incorrect results!", pre_post);
 433             fprintf(fplog, "%s\n", strout_buf);
 434             gmx_warning(strout_buf);
 435             break; /* case 0 */
 436
 437         case 1: /* quick memtest */
 438             fprintf(fplog,  "%s-simulation %s GPU memtest in progress...\n", pre_post, test_type);
 439             fprintf(stdout, "\n%s-simulation %s GPU memtest in progress...", pre_post, test_type);
 440             fflush(fplog);
 441             fflush(stdout);
 442             res = do_quick_memtest(devId);
 443             break; /* case 1 */
 444
 445         case 2: /* full memtest */
 446             fprintf(fplog,  "%s-simulation %s memtest in progress...\n", pre_post, test_type);
 447             fprintf(stdout, "\n%s-simulation %s memtest in progress...", pre_post, test_type);
 448             fflush(fplog);
 449             fflush(stdout);
 450             res = do_full_memtest(devId);
 451             break; /* case 2 */
 452
 453         default: /* timed memtest */
 454             fprintf(fplog,  "%s-simulation ~%ds memtest in progress...\n", pre_post, which_test);
 455             fprintf(stdout, "\n%s-simulation ~%ds memtest in progress...", pre_post, which_test);
 456             fflush(fplog);
 457             fflush(stdout);
 458             res = do_timed_memtest(devId, which_test);
 459         }
 460
 461         if (which_test != 0)
 462         {
 463             if (res != 0)
 464             {
 465                 gmx_fatal(FARGS, MEM_ERR_MSG(pre_post));
 466             }
 467             else
 468             {
 469                 fprintf(fplog,  "Memory test completed without errors.\n");
 470                 fflush(fplog);
 471                 fprintf(stdout, "done, no errors detected\n");
 472                 fflush(stdout);
 473             }
 474         }
 475 }
 476
 477 /*!
 478  * \brief Convert Lennard-Jones parameters c12 and c6 to sigma and epsilon.
 479  *
 480  * \param[in] c12
 481  * \param[in] c6
 482  * \param[out] sigma
 483  * \param[out] epsilon
 484  */
 485 static void convert_c_12_6(double c12, double c6, double *sigma, double *epsilon)
 486 {
 487     if (c12 == 0 && c6 == 0)
 488     {
 489         *epsilon    = 0.0;
 490         *sigma      = 1.0;
 491     }
 492     else if (c12 > 0 && c6 > 0)
 493     {
 494         *epsilon    = (c6*c6)/(4.0*c12);
 495         *sigma      = pow(c12/c6, 1.0/6.0);
 496     }
 497     else
 498     {
 499         gmx_fatal(FARGS,"OpenMM only supports c6 > 0 and c12 > 0 or c6 = c12 = 0.");
 500     }
 501 }
 502
 503 /*!
 504  * \brief Does gromacs option checking.
 505  *
 506  * Checks the gromacs mdp options for features unsupported in OpenMM, case in which
 507  * interrupts the execution. It also warns the user about pecularities of OpenMM
 508  * implementations.
 509  * \param[in] fplog         Gromacs log file pointer.
 510  * \param[in] ir            Gromacs input parameters, see ::t_inputrec
 511  * \param[in] top           Gromacs node local topology, \see gmx_localtop_t
 512  * \param[in] state         Gromacs state structure \see ::t_state
 513  * \param[in] mdatoms       Gromacs atom parameters, \see ::t_mdatoms
 514  * \param[in] fr            \see ::t_forcerec
 515  * \param[in] state         Gromacs systems state, \see ::t_state
 516  */
 517 static void checkGmxOptions(FILE* fplog, GmxOpenMMPlatformOptions *opt,
 518                             t_inputrec *ir, gmx_localtop_t *top,
 519                             t_forcerec *fr, t_state *state)
 520 {
 521     char    warn_buf[STRLEN];
 522     int     i, j, natoms;
 523     double  c6, c12;
 524     double  sigma_ij=0, sigma_ji=0, sigma_ii=0, sigma_jj=0, sigma_comb;
 525     double  eps_ij=0, eps_ji=0, eps_ii=0, eps_jj=0, eps_comb;
 526
 527     /* Abort if unsupported critical options are present */
 528
 529     /* Integrator */
 530     if (ir->eI ==  eiMD)
 531     {
 532         gmx_warning( "OpenMM does not support leap-frog, will use velocity-verlet integrator.");
 533     }
 534
 535     if (    (ir->eI !=  eiMD)   &&
 536             (ir->eI !=  eiVV)   &&
 537             (ir->eI !=  eiVVAK) &&
 538             (ir->eI !=  eiSD1)  &&
 539             (ir->eI !=  eiSD2)  &&
 540             (ir->eI !=  eiBD) )
 541     {
 542         gmx_fatal(FARGS, "OpenMM supports only the following integrators: md/md-vv/md-vv-avek, sd/sd1, and bd.");
 543     }
 544
 545     /* Electroctstics */
 546     if (   !(ir->coulombtype == eelPME   ||
 547              EEL_RF(ir->coulombtype)     ||
 548              ir->coulombtype == eelRF    ||
 549              ir->coulombtype == eelEWALD ||
 550              // no-cutoff
 551              (ir->coulombtype == eelCUT && ir->rcoulomb == 0 &&  ir->rvdw == 0) ||
 552              // we could have cut-off combined with GBSA (openmm will use RF)
 553              ir->implicit_solvent == eisGBSA)   )
 554     {
 555         gmx_fatal(FARGS,"OpenMM supports only the following methods for electrostatics: "
 556                 "NoCutoff (i.e. rcoulomb = rvdw = 0 ),Reaction-Field, Ewald or PME.");
 557     }
 558
 559     if (EEL_RF(ir->coulombtype) && ir->epsilon_rf != 0)
 560     {
 561         // openmm has epsilon_rf=inf hard-coded
 562         gmx_warning("OpenMM will use a Reaction-Field epsilon of infinity instead of %g.",ir->epsilon_rf);
 563     }
 564
 565     if (ir->etc != etcNO &&
 566         ir->eI  != eiSD1 &&
 567         ir->eI  != eiSD2 &&
 568         ir->eI  != eiBD )
 569     {
 570         gmx_warning("OpenMM supports only Andersen thermostat with the md/md-vv/md-vv-avek integrators.");
 571     }
 572
 573     if (ir->implicit_solvent == eisGBSA &&
 574         ir->gb_algorithm != egbOBC  )
 575     {
 576         gmx_warning("OpenMM does not support the specified algorithm for Generalized Born, will use OBC instead.");
 577     }
 578
 579     if (ir->opts.ngtc > 1)
 580         gmx_fatal(FARGS,"OpenMM does not support multiple temperature coupling groups.");
 581
 582     if (ir->epc != epcNO)
 583         gmx_warning("OpenMM supports only Monte Carlo barostat for pressure coupling.");
 584
 585     if (ir->opts.annealing[0])
 586         gmx_fatal(FARGS,"OpenMM does not support simulated annealing.");
 587
 588     if (top->idef.il[F_CONSTR].nr > 0 && ir->eConstrAlg != econtSHAKE)
 589         gmx_warning("OpenMM provides contraints as a combination "
 590                     "of SHAKE, SETTLE and CCMA. Accuracy is based on the SHAKE tolerance set "
 591                     "by the \"shake_tol\" option.");
 592
 593     if (ir->nwall != 0)
 594         gmx_fatal(FARGS,"OpenMM does not support walls.");
 595
 596     if (ir->ePull != epullNO)
 597         gmx_fatal(FARGS,"OpenMM does not support pulling.");
 598
 599     /* check for interaction types */
 600     for (i = 0; i < F_EPOT; i++)
 601     {
 602         if (!(i == F_CONSTR ||
 603             i == F_SETTLE   ||
 604             i == F_BONDS    ||
 605             i == F_HARMONIC ||
 606             i == F_UREY_BRADLEY ||
 607             i == F_ANGLES   ||
 608             i == F_PDIHS    ||
 609             i == F_RBDIHS   ||
 610             i == F_PIDIHS   ||
 611             i == F_IDIHS    ||
 612             i == F_LJ14     ||
 613             i == F_GB12     || /* The GB parameters are hardcoded both in */
 614             i == F_GB13     || /* Gromacs and OpenMM */
 615             i == F_GB14   ) &&
 616             top->idef.il[i].nr > 0)
 617         {
 618             gmx_fatal(FARGS, "OpenMM does not support (some) of the provided interaction "
 619                     "type(s) (%s) ", interaction_function[i].longname);
 620         }
 621     }
 622
 623     if (ir->efep != efepNO)
 624         gmx_fatal(FARGS,"OpenMM does not support free energy calculations.");
 625
 626     if (ir->opts.ngacc > 1)
 627         gmx_fatal(FARGS,"OpenMM does not support non-equilibrium MD (accelerated groups).");
 628
 629     if (IR_ELEC_FIELD(*ir))
 630         gmx_fatal(FARGS,"OpenMM does not support electric fields.");
 631
 632     if (ir->bQMMM)
 633         gmx_fatal(FARGS,"OpenMM does not support QMMM calculations.");
 634
 635     if (ir->rcoulomb != ir->rvdw)
 636         gmx_fatal(FARGS,"OpenMM uses a single cutoff for both Coulomb "
 637                   "and VdW interactions. Please set rcoulomb equal to rvdw.");
 638
 639     if (EEL_FULL(ir->coulombtype))
 640     {
 641         if (ir->ewald_geometry == eewg3DC)
 642             gmx_fatal(FARGS,"OpenMM supports only Ewald 3D geometry.");
 643         if (ir->epsilon_surface != 0)
 644             gmx_fatal(FARGS,"OpenMM does not support dipole correction in Ewald summation.");
 645     }
 646
 647     if (TRICLINIC(state->box))
 648     {
 649         gmx_fatal(FARGS,"OpenMM does not support triclinic unit cells.");
 650     }
 651
 652     /* XXX this is just debugging code to disable the combination rule check */
 653     if ( isStringEqNCase(opt->getOptionValue("check-combrule"), "yes") )
 654     {
 655     /* As OpenMM by default uses hardcoded combination rules
 656        sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j)
 657        we need to check whether the force field params obey this
 658        and if not, we can't use this force field so we exit
 659        grace-fatal-fully. */
 660     real *nbfp = fr->nbfp;
 661     natoms = fr->ntype;
 662     if (debug)
 663     {
 664         fprintf(debug, ">> Atom parameters: <<\n%10s%5s %5s %5s %5s COMB\n",
 665                 "", "i-j", "j-i", "i-i", "j-j");
 666     }
 667     /* loop over all i-j atom pairs and verify if
 668        sigma_ij = sigma_ji = sigma_comb and eps_ij = eps_ji = eps_comb */
 669     for (i = 0; i < natoms; i++)
 670     {
 671         /* i-i */
 672         c12 = C12(nbfp, natoms, i, i);
 673         c6  = C6(nbfp,  natoms, i, i);
 674         convert_c_12_6(c12, c6, &sigma_ii, &eps_ii);
 675
 676         for (j = 0; j < i; j++)
 677         {
 678             /* i-j */
 679             c12 = C12(nbfp, natoms, i, j);
 680             c6  = C6(nbfp,  natoms, i, j);
 681             convert_c_12_6(c12, c6, &sigma_ij, &eps_ij);
 682             /* j-i */
 683             c12 = C12(nbfp, natoms, j, i);
 684             c6  = C6(nbfp,  natoms, j, i);
 685             convert_c_12_6(c12, c6, &sigma_ji, &eps_ji);
 686             /* j-j */
 687             c12 = C12(nbfp, natoms, j, j);
 688             c6  = C6(nbfp,  natoms, j, j);
 689             convert_c_12_6(c12, c6, &sigma_jj, &eps_jj);
 690             /* OpenMM hardcoded combination rules */
 691             sigma_comb = COMBRULE_SIGMA(sigma_ii, sigma_jj);
 692             eps_comb = COMBRULE_EPS(eps_ii, eps_jj);
 693
 694             if (debug)
 695             {
 696                 fprintf(debug, "i=%-3d j=%-3d", i, j);
 697                 fprintf(debug, "%-11s", "sigma");
 698                 fprintf(debug, "%5.3f %5.3f %5.3f %5.3f %5.3f\n",
 699                         sigma_ij, sigma_ji, sigma_ii, sigma_jj, sigma_comb);
 700                 fprintf(debug, "%11s%-11s", "", "epsilon");
 701                 fprintf(debug, "%5.3f %5.3f %5.3f %5.3f %5.3f\n",
 702                         eps_ij, eps_ji, eps_ii, eps_jj, eps_comb);
 703             }
 704
 705             /* check the values against the rule used by omm */
 706             if((fabs(eps_ij) > COMBRULE_CHK_TOL &&
 707                 fabs(eps_ji) > COMBRULE_CHK_TOL) &&
 708                (fabs(sigma_comb - sigma_ij) > COMBRULE_CHK_TOL ||
 709                fabs(sigma_comb - sigma_ji) > COMBRULE_CHK_TOL ||
 710                fabs(eps_comb - eps_ij) > COMBRULE_CHK_TOL ||
 711                fabs(eps_comb - eps_ji) > COMBRULE_CHK_TOL ))
 712             {
 713                 gmx_fatal(FARGS,
 714                         "The combination rules of the used force-field do not "
 715                         "match the one supported by OpenMM:  "
 716                         "sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j). "
 717                         "Switch to a force-field that uses these rules in order to "
 718                         "simulate this system using OpenMM.\n");
 719             }
 720         }
 721     }
 722     if (debug) { fprintf(debug, ">><<\n\n"); }
 723
 724     /* if we got here, log that everything is fine */
 725     if (debug)
 726     {
 727         fprintf(debug, ">> The combination rule of the used force matches the one used by OpenMM.\n");
 728     }
 729     fprintf(fplog, "The combination rule of the used force field matches the one used by OpenMM.\n");
 730
 731     } /* if (are we checking the combination rules) ... */
 732 }
 733
 734
 735 /*!
 736  * \brief Initialize OpenMM, run sanity/consistency checks, and return a pointer to
 737  * the OpenMMData.
 738  *
 739  * Various gromacs data structures are passed that contain the parameters, state and
 740  * other porperties of the system to simulate. These serve as input for initializing
 741  * OpenMM. Besides, a set of misc action are taken:
 742  *  - OpenMM plugins are loaded;
 743  *  - platform options in \p platformOptStr are parsed and checked;
 744  *  - Gromacs parameters are checked for OpenMM support and consistency;
 745  *  - after the OpenMM is initialized memtest executed in the same GPU context.
 746  *
 747  * \param[in] fplog             Gromacs log file handler.
 748  * \param[in] platformOptStr    Platform option string.
 749  * \param[in] ir                The Gromacs input parameters, see ::t_inputrec
 750  * \param[in] top_global        Gromacs system toppology, \see ::gmx_mtop_t
 751  * \param[in] top               Gromacs node local topology, \see gmx_localtop_t
 752  * \param[in] mdatoms           Gromacs atom parameters, \see ::t_mdatoms
 753  * \param[in] fr                \see ::t_forcerec
 754  * \param[in] state             Gromacs systems state, \see ::t_state
 755  * \returns                     Pointer to a
 756  *
 757  */
 758 void* openmm_init(FILE *fplog, const char *platformOptStr,
 759                   t_inputrec *ir,
 760                   gmx_mtop_t *top_global, gmx_localtop_t *top,
 761                   t_mdatoms *mdatoms, t_forcerec *fr, t_state *state)
 762 {
 763
 764     char warn_buf[STRLEN];
 765     static gmx_bool hasLoadedPlugins = false;
 766     string usedPluginDir;
 767     int devId;
 768
 769     try
 770     {
 771         if (!hasLoadedPlugins)
 772         {
 773             vector<string> loadedPlugins;
 774             /*  Look for OpenMM plugins at various locations (listed in order of priority):
 775                 - on the path in OPENMM_PLUGIN_DIR environment variable if this is specified
 776                 - on the path in the OPENMM_PLUGIN_DIR macro that is set by the build script
 777                 - at the default location assumed by OpenMM
 778             */
 779             /* env var */
 780             char *pluginDir = getenv("OPENMM_PLUGIN_DIR");
 781             trim(pluginDir);
 782             /* no env var or empty */
 783             if (pluginDir != NULL && *pluginDir != '\0')
 784             {
 785                 loadedPlugins = Platform::loadPluginsFromDirectory(pluginDir);
 786                 if (!loadedPlugins.empty())
 787                 {
 788                     hasLoadedPlugins = true;
 789                     usedPluginDir = pluginDir;
 790                 }
 791                 else
 792                 {
 793                     gmx_fatal(FARGS, "The directory provided in the OPENMM_PLUGIN_DIR environment variable "
 794                               "(%s) does not contain valid OpenMM plugins. Check your OpenMM installation!",
 795                               pluginDir);
 796                 }
 797             }
 798
 799             /* macro set at build time  */
 800 #ifdef OPENMM_PLUGIN_DIR
 801             if (!hasLoadedPlugins)
 802             {
 803                 loadedPlugins = Platform::loadPluginsFromDirectory(OPENMM_PLUGIN_DIR);
 804                 if (!loadedPlugins.empty())
 805                 {
 806                     hasLoadedPlugins = true;
 807                     usedPluginDir = OPENMM_PLUGIN_DIR;
 808                 }
 809             }
 810 #endif
 811             /* default loocation */
 812             if (!hasLoadedPlugins)
 813             {
 814                 loadedPlugins = Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
 815                 if (!loadedPlugins.empty())
 816                 {
 817                     hasLoadedPlugins = true;
 818                     usedPluginDir = Platform::getDefaultPluginsDirectory();
 819                 }
 820             }
 821
 822             /* if there are still no plugins loaded there won't be any */
 823             if (!hasLoadedPlugins)
 824             {
 825                 gmx_fatal(FARGS, "No OpenMM plugins were found! You can provide the"
 826                           " plugin directory in the OPENMM_PLUGIN_DIR environment variable.", pluginDir);
 827             }
 828
 829             fprintf(fplog, "\nOpenMM plugins loaded from directory %s:\t", usedPluginDir.c_str());
 830             for (int i = 0; i < (int)loadedPlugins.size(); i++)
 831             {
 832                 fprintf(fplog, "%s, ", loadedPlugins[i].c_str());
 833             }
 834             fprintf(fplog, "\n");
 835         }
 836
 837         /* parse option string */
 838         GmxOpenMMPlatformOptions *opt = new GmxOpenMMPlatformOptions(platformOptStr);
 839         devId = atoi(opt->getOptionValue("deviceid").c_str());
 840
 841         if (debug)
 842         {
 843             opt->print();
 844         }
 845
 846         /* check wheter Gromacs options compatibility with OpenMM */
 847         checkGmxOptions(fplog, opt, ir, top, fr, state);
 848
 849         /* Create the system. */
 850         const t_idef& idef = top->idef;
 851         const int numAtoms = top_global->natoms;
 852         const int numConstraints = idef.il[F_CONSTR].nr/3;
 853         const int numSettle = idef.il[F_SETTLE].nr/2;
 854         const int numBonds = idef.il[F_BONDS].nr/3;
 855         const int numHarmonic = idef.il[F_HARMONIC].nr/3;
 856         const int numUB = idef.il[F_UREY_BRADLEY].nr/4;
 857         const int numAngles = idef.il[F_ANGLES].nr/4;
 858         const int numPeriodic = idef.il[F_PDIHS].nr/5;
 859         const int numPeriodicImproper = idef.il[F_PIDIHS].nr/5;
 860         const int numRB = idef.il[F_RBDIHS].nr/5;
 861         const int numImproperDih = idef.il[F_IDIHS].nr/5;
 862         const int num14 = idef.il[F_LJ14].nr/3;
 863         System* sys = new System();
 864         if (ir->nstcomm > 0)
 865             sys->addForce(new CMMotionRemover(ir->nstcomm));
 866
 867         /* Set bonded force field terms. */
 868
 869                 /*
 870                  * CUDA platform currently doesn't support more than one
 871                  * instance of a force object, so we pack all forces that
 872                  * use the same form into one.
 873                 */
 874
 875         const int* bondAtoms = (int*) idef.il[F_BONDS].iatoms;
 876         HarmonicBondForce* bondForce = new HarmonicBondForce();
 877         sys->addForce(bondForce);
 878         int offset = 0;
 879         for (int i = 0; i < numBonds; ++i)
 880         {
 881             int type = bondAtoms[offset++];
 882             int atom1 = bondAtoms[offset++];
 883             int atom2 = bondAtoms[offset++];
 884             bondForce->addBond(atom1, atom2,
 885                                idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
 886         }
 887
 888         const int* harmonicAtoms = (int*) idef.il[F_HARMONIC].iatoms;
 889         offset = 0;
 890         for (int i = 0; i < numHarmonic; ++i)
 891         {
 892             int type = harmonicAtoms[offset++];
 893             int atom1 = harmonicAtoms[offset++];
 894             int atom2 = harmonicAtoms[offset++];
 895             bondForce->addBond(atom1, atom2,
 896                                idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
 897         }
 898
 899                 /* Set the angle force field terms */
 900         const int* angleAtoms = (int*) idef.il[F_ANGLES].iatoms;
 901         HarmonicAngleForce* angleForce = new HarmonicAngleForce();
 902         sys->addForce(angleForce);
 903         offset = 0;
 904         for (int i = 0; i < numAngles; ++i)
 905         {
 906             int type = angleAtoms[offset++];
 907             int atom1 = angleAtoms[offset++];
 908             int atom2 = angleAtoms[offset++];
 909             int atom3 = angleAtoms[offset++];
 910             angleForce->addAngle(atom1, atom2, atom3,
 911                     idef.iparams[type].harmonic.rA*M_PI/180.0, idef.iparams[type].harmonic.krA);
 912         }
 913
 914         /* Urey-Bradley includes both the angle and bond potential for 1-3 interactions */
 915         const int* ubAtoms = (int*) idef.il[F_UREY_BRADLEY].iatoms;
 916                 /* HarmonicBondForce* ubBondForce = new HarmonicBondForce(); */
 917                 /*  HarmonicAngleForce* ubAngleForce = new HarmonicAngleForce(); */
 918         /* sys->addForce(ubBondForce); */
 919         /* sys->addForce(ubAngleForce); */
 920         offset = 0;
 921         for (int i = 0; i < numUB; ++i)
 922         {
 923             int type = ubAtoms[offset++];
 924             int atom1 = ubAtoms[offset++];
 925             int atom2 = ubAtoms[offset++];
 926             int atom3 = ubAtoms[offset++];
 927             /* ubBondForce->addBond(atom1, atom3, */
 928             bondForce->addBond(atom1, atom3,
 929                                idef.iparams[type].u_b.r13A, idef.iparams[type].u_b.kUBA);
 930             /* ubAngleForce->addAngle(atom1, atom2, atom3, */
 931             angleForce->addAngle(atom1, atom2, atom3,
 932                     idef.iparams[type].u_b.thetaA*M_PI/180.0, idef.iparams[type].u_b.kthetaA);
 933         }
 934
 935                 /* Set proper dihedral terms */
 936         const int* periodicAtoms = (int*) idef.il[F_PDIHS].iatoms;
 937         PeriodicTorsionForce* periodicForce = new PeriodicTorsionForce();
 938         sys->addForce(periodicForce);
 939         offset = 0;
 940         for (int i = 0; i < numPeriodic; ++i)
 941         {
 942             int type = periodicAtoms[offset++];
 943             int atom1 = periodicAtoms[offset++];
 944             int atom2 = periodicAtoms[offset++];
 945             int atom3 = periodicAtoms[offset++];
 946             int atom4 = periodicAtoms[offset++];
 947             periodicForce->addTorsion(atom1, atom2, atom3, atom4,
 948                                       idef.iparams[type].pdihs.mult,
 949                                       idef.iparams[type].pdihs.phiA*M_PI/180.0,
 950                                       idef.iparams[type].pdihs.cpA);
 951         }
 952
 953                 /* Set improper dihedral terms that are represented by a periodic function (as in AMBER FF) */
 954         const int* periodicImproperAtoms = (int*) idef.il[F_PIDIHS].iatoms;
 955         /* PeriodicTorsionForce* periodicImproperForce = new PeriodicTorsionForce(); */
 956         /* sys->addForce(periodicImproperForce); */
 957         offset = 0;
 958         for (int i = 0; i < numPeriodicImproper; ++i)
 959         {
 960             int type = periodicImproperAtoms[offset++];
 961             int atom1 = periodicImproperAtoms[offset++];
 962             int atom2 = periodicImproperAtoms[offset++];
 963             int atom3 = periodicImproperAtoms[offset++];
 964             int atom4 = periodicImproperAtoms[offset++];
 965             /* periodicImproperForce->addTorsion(atom1, atom2, atom3, atom4, */
 966             periodicForce->addTorsion(atom1, atom2, atom3, atom4,
 967                                       idef.iparams[type].pdihs.mult,
 968                                       idef.iparams[type].pdihs.phiA*M_PI/180.0,
 969                                       idef.iparams[type].pdihs.cpA);
 970         }
 971
 972         /* Ryckaert-Bellemans dihedrals */
 973         const int* rbAtoms = (int*) idef.il[F_RBDIHS].iatoms;
 974         RBTorsionForce* rbForce = new RBTorsionForce();
 975         sys->addForce(rbForce);
 976         offset = 0;
 977         for (int i = 0; i < numRB; ++i)
 978         {
 979             int type = rbAtoms[offset++];
 980             int atom1 = rbAtoms[offset++];
 981             int atom2 = rbAtoms[offset++];
 982             int atom3 = rbAtoms[offset++];
 983             int atom4 = rbAtoms[offset++];
 984             rbForce->addTorsion(atom1, atom2, atom3, atom4,
 985                                 idef.iparams[type].rbdihs.rbcA[0], idef.iparams[type].rbdihs.rbcA[1],
 986                                 idef.iparams[type].rbdihs.rbcA[2], idef.iparams[type].rbdihs.rbcA[3],
 987                                 idef.iparams[type].rbdihs.rbcA[4], idef.iparams[type].rbdihs.rbcA[5]);
 988         }
 989
 990                 /* Set improper dihedral terms (as in CHARMM FF) */
 991         const int* improperDihAtoms = (int*) idef.il[F_IDIHS].iatoms;
 992                 CustomTorsionForce* improperDihForce = new CustomTorsionForce("2.0*k*asin(sin((theta-theta0)/2))^2");
 993         sys->addForce(improperDihForce);
 994                 improperDihForce->addPerTorsionParameter("k");
 995                 improperDihForce->addPerTorsionParameter("theta0");
 996                 vector<double> improperDihParameters(2);
 997         offset = 0;
 998         for (int i = 0; i < numImproperDih; ++i)
 999         {
1000             int type = improperDihAtoms[offset++];
1001             int atom1 = improperDihAtoms[offset++];
1002             int atom2 = improperDihAtoms[offset++];
1003             int atom3 = improperDihAtoms[offset++];
1004             int atom4 = improperDihAtoms[offset++];
1005                         improperDihParameters[0] = idef.iparams[type].harmonic.krA;
1006                         improperDihParameters[1] = idef.iparams[type].harmonic.rA*M_PI/180.0;
1007             improperDihForce->addTorsion(atom1, atom2, atom3, atom4,
1008                                 improperDihParameters);
1009         }
1010
1011         /* Set nonbonded parameters and masses. */
1012         int ntypes = fr->ntype;
1013         int* types = mdatoms->typeA;
1014         real* nbfp = fr->nbfp;
1015         real* charges = mdatoms->chargeA;
1016         real* masses = mdatoms->massT;
1017         NonbondedForce* nonbondedForce = new NonbondedForce();
1018         sys->addForce(nonbondedForce);
1019
1020         switch (ir->ePBC)
1021         {
1022         case epbcNONE:
1023             if (ir->rcoulomb == 0)
1024             {
1025                 nonbondedForce->setNonbondedMethod(NonbondedForce::NoCutoff);
1026             }
1027             else
1028             {
1029                 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
1030             }
1031             break;
1032         case epbcXYZ:
1033             switch (ir->coulombtype)
1034             {
1035             case eelCUT:
1036             case eelRF:
1037             case eelGRF:
1038             case eelRF_NEC:
1039             case eelRF_ZERO:
1040                 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
1041                 break;
1042
1043             case eelEWALD:
1044                 nonbondedForce->setNonbondedMethod(NonbondedForce::Ewald);
1045                 break;
1046
1047             case eelPME:
1048                 nonbondedForce->setNonbondedMethod(NonbondedForce::PME);
1049                 break;
1050
1051             default:
1052                 gmx_fatal(FARGS,"Internal error: you should not see this message, it means that the"
1053                           "electrosatics option check failed. Please report this error!");
1054             }
1055             sys->setDefaultPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0),
1056                                        Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));
1057             nonbondedForce->setCutoffDistance(ir->rcoulomb);
1058
1059             break;
1060         default:
1061             gmx_fatal(FARGS,"OpenMM supports only full periodic boundary conditions "
1062                               "(pbc = xyz), or none (pbc = no).");
1063         }
1064
1065
1066         /* Fix for PME and Ewald error tolerance
1067          *
1068                  *  OpenMM uses approximate formulas to calculate the Ewald parameter:
1069                  *  alpha = (1.0/cutoff)*sqrt(-log(2.0*tolerlance));
1070                  *  and the grid spacing for PME:
1071                  *  gridX = ceil(2*alpha*box[0][0]/3*(pow(tol, 0.2)))
1072                  *  gridY = ceil(2*alpha*box[1][1]/3*(pow(tol, 0.2)));
1073                  *  gridZ = ceil(2*alpha*box[2][2]/3*(pow(tol, 0.2)));
1074                  *
1075                  *
1076                  *  If the default ewald_rtol=1e-5 is used we silently adjust the value to the
1077                  *  OpenMM default of 5e-4 otherwise a warning is issued about the action taken.
1078                  *
1079                 */
1080         double corr_ewald_rtol = 50.0 * ir->ewald_rtol;
1081         if ((ir->ePBC == epbcXYZ) &&
1082             (ir->coulombtype == eelEWALD || ir->coulombtype == eelPME))
1083         {
1084             if (debug)
1085             {
1086                 fprintf(debug, ">> ewald_rtol = %e (corrected = %e) \n",
1087                     ir->ewald_rtol, corr_ewald_rtol);
1088             }
1089
1090             if (fabs(ir->ewald_rtol - 1e-5) > 1e-10)
1091             {
1092                 gmx_warning("OpenMM uses the ewald_rtol parameter with approximate formulas "
1093                         "to calculate the alpha and grid spacing parameters of the Ewald "
1094                         "and PME methods. This tolerance need to be corrected in order to get "
1095                         "settings close to the ones used in GROMACS. Although the internal correction "
1096                         "should work for any reasonable value of ewald_rtol, using values other than "
1097                         "the default 1e-5 might cause incorrect behavior.");
1098
1099                 if (corr_ewald_rtol > 1)
1100                 {
1101                     gmx_fatal(FARGS, "The ewald_rtol accuracy term is >1 after the "
1102                             "adjustment for OpenMM (%e)", corr_ewald_rtol);
1103                 }
1104             }
1105             nonbondedForce->setEwaldErrorTolerance(corr_ewald_rtol);
1106         }
1107
1108         for (int i = 0; i < numAtoms; ++i)
1109         {
1110             double c12 = nbfp[types[i]*2*ntypes+types[i]*2+1];
1111             double c6 = nbfp[types[i]*2*ntypes+types[i]*2];
1112             double sigma=0.0, epsilon=0.0;
1113             convert_c_12_6(c12, c6, &sigma, &epsilon);
1114             nonbondedForce->addParticle(charges[i], sigma, epsilon);
1115             sys->addParticle(masses[i]);
1116         }
1117
1118         // Build a table of all exclusions.
1119         vector<set<int> > exclusions(numAtoms);
1120         for (int i = 0; i < numAtoms; i++)
1121         {
1122             int start = top->excls.index[i];
1123             int end = top->excls.index[i+1];
1124             for (int j = start; j < end; j++)
1125                 exclusions[i].insert(top->excls.a[j]);
1126         }
1127
1128         // Record the 1-4 interactions, and remove them from the list of exclusions.
1129         const int* nb14Atoms = (int*) idef.il[F_LJ14].iatoms;
1130         offset = 0;
1131         for (int i = 0; i < num14; ++i)
1132         {
1133             int type = nb14Atoms[offset++];
1134             int atom1 = nb14Atoms[offset++];
1135             int atom2 = nb14Atoms[offset++];
1136             double sigma=0, epsilon=0;
1137             convert_c_12_6(idef.iparams[type].lj14.c12A,
1138                     idef.iparams[type].lj14.c6A,
1139                     &sigma, &epsilon);
1140             nonbondedForce->addException(atom1, atom2,
1141                                          fr->fudgeQQ*charges[atom1]*charges[atom2], sigma, epsilon);
1142             exclusions[atom1].erase(atom2);
1143             exclusions[atom2].erase(atom1);
1144         }
1145
1146         // Record exclusions.
1147         for (int i = 0; i < numAtoms; i++)
1148         {
1149             for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
1150             {
1151                 if (i < *iter)
1152                 {
1153                     nonbondedForce->addException(i, *iter, 0.0, 1.0, 0.0);
1154                 }
1155             }
1156         }
1157
1158         // Add GBSA if needed.
1159         if (ir->implicit_solvent == eisGBSA)
1160         {
1161             gmx_warning("The OBC scale factors alpha, beta and gamma are hardcoded in OpenMM with the default Gromacs values.");
1162             t_atoms atoms       = gmx_mtop_global_atoms(top_global);
1163             GBSAOBCForce* gbsa  = new GBSAOBCForce();
1164
1165             sys->addForce(gbsa);
1166             gbsa->setSoluteDielectric(ir->epsilon_r);
1167             gbsa->setSolventDielectric(ir->gb_epsilon_solvent);
1168             gbsa->setCutoffDistance(nonbondedForce->getCutoffDistance());
1169             if (nonbondedForce->getNonbondedMethod() == NonbondedForce::NoCutoff)
1170                 gbsa->setNonbondedMethod(GBSAOBCForce::NoCutoff);
1171             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
1172                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
1173             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffPeriodic)
1174                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
1175             else
1176                 gmx_fatal(FARGS,"OpenMM supports only Reaction-Field electrostatics with OBC/GBSA.");
1177
1178             for (int i = 0; i < numAtoms; ++i)
1179             {
1180                 gbsa->addParticle(charges[i],
1181                                   top_global->atomtypes.gb_radius[atoms.atom[i].type],
1182                                   top_global->atomtypes.S_hct[atoms.atom[i].type]);
1183             }
1184             free_t_atoms(&atoms, FALSE);
1185         }
1186
1187         // Set constraints.
1188         const int* constraintAtoms = (int*) idef.il[F_CONSTR].iatoms;
1189         offset = 0;
1190         for (int i = 0; i < numConstraints; ++i)
1191         {
1192             int type = constraintAtoms[offset++];
1193             int atom1 = constraintAtoms[offset++];
1194             int atom2 = constraintAtoms[offset++];
1195             sys->addConstraint(atom1, atom2, idef.iparams[type].constr.dA);
1196         }
1197         const int* settleAtoms = (int*) idef.il[F_SETTLE].iatoms;
1198         offset = 0;
1199         for (int i = 0; i < numSettle; ++i)
1200         {
1201             int type = settleAtoms[offset++];
1202             int oxygen = settleAtoms[offset++];
1203             sys->addConstraint(oxygen, oxygen+1, idef.iparams[type].settle.doh);
1204             sys->addConstraint(oxygen, oxygen+2, idef.iparams[type].settle.doh);
1205             sys->addConstraint(oxygen+1, oxygen+2, idef.iparams[type].settle.dhh);
1206         }
1207
1208         // Create an integrator for simulating the system.
1209         double friction = (ir->opts.tau_t[0] == 0.0 ? 0.0 : 1.0/ir->opts.tau_t[0]);
1210         Integrator* integ;
1211         if (ir->eI == eiBD)
1212         {
1213             integ = new BrownianIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1214             static_cast<BrownianIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1215         }
1216         else if (EI_SD(ir->eI))
1217         {
1218             integ = new LangevinIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1219             static_cast<LangevinIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1220         }
1221         else
1222         {
1223             integ = new VerletIntegrator(ir->delta_t);
1224             if ( ir->etc != etcNO)
1225             {
1226                 AndersenThermostat* thermostat = new AndersenThermostat(ir->opts.ref_t[0], friction);
1227                 sys->addForce(thermostat);
1228             }
1229         }
1230
1231                 // Add pressure coupling
1232         if (ir->epc != epcNO)
1233                 {
1234           // convert gromacs pressure tensor to a scalar
1235           double pressure = (ir->ref_p[0][0] + ir->ref_p[1][1] + ir->ref_p[2][2]) / 3.0;
1236           int frequency = int(ir->tau_p / ir->delta_t); // update frequency in time steps
1237           if (frequency < 1) frequency = 1;
1238           double temperature = ir->opts.ref_t[0]; // in kelvin
1239           sys->addForce(new MonteCarloBarostat(pressure, temperature, frequency));
1240                 }
1241
1242         integ->setConstraintTolerance(ir->shake_tol);
1243
1244         // Create a context and initialize it.
1245         Context* context = NULL;
1246
1247         /*
1248         OpenMM could automatically select the "best" GPU, however we're not't
1249         going to let it do that for now, as the current algorithm is very rudimentary
1250         and we anyway support only CUDA.
1251         if (platformOptStr == NULL || platformOptStr == "")
1252         {
1253             context = new Context(*sys, *integ);
1254         }
1255         else
1256         */
1257         {
1258             /* which platform should we use */
1259             for (int i = 0; i < (int)Platform::getNumPlatforms() && context == NULL; i++)
1260             {
1261                 if (isStringEqNCase(opt->getOptionValue("platform"), Platform::getPlatform(i).getName()))
1262                 {
1263                     Platform& platform = Platform::getPlatform(i);
1264                     // set standard properties
1265                     platform.setPropertyDefaultValue("CudaDevice", opt->getOptionValue("deviceid"));
1266                     // TODO add extra properties
1267                     context = new Context(*sys, *integ, platform);
1268                 }
1269             }
1270             if (context == NULL)
1271             {
1272                 gmx_fatal(FARGS, "The requested platform \"%s\" could not be found.",
1273                         opt->getOptionValue("platform").c_str());
1274             }
1275         }
1276
1277         Platform& platform = context->getPlatform();
1278         fprintf(fplog, "Gromacs will use the OpenMM platform: %s\n", platform.getName().c_str());
1279
1280         const vector<string>& properties = platform.getPropertyNames();
1281         if (debug)
1282         {
1283             for (int i = 0; i < (int)properties.size(); i++)
1284             {
1285                 fprintf(debug, ">> %s: %s\n", properties[i].c_str(),
1286                         platform.getPropertyValue(*context, properties[i]).c_str());
1287             }
1288         }
1289
1290         /* only for CUDA */
1291         if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1292         {
1293             int tmp;
1294             if (!from_string<int>(tmp, platform.getPropertyValue(*context, "CudaDevice"), std::dec))
1295             {
1296                 gmx_fatal(FARGS, "Internal error: couldn't determine the device selected by OpenMM");
1297
1298             }
1299
1300             /* For now this is just to double-check if OpenMM selected the GPU we wanted,
1301             but when we'll let OpenMM select the GPU automatically, it will query the deviceId.
1302             */
1303             if (tmp != devId)
1304             {
1305                 gmx_fatal(FARGS, "Internal error: OpenMM is using device #%d"
1306                         "while initialized for device #%d", tmp, devId);
1307             }
1308
1309             /* check GPU compatibility */
1310             char gpuname[STRLEN];
1311             devId = atoi(opt->getOptionValue("deviceid").c_str());
1312             if (!is_gmx_openmm_supported_gpu(-1, gpuname))
1313             {
1314                 if (!gmx_strcasecmp(opt->getOptionValue("force-device").c_str(), "yes"))
1315                 {
1316                     sprintf(warn_buf, "Non-supported GPU selected (#%d, %s), forced continuing."
1317                             "Note, that the simulation can be slow or it migth even crash.",
1318                             devId, gpuname);
1319                     fprintf(fplog, "%s\n", warn_buf);
1320                     gmx_warning(warn_buf);
1321                 }
1322                 else
1323                 {
1324                     gmx_fatal(FARGS, "The selected GPU (#%d, %s) is not supported by Gromacs! "
1325                               "Most probably you have a low-end GPU which would not perform well, "
1326                               "or new hardware that has not been tested with the current release. "
1327                               "If you still want to try using the device, use the force-device=yes option.",
1328                               devId, gpuname);
1329                 }
1330             }
1331             else
1332             {
1333                 fprintf(fplog, "Gromacs will run on the GPU #%d (%s).\n", devId, gpuname);
1334             }
1335         }
1336
1337         /* only for CUDA */
1338         if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1339         {
1340             /* pre-simulation memtest */
1341             runMemtest(fplog, -1, "Pre", opt);
1342         }
1343
1344         vector<Vec3> pos(numAtoms);
1345         vector<Vec3> vel(numAtoms);
1346         for (int i = 0; i < numAtoms; ++i)
1347         {
1348             pos[i] = Vec3(state->x[i][0], state->x[i][1], state->x[i][2]);
1349             vel[i] = Vec3(state->v[i][0], state->v[i][1], state->v[i][2]);
1350         }
1351         context->setPositions(pos);
1352         context->setVelocities(vel);
1353
1354         // Return a structure containing the system, integrator, and context.
1355         OpenMMData* data = new OpenMMData();
1356         data->system = sys;
1357         data->integrator = integ;
1358         data->context = context;
1359         data->removeCM = (ir->nstcomm > 0);
1360         data->platformOpt = opt;
1361         return data;
1362     }
1363     catch (std::exception& e)
1364     {
1365         gmx_fatal(FARGS, "OpenMM exception caught while initializating: %s", e.what());
1366     }
1367     return NULL; /* just to avoid warnings */
1368 }
1369
1370 /*!
1371  * \brief Integrate one step.
1372  *
1373  * \param[in] data  OpenMMData object created by openmm_init().
1374  */
1375 void openmm_take_one_step(void* data)
1376 {
1377     // static int step = 0; printf("----> taking step #%d\n", step++);
1378     try
1379     {
1380         static_cast<OpenMMData*>(data)->integrator->step(1);
1381     }
1382     catch (std::exception& e)
1383     {
1384         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1385     }
1386 }
1387
1388 /*!
1389  * \brief Integrate n steps.
1390  *
1391  * \param[in] data  OpenMMData object created by openmm_init().
1392  */
1393 void openmm_take_steps(void* data, int nstep)
1394 {
1395     try
1396     {
1397         static_cast<OpenMMData*>(data)->integrator->step(nstep);
1398     }
1399     catch (std::exception& e)
1400     {
1401         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1402     }
1403 }
1404
1405 /*!
1406  * \brief Clean up the data structures cretead for OpenMM.
1407  *
1408  * \param[in] log   Log file pointer.
1409  * \param[in] data  OpenMMData object created by openmm_init().
1410  */
1411 void openmm_cleanup(FILE* fplog, void* data)
1412 {
1413     OpenMMData* d = static_cast<OpenMMData*>(data);
1414     /* only for CUDA */
1415     if (isStringEqNCase(d->platformOpt->getOptionValue("platform"), "CUDA"))
1416     {
1417         /* post-simulation memtest */
1418         runMemtest(fplog, -1, "Post", d->platformOpt);
1419     }
1420     delete d->system;
1421     delete d->integrator;
1422     delete d->context;
1423     delete d->platformOpt;
1424     delete d;
1425 }
1426
1427 /*!
1428  * \brief Copy the current state information from OpenMM into the Gromacs data structures.
1429  *
1430  * This function results in the requested proprties to be copied from the
1431  * GPU to host. As this represents a bottleneck, the frequency of pulling data
1432  * should be minimized.
1433  *
1434  * \param[in]   data        OpenMMData object created by openmm_init().
1435  * \param[out]  time        Simulation time for which the state was created.
1436  * \param[out]  state       State of the system: coordinates and velocities.
1437  * \param[out]  f           Forces.
1438  * \param[out]  enerd       Energies.
1439  * \param[in]   includePos  True if coordinates are requested.
1440  * \param[in]   includeVel  True if velocities are requested.
1441  * \param[in]   includeForce True if forces are requested.
1442  * \param[in]   includeEnergy True if energies are requested.
1443  */
1444 void openmm_copy_state(void *data,
1445                        t_state *state, double *time,
1446                        rvec f[], gmx_enerdata_t *enerd,
1447                        gmx_bool includePos, gmx_bool includeVel, gmx_bool includeForce, gmx_bool includeEnergy)
1448 {
1449     int types = 0;
1450     if (includePos)
1451         types += State::Positions;
1452     if (includeVel)
1453         types += State::Velocities;
1454     if (includeForce)
1455         types += State::Forces;
1456     if (includeEnergy)
1457         types += State::Energy;
1458     if (types == 0)
1459         return;
1460     try
1461     {
1462         State currentState = static_cast<OpenMMData*>(data)->context->getState(types);
1463         int numAtoms =  static_cast<OpenMMData*>(data)->system->getNumParticles();
1464         if (includePos)
1465         {
1466             for (int i = 0; i < numAtoms; i++)
1467             {
1468                 Vec3 x = currentState.getPositions()[i];
1469                 state->x[i][0] = x[0];
1470                 state->x[i][1] = x[1];
1471                 state->x[i][2] = x[2];
1472             }
1473         }
1474         if (includeVel)
1475         {
1476             for (int i = 0; i < numAtoms; i++)
1477             {
1478                 Vec3 v = currentState.getVelocities()[i];
1479                 state->v[i][0] = v[0];
1480                 state->v[i][1] = v[1];
1481                 state->v[i][2] = v[2];
1482             }
1483         }
1484         if (includeForce)
1485         {
1486             for (int i = 0; i < numAtoms; i++)
1487             {
1488                 Vec3 force = currentState.getForces()[i];
1489                 f[i][0] = force[0];
1490                 f[i][1] = force[1];
1491                 f[i][2] = force[2];
1492             }
1493         }
1494         if (includeEnergy)
1495         {
1496             int numConstraints = static_cast<OpenMMData*>(data)->system->getNumConstraints();
1497             int dof = 3*numAtoms-numConstraints;
1498             if (static_cast<OpenMMData*>(data)->removeCM)
1499                 dof -= 3;
1500             enerd->term[F_EPOT] = currentState.getPotentialEnergy();
1501             enerd->term[F_EKIN] = currentState.getKineticEnergy();
1502             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1503             enerd->term[F_TEMP] = 2.0*enerd->term[F_EKIN]/dof/BOLTZ;
1504         }
1505         *time = currentState.getTime();
1506     }
1507     catch (std::exception& e)
1508     {
1509         gmx_fatal(FARGS, "OpenMM exception caught while retrieving state information: %s", e.what());
1510     }
1511 }