1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
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
56 #include "gmx_fatal.h"
61 #include "gmx_gpu_utils.h"
62 #include "mtop_util.h"
64 #include "openmm_wrapper.h"
66 using namespace OpenMM;
69 #define MEM_ERR_MSG(str) \
70 "The %s-simulation GPU memory test detected errors. As memory errors would cause incorrect " \
71 "simulation results, gromacs has aborted execution.\n Make sure that your GPU's memory is not " \
72 "overclocked and that the device is properly cooled.\n", (str)
75 #define COMBRULE_CHK_TOL 1e-6
76 #define COMBRULE_SIGMA(sig1, sig2) (((sig1) + (sig2))/2)
77 #define COMBRULE_EPS(eps1, eps2) (sqrt((eps1) * (eps2)))
80 * \brief Convert string to integer type.
81 * \param[in] s String to convert from.
82 * \param[in] f Basefield format flag that takes any of the following I/O
83 * manipulators: dec, hex, oct.
84 * \param[out] t Destination variable to convert to.
87 static bool from_string(T& t, const string& s, ios_base& (*f)(ios_base&))
90 return !(iss >> f >> t).fail();
94 * \brief Split string around a given delimiter.
95 * \param[in] s String to split.
96 * \param[in] delim Delimiter character.
97 * \returns Vector of strings found in \p s.
99 static vector<string> split(const string &s, char delim)
101 vector<string> elems;
104 while (getline(ss, item, delim))
106 if (item.length() != 0)
107 elems.push_back(item);
113 * \brief Split a string of the form "option=value" into "option" and "value" strings.
114 * This string corresponds to one option and the associated value from the option list
115 * in the mdrun -device argument.
117 * \param[in] s A string containing an "option=value" pair that needs to be split up.
118 * \param[out] opt The name of the option.
119 * \param[out] val Value of the option.
121 static void splitOptionValue(const string &s, string &opt, string &val)
123 size_t eqPos = s.find('=');
124 if (eqPos != string::npos)
126 opt = s.substr(0, eqPos);
127 if (eqPos != s.length()) val = s.substr(eqPos+1);
132 * \brief Compare two strings ignoring case.
133 * This function is in fact a wrapper around the gromacs function gmx_strncasecmp().
134 * \param[in] s1 String.
135 * \param[in] s2 String.
136 * \returns Similarly to the C function strncasecmp(), the return value is an
137 integer less than, equal to, or greater than 0 if \p s1 less than,
138 identical to, or greater than \p s2.
140 static bool isStringEqNCase(const string s1, const string s2)
142 return (gmx_strncasecmp(s1.c_str(), s2.c_str(), max(s1.length(), s2.length())) == 0);
146 * \brief Convert string to upper case.
148 * \param[in] s String to convert to uppercase.
149 * \returns The given string converted to uppercase.
151 static string toUpper(const string &s)
154 std::transform(stmp.begin(), stmp.end(), stmp.begin(), static_cast < int(*)(int) > (toupper));
159 \name Sizes of constant device option arrays GmxOpenMMPlatformOptions#platforms,
160 GmxOpenMMPlatformOptions#memtests, GmxOpenMMPlatformOptions#deviceid,
161 GmxOpenMMPlatformOptions#force_dev. */
163 #define SIZEOF_PLATFORMS 2 // 2
164 #define SIZEOF_MEMTESTS 3
165 #define SIZEOF_DEVICEIDS 1
166 #define SIZEOF_FORCE_DEV 2
168 #define SIZEOF_CHECK_COMBRULE 2
171 /*! Possible platform options in the mdrun -device option. */
172 static const char *devOptStrings[] = { "platform", "deviceid", "memtest", "force-device", "check-combrule" };
174 /*! Enumerated platform options in the mdrun -device option. */
184 * \brief Class to extract and manage the platform options in the mdrun -device option.
187 class GmxOpenMMPlatformOptions
190 GmxOpenMMPlatformOptions(const char *opt);
191 ~GmxOpenMMPlatformOptions() { options.clear(); }
192 string getOptionValue(const string &opt);
193 void remOption(const string &opt);
196 void setOption(const string &opt, const string &val);
198 map<string, string> options; /*!< Data structure to store the option (name, value) pairs. */
200 static const char * const platforms[SIZEOF_PLATFORMS]; /*!< Available OpenMM platforms; size #SIZEOF_PLATFORMS */
201 static const char * const memtests[SIZEOF_MEMTESTS]; /*!< Available types of memory tests, also valid
202 any positive integer >=15; size #SIZEOF_MEMTESTS */
203 static const char * const deviceid[SIZEOF_DEVICEIDS]; /*!< Possible values for deviceid option;
204 also valid any positive integer; size #SIZEOF_DEVICEIDS */
205 static const char * const force_dev[SIZEOF_FORCE_DEV]; /*!< Possible values for for force-device option;
206 size #SIZEOF_FORCE_DEV */
207 static const char * const check_combrule[SIZEOF_CHECK_COMBRULE]; /* XXX temporary debug feature to
208 turn off combination rule check */
211 const char * const GmxOpenMMPlatformOptions::platforms[SIZEOF_PLATFORMS]
212 = {"CUDA", "Reference"};
213 //= { "Reference", "CUDA" /*,"OpenCL"*/ };
214 const char * const GmxOpenMMPlatformOptions::memtests[SIZEOF_MEMTESTS]
215 = { "15", "full", "off" };
216 const char * const GmxOpenMMPlatformOptions::deviceid[SIZEOF_DEVICEIDS]
218 const char * const GmxOpenMMPlatformOptions::force_dev[SIZEOF_FORCE_DEV]
220 const char * const GmxOpenMMPlatformOptions::check_combrule[SIZEOF_CHECK_COMBRULE]
225 * Takes the option list, parses it, checks the options and their values for validity.
226 * When certain options are not provided by the user, as default value the first item
227 * of the respective constant array is taken (GmxOpenMMPlatformOptions#platforms,
228 * GmxOpenMMPlatformOptions#memtests, GmxOpenMMPlatformOptions#deviceid,
229 * GmxOpenMMPlatformOptions#force_dev).
230 * \param[in] optionString Option list part of the mdrun -device parameter.
232 GmxOpenMMPlatformOptions::GmxOpenMMPlatformOptions(const char *optionString)
234 // set default values
235 setOption("platform", platforms[0]);
236 setOption("memtest", memtests[0]);
237 setOption("deviceid", deviceid[0]);
238 setOption("force-device", force_dev[0]);
239 setOption("check-combrule", check_combrule[0]);
241 string opt(optionString);
243 // remove all whitespaces
244 opt.erase(remove_if(opt.begin(), opt.end(), ::isspace), opt.end());
245 // tokenize around ","-s
246 vector<string> tokens = split(opt, ',');
248 for (vector<string>::iterator it = tokens.begin(); it != tokens.end(); ++it)
250 string opt = "", val = "";
251 splitOptionValue(*it, opt, val);
253 if (isStringEqNCase(opt, "platform"))
255 /* no check, this will fail if platform does not exist when we try to set it */
260 if (isStringEqNCase(opt, "memtest"))
262 /* the value has to be an integer >15(s) or "full" OR "off" */
263 if (!isStringEqNCase(val, "full") && !isStringEqNCase(val, "off"))
266 if (!from_string<int>(secs, val, std::dec))
268 gmx_fatal(FARGS, "Invalid value for option memtest option: \"%s\"!", val.c_str());
272 gmx_fatal(FARGS, "Incorrect value for memtest option (%d). "
273 "Memtest needs to run for at least 15s!", secs);
280 if (isStringEqNCase(opt, "deviceid"))
283 if (!from_string<int>(id, val, std::dec) )
285 gmx_fatal(FARGS, "Invalid device id: \"%s\"!", val.c_str());
291 if (isStringEqNCase(opt, "force-device"))
294 if (!isStringEqNCase(val, "yes") && !isStringEqNCase(val, "no"))
296 gmx_fatal(FARGS, "Invalid OpenMM force option: \"%s\"!", val.c_str());
302 if (isStringEqNCase(opt, "check-combrule"))
305 if (!isStringEqNCase(val, "yes") && !isStringEqNCase(val, "no"))
307 gmx_fatal(FARGS, "Invalid OpenMM force option: \"%s\"!", val.c_str());
314 // if we got till here something went wrong
315 gmx_fatal(FARGS, "Invalid OpenMM platform option: \"%s\"!", (*it).c_str());
321 * \brief Getter function.
322 * \param[in] opt Name of the option.
323 * \returns Returns the value associated to an option.
325 string GmxOpenMMPlatformOptions::getOptionValue(const string &opt)
327 map<string, string> :: const_iterator it = options.find(toUpper(opt));
328 if (it != options.end())
339 * \brief Setter function - private, only used from contructor.
340 * \param[in] opt Name of the option.
341 * \param[in] val Value for the option.
343 void GmxOpenMMPlatformOptions::setOption(const string &opt, const string &val)
345 options[toUpper(opt)] = val;
349 * \brief Removes an option with its value from the map structure. If the option
350 * does not exist, returns without any action.
351 * \param[in] opt Name of the option.
353 void GmxOpenMMPlatformOptions::remOption(const string &opt)
355 options.erase(toUpper(opt));
359 * \brief Print option-value pairs to a file (debugging function).
361 void GmxOpenMMPlatformOptions::print()
363 cout << ">> Platform options: " << endl
364 << ">> platform = " << getOptionValue("platform") << endl
365 << ">> deviceID = " << getOptionValue("deviceid") << endl
366 << ">> memtest = " << getOptionValue("memtest") << endl
367 << ">> force-device = " << getOptionValue("force-device") << endl;
371 * \brief Container for OpenMM related data structures that represent the bridge
372 * between the Gromacs data-structures and the OpenMM library and is but it's
373 * only passed through the API functions as void to disable direct access.
378 System* system; /*! The system to simulate. */
379 Context* context; /*! The OpenMM context in which the simulation is carried out. */
380 Integrator* integrator; /*! The integrator used in the simulation. */
381 bool removeCM; /*! If \true remove venter of motion, false otherwise. */
382 GmxOpenMMPlatformOptions *platformOpt; /*! Platform options. */
386 * \brief Runs memtest on the GPU that has alreaby been initialized by OpenMM.
387 * \param[in] fplog Pointer to gromacs log file.
388 * \param[in] devId Device id of the GPU to run the test on.
389 Note: as OpenMM previously creates the context,for now this is always -1.
390 * \param[in] pre_post Contains either "Pre" or "Post" just to be able to differentiate in
391 * stdout messages/log between memtest carried out before and after simulation.
392 * \param[in] opt Pointer to platform options object.
394 static void runMemtest(FILE* fplog, int devId, const char* pre_post, GmxOpenMMPlatformOptions *opt)
396 char strout_buf[STRLEN];
399 string s = opt->getOptionValue("memtest");
400 const char *test_type = s.c_str();
402 if (!gmx_strcasecmp(test_type, "off"))
408 if (!gmx_strcasecmp(test_type, "full"))
414 from_string<int>(which_test, test_type, std::dec);
420 gmx_fatal(FARGS, "Amount of seconds for memetest is negative (%d). ", which_test);
425 case 0: /* no memtest */
426 sprintf(strout_buf, "%s-simulation GPU memtest skipped. Note, that faulty memory can cause "
427 "incorrect results!", pre_post);
428 fprintf(fplog, "%s\n", strout_buf);
429 gmx_warning(strout_buf);
432 case 1: /* quick memtest */
433 fprintf(fplog, "%s-simulation %s GPU memtest in progress...\n", pre_post, test_type);
434 fprintf(stdout, "\n%s-simulation %s GPU memtest in progress...", pre_post, test_type);
437 res = do_quick_memtest(devId);
440 case 2: /* full memtest */
441 fprintf(fplog, "%s-simulation %s memtest in progress...\n", pre_post, test_type);
442 fprintf(stdout, "\n%s-simulation %s memtest in progress...", pre_post, test_type);
445 res = do_full_memtest(devId);
448 default: /* timed memtest */
449 fprintf(fplog, "%s-simulation ~%ds memtest in progress...\n", pre_post, which_test);
450 fprintf(stdout, "\n%s-simulation ~%ds memtest in progress...", pre_post, which_test);
453 res = do_timed_memtest(devId, which_test);
460 gmx_fatal(FARGS, MEM_ERR_MSG(pre_post));
464 fprintf(fplog, "Memory test completed without errors.\n");
466 fprintf(stdout, "done, no errors detected\n");
473 * \brief Convert Lennard-Jones parameters c12 and c6 to sigma and epsilon.
478 * \param[out] epsilon
480 static void convert_c_12_6(double c12, double c6, double *sigma, double *epsilon)
482 if (c12 == 0 && c6 == 0)
487 else if (c12 > 0 && c6 > 0)
489 *epsilon = (c6*c6)/(4.0*c12);
490 *sigma = pow(c12/c6, 1.0/6.0);
494 gmx_fatal(FARGS,"OpenMM only supports c6 > 0 and c12 > 0 or c6 = c12 = 0.");
499 * \brief Does gromacs option checking.
501 * Checks the gromacs mdp options for features unsupported in OpenMM, case in which
502 * interrupts the execution. It also warns the user about pecularities of OpenMM
504 * \param[in] fplog Gromacs log file pointer.
505 * \param[in] ir Gromacs input parameters, see ::t_inputrec
506 * \param[in] top Gromacs node local topology, \see gmx_localtop_t
507 * \param[in] state Gromacs state structure \see ::t_state
508 * \param[in] mdatoms Gromacs atom parameters, \see ::t_mdatoms
509 * \param[in] fr \see ::t_forcerec
510 * \param[in] state Gromacs systems state, \see ::t_state
512 static void checkGmxOptions(FILE* fplog, GmxOpenMMPlatformOptions *opt,
513 t_inputrec *ir, gmx_localtop_t *top,
514 t_forcerec *fr, t_state *state)
516 char warn_buf[STRLEN];
519 double sigma_ij=0, sigma_ji=0, sigma_ii=0, sigma_jj=0, sigma_comb;
520 double eps_ij=0, eps_ji=0, eps_ii=0, eps_jj=0, eps_comb;
522 /* Abort if unsupported critical options are present */
527 gmx_warning( "OpenMM does not support leap-frog, will use velocity-verlet integrator.");
530 if ( (ir->eI != eiMD) &&
532 (ir->eI != eiVVAK) &&
537 gmx_fatal(FARGS, "OpenMM supports only the following integrators: md/md-vv/md-vv-avek, sd/sd1, and bd.");
541 if ( !(ir->coulombtype == eelPME ||
542 EEL_RF(ir->coulombtype) ||
543 ir->coulombtype == eelRF ||
544 ir->coulombtype == eelEWALD ||
546 (ir->coulombtype == eelCUT && ir->rcoulomb == 0 && ir->rvdw == 0) ||
547 // we could have cut-off combined with GBSA (openmm will use RF)
548 ir->implicit_solvent == eisGBSA) )
550 gmx_fatal(FARGS,"OpenMM supports only the following methods for electrostatics: "
551 "NoCutoff (i.e. rcoulomb = rvdw = 0 ),Reaction-Field, Ewald or PME.");
554 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf != 0)
556 // openmm has epsilon_rf=inf hard-coded
557 gmx_warning("OpenMM will use a Reaction-Field epsilon of infinity instead of %g.",ir->epsilon_rf);
560 if (ir->etc != etcNO &&
565 gmx_warning("OpenMM supports only Andersen thermostat with the md/md-vv/md-vv-avek integrators.");
568 if (ir->implicit_solvent == eisGBSA &&
569 ir->gb_algorithm != egbOBC )
571 gmx_warning("OpenMM does not support the specified algorithm for Generalized Born, will use OBC instead.");
574 if (ir->opts.ngtc > 1)
575 gmx_fatal(FARGS,"OpenMM does not support multiple temperature coupling groups.");
577 if (ir->epc != epcNO)
578 gmx_warning("OpenMM supports only Monte Carlo barostat for pressure coupling.");
580 if (ir->opts.annealing[0])
581 gmx_fatal(FARGS,"OpenMM does not support simulated annealing.");
583 if (top->idef.il[F_CONSTR].nr > 0 && ir->eConstrAlg != econtSHAKE)
584 gmx_warning("OpenMM provides contraints as a combination "
585 "of SHAKE, SETTLE and CCMA. Accuracy is based on the SHAKE tolerance set "
586 "by the \"shake_tol\" option.");
589 gmx_fatal(FARGS,"OpenMM does not support walls.");
591 if (ir->ePull != epullNO)
592 gmx_fatal(FARGS,"OpenMM does not support pulling.");
594 /* check for interaction types */
595 for (i = 0; i < F_EPOT; i++)
597 if (!(i == F_CONSTR ||
600 i == F_UREY_BRADLEY ||
607 i == F_GB12 || /* The GB parameters are hardcoded both in */
608 i == F_GB13 || /* Gromacs and OpenMM */
610 top->idef.il[i].nr > 0)
612 gmx_fatal(FARGS, "OpenMM does not support (some) of the provided interaction "
613 "type(s) (%s) ", interaction_function[i].longname);
617 if (ir->efep != efepNO)
618 gmx_fatal(FARGS,"OpenMM does not support free energy calculations.");
620 if (ir->opts.ngacc > 1)
621 gmx_fatal(FARGS,"OpenMM does not support non-equilibrium MD (accelerated groups).");
623 if (IR_ELEC_FIELD(*ir))
624 gmx_fatal(FARGS,"OpenMM does not support electric fields.");
627 gmx_fatal(FARGS,"OpenMM does not support QMMM calculations.");
629 if (ir->rcoulomb != ir->rvdw)
630 gmx_fatal(FARGS,"OpenMM uses a single cutoff for both Coulomb "
631 "and VdW interactions. Please set rcoulomb equal to rvdw.");
633 if (EEL_FULL(ir->coulombtype))
635 if (ir->ewald_geometry == eewg3DC)
636 gmx_fatal(FARGS,"OpenMM supports only Ewald 3D geometry.");
637 if (ir->epsilon_surface != 0)
638 gmx_fatal(FARGS,"OpenMM does not support dipole correction in Ewald summation.");
641 if (TRICLINIC(state->box))
643 gmx_fatal(FARGS,"OpenMM does not support triclinic unit cells.");
646 /* XXX this is just debugging code to disable the combination rule check */
647 if ( isStringEqNCase(opt->getOptionValue("check-combrule"), "yes") )
649 /* As OpenMM by default uses hardcoded combination rules
650 sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j)
651 we need to check whether the force field params obey this
652 and if not, we can't use this force field so we exit
653 grace-fatal-fully. */
654 real *nbfp = fr->nbfp;
658 fprintf(debug, ">> Atom parameters: <<\n%10s%5s %5s %5s %5s COMB\n",
659 "", "i-j", "j-i", "i-i", "j-j");
661 /* loop over all i-j atom pairs and verify if
662 sigma_ij = sigma_ji = sigma_comb and eps_ij = eps_ji = eps_comb */
663 for (i = 0; i < natoms; i++)
666 c12 = C12(nbfp, natoms, i, i);
667 c6 = C6(nbfp, natoms, i, i);
668 convert_c_12_6(c12, c6, &sigma_ii, &eps_ii);
670 for (j = 0; j < i; j++)
673 c12 = C12(nbfp, natoms, i, j);
674 c6 = C6(nbfp, natoms, i, j);
675 convert_c_12_6(c12, c6, &sigma_ij, &eps_ij);
677 c12 = C12(nbfp, natoms, j, i);
678 c6 = C6(nbfp, natoms, j, i);
679 convert_c_12_6(c12, c6, &sigma_ji, &eps_ji);
681 c12 = C12(nbfp, natoms, j, j);
682 c6 = C6(nbfp, natoms, j, j);
683 convert_c_12_6(c12, c6, &sigma_jj, &eps_jj);
684 /* OpenMM hardcoded combination rules */
685 sigma_comb = COMBRULE_SIGMA(sigma_ii, sigma_jj);
686 eps_comb = COMBRULE_EPS(eps_ii, eps_jj);
690 fprintf(debug, "i=%-3d j=%-3d", i, j);
691 fprintf(debug, "%-11s", "sigma");
692 fprintf(debug, "%5.3f %5.3f %5.3f %5.3f %5.3f\n",
693 sigma_ij, sigma_ji, sigma_ii, sigma_jj, sigma_comb);
694 fprintf(debug, "%11s%-11s", "", "epsilon");
695 fprintf(debug, "%5.3f %5.3f %5.3f %5.3f %5.3f\n",
696 eps_ij, eps_ji, eps_ii, eps_jj, eps_comb);
699 /* check the values against the rule used by omm */
700 if((fabs(eps_ij) > COMBRULE_CHK_TOL &&
701 fabs(eps_ji) > COMBRULE_CHK_TOL) &&
702 (fabs(sigma_comb - sigma_ij) > COMBRULE_CHK_TOL ||
703 fabs(sigma_comb - sigma_ji) > COMBRULE_CHK_TOL ||
704 fabs(eps_comb - eps_ij) > COMBRULE_CHK_TOL ||
705 fabs(eps_comb - eps_ji) > COMBRULE_CHK_TOL ))
708 "The combination rules of the used force-field do not "
709 "match the one supported by OpenMM: "
710 "sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j). "
711 "Switch to a force-field that uses these rules in order to "
712 "simulate this system using OpenMM.\n");
716 if (debug) { fprintf(debug, ">><<\n\n"); }
718 /* if we got here, log that everything is fine */
721 fprintf(debug, ">> The combination rule of the used force matches the one used by OpenMM.\n");
723 fprintf(fplog, "The combination rule of the used force field matches the one used by OpenMM.\n");
725 } /* if (are we checking the combination rules) ... */
730 * \brief Initialize OpenMM, run sanity/consistency checks, and return a pointer to
733 * Various gromacs data structures are passed that contain the parameters, state and
734 * other porperties of the system to simulate. These serve as input for initializing
735 * OpenMM. Besides, a set of misc action are taken:
736 * - OpenMM plugins are loaded;
737 * - platform options in \p platformOptStr are parsed and checked;
738 * - Gromacs parameters are checked for OpenMM support and consistency;
739 * - after the OpenMM is initialized memtest executed in the same GPU context.
741 * \param[in] fplog Gromacs log file handler.
742 * \param[in] platformOptStr Platform option string.
743 * \param[in] ir The Gromacs input parameters, see ::t_inputrec
744 * \param[in] top_global Gromacs system toppology, \see ::gmx_mtop_t
745 * \param[in] top Gromacs node local topology, \see gmx_localtop_t
746 * \param[in] mdatoms Gromacs atom parameters, \see ::t_mdatoms
747 * \param[in] fr \see ::t_forcerec
748 * \param[in] state Gromacs systems state, \see ::t_state
749 * \returns Pointer to a
752 void* openmm_init(FILE *fplog, const char *platformOptStr,
754 gmx_mtop_t *top_global, gmx_localtop_t *top,
755 t_mdatoms *mdatoms, t_forcerec *fr, t_state *state)
758 char warn_buf[STRLEN];
759 static bool hasLoadedPlugins = false;
760 string usedPluginDir;
765 if (!hasLoadedPlugins)
767 vector<string> loadedPlugins;
768 /* Look for OpenMM plugins at various locations (listed in order of priority):
769 - on the path in OPENMM_PLUGIN_DIR environment variable if this is specified
770 - on the path in the OPENMM_PLUGIN_DIR macro that is set by the build script
771 - at the default location assumed by OpenMM
774 char *pluginDir = getenv("OPENMM_PLUGIN_DIR");
776 /* no env var or empty */
777 if (pluginDir != NULL && *pluginDir != '\0')
779 loadedPlugins = Platform::loadPluginsFromDirectory(pluginDir);
780 if (loadedPlugins.size() > 0)
782 hasLoadedPlugins = true;
783 usedPluginDir = pluginDir;
787 gmx_fatal(FARGS, "The directory provided in the OPENMM_PLUGIN_DIR environment variable "
788 "(%s) does not contain valid OpenMM plugins. Check your OpenMM installation!",
793 /* macro set at build time */
794 #ifdef OpenMM_PLUGIN_DIR
795 if (!hasLoadedPlugins)
797 loadedPlugins = Platform::loadPluginsFromDirectory(OPENMM_PLUGIN_DIR);
798 if (loadedPlugins.size() > 0)
800 hasLoadedPlugins = true;
801 usedPluginDir = OPENMM_PLUGIN_DIR;
805 /* default loocation */
806 if (!hasLoadedPlugins)
808 loadedPlugins = Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
809 if (loadedPlugins.size() > 0)
811 hasLoadedPlugins = true;
812 usedPluginDir = Platform::getDefaultPluginsDirectory();
816 /* if there are still no plugins loaded there won't be any */
817 if (!hasLoadedPlugins)
819 gmx_fatal(FARGS, "No OpenMM plugins were found! You can provide the"
820 " plugin directory in the OPENMM_PLUGIN_DIR environment variable.", pluginDir);
823 fprintf(fplog, "\nOpenMM plugins loaded from directory %s:\t", usedPluginDir.c_str());
824 for (int i = 0; i < (int)loadedPlugins.size(); i++)
826 fprintf(fplog, "%s, ", loadedPlugins[i].c_str());
828 fprintf(fplog, "\n");
831 /* parse option string */
832 GmxOpenMMPlatformOptions *opt = new GmxOpenMMPlatformOptions(platformOptStr);
833 devId = atoi(opt->getOptionValue("deviceid").c_str());
840 /* check wheter Gromacs options compatibility with OpenMM */
841 checkGmxOptions(fplog, opt, ir, top, fr, state);
843 // Create the system.
844 const t_idef& idef = top->idef;
845 const int numAtoms = top_global->natoms;
846 const int numConstraints = idef.il[F_CONSTR].nr/3;
847 const int numSettle = idef.il[F_SETTLE].nr/2;
848 const int numBonds = idef.il[F_BONDS].nr/3;
849 const int numUB = idef.il[F_UREY_BRADLEY].nr/4;
850 const int numAngles = idef.il[F_ANGLES].nr/4;
851 const int numPeriodic = idef.il[F_PDIHS].nr/5;
852 const int numPeriodicImproper = idef.il[F_PIDIHS].nr/5;
853 const int numRB = idef.il[F_RBDIHS].nr/5;
854 const int numImproperDih = idef.il[F_IDIHS].nr/5;
855 const int num14 = idef.il[F_LJ14].nr/3;
856 System* sys = new System();
858 sys->addForce(new CMMotionRemover(ir->nstcomm));
860 /* Set bonded force field terms. */
861 const int* bondAtoms = (int*) idef.il[F_BONDS].iatoms;
862 HarmonicBondForce* bondForce = new HarmonicBondForce();
863 sys->addForce(bondForce);
865 for (int i = 0; i < numBonds; ++i)
867 int type = bondAtoms[offset++];
868 int atom1 = bondAtoms[offset++];
869 int atom2 = bondAtoms[offset++];
870 bondForce->addBond(atom1, atom2,
871 idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
874 /* Urey-Bradley includes both the angle and bond potential for 1-3 interactions */
875 const int* ubAtoms = (int*) idef.il[F_UREY_BRADLEY].iatoms;
876 HarmonicBondForce* ubBondForce = new HarmonicBondForce();
877 HarmonicAngleForce* ubAngleForce = new HarmonicAngleForce();
878 sys->addForce(ubBondForce);
879 sys->addForce(ubAngleForce);
881 for (int i = 0; i < numUB; ++i)
883 int type = ubAtoms[offset++];
884 int atom1 = ubAtoms[offset++];
885 int atom2 = ubAtoms[offset++];
886 int atom3 = ubAtoms[offset++];
887 ubBondForce->addBond(atom1, atom3,
888 idef.iparams[type].u_b.r13, idef.iparams[type].u_b.kUB);
889 ubAngleForce->addAngle(atom1, atom2, atom3,
890 idef.iparams[type].u_b.theta*M_PI/180.0, idef.iparams[type].u_b.ktheta);
893 /* Set the angle force field terms */
894 const int* angleAtoms = (int*) idef.il[F_ANGLES].iatoms;
895 HarmonicAngleForce* angleForce = new HarmonicAngleForce();
896 sys->addForce(angleForce);
898 for (int i = 0; i < numAngles; ++i)
900 int type = angleAtoms[offset++];
901 int atom1 = angleAtoms[offset++];
902 int atom2 = angleAtoms[offset++];
903 int atom3 = angleAtoms[offset++];
904 angleForce->addAngle(atom1, atom2, atom3,
905 idef.iparams[type].harmonic.rA*M_PI/180.0, idef.iparams[type].harmonic.krA);
908 /* Set proper dihedral terms */
909 const int* periodicAtoms = (int*) idef.il[F_PDIHS].iatoms;
910 PeriodicTorsionForce* periodicForce = new PeriodicTorsionForce();
911 sys->addForce(periodicForce);
913 for (int i = 0; i < numPeriodic; ++i)
915 int type = periodicAtoms[offset++];
916 int atom1 = periodicAtoms[offset++];
917 int atom2 = periodicAtoms[offset++];
918 int atom3 = periodicAtoms[offset++];
919 int atom4 = periodicAtoms[offset++];
920 periodicForce->addTorsion(atom1, atom2, atom3, atom4,
921 idef.iparams[type].pdihs.mult,
922 idef.iparams[type].pdihs.phiA*M_PI/180.0,
923 idef.iparams[type].pdihs.cpA);
926 /* Set improper dihedral terms that are represented by a periodic function (as in AMBER FF) */
927 const int* periodicImproperAtoms = (int*) idef.il[F_PIDIHS].iatoms;
928 PeriodicTorsionForce* periodicImproperForce = new PeriodicTorsionForce();
929 sys->addForce(periodicImproperForce);
931 for (int i = 0; i < numPeriodicImproper; ++i)
933 int type = periodicImproperAtoms[offset++];
934 int atom1 = periodicImproperAtoms[offset++];
935 int atom2 = periodicImproperAtoms[offset++];
936 int atom3 = periodicImproperAtoms[offset++];
937 int atom4 = periodicImproperAtoms[offset++];
938 periodicImproperForce->addTorsion(atom1, atom2, atom3, atom4,
939 idef.iparams[type].pdihs.mult,
940 idef.iparams[type].pdihs.phiA*M_PI/180.0,
941 idef.iparams[type].pdihs.cpA);
944 /* Ryckaert-Bellemans dihedrals */
945 const int* rbAtoms = (int*) idef.il[F_RBDIHS].iatoms;
946 RBTorsionForce* rbForce = new RBTorsionForce();
947 sys->addForce(rbForce);
949 for (int i = 0; i < numRB; ++i)
951 int type = rbAtoms[offset++];
952 int atom1 = rbAtoms[offset++];
953 int atom2 = rbAtoms[offset++];
954 int atom3 = rbAtoms[offset++];
955 int atom4 = rbAtoms[offset++];
956 rbForce->addTorsion(atom1, atom2, atom3, atom4,
957 idef.iparams[type].rbdihs.rbcA[0], idef.iparams[type].rbdihs.rbcA[1],
958 idef.iparams[type].rbdihs.rbcA[2], idef.iparams[type].rbdihs.rbcA[3],
959 idef.iparams[type].rbdihs.rbcA[4], idef.iparams[type].rbdihs.rbcA[5]);
962 /* Set improper dihedral terms (as in CHARMM FF) */
963 const int* improperDihAtoms = (int*) idef.il[F_IDIHS].iatoms;
964 CustomTorsionForce* improperDihForce = new CustomTorsionForce("2.0*k*asin(sin((theta-theta0)/2))^2");
965 sys->addForce(improperDihForce);
966 improperDihForce->addPerTorsionParameter("k");
967 improperDihForce->addPerTorsionParameter("theta0");
968 vector<double> improperDihParameters(2);
970 for (int i = 0; i < numImproperDih; ++i)
972 int type = improperDihAtoms[offset++];
973 int atom1 = improperDihAtoms[offset++];
974 int atom2 = improperDihAtoms[offset++];
975 int atom3 = improperDihAtoms[offset++];
976 int atom4 = improperDihAtoms[offset++];
977 improperDihParameters[0] = idef.iparams[type].harmonic.krA;
978 improperDihParameters[1] = idef.iparams[type].harmonic.rA*M_PI/180.0;
979 improperDihForce->addTorsion(atom1, atom2, atom3, atom4,
980 improperDihParameters);
983 /* Set nonbonded parameters and masses. */
984 int ntypes = fr->ntype;
985 int* types = mdatoms->typeA;
986 real* nbfp = fr->nbfp;
987 real* charges = mdatoms->chargeA;
988 real* masses = mdatoms->massT;
989 NonbondedForce* nonbondedForce = new NonbondedForce();
990 sys->addForce(nonbondedForce);
995 if (ir->rcoulomb == 0)
997 nonbondedForce->setNonbondedMethod(NonbondedForce::NoCutoff);
1001 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
1005 switch (ir->coulombtype)
1011 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
1015 nonbondedForce->setNonbondedMethod(NonbondedForce::Ewald);
1019 nonbondedForce->setNonbondedMethod(NonbondedForce::PME);
1023 gmx_fatal(FARGS,"Internal error: you should not see this message, it means that the"
1024 "electrosatics option check failed. Please report this error!");
1026 sys->setDefaultPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0),
1027 Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));
1028 nonbondedForce->setCutoffDistance(ir->rcoulomb);
1032 gmx_fatal(FARGS,"OpenMM supports only full periodic boundary conditions "
1033 "(pbc = xyz), or none (pbc = no).");
1037 /* Fix for PME and Ewald error tolerance
1039 * OpenMM uses approximate formulas to calculate the Ewald parameter:
1040 * alpha = (1.0/cutoff)*sqrt(-log(2.0*tolerlance));
1041 * and the grid spacing for PME:
1042 * gridX = ceil(2*alpha*box[0][0]/3*(pow(tol, 0.2)))
1043 * gridY = ceil(2*alpha*box[1][1]/3*(pow(tol, 0.2)));
1044 * gridZ = ceil(2*alpha*box[2][2]/3*(pow(tol, 0.2)));
1047 * If the default ewald_rtol=1e-5 is used we silently adjust the value to the
1048 * OpenMM default of 5e-4 otherwise a warning is issued about the action taken.
1051 double corr_ewald_rtol = 50.0 * ir->ewald_rtol;
1052 if ((ir->ePBC == epbcXYZ) &&
1053 (ir->coulombtype == eelEWALD || ir->coulombtype == eelPME))
1057 fprintf(debug, ">> ewald_rtol = %e (corrected = %e) \n",
1058 ir->ewald_rtol, corr_ewald_rtol);
1061 if (fabs(ir->ewald_rtol - 1e-5) > 1e-10)
1063 gmx_warning("OpenMM uses the ewald_rtol parameter with approximate formulas "
1064 "to calculate the alpha and grid spacing parameters of the Ewald "
1065 "and PME methods. This tolerance need to be corrected in order to get "
1066 "settings close to the ones used in GROMACS. Although the internal correction "
1067 "should work for any reasonable value of ewald_rtol, using values other than "
1068 "the default 1e-5 might cause incorrect behavior.");
1070 if (corr_ewald_rtol > 1)
1072 gmx_fatal(FARGS, "The ewald_rtol accuracy term is >1 after the "
1073 "adjustment for OpenMM (%e)", corr_ewald_rtol);
1076 nonbondedForce->setEwaldErrorTolerance(corr_ewald_rtol);
1079 for (int i = 0; i < numAtoms; ++i)
1081 double c12 = nbfp[types[i]*2*ntypes+types[i]*2+1];
1082 double c6 = nbfp[types[i]*2*ntypes+types[i]*2];
1083 double sigma=0.0, epsilon=0.0;
1084 convert_c_12_6(c12, c6, &sigma, &epsilon);
1085 nonbondedForce->addParticle(charges[i], sigma, epsilon);
1086 sys->addParticle(masses[i]);
1089 // Build a table of all exclusions.
1090 vector<set<int> > exclusions(numAtoms);
1091 for (int i = 0; i < numAtoms; i++)
1093 int start = top->excls.index[i];
1094 int end = top->excls.index[i+1];
1095 for (int j = start; j < end; j++)
1096 exclusions[i].insert(top->excls.a[j]);
1099 // Record the 1-4 interactions, and remove them from the list of exclusions.
1100 const int* nb14Atoms = (int*) idef.il[F_LJ14].iatoms;
1102 for (int i = 0; i < num14; ++i)
1104 int type = nb14Atoms[offset++];
1105 int atom1 = nb14Atoms[offset++];
1106 int atom2 = nb14Atoms[offset++];
1107 double sigma=0, epsilon=0;
1108 convert_c_12_6(idef.iparams[type].lj14.c12A,
1109 idef.iparams[type].lj14.c6A,
1111 nonbondedForce->addException(atom1, atom2,
1112 fr->fudgeQQ*charges[atom1]*charges[atom2], sigma, epsilon);
1113 exclusions[atom1].erase(atom2);
1114 exclusions[atom2].erase(atom1);
1117 // Record exclusions.
1118 for (int i = 0; i < numAtoms; i++)
1120 for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
1124 nonbondedForce->addException(i, *iter, 0.0, 1.0, 0.0);
1129 // Add GBSA if needed.
1130 if (ir->implicit_solvent == eisGBSA)
1132 gmx_warning("The OBC scale factors alpha, beta and gamma are hardcoded in OpenMM with the default Gromacs values.");
1133 t_atoms atoms = gmx_mtop_global_atoms(top_global);
1134 GBSAOBCForce* gbsa = new GBSAOBCForce();
1136 sys->addForce(gbsa);
1137 gbsa->setSoluteDielectric(ir->epsilon_r);
1138 gbsa->setSolventDielectric(ir->gb_epsilon_solvent);
1139 gbsa->setCutoffDistance(nonbondedForce->getCutoffDistance());
1140 if (nonbondedForce->getNonbondedMethod() == NonbondedForce::NoCutoff)
1141 gbsa->setNonbondedMethod(GBSAOBCForce::NoCutoff);
1142 else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
1143 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
1144 else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffPeriodic)
1145 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
1147 gmx_fatal(FARGS,"OpenMM supports only Reaction-Field electrostatics with OBC/GBSA.");
1149 for (int i = 0; i < numAtoms; ++i)
1151 gbsa->addParticle(charges[i],
1152 top_global->atomtypes.gb_radius[atoms.atom[i].type],
1153 top_global->atomtypes.S_hct[atoms.atom[i].type]);
1155 free_t_atoms(&atoms, FALSE);
1159 const int* constraintAtoms = (int*) idef.il[F_CONSTR].iatoms;
1161 for (int i = 0; i < numConstraints; ++i)
1163 int type = constraintAtoms[offset++];
1164 int atom1 = constraintAtoms[offset++];
1165 int atom2 = constraintAtoms[offset++];
1166 sys->addConstraint(atom1, atom2, idef.iparams[type].constr.dA);
1168 const int* settleAtoms = (int*) idef.il[F_SETTLE].iatoms;
1170 for (int i = 0; i < numSettle; ++i)
1172 int type = settleAtoms[offset++];
1173 int oxygen = settleAtoms[offset++];
1174 sys->addConstraint(oxygen, oxygen+1, idef.iparams[type].settle.doh);
1175 sys->addConstraint(oxygen, oxygen+2, idef.iparams[type].settle.doh);
1176 sys->addConstraint(oxygen+1, oxygen+2, idef.iparams[type].settle.dhh);
1179 // Create an integrator for simulating the system.
1180 double friction = (ir->opts.tau_t[0] == 0.0 ? 0.0 : 1.0/ir->opts.tau_t[0]);
1184 integ = new BrownianIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1185 static_cast<BrownianIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1187 else if (EI_SD(ir->eI))
1189 integ = new LangevinIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1190 static_cast<LangevinIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1194 integ = new VerletIntegrator(ir->delta_t);
1195 if ( ir->etc != etcNO)
1197 AndersenThermostat* thermostat = new AndersenThermostat(ir->opts.ref_t[0], friction);
1198 sys->addForce(thermostat);
1202 // Add pressure coupling
1203 if (ir->epc != epcNO)
1205 // convert gromacs pressure tensor to a scalar
1206 double pressure = (ir->ref_p[0][0] + ir->ref_p[1][1] + ir->ref_p[2][2]) / 3.0;
1207 int frequency = int(ir->tau_p / ir->delta_t); // update frequency in time steps
1208 if (frequency < 1) frequency = 1;
1209 double temperature = ir->opts.ref_t[0]; // in kelvin
1210 sys->addForce(new MonteCarloBarostat(pressure, temperature, frequency));
1213 integ->setConstraintTolerance(ir->shake_tol);
1215 // Create a context and initialize it.
1216 Context* context = NULL;
1219 OpenMM could automatically select the "best" GPU, however we're not't
1220 going to let it do that for now, as the current algorithm is very rudimentary
1221 and we anyway support only CUDA.
1222 if (platformOptStr == NULL || platformOptStr == "")
1224 context = new Context(*sys, *integ);
1229 /* which platform should we use */
1230 for (int i = 0; i < (int)Platform::getNumPlatforms() && context == NULL; i++)
1232 if (isStringEqNCase(opt->getOptionValue("platform"), Platform::getPlatform(i).getName()))
1234 Platform& platform = Platform::getPlatform(i);
1235 // set standard properties
1236 platform.setPropertyDefaultValue("CudaDevice", opt->getOptionValue("deviceid"));
1237 // TODO add extra properties
1238 context = new Context(*sys, *integ, platform);
1241 if (context == NULL)
1243 gmx_fatal(FARGS, "The requested platform \"%s\" could not be found.",
1244 opt->getOptionValue("platform").c_str());
1248 Platform& platform = context->getPlatform();
1249 fprintf(fplog, "Gromacs will use the OpenMM platform: %s\n", platform.getName().c_str());
1251 const vector<string>& properties = platform.getPropertyNames();
1254 for (int i = 0; i < (int)properties.size(); i++)
1256 fprintf(debug, ">> %s: %s\n", properties[i].c_str(),
1257 platform.getPropertyValue(*context, properties[i]).c_str());
1262 if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1265 if (!from_string<int>(tmp, platform.getPropertyValue(*context, "CudaDevice"), std::dec))
1267 gmx_fatal(FARGS, "Internal error: couldn't determine the device selected by OpenMM");
1271 /* For now this is just to double-check if OpenMM selected the GPU we wanted,
1272 but when we'll let OpenMM select the GPU automatically, it will query the devideId.
1276 gmx_fatal(FARGS, "Internal error: OpenMM is using device #%d"
1277 "while initialized for device #%d", tmp, devId);
1279 cout << ">>>>> OpenMM devId=" << tmp << endl;
1281 /* check GPU compatibility */
1282 char gpuname[STRLEN];
1283 devId = atoi(opt->getOptionValue("deviceid").c_str());
1284 if (!is_supported_cuda_gpu(-1, gpuname))
1286 if (!gmx_strcasecmp(opt->getOptionValue("force-device").c_str(), "yes"))
1288 sprintf(warn_buf, "Non-supported GPU selected (#%d, %s), forced continuing."
1289 "Note, that the simulation can be slow or it migth even crash.",
1291 fprintf(fplog, "%s\n", warn_buf);
1292 gmx_warning(warn_buf);
1296 gmx_fatal(FARGS, "The selected GPU (#%d, %s) is not supported by Gromacs! "
1297 "Most probably you have a low-end GPU which would not perform well, "
1298 "or new hardware that has not been tested with the current release. "
1299 "If you still want to try using the device, use the force-device=yes option.",
1305 fprintf(fplog, "Gromacs will run on the GPU #%d (%s).\n", devId, gpuname);
1310 if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1312 /* pre-simulation memtest */
1313 runMemtest(fplog, -1, "Pre", opt);
1316 vector<Vec3> pos(numAtoms);
1317 vector<Vec3> vel(numAtoms);
1318 for (int i = 0; i < numAtoms; ++i)
1320 pos[i] = Vec3(state->x[i][0], state->x[i][1], state->x[i][2]);
1321 vel[i] = Vec3(state->v[i][0], state->v[i][1], state->v[i][2]);
1323 context->setPositions(pos);
1324 context->setVelocities(vel);
1326 // Return a structure containing the system, integrator, and context.
1327 OpenMMData* data = new OpenMMData();
1329 data->integrator = integ;
1330 data->context = context;
1331 data->removeCM = (ir->nstcomm > 0);
1332 data->platformOpt = opt;
1335 catch (std::exception& e)
1337 gmx_fatal(FARGS, "OpenMM exception caught while initializating: %s", e.what());
1339 return NULL; /* just to avoid warnings */
1343 * \brief Integrate one step.
1345 * \param[in] data OpenMMData object created by openmm_init().
1347 void openmm_take_one_step(void* data)
1349 // static int step = 0; printf("----> taking step #%d\n", step++);
1352 static_cast<OpenMMData*>(data)->integrator->step(1);
1354 catch (std::exception& e)
1356 gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1361 * \brief Integrate n steps.
1363 * \param[in] data OpenMMData object created by openmm_init().
1365 void openmm_take_steps(void* data, int nstep)
1369 static_cast<OpenMMData*>(data)->integrator->step(nstep);
1371 catch (std::exception& e)
1373 gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1378 * \brief Clean up the data structures cretead for OpenMM.
1380 * \param[in] log Log file pointer.
1381 * \param[in] data OpenMMData object created by openmm_init().
1383 void openmm_cleanup(FILE* fplog, void* data)
1385 OpenMMData* d = static_cast<OpenMMData*>(data);
1387 if (isStringEqNCase(d->platformOpt->getOptionValue("platform"), "CUDA"))
1389 /* post-simulation memtest */
1390 runMemtest(fplog, -1, "Post", d->platformOpt);
1393 delete d->integrator;
1395 delete d->platformOpt;
1400 * \brief Copy the current state information from OpenMM into the Gromacs data structures.
1402 * This function results in the requested proprties to be copied from the
1403 * GPU to host. As this represents a bottleneck, the frequency of pulling data
1404 * should be minimized.
1406 * \param[in] data OpenMMData object created by openmm_init().
1407 * \param[out] time Simulation time for which the state was created.
1408 * \param[out] state State of the system: coordinates and velocities.
1409 * \param[out] f Forces.
1410 * \param[out] enerd Energies.
1411 * \param[in] includePos True if coordinates are requested.
1412 * \param[in] includeVel True if velocities are requested.
1413 * \param[in] includeForce True if forces are requested.
1414 * \param[in] includeEnergy True if energies are requested.
1416 void openmm_copy_state(void *data,
1417 t_state *state, double *time,
1418 rvec f[], gmx_enerdata_t *enerd,
1419 bool includePos, bool includeVel, bool includeForce, bool includeEnergy)
1423 types += State::Positions;
1425 types += State::Velocities;
1427 types += State::Forces;
1429 types += State::Energy;
1434 State currentState = static_cast<OpenMMData*>(data)->context->getState(types);
1435 int numAtoms = static_cast<OpenMMData*>(data)->system->getNumParticles();
1438 for (int i = 0; i < numAtoms; i++)
1440 Vec3 x = currentState.getPositions()[i];
1441 state->x[i][0] = x[0];
1442 state->x[i][1] = x[1];
1443 state->x[i][2] = x[2];
1448 for (int i = 0; i < numAtoms; i++)
1450 Vec3 v = currentState.getVelocities()[i];
1451 state->v[i][0] = v[0];
1452 state->v[i][1] = v[1];
1453 state->v[i][2] = v[2];
1458 for (int i = 0; i < numAtoms; i++)
1460 Vec3 force = currentState.getForces()[i];
1468 int numConstraints = static_cast<OpenMMData*>(data)->system->getNumConstraints();
1469 int dof = 3*numAtoms-numConstraints;
1470 if (static_cast<OpenMMData*>(data)->removeCM)
1472 enerd->term[F_EPOT] = currentState.getPotentialEnergy();
1473 enerd->term[F_EKIN] = currentState.getKineticEnergy();
1474 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1475 enerd->term[F_TEMP] = 2.0*enerd->term[F_EKIN]/dof/BOLTZ;
1477 *time = currentState.getTime();
1479 catch (std::exception& e)
1481 gmx_fatal(FARGS, "OpenMM exception caught while retrieving state information: %s", e.what());