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 1 // 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]
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);
873 // Urey-Bradley includes both the angle and bond potential for 1-3 interactions
874 const int* ubAtoms = (int*) idef.il[F_UREY_BRADLEY].iatoms;
875 HarmonicBondForce* ubBondForce = new HarmonicBondForce();
876 HarmonicAngleForce* ubAngleForce = new HarmonicAngleForce();
877 sys->addForce(ubBondForce);
878 sys->addForce(ubAngleForce);
880 for (int i = 0; i < numUB; ++i)
882 int type = ubAtoms[offset++];
883 int atom1 = ubAtoms[offset++];
884 int atom2 = ubAtoms[offset++];
885 int atom3 = ubAtoms[offset++];
886 ubBondForce->addBond(atom1, atom3,
887 idef.iparams[type].u_b.r13, idef.iparams[type].u_b.kUB);
888 ubAngleForce->addAngle(atom1, atom2, atom3,
889 idef.iparams[type].u_b.theta*M_PI/180.0, idef.iparams[type].u_b.ktheta);
891 const int* angleAtoms = (int*) idef.il[F_ANGLES].iatoms;
892 HarmonicAngleForce* angleForce = new HarmonicAngleForce();
893 sys->addForce(angleForce);
895 for (int i = 0; i < numAngles; ++i)
897 int type = angleAtoms[offset++];
898 int atom1 = angleAtoms[offset++];
899 int atom2 = angleAtoms[offset++];
900 int atom3 = angleAtoms[offset++];
901 angleForce->addAngle(atom1, atom2, atom3,
902 idef.iparams[type].harmonic.rA*M_PI/180.0, idef.iparams[type].harmonic.krA);
904 const int* periodicAtoms = (int*) idef.il[F_PDIHS].iatoms;
905 PeriodicTorsionForce* periodicForce = new PeriodicTorsionForce();
906 sys->addForce(periodicForce);
908 for (int i = 0; i < numPeriodic; ++i)
910 int type = periodicAtoms[offset++];
911 int atom1 = periodicAtoms[offset++];
912 int atom2 = periodicAtoms[offset++];
913 int atom3 = periodicAtoms[offset++];
914 int atom4 = periodicAtoms[offset++];
915 periodicForce->addTorsion(atom1, atom2, atom3, atom4,
916 idef.iparams[type].pdihs.mult,
917 idef.iparams[type].pdihs.phiA*M_PI/180.0,
918 idef.iparams[type].pdihs.cpA);
920 const int* periodicImproperAtoms = (int*) idef.il[F_PIDIHS].iatoms;
921 PeriodicTorsionForce* periodicImproperForce = new PeriodicTorsionForce();
922 sys->addForce(periodicImproperForce);
924 for (int i = 0; i < numPeriodicImproper; ++i)
926 int type = periodicImproperAtoms[offset++];
927 int atom1 = periodicImproperAtoms[offset++];
928 int atom2 = periodicImproperAtoms[offset++];
929 int atom3 = periodicImproperAtoms[offset++];
930 int atom4 = periodicImproperAtoms[offset++];
931 periodicImproperForce->addTorsion(atom1, atom2, atom3, atom4,
932 idef.iparams[type].pdihs.mult,
933 idef.iparams[type].pdihs.phiA*M_PI/180.0,
934 idef.iparams[type].pdihs.cpA);
936 const int* rbAtoms = (int*) idef.il[F_RBDIHS].iatoms;
937 RBTorsionForce* rbForce = new RBTorsionForce();
938 sys->addForce(rbForce);
940 for (int i = 0; i < numRB; ++i)
942 int type = rbAtoms[offset++];
943 int atom1 = rbAtoms[offset++];
944 int atom2 = rbAtoms[offset++];
945 int atom3 = rbAtoms[offset++];
946 int atom4 = rbAtoms[offset++];
947 rbForce->addTorsion(atom1, atom2, atom3, atom4,
948 idef.iparams[type].rbdihs.rbcA[0], idef.iparams[type].rbdihs.rbcA[1],
949 idef.iparams[type].rbdihs.rbcA[2], idef.iparams[type].rbdihs.rbcA[3],
950 idef.iparams[type].rbdihs.rbcA[4], idef.iparams[type].rbdihs.rbcA[5]);
953 const int* improperDihAtoms = (int*) idef.il[F_IDIHS].iatoms;
954 CustomTorsionForce* improperDihForce = new CustomTorsionForce("2.0*k*asin(sin((theta-theta0)/2))^2");
955 sys->addForce(improperDihForce);
956 improperDihForce->addPerTorsionParameter("k");
957 improperDihForce->addPerTorsionParameter("theta0");
958 vector<double> improperDihParameters(2);
960 for (int i = 0; i < numImproperDih; ++i)
962 int type = improperDihAtoms[offset++];
963 int atom1 = improperDihAtoms[offset++];
964 int atom2 = improperDihAtoms[offset++];
965 int atom3 = improperDihAtoms[offset++];
966 int atom4 = improperDihAtoms[offset++];
967 improperDihParameters[0] = idef.iparams[type].harmonic.krA;
968 improperDihParameters[1] = idef.iparams[type].harmonic.rA*M_PI/180.0;
969 improperDihForce->addTorsion(atom1, atom2, atom3, atom4,
970 improperDihParameters);
973 // Set nonbonded parameters and masses.
974 int ntypes = fr->ntype;
975 int* types = mdatoms->typeA;
976 real* nbfp = fr->nbfp;
977 real* charges = mdatoms->chargeA;
978 real* masses = mdatoms->massT;
979 NonbondedForce* nonbondedForce = new NonbondedForce();
980 sys->addForce(nonbondedForce);
985 if (ir->rcoulomb == 0)
987 nonbondedForce->setNonbondedMethod(NonbondedForce::NoCutoff);
991 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
995 switch (ir->coulombtype)
1001 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
1005 nonbondedForce->setNonbondedMethod(NonbondedForce::Ewald);
1009 nonbondedForce->setNonbondedMethod(NonbondedForce::PME);
1013 gmx_fatal(FARGS,"Internal error: you should not see this message, it that the"
1014 "electrosatics option check failed. Please report this error!");
1016 sys->setDefaultPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0),
1017 Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));
1018 nonbondedForce->setCutoffDistance(ir->rcoulomb);
1022 gmx_fatal(FARGS,"OpenMM supports only full periodic boundary conditions "
1023 "(pbc = xyz), or none (pbc = no).");
1027 /* Fix for PME and Ewald error tolerance
1029 * OpenMM uses approximate formulas to calculate the Ewald parameter:
1030 * alpha = (1.0/cutoff)*sqrt(-log(2.0*tolerlance));
1031 * and the grid spacing for PME:
1032 * gridX = ceil(2*alpha*box[0][0]/3*(pow(tol, 0.2)))
1033 * gridY = ceil(2*alpha*box[1][1]/3*(pow(tol, 0.2)));
1034 * gridZ = ceil(2*alpha*box[2][2]/3*(pow(tol, 0.2)));
1037 * If the default ewald_rtol=1e-5 is used we silently adjust the value to the
1038 * OpenMM default of 5e-4 otherwise a warning is issued about the action taken.
1041 double corr_ewald_rtol = 50.0 * ir->ewald_rtol;
1042 if ((ir->ePBC == epbcXYZ) &&
1043 (ir->coulombtype == eelEWALD || ir->coulombtype == eelPME))
1047 fprintf(debug, ">> ewald_rtol = %e (corrected = %e) \n",
1048 ir->ewald_rtol, corr_ewald_rtol);
1051 if (fabs(ir->ewald_rtol - 1e-5) > 1e-10)
1053 gmx_warning("OpenMM uses the ewald_rtol parameter with approximate formulas "
1054 "to calculate the alpha and grid spacing parameters of the Ewald "
1055 "and PME methods. This tolerance need to be corrected in order to get "
1056 "settings close to the ones used in GROMACS. Although the internal correction "
1057 "should work for any reasonable value of ewald_rtol, using values other than "
1058 "the default 1e-5 might cause incorrect behavior.");
1060 if (corr_ewald_rtol > 1)
1062 gmx_fatal(FARGS, "The ewald_rtol accuracy term is >1 after the "
1063 "adjustment for OpenMM (%e)", corr_ewald_rtol);
1066 nonbondedForce->setEwaldErrorTolerance(corr_ewald_rtol);
1069 for (int i = 0; i < numAtoms; ++i)
1071 double c12 = nbfp[types[i]*2*ntypes+types[i]*2+1];
1072 double c6 = nbfp[types[i]*2*ntypes+types[i]*2];
1073 double sigma=0.0, epsilon=0.0;
1074 convert_c_12_6(c12, c6, &sigma, &epsilon);
1075 nonbondedForce->addParticle(charges[i], sigma, epsilon);
1076 sys->addParticle(masses[i]);
1079 // Build a table of all exclusions.
1080 vector<set<int> > exclusions(numAtoms);
1081 for (int i = 0; i < numAtoms; i++)
1083 int start = top->excls.index[i];
1084 int end = top->excls.index[i+1];
1085 for (int j = start; j < end; j++)
1086 exclusions[i].insert(top->excls.a[j]);
1089 // Record the 1-4 interactions, and remove them from the list of exclusions.
1090 const int* nb14Atoms = (int*) idef.il[F_LJ14].iatoms;
1092 for (int i = 0; i < num14; ++i)
1094 int type = nb14Atoms[offset++];
1095 int atom1 = nb14Atoms[offset++];
1096 int atom2 = nb14Atoms[offset++];
1097 double sigma=0, epsilon=0;
1098 convert_c_12_6(idef.iparams[type].lj14.c12A,
1099 idef.iparams[type].lj14.c6A,
1101 nonbondedForce->addException(atom1, atom2,
1102 fr->fudgeQQ*charges[atom1]*charges[atom2], sigma, epsilon);
1103 exclusions[atom1].erase(atom2);
1104 exclusions[atom2].erase(atom1);
1107 // Record exclusions.
1108 for (int i = 0; i < numAtoms; i++)
1110 for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
1114 nonbondedForce->addException(i, *iter, 0.0, 1.0, 0.0);
1119 // Add GBSA if needed.
1120 if (ir->implicit_solvent == eisGBSA)
1122 gmx_warning("The OBC scale factors alpha, beta and gamma are hardcoded in OpenMM with the default Gromacs values.");
1123 t_atoms atoms = gmx_mtop_global_atoms(top_global);
1124 GBSAOBCForce* gbsa = new GBSAOBCForce();
1126 sys->addForce(gbsa);
1127 gbsa->setSoluteDielectric(ir->epsilon_r);
1128 gbsa->setSolventDielectric(ir->gb_epsilon_solvent);
1129 gbsa->setCutoffDistance(nonbondedForce->getCutoffDistance());
1130 if (nonbondedForce->getNonbondedMethod() == NonbondedForce::NoCutoff)
1131 gbsa->setNonbondedMethod(GBSAOBCForce::NoCutoff);
1132 else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
1133 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
1134 else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffPeriodic)
1135 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
1137 gmx_fatal(FARGS,"OpenMM supports only Reaction-Field electrostatics with OBC/GBSA.");
1139 for (int i = 0; i < numAtoms; ++i)
1141 gbsa->addParticle(charges[i],
1142 top_global->atomtypes.gb_radius[atoms.atom[i].type],
1143 top_global->atomtypes.S_hct[atoms.atom[i].type]);
1145 free_t_atoms(&atoms, FALSE);
1149 const int* constraintAtoms = (int*) idef.il[F_CONSTR].iatoms;
1151 for (int i = 0; i < numConstraints; ++i)
1153 int type = constraintAtoms[offset++];
1154 int atom1 = constraintAtoms[offset++];
1155 int atom2 = constraintAtoms[offset++];
1156 sys->addConstraint(atom1, atom2, idef.iparams[type].constr.dA);
1158 const int* settleAtoms = (int*) idef.il[F_SETTLE].iatoms;
1160 for (int i = 0; i < numSettle; ++i)
1162 int type = settleAtoms[offset++];
1163 int oxygen = settleAtoms[offset++];
1164 sys->addConstraint(oxygen, oxygen+1, idef.iparams[type].settle.doh);
1165 sys->addConstraint(oxygen, oxygen+2, idef.iparams[type].settle.doh);
1166 sys->addConstraint(oxygen+1, oxygen+2, idef.iparams[type].settle.dhh);
1169 // Create an integrator for simulating the system.
1170 double friction = (ir->opts.tau_t[0] == 0.0 ? 0.0 : 1.0/ir->opts.tau_t[0]);
1174 integ = new BrownianIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1175 static_cast<BrownianIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1177 else if (EI_SD(ir->eI))
1179 integ = new LangevinIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1180 static_cast<LangevinIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1184 integ = new VerletIntegrator(ir->delta_t);
1185 if ( ir->etc != etcNO)
1187 AndersenThermostat* thermostat = new AndersenThermostat(ir->opts.ref_t[0], friction);
1188 sys->addForce(thermostat);
1192 // Add pressure coupling
1193 if (ir->epc != epcNO)
1195 // convert gromacs pressure tensor to a scalar
1196 double pressure = (ir->ref_p[0][0] + ir->ref_p[1][1] + ir->ref_p[2][2]) / 3.0;
1197 int frequency = int(ir->tau_p / ir->delta_t); // update frequency in time steps
1198 if (frequency < 1) frequency = 1;
1199 double temperature = ir->opts.ref_t[0]; // in kelvin
1200 sys->addForce(new MonteCarloBarostat(pressure, temperature, frequency));
1203 integ->setConstraintTolerance(ir->shake_tol);
1205 // Create a context and initialize it.
1206 Context* context = NULL;
1209 OpenMM could automatically select the "best" GPU, however we're not't
1210 going to let it do that for now, as the current algorithm is very rudimentary
1211 and we anyway support only CUDA.
1212 if (platformOptStr == NULL || platformOptStr == "")
1214 context = new Context(*sys, *integ);
1219 /* which platform should we use */
1220 for (int i = 0; i < (int)Platform::getNumPlatforms() && context == NULL; i++)
1222 if (isStringEqNCase(opt->getOptionValue("platform"), Platform::getPlatform(i).getName()))
1224 Platform& platform = Platform::getPlatform(i);
1225 // set standard properties
1226 platform.setPropertyDefaultValue("CudaDevice", opt->getOptionValue("deviceid"));
1227 // TODO add extra properties
1228 context = new Context(*sys, *integ, platform);
1231 if (context == NULL)
1233 gmx_fatal(FARGS, "The requested platform \"%s\" could not be found.",
1234 opt->getOptionValue("platform").c_str());
1238 Platform& platform = context->getPlatform();
1239 fprintf(fplog, "Gromacs will use the OpenMM platform: %s\n", platform.getName().c_str());
1241 const vector<string>& properties = platform.getPropertyNames();
1244 for (int i = 0; i < (int)properties.size(); i++)
1246 fprintf(debug, ">> %s: %s\n", properties[i].c_str(),
1247 platform.getPropertyValue(*context, properties[i]).c_str());
1252 if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1255 if (!from_string<int>(tmp, platform.getPropertyValue(*context, "CudaDevice"), std::dec))
1257 gmx_fatal(FARGS, "Internal error: couldn't determine the device selected by OpenMM");
1261 /* For now this is just to double-check if OpenMM selected the GPU we wanted,
1262 but when we'll let OpenMM select the GPU automatically, it will query the devideId.
1266 gmx_fatal(FARGS, "Internal error: OpenMM is using device #%d"
1267 "while initialized for device #%d", tmp, devId);
1269 cout << ">>>>> OpenMM devId=" << tmp << endl;
1271 /* check GPU compatibility */
1272 char gpuname[STRLEN];
1273 devId = atoi(opt->getOptionValue("deviceid").c_str());
1274 if (!is_supported_cuda_gpu(-1, gpuname))
1276 if (!gmx_strcasecmp(opt->getOptionValue("force-device").c_str(), "yes"))
1278 sprintf(warn_buf, "Non-supported GPU selected (#%d, %s), forced continuing."
1279 "Note, that the simulation can be slow or it migth even crash.",
1281 fprintf(fplog, "%s\n", warn_buf);
1282 gmx_warning(warn_buf);
1286 gmx_fatal(FARGS, "The selected GPU (#%d, %s) is not supported by Gromacs! "
1287 "Most probably you have a low-end GPU which would not perform well, "
1288 "or new hardware that has not been tested with the current release. "
1289 "If you still want to try using the device, use the force-device=yes option.",
1295 fprintf(fplog, "Gromacs will run on the GPU #%d (%s).\n", devId, gpuname);
1300 if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1302 /* pre-simulation memtest */
1303 runMemtest(fplog, -1, "Pre", opt);
1306 vector<Vec3> pos(numAtoms);
1307 vector<Vec3> vel(numAtoms);
1308 for (int i = 0; i < numAtoms; ++i)
1310 pos[i] = Vec3(state->x[i][0], state->x[i][1], state->x[i][2]);
1311 vel[i] = Vec3(state->v[i][0], state->v[i][1], state->v[i][2]);
1313 context->setPositions(pos);
1314 context->setVelocities(vel);
1316 // Return a structure containing the system, integrator, and context.
1317 OpenMMData* data = new OpenMMData();
1319 data->integrator = integ;
1320 data->context = context;
1321 data->removeCM = (ir->nstcomm > 0);
1322 data->platformOpt = opt;
1325 catch (std::exception& e)
1327 gmx_fatal(FARGS, "OpenMM exception caught while initializating: %s", e.what());
1329 return NULL; /* just to avoid warnings */
1333 * \brief Integrate one step.
1335 * \param[in] data OpenMMData object created by openmm_init().
1337 void openmm_take_one_step(void* data)
1339 // static int step = 0; printf("----> taking step #%d\n", step++);
1342 static_cast<OpenMMData*>(data)->integrator->step(1);
1344 catch (std::exception& e)
1346 gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1351 * \brief Integrate n steps.
1353 * \param[in] data OpenMMData object created by openmm_init().
1355 void openmm_take_steps(void* data, int nstep)
1359 static_cast<OpenMMData*>(data)->integrator->step(nstep);
1361 catch (std::exception& e)
1363 gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1368 * \brief Clean up the data structures cretead for OpenMM.
1370 * \param[in] log Log file pointer.
1371 * \param[in] data OpenMMData object created by openmm_init().
1373 void openmm_cleanup(FILE* fplog, void* data)
1375 OpenMMData* d = static_cast<OpenMMData*>(data);
1377 if (isStringEqNCase(d->platformOpt->getOptionValue("platform"), "CUDA"))
1379 /* post-simulation memtest */
1380 runMemtest(fplog, -1, "Post", d->platformOpt);
1383 delete d->integrator;
1385 delete d->platformOpt;
1390 * \brief Copy the current state information from OpenMM into the Gromacs data structures.
1392 * This function results in the requested proprties to be copied from the
1393 * GPU to host. As this represents a bottleneck, the frequency of pulling data
1394 * should be minimized.
1396 * \param[in] data OpenMMData object created by openmm_init().
1397 * \param[out] time Simulation time for which the state was created.
1398 * \param[out] state State of the system: coordinates and velocities.
1399 * \param[out] f Forces.
1400 * \param[out] enerd Energies.
1401 * \param[in] includePos True if coordinates are requested.
1402 * \param[in] includeVel True if velocities are requested.
1403 * \param[in] includeForce True if forces are requested.
1404 * \param[in] includeEnergy True if energies are requested.
1406 void openmm_copy_state(void *data,
1407 t_state *state, double *time,
1408 rvec f[], gmx_enerdata_t *enerd,
1409 bool includePos, bool includeVel, bool includeForce, bool includeEnergy)
1413 types += State::Positions;
1415 types += State::Velocities;
1417 types += State::Forces;
1419 types += State::Energy;
1424 State currentState = static_cast<OpenMMData*>(data)->context->getState(types);
1425 int numAtoms = static_cast<OpenMMData*>(data)->system->getNumParticles();
1428 for (int i = 0; i < numAtoms; i++)
1430 Vec3 x = currentState.getPositions()[i];
1431 state->x[i][0] = x[0];
1432 state->x[i][1] = x[1];
1433 state->x[i][2] = x[2];
1438 for (int i = 0; i < numAtoms; i++)
1440 Vec3 v = currentState.getVelocities()[i];
1441 state->v[i][0] = v[0];
1442 state->v[i][1] = v[1];
1443 state->v[i][2] = v[2];
1448 for (int i = 0; i < numAtoms; i++)
1450 Vec3 force = currentState.getForces()[i];
1458 int numConstraints = static_cast<OpenMMData*>(data)->system->getNumConstraints();
1459 int dof = 3*numAtoms-numConstraints;
1460 if (static_cast<OpenMMData*>(data)->removeCM)
1462 enerd->term[F_EPOT] = currentState.getPotentialEnergy();
1463 enerd->term[F_EKIN] = currentState.getKineticEnergy();
1464 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1465 enerd->term[F_TEMP] = 2.0*enerd->term[F_EKIN]/dof/BOLTZ;
1467 *time = currentState.getTime();
1469 catch (std::exception& e)
1471 gmx_fatal(FARGS, "OpenMM exception caught while retrieving state information: %s", e.what());