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, 2013, 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 Lesser 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
46 #include <types/simple.h>
61 #include "gmx_fatal.h"
66 #include "openmm_gpu_utils.h"
67 #include "gpu_utils.h"
68 #include "mtop_util.h"
70 #include "openmm_wrapper.h"
72 using namespace OpenMM;
75 #define MEM_ERR_MSG(str) \
76 "The %s-simulation GPU memory test detected errors. As memory errors would cause incorrect " \
77 "simulation results, gromacs has aborted execution.\n Make sure that your GPU's memory is not " \
78 "overclocked and that the device is properly cooled.\n", (str)
81 #define COMBRULE_CHK_TOL 1e-6
82 #define COMBRULE_SIGMA(sig1, sig2) (((sig1) + (sig2))/2)
83 #define COMBRULE_EPS(eps1, eps2) (sqrt((eps1) * (eps2)))
86 * \brief Convert string to integer type.
87 * \param[in] s String to convert from.
88 * \param[in] f Basefield format flag that takes any of the following I/O
89 * manipulators: dec, hex, oct.
90 * \param[out] t Destination variable to convert to.
93 static gmx_bool from_string(T& t, const string& s, ios_base& (*f)(ios_base&))
96 return !(iss >> f >> t).fail();
100 * \brief Split string around a given delimiter.
101 * \param[in] s String to split.
102 * \param[in] delim Delimiter character.
103 * \returns Vector of strings found in \p s.
105 static vector<string> split(const string &s, char delim)
107 vector<string> elems;
110 while (getline(ss, item, delim))
112 if (item.length() != 0)
113 elems.push_back(item);
119 * \brief Split a string of the form "option=value" into "option" and "value" strings.
120 * This string corresponds to one option and the associated value from the option list
121 * in the mdrun -device argument.
123 * \param[in] s A string containing an "option=value" pair that needs to be split up.
124 * \param[out] opt The name of the option.
125 * \param[out] val Value of the option.
127 static void splitOptionValue(const string &s, string &opt, string &val)
129 size_t eqPos = s.find('=');
130 if (eqPos != string::npos)
132 opt = s.substr(0, eqPos);
133 if (eqPos != s.length()) val = s.substr(eqPos+1);
138 * \brief Compare two strings ignoring case.
139 * This function is in fact a wrapper around the gromacs function gmx_strncasecmp().
140 * \param[in] s1 String.
141 * \param[in] s2 String.
142 * \returns Similarly to the C function strncasecmp(), the return value is an
143 integer less than, equal to, or greater than 0 if \p s1 less than,
144 identical to, or greater than \p s2.
146 static gmx_bool isStringEqNCase(const string& s1, const string& s2)
148 return (gmx_strncasecmp(s1.c_str(), s2.c_str(), max(s1.length(), s2.length())) == 0);
152 * \brief Convert string to upper case.
154 * \param[in] s String to convert to uppercase.
155 * \returns The given string converted to uppercase.
157 static string toUpper(const string &s)
160 std::transform(stmp.begin(), stmp.end(), stmp.begin(), static_cast < int(*)(int) > (toupper));
165 \name Sizes of constant device option arrays GmxOpenMMPlatformOptions#platforms,
166 GmxOpenMMPlatformOptions#memtests, GmxOpenMMPlatformOptions#deviceid,
167 GmxOpenMMPlatformOptions#force_dev. */
169 #define SIZEOF_PLATFORMS 2 // 2
170 #define SIZEOF_MEMTESTS 3
171 #define SIZEOF_DEVICEIDS 1
172 #define SIZEOF_FORCE_DEV 2
174 #define SIZEOF_CHECK_COMBRULE 2
177 /*! Possible platform options in the mdrun -device option. */
178 static const char *devOptStrings[] = { "platform", "deviceid", "memtest", "force-device", "check-combrule" };
180 /*! Enumerated platform options in the mdrun -device option. */
190 * \brief Class to extract and manage the platform options in the mdrun -device option.
193 class GmxOpenMMPlatformOptions
196 GmxOpenMMPlatformOptions(const char *opt);
197 ~GmxOpenMMPlatformOptions() { options.clear(); }
198 string getOptionValue(const string &opt);
199 void remOption(const string &opt);
202 void setOption(const string &opt, const string &val);
204 map<string, string> options; /*!< Data structure to store the option (name, value) pairs. */
206 static const char * const platforms[SIZEOF_PLATFORMS]; /*!< Available OpenMM platforms; size #SIZEOF_PLATFORMS */
207 static const char * const memtests[SIZEOF_MEMTESTS]; /*!< Available types of memory tests, also valid
208 any positive integer >=15; size #SIZEOF_MEMTESTS */
209 static const char * const deviceid[SIZEOF_DEVICEIDS]; /*!< Possible values for deviceid option;
210 also valid any positive integer; size #SIZEOF_DEVICEIDS */
211 static const char * const force_dev[SIZEOF_FORCE_DEV]; /*!< Possible values for for force-device option;
212 size #SIZEOF_FORCE_DEV */
213 static const char * const check_combrule[SIZEOF_CHECK_COMBRULE]; /* XXX temporary debug feature to
214 turn off combination rule check */
217 const char * const GmxOpenMMPlatformOptions::platforms[SIZEOF_PLATFORMS]
218 = {"CUDA", "Reference"};
219 //= { "Reference", "CUDA" /*,"OpenCL"*/ };
220 const char * const GmxOpenMMPlatformOptions::memtests[SIZEOF_MEMTESTS]
221 = { "15", "full", "off" };
222 const char * const GmxOpenMMPlatformOptions::deviceid[SIZEOF_DEVICEIDS]
224 const char * const GmxOpenMMPlatformOptions::force_dev[SIZEOF_FORCE_DEV]
226 const char * const GmxOpenMMPlatformOptions::check_combrule[SIZEOF_CHECK_COMBRULE]
231 * Takes the option list, parses it, checks the options and their values for validity.
232 * When certain options are not provided by the user, as default value the first item
233 * of the respective constant array is taken (GmxOpenMMPlatformOptions#platforms,
234 * GmxOpenMMPlatformOptions#memtests, GmxOpenMMPlatformOptions#deviceid,
235 * GmxOpenMMPlatformOptions#force_dev).
236 * \param[in] optionString Option list part of the mdrun -device parameter.
238 GmxOpenMMPlatformOptions::GmxOpenMMPlatformOptions(const char *optionString)
240 // set default values
241 setOption("platform", platforms[0]);
242 setOption("memtest", memtests[0]);
243 setOption("deviceid", deviceid[0]);
244 setOption("force-device", force_dev[0]);
245 setOption("check-combrule", check_combrule[0]);
247 string opt(optionString);
249 // remove all whitespaces
250 opt.erase(remove_if(opt.begin(), opt.end(), ::isspace), opt.end());
251 // tokenize around ","-s
252 vector<string> tokens = split(opt, ',');
254 for (vector<string>::iterator it = tokens.begin(); it != tokens.end(); ++it)
256 string opt = "", val = "";
257 splitOptionValue(*it, opt, val);
259 if (isStringEqNCase(opt, "platform"))
261 /* no check, this will fail if platform does not exist when we try to set it */
266 if (isStringEqNCase(opt, "memtest"))
268 /* the value has to be an integer >15(s) or "full" OR "off" */
269 if (!isStringEqNCase(val, "full") && !isStringEqNCase(val, "off"))
272 if (!from_string<int>(secs, val, std::dec))
274 gmx_fatal(FARGS, "Invalid value for option memtest option: \"%s\"!", val.c_str());
278 gmx_fatal(FARGS, "Incorrect value for memtest option (%d). "
279 "Memtest needs to run for at least 15s!", secs);
286 if (isStringEqNCase(opt, "deviceid"))
289 if (!from_string<int>(id, val, std::dec) )
291 gmx_fatal(FARGS, "Invalid device id: \"%s\"!", val.c_str());
297 if (isStringEqNCase(opt, "force-device"))
300 if (!isStringEqNCase(val, "yes") && !isStringEqNCase(val, "no"))
302 gmx_fatal(FARGS, "Invalid OpenMM force option: \"%s\"!", val.c_str());
308 if (isStringEqNCase(opt, "check-combrule"))
311 if (!isStringEqNCase(val, "yes") && !isStringEqNCase(val, "no"))
313 gmx_fatal(FARGS, "Invalid OpenMM force option: \"%s\"!", val.c_str());
320 // if we got till here something went wrong
321 gmx_fatal(FARGS, "Invalid OpenMM platform option: \"%s\"!", (*it).c_str());
327 * \brief Getter function.
328 * \param[in] opt Name of the option.
329 * \returns Returns the value associated to an option.
331 string GmxOpenMMPlatformOptions::getOptionValue(const string &opt)
333 map<string, string> :: const_iterator it = options.find(toUpper(opt));
334 if (it != options.end())
345 * \brief Setter function - private, only used from contructor.
346 * \param[in] opt Name of the option.
347 * \param[in] val Value for the option.
349 void GmxOpenMMPlatformOptions::setOption(const string &opt, const string &val)
351 options[toUpper(opt)] = val;
355 * \brief Removes an option with its value from the map structure. If the option
356 * does not exist, returns without any action.
357 * \param[in] opt Name of the option.
359 void GmxOpenMMPlatformOptions::remOption(const string &opt)
361 options.erase(toUpper(opt));
365 * \brief Print option-value pairs to a file (debugging function).
367 void GmxOpenMMPlatformOptions::print()
369 cout << ">> Platform options: " << endl
370 << ">> platform = " << getOptionValue("platform") << endl
371 << ">> deviceID = " << getOptionValue("deviceid") << endl
372 << ">> memtest = " << getOptionValue("memtest") << endl
373 << ">> force-device = " << getOptionValue("force-device") << endl;
377 * \brief Container for OpenMM related data structures that represent the bridge
378 * between the Gromacs data-structures and the OpenMM library and is but it's
379 * only passed through the API functions as void to disable direct access.
384 System* system; /*! The system to simulate. */
385 Context* context; /*! The OpenMM context in which the simulation is carried out. */
386 Integrator* integrator; /*! The integrator used in the simulation. */
387 gmx_bool removeCM; /*! If \true remove venter of motion, false otherwise. */
388 GmxOpenMMPlatformOptions *platformOpt; /*! Platform options. */
392 * \brief Runs memtest on the GPU that has alreaby been initialized by OpenMM.
393 * \param[in] fplog Pointer to gromacs log file.
394 * \param[in] devId Device id of the GPU to run the test on.
395 Note: as OpenMM previously creates the context,for now this is always -1.
396 * \param[in] pre_post Contains either "Pre" or "Post" just to be able to differentiate in
397 * stdout messages/log between memtest carried out before and after simulation.
398 * \param[in] opt Pointer to platform options object.
400 static void runMemtest(FILE* fplog, int devId, const char* pre_post, GmxOpenMMPlatformOptions *opt)
402 char strout_buf[STRLEN];
405 string s = opt->getOptionValue("memtest");
406 const char *test_type = s.c_str();
408 if (!gmx_strcasecmp(test_type, "off"))
414 if (!gmx_strcasecmp(test_type, "full"))
420 from_string<int>(which_test, test_type, std::dec);
426 gmx_fatal(FARGS, "Amount of seconds for memetest is negative (%d). ", which_test);
431 case 0: /* no memtest */
432 sprintf(strout_buf, "%s-simulation GPU memtest skipped. Note, that faulty memory can cause "
433 "incorrect results!", pre_post);
434 fprintf(fplog, "%s\n", strout_buf);
435 gmx_warning(strout_buf);
438 case 1: /* quick memtest */
439 fprintf(fplog, "%s-simulation %s GPU memtest in progress...\n", pre_post, test_type);
440 fprintf(stdout, "\n%s-simulation %s GPU memtest in progress...", pre_post, test_type);
443 res = do_quick_memtest(devId);
446 case 2: /* full memtest */
447 fprintf(fplog, "%s-simulation %s memtest in progress...\n", pre_post, test_type);
448 fprintf(stdout, "\n%s-simulation %s memtest in progress...", pre_post, test_type);
451 res = do_full_memtest(devId);
454 default: /* timed memtest */
455 fprintf(fplog, "%s-simulation ~%ds memtest in progress...\n", pre_post, which_test);
456 fprintf(stdout, "\n%s-simulation ~%ds memtest in progress...", pre_post, which_test);
459 res = do_timed_memtest(devId, which_test);
466 gmx_fatal(FARGS, MEM_ERR_MSG(pre_post));
470 fprintf(fplog, "Memory test completed without errors.\n");
472 fprintf(stdout, "done, no errors detected\n");
479 * \brief Convert Lennard-Jones parameters c12 and c6 to sigma and epsilon.
484 * \param[out] epsilon
486 static void convert_c_12_6(double c12, double c6, double *sigma, double *epsilon)
488 if (c12 == 0 && c6 == 0)
493 else if (c12 > 0 && c6 > 0)
495 *epsilon = (c6*c6)/(4.0*c12);
496 *sigma = pow(c12/c6, 1.0/6.0);
500 gmx_fatal(FARGS,"OpenMM only supports c6 > 0 and c12 > 0 or c6 = c12 = 0.");
505 * \brief Does gromacs option checking.
507 * Checks the gromacs mdp options for features unsupported in OpenMM, case in which
508 * interrupts the execution. It also warns the user about pecularities of OpenMM
510 * \param[in] fplog Gromacs log file pointer.
511 * \param[in] ir Gromacs input parameters, see ::t_inputrec
512 * \param[in] top Gromacs node local topology, \see gmx_localtop_t
513 * \param[in] state Gromacs state structure \see ::t_state
514 * \param[in] mdatoms Gromacs atom parameters, \see ::t_mdatoms
515 * \param[in] fr \see ::t_forcerec
516 * \param[in] state Gromacs systems state, \see ::t_state
518 static void checkGmxOptions(FILE* fplog, GmxOpenMMPlatformOptions *opt,
519 t_inputrec *ir, gmx_localtop_t *top,
520 t_forcerec *fr, t_state *state)
522 char warn_buf[STRLEN];
525 double sigma_ij=0, sigma_ji=0, sigma_ii=0, sigma_jj=0, sigma_comb;
526 double eps_ij=0, eps_ji=0, eps_ii=0, eps_jj=0, eps_comb;
528 /* Abort if unsupported critical options are present */
533 gmx_warning( "OpenMM does not support leap-frog, will use velocity-verlet integrator.");
536 if ( (ir->eI != eiMD) &&
538 (ir->eI != eiVVAK) &&
543 gmx_fatal(FARGS, "OpenMM supports only the following integrators: md/md-vv/md-vv-avek, sd/sd1, and bd.");
547 if ( !(ir->coulombtype == eelPME ||
548 EEL_RF(ir->coulombtype) ||
549 ir->coulombtype == eelRF ||
550 ir->coulombtype == eelEWALD ||
552 (ir->coulombtype == eelCUT && ir->rcoulomb == 0 && ir->rvdw == 0) ||
553 // we could have cut-off combined with GBSA (openmm will use RF)
554 ir->implicit_solvent == eisGBSA) )
556 gmx_fatal(FARGS,"OpenMM supports only the following methods for electrostatics: "
557 "NoCutoff (i.e. rcoulomb = rvdw = 0 ),Reaction-Field, Ewald or PME.");
560 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf != 0)
562 // openmm has epsilon_rf=inf hard-coded
563 gmx_warning("OpenMM will use a Reaction-Field epsilon of infinity instead of %g.",ir->epsilon_rf);
566 if (ir->etc != etcNO &&
571 gmx_warning("OpenMM supports only Andersen thermostat with the md/md-vv/md-vv-avek integrators.");
574 if (ir->implicit_solvent == eisGBSA &&
575 ir->gb_algorithm != egbOBC )
577 gmx_warning("OpenMM does not support the specified algorithm for Generalized Born, will use OBC instead.");
580 if (ir->opts.ngtc > 1)
581 gmx_fatal(FARGS,"OpenMM does not support multiple temperature coupling groups.");
583 if (ir->epc != epcNO)
584 gmx_warning("OpenMM supports only Monte Carlo barostat for pressure coupling.");
586 if (ir->opts.annealing[0])
587 gmx_fatal(FARGS,"OpenMM does not support simulated annealing.");
589 if (top->idef.il[F_CONSTR].nr > 0 && ir->eConstrAlg != econtSHAKE)
590 gmx_warning("OpenMM provides contraints as a combination "
591 "of SHAKE, SETTLE and CCMA. Accuracy is based on the SHAKE tolerance set "
592 "by the \"shake_tol\" option.");
595 gmx_fatal(FARGS,"OpenMM does not support walls.");
597 if (ir->ePull != epullNO)
598 gmx_fatal(FARGS,"OpenMM does not support pulling.");
600 /* check for interaction types */
601 for (i = 0; i < F_EPOT; i++)
603 if (!(i == F_CONSTR ||
607 i == F_UREY_BRADLEY ||
614 i == F_GB12 || /* The GB parameters are hardcoded both in */
615 i == F_GB13 || /* Gromacs and OpenMM */
617 top->idef.il[i].nr > 0)
619 gmx_fatal(FARGS, "OpenMM does not support (some) of the provided interaction "
620 "type(s) (%s) ", interaction_function[i].longname);
624 if (ir->efep != efepNO)
625 gmx_fatal(FARGS,"OpenMM does not support free energy calculations.");
627 if (ir->opts.ngacc > 1)
628 gmx_fatal(FARGS,"OpenMM does not support non-equilibrium MD (accelerated groups).");
630 if (IR_ELEC_FIELD(*ir))
631 gmx_fatal(FARGS,"OpenMM does not support electric fields.");
634 gmx_fatal(FARGS,"OpenMM does not support QMMM calculations.");
636 if (ir->rcoulomb != ir->rvdw)
637 gmx_fatal(FARGS,"OpenMM uses a single cutoff for both Coulomb "
638 "and VdW interactions. Please set rcoulomb equal to rvdw.");
640 if (EEL_FULL(ir->coulombtype))
642 if (ir->ewald_geometry == eewg3DC)
643 gmx_fatal(FARGS,"OpenMM supports only Ewald 3D geometry.");
644 if (ir->epsilon_surface != 0)
645 gmx_fatal(FARGS,"OpenMM does not support dipole correction in Ewald summation.");
648 if (TRICLINIC(state->box))
650 gmx_fatal(FARGS,"OpenMM does not support triclinic unit cells.");
653 /* XXX this is just debugging code to disable the combination rule check */
654 if ( isStringEqNCase(opt->getOptionValue("check-combrule"), "yes") )
656 /* As OpenMM by default uses hardcoded combination rules
657 sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j)
658 we need to check whether the force field params obey this
659 and if not, we can't use this force field so we exit
660 grace-fatal-fully. */
661 real *nbfp = fr->nbfp;
665 fprintf(debug, ">> Atom parameters: <<\n%10s%5s %5s %5s %5s COMB\n",
666 "", "i-j", "j-i", "i-i", "j-j");
668 /* loop over all i-j atom pairs and verify if
669 sigma_ij = sigma_ji = sigma_comb and eps_ij = eps_ji = eps_comb */
670 for (i = 0; i < natoms; i++)
673 /* nbfp now includes the 6.0/12.0 prefactors to save flops in kernels */
674 c12 = C12(nbfp, natoms, i, i)/12.0;
675 c6 = C6(nbfp, natoms, i, i)/6.0;
676 convert_c_12_6(c12, c6, &sigma_ii, &eps_ii);
678 for (j = 0; j < i; j++)
681 c12 = C12(nbfp, natoms, i, j)/12.0;
682 c6 = C6(nbfp, natoms, i, j)/6.0;
683 convert_c_12_6(c12, c6, &sigma_ij, &eps_ij);
685 c12 = C12(nbfp, natoms, j, i)/12.0;
686 c6 = C6(nbfp, natoms, j, i)/6.0;
687 convert_c_12_6(c12, c6, &sigma_ji, &eps_ji);
689 c12 = C12(nbfp, natoms, j, j)/12.0;
690 c6 = C6(nbfp, natoms, j, j)/6.0;
691 convert_c_12_6(c12, c6, &sigma_jj, &eps_jj);
692 /* OpenMM hardcoded combination rules */
693 sigma_comb = COMBRULE_SIGMA(sigma_ii, sigma_jj);
694 eps_comb = COMBRULE_EPS(eps_ii, eps_jj);
698 fprintf(debug, "i=%-3d j=%-3d", i, j);
699 fprintf(debug, "%-11s", "sigma");
700 fprintf(debug, "%5.3f %5.3f %5.3f %5.3f %5.3f\n",
701 sigma_ij, sigma_ji, sigma_ii, sigma_jj, sigma_comb);
702 fprintf(debug, "%11s%-11s", "", "epsilon");
703 fprintf(debug, "%5.3f %5.3f %5.3f %5.3f %5.3f\n",
704 eps_ij, eps_ji, eps_ii, eps_jj, eps_comb);
707 /* check the values against the rule used by omm */
708 if((fabs(eps_ij) > COMBRULE_CHK_TOL &&
709 fabs(eps_ji) > COMBRULE_CHK_TOL) &&
710 (fabs(sigma_comb - sigma_ij) > COMBRULE_CHK_TOL ||
711 fabs(sigma_comb - sigma_ji) > COMBRULE_CHK_TOL ||
712 fabs(eps_comb - eps_ij) > COMBRULE_CHK_TOL ||
713 fabs(eps_comb - eps_ji) > COMBRULE_CHK_TOL ))
716 "The combination rules of the used force-field do not "
717 "match the one supported by OpenMM: "
718 "sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j). "
719 "Switch to a force-field that uses these rules in order to "
720 "simulate this system using OpenMM.\n");
724 if (debug) { fprintf(debug, ">><<\n\n"); }
726 /* if we got here, log that everything is fine */
729 fprintf(debug, ">> The combination rule of the used force matches the one used by OpenMM.\n");
731 fprintf(fplog, "The combination rule of the used force field matches the one used by OpenMM.\n");
733 } /* if (are we checking the combination rules) ... */
738 * \brief Initialize OpenMM, run sanity/consistency checks, and return a pointer to
741 * Various gromacs data structures are passed that contain the parameters, state and
742 * other porperties of the system to simulate. These serve as input for initializing
743 * OpenMM. Besides, a set of misc action are taken:
744 * - OpenMM plugins are loaded;
745 * - platform options in \p platformOptStr are parsed and checked;
746 * - Gromacs parameters are checked for OpenMM support and consistency;
747 * - after the OpenMM is initialized memtest executed in the same GPU context.
749 * \param[in] fplog Gromacs log file handler.
750 * \param[in] platformOptStr Platform option string.
751 * \param[in] ir The Gromacs input parameters, see ::t_inputrec
752 * \param[in] top_global Gromacs system toppology, \see ::gmx_mtop_t
753 * \param[in] top Gromacs node local topology, \see gmx_localtop_t
754 * \param[in] mdatoms Gromacs atom parameters, \see ::t_mdatoms
755 * \param[in] fr \see ::t_forcerec
756 * \param[in] state Gromacs systems state, \see ::t_state
757 * \returns Pointer to a
760 void* openmm_init(FILE *fplog, const char *platformOptStr,
762 gmx_mtop_t *top_global, gmx_localtop_t *top,
763 t_mdatoms *mdatoms, t_forcerec *fr, t_state *state)
766 char warn_buf[STRLEN];
767 static gmx_bool hasLoadedPlugins = false;
768 string usedPluginDir;
773 if (!hasLoadedPlugins)
775 vector<string> loadedPlugins;
776 /* Look for OpenMM plugins at various locations (listed in order of priority):
777 - on the path in OPENMM_PLUGIN_DIR environment variable if this is specified
778 - on the path in the OPENMM_PLUGIN_DIR macro that is set by the build script
779 - at the default location assumed by OpenMM
782 char *pluginDir = getenv("OPENMM_PLUGIN_DIR");
784 /* no env var or empty */
785 if (pluginDir != NULL && *pluginDir != '\0')
787 loadedPlugins = Platform::loadPluginsFromDirectory(pluginDir);
788 if (!loadedPlugins.empty())
790 hasLoadedPlugins = true;
791 usedPluginDir = pluginDir;
795 gmx_fatal(FARGS, "The directory provided in the OPENMM_PLUGIN_DIR environment variable "
796 "(%s) does not contain valid OpenMM plugins. Check your OpenMM installation!",
801 /* macro set at build time */
802 #ifdef OPENMM_PLUGIN_DIR
803 if (!hasLoadedPlugins)
805 loadedPlugins = Platform::loadPluginsFromDirectory(OPENMM_PLUGIN_DIR);
806 if (!loadedPlugins.empty())
808 hasLoadedPlugins = true;
809 usedPluginDir = OPENMM_PLUGIN_DIR;
813 /* default loocation */
814 if (!hasLoadedPlugins)
816 loadedPlugins = Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
817 if (!loadedPlugins.empty())
819 hasLoadedPlugins = true;
820 usedPluginDir = Platform::getDefaultPluginsDirectory();
824 /* if there are still no plugins loaded there won't be any */
825 if (!hasLoadedPlugins)
827 gmx_fatal(FARGS, "No OpenMM plugins were found! You can provide the"
828 " plugin directory in the OPENMM_PLUGIN_DIR environment variable.", pluginDir);
831 fprintf(fplog, "\nOpenMM plugins loaded from directory %s:\t", usedPluginDir.c_str());
832 for (int i = 0; i < (int)loadedPlugins.size(); i++)
834 fprintf(fplog, "%s, ", loadedPlugins[i].c_str());
836 fprintf(fplog, "\n");
839 /* parse option string */
840 GmxOpenMMPlatformOptions *opt = new GmxOpenMMPlatformOptions(platformOptStr);
841 devId = atoi(opt->getOptionValue("deviceid").c_str());
848 /* check wheter Gromacs options compatibility with OpenMM */
849 checkGmxOptions(fplog, opt, ir, top, fr, state);
851 /* Create the system. */
852 const t_idef& idef = top->idef;
853 const int numAtoms = top_global->natoms;
854 const int numConstraints = idef.il[F_CONSTR].nr/3;
855 const int numSettle = idef.il[F_SETTLE].nr/2;
856 const int numBonds = idef.il[F_BONDS].nr/3;
857 const int numHarmonic = idef.il[F_HARMONIC].nr/3;
858 const int numUB = idef.il[F_UREY_BRADLEY].nr/4;
859 const int numAngles = idef.il[F_ANGLES].nr/4;
860 const int numPeriodic = idef.il[F_PDIHS].nr/5;
861 const int numPeriodicImproper = idef.il[F_PIDIHS].nr/5;
862 const int numRB = idef.il[F_RBDIHS].nr/5;
863 const int numImproperDih = idef.il[F_IDIHS].nr/5;
864 const int num14 = idef.il[F_LJ14].nr/3;
865 System* sys = new System();
867 sys->addForce(new CMMotionRemover(ir->nstcomm));
869 /* Set bonded force field terms. */
872 * CUDA platform currently doesn't support more than one
873 * instance of a force object, so we pack all forces that
874 * use the same form into one.
877 const int* bondAtoms = (int*) idef.il[F_BONDS].iatoms;
878 HarmonicBondForce* bondForce = new HarmonicBondForce();
879 sys->addForce(bondForce);
881 for (int i = 0; i < numBonds; ++i)
883 int type = bondAtoms[offset++];
884 int atom1 = bondAtoms[offset++];
885 int atom2 = bondAtoms[offset++];
886 bondForce->addBond(atom1, atom2,
887 idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
890 const int* harmonicAtoms = (int*) idef.il[F_HARMONIC].iatoms;
892 for (int i = 0; i < numHarmonic; ++i)
894 int type = harmonicAtoms[offset++];
895 int atom1 = harmonicAtoms[offset++];
896 int atom2 = harmonicAtoms[offset++];
897 bondForce->addBond(atom1, atom2,
898 idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
901 /* Set the angle force field terms */
902 const int* angleAtoms = (int*) idef.il[F_ANGLES].iatoms;
903 HarmonicAngleForce* angleForce = new HarmonicAngleForce();
904 sys->addForce(angleForce);
906 for (int i = 0; i < numAngles; ++i)
908 int type = angleAtoms[offset++];
909 int atom1 = angleAtoms[offset++];
910 int atom2 = angleAtoms[offset++];
911 int atom3 = angleAtoms[offset++];
912 angleForce->addAngle(atom1, atom2, atom3,
913 idef.iparams[type].harmonic.rA*M_PI/180.0, idef.iparams[type].harmonic.krA);
916 /* Urey-Bradley includes both the angle and bond potential for 1-3 interactions */
917 const int* ubAtoms = (int*) idef.il[F_UREY_BRADLEY].iatoms;
918 /* HarmonicBondForce* ubBondForce = new HarmonicBondForce(); */
919 /* HarmonicAngleForce* ubAngleForce = new HarmonicAngleForce(); */
920 /* sys->addForce(ubBondForce); */
921 /* sys->addForce(ubAngleForce); */
923 for (int i = 0; i < numUB; ++i)
925 int type = ubAtoms[offset++];
926 int atom1 = ubAtoms[offset++];
927 int atom2 = ubAtoms[offset++];
928 int atom3 = ubAtoms[offset++];
929 /* ubBondForce->addBond(atom1, atom3, */
930 bondForce->addBond(atom1, atom3,
931 idef.iparams[type].u_b.r13A, idef.iparams[type].u_b.kUBA);
932 /* ubAngleForce->addAngle(atom1, atom2, atom3, */
933 angleForce->addAngle(atom1, atom2, atom3,
934 idef.iparams[type].u_b.thetaA*M_PI/180.0, idef.iparams[type].u_b.kthetaA);
937 /* Set proper dihedral terms */
938 const int* periodicAtoms = (int*) idef.il[F_PDIHS].iatoms;
939 PeriodicTorsionForce* periodicForce = new PeriodicTorsionForce();
940 sys->addForce(periodicForce);
942 for (int i = 0; i < numPeriodic; ++i)
944 int type = periodicAtoms[offset++];
945 int atom1 = periodicAtoms[offset++];
946 int atom2 = periodicAtoms[offset++];
947 int atom3 = periodicAtoms[offset++];
948 int atom4 = periodicAtoms[offset++];
949 periodicForce->addTorsion(atom1, atom2, atom3, atom4,
950 idef.iparams[type].pdihs.mult,
951 idef.iparams[type].pdihs.phiA*M_PI/180.0,
952 idef.iparams[type].pdihs.cpA);
955 /* Set improper dihedral terms that are represented by a periodic function (as in AMBER FF) */
956 const int* periodicImproperAtoms = (int*) idef.il[F_PIDIHS].iatoms;
957 /* PeriodicTorsionForce* periodicImproperForce = new PeriodicTorsionForce(); */
958 /* sys->addForce(periodicImproperForce); */
960 for (int i = 0; i < numPeriodicImproper; ++i)
962 int type = periodicImproperAtoms[offset++];
963 int atom1 = periodicImproperAtoms[offset++];
964 int atom2 = periodicImproperAtoms[offset++];
965 int atom3 = periodicImproperAtoms[offset++];
966 int atom4 = periodicImproperAtoms[offset++];
967 /* periodicImproperForce->addTorsion(atom1, atom2, atom3, atom4, */
968 periodicForce->addTorsion(atom1, atom2, atom3, atom4,
969 idef.iparams[type].pdihs.mult,
970 idef.iparams[type].pdihs.phiA*M_PI/180.0,
971 idef.iparams[type].pdihs.cpA);
974 /* Ryckaert-Bellemans dihedrals */
975 const int* rbAtoms = (int*) idef.il[F_RBDIHS].iatoms;
976 RBTorsionForce* rbForce = new RBTorsionForce();
977 sys->addForce(rbForce);
979 for (int i = 0; i < numRB; ++i)
981 int type = rbAtoms[offset++];
982 int atom1 = rbAtoms[offset++];
983 int atom2 = rbAtoms[offset++];
984 int atom3 = rbAtoms[offset++];
985 int atom4 = rbAtoms[offset++];
986 rbForce->addTorsion(atom1, atom2, atom3, atom4,
987 idef.iparams[type].rbdihs.rbcA[0], idef.iparams[type].rbdihs.rbcA[1],
988 idef.iparams[type].rbdihs.rbcA[2], idef.iparams[type].rbdihs.rbcA[3],
989 idef.iparams[type].rbdihs.rbcA[4], idef.iparams[type].rbdihs.rbcA[5]);
992 /* Set improper dihedral terms (as in CHARMM FF) */
993 const int* improperDihAtoms = (int*) idef.il[F_IDIHS].iatoms;
994 CustomTorsionForce* improperDihForce = new CustomTorsionForce("2.0*k*asin(sin((theta-theta0)/2))^2");
995 sys->addForce(improperDihForce);
996 improperDihForce->addPerTorsionParameter("k");
997 improperDihForce->addPerTorsionParameter("theta0");
998 vector<double> improperDihParameters(2);
1000 for (int i = 0; i < numImproperDih; ++i)
1002 int type = improperDihAtoms[offset++];
1003 int atom1 = improperDihAtoms[offset++];
1004 int atom2 = improperDihAtoms[offset++];
1005 int atom3 = improperDihAtoms[offset++];
1006 int atom4 = improperDihAtoms[offset++];
1007 improperDihParameters[0] = idef.iparams[type].harmonic.krA;
1008 improperDihParameters[1] = idef.iparams[type].harmonic.rA*M_PI/180.0;
1009 improperDihForce->addTorsion(atom1, atom2, atom3, atom4,
1010 improperDihParameters);
1013 /* Set nonbonded parameters and masses. */
1014 int ntypes = fr->ntype;
1015 int* types = mdatoms->typeA;
1016 real* nbfp = fr->nbfp;
1017 real* charges = mdatoms->chargeA;
1018 real* masses = mdatoms->massT;
1019 NonbondedForce* nonbondedForce = new NonbondedForce();
1020 sys->addForce(nonbondedForce);
1025 if (ir->rcoulomb == 0)
1027 nonbondedForce->setNonbondedMethod(NonbondedForce::NoCutoff);
1031 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
1035 switch (ir->coulombtype)
1042 nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
1046 nonbondedForce->setNonbondedMethod(NonbondedForce::Ewald);
1050 nonbondedForce->setNonbondedMethod(NonbondedForce::PME);
1054 gmx_fatal(FARGS,"Internal error: you should not see this message, it means that the"
1055 "electrosatics option check failed. Please report this error!");
1057 sys->setDefaultPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0),
1058 Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));
1059 nonbondedForce->setCutoffDistance(ir->rcoulomb);
1063 gmx_fatal(FARGS,"OpenMM supports only full periodic boundary conditions "
1064 "(pbc = xyz), or none (pbc = no).");
1068 /* Fix for PME and Ewald error tolerance
1070 * OpenMM uses approximate formulas to calculate the Ewald parameter:
1071 * alpha = (1.0/cutoff)*sqrt(-log(2.0*tolerlance));
1072 * and the grid spacing for PME:
1073 * gridX = ceil(2*alpha*box[0][0]/3*(pow(tol, 0.2)))
1074 * gridY = ceil(2*alpha*box[1][1]/3*(pow(tol, 0.2)));
1075 * gridZ = ceil(2*alpha*box[2][2]/3*(pow(tol, 0.2)));
1078 * If the default ewald_rtol=1e-5 is used we silently adjust the value to the
1079 * OpenMM default of 5e-4 otherwise a warning is issued about the action taken.
1082 double corr_ewald_rtol = 50.0 * ir->ewald_rtol;
1083 if ((ir->ePBC == epbcXYZ) &&
1084 (ir->coulombtype == eelEWALD || ir->coulombtype == eelPME))
1088 fprintf(debug, ">> ewald_rtol = %e (corrected = %e) \n",
1089 ir->ewald_rtol, corr_ewald_rtol);
1092 if (fabs(ir->ewald_rtol - 1e-5) > 1e-10)
1094 gmx_warning("OpenMM uses the ewald_rtol parameter with approximate formulas "
1095 "to calculate the alpha and grid spacing parameters of the Ewald "
1096 "and PME methods. This tolerance need to be corrected in order to get "
1097 "settings close to the ones used in GROMACS. Although the internal correction "
1098 "should work for any reasonable value of ewald_rtol, using values other than "
1099 "the default 1e-5 might cause incorrect behavior.");
1101 if (corr_ewald_rtol > 1)
1103 gmx_fatal(FARGS, "The ewald_rtol accuracy term is >1 after the "
1104 "adjustment for OpenMM (%e)", corr_ewald_rtol);
1107 nonbondedForce->setEwaldErrorTolerance(corr_ewald_rtol);
1110 for (int i = 0; i < numAtoms; ++i)
1112 /* nbfp now includes the 6.0/12.0 derivative prefactors to save flops in kernels*/
1113 double c12 = nbfp[types[i]*2*ntypes+types[i]*2+1]/12.0;
1114 double c6 = nbfp[types[i]*2*ntypes+types[i]*2]/6.0;
1115 double sigma=0.0, epsilon=0.0;
1116 convert_c_12_6(c12, c6, &sigma, &epsilon);
1117 nonbondedForce->addParticle(charges[i], sigma, epsilon);
1118 sys->addParticle(masses[i]);
1121 // Build a table of all exclusions.
1122 vector<set<int> > exclusions(numAtoms);
1123 for (int i = 0; i < numAtoms; i++)
1125 int start = top->excls.index[i];
1126 int end = top->excls.index[i+1];
1127 for (int j = start; j < end; j++)
1128 exclusions[i].insert(top->excls.a[j]);
1131 // Record the 1-4 interactions, and remove them from the list of exclusions.
1132 const int* nb14Atoms = (int*) idef.il[F_LJ14].iatoms;
1134 for (int i = 0; i < num14; ++i)
1136 int type = nb14Atoms[offset++];
1137 int atom1 = nb14Atoms[offset++];
1138 int atom2 = nb14Atoms[offset++];
1139 double sigma=0, epsilon=0;
1140 convert_c_12_6(idef.iparams[type].lj14.c12A,
1141 idef.iparams[type].lj14.c6A,
1143 nonbondedForce->addException(atom1, atom2,
1144 fr->fudgeQQ*charges[atom1]*charges[atom2], sigma, epsilon);
1145 exclusions[atom1].erase(atom2);
1146 exclusions[atom2].erase(atom1);
1149 // Record exclusions.
1150 for (int i = 0; i < numAtoms; i++)
1152 for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
1156 nonbondedForce->addException(i, *iter, 0.0, 1.0, 0.0);
1161 // Add GBSA if needed.
1162 if (ir->implicit_solvent == eisGBSA)
1164 gmx_warning("The OBC scale factors alpha, beta and gamma are hardcoded in OpenMM with the default Gromacs values.");
1165 t_atoms atoms = gmx_mtop_global_atoms(top_global);
1166 GBSAOBCForce* gbsa = new GBSAOBCForce();
1168 sys->addForce(gbsa);
1169 gbsa->setSoluteDielectric(ir->epsilon_r);
1170 gbsa->setSolventDielectric(ir->gb_epsilon_solvent);
1171 gbsa->setCutoffDistance(nonbondedForce->getCutoffDistance());
1172 if (nonbondedForce->getNonbondedMethod() == NonbondedForce::NoCutoff)
1173 gbsa->setNonbondedMethod(GBSAOBCForce::NoCutoff);
1174 else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
1175 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
1176 else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffPeriodic)
1177 gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
1179 gmx_fatal(FARGS,"OpenMM supports only Reaction-Field electrostatics with OBC/GBSA.");
1181 for (int i = 0; i < numAtoms; ++i)
1183 gbsa->addParticle(charges[i],
1184 top_global->atomtypes.gb_radius[atoms.atom[i].type],
1185 top_global->atomtypes.S_hct[atoms.atom[i].type]);
1187 free_t_atoms(&atoms, FALSE);
1191 const int* constraintAtoms = (int*) idef.il[F_CONSTR].iatoms;
1193 for (int i = 0; i < numConstraints; ++i)
1195 int type = constraintAtoms[offset++];
1196 int atom1 = constraintAtoms[offset++];
1197 int atom2 = constraintAtoms[offset++];
1198 sys->addConstraint(atom1, atom2, idef.iparams[type].constr.dA);
1200 const int* settleAtoms = (int*) idef.il[F_SETTLE].iatoms;
1202 for (int i = 0; i < numSettle; ++i)
1204 int type = settleAtoms[offset++];
1205 int oxygen = settleAtoms[offset++];
1206 sys->addConstraint(oxygen, oxygen+1, idef.iparams[type].settle.doh);
1207 sys->addConstraint(oxygen, oxygen+2, idef.iparams[type].settle.doh);
1208 sys->addConstraint(oxygen+1, oxygen+2, idef.iparams[type].settle.dhh);
1211 // Create an integrator for simulating the system.
1212 double friction = (ir->opts.tau_t[0] == 0.0 ? 0.0 : 1.0/ir->opts.tau_t[0]);
1216 integ = new BrownianIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1217 static_cast<BrownianIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1219 else if (EI_SD(ir->eI))
1221 integ = new LangevinIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
1222 static_cast<LangevinIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed);
1226 integ = new VerletIntegrator(ir->delta_t);
1227 if ( ir->etc != etcNO)
1229 AndersenThermostat* thermostat = new AndersenThermostat(ir->opts.ref_t[0], friction);
1230 sys->addForce(thermostat);
1234 // Add pressure coupling
1235 if (ir->epc != epcNO)
1237 // convert gromacs pressure tensor to a scalar
1238 double pressure = (ir->ref_p[0][0] + ir->ref_p[1][1] + ir->ref_p[2][2]) / 3.0;
1239 int frequency = int(ir->tau_p / ir->delta_t); // update frequency in time steps
1240 if (frequency < 1) frequency = 1;
1241 double temperature = ir->opts.ref_t[0]; // in kelvin
1242 sys->addForce(new MonteCarloBarostat(pressure, temperature, frequency));
1245 integ->setConstraintTolerance(ir->shake_tol);
1247 // Create a context and initialize it.
1248 Context* context = NULL;
1251 OpenMM could automatically select the "best" GPU, however we're not't
1252 going to let it do that for now, as the current algorithm is very rudimentary
1253 and we anyway support only CUDA.
1254 if (platformOptStr == NULL || platformOptStr == "")
1256 context = new Context(*sys, *integ);
1261 /* which platform should we use */
1262 for (int i = 0; i < (int)Platform::getNumPlatforms() && context == NULL; i++)
1264 if (isStringEqNCase(opt->getOptionValue("platform"), Platform::getPlatform(i).getName()))
1266 Platform& platform = Platform::getPlatform(i);
1267 // set standard properties
1268 platform.setPropertyDefaultValue("CudaDevice", opt->getOptionValue("deviceid"));
1269 // TODO add extra properties
1270 context = new Context(*sys, *integ, platform);
1273 if (context == NULL)
1275 gmx_fatal(FARGS, "The requested platform \"%s\" could not be found.",
1276 opt->getOptionValue("platform").c_str());
1280 Platform& platform = context->getPlatform();
1281 fprintf(fplog, "Gromacs will use the OpenMM platform: %s\n", platform.getName().c_str());
1283 const vector<string>& properties = platform.getPropertyNames();
1286 for (int i = 0; i < (int)properties.size(); i++)
1288 fprintf(debug, ">> %s: %s\n", properties[i].c_str(),
1289 platform.getPropertyValue(*context, properties[i]).c_str());
1294 if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1297 if (!from_string<int>(tmp, platform.getPropertyValue(*context, "CudaDevice"), std::dec))
1299 gmx_fatal(FARGS, "Internal error: couldn't determine the device selected by OpenMM");
1303 /* For now this is just to double-check if OpenMM selected the GPU we wanted,
1304 but when we'll let OpenMM select the GPU automatically, it will query the deviceId.
1308 gmx_fatal(FARGS, "Internal error: OpenMM is using device #%d"
1309 "while initialized for device #%d", tmp, devId);
1312 /* check GPU compatibility */
1313 char gpuname[STRLEN];
1314 devId = atoi(opt->getOptionValue("deviceid").c_str());
1315 if (!is_gmx_openmm_supported_gpu(-1, gpuname))
1317 if (!gmx_strcasecmp(opt->getOptionValue("force-device").c_str(), "yes"))
1319 sprintf(warn_buf, "Non-supported GPU selected (#%d, %s), forced continuing."
1320 "Note, that the simulation can be slow or it migth even crash.",
1322 fprintf(fplog, "%s\n", warn_buf);
1323 gmx_warning(warn_buf);
1327 gmx_fatal(FARGS, "The selected GPU (#%d, %s) is not supported by Gromacs! "
1328 "Most probably you have a low-end GPU which would not perform well, "
1329 "or new hardware that has not been tested with the current release. "
1330 "If you still want to try using the device, use the force-device=yes option.",
1336 fprintf(fplog, "Gromacs will run on the GPU #%d (%s).\n", devId, gpuname);
1341 if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
1343 /* pre-simulation memtest */
1344 runMemtest(fplog, -1, "Pre", opt);
1347 vector<Vec3> pos(numAtoms);
1348 vector<Vec3> vel(numAtoms);
1349 for (int i = 0; i < numAtoms; ++i)
1351 pos[i] = Vec3(state->x[i][0], state->x[i][1], state->x[i][2]);
1352 vel[i] = Vec3(state->v[i][0], state->v[i][1], state->v[i][2]);
1354 context->setPositions(pos);
1355 context->setVelocities(vel);
1357 // Return a structure containing the system, integrator, and context.
1358 OpenMMData* data = new OpenMMData();
1360 data->integrator = integ;
1361 data->context = context;
1362 data->removeCM = (ir->nstcomm > 0);
1363 data->platformOpt = opt;
1366 catch (std::exception& e)
1368 gmx_fatal(FARGS, "OpenMM exception caught while initializating: %s", e.what());
1370 return NULL; /* just to avoid warnings */
1374 * \brief Integrate one step.
1376 * \param[in] data OpenMMData object created by openmm_init().
1378 void openmm_take_one_step(void* data)
1380 // static int step = 0; printf("----> taking step #%d\n", step++);
1383 static_cast<OpenMMData*>(data)->integrator->step(1);
1385 catch (std::exception& e)
1387 gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1392 * \brief Integrate n steps.
1394 * \param[in] data OpenMMData object created by openmm_init().
1396 void openmm_take_steps(void* data, int nstep)
1400 static_cast<OpenMMData*>(data)->integrator->step(nstep);
1402 catch (std::exception& e)
1404 gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
1409 * \brief Clean up the data structures cretead for OpenMM.
1411 * \param[in] log Log file pointer.
1412 * \param[in] data OpenMMData object created by openmm_init().
1414 void openmm_cleanup(FILE* fplog, void* data)
1416 OpenMMData* d = static_cast<OpenMMData*>(data);
1418 if (isStringEqNCase(d->platformOpt->getOptionValue("platform"), "CUDA"))
1420 /* post-simulation memtest */
1421 runMemtest(fplog, -1, "Post", d->platformOpt);
1424 delete d->integrator;
1426 delete d->platformOpt;
1431 * \brief Copy the current state information from OpenMM into the Gromacs data structures.
1433 * This function results in the requested proprties to be copied from the
1434 * GPU to host. As this represents a bottleneck, the frequency of pulling data
1435 * should be minimized.
1437 * \param[in] data OpenMMData object created by openmm_init().
1438 * \param[out] time Simulation time for which the state was created.
1439 * \param[out] state State of the system: coordinates and velocities.
1440 * \param[out] f Forces.
1441 * \param[out] enerd Energies.
1442 * \param[in] includePos True if coordinates are requested.
1443 * \param[in] includeVel True if velocities are requested.
1444 * \param[in] includeForce True if forces are requested.
1445 * \param[in] includeEnergy True if energies are requested.
1447 void openmm_copy_state(void *data,
1448 t_state *state, double *time,
1449 rvec f[], gmx_enerdata_t *enerd,
1450 gmx_bool includePos, gmx_bool includeVel, gmx_bool includeForce, gmx_bool includeEnergy)
1454 types += State::Positions;
1456 types += State::Velocities;
1458 types += State::Forces;
1460 types += State::Energy;
1465 State currentState = static_cast<OpenMMData*>(data)->context->getState(types);
1466 int numAtoms = static_cast<OpenMMData*>(data)->system->getNumParticles();
1469 for (int i = 0; i < numAtoms; i++)
1471 Vec3 x = currentState.getPositions()[i];
1472 state->x[i][0] = x[0];
1473 state->x[i][1] = x[1];
1474 state->x[i][2] = x[2];
1479 for (int i = 0; i < numAtoms; i++)
1481 Vec3 v = currentState.getVelocities()[i];
1482 state->v[i][0] = v[0];
1483 state->v[i][1] = v[1];
1484 state->v[i][2] = v[2];
1489 for (int i = 0; i < numAtoms; i++)
1491 Vec3 force = currentState.getForces()[i];
1499 int numConstraints = static_cast<OpenMMData*>(data)->system->getNumConstraints();
1500 int dof = 3*numAtoms-numConstraints;
1501 if (static_cast<OpenMMData*>(data)->removeCM)
1503 enerd->term[F_EPOT] = currentState.getPotentialEnergy();
1504 enerd->term[F_EKIN] = currentState.getKineticEnergy();
1505 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1506 enerd->term[F_TEMP] = 2.0*enerd->term[F_EKIN]/dof/BOLTZ;
1508 *time = currentState.getTime();
1510 catch (std::exception& e)
1512 gmx_fatal(FARGS, "OpenMM exception caught while retrieving state information: %s", e.what());