src/programs/mdrun/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 center of mass 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     int     i, j, natoms;
 522     double  c6, c12;
 523     double  sigma_ij=0, sigma_ji=0, sigma_ii=0, sigma_jj=0, sigma_comb;
 524     double  eps_ij=0, eps_ji=0, eps_ii=0, eps_jj=0, eps_comb;
 525
 526     /* Abort if unsupported critical options are present */
 527
 528     /* Integrator */
 529     if (ir->eI ==  eiMD)
 530     {
 531         gmx_warning( "OpenMM does not support leap-frog, will use velocity-verlet integrator.");
 532     }
 533
 534     if (    (ir->eI !=  eiMD)   &&
 535             (ir->eI !=  eiVV)   &&
 536             (ir->eI !=  eiVVAK) &&
 537             (ir->eI !=  eiSD1)  &&
 538             (ir->eI !=  eiSD2)  &&
 539             (ir->eI !=  eiBD) )
 540     {
 541         gmx_fatal(FARGS, "OpenMM supports only the following integrators: md/md-vv/md-vv-avek, sd/sd1, and bd.");
 542     }
 543
 544     /* Electroctstics */
 545     if (   !(ir->coulombtype == eelPME   ||
 546              EEL_RF(ir->coulombtype)     ||
 547              ir->coulombtype == eelRF    ||
 548              ir->coulombtype == eelEWALD ||
 549              // no-cutoff
 550              (ir->coulombtype == eelCUT && ir->rcoulomb == 0 &&  ir->rvdw == 0) ||
 551              // we could have cut-off combined with GBSA (openmm will use RF)
 552              ir->implicit_solvent == eisGBSA)   )
 553     {
 554         gmx_fatal(FARGS,"OpenMM supports only the following methods for electrostatics: "
 555                 "NoCutoff (i.e. rcoulomb = rvdw = 0 ),Reaction-Field, Ewald or PME.");
 556     }
 557
 558     if (EEL_RF(ir->coulombtype) && ir->epsilon_rf != 0)
 559     {
 560         // openmm has epsilon_rf=inf hard-coded
 561         gmx_warning("OpenMM will use a Reaction-Field epsilon of infinity instead of %g.",ir->epsilon_rf);
 562     }
 563
 564     if (ir->etc != etcNO &&
 565         ir->eI  != eiSD1 &&
 566         ir->eI  != eiSD2 &&
 567         ir->eI  != eiBD )
 568     {
 569         gmx_warning("OpenMM supports only Andersen thermostat with the md/md-vv/md-vv-avek integrators.");
 570     }
 571
 572     if (ir->implicit_solvent == eisGBSA &&
 573         ir->gb_algorithm != egbOBC  )
 574     {
 575         gmx_warning("OpenMM does not support the specified algorithm for Generalized Born, will use OBC instead.");
 576     }
 577
 578     if (ir->opts.ngtc > 1)
 579         gmx_fatal(FARGS,"OpenMM does not support multiple temperature coupling groups.");
 580
 581     if (ir->epc != epcNO)
 582         gmx_warning("OpenMM supports only Monte Carlo barostat for pressure coupling.");
 583
 584     if (ir->opts.annealing[0])
 585         gmx_fatal(FARGS,"OpenMM does not support simulated annealing.");
 586
 587     if (top->idef.il[F_CONSTR].nr > 0 && ir->eConstrAlg != econtSHAKE)
 588         gmx_warning("OpenMM provides contraints as a combination "
 589                     "of SHAKE, SETTLE and CCMA. Accuracy is based on the SHAKE tolerance set "
 590                     "by the \"shake_tol\" option.");
 591
 592     if (ir->nwall != 0)
 593         gmx_fatal(FARGS,"OpenMM does not support walls.");
 594
 595     if (ir->ePull != epullNO)
 596         gmx_fatal(FARGS,"OpenMM does not support pulling.");
 597
 598     /* check for interaction types */
 599     for (i = 0; i < F_EPOT; i++)
 600     {
 601         if (!(i == F_CONSTR ||
 602             i == F_SETTLE   ||
 603             i == F_BONDS    ||
 604             i == F_HARMONIC ||
 605             i == F_UREY_BRADLEY ||
 606             i == F_ANGLES   ||
 607             i == F_PDIHS    ||
 608             i == F_RBDIHS   ||
 609             i == F_PIDIHS   ||
 610             i == F_IDIHS    ||
 611             i == F_LJ14     ||
 612             i == F_GB12     || /* The GB parameters are hardcoded both in */
 613             i == F_GB13     || /* Gromacs and OpenMM */
 614             i == F_GB14   ) &&
 615             top->idef.il[i].nr > 0)
 616         {
 617             gmx_fatal(FARGS, "OpenMM does not support (some) of the provided interaction "
 618                     "type(s) (%s) ", interaction_function[i].longname);
 619         }
 620     }
 621
 622     if (ir->efep != efepNO)
 623         gmx_fatal(FARGS,"OpenMM does not support free energy calculations.");
 624
 625     if (ir->opts.ngacc > 1)
 626         gmx_fatal(FARGS,"OpenMM does not support non-equilibrium MD (accelerated groups).");
 627
 628     if (IR_ELEC_FIELD(*ir))
 629         gmx_fatal(FARGS,"OpenMM does not support electric fields.");
 630
 631     if (ir->bQMMM)
 632         gmx_fatal(FARGS,"OpenMM does not support QMMM calculations.");
 633
 634     if (ir->rcoulomb != ir->rvdw)
 635         gmx_fatal(FARGS,"OpenMM uses a single cutoff for both Coulomb "
 636                   "and VdW interactions. Please set rcoulomb equal to rvdw.");
 637
 638     if (EEL_FULL(ir->coulombtype))
 639     {
 640         if (ir->ewald_geometry == eewg3DC)
 641             gmx_fatal(FARGS,"OpenMM supports only Ewald 3D geometry.");
 642         if (ir->epsilon_surface != 0)
 643             gmx_fatal(FARGS,"OpenMM does not support dipole correction in Ewald summation.");
 644     }
 645
 646     if (TRICLINIC(state->box))
 647     {
 648         gmx_fatal(FARGS,"OpenMM does not support triclinic unit cells.");
 649     }
 650
 651     /* XXX this is just debugging code to disable the combination rule check */
 652     if ( isStringEqNCase(opt->getOptionValue("check-combrule"), "yes") )
 653     {
 654     /* As OpenMM by default uses hardcoded combination rules
 655        sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j)
 656        we need to check whether the force field params obey this
 657        and if not, we can't use this force field so we exit
 658        grace-fatal-fully. */
 659     real *nbfp = fr->nbfp;
 660     natoms = fr->ntype;
 661     if (debug)
 662     {
 663         fprintf(debug, ">> Atom parameters: <<\n%10s%5s %5s %5s %5s COMB\n",
 664                 "", "i-j", "j-i", "i-i", "j-j");
 665     }
 666     /* loop over all i-j atom pairs and verify if
 667        sigma_ij = sigma_ji = sigma_comb and eps_ij = eps_ji = eps_comb */
 668     for (i = 0; i < natoms; i++)
 669     {
 670         /* i-i */
 671         /* nbfp now includes the 6.0/12.0 prefactors to save flops in kernels */
 672         c12 = C12(nbfp, natoms, i, i)/12.0;
 673         c6  = C6(nbfp,  natoms, i, i)/6.0;
 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)/12.0;
 680             c6  = C6(nbfp,  natoms, i, j)/6.0;
 681             convert_c_12_6(c12, c6, &sigma_ij, &eps_ij);
 682             /* j-i */
 683             c12 = C12(nbfp, natoms, j, i)/12.0;
 684             c6  = C6(nbfp,  natoms, j, i)/6.0;
 685             convert_c_12_6(c12, c6, &sigma_ji, &eps_ji);
 686             /* j-j */
 687             c12 = C12(nbfp, natoms, j, j)/12.0;
 688             c6  = C6(nbfp,  natoms, j, j)/6.0;
 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             /* nbfp now includes the 6.0/12.0 derivative prefactors to save flops in kernels*/
1111             double c12 = nbfp[types[i]*2*ntypes+types[i]*2+1]/12.0;
1112             double c6 = nbfp[types[i]*2*ntypes+types[i]*2]/6.0;
1113             double sigma=0.0, epsilon=0.0;
1114             convert_c_12_6(c12, c6, &sigma, &epsilon);
1115             nonbondedForce->addParticle(charges[i], sigma, epsilon);
1116             sys->addParticle(masses[i]);
1117         }
1118
1119         // Build a table of all exclusions.
1120         vector<set<int> > exclusions(numAtoms);
1121         for (int i = 0; i < numAtoms; i++)
1122         {
1123             int start = top->excls.index[i];
1124             int end = top->excls.index[i+1];
1125             for (int j = start; j < end; j++)
1126                 exclusions[i].insert(top->excls.a[j]);
1127         }
1128
1129         // Record the 1-4 interactions, and remove them from the list of exclusions.
1130         const int* nb14Atoms = (int*) idef.il[F_LJ14].iatoms;
1131         offset = 0;
1132         for (int i = 0; i < num14; ++i)
1133         {
1134             int type = nb14Atoms[offset++];
1135             int atom1 = nb14Atoms[offset++];
1136             int atom2 = nb14Atoms[offset++];
1137             double sigma=0, epsilon=0;
1138             convert_c_12_6(idef.iparams[type].lj14.c12A,
1139                     idef.iparams[type].lj14.c6A,
1140                     &sigma, &epsilon);
1141             nonbondedForce->addException(atom1, atom2,
1142                                          fr->fudgeQQ*charges[atom1]*charges[atom2], sigma, epsilon);
1143             exclusions[atom1].erase(atom2);
1144             exclusions[atom2].erase(atom1);
1145         }
1146
1147         // Record exclusions.
1148         for (int i = 0; i < numAtoms; i++)
1149         {
1150             for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
1151             {
1152                 if (i < *iter)
1153                 {
1154                     nonbondedForce->addException(i, *iter, 0.0, 1.0, 0.0);
1155                 }
1156             }
1157         }
1158
1159         // Add GBSA if needed.
1160         if (ir->implicit_solvent == eisGBSA)
1161         {
1162             gmx_warning("The OBC scale factors alpha, beta and gamma are hardcoded in OpenMM with the default Gromacs values.");
1163             t_atoms atoms       = gmx_mtop_global_atoms(top_global);
1164             GBSAOBCForce* gbsa  = new GBSAOBCForce();
1165
1166             sys->addForce(gbsa);
1167             gbsa->setSoluteDielectric(ir->epsilon_r);
1168             gbsa->setSolventDielectric(ir->gb_epsilon_solvent);
1169             gbsa->setCutoffDistance(nonbondedForce->getCutoffDistance());
1170             if (nonbondedForce->getNonbondedMethod() == NonbondedForce::NoCutoff)
1171                 gbsa->setNonbondedMethod(GBSAOBCForce::NoCutoff);
1172             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
1173                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
1174             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffPeriodic)
1175                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
1176             else
1177                 gmx_fatal(FARGS,"OpenMM supports only Reaction-Field electrostatics with OBC/GBSA.");
1178
1179             for (int i = 0; i < numAtoms; ++i)
1180             {
1181                 gbsa->addParticle(charges[i],
1182                                   top_global->atomtypes.gb_radius[atoms.atom[i].type],
1183                                   top_global->atomtypes.S_hct[atoms.atom[i].type]);
1184             }
1185             free_t_atoms(&atoms, FALSE);
1186         }
1187
1188         // Set constraints.
1189         const int* constraintAtoms = (int*) idef.il[F_CONSTR].iatoms;
1190         offset = 0;
1191         for (int i = 0; i < numConstraints; ++i)
1192         {
1193             int type = constraintAtoms[offset++];
1194             int atom1 = constraintAtoms[offset++];
1195             int atom2 = constraintAtoms[offset++];
1196             sys->addConstraint(atom1, atom2, idef.iparams[type].constr.dA);
1197         }
1198         const int* settleAtoms = (int*) idef.il[F_SETTLE].iatoms;
1199         offset = 0;
1200         for (int i = 0; i < numSettle; ++i)
1201         {
1202             int type = settleAtoms[offset++];
1203             int oxygen = settleAtoms[offset++];
1204             sys->addConstraint(oxygen, oxygen+1, idef.iparams[type].settle.doh);
1205             sys->addConstraint(oxygen, oxygen+2, idef.iparams[type].settle.doh);
1206             sys->addConstraint(oxygen+1, oxygen+2, idef.iparams[type].settle.dhh);
1207         }
1208
1209         // Create an integrator for simulating the system.
1210         double friction = (ir->opts.tau_t[0] == 0.0 ? 0.0 : 1.0/ir->opts.tau_t[0]);
1211         Integrator* integ;
1212         if (ir->eI == eiBD)
1213         {
1214             integ = new BrownianIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1215             static_cast<BrownianIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1216         }
1217         else if (EI_SD(ir->eI))
1218         {
1219             integ = new LangevinIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1220             static_cast<LangevinIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1221         }
1222         else
1223         {
1224             integ = new VerletIntegrator(ir->delta_t);
1225             if ( ir->etc != etcNO)
1226             {
1227                 AndersenThermostat* thermostat = new AndersenThermostat(ir->opts.ref_t[0], friction);
1228                 sys->addForce(thermostat);
1229             }
1230         }
1231
1232                 // Add pressure coupling
1233         if (ir->epc != epcNO)
1234                 {
1235           // convert gromacs pressure tensor to a scalar
1236           double pressure = (ir->ref_p[0][0] + ir->ref_p[1][1] + ir->ref_p[2][2]) / 3.0;
1237           int frequency = int(ir->tau_p / ir->delta_t); // update frequency in time steps
1238           if (frequency < 1) frequency = 1;
1239           double temperature = ir->opts.ref_t[0]; // in kelvin
1240           sys->addForce(new MonteCarloBarostat(pressure, temperature, frequency));
1241                 }
1242
1243         integ->setConstraintTolerance(ir->shake_tol);
1244
1245         // Create a context and initialize it.
1246         Context* context = NULL;
1247
1248         /*
1249         OpenMM could automatically select the "best" GPU, however we're not't
1250         going to let it do that for now, as the current algorithm is very rudimentary
1251         and we anyway support only CUDA.
1252         if (platformOptStr == NULL || platformOptStr == "")
1253         {
1254             context = new Context(*sys, *integ);
1255         }
1256         else
1257         */
1258         {
1259             /* which platform should we use */
1260             for (int i = 0; i < (int)Platform::getNumPlatforms() && context == NULL; i++)
1261             {
1262                 if (isStringEqNCase(opt->getOptionValue("platform"), Platform::getPlatform(i).getName()))
1263                 {
1264                     Platform& platform = Platform::getPlatform(i);
1265                     // set standard properties
1266                     platform.setPropertyDefaultValue("CudaDevice", opt->getOptionValue("deviceid"));
1267                     // TODO add extra properties
1268                     context = new Context(*sys, *integ, platform);
1269                 }
1270             }
1271             if (context == NULL)
1272             {
1273                 gmx_fatal(FARGS, "The requested platform \"%s\" could not be found.",
1274                         opt->getOptionValue("platform").c_str());
1275             }
1276         }
1277
1278         Platform& platform = context->getPlatform();
1279         fprintf(fplog, "Gromacs will use the OpenMM platform: %s\n", platform.getName().c_str());
1280
1281         const vector<string>& properties = platform.getPropertyNames();
1282         if (debug)
1283         {
1284             for (int i = 0; i < (int)properties.size(); i++)
1285             {
1286                 fprintf(debug, ">> %s: %s\n", properties[i].c_str(),
1287                         platform.getPropertyValue(*context, properties[i]).c_str());
1288             }
1289         }
1290
1291         /* only for CUDA */
1292         if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1293         {
1294             int tmp;
1295             if (!from_string<int>(tmp, platform.getPropertyValue(*context, "CudaDevice"), std::dec))
1296             {
1297                 gmx_fatal(FARGS, "Internal error: couldn't determine the device selected by OpenMM");
1298
1299             }
1300
1301             /* For now this is just to double-check if OpenMM selected the GPU we wanted,
1302             but when we'll let OpenMM select the GPU automatically, it will query the deviceId.
1303             */
1304             if (tmp != devId)
1305             {
1306                 gmx_fatal(FARGS, "Internal error: OpenMM is using device #%d"
1307                         "while initialized for device #%d", tmp, devId);
1308             }
1309
1310             /* check GPU compatibility */
1311             char gpuname[STRLEN];
1312             devId = atoi(opt->getOptionValue("deviceid").c_str());
1313             if (!is_gmx_openmm_supported_gpu(-1, gpuname))
1314             {
1315                 if (!gmx_strcasecmp(opt->getOptionValue("force-device").c_str(), "yes"))
1316                 {
1317                     sprintf(warn_buf, "Non-supported GPU selected (#%d, %s), forced continuing."
1318                             "Note, that the simulation can be slow or it migth even crash.",
1319                             devId, gpuname);
1320                     fprintf(fplog, "%s\n", warn_buf);
1321                     gmx_warning(warn_buf);
1322                 }
1323                 else
1324                 {
1325                     gmx_fatal(FARGS, "The selected GPU (#%d, %s) is not supported by Gromacs! "
1326                               "Most probably you have a low-end GPU which would not perform well, "
1327                               "or new hardware that has not been tested with the current release. "
1328                               "If you still want to try using the device, use the force-device=yes option.",
1329                               devId, gpuname);
1330                 }
1331             }
1332             else
1333             {
1334                 fprintf(fplog, "Gromacs will run on the GPU #%d (%s).\n", devId, gpuname);
1335             }
1336         }
1337
1338         /* only for CUDA */
1339         if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1340         {
1341             /* pre-simulation memtest */
1342             runMemtest(fplog, -1, "Pre", opt);
1343         }
1344
1345         vector<Vec3> pos(numAtoms);
1346         vector<Vec3> vel(numAtoms);
1347         for (int i = 0; i < numAtoms; ++i)
1348         {
1349             pos[i] = Vec3(state->x[i][0], state->x[i][1], state->x[i][2]);
1350             vel[i] = Vec3(state->v[i][0], state->v[i][1], state->v[i][2]);
1351         }
1352         context->setPositions(pos);
1353         context->setVelocities(vel);
1354
1355         // Return a structure containing the system, integrator, and context.
1356         OpenMMData* data = new OpenMMData();
1357         data->system = sys;
1358         data->integrator = integ;
1359         data->context = context;
1360         data->removeCM = (ir->nstcomm > 0);
1361         data->platformOpt = opt;
1362         return data;
1363     }
1364     catch (std::exception& e)
1365     {
1366         gmx_fatal(FARGS, "OpenMM exception caught while initializating: %s", e.what());
1367     }
1368     return NULL; /* just to avoid warnings */
1369 }
1370
1371 /*!
1372  * \brief Integrate one step.
1373  *
1374  * \param[in] data  OpenMMData object created by openmm_init().
1375  */
1376 void openmm_take_one_step(void* data)
1377 {
1378     // static int step = 0; printf("----> taking step #%d\n", step++);
1379     try
1380     {
1381         static_cast<OpenMMData*>(data)->integrator->step(1);
1382     }
1383     catch (std::exception& e)
1384     {
1385         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1386     }
1387 }
1388
1389 /*!
1390  * \brief Integrate n steps.
1391  *
1392  * \param[in] data  OpenMMData object created by openmm_init().
1393  */
1394 void openmm_take_steps(void* data, int nstep)
1395 {
1396     try
1397     {
1398         static_cast<OpenMMData*>(data)->integrator->step(nstep);
1399     }
1400     catch (std::exception& e)
1401     {
1402         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1403     }
1404 }
1405
1406 /*!
1407  * \brief Clean up the data structures cretead for OpenMM.
1408  *
1409  * \param[in] log   Log file pointer.
1410  * \param[in] data  OpenMMData object created by openmm_init().
1411  */
1412 void openmm_cleanup(FILE* fplog, void* data)
1413 {
1414     OpenMMData* d = static_cast<OpenMMData*>(data);
1415     /* only for CUDA */
1416     if (isStringEqNCase(d->platformOpt->getOptionValue("platform"), "CUDA"))
1417     {
1418         /* post-simulation memtest */
1419         runMemtest(fplog, -1, "Post", d->platformOpt);
1420     }
1421     delete d->system;
1422     delete d->integrator;
1423     delete d->context;
1424     delete d->platformOpt;
1425     delete d;
1426 }
1427
1428 /*!
1429  * \brief Copy the current state information from OpenMM into the Gromacs data structures.
1430  *
1431  * This function results in the requested proprties to be copied from the
1432  * GPU to host. As this represents a bottleneck, the frequency of pulling data
1433  * should be minimized.
1434  *
1435  * \param[in]   data        OpenMMData object created by openmm_init().
1436  * \param[out]  time        Simulation time for which the state was created.
1437  * \param[out]  state       State of the system: coordinates and velocities.
1438  * \param[out]  f           Forces.
1439  * \param[out]  enerd       Energies.
1440  * \param[in]   includePos  True if coordinates are requested.
1441  * \param[in]   includeVel  True if velocities are requested.
1442  * \param[in]   includeForce True if forces are requested.
1443  * \param[in]   includeEnergy True if energies are requested.
1444  */
1445 void openmm_copy_state(void *data,
1446                        t_state *state, double *time,
1447                        rvec f[], gmx_enerdata_t *enerd,
1448                        gmx_bool includePos, gmx_bool includeVel, gmx_bool includeForce, gmx_bool includeEnergy)
1449 {
1450     int types = 0;
1451     if (includePos)
1452         types += State::Positions;
1453     if (includeVel)
1454         types += State::Velocities;
1455     if (includeForce)
1456         types += State::Forces;
1457     if (includeEnergy)
1458         types += State::Energy;
1459     if (types == 0)
1460         return;
1461     try
1462     {
1463         State currentState = static_cast<OpenMMData*>(data)->context->getState(types);
1464         int numAtoms =  static_cast<OpenMMData*>(data)->system->getNumParticles();
1465         if (includePos)
1466         {
1467             for (int i = 0; i < numAtoms; i++)
1468             {
1469                 Vec3 x = currentState.getPositions()[i];
1470                 state->x[i][0] = x[0];
1471                 state->x[i][1] = x[1];
1472                 state->x[i][2] = x[2];
1473             }
1474         }
1475         if (includeVel)
1476         {
1477             for (int i = 0; i < numAtoms; i++)
1478             {
1479                 Vec3 v = currentState.getVelocities()[i];
1480                 state->v[i][0] = v[0];
1481                 state->v[i][1] = v[1];
1482                 state->v[i][2] = v[2];
1483             }
1484         }
1485         if (includeForce)
1486         {
1487             for (int i = 0; i < numAtoms; i++)
1488             {
1489                 Vec3 force = currentState.getForces()[i];
1490                 f[i][0] = force[0];
1491                 f[i][1] = force[1];
1492                 f[i][2] = force[2];
1493             }
1494         }
1495         if (includeEnergy)
1496         {
1497             int numConstraints = static_cast<OpenMMData*>(data)->system->getNumConstraints();
1498             int dof = 3*numAtoms-numConstraints;
1499             if (static_cast<OpenMMData*>(data)->removeCM)
1500                 dof -= 3;
1501             enerd->term[F_EPOT] = currentState.getPotentialEnergy();
1502             enerd->term[F_EKIN] = currentState.getKineticEnergy();
1503             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1504             enerd->term[F_TEMP] = 2.0*enerd->term[F_EKIN]/dof/BOLTZ;
1505         }
1506         *time = currentState.getTime();
1507     }
1508     catch (std::exception& e)
1509     {
1510         gmx_fatal(FARGS, "OpenMM exception caught while retrieving state information: %s", e.what());
1511     }
1512 }