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
46 #include <types/simple.h>
61 #include "gmx_fatal.h"
66 #include "gpu_utils.h"
67 #include "mtop_util.h"
69 #include "openmm_wrapper.h"
71 using namespace OpenMM
;
74 #define MEM_ERR_MSG(str) \
75 "The %s-simulation GPU memory test detected errors. As memory errors would cause incorrect " \
76 "simulation results, gromacs has aborted execution.\n Make sure that your GPU's memory is not " \
77 "overclocked and that the device is properly cooled.\n", (str)
80 #define COMBRULE_CHK_TOL 1e-6
81 #define COMBRULE_SIGMA(sig1, sig2) (((sig1) + (sig2))/2)
82 #define COMBRULE_EPS(eps1, eps2) (sqrt((eps1) * (eps2)))
85 * \brief Convert string to integer type.
86 * \param[in] s String to convert from.
87 * \param[in] f Basefield format flag that takes any of the following I/O
88 * manipulators: dec, hex, oct.
89 * \param[out] t Destination variable to convert to.
92 static gmx_bool
from_string(T
& t
, const string
& s
, ios_base
& (*f
)(ios_base
&))
95 return !(iss
>> f
>> t
).fail();
99 * \brief Split string around a given delimiter.
100 * \param[in] s String to split.
101 * \param[in] delim Delimiter character.
102 * \returns Vector of strings found in \p s.
104 static vector
<string
> split(const string
&s
, char delim
)
106 vector
<string
> elems
;
109 while (getline(ss
, item
, delim
))
111 if (item
.length() != 0)
112 elems
.push_back(item
);
118 * \brief Split a string of the form "option=value" into "option" and "value" strings.
119 * This string corresponds to one option and the associated value from the option list
120 * in the mdrun -device argument.
122 * \param[in] s A string containing an "option=value" pair that needs to be split up.
123 * \param[out] opt The name of the option.
124 * \param[out] val Value of the option.
126 static void splitOptionValue(const string
&s
, string
&opt
, string
&val
)
128 size_t eqPos
= s
.find('=');
129 if (eqPos
!= string::npos
)
131 opt
= s
.substr(0, eqPos
);
132 if (eqPos
!= s
.length()) val
= s
.substr(eqPos
+1);
137 * \brief Compare two strings ignoring case.
138 * This function is in fact a wrapper around the gromacs function gmx_strncasecmp().
139 * \param[in] s1 String.
140 * \param[in] s2 String.
141 * \returns Similarly to the C function strncasecmp(), the return value is an
142 integer less than, equal to, or greater than 0 if \p s1 less than,
143 identical to, or greater than \p s2.
145 static gmx_bool
isStringEqNCase(const string
& s1
, const string
& s2
)
147 return (gmx_strncasecmp(s1
.c_str(), s2
.c_str(), max(s1
.length(), s2
.length())) == 0);
151 * \brief Convert string to upper case.
153 * \param[in] s String to convert to uppercase.
154 * \returns The given string converted to uppercase.
156 static string
toUpper(const string
&s
)
159 std::transform(stmp
.begin(), stmp
.end(), stmp
.begin(), static_cast < int(*)(int) > (toupper
));
164 \name Sizes of constant device option arrays GmxOpenMMPlatformOptions#platforms,
165 GmxOpenMMPlatformOptions#memtests, GmxOpenMMPlatformOptions#deviceid,
166 GmxOpenMMPlatformOptions#force_dev. */
168 #define SIZEOF_PLATFORMS 2 // 2
169 #define SIZEOF_MEMTESTS 3
170 #define SIZEOF_DEVICEIDS 1
171 #define SIZEOF_FORCE_DEV 2
173 #define SIZEOF_CHECK_COMBRULE 2
176 /*! Possible platform options in the mdrun -device option. */
177 static const char *devOptStrings
[] = { "platform", "deviceid", "memtest", "force-device", "check-combrule" };
179 /*! Enumerated platform options in the mdrun -device option. */
189 * \brief Class to extract and manage the platform options in the mdrun -device option.
192 class GmxOpenMMPlatformOptions
195 GmxOpenMMPlatformOptions(const char *opt
);
196 ~GmxOpenMMPlatformOptions() { options
.clear(); }
197 string
getOptionValue(const string
&opt
);
198 void remOption(const string
&opt
);
201 void setOption(const string
&opt
, const string
&val
);
203 map
<string
, string
> options
; /*!< Data structure to store the option (name, value) pairs. */
205 static const char * const platforms
[SIZEOF_PLATFORMS
]; /*!< Available OpenMM platforms; size #SIZEOF_PLATFORMS */
206 static const char * const memtests
[SIZEOF_MEMTESTS
]; /*!< Available types of memory tests, also valid
207 any positive integer >=15; size #SIZEOF_MEMTESTS */
208 static const char * const deviceid
[SIZEOF_DEVICEIDS
]; /*!< Possible values for deviceid option;
209 also valid any positive integer; size #SIZEOF_DEVICEIDS */
210 static const char * const force_dev
[SIZEOF_FORCE_DEV
]; /*!< Possible values for for force-device option;
211 size #SIZEOF_FORCE_DEV */
212 static const char * const check_combrule
[SIZEOF_CHECK_COMBRULE
]; /* XXX temporary debug feature to
213 turn off combination rule check */
216 const char * const GmxOpenMMPlatformOptions::platforms
[SIZEOF_PLATFORMS
]
217 = {"CUDA", "Reference"};
218 //= { "Reference", "CUDA" /*,"OpenCL"*/ };
219 const char * const GmxOpenMMPlatformOptions::memtests
[SIZEOF_MEMTESTS
]
220 = { "15", "full", "off" };
221 const char * const GmxOpenMMPlatformOptions::deviceid
[SIZEOF_DEVICEIDS
]
223 const char * const GmxOpenMMPlatformOptions::force_dev
[SIZEOF_FORCE_DEV
]
225 const char * const GmxOpenMMPlatformOptions::check_combrule
[SIZEOF_CHECK_COMBRULE
]
230 * Takes the option list, parses it, checks the options and their values for validity.
231 * When certain options are not provided by the user, as default value the first item
232 * of the respective constant array is taken (GmxOpenMMPlatformOptions#platforms,
233 * GmxOpenMMPlatformOptions#memtests, GmxOpenMMPlatformOptions#deviceid,
234 * GmxOpenMMPlatformOptions#force_dev).
235 * \param[in] optionString Option list part of the mdrun -device parameter.
237 GmxOpenMMPlatformOptions::GmxOpenMMPlatformOptions(const char *optionString
)
239 // set default values
240 setOption("platform", platforms
[0]);
241 setOption("memtest", memtests
[0]);
242 setOption("deviceid", deviceid
[0]);
243 setOption("force-device", force_dev
[0]);
244 setOption("check-combrule", check_combrule
[0]);
246 string
opt(optionString
);
248 // remove all whitespaces
249 opt
.erase(remove_if(opt
.begin(), opt
.end(), ::isspace
), opt
.end());
250 // tokenize around ","-s
251 vector
<string
> tokens
= split(opt
, ',');
253 for (vector
<string
>::iterator it
= tokens
.begin(); it
!= tokens
.end(); ++it
)
255 string opt
= "", val
= "";
256 splitOptionValue(*it
, opt
, val
);
258 if (isStringEqNCase(opt
, "platform"))
260 /* no check, this will fail if platform does not exist when we try to set it */
265 if (isStringEqNCase(opt
, "memtest"))
267 /* the value has to be an integer >15(s) or "full" OR "off" */
268 if (!isStringEqNCase(val
, "full") && !isStringEqNCase(val
, "off"))
271 if (!from_string
<int>(secs
, val
, std::dec
))
273 gmx_fatal(FARGS
, "Invalid value for option memtest option: \"%s\"!", val
.c_str());
277 gmx_fatal(FARGS
, "Incorrect value for memtest option (%d). "
278 "Memtest needs to run for at least 15s!", secs
);
285 if (isStringEqNCase(opt
, "deviceid"))
288 if (!from_string
<int>(id
, val
, std::dec
) )
290 gmx_fatal(FARGS
, "Invalid device id: \"%s\"!", val
.c_str());
296 if (isStringEqNCase(opt
, "force-device"))
299 if (!isStringEqNCase(val
, "yes") && !isStringEqNCase(val
, "no"))
301 gmx_fatal(FARGS
, "Invalid OpenMM force option: \"%s\"!", val
.c_str());
307 if (isStringEqNCase(opt
, "check-combrule"))
310 if (!isStringEqNCase(val
, "yes") && !isStringEqNCase(val
, "no"))
312 gmx_fatal(FARGS
, "Invalid OpenMM force option: \"%s\"!", val
.c_str());
319 // if we got till here something went wrong
320 gmx_fatal(FARGS
, "Invalid OpenMM platform option: \"%s\"!", (*it
).c_str());
326 * \brief Getter function.
327 * \param[in] opt Name of the option.
328 * \returns Returns the value associated to an option.
330 string
GmxOpenMMPlatformOptions::getOptionValue(const string
&opt
)
332 map
<string
, string
> :: const_iterator it
= options
.find(toUpper(opt
));
333 if (it
!= options
.end())
344 * \brief Setter function - private, only used from contructor.
345 * \param[in] opt Name of the option.
346 * \param[in] val Value for the option.
348 void GmxOpenMMPlatformOptions::setOption(const string
&opt
, const string
&val
)
350 options
[toUpper(opt
)] = val
;
354 * \brief Removes an option with its value from the map structure. If the option
355 * does not exist, returns without any action.
356 * \param[in] opt Name of the option.
358 void GmxOpenMMPlatformOptions::remOption(const string
&opt
)
360 options
.erase(toUpper(opt
));
364 * \brief Print option-value pairs to a file (debugging function).
366 void GmxOpenMMPlatformOptions::print()
368 cout
<< ">> Platform options: " << endl
369 << ">> platform = " << getOptionValue("platform") << endl
370 << ">> deviceID = " << getOptionValue("deviceid") << endl
371 << ">> memtest = " << getOptionValue("memtest") << endl
372 << ">> force-device = " << getOptionValue("force-device") << endl
;
376 * \brief Container for OpenMM related data structures that represent the bridge
377 * between the Gromacs data-structures and the OpenMM library and is but it's
378 * only passed through the API functions as void to disable direct access.
383 System
* system
; /*! The system to simulate. */
384 Context
* context
; /*! The OpenMM context in which the simulation is carried out. */
385 Integrator
* integrator
; /*! The integrator used in the simulation. */
386 gmx_bool removeCM
; /*! If \true remove venter of motion, false otherwise. */
387 GmxOpenMMPlatformOptions
*platformOpt
; /*! Platform options. */
391 * \brief Runs memtest on the GPU that has alreaby been initialized by OpenMM.
392 * \param[in] fplog Pointer to gromacs log file.
393 * \param[in] devId Device id of the GPU to run the test on.
394 Note: as OpenMM previously creates the context,for now this is always -1.
395 * \param[in] pre_post Contains either "Pre" or "Post" just to be able to differentiate in
396 * stdout messages/log between memtest carried out before and after simulation.
397 * \param[in] opt Pointer to platform options object.
399 static void runMemtest(FILE* fplog
, int devId
, const char* pre_post
, GmxOpenMMPlatformOptions
*opt
)
401 char strout_buf
[STRLEN
];
404 string s
= opt
->getOptionValue("memtest");
405 const char *test_type
= s
.c_str();
407 if (!gmx_strcasecmp(test_type
, "off"))
413 if (!gmx_strcasecmp(test_type
, "full"))
419 from_string
<int>(which_test
, test_type
, std::dec
);
425 gmx_fatal(FARGS
, "Amount of seconds for memetest is negative (%d). ", which_test
);
430 case 0: /* no memtest */
431 sprintf(strout_buf
, "%s-simulation GPU memtest skipped. Note, that faulty memory can cause "
432 "incorrect results!", pre_post
);
433 fprintf(fplog
, "%s\n", strout_buf
);
434 gmx_warning(strout_buf
);
437 case 1: /* quick memtest */
438 fprintf(fplog
, "%s-simulation %s GPU memtest in progress...\n", pre_post
, test_type
);
439 fprintf(stdout
, "\n%s-simulation %s GPU memtest in progress...", pre_post
, test_type
);
442 res
= do_quick_memtest(devId
);
445 case 2: /* full memtest */
446 fprintf(fplog
, "%s-simulation %s memtest in progress...\n", pre_post
, test_type
);
447 fprintf(stdout
, "\n%s-simulation %s memtest in progress...", pre_post
, test_type
);
450 res
= do_full_memtest(devId
);
453 default: /* timed memtest */
454 fprintf(fplog
, "%s-simulation ~%ds memtest in progress...\n", pre_post
, which_test
);
455 fprintf(stdout
, "\n%s-simulation ~%ds memtest in progress...", pre_post
, which_test
);
458 res
= do_timed_memtest(devId
, which_test
);
465 gmx_fatal(FARGS
, MEM_ERR_MSG(pre_post
));
469 fprintf(fplog
, "Memory test completed without errors.\n");
471 fprintf(stdout
, "done, no errors detected\n");
478 * \brief Convert Lennard-Jones parameters c12 and c6 to sigma and epsilon.
483 * \param[out] epsilon
485 static void convert_c_12_6(double c12
, double c6
, double *sigma
, double *epsilon
)
487 if (c12
== 0 && c6
== 0)
492 else if (c12
> 0 && c6
> 0)
494 *epsilon
= (c6
*c6
)/(4.0*c12
);
495 *sigma
= pow(c12
/c6
, 1.0/6.0);
499 gmx_fatal(FARGS
,"OpenMM only supports c6 > 0 and c12 > 0 or c6 = c12 = 0.");
504 * \brief Does gromacs option checking.
506 * Checks the gromacs mdp options for features unsupported in OpenMM, case in which
507 * interrupts the execution. It also warns the user about pecularities of OpenMM
509 * \param[in] fplog Gromacs log file pointer.
510 * \param[in] ir Gromacs input parameters, see ::t_inputrec
511 * \param[in] top Gromacs node local topology, \see gmx_localtop_t
512 * \param[in] state Gromacs state structure \see ::t_state
513 * \param[in] mdatoms Gromacs atom parameters, \see ::t_mdatoms
514 * \param[in] fr \see ::t_forcerec
515 * \param[in] state Gromacs systems state, \see ::t_state
517 static void checkGmxOptions(FILE* fplog
, GmxOpenMMPlatformOptions
*opt
,
518 t_inputrec
*ir
, gmx_localtop_t
*top
,
519 t_forcerec
*fr
, t_state
*state
)
521 char warn_buf
[STRLEN
];
524 double sigma_ij
=0, sigma_ji
=0, sigma_ii
=0, sigma_jj
=0, sigma_comb
;
525 double eps_ij
=0, eps_ji
=0, eps_ii
=0, eps_jj
=0, eps_comb
;
527 /* Abort if unsupported critical options are present */
532 gmx_warning( "OpenMM does not support leap-frog, will use velocity-verlet integrator.");
535 if ( (ir
->eI
!= eiMD
) &&
537 (ir
->eI
!= eiVVAK
) &&
542 gmx_fatal(FARGS
, "OpenMM supports only the following integrators: md/md-vv/md-vv-avek, sd/sd1, and bd.");
546 if ( !(ir
->coulombtype
== eelPME
||
547 EEL_RF(ir
->coulombtype
) ||
548 ir
->coulombtype
== eelRF
||
549 ir
->coulombtype
== eelEWALD
||
551 (ir
->coulombtype
== eelCUT
&& ir
->rcoulomb
== 0 && ir
->rvdw
== 0) ||
552 // we could have cut-off combined with GBSA (openmm will use RF)
553 ir
->implicit_solvent
== eisGBSA
) )
555 gmx_fatal(FARGS
,"OpenMM supports only the following methods for electrostatics: "
556 "NoCutoff (i.e. rcoulomb = rvdw = 0 ),Reaction-Field, Ewald or PME.");
559 if (EEL_RF(ir
->coulombtype
) && ir
->epsilon_rf
!= 0)
561 // openmm has epsilon_rf=inf hard-coded
562 gmx_warning("OpenMM will use a Reaction-Field epsilon of infinity instead of %g.",ir
->epsilon_rf
);
565 if (ir
->etc
!= etcNO
&&
570 gmx_warning("OpenMM supports only Andersen thermostat with the md/md-vv/md-vv-avek integrators.");
573 if (ir
->implicit_solvent
== eisGBSA
&&
574 ir
->gb_algorithm
!= egbOBC
)
576 gmx_warning("OpenMM does not support the specified algorithm for Generalized Born, will use OBC instead.");
579 if (ir
->opts
.ngtc
> 1)
580 gmx_fatal(FARGS
,"OpenMM does not support multiple temperature coupling groups.");
582 if (ir
->epc
!= epcNO
)
583 gmx_warning("OpenMM supports only Monte Carlo barostat for pressure coupling.");
585 if (ir
->opts
.annealing
[0])
586 gmx_fatal(FARGS
,"OpenMM does not support simulated annealing.");
588 if (top
->idef
.il
[F_CONSTR
].nr
> 0 && ir
->eConstrAlg
!= econtSHAKE
)
589 gmx_warning("OpenMM provides contraints as a combination "
590 "of SHAKE, SETTLE and CCMA. Accuracy is based on the SHAKE tolerance set "
591 "by the \"shake_tol\" option.");
594 gmx_fatal(FARGS
,"OpenMM does not support walls.");
596 if (ir
->ePull
!= epullNO
)
597 gmx_fatal(FARGS
,"OpenMM does not support pulling.");
599 /* check for interaction types */
600 for (i
= 0; i
< F_EPOT
; i
++)
602 if (!(i
== F_CONSTR
||
606 i
== F_UREY_BRADLEY
||
613 i
== F_GB12
|| /* The GB parameters are hardcoded both in */
614 i
== F_GB13
|| /* Gromacs and OpenMM */
616 top
->idef
.il
[i
].nr
> 0)
618 gmx_fatal(FARGS
, "OpenMM does not support (some) of the provided interaction "
619 "type(s) (%s) ", interaction_function
[i
].longname
);
623 if (ir
->efep
!= efepNO
)
624 gmx_fatal(FARGS
,"OpenMM does not support free energy calculations.");
626 if (ir
->opts
.ngacc
> 1)
627 gmx_fatal(FARGS
,"OpenMM does not support non-equilibrium MD (accelerated groups).");
629 if (IR_ELEC_FIELD(*ir
))
630 gmx_fatal(FARGS
,"OpenMM does not support electric fields.");
633 gmx_fatal(FARGS
,"OpenMM does not support QMMM calculations.");
635 if (ir
->rcoulomb
!= ir
->rvdw
)
636 gmx_fatal(FARGS
,"OpenMM uses a single cutoff for both Coulomb "
637 "and VdW interactions. Please set rcoulomb equal to rvdw.");
639 if (EEL_FULL(ir
->coulombtype
))
641 if (ir
->ewald_geometry
== eewg3DC
)
642 gmx_fatal(FARGS
,"OpenMM supports only Ewald 3D geometry.");
643 if (ir
->epsilon_surface
!= 0)
644 gmx_fatal(FARGS
,"OpenMM does not support dipole correction in Ewald summation.");
647 if (TRICLINIC(state
->box
))
649 gmx_fatal(FARGS
,"OpenMM does not support triclinic unit cells.");
652 /* XXX this is just debugging code to disable the combination rule check */
653 if ( isStringEqNCase(opt
->getOptionValue("check-combrule"), "yes") )
655 /* As OpenMM by default uses hardcoded combination rules
656 sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j)
657 we need to check whether the force field params obey this
658 and if not, we can't use this force field so we exit
659 grace-fatal-fully. */
660 real
*nbfp
= fr
->nbfp
;
664 fprintf(debug
, ">> Atom parameters: <<\n%10s%5s %5s %5s %5s COMB\n",
665 "", "i-j", "j-i", "i-i", "j-j");
667 /* loop over all i-j atom pairs and verify if
668 sigma_ij = sigma_ji = sigma_comb and eps_ij = eps_ji = eps_comb */
669 for (i
= 0; i
< natoms
; i
++)
672 c12
= C12(nbfp
, natoms
, i
, i
);
673 c6
= C6(nbfp
, natoms
, i
, i
);
674 convert_c_12_6(c12
, c6
, &sigma_ii
, &eps_ii
);
676 for (j
= 0; j
< i
; j
++)
679 c12
= C12(nbfp
, natoms
, i
, j
);
680 c6
= C6(nbfp
, natoms
, i
, j
);
681 convert_c_12_6(c12
, c6
, &sigma_ij
, &eps_ij
);
683 c12
= C12(nbfp
, natoms
, j
, i
);
684 c6
= C6(nbfp
, natoms
, j
, i
);
685 convert_c_12_6(c12
, c6
, &sigma_ji
, &eps_ji
);
687 c12
= C12(nbfp
, natoms
, j
, j
);
688 c6
= C6(nbfp
, natoms
, j
, j
);
689 convert_c_12_6(c12
, c6
, &sigma_jj
, &eps_jj
);
690 /* OpenMM hardcoded combination rules */
691 sigma_comb
= COMBRULE_SIGMA(sigma_ii
, sigma_jj
);
692 eps_comb
= COMBRULE_EPS(eps_ii
, eps_jj
);
696 fprintf(debug
, "i=%-3d j=%-3d", i
, j
);
697 fprintf(debug
, "%-11s", "sigma");
698 fprintf(debug
, "%5.3f %5.3f %5.3f %5.3f %5.3f\n",
699 sigma_ij
, sigma_ji
, sigma_ii
, sigma_jj
, sigma_comb
);
700 fprintf(debug
, "%11s%-11s", "", "epsilon");
701 fprintf(debug
, "%5.3f %5.3f %5.3f %5.3f %5.3f\n",
702 eps_ij
, eps_ji
, eps_ii
, eps_jj
, eps_comb
);
705 /* check the values against the rule used by omm */
706 if((fabs(eps_ij
) > COMBRULE_CHK_TOL
&&
707 fabs(eps_ji
) > COMBRULE_CHK_TOL
) &&
708 (fabs(sigma_comb
- sigma_ij
) > COMBRULE_CHK_TOL
||
709 fabs(sigma_comb
- sigma_ji
) > COMBRULE_CHK_TOL
||
710 fabs(eps_comb
- eps_ij
) > COMBRULE_CHK_TOL
||
711 fabs(eps_comb
- eps_ji
) > COMBRULE_CHK_TOL
))
714 "The combination rules of the used force-field do not "
715 "match the one supported by OpenMM: "
716 "sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j). "
717 "Switch to a force-field that uses these rules in order to "
718 "simulate this system using OpenMM.\n");
722 if (debug
) { fprintf(debug
, ">><<\n\n"); }
724 /* if we got here, log that everything is fine */
727 fprintf(debug
, ">> The combination rule of the used force matches the one used by OpenMM.\n");
729 fprintf(fplog
, "The combination rule of the used force field matches the one used by OpenMM.\n");
731 } /* if (are we checking the combination rules) ... */
736 * \brief Initialize OpenMM, run sanity/consistency checks, and return a pointer to
739 * Various gromacs data structures are passed that contain the parameters, state and
740 * other porperties of the system to simulate. These serve as input for initializing
741 * OpenMM. Besides, a set of misc action are taken:
742 * - OpenMM plugins are loaded;
743 * - platform options in \p platformOptStr are parsed and checked;
744 * - Gromacs parameters are checked for OpenMM support and consistency;
745 * - after the OpenMM is initialized memtest executed in the same GPU context.
747 * \param[in] fplog Gromacs log file handler.
748 * \param[in] platformOptStr Platform option string.
749 * \param[in] ir The Gromacs input parameters, see ::t_inputrec
750 * \param[in] top_global Gromacs system toppology, \see ::gmx_mtop_t
751 * \param[in] top Gromacs node local topology, \see gmx_localtop_t
752 * \param[in] mdatoms Gromacs atom parameters, \see ::t_mdatoms
753 * \param[in] fr \see ::t_forcerec
754 * \param[in] state Gromacs systems state, \see ::t_state
755 * \returns Pointer to a
758 void* openmm_init(FILE *fplog
, const char *platformOptStr
,
760 gmx_mtop_t
*top_global
, gmx_localtop_t
*top
,
761 t_mdatoms
*mdatoms
, t_forcerec
*fr
, t_state
*state
)
764 char warn_buf
[STRLEN
];
765 static gmx_bool hasLoadedPlugins
= false;
766 string usedPluginDir
;
771 if (!hasLoadedPlugins
)
773 vector
<string
> loadedPlugins
;
774 /* Look for OpenMM plugins at various locations (listed in order of priority):
775 - on the path in OPENMM_PLUGIN_DIR environment variable if this is specified
776 - on the path in the OPENMM_PLUGIN_DIR macro that is set by the build script
777 - at the default location assumed by OpenMM
780 char *pluginDir
= getenv("OPENMM_PLUGIN_DIR");
782 /* no env var or empty */
783 if (pluginDir
!= NULL
&& *pluginDir
!= '\0')
785 loadedPlugins
= Platform::loadPluginsFromDirectory(pluginDir
);
786 if (!loadedPlugins
.empty())
788 hasLoadedPlugins
= true;
789 usedPluginDir
= pluginDir
;
793 gmx_fatal(FARGS
, "The directory provided in the OPENMM_PLUGIN_DIR environment variable "
794 "(%s) does not contain valid OpenMM plugins. Check your OpenMM installation!",
799 /* macro set at build time */
800 #ifdef OPENMM_PLUGIN_DIR
801 if (!hasLoadedPlugins
)
803 loadedPlugins
= Platform::loadPluginsFromDirectory(OPENMM_PLUGIN_DIR
);
804 if (!loadedPlugins
.empty())
806 hasLoadedPlugins
= true;
807 usedPluginDir
= OPENMM_PLUGIN_DIR
;
811 /* default loocation */
812 if (!hasLoadedPlugins
)
814 loadedPlugins
= Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
815 if (!loadedPlugins
.empty())
817 hasLoadedPlugins
= true;
818 usedPluginDir
= Platform::getDefaultPluginsDirectory();
822 /* if there are still no plugins loaded there won't be any */
823 if (!hasLoadedPlugins
)
825 gmx_fatal(FARGS
, "No OpenMM plugins were found! You can provide the"
826 " plugin directory in the OPENMM_PLUGIN_DIR environment variable.", pluginDir
);
829 fprintf(fplog
, "\nOpenMM plugins loaded from directory %s:\t", usedPluginDir
.c_str());
830 for (int i
= 0; i
< (int)loadedPlugins
.size(); i
++)
832 fprintf(fplog
, "%s, ", loadedPlugins
[i
].c_str());
834 fprintf(fplog
, "\n");
837 /* parse option string */
838 GmxOpenMMPlatformOptions
*opt
= new GmxOpenMMPlatformOptions(platformOptStr
);
839 devId
= atoi(opt
->getOptionValue("deviceid").c_str());
846 /* check wheter Gromacs options compatibility with OpenMM */
847 checkGmxOptions(fplog
, opt
, ir
, top
, fr
, state
);
849 /* Create the system. */
850 const t_idef
& idef
= top
->idef
;
851 const int numAtoms
= top_global
->natoms
;
852 const int numConstraints
= idef
.il
[F_CONSTR
].nr
/3;
853 const int numSettle
= idef
.il
[F_SETTLE
].nr
/2;
854 const int numBonds
= idef
.il
[F_BONDS
].nr
/3;
855 const int numHarmonic
= idef
.il
[F_HARMONIC
].nr
/3;
856 const int numUB
= idef
.il
[F_UREY_BRADLEY
].nr
/4;
857 const int numAngles
= idef
.il
[F_ANGLES
].nr
/4;
858 const int numPeriodic
= idef
.il
[F_PDIHS
].nr
/5;
859 const int numPeriodicImproper
= idef
.il
[F_PIDIHS
].nr
/5;
860 const int numRB
= idef
.il
[F_RBDIHS
].nr
/5;
861 const int numImproperDih
= idef
.il
[F_IDIHS
].nr
/5;
862 const int num14
= idef
.il
[F_LJ14
].nr
/3;
863 System
* sys
= new System();
865 sys
->addForce(new CMMotionRemover(ir
->nstcomm
));
867 /* Set bonded force field terms. */
870 * CUDA platform currently doesn't support more than one
871 * instance of a force object, so we pack all forces that
872 * use the same form into one.
875 const int* bondAtoms
= (int*) idef
.il
[F_BONDS
].iatoms
;
876 HarmonicBondForce
* bondForce
= new HarmonicBondForce();
877 sys
->addForce(bondForce
);
879 for (int i
= 0; i
< numBonds
; ++i
)
881 int type
= bondAtoms
[offset
++];
882 int atom1
= bondAtoms
[offset
++];
883 int atom2
= bondAtoms
[offset
++];
884 bondForce
->addBond(atom1
, atom2
,
885 idef
.iparams
[type
].harmonic
.rA
, idef
.iparams
[type
].harmonic
.krA
);
888 const int* harmonicAtoms
= (int*) idef
.il
[F_HARMONIC
].iatoms
;
890 for (int i
= 0; i
< numHarmonic
; ++i
)
892 int type
= harmonicAtoms
[offset
++];
893 int atom1
= harmonicAtoms
[offset
++];
894 int atom2
= harmonicAtoms
[offset
++];
895 bondForce
->addBond(atom1
, atom2
,
896 idef
.iparams
[type
].harmonic
.rA
, idef
.iparams
[type
].harmonic
.krA
);
899 /* Set the angle force field terms */
900 const int* angleAtoms
= (int*) idef
.il
[F_ANGLES
].iatoms
;
901 HarmonicAngleForce
* angleForce
= new HarmonicAngleForce();
902 sys
->addForce(angleForce
);
904 for (int i
= 0; i
< numAngles
; ++i
)
906 int type
= angleAtoms
[offset
++];
907 int atom1
= angleAtoms
[offset
++];
908 int atom2
= angleAtoms
[offset
++];
909 int atom3
= angleAtoms
[offset
++];
910 angleForce
->addAngle(atom1
, atom2
, atom3
,
911 idef
.iparams
[type
].harmonic
.rA
*M_PI
/180.0, idef
.iparams
[type
].harmonic
.krA
);
914 /* Urey-Bradley includes both the angle and bond potential for 1-3 interactions */
915 const int* ubAtoms
= (int*) idef
.il
[F_UREY_BRADLEY
].iatoms
;
916 /* HarmonicBondForce* ubBondForce = new HarmonicBondForce(); */
917 /* HarmonicAngleForce* ubAngleForce = new HarmonicAngleForce(); */
918 /* sys->addForce(ubBondForce); */
919 /* sys->addForce(ubAngleForce); */
921 for (int i
= 0; i
< numUB
; ++i
)
923 int type
= ubAtoms
[offset
++];
924 int atom1
= ubAtoms
[offset
++];
925 int atom2
= ubAtoms
[offset
++];
926 int atom3
= ubAtoms
[offset
++];
927 /* ubBondForce->addBond(atom1, atom3, */
928 bondForce
->addBond(atom1
, atom3
,
929 idef
.iparams
[type
].u_b
.r13A
, idef
.iparams
[type
].u_b
.kUBA
);
930 /* ubAngleForce->addAngle(atom1, atom2, atom3, */
931 angleForce
->addAngle(atom1
, atom2
, atom3
,
932 idef
.iparams
[type
].u_b
.thetaA
*M_PI
/180.0, idef
.iparams
[type
].u_b
.kthetaA
);
935 /* Set proper dihedral terms */
936 const int* periodicAtoms
= (int*) idef
.il
[F_PDIHS
].iatoms
;
937 PeriodicTorsionForce
* periodicForce
= new PeriodicTorsionForce();
938 sys
->addForce(periodicForce
);
940 for (int i
= 0; i
< numPeriodic
; ++i
)
942 int type
= periodicAtoms
[offset
++];
943 int atom1
= periodicAtoms
[offset
++];
944 int atom2
= periodicAtoms
[offset
++];
945 int atom3
= periodicAtoms
[offset
++];
946 int atom4
= periodicAtoms
[offset
++];
947 periodicForce
->addTorsion(atom1
, atom2
, atom3
, atom4
,
948 idef
.iparams
[type
].pdihs
.mult
,
949 idef
.iparams
[type
].pdihs
.phiA
*M_PI
/180.0,
950 idef
.iparams
[type
].pdihs
.cpA
);
953 /* Set improper dihedral terms that are represented by a periodic function (as in AMBER FF) */
954 const int* periodicImproperAtoms
= (int*) idef
.il
[F_PIDIHS
].iatoms
;
955 /* PeriodicTorsionForce* periodicImproperForce = new PeriodicTorsionForce(); */
956 /* sys->addForce(periodicImproperForce); */
958 for (int i
= 0; i
< numPeriodicImproper
; ++i
)
960 int type
= periodicImproperAtoms
[offset
++];
961 int atom1
= periodicImproperAtoms
[offset
++];
962 int atom2
= periodicImproperAtoms
[offset
++];
963 int atom3
= periodicImproperAtoms
[offset
++];
964 int atom4
= periodicImproperAtoms
[offset
++];
965 /* periodicImproperForce->addTorsion(atom1, atom2, atom3, atom4, */
966 periodicForce
->addTorsion(atom1
, atom2
, atom3
, atom4
,
967 idef
.iparams
[type
].pdihs
.mult
,
968 idef
.iparams
[type
].pdihs
.phiA
*M_PI
/180.0,
969 idef
.iparams
[type
].pdihs
.cpA
);
972 /* Ryckaert-Bellemans dihedrals */
973 const int* rbAtoms
= (int*) idef
.il
[F_RBDIHS
].iatoms
;
974 RBTorsionForce
* rbForce
= new RBTorsionForce();
975 sys
->addForce(rbForce
);
977 for (int i
= 0; i
< numRB
; ++i
)
979 int type
= rbAtoms
[offset
++];
980 int atom1
= rbAtoms
[offset
++];
981 int atom2
= rbAtoms
[offset
++];
982 int atom3
= rbAtoms
[offset
++];
983 int atom4
= rbAtoms
[offset
++];
984 rbForce
->addTorsion(atom1
, atom2
, atom3
, atom4
,
985 idef
.iparams
[type
].rbdihs
.rbcA
[0], idef
.iparams
[type
].rbdihs
.rbcA
[1],
986 idef
.iparams
[type
].rbdihs
.rbcA
[2], idef
.iparams
[type
].rbdihs
.rbcA
[3],
987 idef
.iparams
[type
].rbdihs
.rbcA
[4], idef
.iparams
[type
].rbdihs
.rbcA
[5]);
990 /* Set improper dihedral terms (as in CHARMM FF) */
991 const int* improperDihAtoms
= (int*) idef
.il
[F_IDIHS
].iatoms
;
992 CustomTorsionForce
* improperDihForce
= new CustomTorsionForce("2.0*k*asin(sin((theta-theta0)/2))^2");
993 sys
->addForce(improperDihForce
);
994 improperDihForce
->addPerTorsionParameter("k");
995 improperDihForce
->addPerTorsionParameter("theta0");
996 vector
<double> improperDihParameters(2);
998 for (int i
= 0; i
< numImproperDih
; ++i
)
1000 int type
= improperDihAtoms
[offset
++];
1001 int atom1
= improperDihAtoms
[offset
++];
1002 int atom2
= improperDihAtoms
[offset
++];
1003 int atom3
= improperDihAtoms
[offset
++];
1004 int atom4
= improperDihAtoms
[offset
++];
1005 improperDihParameters
[0] = idef
.iparams
[type
].harmonic
.krA
;
1006 improperDihParameters
[1] = idef
.iparams
[type
].harmonic
.rA
*M_PI
/180.0;
1007 improperDihForce
->addTorsion(atom1
, atom2
, atom3
, atom4
,
1008 improperDihParameters
);
1011 /* Set nonbonded parameters and masses. */
1012 int ntypes
= fr
->ntype
;
1013 int* types
= mdatoms
->typeA
;
1014 real
* nbfp
= fr
->nbfp
;
1015 real
* charges
= mdatoms
->chargeA
;
1016 real
* masses
= mdatoms
->massT
;
1017 NonbondedForce
* nonbondedForce
= new NonbondedForce();
1018 sys
->addForce(nonbondedForce
);
1023 if (ir
->rcoulomb
== 0)
1025 nonbondedForce
->setNonbondedMethod(NonbondedForce::NoCutoff
);
1029 nonbondedForce
->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic
);
1033 switch (ir
->coulombtype
)
1040 nonbondedForce
->setNonbondedMethod(NonbondedForce::CutoffPeriodic
);
1044 nonbondedForce
->setNonbondedMethod(NonbondedForce::Ewald
);
1048 nonbondedForce
->setNonbondedMethod(NonbondedForce::PME
);
1052 gmx_fatal(FARGS
,"Internal error: you should not see this message, it means that the"
1053 "electrosatics option check failed. Please report this error!");
1055 sys
->setDefaultPeriodicBoxVectors(Vec3(state
->box
[0][0], 0, 0),
1056 Vec3(0, state
->box
[1][1], 0), Vec3(0, 0, state
->box
[2][2]));
1057 nonbondedForce
->setCutoffDistance(ir
->rcoulomb
);
1061 gmx_fatal(FARGS
,"OpenMM supports only full periodic boundary conditions "
1062 "(pbc = xyz), or none (pbc = no).");
1066 /* Fix for PME and Ewald error tolerance
1068 * OpenMM uses approximate formulas to calculate the Ewald parameter:
1069 * alpha = (1.0/cutoff)*sqrt(-log(2.0*tolerlance));
1070 * and the grid spacing for PME:
1071 * gridX = ceil(2*alpha*box[0][0]/3*(pow(tol, 0.2)))
1072 * gridY = ceil(2*alpha*box[1][1]/3*(pow(tol, 0.2)));
1073 * gridZ = ceil(2*alpha*box[2][2]/3*(pow(tol, 0.2)));
1076 * If the default ewald_rtol=1e-5 is used we silently adjust the value to the
1077 * OpenMM default of 5e-4 otherwise a warning is issued about the action taken.
1080 double corr_ewald_rtol
= 50.0 * ir
->ewald_rtol
;
1081 if ((ir
->ePBC
== epbcXYZ
) &&
1082 (ir
->coulombtype
== eelEWALD
|| ir
->coulombtype
== eelPME
))
1086 fprintf(debug
, ">> ewald_rtol = %e (corrected = %e) \n",
1087 ir
->ewald_rtol
, corr_ewald_rtol
);
1090 if (fabs(ir
->ewald_rtol
- 1e-5) > 1e-10)
1092 gmx_warning("OpenMM uses the ewald_rtol parameter with approximate formulas "
1093 "to calculate the alpha and grid spacing parameters of the Ewald "
1094 "and PME methods. This tolerance need to be corrected in order to get "
1095 "settings close to the ones used in GROMACS. Although the internal correction "
1096 "should work for any reasonable value of ewald_rtol, using values other than "
1097 "the default 1e-5 might cause incorrect behavior.");
1099 if (corr_ewald_rtol
> 1)
1101 gmx_fatal(FARGS
, "The ewald_rtol accuracy term is >1 after the "
1102 "adjustment for OpenMM (%e)", corr_ewald_rtol
);
1105 nonbondedForce
->setEwaldErrorTolerance(corr_ewald_rtol
);
1108 for (int i
= 0; i
< numAtoms
; ++i
)
1110 double c12
= nbfp
[types
[i
]*2*ntypes
+types
[i
]*2+1];
1111 double c6
= nbfp
[types
[i
]*2*ntypes
+types
[i
]*2];
1112 double sigma
=0.0, epsilon
=0.0;
1113 convert_c_12_6(c12
, c6
, &sigma
, &epsilon
);
1114 nonbondedForce
->addParticle(charges
[i
], sigma
, epsilon
);
1115 sys
->addParticle(masses
[i
]);
1118 // Build a table of all exclusions.
1119 vector
<set
<int> > exclusions(numAtoms
);
1120 for (int i
= 0; i
< numAtoms
; i
++)
1122 int start
= top
->excls
.index
[i
];
1123 int end
= top
->excls
.index
[i
+1];
1124 for (int j
= start
; j
< end
; j
++)
1125 exclusions
[i
].insert(top
->excls
.a
[j
]);
1128 // Record the 1-4 interactions, and remove them from the list of exclusions.
1129 const int* nb14Atoms
= (int*) idef
.il
[F_LJ14
].iatoms
;
1131 for (int i
= 0; i
< num14
; ++i
)
1133 int type
= nb14Atoms
[offset
++];
1134 int atom1
= nb14Atoms
[offset
++];
1135 int atom2
= nb14Atoms
[offset
++];
1136 double sigma
=0, epsilon
=0;
1137 convert_c_12_6(idef
.iparams
[type
].lj14
.c12A
,
1138 idef
.iparams
[type
].lj14
.c6A
,
1140 nonbondedForce
->addException(atom1
, atom2
,
1141 fr
->fudgeQQ
*charges
[atom1
]*charges
[atom2
], sigma
, epsilon
);
1142 exclusions
[atom1
].erase(atom2
);
1143 exclusions
[atom2
].erase(atom1
);
1146 // Record exclusions.
1147 for (int i
= 0; i
< numAtoms
; i
++)
1149 for (set
<int>::const_iterator iter
= exclusions
[i
].begin(); iter
!= exclusions
[i
].end(); ++iter
)
1153 nonbondedForce
->addException(i
, *iter
, 0.0, 1.0, 0.0);
1158 // Add GBSA if needed.
1159 if (ir
->implicit_solvent
== eisGBSA
)
1161 gmx_warning("The OBC scale factors alpha, beta and gamma are hardcoded in OpenMM with the default Gromacs values.");
1162 t_atoms atoms
= gmx_mtop_global_atoms(top_global
);
1163 GBSAOBCForce
* gbsa
= new GBSAOBCForce();
1165 sys
->addForce(gbsa
);
1166 gbsa
->setSoluteDielectric(ir
->epsilon_r
);
1167 gbsa
->setSolventDielectric(ir
->gb_epsilon_solvent
);
1168 gbsa
->setCutoffDistance(nonbondedForce
->getCutoffDistance());
1169 if (nonbondedForce
->getNonbondedMethod() == NonbondedForce::NoCutoff
)
1170 gbsa
->setNonbondedMethod(GBSAOBCForce::NoCutoff
);
1171 else if (nonbondedForce
->getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic
)
1172 gbsa
->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic
);
1173 else if (nonbondedForce
->getNonbondedMethod() == NonbondedForce::CutoffPeriodic
)
1174 gbsa
->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic
);
1176 gmx_fatal(FARGS
,"OpenMM supports only Reaction-Field electrostatics with OBC/GBSA.");
1178 for (int i
= 0; i
< numAtoms
; ++i
)
1180 gbsa
->addParticle(charges
[i
],
1181 top_global
->atomtypes
.gb_radius
[atoms
.atom
[i
].type
],
1182 top_global
->atomtypes
.S_hct
[atoms
.atom
[i
].type
]);
1184 free_t_atoms(&atoms
, FALSE
);
1188 const int* constraintAtoms
= (int*) idef
.il
[F_CONSTR
].iatoms
;
1190 for (int i
= 0; i
< numConstraints
; ++i
)
1192 int type
= constraintAtoms
[offset
++];
1193 int atom1
= constraintAtoms
[offset
++];
1194 int atom2
= constraintAtoms
[offset
++];
1195 sys
->addConstraint(atom1
, atom2
, idef
.iparams
[type
].constr
.dA
);
1197 const int* settleAtoms
= (int*) idef
.il
[F_SETTLE
].iatoms
;
1199 for (int i
= 0; i
< numSettle
; ++i
)
1201 int type
= settleAtoms
[offset
++];
1202 int oxygen
= settleAtoms
[offset
++];
1203 sys
->addConstraint(oxygen
, oxygen
+1, idef
.iparams
[type
].settle
.doh
);
1204 sys
->addConstraint(oxygen
, oxygen
+2, idef
.iparams
[type
].settle
.doh
);
1205 sys
->addConstraint(oxygen
+1, oxygen
+2, idef
.iparams
[type
].settle
.dhh
);
1208 // Create an integrator for simulating the system.
1209 double friction
= (ir
->opts
.tau_t
[0] == 0.0 ? 0.0 : 1.0/ir
->opts
.tau_t
[0]);
1213 integ
= new BrownianIntegrator(ir
->opts
.ref_t
[0], friction
, ir
->delta_t
);
1214 static_cast<BrownianIntegrator
*>(integ
)->setRandomNumberSeed(ir
->ld_seed
);
1216 else if (EI_SD(ir
->eI
))
1218 integ
= new LangevinIntegrator(ir
->opts
.ref_t
[0], friction
, ir
->delta_t
);
1219 static_cast<LangevinIntegrator
*>(integ
)->setRandomNumberSeed(ir
->ld_seed
);
1223 integ
= new VerletIntegrator(ir
->delta_t
);
1224 if ( ir
->etc
!= etcNO
)
1226 AndersenThermostat
* thermostat
= new AndersenThermostat(ir
->opts
.ref_t
[0], friction
);
1227 sys
->addForce(thermostat
);
1231 // Add pressure coupling
1232 if (ir
->epc
!= epcNO
)
1234 // convert gromacs pressure tensor to a scalar
1235 double pressure
= (ir
->ref_p
[0][0] + ir
->ref_p
[1][1] + ir
->ref_p
[2][2]) / 3.0;
1236 int frequency
= int(ir
->tau_p
/ ir
->delta_t
); // update frequency in time steps
1237 if (frequency
< 1) frequency
= 1;
1238 double temperature
= ir
->opts
.ref_t
[0]; // in kelvin
1239 sys
->addForce(new MonteCarloBarostat(pressure
, temperature
, frequency
));
1242 integ
->setConstraintTolerance(ir
->shake_tol
);
1244 // Create a context and initialize it.
1245 Context
* context
= NULL
;
1248 OpenMM could automatically select the "best" GPU, however we're not't
1249 going to let it do that for now, as the current algorithm is very rudimentary
1250 and we anyway support only CUDA.
1251 if (platformOptStr == NULL || platformOptStr == "")
1253 context = new Context(*sys, *integ);
1258 /* which platform should we use */
1259 for (int i
= 0; i
< (int)Platform::getNumPlatforms() && context
== NULL
; i
++)
1261 if (isStringEqNCase(opt
->getOptionValue("platform"), Platform::getPlatform(i
).getName()))
1263 Platform
& platform
= Platform::getPlatform(i
);
1264 // set standard properties
1265 platform
.setPropertyDefaultValue("CudaDevice", opt
->getOptionValue("deviceid"));
1266 // TODO add extra properties
1267 context
= new Context(*sys
, *integ
, platform
);
1270 if (context
== NULL
)
1272 gmx_fatal(FARGS
, "The requested platform \"%s\" could not be found.",
1273 opt
->getOptionValue("platform").c_str());
1277 Platform
& platform
= context
->getPlatform();
1278 fprintf(fplog
, "Gromacs will use the OpenMM platform: %s\n", platform
.getName().c_str());
1280 const vector
<string
>& properties
= platform
.getPropertyNames();
1283 for (int i
= 0; i
< (int)properties
.size(); i
++)
1285 fprintf(debug
, ">> %s: %s\n", properties
[i
].c_str(),
1286 platform
.getPropertyValue(*context
, properties
[i
]).c_str());
1291 if (isStringEqNCase(opt
->getOptionValue("platform"), "CUDA"))
1294 if (!from_string
<int>(tmp
, platform
.getPropertyValue(*context
, "CudaDevice"), std::dec
))
1296 gmx_fatal(FARGS
, "Internal error: couldn't determine the device selected by OpenMM");
1300 /* For now this is just to double-check if OpenMM selected the GPU we wanted,
1301 but when we'll let OpenMM select the GPU automatically, it will query the deviceId.
1305 gmx_fatal(FARGS
, "Internal error: OpenMM is using device #%d"
1306 "while initialized for device #%d", tmp
, devId
);
1309 /* check GPU compatibility */
1310 char gpuname
[STRLEN
];
1311 devId
= atoi(opt
->getOptionValue("deviceid").c_str());
1312 if (!is_gmx_openmm_supported_gpu(-1, gpuname
))
1314 if (!gmx_strcasecmp(opt
->getOptionValue("force-device").c_str(), "yes"))
1316 sprintf(warn_buf
, "Non-supported GPU selected (#%d, %s), forced continuing."
1317 "Note, that the simulation can be slow or it migth even crash.",
1319 fprintf(fplog
, "%s\n", warn_buf
);
1320 gmx_warning(warn_buf
);
1324 gmx_fatal(FARGS
, "The selected GPU (#%d, %s) is not supported by Gromacs! "
1325 "Most probably you have a low-end GPU which would not perform well, "
1326 "or new hardware that has not been tested with the current release. "
1327 "If you still want to try using the device, use the force-device=yes option.",
1333 fprintf(fplog
, "Gromacs will run on the GPU #%d (%s).\n", devId
, gpuname
);
1338 if (isStringEqNCase(opt
->getOptionValue("platform"), "CUDA"))
1340 /* pre-simulation memtest */
1341 runMemtest(fplog
, -1, "Pre", opt
);
1344 vector
<Vec3
> pos(numAtoms
);
1345 vector
<Vec3
> vel(numAtoms
);
1346 for (int i
= 0; i
< numAtoms
; ++i
)
1348 pos
[i
] = Vec3(state
->x
[i
][0], state
->x
[i
][1], state
->x
[i
][2]);
1349 vel
[i
] = Vec3(state
->v
[i
][0], state
->v
[i
][1], state
->v
[i
][2]);
1351 context
->setPositions(pos
);
1352 context
->setVelocities(vel
);
1354 // Return a structure containing the system, integrator, and context.
1355 OpenMMData
* data
= new OpenMMData();
1357 data
->integrator
= integ
;
1358 data
->context
= context
;
1359 data
->removeCM
= (ir
->nstcomm
> 0);
1360 data
->platformOpt
= opt
;
1363 catch (std::exception
& e
)
1365 gmx_fatal(FARGS
, "OpenMM exception caught while initializating: %s", e
.what());
1367 return NULL
; /* just to avoid warnings */
1371 * \brief Integrate one step.
1373 * \param[in] data OpenMMData object created by openmm_init().
1375 void openmm_take_one_step(void* data
)
1377 // static int step = 0; printf("----> taking step #%d\n", step++);
1380 static_cast<OpenMMData
*>(data
)->integrator
->step(1);
1382 catch (std::exception
& e
)
1384 gmx_fatal(FARGS
, "OpenMM exception caught while taking a step: %s", e
.what());
1389 * \brief Integrate n steps.
1391 * \param[in] data OpenMMData object created by openmm_init().
1393 void openmm_take_steps(void* data
, int nstep
)
1397 static_cast<OpenMMData
*>(data
)->integrator
->step(nstep
);
1399 catch (std::exception
& e
)
1401 gmx_fatal(FARGS
, "OpenMM exception caught while taking a step: %s", e
.what());
1406 * \brief Clean up the data structures cretead for OpenMM.
1408 * \param[in] log Log file pointer.
1409 * \param[in] data OpenMMData object created by openmm_init().
1411 void openmm_cleanup(FILE* fplog
, void* data
)
1413 OpenMMData
* d
= static_cast<OpenMMData
*>(data
);
1415 if (isStringEqNCase(d
->platformOpt
->getOptionValue("platform"), "CUDA"))
1417 /* post-simulation memtest */
1418 runMemtest(fplog
, -1, "Post", d
->platformOpt
);
1421 delete d
->integrator
;
1423 delete d
->platformOpt
;
1428 * \brief Copy the current state information from OpenMM into the Gromacs data structures.
1430 * This function results in the requested proprties to be copied from the
1431 * GPU to host. As this represents a bottleneck, the frequency of pulling data
1432 * should be minimized.
1434 * \param[in] data OpenMMData object created by openmm_init().
1435 * \param[out] time Simulation time for which the state was created.
1436 * \param[out] state State of the system: coordinates and velocities.
1437 * \param[out] f Forces.
1438 * \param[out] enerd Energies.
1439 * \param[in] includePos True if coordinates are requested.
1440 * \param[in] includeVel True if velocities are requested.
1441 * \param[in] includeForce True if forces are requested.
1442 * \param[in] includeEnergy True if energies are requested.
1444 void openmm_copy_state(void *data
,
1445 t_state
*state
, double *time
,
1446 rvec f
[], gmx_enerdata_t
*enerd
,
1447 gmx_bool includePos
, gmx_bool includeVel
, gmx_bool includeForce
, gmx_bool includeEnergy
)
1451 types
+= State::Positions
;
1453 types
+= State::Velocities
;
1455 types
+= State::Forces
;
1457 types
+= State::Energy
;
1462 State currentState
= static_cast<OpenMMData
*>(data
)->context
->getState(types
);
1463 int numAtoms
= static_cast<OpenMMData
*>(data
)->system
->getNumParticles();
1466 for (int i
= 0; i
< numAtoms
; i
++)
1468 Vec3 x
= currentState
.getPositions()[i
];
1469 state
->x
[i
][0] = x
[0];
1470 state
->x
[i
][1] = x
[1];
1471 state
->x
[i
][2] = x
[2];
1476 for (int i
= 0; i
< numAtoms
; i
++)
1478 Vec3 v
= currentState
.getVelocities()[i
];
1479 state
->v
[i
][0] = v
[0];
1480 state
->v
[i
][1] = v
[1];
1481 state
->v
[i
][2] = v
[2];
1486 for (int i
= 0; i
< numAtoms
; i
++)
1488 Vec3 force
= currentState
.getForces()[i
];
1496 int numConstraints
= static_cast<OpenMMData
*>(data
)->system
->getNumConstraints();
1497 int dof
= 3*numAtoms
-numConstraints
;
1498 if (static_cast<OpenMMData
*>(data
)->removeCM
)
1500 enerd
->term
[F_EPOT
] = currentState
.getPotentialEnergy();
1501 enerd
->term
[F_EKIN
] = currentState
.getKineticEnergy();
1502 enerd
->term
[F_ETOT
] = enerd
->term
[F_EPOT
] + enerd
->term
[F_EKIN
];
1503 enerd
->term
[F_TEMP
] = 2.0*enerd
->term
[F_EKIN
]/dof
/BOLTZ
;
1505 *time
= currentState
.getTime();
1507 catch (std::exception
& e
)
1509 gmx_fatal(FARGS
, "OpenMM exception caught while retrieving state information: %s", e
.what());