src/kernel/openmm_wrapper.cpp

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2010, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35
  36 /*
  37  * Note, that parts of this source code originate from the Simtk release
  38  * of OpenMM accelerated Gromacs, for more details see:
  39  * https://simtk.org/project/xml/downloads.xml?group_id=161#package_id600
  40  */
  41
  42 #ifdef HAVE_CONFIG_H
  43 #include <config.h>
  44 #endif
  45
  46 #include <types/simple.h>
  47 #include <cmath>
  48 #include <set>
  49 #include <iostream>
  50 #include <sstream>
  51 #include <fstream>
  52 #include <map>
  53 #include <vector>
  54 #include <cctype>
  55 #include <algorithm>
  56
  57 using namespace std;
  58
  59 #include "OpenMM.h"
  60
  61 #include "gmx_fatal.h"
  62 #include "typedefs.h"
  63 #include "mdrun.h"
  64 #include "physics.h"
  65 #include "string2.h"
  66 #include "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 venter of motion, false otherwise. */
 387     GmxOpenMMPlatformOptions *platformOpt; /*! Platform options. */
 388 };
 389
 390 /*!
 391  *  \brief Runs memtest on the GPU that has alreaby been initialized by OpenMM.
 392  *  \param[in] fplog    Pointer to gromacs log file.
 393  *  \param[in] devId    Device id of the GPU to run the test on.
 394                         Note: as OpenMM previously creates the context,for now this is always -1.
 395  *  \param[in] pre_post Contains either "Pre" or "Post" just to be able to differentiate in
 396  *                      stdout messages/log between memtest carried out before and after simulation.
 397  *  \param[in] opt      Pointer to platform options object.
 398  */
 399 static void runMemtest(FILE* fplog, int devId, const char* pre_post, GmxOpenMMPlatformOptions *opt)
 400 {
 401     char        strout_buf[STRLEN];
 402     int         which_test;
 403     int         res = 0;
 404     string      s = opt->getOptionValue("memtest");
 405     const char  *test_type = s.c_str();
 406
 407     if (!gmx_strcasecmp(test_type, "off"))
 408     {
 409         which_test = 0;
 410     }
 411     else
 412     {
 413         if (!gmx_strcasecmp(test_type, "full"))
 414         {
 415             which_test = 2;
 416         }
 417         else
 418         {
 419             from_string<int>(which_test, test_type, std::dec);
 420         }
 421     }
 422
 423     if (which_test < 0)
 424     {
 425         gmx_fatal(FARGS, "Amount of seconds for memetest is negative (%d). ", which_test);
 426     }
 427
 428     switch (which_test)
 429     {
 430         case 0: /* no memtest */
 431             sprintf(strout_buf, "%s-simulation GPU memtest skipped. Note, that faulty memory can cause "
 432                 "incorrect results!", pre_post);
 433             fprintf(fplog, "%s\n", strout_buf);
 434             gmx_warning(strout_buf);
 435             break; /* case 0 */
 436
 437         case 1: /* quick memtest */
 438             fprintf(fplog,  "%s-simulation %s GPU memtest in progress...\n", pre_post, test_type);
 439             fprintf(stdout, "\n%s-simulation %s GPU memtest in progress...", pre_post, test_type);
 440             fflush(fplog);
 441             fflush(stdout);
 442             res = do_quick_memtest(devId);
 443             break; /* case 1 */
 444
 445         case 2: /* full memtest */
 446             fprintf(fplog,  "%s-simulation %s memtest in progress...\n", pre_post, test_type);
 447             fprintf(stdout, "\n%s-simulation %s memtest in progress...", pre_post, test_type);
 448             fflush(fplog);
 449             fflush(stdout);
 450             res = do_full_memtest(devId);
 451             break; /* case 2 */
 452
 453         default: /* timed memtest */
 454             fprintf(fplog,  "%s-simulation ~%ds memtest in progress...\n", pre_post, which_test);
 455             fprintf(stdout, "\n%s-simulation ~%ds memtest in progress...", pre_post, which_test);
 456             fflush(fplog);
 457             fflush(stdout);
 458             res = do_timed_memtest(devId, which_test);
 459         }
 460
 461         if (which_test != 0)
 462         {
 463             if (res != 0)
 464             {
 465                 gmx_fatal(FARGS, MEM_ERR_MSG(pre_post));
 466             }
 467             else
 468             {
 469                 fprintf(fplog,  "Memory test completed without errors.\n");
 470                 fflush(fplog);
 471                 fprintf(stdout, "done, no errors detected\n");
 472                 fflush(stdout);
 473             }
 474         }
 475 }
 476
 477 /*!
 478  * \brief Convert Lennard-Jones parameters c12 and c6 to sigma and epsilon.
 479  *
 480  * \param[in] c12
 481  * \param[in] c6
 482  * \param[out] sigma
 483  * \param[out] epsilon
 484  */
 485 static void convert_c_12_6(double c12, double c6, double *sigma, double *epsilon)
 486 {
 487     if (c12 == 0 && c6 == 0)
 488     {
 489         *epsilon    = 0.0;
 490         *sigma      = 1.0;
 491     }
 492     else if (c12 > 0 && c6 > 0)
 493     {
 494         *epsilon    = (c6*c6)/(4.0*c12);
 495         *sigma      = pow(c12/c6, 1.0/6.0);
 496     }
 497     else
 498     {
 499         gmx_fatal(FARGS,"OpenMM only supports c6 > 0 and c12 > 0 or c6 = c12 = 0.");
 500     }
 501 }
 502
 503 /*!
 504  * \brief Does gromacs option checking.
 505  *
 506  * Checks the gromacs mdp options for features unsupported in OpenMM, case in which
 507  * interrupts the execution. It also warns the user about pecularities of OpenMM
 508  * implementations.
 509  * \param[in] fplog         Gromacs log file pointer.
 510  * \param[in] ir            Gromacs input parameters, see ::t_inputrec
 511  * \param[in] top           Gromacs node local topology, \see gmx_localtop_t
 512  * \param[in] state         Gromacs state structure \see ::t_state
 513  * \param[in] mdatoms       Gromacs atom parameters, \see ::t_mdatoms
 514  * \param[in] fr            \see ::t_forcerec
 515  * \param[in] state         Gromacs systems state, \see ::t_state
 516  */
 517 static void checkGmxOptions(FILE* fplog, GmxOpenMMPlatformOptions *opt,
 518                             t_inputrec *ir, gmx_localtop_t *top,
 519                             t_forcerec *fr, t_state *state)
 520 {
 521     char    warn_buf[STRLEN];
 522     int     i, j, natoms;
 523     double  c6, c12;
 524     double  sigma_ij=0, sigma_ji=0, sigma_ii=0, sigma_jj=0, sigma_comb;
 525     double  eps_ij=0, eps_ji=0, eps_ii=0, eps_jj=0, eps_comb;
 526
 527     /* Abort if unsupported critical options are present */
 528
 529     /* Integrator */
 530     if (ir->eI ==  eiMD)
 531     {
 532         gmx_warning( "OpenMM does not support leap-frog, will use velocity-verlet integrator.");
 533     }
 534
 535     if (    (ir->eI !=  eiMD)   &&
 536             (ir->eI !=  eiVV)   &&
 537             (ir->eI !=  eiVVAK) &&
 538             (ir->eI !=  eiSD1)  &&
 539             (ir->eI !=  eiSD2)  &&
 540             (ir->eI !=  eiBD) )
 541     {
 542         gmx_fatal(FARGS, "OpenMM supports only the following integrators: md/md-vv/md-vv-avek, sd/sd1, and bd.");
 543     }
 544
 545     /* Electroctstics */
 546     if (   !(ir->coulombtype == eelPME   ||
 547              EEL_RF(ir->coulombtype)     ||
 548              ir->coulombtype == eelRF    ||
 549              ir->coulombtype == eelEWALD ||
 550              // no-cutoff
 551              (ir->coulombtype == eelCUT && ir->rcoulomb == 0 &&  ir->rvdw == 0) ||
 552              // we could have cut-off combined with GBSA (openmm will use RF)
 553              ir->implicit_solvent == eisGBSA)   )
 554     {
 555         gmx_fatal(FARGS,"OpenMM supports only the following methods for electrostatics: "
 556                 "NoCutoff (i.e. rcoulomb = rvdw = 0 ),Reaction-Field, Ewald or PME.");
 557     }
 558
 559     if (EEL_RF(ir->coulombtype) && ir->epsilon_rf != 0)
 560     {
 561         // openmm has epsilon_rf=inf hard-coded
 562         gmx_warning("OpenMM will use a Reaction-Field epsilon of infinity instead of %g.",ir->epsilon_rf);
 563     }
 564
 565     if (ir->etc != etcNO &&
 566         ir->eI  != eiSD1 &&
 567         ir->eI  != eiSD2 &&
 568         ir->eI  != eiBD )
 569     {
 570         gmx_warning("OpenMM supports only Andersen thermostat with the md/md-vv/md-vv-avek integrators.");
 571     }
 572
 573     if (ir->implicit_solvent == eisGBSA &&
 574         ir->gb_algorithm != egbOBC  )
 575     {
 576         gmx_warning("OpenMM does not support the specified algorithm for Generalized Born, will use OBC instead.");
 577     }
 578
 579     if (ir->opts.ngtc > 1)
 580         gmx_fatal(FARGS,"OpenMM does not support multiple temperature coupling groups.");
 581
 582     if (ir->epc != epcNO)
 583         gmx_warning("OpenMM supports only Monte Carlo barostat for pressure coupling.");
 584
 585     if (ir->opts.annealing[0])
 586         gmx_fatal(FARGS,"OpenMM does not support simulated annealing.");
 587
 588     if (top->idef.il[F_CONSTR].nr > 0 && ir->eConstrAlg != econtSHAKE)
 589         gmx_warning("OpenMM provides contraints as a combination "
 590                     "of SHAKE, SETTLE and CCMA. Accuracy is based on the SHAKE tolerance set "
 591                     "by the \"shake_tol\" option.");
 592
 593     if (ir->nwall != 0)
 594         gmx_fatal(FARGS,"OpenMM does not support walls.");
 595
 596     if (ir->ePull != epullNO)
 597         gmx_fatal(FARGS,"OpenMM does not support pulling.");
 598
 599     /* check for interaction types */
 600     for (i = 0; i < F_EPOT; i++)
 601     {
 602         if (!(i == F_CONSTR ||
 603             i == F_SETTLE   ||
 604             i == F_BONDS    ||
 605             i == F_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.size() > 0)
 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.size() > 0)
 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.size() > 0)
 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 numUB = idef.il[F_UREY_BRADLEY].nr/4;
 855         const int numAngles = idef.il[F_ANGLES].nr/4;
 856         const int numPeriodic = idef.il[F_PDIHS].nr/5;
 857         const int numPeriodicImproper = idef.il[F_PIDIHS].nr/5;
 858         const int numRB = idef.il[F_RBDIHS].nr/5;
 859         const int numImproperDih = idef.il[F_IDIHS].nr/5;
 860         const int num14 = idef.il[F_LJ14].nr/3;
 861         System* sys = new System();
 862         if (ir->nstcomm > 0)
 863             sys->addForce(new CMMotionRemover(ir->nstcomm));
 864
 865         /* Set bonded force field terms. */
 866
 867                 /*
 868                  * CUDA platform currently doesn't support more than one
 869                  * instance of a force object, so we pack all forces that
 870                  * use the same form into one.
 871                 */
 872
 873         const int* bondAtoms = (int*) idef.il[F_BONDS].iatoms;
 874         HarmonicBondForce* bondForce = new HarmonicBondForce();
 875         sys->addForce(bondForce);
 876         int offset = 0;
 877         for (int i = 0; i < numBonds; ++i)
 878         {
 879             int type = bondAtoms[offset++];
 880             int atom1 = bondAtoms[offset++];
 881             int atom2 = bondAtoms[offset++];
 882             bondForce->addBond(atom1, atom2,
 883                                idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
 884         }
 885
 886                 /* Set the angle force field terms */
 887         const int* angleAtoms = (int*) idef.il[F_ANGLES].iatoms;
 888         HarmonicAngleForce* angleForce = new HarmonicAngleForce();
 889         sys->addForce(angleForce);
 890         offset = 0;
 891         for (int i = 0; i < numAngles; ++i)
 892         {
 893             int type = angleAtoms[offset++];
 894             int atom1 = angleAtoms[offset++];
 895             int atom2 = angleAtoms[offset++];
 896             int atom3 = angleAtoms[offset++];
 897             angleForce->addAngle(atom1, atom2, atom3,
 898                     idef.iparams[type].harmonic.rA*M_PI/180.0, idef.iparams[type].harmonic.krA);
 899         }
 900
 901         /* Urey-Bradley includes both the angle and bond potential for 1-3 interactions */
 902         const int* ubAtoms = (int*) idef.il[F_UREY_BRADLEY].iatoms;
 903                 /* HarmonicBondForce* ubBondForce = new HarmonicBondForce(); */
 904                 /*  HarmonicAngleForce* ubAngleForce = new HarmonicAngleForce(); */
 905         /* sys->addForce(ubBondForce); */
 906         /* sys->addForce(ubAngleForce); */
 907         offset = 0;
 908         for (int i = 0; i < numUB; ++i)
 909         {
 910             int type = ubAtoms[offset++];
 911             int atom1 = ubAtoms[offset++];
 912             int atom2 = ubAtoms[offset++];
 913             int atom3 = ubAtoms[offset++];
 914             /* ubBondForce->addBond(atom1, atom3, */
 915             bondForce->addBond(atom1, atom3,
 916                                idef.iparams[type].u_b.r13, idef.iparams[type].u_b.kUB);
 917             /* ubAngleForce->addAngle(atom1, atom2, atom3, */
 918             angleForce->addAngle(atom1, atom2, atom3,
 919                     idef.iparams[type].u_b.theta*M_PI/180.0, idef.iparams[type].u_b.ktheta);
 920         }
 921
 922                 /* Set proper dihedral terms */
 923         const int* periodicAtoms = (int*) idef.il[F_PDIHS].iatoms;
 924         PeriodicTorsionForce* periodicForce = new PeriodicTorsionForce();
 925         sys->addForce(periodicForce);
 926         offset = 0;
 927         for (int i = 0; i < numPeriodic; ++i)
 928         {
 929             int type = periodicAtoms[offset++];
 930             int atom1 = periodicAtoms[offset++];
 931             int atom2 = periodicAtoms[offset++];
 932             int atom3 = periodicAtoms[offset++];
 933             int atom4 = periodicAtoms[offset++];
 934             periodicForce->addTorsion(atom1, atom2, atom3, atom4,
 935                                       idef.iparams[type].pdihs.mult,
 936                                       idef.iparams[type].pdihs.phiA*M_PI/180.0,
 937                                       idef.iparams[type].pdihs.cpA);
 938         }
 939
 940                 /* Set improper dihedral terms that are represented by a periodic function (as in AMBER FF) */
 941         const int* periodicImproperAtoms = (int*) idef.il[F_PIDIHS].iatoms;
 942         /* PeriodicTorsionForce* periodicImproperForce = new PeriodicTorsionForce(); */
 943         /* sys->addForce(periodicImproperForce); */
 944         offset = 0;
 945         for (int i = 0; i < numPeriodicImproper; ++i)
 946         {
 947             int type = periodicImproperAtoms[offset++];
 948             int atom1 = periodicImproperAtoms[offset++];
 949             int atom2 = periodicImproperAtoms[offset++];
 950             int atom3 = periodicImproperAtoms[offset++];
 951             int atom4 = periodicImproperAtoms[offset++];
 952             /* periodicImproperForce->addTorsion(atom1, atom2, atom3, atom4, */
 953             periodicForce->addTorsion(atom1, atom2, atom3, atom4,
 954                                       idef.iparams[type].pdihs.mult,
 955                                       idef.iparams[type].pdihs.phiA*M_PI/180.0,
 956                                       idef.iparams[type].pdihs.cpA);
 957         }
 958
 959         /* Ryckaert-Bellemans dihedrals */
 960         const int* rbAtoms = (int*) idef.il[F_RBDIHS].iatoms;
 961         RBTorsionForce* rbForce = new RBTorsionForce();
 962         sys->addForce(rbForce);
 963         offset = 0;
 964         for (int i = 0; i < numRB; ++i)
 965         {
 966             int type = rbAtoms[offset++];
 967             int atom1 = rbAtoms[offset++];
 968             int atom2 = rbAtoms[offset++];
 969             int atom3 = rbAtoms[offset++];
 970             int atom4 = rbAtoms[offset++];
 971             rbForce->addTorsion(atom1, atom2, atom3, atom4,
 972                                 idef.iparams[type].rbdihs.rbcA[0], idef.iparams[type].rbdihs.rbcA[1],
 973                                 idef.iparams[type].rbdihs.rbcA[2], idef.iparams[type].rbdihs.rbcA[3],
 974                                 idef.iparams[type].rbdihs.rbcA[4], idef.iparams[type].rbdihs.rbcA[5]);
 975         }
 976
 977                 /* Set improper dihedral terms (as in CHARMM FF) */
 978         const int* improperDihAtoms = (int*) idef.il[F_IDIHS].iatoms;
 979                 CustomTorsionForce* improperDihForce = new CustomTorsionForce("2.0*k*asin(sin((theta-theta0)/2))^2");
 980         sys->addForce(improperDihForce);
 981                 improperDihForce->addPerTorsionParameter("k");
 982                 improperDihForce->addPerTorsionParameter("theta0");
 983                 vector<double> improperDihParameters(2);
 984         offset = 0;
 985         for (int i = 0; i < numImproperDih; ++i)
 986         {
 987             int type = improperDihAtoms[offset++];
 988             int atom1 = improperDihAtoms[offset++];
 989             int atom2 = improperDihAtoms[offset++];
 990             int atom3 = improperDihAtoms[offset++];
 991             int atom4 = improperDihAtoms[offset++];
 992                         improperDihParameters[0] = idef.iparams[type].harmonic.krA;
 993                         improperDihParameters[1] = idef.iparams[type].harmonic.rA*M_PI/180.0;
 994             improperDihForce->addTorsion(atom1, atom2, atom3, atom4,
 995                                 improperDihParameters);
 996         }
 997
 998         /* Set nonbonded parameters and masses. */
 999         int ntypes = fr->ntype;
1000         int* types = mdatoms->typeA;
1001         real* nbfp = fr->nbfp;
1002         real* charges = mdatoms->chargeA;
1003         real* masses = mdatoms->massT;
1004         NonbondedForce* nonbondedForce = new NonbondedForce();
1005         sys->addForce(nonbondedForce);
1006
1007         switch (ir->ePBC)
1008         {
1009         case epbcNONE:
1010             if (ir->rcoulomb == 0)
1011             {
1012                 nonbondedForce->setNonbondedMethod(NonbondedForce::NoCutoff);
1013             }
1014             else
1015             {
1016                 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
1017             }
1018             break;
1019         case epbcXYZ:
1020             switch (ir->coulombtype)
1021             {
1022             case eelRF:
1023             case eelGRF:
1024             case eelRF_NEC:
1025             case eelRF_ZERO:
1026                 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
1027                 break;
1028
1029             case eelEWALD:
1030                 nonbondedForce->setNonbondedMethod(NonbondedForce::Ewald);
1031                 break;
1032
1033             case eelPME:
1034                 nonbondedForce->setNonbondedMethod(NonbondedForce::PME);
1035                 break;
1036
1037             default:
1038                 gmx_fatal(FARGS,"Internal error: you should not see this message, it means that the"
1039                           "electrosatics option check failed. Please report this error!");
1040             }
1041             sys->setDefaultPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0),
1042                                        Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));
1043             nonbondedForce->setCutoffDistance(ir->rcoulomb);
1044
1045             break;
1046         default:
1047             gmx_fatal(FARGS,"OpenMM supports only full periodic boundary conditions "
1048                               "(pbc = xyz), or none (pbc = no).");
1049         }
1050
1051
1052         /* Fix for PME and Ewald error tolerance
1053          *
1054                  *  OpenMM uses approximate formulas to calculate the Ewald parameter:
1055                  *  alpha = (1.0/cutoff)*sqrt(-log(2.0*tolerlance));
1056                  *  and the grid spacing for PME:
1057                  *  gridX = ceil(2*alpha*box[0][0]/3*(pow(tol, 0.2)))
1058                  *  gridY = ceil(2*alpha*box[1][1]/3*(pow(tol, 0.2)));
1059                  *  gridZ = ceil(2*alpha*box[2][2]/3*(pow(tol, 0.2)));
1060                  *
1061                  *
1062                  *  If the default ewald_rtol=1e-5 is used we silently adjust the value to the
1063                  *  OpenMM default of 5e-4 otherwise a warning is issued about the action taken.
1064                  *
1065                 */
1066         double corr_ewald_rtol = 50.0 * ir->ewald_rtol;
1067         if ((ir->ePBC == epbcXYZ) &&
1068             (ir->coulombtype == eelEWALD || ir->coulombtype == eelPME))
1069         {
1070             if (debug)
1071             {
1072                 fprintf(debug, ">> ewald_rtol = %e (corrected = %e) \n",
1073                     ir->ewald_rtol, corr_ewald_rtol);
1074             }
1075
1076             if (fabs(ir->ewald_rtol - 1e-5) > 1e-10)
1077             {
1078                 gmx_warning("OpenMM uses the ewald_rtol parameter with approximate formulas "
1079                         "to calculate the alpha and grid spacing parameters of the Ewald "
1080                         "and PME methods. This tolerance need to be corrected in order to get "
1081                         "settings close to the ones used in GROMACS. Although the internal correction "
1082                         "should work for any reasonable value of ewald_rtol, using values other than "
1083                         "the default 1e-5 might cause incorrect behavior.");
1084
1085                 if (corr_ewald_rtol > 1)
1086                 {
1087                     gmx_fatal(FARGS, "The ewald_rtol accuracy term is >1 after the "
1088                             "adjustment for OpenMM (%e)", corr_ewald_rtol);
1089                 }
1090             }
1091             nonbondedForce->setEwaldErrorTolerance(corr_ewald_rtol);
1092         }
1093
1094         for (int i = 0; i < numAtoms; ++i)
1095         {
1096             double c12 = nbfp[types[i]*2*ntypes+types[i]*2+1];
1097             double c6 = nbfp[types[i]*2*ntypes+types[i]*2];
1098             double sigma=0.0, epsilon=0.0;
1099             convert_c_12_6(c12, c6, &sigma, &epsilon);
1100             nonbondedForce->addParticle(charges[i], sigma, epsilon);
1101             sys->addParticle(masses[i]);
1102         }
1103
1104         // Build a table of all exclusions.
1105         vector<set<int> > exclusions(numAtoms);
1106         for (int i = 0; i < numAtoms; i++)
1107         {
1108             int start = top->excls.index[i];
1109             int end = top->excls.index[i+1];
1110             for (int j = start; j < end; j++)
1111                 exclusions[i].insert(top->excls.a[j]);
1112         }
1113
1114         // Record the 1-4 interactions, and remove them from the list of exclusions.
1115         const int* nb14Atoms = (int*) idef.il[F_LJ14].iatoms;
1116         offset = 0;
1117         for (int i = 0; i < num14; ++i)
1118         {
1119             int type = nb14Atoms[offset++];
1120             int atom1 = nb14Atoms[offset++];
1121             int atom2 = nb14Atoms[offset++];
1122             double sigma=0, epsilon=0;
1123             convert_c_12_6(idef.iparams[type].lj14.c12A,
1124                     idef.iparams[type].lj14.c6A,
1125                     &sigma, &epsilon);
1126             nonbondedForce->addException(atom1, atom2,
1127                                          fr->fudgeQQ*charges[atom1]*charges[atom2], sigma, epsilon);
1128             exclusions[atom1].erase(atom2);
1129             exclusions[atom2].erase(atom1);
1130         }
1131
1132         // Record exclusions.
1133         for (int i = 0; i < numAtoms; i++)
1134         {
1135             for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
1136             {
1137                 if (i < *iter)
1138                 {
1139                     nonbondedForce->addException(i, *iter, 0.0, 1.0, 0.0);
1140                 }
1141             }
1142         }
1143
1144         // Add GBSA if needed.
1145         if (ir->implicit_solvent == eisGBSA)
1146         {
1147             gmx_warning("The OBC scale factors alpha, beta and gamma are hardcoded in OpenMM with the default Gromacs values.");
1148             t_atoms atoms       = gmx_mtop_global_atoms(top_global);
1149             GBSAOBCForce* gbsa  = new GBSAOBCForce();
1150
1151             sys->addForce(gbsa);
1152             gbsa->setSoluteDielectric(ir->epsilon_r);
1153             gbsa->setSolventDielectric(ir->gb_epsilon_solvent);
1154             gbsa->setCutoffDistance(nonbondedForce->getCutoffDistance());
1155             if (nonbondedForce->getNonbondedMethod() == NonbondedForce::NoCutoff)
1156                 gbsa->setNonbondedMethod(GBSAOBCForce::NoCutoff);
1157             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
1158                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
1159             else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffPeriodic)
1160                 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
1161             else
1162                 gmx_fatal(FARGS,"OpenMM supports only Reaction-Field electrostatics with OBC/GBSA.");
1163
1164             for (int i = 0; i < numAtoms; ++i)
1165             {
1166                 gbsa->addParticle(charges[i],
1167                                   top_global->atomtypes.gb_radius[atoms.atom[i].type],
1168                                   top_global->atomtypes.S_hct[atoms.atom[i].type]);
1169             }
1170             free_t_atoms(&atoms, FALSE);
1171         }
1172
1173         // Set constraints.
1174         const int* constraintAtoms = (int*) idef.il[F_CONSTR].iatoms;
1175         offset = 0;
1176         for (int i = 0; i < numConstraints; ++i)
1177         {
1178             int type = constraintAtoms[offset++];
1179             int atom1 = constraintAtoms[offset++];
1180             int atom2 = constraintAtoms[offset++];
1181             sys->addConstraint(atom1, atom2, idef.iparams[type].constr.dA);
1182         }
1183         const int* settleAtoms = (int*) idef.il[F_SETTLE].iatoms;
1184         offset = 0;
1185         for (int i = 0; i < numSettle; ++i)
1186         {
1187             int type = settleAtoms[offset++];
1188             int oxygen = settleAtoms[offset++];
1189             sys->addConstraint(oxygen, oxygen+1, idef.iparams[type].settle.doh);
1190             sys->addConstraint(oxygen, oxygen+2, idef.iparams[type].settle.doh);
1191             sys->addConstraint(oxygen+1, oxygen+2, idef.iparams[type].settle.dhh);
1192         }
1193
1194         // Create an integrator for simulating the system.
1195         double friction = (ir->opts.tau_t[0] == 0.0 ? 0.0 : 1.0/ir->opts.tau_t[0]);
1196         Integrator* integ;
1197         if (ir->eI == eiBD)
1198         {
1199             integ = new BrownianIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1200             static_cast<BrownianIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1201         }
1202         else if (EI_SD(ir->eI))
1203         {
1204             integ = new LangevinIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1205             static_cast<LangevinIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1206         }
1207         else
1208         {
1209             integ = new VerletIntegrator(ir->delta_t);
1210             if ( ir->etc != etcNO)
1211             {
1212                 AndersenThermostat* thermostat = new AndersenThermostat(ir->opts.ref_t[0], friction);
1213                 sys->addForce(thermostat);
1214             }
1215         }
1216
1217                 // Add pressure coupling
1218         if (ir->epc != epcNO)
1219                 {
1220           // convert gromacs pressure tensor to a scalar
1221           double pressure = (ir->ref_p[0][0] + ir->ref_p[1][1] + ir->ref_p[2][2]) / 3.0;
1222           int frequency = int(ir->tau_p / ir->delta_t); // update frequency in time steps
1223           if (frequency < 1) frequency = 1;
1224           double temperature = ir->opts.ref_t[0]; // in kelvin
1225           sys->addForce(new MonteCarloBarostat(pressure, temperature, frequency));
1226                 }
1227
1228         integ->setConstraintTolerance(ir->shake_tol);
1229
1230         // Create a context and initialize it.
1231         Context* context = NULL;
1232
1233         /*
1234         OpenMM could automatically select the "best" GPU, however we're not't
1235         going to let it do that for now, as the current algorithm is very rudimentary
1236         and we anyway support only CUDA.
1237         if (platformOptStr == NULL || platformOptStr == "")
1238         {
1239             context = new Context(*sys, *integ);
1240         }
1241         else
1242         */
1243         {
1244             /* which platform should we use */
1245             for (int i = 0; i < (int)Platform::getNumPlatforms() && context == NULL; i++)
1246             {
1247                 if (isStringEqNCase(opt->getOptionValue("platform"), Platform::getPlatform(i).getName()))
1248                 {
1249                     Platform& platform = Platform::getPlatform(i);
1250                     // set standard properties
1251                     platform.setPropertyDefaultValue("CudaDevice", opt->getOptionValue("deviceid"));
1252                     // TODO add extra properties
1253                     context = new Context(*sys, *integ, platform);
1254                 }
1255             }
1256             if (context == NULL)
1257             {
1258                 gmx_fatal(FARGS, "The requested platform \"%s\" could not be found.",
1259                         opt->getOptionValue("platform").c_str());
1260             }
1261         }
1262
1263         Platform& platform = context->getPlatform();
1264         fprintf(fplog, "Gromacs will use the OpenMM platform: %s\n", platform.getName().c_str());
1265
1266         const vector<string>& properties = platform.getPropertyNames();
1267         if (debug)
1268         {
1269             for (int i = 0; i < (int)properties.size(); i++)
1270             {
1271                 fprintf(debug, ">> %s: %s\n", properties[i].c_str(),
1272                         platform.getPropertyValue(*context, properties[i]).c_str());
1273             }
1274         }
1275
1276         /* only for CUDA */
1277         if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1278         {
1279             int tmp;
1280             if (!from_string<int>(tmp, platform.getPropertyValue(*context, "CudaDevice"), std::dec))
1281             {
1282                 gmx_fatal(FARGS, "Internal error: couldn't determine the device selected by OpenMM");
1283
1284             }
1285
1286             /* For now this is just to double-check if OpenMM selected the GPU we wanted,
1287             but when we'll let OpenMM select the GPU automatically, it will query the deviceId.
1288             */
1289             if (tmp != devId)
1290             {
1291                 gmx_fatal(FARGS, "Internal error: OpenMM is using device #%d"
1292                         "while initialized for device #%d", tmp, devId);
1293             }
1294
1295             /* check GPU compatibility */
1296             char gpuname[STRLEN];
1297             devId = atoi(opt->getOptionValue("deviceid").c_str());
1298             if (!is_supported_cuda_gpu(-1, gpuname))
1299             {
1300                 if (!gmx_strcasecmp(opt->getOptionValue("force-device").c_str(), "yes"))
1301                 {
1302                     sprintf(warn_buf, "Non-supported GPU selected (#%d, %s), forced continuing."
1303                             "Note, that the simulation can be slow or it migth even crash.",
1304                             devId, gpuname);
1305                     fprintf(fplog, "%s\n", warn_buf);
1306                     gmx_warning(warn_buf);
1307                 }
1308                 else
1309                 {
1310                     gmx_fatal(FARGS, "The selected GPU (#%d, %s) is not supported by Gromacs! "
1311                               "Most probably you have a low-end GPU which would not perform well, "
1312                               "or new hardware that has not been tested with the current release. "
1313                               "If you still want to try using the device, use the force-device=yes option.",
1314                               devId, gpuname);
1315                 }
1316             }
1317             else
1318             {
1319                 fprintf(fplog, "Gromacs will run on the GPU #%d (%s).\n", devId, gpuname);
1320             }
1321         }
1322
1323         /* only for CUDA */
1324         if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1325         {
1326             /* pre-simulation memtest */
1327             runMemtest(fplog, -1, "Pre", opt);
1328         }
1329
1330         vector<Vec3> pos(numAtoms);
1331         vector<Vec3> vel(numAtoms);
1332         for (int i = 0; i < numAtoms; ++i)
1333         {
1334             pos[i] = Vec3(state->x[i][0], state->x[i][1], state->x[i][2]);
1335             vel[i] = Vec3(state->v[i][0], state->v[i][1], state->v[i][2]);
1336         }
1337         context->setPositions(pos);
1338         context->setVelocities(vel);
1339
1340         // Return a structure containing the system, integrator, and context.
1341         OpenMMData* data = new OpenMMData();
1342         data->system = sys;
1343         data->integrator = integ;
1344         data->context = context;
1345         data->removeCM = (ir->nstcomm > 0);
1346         data->platformOpt = opt;
1347         return data;
1348     }
1349     catch (std::exception& e)
1350     {
1351         gmx_fatal(FARGS, "OpenMM exception caught while initializating: %s", e.what());
1352     }
1353     return NULL; /* just to avoid warnings */
1354 }
1355
1356 /*!
1357  * \brief Integrate one step.
1358  *
1359  * \param[in] data  OpenMMData object created by openmm_init().
1360  */
1361 void openmm_take_one_step(void* data)
1362 {
1363     // static int step = 0; printf("----> taking step #%d\n", step++);
1364     try
1365     {
1366         static_cast<OpenMMData*>(data)->integrator->step(1);
1367     }
1368     catch (std::exception& e)
1369     {
1370         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1371     }
1372 }
1373
1374 /*!
1375  * \brief Integrate n steps.
1376  *
1377  * \param[in] data  OpenMMData object created by openmm_init().
1378  */
1379 void openmm_take_steps(void* data, int nstep)
1380 {
1381     try
1382     {
1383         static_cast<OpenMMData*>(data)->integrator->step(nstep);
1384     }
1385     catch (std::exception& e)
1386     {
1387         gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1388     }
1389 }
1390
1391 /*!
1392  * \brief Clean up the data structures cretead for OpenMM.
1393  *
1394  * \param[in] log   Log file pointer.
1395  * \param[in] data  OpenMMData object created by openmm_init().
1396  */
1397 void openmm_cleanup(FILE* fplog, void* data)
1398 {
1399     OpenMMData* d = static_cast<OpenMMData*>(data);
1400     /* only for CUDA */
1401     if (isStringEqNCase(d->platformOpt->getOptionValue("platform"), "CUDA"))
1402     {
1403         /* post-simulation memtest */
1404         runMemtest(fplog, -1, "Post", d->platformOpt);
1405     }
1406     delete d->system;
1407     delete d->integrator;
1408     delete d->context;
1409     delete d->platformOpt;
1410     delete d;
1411 }
1412
1413 /*!
1414  * \brief Copy the current state information from OpenMM into the Gromacs data structures.
1415  *
1416  * This function results in the requested proprties to be copied from the
1417  * GPU to host. As this represents a bottleneck, the frequency of pulling data
1418  * should be minimized.
1419  *
1420  * \param[in]   data        OpenMMData object created by openmm_init().
1421  * \param[out]  time        Simulation time for which the state was created.
1422  * \param[out]  state       State of the system: coordinates and velocities.
1423  * \param[out]  f           Forces.
1424  * \param[out]  enerd       Energies.
1425  * \param[in]   includePos  True if coordinates are requested.
1426  * \param[in]   includeVel  True if velocities are requested.
1427  * \param[in]   includeForce True if forces are requested.
1428  * \param[in]   includeEnergy True if energies are requested.
1429  */
1430 void openmm_copy_state(void *data,
1431                        t_state *state, double *time,
1432                        rvec f[], gmx_enerdata_t *enerd,
1433                        gmx_bool includePos, gmx_bool includeVel, gmx_bool includeForce, gmx_bool includeEnergy)
1434 {
1435     int types = 0;
1436     if (includePos)
1437         types += State::Positions;
1438     if (includeVel)
1439         types += State::Velocities;
1440     if (includeForce)
1441         types += State::Forces;
1442     if (includeEnergy)
1443         types += State::Energy;
1444     if (types == 0)
1445         return;
1446     try
1447     {
1448         State currentState = static_cast<OpenMMData*>(data)->context->getState(types);
1449         int numAtoms =  static_cast<OpenMMData*>(data)->system->getNumParticles();
1450         if (includePos)
1451         {
1452             for (int i = 0; i < numAtoms; i++)
1453             {
1454                 Vec3 x = currentState.getPositions()[i];
1455                 state->x[i][0] = x[0];
1456                 state->x[i][1] = x[1];
1457                 state->x[i][2] = x[2];
1458             }
1459         }
1460         if (includeVel)
1461         {
1462             for (int i = 0; i < numAtoms; i++)
1463             {
1464                 Vec3 v = currentState.getVelocities()[i];
1465                 state->v[i][0] = v[0];
1466                 state->v[i][1] = v[1];
1467                 state->v[i][2] = v[2];
1468             }
1469         }
1470         if (includeForce)
1471         {
1472             for (int i = 0; i < numAtoms; i++)
1473             {
1474                 Vec3 force = currentState.getForces()[i];
1475                 f[i][0] = force[0];
1476                 f[i][1] = force[1];
1477                 f[i][2] = force[2];
1478             }
1479         }
1480         if (includeEnergy)
1481         {
1482             int numConstraints = static_cast<OpenMMData*>(data)->system->getNumConstraints();
1483             int dof = 3*numAtoms-numConstraints;
1484             if (static_cast<OpenMMData*>(data)->removeCM)
1485                 dof -= 3;
1486             enerd->term[F_EPOT] = currentState.getPotentialEnergy();
1487             enerd->term[F_EKIN] = currentState.getKineticEnergy();
1488             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1489             enerd->term[F_TEMP] = 2.0*enerd->term[F_EKIN]/dof/BOLTZ;
1490         }
1491         *time = currentState.getTime();
1492     }
1493     catch (std::exception& e)
1494     {
1495         gmx_fatal(FARGS, "OpenMM exception caught while retrieving state information: %s", e.what());
1496     }
1497 }