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 "gmx_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         c12 = C12(nbfp, natoms, i, i);
 672         c6  = C6(nbfp,  natoms, i, i);
 673         convert_c_12_6(c12, c6, &sigma_ii, &eps_ii);
 674
 675         for (j = 0; j < i; j++)
 676         {
 677             /* i-j */
 678             c12 = C12(nbfp, natoms, i, j);
 679             c6  = C6(nbfp,  natoms, i, j);
 680             convert_c_12_6(c12, c6, &sigma_ij, &eps_ij);
 681             /* j-i */
 682             c12 = C12(nbfp, natoms, j, i);
 683             c6  = C6(nbfp,  natoms, j, i);
 684             convert_c_12_6(c12, c6, &sigma_ji, &eps_ji);
 685             /* j-j */
 686             c12 = C12(nbfp, natoms, j, j);
 687             c6  = C6(nbfp,  natoms, j, j);
 688             convert_c_12_6(c12, c6, &sigma_jj, &eps_jj);
 689             /* OpenMM hardcoded combination rules */
 690             sigma_comb = COMBRULE_SIGMA(sigma_ii, sigma_jj);
 691             eps_comb = COMBRULE_EPS(eps_ii, eps_jj);
 692
 693             if (debug)
 694             {
 695                 fprintf(debug, "i=%-3d j=%-3d", i, j);
 696                 fprintf(debug, "%-11s", "sigma");
 697                 fprintf(debug, "%5.3f %5.3f %5.3f %5.3f %5.3f\n",
 698                         sigma_ij, sigma_ji, sigma_ii, sigma_jj, sigma_comb);
 699                 fprintf(debug, "%11s%-11s", "", "epsilon");
 700                 fprintf(debug, "%5.3f %5.3f %5.3f %5.3f %5.3f\n",
 701                         eps_ij, eps_ji, eps_ii, eps_jj, eps_comb);
 702             }
 703
 704             /* check the values against the rule used by omm */
 705             if((fabs(eps_ij) > COMBRULE_CHK_TOL &&
 706                 fabs(eps_ji) > COMBRULE_CHK_TOL) &&
 707                (fabs(sigma_comb - sigma_ij) > COMBRULE_CHK_TOL ||
 708                fabs(sigma_comb - sigma_ji) > COMBRULE_CHK_TOL ||
 709                fabs(eps_comb - eps_ij) > COMBRULE_CHK_TOL ||
 710                fabs(eps_comb - eps_ji) > COMBRULE_CHK_TOL ))
 711             {
 712                 gmx_fatal(FARGS,
 713                         "The combination rules of the used force-field do not "
 714                         "match the one supported by OpenMM:  "
 715                         "sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j). "
 716                         "Switch to a force-field that uses these rules in order to "
 717                         "simulate this system using OpenMM.\n");
 718             }
 719         }
 720     }
 721     if (debug) { fprintf(debug, ">><<\n\n"); }
 722
 723     /* if we got here, log that everything is fine */
 724     if (debug)
 725     {
 726         fprintf(debug, ">> The combination rule of the used force matches the one used by OpenMM.\n");
 727     }
 728     fprintf(fplog, "The combination rule of the used force field matches the one used by OpenMM.\n");
 729
 730     } /* if (are we checking the combination rules) ... */
 731 }
 732
 733
 734 /*!
 735  * \brief Initialize OpenMM, run sanity/consistency checks, and return a pointer to
 736  * the OpenMMData.
 737  *
 738  * Various gromacs data structures are passed that contain the parameters, state and
 739  * other porperties of the system to simulate. These serve as input for initializing
 740  * OpenMM. Besides, a set of misc action are taken:
 741  *  - OpenMM plugins are loaded;
 742  *  - platform options in \p platformOptStr are parsed and checked;
 743  *  - Gromacs parameters are checked for OpenMM support and consistency;
 744  *  - after the OpenMM is initialized memtest executed in the same GPU context.
 745  *
 746  * \param[in] fplog             Gromacs log file handler.
 747  * \param[in] platformOptStr    Platform option string.
 748  * \param[in] ir                The Gromacs input parameters, see ::t_inputrec
 749  * \param[in] top_global        Gromacs system toppology, \see ::gmx_mtop_t
 750  * \param[in] top               Gromacs node local topology, \see gmx_localtop_t
 751  * \param[in] mdatoms           Gromacs atom parameters, \see ::t_mdatoms
 752  * \param[in] fr                \see ::t_forcerec
 753  * \param[in] state             Gromacs systems state, \see ::t_state
 754  * \returns                     Pointer to a
 755  *
 756  */
 757 void* openmm_init(FILE *fplog, const char *platformOptStr,
 758                   t_inputrec *ir,
 759                   gmx_mtop_t *top_global, gmx_localtop_t *top,
 760                   t_mdatoms *mdatoms, t_forcerec *fr, t_state *state)
 761 {
 762
 763     char warn_buf[STRLEN];
 764     static gmx_bool hasLoadedPlugins = false;
 765     string usedPluginDir;
 766     int devId;
 767
 768     try
 769     {
 770         if (!hasLoadedPlugins)
 771         {
 772             vector<string> loadedPlugins;
 773             /*  Look for OpenMM plugins at various locations (listed in order of priority):
 774                 - on the path in OPENMM_PLUGIN_DIR environment variable if this is specified
 775                 - on the path in the OPENMM_PLUGIN_DIR macro that is set by the build script
 776                 - at the default location assumed by OpenMM
 777             */
 778             /* env var */
 779             char *pluginDir = getenv("OPENMM_PLUGIN_DIR");
 780             trim(pluginDir);
 781             /* no env var or empty */
 782             if (pluginDir != NULL && *pluginDir != '\0')
 783             {
 784                 loadedPlugins = Platform::loadPluginsFromDirectory(pluginDir);
 785                 if (!loadedPlugins.empty())
 786                 {
 787                     hasLoadedPlugins = true;
 788                     usedPluginDir = pluginDir;
 789                 }
 790                 else
 791                 {
 792                     gmx_fatal(FARGS, "The directory provided in the OPENMM_PLUGIN_DIR environment variable "
 793                               "(%s) does not contain valid OpenMM plugins. Check your OpenMM installation!",
 794                               pluginDir);
 795                 }
 796             }
 797
 798             /* macro set at build time  */
 799 #ifdef OPENMM_PLUGIN_DIR
 800             if (!hasLoadedPlugins)
 801             {
 802                 loadedPlugins = Platform::loadPluginsFromDirectory(OPENMM_PLUGIN_DIR);
 803                 if (!loadedPlugins.empty())
 804                 {
 805                     hasLoadedPlugins = true;
 806                     usedPluginDir = OPENMM_PLUGIN_DIR;
 807                 }
 808             }
 809 #endif
 810             /* default loocation */
 811             if (!hasLoadedPlugins)
 812             {
 813                 loadedPlugins = Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
 814                 if (!loadedPlugins.empty())
 815                 {
 816                     hasLoadedPlugins = true;
 817                     usedPluginDir = Platform::getDefaultPluginsDirectory();
 818                 }
 819             }
 820
 821             /* if there are still no plugins loaded there won't be any */
 822             if (!hasLoadedPlugins)
 823             {
 824                 gmx_fatal(FARGS, "No OpenMM plugins were found! You can provide the"
 825                           " plugin directory in the OPENMM_PLUGIN_DIR environment variable.", pluginDir);
 826             }
 827
 828             fprintf(fplog, "\nOpenMM plugins loaded from directory %s:\t", usedPluginDir.c_str());
 829             for (int i = 0; i < (int)loadedPlugins.size(); i++)
 830             {
 831                 fprintf(fplog, "%s, ", loadedPlugins[i].c_str());
 832             }
 833             fprintf(fplog, "\n");
 834         }
 835
 836         /* parse option string */
 837         GmxOpenMMPlatformOptions *opt = new GmxOpenMMPlatformOptions(platformOptStr);
 838         devId = atoi(opt->getOptionValue("deviceid").c_str());
 839
 840         if (debug)
 841         {
 842             opt->print();
 843         }
 844
 845         /* check wheter Gromacs options compatibility with OpenMM */
 846         checkGmxOptions(fplog, opt, ir, top, fr, state);
 847
 848         /* Create the system. */
 849         const t_idef& idef = top->idef;
 850         const int numAtoms = top_global->natoms;
 851         const int numConstraints = idef.il[F_CONSTR].nr/3;
 852         const int numSettle = idef.il[F_SETTLE].nr/2;
 853         const int numBonds = idef.il[F_BONDS].nr/3;
 854         const int numHarmonic = idef.il[F_HARMONIC].nr/3;
 855         const int numUB = idef.il[F_UREY_BRADLEY].nr/4;
 856         const int numAngles = idef.il[F_ANGLES].nr/4;
 857         const int numPeriodic = idef.il[F_PDIHS].nr/5;
 858         const int numPeriodicImproper = idef.il[F_PIDIHS].nr/5;
 859         const int numRB = idef.il[F_RBDIHS].nr/5;
 860         const int numImproperDih = idef.il[F_IDIHS].nr/5;
 861         const int num14 = idef.il[F_LJ14].nr/3;
 862         System* sys = new System();
 863         if (ir->nstcomm > 0)
 864             sys->addForce(new CMMotionRemover(ir->nstcomm));
 865
 866         /* Set bonded force field terms. */
 867
 868                 /*
 869                  * CUDA platform currently doesn't support more than one
 870                  * instance of a force object, so we pack all forces that
 871                  * use the same form into one.
 872                 */
 873
 874         const int* bondAtoms = (int*) idef.il[F_BONDS].iatoms;
 875         HarmonicBondForce* bondForce = new HarmonicBondForce();
 876         sys->addForce(bondForce);
 877         int offset = 0;
 878         for (int i = 0; i < numBonds; ++i)
 879         {
 880             int type = bondAtoms[offset++];
 881             int atom1 = bondAtoms[offset++];
 882             int atom2 = bondAtoms[offset++];
 883             bondForce->addBond(atom1, atom2,
 884                                idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
 885         }
 886
 887         const int* harmonicAtoms = (int*) idef.il[F_HARMONIC].iatoms;
 888         offset = 0;
 889         for (int i = 0; i < numHarmonic; ++i)
 890         {
 891             int type = harmonicAtoms[offset++];
 892             int atom1 = harmonicAtoms[offset++];
 893             int atom2 = harmonicAtoms[offset++];
 894             bondForce->addBond(atom1, atom2,
 895                                idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
 896         }
 897
 898                 /* Set the angle force field terms */
 899         const int* angleAtoms = (int*) idef.il[F_ANGLES].iatoms;
 900         HarmonicAngleForce* angleForce = new HarmonicAngleForce();
 901         sys->addForce(angleForce);
 902         offset = 0;
 903         for (int i = 0; i < numAngles; ++i)
 904         {
 905             int type = angleAtoms[offset++];
 906             int atom1 = angleAtoms[offset++];
 907             int atom2 = angleAtoms[offset++];
 908             int atom3 = angleAtoms[offset++];
 909             angleForce->addAngle(atom1, atom2, atom3,
 910                     idef.iparams[type].harmonic.rA*M_PI/180.0, idef.iparams[type].harmonic.krA);
 911         }
 912
 913         /* Urey-Bradley includes both the angle and bond potential for 1-3 interactions */
 914         const int* ubAtoms = (int*) idef.il[F_UREY_BRADLEY].iatoms;
 915                 /* HarmonicBondForce* ubBondForce = new HarmonicBondForce(); */
 916                 /*  HarmonicAngleForce* ubAngleForce = new HarmonicAngleForce(); */
 917         /* sys->addForce(ubBondForce); */
 918         /* sys->addForce(ubAngleForce); */
 919         offset = 0;
 920         for (int i = 0; i < numUB; ++i)
 921         {
 922             int type = ubAtoms[offset++];
 923             int atom1 = ubAtoms[offset++];
 924             int atom2 = ubAtoms[offset++];
 925             int atom3 = ubAtoms[offset++];
 926             /* ubBondForce->addBond(atom1, atom3, */
 927             bondForce->addBond(atom1, atom3,
 928                                idef.iparams[type].u_b.r13, idef.iparams[type].u_b.kUB);
 929             /* ubAngleForce->addAngle(atom1, atom2, atom3, */
 930             angleForce->addAngle(atom1, atom2, atom3,
 931                     idef.iparams[type].u_b.theta*M_PI/180.0, idef.iparams[type].u_b.ktheta);
 932         }
 933
 934                 /* Set proper dihedral terms */
 935         const int* periodicAtoms = (int*) idef.il[F_PDIHS].iatoms;
 936         PeriodicTorsionForce* periodicForce = new PeriodicTorsionForce();
 937         sys->addForce(periodicForce);
 938         offset = 0;
 939         for (int i = 0; i < numPeriodic; ++i)
 940         {
 941             int type = periodicAtoms[offset++];
 942             int atom1 = periodicAtoms[offset++];
 943             int atom2 = periodicAtoms[offset++];
 944             int atom3 = periodicAtoms[offset++];
 945             int atom4 = periodicAtoms[offset++];
 946             periodicForce->addTorsion(atom1, atom2, atom3, atom4,
 947                                       idef.iparams[type].pdihs.mult,
 948                                       idef.iparams[type].pdihs.phiA*M_PI/180.0,
 949                                       idef.iparams[type].pdihs.cpA);
 950         }
 951
 952                 /* Set improper dihedral terms that are represented by a periodic function (as in AMBER FF) */
 953         const int* periodicImproperAtoms = (int*) idef.il[F_PIDIHS].iatoms;
 954         /* PeriodicTorsionForce* periodicImproperForce = new PeriodicTorsionForce(); */
 955         /* sys->addForce(periodicImproperForce); */
 956         offset = 0;
 957         for (int i = 0; i < numPeriodicImproper; ++i)
 958         {
 959             int type = periodicImproperAtoms[offset++];
 960             int atom1 = periodicImproperAtoms[offset++];
 961             int atom2 = periodicImproperAtoms[offset++];
 962             int atom3 = periodicImproperAtoms[offset++];
 963             int atom4 = periodicImproperAtoms[offset++];
 964             /* periodicImproperForce->addTorsion(atom1, atom2, atom3, atom4, */
 965             periodicForce->addTorsion(atom1, atom2, atom3, atom4,
 966                                       idef.iparams[type].pdihs.mult,
 967                                       idef.iparams[type].pdihs.phiA*M_PI/180.0,
 968                                       idef.iparams[type].pdihs.cpA);
 969         }
 970
 971         /* Ryckaert-Bellemans dihedrals */
 972         const int* rbAtoms = (int*) idef.il[F_RBDIHS].iatoms;
 973         RBTorsionForce* rbForce = new RBTorsionForce();
 974         sys->addForce(rbForce);
 975         offset = 0;
 976         for (int i = 0; i < numRB; ++i)
 977         {
 978             int type = rbAtoms[offset++];
 979             int atom1 = rbAtoms[offset++];
 980             int atom2 = rbAtoms[offset++];
 981             int atom3 = rbAtoms[offset++];
 982             int atom4 = rbAtoms[offset++];
 983             rbForce->addTorsion(atom1, atom2, atom3, atom4,
 984                                 idef.iparams[type].rbdihs.rbcA[0], idef.iparams[type].rbdihs.rbcA[1],
 985                                 idef.iparams[type].rbdihs.rbcA[2], idef.iparams[type].rbdihs.rbcA[3],
 986                                 idef.iparams[type].rbdihs.rbcA[4], idef.iparams[type].rbdihs.rbcA[5]);
 987         }
 988
 989                 /* Set improper dihedral terms (as in CHARMM FF) */
 990         const int* improperDihAtoms = (int*) idef.il[F_IDIHS].iatoms;
 991                 CustomTorsionForce* improperDihForce = new CustomTorsionForce("2.0*k*asin(sin((theta-theta0)/2))^2");
 992         sys->addForce(improperDihForce);
 993                 improperDihForce->addPerTorsionParameter("k");
 994                 improperDihForce->addPerTorsionParameter("theta0");
 995                 vector<double> improperDihParameters(2);
 996         offset = 0;
 997         for (int i = 0; i < numImproperDih; ++i)
 998         {
 999             int type = improperDihAtoms[offset++];
1000             int atom1 = improperDihAtoms[offset++];
1001             int atom2 = improperDihAtoms[offset++];
1002             int atom3 = improperDihAtoms[offset++];
1003             int atom4 = improperDihAtoms[offset++];
1004                         improperDihParameters[0] = idef.iparams[type].harmonic.krA;
1005                         improperDihParameters[1] = idef.iparams[type].harmonic.rA*M_PI/180.0;
1006             improperDihForce->addTorsion(atom1, atom2, atom3, atom4,
1007                                 improperDihParameters);
1008         }
1009
1010         /* Set nonbonded parameters and masses. */
1011         int ntypes = fr->ntype;
1012         int* types = mdatoms->typeA;
1013         real* nbfp = fr->nbfp;
1014         real* charges = mdatoms->chargeA;
1015         real* masses = mdatoms->massT;
1016         NonbondedForce* nonbondedForce = new NonbondedForce();
1017         sys->addForce(nonbondedForce);
1018
1019         switch (ir->ePBC)
1020         {
1021         case epbcNONE:
1022             if (ir->rcoulomb == 0)
1023             {
1024                 nonbondedForce->setNonbondedMethod(NonbondedForce::NoCutoff);
1025             }
1026             else
1027             {
1028                 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
1029             }
1030             break;
1031         case epbcXYZ:
1032             switch (ir->coulombtype)
1033             {
1034             case eelCUT:
1035             case eelRF:
1036             case eelGRF:
1037             case eelRF_NEC:
1038             case eelRF_ZERO:
1039                 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
1040                 break;
1041
1042             case eelEWALD:
1043                 nonbondedForce->setNonbondedMethod(NonbondedForce::Ewald);
1044                 break;
1045
1046             case eelPME:
1047                 nonbondedForce->setNonbondedMethod(NonbondedForce::PME);
1048                 break;
1049
1050             default:
1051                 gmx_fatal(FARGS,"Internal error: you should not see this message, it means that the"
1052                           "electrosatics option check failed. Please report this error!");
1053             }
1054             sys->setDefaultPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0),
1055                                        Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));
1056             nonbondedForce->setCutoffDistance(ir->rcoulomb);
1057
1058             break;
1059         default:
1060             gmx_fatal(FARGS,"OpenMM supports only full periodic boundary conditions "
1061                               "(pbc = xyz), or none (pbc = no).");
1062         }
1063
1064
1065         /* Fix for PME and Ewald error tolerance
1066          *
1067                  *  OpenMM uses approximate formulas to calculate the Ewald parameter:
1068                  *  alpha = (1.0/cutoff)*sqrt(-log(2.0*tolerlance));
1069                  *  and the grid spacing for PME:
1070                  *  gridX = ceil(2*alpha*box[0][0]/3*(pow(tol, 0.2)))
1071                  *  gridY = ceil(2*alpha*box[1][1]/3*(pow(tol, 0.2)));
1072                  *  gridZ = ceil(2*alpha*box[2][2]/3*(pow(tol, 0.2)));
1073                  *
1074                  *
1075                  *  If the default ewald_rtol=1e-5 is used we silently adjust the value to the
1076                  *  OpenMM default of 5e-4 otherwise a warning is issued about the action taken.
1077                  *
1078                 */
1079         double corr_ewald_rtol = 50.0 * ir->ewald_rtol;
1080         if ((ir->ePBC == epbcXYZ) &&
1081             (ir->coulombtype == eelEWALD || ir->coulombtype == eelPME))
1082         {
1083             if (debug)
1084             {
1085                 fprintf(debug, ">> ewald_rtol = %e (corrected = %e) \n",
1086                     ir->ewald_rtol, corr_ewald_rtol);
1087             }
1088
1089             if (fabs(ir->ewald_rtol - 1e-5) > 1e-10)
1090             {
1091                 gmx_warning("OpenMM uses the ewald_rtol parameter with approximate formulas "
1092                         "to calculate the alpha and grid spacing parameters of the Ewald "
1093                         "and PME methods. This tolerance need to be corrected in order to get "
1094                         "settings close to the ones used in GROMACS. Although the internal correction "
1095                         "should work for any reasonable value of ewald_rtol, using values other than "
1096                         "the default 1e-5 might cause incorrect behavior.");
1097
1098                 if (corr_ewald_rtol > 1)
1099                 {
1100                     gmx_fatal(FARGS, "The ewald_rtol accuracy term is >1 after the "
1101                             "adjustment for OpenMM (%e)", corr_ewald_rtol);
1102                 }
1103             }
1104             nonbondedForce->setEwaldErrorTolerance(corr_ewald_rtol);
1105         }
1106
1107         for (int i = 0; i < numAtoms; ++i)
1108         {
1109             double c12 = nbfp[types[i]*2*ntypes+types[i]*2+1];
1110             double c6 = nbfp[types[i]*2*ntypes+types[i]*2];
1111             double sigma=0.0, epsilon=0.0;
1112             convert_c_12_6(c12, c6, &sigma, &epsilon);
1113             nonbondedForce->addParticle(charges[i], sigma, epsilon);
1114             sys->addParticle(masses[i]);
1115         }
1116
1117         // Build a table of all exclusions.
1118         vector<set<int> > exclusions(numAtoms);
1119         for (int i = 0; i < numAtoms; i++)
1120         {
1121             int start = top->excls.index[i];
1122             int end = top->excls.index[i+1];
1123             for (int j = start; j < end; j++)
1124                 exclusions[i].insert(top->excls.a[j]);
1125         }
1126
1127         // Record the 1-4 interactions, and remove them from the list of exclusions.
1128         const int* nb14Atoms = (int*) idef.il[F_LJ14].iatoms;
1129         offset = 0;
1130         for (int i = 0; i < num14; ++i)
1131         {
1132             int type = nb14Atoms[offset++];
1133             int atom1 = nb14Atoms[offset++];
1134             int atom2 = nb14Atoms[offset++];
1135             double sigma=0, epsilon=0;
1136             convert_c_12_6(idef.iparams[type].lj14.c12A,
1137                     idef.iparams[type].lj14.c6A,
1138                     &sigma, &epsilon);
1139             nonbondedForce->addException(atom1, atom2,
1140                                          fr->fudgeQQ*charges[atom1]*charges[atom2], sigma, epsilon);
1141             exclusions[atom1].erase(atom2);
1142             exclusions[atom2].erase(atom1);
1143         }
1144
1145         // Record exclusions.
1146         for (int i = 0; i < numAtoms; i++)
1147         {
1148             for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
1149             {
1150                 if (i < *iter)
1151                 {
1152                     nonbondedForce->addException(i, *iter, 0.0, 1.0, 0.0);
1153                 }
1154             }
1155         }
1156
1157         // Add GBSA if needed.
1158         if (ir->implicit_solvent == eisGBSA)
1159         {
1160             gmx_warning("The OBC scale factors alpha, beta and gamma are hardcoded in OpenMM with the default Gromacs values.");
1161             t_atoms atoms       = gmx_mtop_global_atoms(top_global);
1162             GBSAOBCForce* gbsa  = new GBSAOBCForce();
1163
1164             sys->addForce(gbsa);
1165             gbsa->setSoluteDielectric(ir->epsilon_r);
1166             gbsa->setSolventDielectric(ir->gb_epsilon_solvent);
1167             gbsa->setCutoffDistance(nonbondedForce->getCutoffDistance());
1168             if (nonbondedForce->getNonbondedMethod() == NonbondedForce::NoCutoff)
1169                 gbsa->setNonbondedMethod(GBSAOBCForce::NoCutoff);
1170             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
1171                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
1172             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffPeriodic)
1173                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
1174             else
1175                 gmx_fatal(FARGS,"OpenMM supports only Reaction-Field electrostatics with OBC/GBSA.");
1176
1177             for (int i = 0; i < numAtoms; ++i)
1178             {
1179                 gbsa->addParticle(charges[i],
1180                                   top_global->atomtypes.gb_radius[atoms.atom[i].type],
1181                                   top_global->atomtypes.S_hct[atoms.atom[i].type]);
1182             }
1183             free_t_atoms(&atoms, FALSE);
1184         }
1185
1186         // Set constraints.
1187         const int* constraintAtoms = (int*) idef.il[F_CONSTR].iatoms;
1188         offset = 0;
1189         for (int i = 0; i < numConstraints; ++i)
1190         {
1191             int type = constraintAtoms[offset++];
1192             int atom1 = constraintAtoms[offset++];
1193             int atom2 = constraintAtoms[offset++];
1194             sys->addConstraint(atom1, atom2, idef.iparams[type].constr.dA);
1195         }
1196         const int* settleAtoms = (int*) idef.il[F_SETTLE].iatoms;
1197         offset = 0;
1198         for (int i = 0; i < numSettle; ++i)
1199         {
1200             int type = settleAtoms[offset++];
1201             int oxygen = settleAtoms[offset++];
1202             sys->addConstraint(oxygen, oxygen+1, idef.iparams[type].settle.doh);
1203             sys->addConstraint(oxygen, oxygen+2, idef.iparams[type].settle.doh);
1204             sys->addConstraint(oxygen+1, oxygen+2, idef.iparams[type].settle.dhh);
1205         }
1206
1207         // Create an integrator for simulating the system.
1208         double friction = (ir->opts.tau_t[0] == 0.0 ? 0.0 : 1.0/ir->opts.tau_t[0]);
1209         Integrator* integ;
1210         if (ir->eI == eiBD)
1211         {
1212             integ = new BrownianIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1213             static_cast<BrownianIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1214         }
1215         else if (EI_SD(ir->eI))
1216         {
1217             integ = new LangevinIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1218             static_cast<LangevinIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1219         }
1220         else
1221         {
1222             integ = new VerletIntegrator(ir->delta_t);
1223             if ( ir->etc != etcNO)
1224             {
1225                 AndersenThermostat* thermostat = new AndersenThermostat(ir->opts.ref_t[0], friction);
1226                 sys->addForce(thermostat);
1227             }
1228         }
1229
1230                 // Add pressure coupling
1231         if (ir->epc != epcNO)
1232                 {
1233           // convert gromacs pressure tensor to a scalar
1234           double pressure = (ir->ref_p[0][0] + ir->ref_p[1][1] + ir->ref_p[2][2]) / 3.0;
1235           int frequency = int(ir->tau_p / ir->delta_t); // update frequency in time steps
1236           if (frequency < 1) frequency = 1;
1237           double temperature = ir->opts.ref_t[0]; // in kelvin
1238           sys->addForce(new MonteCarloBarostat(pressure, temperature, frequency));
1239                 }
1240
1241         integ->setConstraintTolerance(ir->shake_tol);
1242
1243         // Create a context and initialize it.
1244         Context* context = NULL;
1245
1246         /*
1247         OpenMM could automatically select the "best" GPU, however we're not't
1248         going to let it do that for now, as the current algorithm is very rudimentary
1249         and we anyway support only CUDA.
1250         if (platformOptStr == NULL || platformOptStr == "")
1251         {
1252             context = new Context(*sys, *integ);
1253         }
1254         else
1255         */
1256         {
1257             /* which platform should we use */
1258             for (int i = 0; i < (int)Platform::getNumPlatforms() && context == NULL; i++)
1259             {
1260                 if (isStringEqNCase(opt->getOptionValue("platform"), Platform::getPlatform(i).getName()))
1261                 {
1262                     Platform& platform = Platform::getPlatform(i);
1263                     // set standard properties
1264                     platform.setPropertyDefaultValue("CudaDevice", opt->getOptionValue("deviceid"));
1265                     // TODO add extra properties
1266                     context = new Context(*sys, *integ, platform);
1267                 }
1268             }
1269             if (context == NULL)
1270             {
1271                 gmx_fatal(FARGS, "The requested platform \"%s\" could not be found.",
1272                         opt->getOptionValue("platform").c_str());
1273             }
1274         }
1275
1276         Platform& platform = context->getPlatform();
1277         fprintf(fplog, "Gromacs will use the OpenMM platform: %s\n", platform.getName().c_str());
1278
1279         const vector<string>& properties = platform.getPropertyNames();
1280         if (debug)
1281         {
1282             for (int i = 0; i < (int)properties.size(); i++)
1283             {
1284                 fprintf(debug, ">> %s: %s\n", properties[i].c_str(),
1285                         platform.getPropertyValue(*context, properties[i]).c_str());
1286             }
1287         }
1288
1289         /* only for CUDA */
1290         if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1291         {
1292             int tmp;
1293             if (!from_string<int>(tmp, platform.getPropertyValue(*context, "CudaDevice"), std::dec))
1294             {
1295                 gmx_fatal(FARGS, "Internal error: couldn't determine the device selected by OpenMM");
1296
1297             }
1298
1299             /* For now this is just to double-check if OpenMM selected the GPU we wanted,
1300             but when we'll let OpenMM select the GPU automatically, it will query the deviceId.
1301             */
1302             if (tmp != devId)
1303             {
1304                 gmx_fatal(FARGS, "Internal error: OpenMM is using device #%d"
1305                         "while initialized for device #%d", tmp, devId);
1306             }
1307
1308             /* check GPU compatibility */
1309             char gpuname[STRLEN];
1310             devId = atoi(opt->getOptionValue("deviceid").c_str());
1311             if (!is_supported_cuda_gpu(-1, gpuname))
1312             {
1313                 if (!gmx_strcasecmp(opt->getOptionValue("force-device").c_str(), "yes"))
1314                 {
1315                     sprintf(warn_buf, "Non-supported GPU selected (#%d, %s), forced continuing."
1316                             "Note, that the simulation can be slow or it migth even crash.",
1317                             devId, gpuname);
1318                     fprintf(fplog, "%s\n", warn_buf);
1319                     gmx_warning(warn_buf);
1320                 }
1321                 else
1322                 {
1323                     gmx_fatal(FARGS, "The selected GPU (#%d, %s) is not supported by Gromacs! "
1324                               "Most probably you have a low-end GPU which would not perform well, "
1325                               "or new hardware that has not been tested with the current release. "
1326                               "If you still want to try using the device, use the force-device=yes option.",
1327                               devId, gpuname);
1328                 }
1329             }
1330             else
1331             {
1332                 fprintf(fplog, "Gromacs will run on the GPU #%d (%s).\n", devId, gpuname);
1333             }
1334         }
1335
1336         /* only for CUDA */
1337         if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1338         {
1339             /* pre-simulation memtest */
1340             runMemtest(fplog, -1, "Pre", opt);
1341         }
1342
1343         vector<Vec3> pos(numAtoms);
1344         vector<Vec3> vel(numAtoms);
1345         for (int i = 0; i < numAtoms; ++i)
1346         {
1347             pos[i] = Vec3(state->x[i][0], state->x[i][1], state->x[i][2]);
1348             vel[i] = Vec3(state->v[i][0], state->v[i][1], state->v[i][2]);
1349         }
1350         context->setPositions(pos);
1351         context->setVelocities(vel);
1352
1353         // Return a structure containing the system, integrator, and context.
1354         OpenMMData* data = new OpenMMData();
1355         data->system = sys;
1356         data->integrator = integ;
1357         data->context = context;
1358         data->removeCM = (ir->nstcomm > 0);
1359         data->platformOpt = opt;
1360         return data;
1361     }
1362     catch (std::exception& e)
1363     {
1364         gmx_fatal(FARGS, "OpenMM exception caught while initializating: %s", e.what());
1365     }
1366     return NULL; /* just to avoid warnings */
1367 }
1368
1369 /*!
1370  * \brief Integrate one step.
1371  *
1372  * \param[in] data  OpenMMData object created by openmm_init().
1373  */
1374 void openmm_take_one_step(void* data)
1375 {
1376     // static int step = 0; printf("----> taking step #%d\n", step++);
1377     try
1378     {
1379         static_cast<OpenMMData*>(data)->integrator->step(1);
1380     }
1381     catch (std::exception& e)
1382     {
1383         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1384     }
1385 }
1386
1387 /*!
1388  * \brief Integrate n steps.
1389  *
1390  * \param[in] data  OpenMMData object created by openmm_init().
1391  */
1392 void openmm_take_steps(void* data, int nstep)
1393 {
1394     try
1395     {
1396         static_cast<OpenMMData*>(data)->integrator->step(nstep);
1397     }
1398     catch (std::exception& e)
1399     {
1400         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1401     }
1402 }
1403
1404 /*!
1405  * \brief Clean up the data structures cretead for OpenMM.
1406  *
1407  * \param[in] log   Log file pointer.
1408  * \param[in] data  OpenMMData object created by openmm_init().
1409  */
1410 void openmm_cleanup(FILE* fplog, void* data)
1411 {
1412     OpenMMData* d = static_cast<OpenMMData*>(data);
1413     /* only for CUDA */
1414     if (isStringEqNCase(d->platformOpt->getOptionValue("platform"), "CUDA"))
1415     {
1416         /* post-simulation memtest */
1417         runMemtest(fplog, -1, "Post", d->platformOpt);
1418     }
1419     delete d->system;
1420     delete d->integrator;
1421     delete d->context;
1422     delete d->platformOpt;
1423     delete d;
1424 }
1425
1426 /*!
1427  * \brief Copy the current state information from OpenMM into the Gromacs data structures.
1428  *
1429  * This function results in the requested proprties to be copied from the
1430  * GPU to host. As this represents a bottleneck, the frequency of pulling data
1431  * should be minimized.
1432  *
1433  * \param[in]   data        OpenMMData object created by openmm_init().
1434  * \param[out]  time        Simulation time for which the state was created.
1435  * \param[out]  state       State of the system: coordinates and velocities.
1436  * \param[out]  f           Forces.
1437  * \param[out]  enerd       Energies.
1438  * \param[in]   includePos  True if coordinates are requested.
1439  * \param[in]   includeVel  True if velocities are requested.
1440  * \param[in]   includeForce True if forces are requested.
1441  * \param[in]   includeEnergy True if energies are requested.
1442  */
1443 void openmm_copy_state(void *data,
1444                        t_state *state, double *time,
1445                        rvec f[], gmx_enerdata_t *enerd,
1446                        gmx_bool includePos, gmx_bool includeVel, gmx_bool includeForce, gmx_bool includeEnergy)
1447 {
1448     int types = 0;
1449     if (includePos)
1450         types += State::Positions;
1451     if (includeVel)
1452         types += State::Velocities;
1453     if (includeForce)
1454         types += State::Forces;
1455     if (includeEnergy)
1456         types += State::Energy;
1457     if (types == 0)
1458         return;
1459     try
1460     {
1461         State currentState = static_cast<OpenMMData*>(data)->context->getState(types);
1462         int numAtoms =  static_cast<OpenMMData*>(data)->system->getNumParticles();
1463         if (includePos)
1464         {
1465             for (int i = 0; i < numAtoms; i++)
1466             {
1467                 Vec3 x = currentState.getPositions()[i];
1468                 state->x[i][0] = x[0];
1469                 state->x[i][1] = x[1];
1470                 state->x[i][2] = x[2];
1471             }
1472         }
1473         if (includeVel)
1474         {
1475             for (int i = 0; i < numAtoms; i++)
1476             {
1477                 Vec3 v = currentState.getVelocities()[i];
1478                 state->v[i][0] = v[0];
1479                 state->v[i][1] = v[1];
1480                 state->v[i][2] = v[2];
1481             }
1482         }
1483         if (includeForce)
1484         {
1485             for (int i = 0; i < numAtoms; i++)
1486             {
1487                 Vec3 force = currentState.getForces()[i];
1488                 f[i][0] = force[0];
1489                 f[i][1] = force[1];
1490                 f[i][2] = force[2];
1491             }
1492         }
1493         if (includeEnergy)
1494         {
1495             int numConstraints = static_cast<OpenMMData*>(data)->system->getNumConstraints();
1496             int dof = 3*numAtoms-numConstraints;
1497             if (static_cast<OpenMMData*>(data)->removeCM)
1498                 dof -= 3;
1499             enerd->term[F_EPOT] = currentState.getPotentialEnergy();
1500             enerd->term[F_EKIN] = currentState.getKineticEnergy();
1501             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1502             enerd->term[F_TEMP] = 2.0*enerd->term[F_EKIN]/dof/BOLTZ;
1503         }
1504         *time = currentState.getTime();
1505     }
1506     catch (std::exception& e)
1507     {
1508         gmx_fatal(FARGS, "OpenMM exception caught while retrieving state information: %s", e.what());
1509     }
1510 }