Remove OpenMM
authorMark Abraham <mark.j.abraham@gmail.com>
Fri, 5 Jul 2013 21:51:40 +0000 (23:51 +0200)
committerGerrit Code Review <gerrit@gerrit.gromacs.org>
Tue, 20 Aug 2013 22:16:03 +0000 (00:16 +0200)
This code is not functional in release-4-6 branch, and nobody has
plans to fix it. Many things are already out of sync in the duplicated
code. This will get much worse as time passes, so simplifying the few
places that still touch the main code base has slight advantages. And
nobody will feel the need to clean up the duplicated code.

Change-Id: I06c2eb71042fd830fd888fbee1b5e3efcae0379d

16 files changed:
CMakeLists.txt
admin/installguide/installguide.tex
src/contrib/BuildMdrunOpenMM.cmake [deleted file]
src/contrib/CMakeLists.txt
src/contrib/FindOpenMM.cmake [deleted file]
src/contrib/md_openmm.c [deleted file]
src/contrib/md_openmm.h [deleted file]
src/contrib/mdrun_openmm.c [deleted file]
src/contrib/openmm_gpu_utils.cu [deleted file]
src/contrib/openmm_gpu_utils.h [deleted file]
src/contrib/openmm_wrapper.cpp [deleted file]
src/contrib/openmm_wrapper.h [deleted file]
src/contrib/runner_openmm.c [deleted file]
src/gromacs/mdlib/mdebin.c
src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_utils.cuh
src/programs/mdrun/CMakeLists.txt

index 0d5d0a59cdb56b9bc3be2db0026bf54314559a75..32f062ee513b78fe8f7144f04f419282aa105c36 100644 (file)
@@ -127,10 +127,6 @@ mark_as_advanced(GMX_FAHCORE)
 # decide on GPU settings based on user-settings and GPU/CUDA detection
 include(gmxManageGPU)
 
-# TODO: move OpenMM to contrib
-option(GMX_OPENMM "Accelerated execution on GPUs through the OpenMM library (rerun cmake after changing to see relevant options)" OFF)
-mark_as_advanced(GMX_OPENMM)
-
 include(gmxDetectAcceleration)
 if(NOT DEFINED GMX_CPU_ACCELERATION)
     if(CMAKE_CROSSCOMPILING)
index aabe66a079ac83947420e58b5b520e7ff57cd768..6881cbf2999405adacced4ce774249275088e664 100644 (file)
@@ -28,7 +28,6 @@
 \newcommand{\threadmpi}{ThreadMPI}
 \newcommand{\openmpi}{OpenMPI}
 \newcommand{\openmp}{OpenMP}
-\newcommand{\openmm}{OpenMM}
 \newcommand{\lammpi}{LAM/MPI}
 \newcommand{\mpich}{MPICH}
 \newcommand{\cmake}{CMake}
@@ -145,17 +144,6 @@ If you wish to use the excellent new native GPU support in \gromacs,
 version is strongly encouraged. \nvidia{} GPUs with at least \nvidia{} compute
 capability 2.0 are required, e.g. Fermi or Kepler cards.
 
-The GPU support from \gromacs{} version 4.5 using \openmm{}
-\url{https://simtk.org/home/openmm} is still contained in the code,
-but in the ``user contributions'' section (\verb+src/contrib+). You
-will need to set
-\verb+-DGMX_OPENMM=on -DGMX_GPU=off -DGMX_MPI=off
--DGMX_THREAD_MPI=off\+ in order to build it. It also requires \cuda{},
-and remains the only hardware-based acceleration available for
-implicit solvent simulations in \gromacs{} at the moment. However, the
-long-term plan is to enable this functionality in core Gromacs, and
-not have the OpenMM interface supported by the \gromacs team.
-
 If you wish to run in parallel on multiple machines across a network,
 you will need to have
 \begin{itemize}
diff --git a/src/contrib/BuildMdrunOpenMM.cmake b/src/contrib/BuildMdrunOpenMM.cmake
deleted file mode 100644 (file)
index 4f82c34..0000000
+++ /dev/null
@@ -1,64 +0,0 @@
-#
-# This file is part of the GROMACS molecular simulation package.
-#
-# Copyright (c) 2012, by the GROMACS development team, led by
-# David van der Spoel, Berk Hess, Erik Lindahl, and including many
-# others, as listed in the AUTHORS file in the top-level source
-# directory and at http://www.gromacs.org.
-#
-# GROMACS is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public License
-# as published by the Free Software Foundation; either version 2.1
-# of the License, or (at your option) any later version.
-#
-# GROMACS is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with GROMACS; if not, see
-# http://www.gnu.org/licenses, or write to the Free Software Foundation,
-# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
-#
-# If you want to redistribute modifications to GROMACS, please
-# consider that scientific software is very special. Version
-# control is crucial - bugs must be traceable. We will be happy to
-# consider code for inclusion in the official distribution, but
-# derived work must not be called official GROMACS. Details are found
-# in the README & COPYING files - if they are missing, get the
-# official version at http://www.gromacs.org.
-#
-# To help us fund GROMACS development, we humbly ask that you cite
-# the research papers on the package. Check out http://www.gromacs.org.
-#
-
-include_directories(${OpenMM_INCLUDE_DIR})
-link_directories(${OpenMM_LIBRARY_DIR}) 
-# with this define no evn.var. is needed with OPENMM_PLUGIN_DIR
-# if the same OpenMM installation is used for running and building 
-add_definitions( -DOPENMM_PLUGIN_DIR="${OpenMM_PLUGIN_DIR}" ) 
-file(TO_CMAKE_PATH ${OpenMM_PLUGIN_DIR} _path)
-add_library(openmm_api_wrapper STATIC ${CMAKE_SOURCE_DIR}/src/contrib/openmm_wrapper.cpp)
-target_link_libraries(openmm_api_wrapper ${OpenMM_LIBRARIES})
-
-list(REMOVE_ITEM MDRUN_SOURCES mdrun.c runner.c)
-list(APPEND MDRUN_SOURCES
-    ${CMAKE_SOURCE_DIR}/src/contrib/md_openmm.c
-    ${CMAKE_SOURCE_DIR}/src/contrib/mdrun_openmm.c
-    ${CMAKE_SOURCE_DIR}/src/contrib/runner_openmm.c
-    )
-
-# this is to circumvent the following MSVC error: 
-# warning LNK4098: defaultlib 'LIBCMT' conflicts with use of other libs
-# fatal error LNK1169: one or more multiply defined symbols found
-if(GMX_OPENMM AND MSVC)
-    set_target_properties(mdrun PROPERTIES LINK_FLAGS "/NODEFAULTLIB:LIBCMT")
-endif()
-
-include_directories(${CMAKE_SOURCE_DIR}/src/gmxlib/gpu_utils)
-
-set_source_files_properties(main.c PROPERTIES LANGUAGE CXX)
-if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
-    set_source_files_properties(main.c PROPERTIES COMPILE_FLAGS "-x c++")
-endif()
index 3dff35583d6196be25ed899c58fac3dcd42cca05..32cb46f44ff8b4ebe4c8e03ca2c1567733d309f9 100644 (file)
@@ -3,82 +3,6 @@ set(CONTRIB_PROGRAMS
 
 )
 
-# Uncomment the next line to build OpenMM:
-#option(GMX_OPENMM "Accelerated execution on GPUs through the OpenMM library (no longer supported" ON)
-
-# At run time, you may need to set the environment variable
-# OPENMM_PLUGIN_DIR=PATH_TO_GROMACS/openmm/lib/plugins
-# to make things work
-
-if(GMX_OPENMM)
-    if(GMX_GPU)
-        message(FATAL_ERROR "The OpenMM build is not compatible with the native GPU build")
-    endif()
-
-    # These won't actually do anything unless they precede the
-    # definition of these options elsewhere. However, they can't be
-    # FORCE-d either.
-    set(GMX_BINARY_SUFFIX "-openmm" CACHE STRING "Suffix to distinguish OpenMM build from normal")
-    set(GMX_LIBS_SUFFIX "_openmm" CACHE STRING "Suffix to distinguish OpenMM build from normal")
-
-#######################################################################
-# Check for options incompatible with OpenMM build                    #
-#######################################################################
-    # we'll use the built-in fft to avoid unnecessary dependencies
-    string(TOUPPER ${GMX_FFT_LIBRARY} GMX_FFT_LIBRARY)
-    if(NOT ${GMX_FFT_LIBRARY} STREQUAL "FFTPACK")
-        message(STATUS "No external FFT libraries needed for the OpenMM build, switching to fftpack!")
-        set(GMX_FFT_LIBRARY "fftpack" CACHE STRING 
-               "No external FFT libraries needed for the OpenMM build, switching to  fftpack!" FORCE)
-    endif()
-    if(GMX_MPI)
-        message(FATAL_ERROR "The OpenMM build is not compatible with MPI!")
-    endif(GMX_MPI)
-    if(GMX_THREAD_MPI)
-        message(STATUS "Thread-MPI not compatible with OpenMM, disabled!")
-        set(GMX_THREAD_MPI OFF CACHE BOOL
-               "Thread-MPI not compatible with OpenMM build, disabled!" FORCE)
-    endif(GMX_THREAD_MPI)
-    if(GMX_OPENMP)
-        message(STATUS "OpenMP multithreading not compatible with OpenMM, disabled")
-        set(GMX_OPENMP OFF CACHE BOOL
-            "OpenMP multithreading not compatible with OpenMM, disabled!" FORCE)
-    endif()
-    if(GMX_SOFTWARE_INVSQRT)
-        set(GMX_SOFTWARE_INVSQRT OFF CACHE STRING 
-                "The OpenMM build does not need GROMACS software 1/sqrt!" FORCE)
-    endif(GMX_SOFTWARE_INVSQRT)
-    string(TOUPPER ${GMX_CPU_ACCELERATION} GMX_CPU_ACCELERATION)
-    if(NOT GMX_CPU_ACCELERATION STREQUAL "NONE")
-        message(STATUS "Switching off CPU-based acceleration, the OpenMM build does not support/need any!")
-        set(GMX_CPU_ACCELERATION "None" CACHE STRING
-            "Switching off CPU-based acceleration, the OpenMM build does not support/need any!" FORCE)
-    endif()
-    if(GMX_FAHCORE)
-        message(FATAL_ERROR "The OpenMM build does not support FAH build!")
-    endif(GMX_FAHCORE)
-    if(GMX_DOUBLE)
-        message(FATAL_ERROR  "The OpenMM-build does not support double precision calculations!")
-    endif()
-
-    set(CUDA_BUILD_EMULATION OFF)
-    find_package(CUDA 3.1 REQUIRED)
-    add_definitions(-DGMX_OPENMM)
-    if(CMAKE_BUILD_TYPE STREQUAL "DEBUG")    
-        set(CUDA_VERBOSE_BUILD ON)
-    endif()
-
-    # mark as advanced the unused variables
-    mark_as_advanced(FORCE GMX_CPU_ACCELERATION GMX_MPI GMX_FFT_LIBRARY 
-        GMX_QMMM_PROGRAM GMX_THREAD_MPI GMX_DOUBLE)
-
-else(GMX_OPENMM)
-
-     mark_as_advanced(CLEAR GMX_CPU_ACCELERATION GMX_MPI GMX_FFT_LIBRARY 
-        GMX_QMMM_PROGRAM GMX_THREAD_MPI GMX_DOUBLE)
-
-endif(GMX_OPENMM)
-
 foreach(PROG ${CONTRIB_PROGRAMS})
         add_executable(${PROG} ${PROG}.c ${NGMX_COMMON_SOURCE})
        set_target_properties(${PROG} PROPERTIES OUTPUT_NAME "${PROG}${GMX_BINARY_SUFFIX}")
diff --git a/src/contrib/FindOpenMM.cmake b/src/contrib/FindOpenMM.cmake
deleted file mode 100644 (file)
index 541b2a4..0000000
+++ /dev/null
@@ -1,104 +0,0 @@
-# Find OpenMM library.
-#
-# Looks for the OpenMM libraries at the default (/usr/local) location 
-# or custom location found in the OPENMM_ROOT_DIR environment variable. 
-#
-# The script defines defines: 
-#  OpenMM_FOUND     
-#  OpenMM_ROOT_DIR
-#  OpenMM_INCLUDE_DIR
-#  OpenMM_LIBRARY_DIR
-#  OpenMM_LIBRARIES      
-#  OpenMM_LIBRARIES_D   - debug version of libraries 
-#  OpenMM_PLUGIN_DIR
-#
-
-# Author: Szilard Pall (pszilard@cbr.su.se)
-
-if(OpenMM_INCLUDE_DIR AND OpenMM_LIBRARY_DIR AND OpenMM_PLUGIN_DIR)
-    set(OpenMM_FIND_QUIETLY)
-endif()
-
-file(TO_CMAKE_PATH "$ENV{OPENMM_ROOT_DIR}" _env_OPENMM_ROOT_DIR)
-
-if(IS_DIRECTORY ${_env_OPENMM_ROOT_DIR})
-    set(OpenMM_ROOT_DIR "${_env_OPENMM_ROOT_DIR}" CACHE PATH "OpenMM installation directory" FORCE)
-else()
-    if(IS_DIRECTORY "/usr/local/openmm")
-        set(OpenMM_ROOT_DIR "/usr/local/openmm" CACHE PATH "OpenMM installation directory" FORCE)
-    endif()
-endif()
-
-find_library(OpenMM_LIBRARIES
-    NAMES OpenMM
-    PATHS "${OpenMM_ROOT_DIR}/lib"
-    CACHE STRING "OpenMM libraries")
-
-find_library(OpenMM_LIBRARIES_D
-    NAMES OpenMM_d
-    PATHS "${OpenMM_ROOT_DIR}/lib"
-    CACHE STRING "OpenMM debug libraries")
-
-if(OpenMM_LIBRARIES_D AND NOT OpenMM_LIBRARIES)
-    set(OpenMM_LIBRARIES ${OpenMM_LIBRARIES_D}
-        CACHE STRING "OpenMM libraries" FORCE)
-    message(WARNING " Only found debug versions of the OpenMM libraries!")
-endif()
-
-get_filename_component(OpenMM_LIBRARY_DIR 
-    ${OpenMM_LIBRARIES} 
-    PATH
-    CACHE STRING "OpenMM library path")
-
-find_path(OpenMM_INCLUDE_DIR 
-    NAMES OpenMM.h 
-    PATHS "${OpenMM_ROOT_DIR}/include" "${OpenMM_LIBRARY_DIR}/../include"
-    CACHE STRING "OpenMM include directory")    
-
-# if we did not manage to set the root dir at the beginning but found the 
-# libs then set the ${OpenMM_LIBRARY_DIR}/.. as root
-if(NOT IS_DIRECTORY ${OpenMM_ROOT_DIR})
-    if (IS_DIRECTORY "${OpenMM_LIBRARY_DIR}/..") # just double-checking
-        get_filename_component(OpenMM_ROOT_DIR 
-            "${OpenMM_LIBRARY_DIR}/.." 
-            ABSOLUTE)
-    endif()   
-endif()
-
-if(NOT IS_DIRECTORY ${OpenMM_ROOT_DIR})
-    message(FATAL_ERROR "Could not find OpenMM! Set the OPENMM_ROOT_DIR environment "
-    "variable to contain the path of the OpenMM installation.")
-endif()
-
-if(NOT IS_DIRECTORY ${OpenMM_LIBRARY_DIR})
-    message(FATAL_ERROR "Can't find OpenMM libraries. Check your OpenMM installation!")
-endif()
-
-# now we can be sure that we have the library dir
-if(IS_DIRECTORY "${OpenMM_LIBRARY_DIR}/plugins")
-    get_filename_component(OpenMM_PLUGIN_DIR
-        "${OpenMM_LIBRARY_DIR}/plugins"
-        ABSOLUTE)
-    set(OpenMM_PLUGIN_DIR ${OpenMM_PLUGIN_DIR} CACHE PATH "OpenMM plugins directory")
-else()
-    message(WARNING "Could not detect the OpenMM plugin directory at the default location (${OpenMM_LIBRARY_DIR}/plugins)."
-            "Check your OpenMM installation or set the OPENMM_PLUGIN_DIR environment variable!")
-endif()
-
-if(NOT OpenMM_INCLUDE_DIR)
-    message(FATAL_ERROR "Can't find OpenMM includes. Check your OpenMM installation!")
-endif()
-
-set(OpenMM_ROOT_DIR ${OpenMM_ROOT_DIR} CACHE PATH "OpenMM installation directory")
-
-include(FindPackageHandleStandardArgs)
-find_package_handle_standard_args(OpenMM DEFAULT_MSG 
-                                    OpenMM_ROOT_DIR
-                                    OpenMM_LIBRARIES 
-                                    OpenMM_LIBRARY_DIR 
-                                    OpenMM_INCLUDE_DIR)
-
-mark_as_advanced(OpenMM_INCLUDE_DIR
-    OpenMM_LIBRARIES
-    OpenMM_LIBRARIES_D
-    OpenMM_LIBRARY_DIR)
diff --git a/src/contrib/md_openmm.c b/src/contrib/md_openmm.c
deleted file mode 100644 (file)
index d5a15e7..0000000
+++ /dev/null
@@ -1,585 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2010, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
- * Copyright (c) 2012, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <signal.h>
-#include <stdlib.h>
-
-#include "typedefs.h"
-#include "smalloc.h"
-#include "sysstuff.h"
-#include "vec.h"
-#include "statutil.h"
-#include "vcm.h"
-#include "mdebin.h"
-#include "nrnb.h"
-#include "calcmu.h"
-#include "index.h"
-#include "vsite.h"
-#include "update.h"
-#include "ns.h"
-#include "trnio.h"
-#include "xtcio.h"
-#include "mdrun.h"
-#include "md_support.h"
-#include "confio.h"
-#include "network.h"
-#include "pull.h"
-#include "xvgr.h"
-#include "physics.h"
-#include "names.h"
-#include "xmdrun.h"
-#include "ionize.h"
-#include "disre.h"
-#include "orires.h"
-#include "pme.h"
-#include "mdatoms.h"
-#include "qmmm.h"
-#include "mpelogging.h"
-#include "domdec.h"
-#include "partdec.h"
-#include "topsort.h"
-#include "coulomb.h"
-#include "constr.h"
-#include "compute_io.h"
-#include "mvdata.h"
-#include "checkpoint.h"
-#include "mtop_util.h"
-#include "sighandler.h"
-#include "genborn.h"
-#include "string2.h"
-#include "copyrite.h"
-#include "membed.h"
-
-#ifdef GMX_THREAD_MPI
-#include "tmpi.h"
-#endif
-
-/* include even when OpenMM not used to force compilation of do_md_openmm */
-#include "openmm_wrapper.h"
-
-double do_md_openmm(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
-                    const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
-                    int nstglobalcomm,
-                    gmx_vsite_t *vsite,gmx_constr_t constr,
-                    int stepout,t_inputrec *ir,
-                    gmx_mtop_t *top_global,
-                    t_fcdata *fcd,
-                    t_state *state_global,
-                    t_mdatoms *mdatoms,
-                    t_nrnb *nrnb,gmx_wallcycle_t wcycle,
-                    gmx_edsam_t ed,t_forcerec *fr,
-                    int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
-                    gmx_membed_t membed,
-                    real cpt_period,real max_hours,
-                    const char *deviceOptions,
-                    unsigned long Flags,
-                    gmx_runtime_t *runtime)
-{
-    gmx_mdoutf_t *outf;
-    gmx_large_int_t step,step_rel;
-    double     run_time;
-    double     t,t0,lam0;
-    gmx_bool       bSimAnn,
-    bFirstStep,bStateFromTPX,bLastStep,bStartingFromCpt;
-    gmx_bool       bInitStep=TRUE;
-    gmx_bool       do_ene,do_log, do_verbose,
-    bX,bV,bF,bCPT;
-    tensor     force_vir,shake_vir,total_vir,pres;
-    int        i,m;
-    int        mdof_flags;
-    rvec       mu_tot;
-    t_vcm      *vcm;
-    int        nchkpt=1;
-    gmx_localtop_t *top;
-    t_mdebin *mdebin;
-    t_state    *state=NULL;
-    rvec       *f_global=NULL;
-    int        n_xtc=-1;
-    rvec       *x_xtc=NULL;
-    gmx_enerdata_t *enerd;
-    rvec       *f=NULL;
-    gmx_global_stat_t gstat;
-    gmx_update_t upd=NULL;
-    t_graph    *graph=NULL;
-    globsig_t   gs;
-
-    gmx_groups_t *groups;
-    gmx_ekindata_t *ekind, *ekind_save;
-    gmx_bool        bAppend;
-    int         a0,a1;
-    matrix      lastbox;
-    real        reset_counters=0,reset_counters_now=0;
-    char        sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
-    int         handled_stop_condition=gmx_stop_cond_none; 
-
-    const char *ommOptions = NULL;
-    void   *openmmData;
-
-#ifdef GMX_DOUBLE
-    /* Checks in cmake should prevent the compilation in double precision
-     * with OpenMM, but just to be sure we check here.
-     */
-    gmx_fatal(FARGS,"Compilation was performed in double precision, but OpenMM only supports single precision. If you want to use to OpenMM, compile in single precision.");
-#endif
-
-    bAppend  = (Flags & MD_APPENDFILES);
-    check_ir_old_tpx_versions(cr,fplog,ir,top_global);
-
-    groups = &top_global->groups;
-
-    /* Initial values */
-    init_md(fplog,cr,ir,oenv,&t,&t0,state_global->lambda,
-            &(state_global->fep_state),&lam0,
-            nrnb,top_global,&upd,
-            nfile,fnm,&outf,&mdebin,
-            force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
-
-    clear_mat(total_vir);
-    clear_mat(pres);
-    /* Energy terms and groups */
-    snew(enerd,1);
-    init_enerdata(top_global->groups.grps[egcENER].nr,ir->fepvals->n_lambda,
-                  enerd);
-    snew(f,top_global->natoms);
-
-    /* Kinetic energy data */
-    snew(ekind,1);
-    init_ekindata(fplog,top_global,&(ir->opts),ekind);
-    /* needed for iteration of constraints */
-    snew(ekind_save,1);
-    init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
-    /* Copy the cos acceleration to the groups struct */
-    ekind->cosacc.cos_accel = ir->cos_accel;
-
-    gstat = global_stat_init(ir);
-    debug_gmx();
-
-    {
-        double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
-        if ((io > 2000) && MASTER(cr))
-            fprintf(stderr,
-                    "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
-                    io);
-    }
-
-    top = gmx_mtop_generate_local_top(top_global,ir);
-
-    a0 = 0;
-    a1 = top_global->natoms;
-
-    state = partdec_init_local_state(cr,state_global);
-    f_global = f;
-
-    atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
-
-    if (vsite)
-    {
-        set_vsite_top(vsite,top,mdatoms,cr);
-    }
-
-    if (ir->ePBC != epbcNONE && !ir->bPeriodicMols)
-    {
-        graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
-    }
-
-    update_mdatoms(mdatoms,state->lambda[efptMASS]);
-
-    if (deviceOptions[0]=='\0')
-    {
-        /* empty options, which should default to OpenMM in this build */
-        ommOptions=deviceOptions;
-    }
-    else
-    {
-        if (gmx_strncasecmp(deviceOptions,"OpenMM",6)!=0)
-        {
-            gmx_fatal(FARGS, "This Gromacs version currently only works with OpenMM. Use -device \"OpenMM:<options>\"");
-        }
-        else
-        {
-            ommOptions=strchr(deviceOptions,':');
-            if (NULL!=ommOptions)
-            {
-                /* Increase the pointer to skip the colon */
-                ommOptions++;
-            }
-        }
-    }
-
-    openmmData = openmm_init(fplog, ommOptions, ir, top_global, top, mdatoms, fr, state);
-    please_cite(fplog,"Friedrichs2009");
-
-    if (MASTER(cr))
-    {
-        /* Update mdebin with energy history if appending to output files */
-        if ( Flags & MD_APPENDFILES )
-        {
-            restore_energyhistory_from_state(mdebin,&state_global->enerhist);
-        }
-        /* Set the initial energy history in state to zero by updating once */
-        update_energyhistory(&state_global->enerhist,mdebin);
-    }
-
-    if (constr)
-    {
-        set_constraints(constr,top,ir,mdatoms,cr);
-    }
-
-    if (!ir->bContinuation)
-    {
-        if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
-        {
-            /* Set the velocities of frozen particles to zero */
-            for (i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
-            {
-                for (m=0; m<DIM; m++)
-                {
-                    if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
-                    {
-                        state->v[i][m] = 0;
-                    }
-                }
-            }
-        }
-
-        if (constr)
-        {
-            /* Constrain the initial coordinates and velocities */
-            do_constrain_first(fplog,constr,ir,mdatoms,state,f,
-                               graph,cr,nrnb,fr,top,shake_vir);
-        }
-        if (vsite)
-        {
-            /* Construct the virtual sites for the initial configuration */
-            construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
-                             top->idef.iparams,top->idef.il,
-                             fr->ePBC,fr->bMolPBC,graph,cr,state->box);
-        }
-    }
-
-    debug_gmx();
-
-    if (MASTER(cr))
-    {
-        char tbuf[20];
-        fprintf(stderr,"starting mdrun '%s'\n",
-                *(top_global->name));
-        if (ir->nsteps >= 0)
-        {
-            sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
-        }
-        else
-        {
-            sprintf(tbuf,"%s","infinite");
-        }
-        if (ir->init_step > 0)
-        {
-            fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
-                    gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
-                    gmx_step_str(ir->init_step,sbuf2),
-                    ir->init_step*ir->delta_t);
-        }
-        else
-        {
-            fprintf(stderr,"%s steps, %s ps.\n",
-                    gmx_step_str(ir->nsteps,sbuf),tbuf);
-        }
-    }
-
-    fprintf(fplog,"\n");
-
-    /* Set and write start time */
-    runtime_start(runtime);
-    print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
-    wallcycle_start(wcycle,ewcRUN);
-    if (fplog)
-        fprintf(fplog,"\n");
-
-    /* safest point to do file checkpointing is here.  More general point would be immediately before integrator call */
-
-    debug_gmx();
-    /***********************************************************
-     *
-     *             Loop over MD steps
-     *
-     ************************************************************/
-
-    /* loop over MD steps or if rerunMD to end of input trajectory */
-    bFirstStep = TRUE;
-    /* Skip the first Nose-Hoover integration when we get the state from tpx */
-    bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
-    bInitStep = bFirstStep && bStateFromTPX;
-    bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
-    bLastStep = FALSE;
-
-    init_global_signals(&gs,cr,ir,repl_ex_nst);
-
-    step = ir->init_step;
-    step_rel = 0;
-
-    while (!bLastStep)
-    {
-        wallcycle_start(wcycle,ewcSTEP);
-
-        GMX_MPE_LOG(ev_timestep1);
-
-        bLastStep = (step_rel == ir->nsteps);
-        t = t0 + step*ir->delta_t;
-
-        if (gs.set[eglsSTOPCOND] != 0)
-        {
-            bLastStep = TRUE;
-        }
-
-        do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
-        do_verbose = bVerbose &&
-                     (step % stepout == 0 || bFirstStep || bLastStep);
-
-        if (MASTER(cr) && do_log)
-        {
-            print_ebin_header(fplog,step,t,state->lambda[efptFEP]);
-        }
-
-        clear_mat(force_vir);
-        GMX_MPE_LOG(ev_timestep2);
-
-        /* We write a checkpoint at this MD step when:
-         * either when we signalled through gs (in OpenMM NS works different),
-         * or at the last step (but not when we do not want confout),
-         * but never at the first step.
-         */
-        bCPT = ((gs.set[eglsCHKPT] ||
-                 (bLastStep && (Flags & MD_CONFOUT))) &&
-                step > ir->init_step );
-        if (bCPT)
-        {
-            gs.set[eglsCHKPT] = 0;
-        }
-
-        /* Now we have the energies and forces corresponding to the
-         * coordinates at time t. We must output all of this before
-         * the update.
-         * for RerunMD t is read from input trajectory
-         */
-        GMX_MPE_LOG(ev_output_start);
-
-        mdof_flags = 0;
-        if (do_per_step(step,ir->nstxout))
-        {
-            mdof_flags |= MDOF_X;
-        }
-        if (do_per_step(step,ir->nstvout))
-        {
-            mdof_flags |= MDOF_V;
-        }
-        if (do_per_step(step,ir->nstfout))
-        {
-            mdof_flags |= MDOF_F;
-        }
-        if (do_per_step(step,ir->nstxtcout))
-        {
-            mdof_flags |= MDOF_XTC;
-        }
-        if (bCPT)
-        {
-            mdof_flags |= MDOF_CPT;
-        };
-        do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
-
-        if (mdof_flags != 0 || do_ene || do_log)
-        {
-            wallcycle_start(wcycle,ewcTRAJ);
-            bF = (mdof_flags & MDOF_F);
-            bX = (mdof_flags & (MDOF_X | MDOF_XTC | MDOF_CPT));
-            bV = (mdof_flags & (MDOF_V | MDOF_CPT));
-
-            openmm_copy_state(openmmData, state, &t, f, enerd, bX, bV, bF, do_ene);
-
-            upd_mdebin(mdebin,FALSE,TRUE,
-                       t,mdatoms->tmass,enerd,state,ir->fepvals,ir->expandedvals,lastbox,
-                       shake_vir,force_vir,total_vir,pres,
-                       ekind,mu_tot,constr);
-            print_ebin(outf->fp_ene,do_ene,FALSE,FALSE,do_log?fplog:NULL,
-                       step,t,
-                       eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
-            write_traj(fplog,cr,outf,mdof_flags,top_global,
-                       step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
-            if (bCPT)
-            {
-                nchkpt++;
-                bCPT = FALSE;
-            }
-            debug_gmx();
-            if (bLastStep && step_rel == ir->nsteps &&
-                    (Flags & MD_CONFOUT) && MASTER(cr))
-            {
-                /* x and v have been collected in write_traj,
-                 * because a checkpoint file will always be written
-                 * at the last step.
-                 */
-                fprintf(stderr,"\nWriting final coordinates.\n");
-                if (ir->ePBC != epbcNONE && !ir->bPeriodicMols)
-                {
-                    /* Make molecules whole only for confout writing */
-                    do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
-                }
-                write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
-                                    *top_global->name,top_global,
-                                    state_global->x,state_global->v,
-                                    ir->ePBC,state->box);
-                debug_gmx();
-            }
-            wallcycle_stop(wcycle,ewcTRAJ);
-        }
-        GMX_MPE_LOG(ev_output_finish);
-
-
-        /* Determine the wallclock run time up till now */
-        run_time = gmx_gettime() - (double)runtime->real;
-
-        /* Check whether everything is still allright */
-        if (((int)gmx_get_stop_condition() > handled_stop_condition)
-#ifdef GMX_THREAD_MPI
-            && MASTER(cr)
-#endif
-            )
-        {
-           /* this is just make gs.sig compatible with the hack 
-               of sending signals around by MPI_Reduce with together with
-               other floats */
-            /* NOTE: this only works for serial code. For code that allows
-               MPI nodes to propagate their condition, see kernel/md.c*/
-            if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
-                gs.set[eglsSTOPCOND]=1;
-            if ( gmx_get_stop_condition() == gmx_stop_cond_next )
-                gs.set[eglsSTOPCOND]=1;
-            /* < 0 means stop at next step, > 0 means stop at next NS step */
-            if (fplog)
-            {
-                fprintf(fplog,
-                        "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
-                        gmx_get_signal_name(),
-                        gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
-                fflush(fplog);
-            }
-            fprintf(stderr,
-                    "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
-                    gmx_get_signal_name(),
-                    gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
-            fflush(stderr);
-            handled_stop_condition=(int)gmx_get_stop_condition();
-        }
-        else if (MASTER(cr) &&
-                 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
-                 gs.set[eglsSTOPCOND] == 0)
-        {
-            /* Signal to terminate the run */
-            gs.set[eglsSTOPCOND] = 1;
-            if (fplog)
-            {
-                fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
-            }
-            fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
-        }
-
-        /* checkpoints */
-        if (MASTER(cr) && (cpt_period >= 0 &&
-                           (cpt_period == 0 ||
-                            run_time >= nchkpt*cpt_period*60.0)) &&
-                gs.set[eglsCHKPT] == 0)
-        {
-            gs.set[eglsCHKPT] = 1;
-        }
-
-        /* Time for performance */
-        if (((step % stepout) == 0) || bLastStep)
-        {
-            runtime_upd_proc(runtime);
-        }
-
-        if (do_per_step(step,ir->nstlog))
-        {
-            if (fflush(fplog) != 0)
-            {
-                gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of disk space?");
-            }
-        }
-
-        /* Remaining runtime */
-        if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
-        {
-            print_time(stderr,runtime,step,ir,cr);
-        }
-
-        bFirstStep = FALSE;
-        bInitStep = FALSE;
-        bStartingFromCpt = FALSE;
-        step++;
-        step_rel++;
-
-        openmm_take_one_step(openmmData);
-    }
-    /* End of main MD loop */
-    debug_gmx();
-
-    /* Stop the time */
-    runtime_end(runtime);
-
-    if (MASTER(cr))
-    {
-        if (ir->nstcalcenergy > 0) 
-        {
-            print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
-                       eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
-        }
-    }
-
-    openmm_cleanup(fplog, openmmData);
-
-    done_mdoutf(outf);
-
-    debug_gmx();
-
-    runtime->nsteps_done = step_rel;
-
-    return 0;
-}
diff --git a/src/contrib/md_openmm.h b/src/contrib/md_openmm.h
deleted file mode 100644 (file)
index d177db7..0000000
+++ /dev/null
@@ -1,59 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2010, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
- * Copyright (c) 2012, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#ifndef _MD_OPENMM_H
-#define _MD_OPENMM_H
-
-double do_md_openmm(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
-             const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
-             int nstglobalcomm,
-             gmx_vsite_t *vsite,gmx_constr_t constr,
-             int stepout,t_inputrec *ir,
-             gmx_mtop_t *top_global,
-             t_fcdata *fcd,
-             t_state *state_global,
-             t_mdatoms *mdatoms,
-             t_nrnb *nrnb,gmx_wallcycle_t wcycle,
-             gmx_edsam_t ed,t_forcerec *fr,
-             int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
-             gmx_membed_t membed,
-             real cpt_period,real max_hours,
-             const char *deviceOptions,
-             unsigned long Flags,
-             gmx_runtime_t *runtime);
-
-#endif /* _MD_OPENMM_H */
diff --git a/src/contrib/mdrun_openmm.c b/src/contrib/mdrun_openmm.c
deleted file mode 100644 (file)
index a258666..0000000
+++ /dev/null
@@ -1,511 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
- * Copyright (c) 2012, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "typedefs.h"
-#include "macros.h"
-#include "copyrite.h"
-#include "main.h"
-#include "statutil.h"
-#include "smalloc.h"
-#include "futil.h"
-#include "smalloc.h"
-#include "edsam.h"
-#include "mdrun.h"
-#include "xmdrun.h"
-#include "checkpoint.h"
-#ifdef GMX_THREAD_MPI
-#include "thread_mpi.h"
-#endif
-
-/* afm stuf */
-#include "pull.h"
-
-int cmain(int argc,char *argv[])
-{
-  const char *desc[] = {
-    "This is an experimental release of GROMACS for accelerated",
-       "Molecular Dynamics simulations on GPU processors. Support is provided",
-       "by the OpenMM library (https://simtk.org/home/openmm).[PAR]",
-       "*Warning*[BR]",
-       "This release is targeted at developers and advanced users and",
-       "care should be taken before production use. The following should be",
-       "noted before using the program:[PAR]",
-       " * The current release runs only on modern nVidia GPU hardware with CUDA support.",
-       "Make sure that the necessary CUDA drivers and libraries for your operating system",
-       "are already installed. The CUDA SDK also should be installed in order to compile",
-       "the program from source (http://www.nvidia.com/object/cuda_home.html).[PAR]",
-       " * Multiple GPU cards are not supported.[PAR]",
-       " * Only a small subset of the GROMACS features and options are supported on the GPUs.",
-       "See below for a detailed list.[PAR]",
-       " * Consumer level GPU cards are known to often have problems with faulty memory.",
-       "It is recommended that a full memory check of the cards is done at least once",
-       "(for example, using the memtest=full option).",
-       "A partial memory check (for example, memtest=15) before and",
-       "after the simulation run would help spot",
-       "problems resulting from processor overheating.[PAR]",
-       " * The maximum size of the simulated systems depends on the available",
-       "GPU memory,for example, a GTX280 with 1GB memory has been tested with systems",
-       "of up to about 100,000 atoms.[PAR]",
-       " * In order to take a full advantage of the GPU platform features, many algorithms",
-       "have been implemented in a very different way than they are on the CPUs.",
-       "Therefore numercal correspondence between properties of the state of",
-       "simulated systems should not be expected. Moreover, the values will likely vary",
-       "when simulations are done on different GPU hardware.[PAR]",
-       " * Frequent retrieval of system state information such as",
-       "trajectory coordinates and energies can greatly influence the performance",
-       "of the program due to slow CPU<->GPU memory transfer speed.[PAR]",
-       " * MD algorithms are complex, and although the Gromacs code is highly tuned for them,",
-       "they often do not translate very well onto the streaming architetures.",
-       "Realistic expectations about the achievable speed-up from test with GTX280:",
-       "For small protein systems in implicit solvent using all-vs-all kernels the acceleration",
-       "can be as high as 20 times, but in most other setups involving cutoffs and PME the",
-       "acceleration is usually only ~4 times relative to a 3GHz CPU.[PAR]",
-       "Supported features:[PAR]",
-       " * Integrators: md/md-vv/md-vv-avek, sd/sd1 and bd.\n",
-       " * Long-range interactions (option coulombtype): Reaction-Field, Ewald, PME, and cut-off (for Implicit Solvent only)\n",
-       " * Temperature control: Supported only with the md/md-vv/md-vv-avek, sd/sd1 and bd integrators.\n",
-       " * Pressure control: Supported.\n",
-       " * Implicit solvent: Supported.\n",
-       "A detailed description can be found on the GROMACS website:\n",
-       "http://www.gromacs.org/gpu[PAR]",
-/* From the original mdrun documentaion */
-    "The [TT]mdrun[tt] program reads the run input file ([TT]-s[tt])",
-    "and distributes the topology over nodes if needed.",
-    "[TT]mdrun[tt] produces at least four output files.",
-    "A single log file ([TT]-g[tt]) is written, unless the option",
-    "[TT]-seppot[tt] is used, in which case each node writes a log file.",
-    "The trajectory file ([TT]-o[tt]), contains coordinates, velocities and",
-    "optionally forces.",
-    "The structure file ([TT]-c[tt]) contains the coordinates and",
-    "velocities of the last step.",
-    "The energy file ([TT]-e[tt]) contains energies, the temperature,",
-    "pressure, etc, a lot of these things are also printed in the log file.",
-    "Optionally coordinates can be written to a compressed trajectory file",
-    "([TT]-x[tt]).[PAR]",
-/* openmm specific information */
-       "Usage with OpenMM:[BR]",
-       "[TT]mdrun -device \"OpenMM:platform=Cuda,memtest=15,deviceid=0,force-device=no\"[tt][PAR]",
-       "Options:[PAR]",
-       "      [TT]platform[tt] = Cuda\t\t:\tThe only available value. OpenCL support will be available in future.\n",
-       "      [TT]memtest[tt] = 15\t\t:\tRun a partial, random GPU memory test for the given amount of seconds. A full test",
-       "(recommended!) can be run with \"memtest=full\". Memory testing can be disabled with \"memtest=off\".\n",
-       "      [TT]deviceid[tt] = 0\t\t:\tSpecify the target device when multiple cards are present.",
-       "Only one card can be used at any given time though.\n",
-       "      [TT]force-device[tt] = no\t\t:\tIf set to \"yes\" [TT]mdrun[tt]  will be forced to execute on",
-       "hardware that is not officially supported. GPU acceleration can also be achieved on older",
-       "but Cuda capable cards, although the simulation might be too slow, and the memory limits too strict.",
-  };
-  t_commrec    *cr;
-  t_filenm fnm[] = {
-    { efTPX, NULL,      NULL,       ffREAD },
-    { efTRN, "-o",      NULL,       ffWRITE },
-    { efXTC, "-x",      NULL,       ffOPTWR },
-    { efCPT, "-cpi",    NULL,       ffOPTRD },
-    { efCPT, "-cpo",    NULL,       ffOPTWR },
-    { efSTO, "-c",      "confout",  ffWRITE },
-    { efEDR, "-e",      "ener",     ffWRITE },
-    { efLOG, "-g",      "md",       ffWRITE },
-    { efXVG, "-dhdl",   "dhdl",     ffOPTWR },
-    { efXVG, "-field",  "field",    ffOPTWR },
-    { efXVG, "-table",  "table",    ffOPTRD },
-    { efXVG, "-tabletf", "tabletf",    ffOPTRD },
-    { efXVG, "-tablep", "tablep",   ffOPTRD },
-    { efXVG, "-tableb", "table",    ffOPTRD },
-    { efTRX, "-rerun",  "rerun",    ffOPTRD },
-    { efXVG, "-tpi",    "tpi",      ffOPTWR },
-    { efXVG, "-tpid",   "tpidist",  ffOPTWR },
-    { efEDI, "-ei",     "sam",      ffOPTRD },
-    { efXVG, "-eo",     "sam",      ffOPTWR },
-    { efGCT, "-j",      "wham",     ffOPTRD },
-    { efGCT, "-jo",     "bam",      ffOPTWR },
-    { efXVG, "-ffout",  "gct",      ffOPTWR },
-    { efXVG, "-devout", "deviatie", ffOPTWR },
-    { efXVG, "-runav",  "runaver",  ffOPTWR },
-    { efXVG, "-px",     "pullx",    ffOPTWR },
-    { efXVG, "-pf",     "pullf",    ffOPTWR },
-    { efXVG, "-ro",     "rotation", ffOPTWR },
-    { efLOG, "-ra",     "rotangles",ffOPTWR },
-    { efLOG, "-rs",     "rotslabs", ffOPTWR },
-    { efLOG, "-rt",     "rottorque",ffOPTWR },
-    { efMTX, "-mtx",    "nm",       ffOPTWR },
-    { efNDX, "-dn",     "dipole",   ffOPTWR },
-    { efRND, "-multidir",NULL,      ffOPTRDMULT},
-    { efDAT, "-membed", "membed",   ffOPTRD },
-    { efTOP, "-mp",     "membed",   ffOPTRD },
-    { efNDX, "-mn",     "membed",   ffOPTRD }
-  };
-#define NFILE asize(fnm)
-
-  /* Command line options ! */
-  gmx_bool bCart        = FALSE;
-  gmx_bool bPPPME       = FALSE;
-  gmx_bool bPartDec     = FALSE;
-  gmx_bool bDDBondCheck = TRUE;
-  gmx_bool bDDBondComm  = TRUE;
-  gmx_bool bTunePME     = TRUE;
-  gmx_bool bTestVerlet  = FALSE;
-  gmx_bool bVerbose     = FALSE;
-  gmx_bool bCompact     = TRUE;
-  gmx_bool bSepPot      = FALSE;
-  gmx_bool bRerunVSite  = FALSE;
-  gmx_bool bIonize      = FALSE;
-  gmx_bool bConfout     = TRUE;
-  gmx_bool bReproducible = FALSE;
-    
-  int  npme=-1;
-  int  nmultisim=0;
-  int  nstglobalcomm=-1;
-  int  repl_ex_nst=0;
-  int  repl_ex_seed=-1;
-  int  repl_ex_nex=0;
-  int  nstepout=100;
-  int  resetstep=-1;
-  int  nsteps=-2; /* the value -2 means that the mdp option will be used */
-  
-  rvec realddxyz={0,0,0};
-  const char *ddno_opt[ddnoNR+1] =
-    { NULL, "interleave", "pp_pme", "cartesian", NULL };
-  const char *dddlb_opt[] =
-    { NULL, "auto", "no", "yes", NULL };
-  const char   *thread_aff_opt[threadaffNR+1] =
-    { NULL, "auto", "no", "yes", NULL };
-  const char *nbpu_opt[] =
-    { NULL, "auto", "cpu", "gpu", "gpu_cpu", NULL };
-  real rdd=0.0,rconstr=0.0,dlb_scale=0.8,pforce=-1;
-  char *ddcsx=NULL,*ddcsy=NULL,*ddcsz=NULL;
-  real cpt_period=15.0,max_hours=-1;
-  gmx_bool bAppendFiles=TRUE;
-  gmx_bool bKeepAndNumCPT=FALSE;
-  gmx_bool bResetCountersHalfWay=FALSE;
-  output_env_t oenv=NULL;
-  const char *deviceOptions = "";
-
-  gmx_hw_opt_t hw_opt={0,0,0,0,TRUE,FALSE,0,NULL};
-
-  t_pargs pa[] = {
-
-    { "-pd",      FALSE, etBOOL,{&bPartDec},
-      "Use particle decompostion" },
-    { "-dd",      FALSE, etRVEC,{&realddxyz},
-      "Domain decomposition grid, 0 is optimize" },
-    { "-ddorder", FALSE, etENUM, {ddno_opt},
-      "DD node order" },
-    { "-npme",    FALSE, etINT, {&npme},
-      "Number of separate nodes to be used for PME, -1 is guess" },
-    { "-nt",      FALSE, etINT, {&hw_opt.nthreads_tot},
-      "Total number of threads to start (0 is guess)" },
-    { "-ntmpi",   FALSE, etINT, {&hw_opt.nthreads_tmpi},
-      "Number of thread-MPI threads to start (0 is guess)" },
-    { "-ntomp",   FALSE, etINT, {&hw_opt.nthreads_omp},
-      "Number of OpenMP threads per MPI process/thread to start (0 is guess)" },
-    { "-ntomp_pme", FALSE, etINT, {&hw_opt.nthreads_omp_pme},
-      "Number of OpenMP threads per MPI process/thread to start (0 is -ntomp)" },
-    { "-pin",     FALSE, etBOOL, {thread_aff_opt},
-      "Pin OpenMP threads to cores" },
-    { "-pinoffset", FALSE, etINT, {&hw_opt.core_pinning_offset},
-      "Core offset for pinning (for running multiple mdrun processes on a single physical node)" },
-    { "-pinstride", FALSE, etINT, {&hw_opt.core_pinning_stride},
-      "Pinning distance in logical cores for threads, use 0 to minimize the number of threads per physical core" },
-    { "-gpu_id",  FALSE, etSTR, {&hw_opt.gpu_id},
-      "List of GPU id's to use" },
-    { "-ddcheck", FALSE, etBOOL, {&bDDBondCheck},
-      "Check for all bonded interactions with DD" },
-    { "-ddbondcomm", FALSE, etBOOL, {&bDDBondComm},
-      "HIDDENUse special bonded atom communication when [TT]-rdd[tt] > cut-off" },
-    { "-rdd",     FALSE, etREAL, {&rdd},
-      "The maximum distance for bonded interactions with DD (nm), 0 is determine from initial coordinates" },
-    { "-rcon",    FALSE, etREAL, {&rconstr},
-      "Maximum distance for P-LINCS (nm), 0 is estimate" },
-    { "-dlb",     FALSE, etENUM, {dddlb_opt},
-      "Dynamic load balancing (with DD)" },
-    { "-dds",     FALSE, etREAL, {&dlb_scale},
-      "Minimum allowed dlb scaling of the DD cell size" },
-    { "-ddcsx",   FALSE, etSTR, {&ddcsx},
-      "HIDDENThe DD cell sizes in x" },
-    { "-ddcsy",   FALSE, etSTR, {&ddcsy},
-      "HIDDENThe DD cell sizes in y" },
-    { "-ddcsz",   FALSE, etSTR, {&ddcsz},
-      "HIDDENThe DD cell sizes in z" },
-    { "-gcom",    FALSE, etINT,{&nstglobalcomm},
-      "Global communication frequency" },
-    { "-nb",      FALSE, etENUM, {&nbpu_opt},
-      "Calculate non-bonded interactions on" },
-    { "-tunepme", FALSE, etBOOL, {&bTunePME},  
-      "Optimize PME load between PP/PME nodes or GPU/CPU" },
-    { "-testverlet", FALSE, etBOOL, {&bTestVerlet},
-      "Test the Verlet non-bonded scheme" },
-    { "-v",       FALSE, etBOOL,{&bVerbose},  
-      "Be loud and noisy" },
-    { "-compact", FALSE, etBOOL,{&bCompact},  
-      "Write a compact log file" },
-    { "-seppot",  FALSE, etBOOL, {&bSepPot},
-      "Write separate V and dVdl terms for each interaction type and node to the log file(s)" },
-    { "-pforce",  FALSE, etREAL, {&pforce},
-      "Print all forces larger than this (kJ/mol nm)" },
-    { "-reprod",  FALSE, etBOOL,{&bReproducible},  
-      "Try to avoid optimizations that affect binary reproducibility" },
-    { "-cpt",     FALSE, etREAL, {&cpt_period},
-      "Checkpoint interval (minutes)" },
-    { "-cpnum",   FALSE, etBOOL, {&bKeepAndNumCPT},
-      "Keep and number checkpoint files" },
-    { "-append",  FALSE, etBOOL, {&bAppendFiles},
-      "Append to previous output files when continuing from checkpoint instead of adding the simulation part number to all file names" },
-    { "-nsteps",  FALSE, etINT, {&nsteps},
-      "Run this number of steps, overrides .mdp file option" },
-    { "-maxh",   FALSE, etREAL, {&max_hours},
-      "Terminate after 0.99 times this time (hours)" },
-    { "-multi",   FALSE, etINT,{&nmultisim}, 
-      "Do multiple simulations in parallel" },
-    { "-replex",  FALSE, etINT, {&repl_ex_nst}, 
-      "Attempt replica exchange periodically with this period (steps)" },
-    { "-nex",  FALSE, etINT, {&repl_ex_nex},
-      "Number of random exchanges to carry out each exchange interval (N^3 is one suggestion).  -nex zero or not specified gives neighbor replica exchange." },
-    { "-reseed",  FALSE, etINT, {&repl_ex_seed}, 
-      "Seed for replica exchange, -1 is generate a seed" },
-    { "-rerunvsite", FALSE, etBOOL, {&bRerunVSite},
-      "HIDDENRecalculate virtual site coordinates with [TT]-rerun[tt]" },
-    { "-ionize",  FALSE, etBOOL,{&bIonize},
-      "Do a simulation including the effect of an X-Ray bombardment on your system" },
-    { "-confout", FALSE, etBOOL, {&bConfout},
-      "HIDDENWrite the last configuration with [TT]-c[tt] and force checkpointing at the last step" },
-    { "-stepout", FALSE, etINT, {&nstepout},
-      "HIDDENFrequency of writing the remaining runtime" },
-    { "-resetstep", FALSE, etINT, {&resetstep},
-      "HIDDENReset cycle counters after these many time steps" },
-    { "-resethway", FALSE, etBOOL, {&bResetCountersHalfWay},
-      "HIDDENReset the cycle counters after half the number of steps or halfway [TT]-maxh[tt]" },
-    { "-device",  FALSE, etSTR, {&deviceOptions},
-      "Device option string" }
-  };
-  gmx_edsam_t  ed;
-  unsigned long Flags, PCA_Flags;
-  ivec     ddxyz;
-  int      dd_node_order;
-  gmx_bool     bAddPart;
-  FILE     *fplog,*fptest;
-  int      sim_part,sim_part_fn;
-  const char *part_suffix=".part";
-  char     suffix[STRLEN];
-  int      rc;
-  char **multidir=NULL;
-
-
-  cr = init_par(&argc,&argv);
-
-  if (MASTER(cr))
-    CopyRight(stderr, argv[0]);
-
-  PCA_Flags = (PCA_CAN_SET_DEFFNM | (MASTER(cr) ? 0 : PCA_QUIET));
-  
-  /* Comment this in to do fexist calls only on master
-   * works not with rerun or tables at the moment
-   * also comment out the version of init_forcerec in md.c 
-   * with NULL instead of opt2fn
-   */
-  /*
-     if (!MASTER(cr))
-     {
-     PCA_Flags |= PCA_NOT_READ_NODE;
-     }
-     */
-
-  parse_common_args(&argc,argv,PCA_Flags, NFILE,fnm,asize(pa),pa,
-                    asize(desc),desc,0,NULL, &oenv);
-
-
-
-  /* we set these early because they might be used in init_multisystem() 
-     Note that there is the potential for npme>nnodes until the number of
-     threads is set later on, if there's thread parallelization. That shouldn't
-     lead to problems. */ 
-  dd_node_order = nenum(ddno_opt);
-  cr->npmenodes = npme;
-
-  /* now check the -multi and -multidir option */
-  if (opt2bSet("-multidir", NFILE, fnm))
-  {
-      int i;
-      if (nmultisim > 0)
-      {
-          gmx_fatal(FARGS, "mdrun -multi and -multidir options are mutually exclusive.");
-      }
-      nmultisim = opt2fns(&multidir, "-multidir", NFILE, fnm);
-  }
-
-
-  if (repl_ex_nst != 0 && nmultisim < 2)
-      gmx_fatal(FARGS,"Need at least two replicas for replica exchange (option -multi)");
-
-  if (repl_ex_nex < 0)
-      gmx_fatal(FARGS,"Replica exchange number of exchanges needs to be positive");
-
-  if (nmultisim > 1) {
-#ifndef GMX_THREAD_MPI
-    gmx_bool bParFn = (multidir == NULL);
-    init_multisystem(cr, nmultisim, multidir, NFILE, fnm, bParFn);
-#else
-    gmx_fatal(FARGS,"mdrun -multi is not supported with the thread library.Please compile GROMACS with MPI support");
-#endif
-  }
-
-  bAddPart = !bAppendFiles;
-
-  /* Check if there is ANY checkpoint file available */        
-  sim_part    = 1;
-  sim_part_fn = sim_part;
-  if (opt2bSet("-cpi",NFILE,fnm))
-  {
-      if (bSepPot && bAppendFiles)
-      {
-          gmx_fatal(FARGS,"Output file appending is not supported with -seppot");
-      }
-
-      bAppendFiles =
-                read_checkpoint_simulation_part(opt2fn_master("-cpi", NFILE,
-                                                              fnm,cr),
-                                                &sim_part_fn,NULL,cr,
-                                                bAppendFiles,NFILE,fnm,
-                                                part_suffix,&bAddPart);
-      if (sim_part_fn==0 && MASTER(cr))
-      {
-          fprintf(stdout,"No previous checkpoint file present, assuming this is a new run.\n");
-      }
-      else
-      {
-          sim_part = sim_part_fn + 1;
-      }
-
-      if (MULTISIM(cr) && MASTER(cr))
-      {
-          check_multi_int(stdout,cr->ms,sim_part,"simulation part", TRUE);
-      }
-  } 
-  else
-  {
-      bAppendFiles = FALSE;
-  }
-
-  if (!bAppendFiles)
-  {
-      sim_part_fn = sim_part;
-  }
-
-  if (bAddPart)
-  {
-      /* Rename all output files (except checkpoint files) */
-      /* create new part name first (zero-filled) */
-      sprintf(suffix,"%s%04d",part_suffix,sim_part_fn);
-
-      add_suffix_to_output_names(fnm,NFILE,suffix);
-      if (MASTER(cr))
-      {
-          fprintf(stdout,"Checkpoint file is from part %d, new output files will be suffixed '%s'.\n",sim_part-1,suffix);
-      }
-  }
-
-  Flags = opt2bSet("-rerun",NFILE,fnm) ? MD_RERUN : 0;
-  Flags = Flags | (bSepPot       ? MD_SEPPOT       : 0);
-  Flags = Flags | (bIonize       ? MD_IONIZE       : 0);
-  Flags = Flags | (bPartDec      ? MD_PARTDEC      : 0);
-  Flags = Flags | (bDDBondCheck  ? MD_DDBONDCHECK  : 0);
-  Flags = Flags | (bDDBondComm   ? MD_DDBONDCOMM   : 0);
-  Flags = Flags | (bTunePME      ? MD_TUNEPME      : 0);
-  Flags = Flags | (bTestVerlet   ? MD_TESTVERLET   : 0);
-  Flags = Flags | (bConfout      ? MD_CONFOUT      : 0);
-  Flags = Flags | (bRerunVSite   ? MD_RERUN_VSITE  : 0);
-  Flags = Flags | (bReproducible ? MD_REPRODUCIBLE : 0);
-  Flags = Flags | (bAppendFiles  ? MD_APPENDFILES  : 0); 
-  Flags = Flags | (opt2parg_bSet("-append", asize(pa),pa) ? MD_APPENDFILESSET : 0); 
-  Flags = Flags | (bKeepAndNumCPT ? MD_KEEPANDNUMCPT : 0); 
-  Flags = Flags | (sim_part>1    ? MD_STARTFROMCPT : 0); 
-  Flags = Flags | (bResetCountersHalfWay ? MD_RESETCOUNTERSHALFWAY : 0);
-
-
-  /* We postpone opening the log file if we are appending, so we can 
-     first truncate the old log file and append to the correct position 
-     there instead.  */
-  if ((MASTER(cr) || bSepPot) && !bAppendFiles) 
-  {
-      gmx_log_open(ftp2fn(efLOG,NFILE,fnm),cr,
-                   !bSepPot,Flags & MD_APPENDFILES,&fplog);
-      CopyRight(fplog,argv[0]);
-      please_cite(fplog,"Hess2008b");
-      please_cite(fplog,"Spoel2005a");
-      please_cite(fplog,"Lindahl2001a");
-      please_cite(fplog,"Berendsen95a");
-  }
-  else if (!MASTER(cr) && bSepPot)
-  {
-      gmx_log_open(ftp2fn(efLOG,NFILE,fnm),cr,!bSepPot,Flags,&fplog);
-  }
-  else
-  {
-      fplog = NULL;
-  }
-
-  ddxyz[XX] = (int)(realddxyz[XX] + 0.5);
-  ddxyz[YY] = (int)(realddxyz[YY] + 0.5);
-  ddxyz[ZZ] = (int)(realddxyz[ZZ] + 0.5);
-
-  rc = mdrunner(&hw_opt, fplog,cr,NFILE,fnm,oenv,bVerbose,bCompact,
-                nstglobalcomm, ddxyz,dd_node_order,rdd,rconstr,
-                dddlb_opt[0],dlb_scale,ddcsx,ddcsy,ddcsz,
-                nbpu_opt[0],
-                nsteps,nstepout,resetstep,
-                nmultisim,repl_ex_nst,repl_ex_nex,repl_ex_seed,
-                pforce, cpt_period,max_hours,deviceOptions,Flags);
-
-  gmx_finalize_par();
-
-  if (MULTIMASTER(cr)) {
-      gmx_thanx(stderr);
-  }
-
-  /* Log file has to be closed in mdrunner if we are appending to it 
-     (fplog not set here) */
-  if (MASTER(cr) && !bAppendFiles) 
-  {
-      gmx_log_close(fplog);
-  }
-
-  return rc;
-}
-
diff --git a/src/contrib/openmm_gpu_utils.cu b/src/contrib/openmm_gpu_utils.cu
deleted file mode 100644 (file)
index cf10dfb..0000000
+++ /dev/null
@@ -1,145 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2010,2012 The GROMACS development team,
- * check out http://www.gromacs.org for more information.
- * Copyright (c) 2012, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-
-#include <stdio.h>
-#include <stdlib.h>
-
-#include "smalloc.h"
-#include "string2.h"
-#include "types/hw_info.h"
-
-#include "openmm_gpu_utils.h"
-#include "../../src/gmxlib/cuda_tools/cudautils.cuh"
-
-// TODO put this list into an external file and include it so that the list is easily accessible
-/*! List of supported GPUs. */
-static const char * const SupportedGPUs[] = {
-    /* GT400 */
-    "Geforce GTX 480",
-    "Geforce GTX 470",
-    "Geforce GTX 465",
-    "Geforce GTX 460",
-
-    "Tesla C2070",
-    "Tesla C2050",
-    "Tesla S2070",
-    "Tesla S2050",
-    "Tesla M2070",
-    "Tesla M2050",
-
-    "Quadro 5000",
-    "Quadro 6000",
-
-    /* GT200 */
-    "Geforce GTX 295",
-    "Geforce GTX 285",
-    "Geforce GTX 280",
-    "Geforce GTX 275",
-    "Geforce GTX 260",
-    "GeForce GTS 250",
-    "GeForce GTS 150",
-
-    "GeForce GTX 285M",
-    "GeForce GTX 280M",
-
-    "Tesla S1070",
-    "Tesla C1060",
-    "Tesla M1060",
-
-    "Quadro FX 5800",
-    "Quadro FX 4800",
-    "Quadro CX",
-    "Quadro Plex 2200 D2",
-    "Quadro Plex 2200 S4",
-
-    /* G90 */
-    "GeForce 9800 G", /* GX2, GTX, GTX+, GT */
-    "GeForce 9800M GTX",
-
-    "Quadro FX 4700",
-    "Quadro Plex 2100 D4"
-};
-
-/*! Number of supported GPUs */
-#define NB_GPUS (sizeof(SupportedGPUs)/sizeof(SupportedGPUs[0]))
-
-FUNC_QUALIFIER
-gmx_bool is_gmx_openmm_supported_gpu(int dev_id, char *gpu_name) FUNC_TERM_INT
-
-/*! 
- * \brief Checks whether the GPU with the given name is supported in Gromacs-OpenMM.
- * 
- * \param[in] gpu_name  the name of the CUDA device
- * \returns             TRUE if the device is supported, otherwise FALSE
- */
-static bool is_gmx_openmm_supported_gpu_name(char *gpuName)
-{
-    size_t i;
-    for (i = 0; i < NB_GPUS; i++)
-    {
-        trim(gpuName);
-        if (gmx_strncasecmp(gpuName, SupportedGPUs[i], strlen(SupportedGPUs[i])) == 0)
-            return 1;
-    }
-    return 0;
-}
-
-/*! \brief Checks whether the GPU with the given device id is supported in Gromacs-OpenMM.
- *
- * \param[in] dev_id    the device id of the GPU or -1 if the device has already been selected
- * \param[out] gpu_name Set to contain the name of the CUDA device, if NULL passed, no device name is set. 
- * \returns             TRUE if the device is supported, otherwise FALSE
- * 
- */
-gmx_bool is_gmx_openmm_supported_gpu(int dev_id, char *gpu_name)
-{
-    cudaDeviceProp dev_prop;
-
-    if (debug) fprintf(debug, "Checking compatibility with device #%d, %s\n", dev_id, gpu_name);
-
-    if (do_sanity_checks(dev_id, &dev_prop) != 0)
-        return -1;
-
-    if (gpu_name != NULL)
-    { 
-        strcpy(gpu_name, dev_prop.name);
-    }
-    return is_gmx_openmm_supported_gpu_name(dev_prop.name);
-}
-
-
diff --git a/src/contrib/openmm_gpu_utils.h b/src/contrib/openmm_gpu_utils.h
deleted file mode 100644 (file)
index c5874d6..0000000
+++ /dev/null
@@ -1,64 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2010, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
- * Copyright (c) 2012, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-
-#ifndef _OPENMM_GPU_UTILS_H_
-#define _OPENMM_GPU_UTILS_H_
-
-#include "types/simple.h"
-#include "types/hw_info.h"
-
-#define FUNC_TERM_INT {return -1;}
-#define FUNC_TERM_VOID {}
-#define FUNC_QUALIFIER static
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-FUNC_QUALIFIER
-gmx_bool is_gmx_openmm_supported_gpu(int dev_id, char *gpu_name) FUNC_TERM_INT
-
-#ifdef __cplusplus
-}
-#endif
-
-#undef FUNC_TERM_INT
-#undef FUNC_TERM_VOID
-#undef FUNC_QUALIFIER
-
-#endif /* _OPENMM_GPU_UTILS_H_ */
diff --git a/src/contrib/openmm_wrapper.cpp b/src/contrib/openmm_wrapper.cpp
deleted file mode 100644 (file)
index 3343aba..0000000
+++ /dev/null
@@ -1,1514 +0,0 @@
-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
- *
- * 
- *                This source code is part of
- * 
- *                 G   R   O   M   A   C   S
- * 
- *          GROningen MAchine for Chemical Simulations
- * 
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2010, 2013, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- * 
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
- * 
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- * 
- * For more info, check our website at http://www.gromacs.org
- * 
- * And Hey:
- * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
- */
-
-/*
- * Note, that parts of this source code originate from the Simtk release 
- * of OpenMM accelerated Gromacs, for more details see: 
- * https://simtk.org/project/xml/downloads.xml?group_id=161#package_id600
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <types/simple.h>
-#include <cmath>
-#include <set>
-#include <iostream>
-#include <sstream>
-#include <fstream>
-#include <map>
-#include <vector>
-#include <cctype>
-#include <algorithm>
-
-using namespace std;
-
-#include "OpenMM.h"
-
-#include "gmx_fatal.h"
-#include "typedefs.h"
-#include "mdrun.h"
-#include "physics.h"
-#include "string2.h"
-#include "openmm_gpu_utils.h"
-#include "gpu_utils.h"
-#include "mtop_util.h"
-
-#include "openmm_wrapper.h"
-
-using namespace OpenMM;
-
-/*! \cond */
-#define MEM_ERR_MSG(str) \
-    "The %s-simulation GPU memory test detected errors. As memory errors would cause incorrect " \
-    "simulation results, gromacs has aborted execution.\n Make sure that your GPU's memory is not " \
-    "overclocked and that the device is properly cooled.\n", (str)
-/*! \endcond */
-
-#define COMBRULE_CHK_TOL            1e-6
-#define COMBRULE_SIGMA(sig1, sig2)  (((sig1) + (sig2))/2)
-#define COMBRULE_EPS(eps1, eps2)    (sqrt((eps1) * (eps2)))
-
-/*! 
- * \brief Convert string to integer type.
- * \param[in]  s    String to convert from.
- * \param[in]  f    Basefield format flag that takes any of the following I/O
- *                  manipulators: dec, hex, oct.
- * \param[out] t    Destination variable to convert to.
- */
-template <class T>
-static gmx_bool from_string(T& t, const string& s, ios_base& (*f)(ios_base&))
-{
-    istringstream iss(s);
-    return !(iss >> f >> t).fail();
-}
-
-/*!
- * \brief Split string around a given delimiter.
- * \param[in] s      String to split.
- * \param[in] delim  Delimiter character.
- * \returns          Vector of strings found in \p s.
- */
-static vector<string> split(const string &s, char delim)
-{
-    vector<string> elems;
-    stringstream ss(s);
-    string item;
-    while (getline(ss, item, delim))
-    {
-        if (item.length() != 0)
-            elems.push_back(item);
-    }
-    return elems;
-}
-
-/*!
- * \brief Split a string of the form "option=value" into "option" and "value" strings.
- * This string corresponds to one option and the associated value from the option list 
- * in the mdrun -device argument.
- *
- * \param[in]  s    A string containing an "option=value" pair that needs to be split up.
- * \param[out] opt  The name of the option.
- * \param[out] val  Value of the option. 
- */
-static void splitOptionValue(const string &s, string &opt, string &val)
-{
-    size_t eqPos = s.find('=');
-    if (eqPos != string::npos)
-    {
-        opt = s.substr(0, eqPos);
-        if (eqPos != s.length())  val = s.substr(eqPos+1);
-    }
-}
-
-/*!
- * \brief Compare two strings ignoring case.
- * This function is in fact a wrapper around the gromacs function gmx_strncasecmp().
- * \param[in] s1 String. 
- * \param[in] s2 String.
- * \returns      Similarly to the C function strncasecmp(), the return value is an  
-                 integer less than, equal to, or greater than 0 if \p s1 less than, 
-                 identical to, or greater than \p s2.
- */
-static gmx_bool isStringEqNCase(const string& s1, const string& s2)
-{
-    return (gmx_strncasecmp(s1.c_str(), s2.c_str(), max(s1.length(), s2.length())) == 0);
-}
-
-/*!
- * \brief Convert string to upper case.
- *
- * \param[in]  s    String to convert to uppercase.
- * \returns         The given string converted to uppercase.
- */
-static string toUpper(const string &s)
-{
-    string stmp(s);
-    std::transform(stmp.begin(), stmp.end(), stmp.begin(), static_cast < int(*)(int) > (toupper));
-    return stmp;
-}
-
-/*! 
-  \name Sizes of constant device option arrays GmxOpenMMPlatformOptions#platforms, 
-  GmxOpenMMPlatformOptions#memtests, GmxOpenMMPlatformOptions#deviceid, 
-  GmxOpenMMPlatformOptions#force_dev.  */
-/* {@ */
-#define SIZEOF_PLATFORMS    2  // 2
-#define SIZEOF_MEMTESTS     3 
-#define SIZEOF_DEVICEIDS    1 
-#define SIZEOF_FORCE_DEV    2 
-
-#define SIZEOF_CHECK_COMBRULE 2
-/* @} */
-
-/*! Possible platform options in the mdrun -device option. */
-static const char *devOptStrings[] = { "platform", "deviceid", "memtest", "force-device", "check-combrule" }; 
-
-/*! Enumerated platform options in the mdrun -device option. */
-enum devOpt
-{
-    PLATFORM     = 0,
-    DEVICEID     = 1,
-    MEMTEST      = 2,
-    FORCE_DEVICE = 3
-};
-
-/*!
- * \brief Class to extract and manage the platform options in the mdrun -device option.
- * 
- */
-class GmxOpenMMPlatformOptions
-{
-public:
-    GmxOpenMMPlatformOptions(const char *opt);
-    ~GmxOpenMMPlatformOptions() { options.clear(); }
-    string getOptionValue(const string &opt);
-    void remOption(const string &opt);
-    void print();
-private:
-    void setOption(const string &opt, const string &val);
-
-    map<string, string> options; /*!< Data structure to store the option (name, value) pairs. */
-
-    static const char * const platforms[SIZEOF_PLATFORMS];  /*!< Available OpenMM platforms; size #SIZEOF_PLATFORMS */
-    static const char * const memtests[SIZEOF_MEMTESTS];    /*!< Available types of memory tests, also valid 
-                                                                 any positive integer >=15; size #SIZEOF_MEMTESTS */
-    static const char * const deviceid[SIZEOF_DEVICEIDS];   /*!< Possible values for deviceid option; 
-                                                                 also valid any positive integer; size #SIZEOF_DEVICEIDS */
-    static const char * const force_dev[SIZEOF_FORCE_DEV];  /*!< Possible values for for force-device option; 
-                                                                 size #SIZEOF_FORCE_DEV */
-    static const char * const check_combrule[SIZEOF_CHECK_COMBRULE]; /* XXX temporary debug feature to 
-                                                                      turn off combination rule check */
-};
-
-const char * const GmxOpenMMPlatformOptions::platforms[SIZEOF_PLATFORMS]
-                    = {"CUDA", "Reference"};
-                    //= { "Reference", "CUDA" /*,"OpenCL"*/ };
-const char * const GmxOpenMMPlatformOptions::memtests[SIZEOF_MEMTESTS]
-                    = { "15", "full", "off" };
-const char * const GmxOpenMMPlatformOptions::deviceid[SIZEOF_DEVICEIDS]
-                    = { "0" };
-const char * const GmxOpenMMPlatformOptions::force_dev[SIZEOF_FORCE_DEV]
-                    = { "no", "yes" };
-const char * const GmxOpenMMPlatformOptions::check_combrule[SIZEOF_CHECK_COMBRULE] 
-                    = { "yes", "no" };
-
-/*!
- * \brief Contructor.
- * Takes the option list, parses it, checks the options and their values for validity.
- * When certain options are not provided by the user, as default value the first item  
- * of the respective constant array is taken (GmxOpenMMPlatformOptions#platforms, 
- * GmxOpenMMPlatformOptions#memtests, GmxOpenMMPlatformOptions#deviceid, 
- * GmxOpenMMPlatformOptions#force_dev). 
- * \param[in] optionString  Option list part of the mdrun -device parameter.
- */
-GmxOpenMMPlatformOptions::GmxOpenMMPlatformOptions(const char *optionString)
-{
-    // set default values
-    setOption("platform",       platforms[0]);
-    setOption("memtest",        memtests[0]);
-    setOption("deviceid",       deviceid[0]);
-    setOption("force-device",   force_dev[0]);
-    setOption("check-combrule", check_combrule[0]);
-
-    string opt(optionString);
-
-    // remove all whitespaces
-    opt.erase(remove_if(opt.begin(), opt.end(), ::isspace), opt.end());
-    // tokenize around ","-s
-    vector<string> tokens = split(opt, ',');
-
-    for (vector<string>::iterator it = tokens.begin(); it != tokens.end(); ++it)
-    {
-        string opt = "", val = "";
-        splitOptionValue(*it, opt, val);
-
-        if (isStringEqNCase(opt, "platform"))
-        {
-            /* no check, this will fail if platform does not exist when we try to set it */
-            setOption(opt, val);
-            continue;
-        }
-
-        if (isStringEqNCase(opt, "memtest"))
-        {
-            /* the value has to be an integer >15(s) or "full" OR "off" */
-            if (!isStringEqNCase(val, "full") && !isStringEqNCase(val, "off")) 
-            {
-                int secs;
-                if (!from_string<int>(secs, val, std::dec))
-                {
-                    gmx_fatal(FARGS, "Invalid value for option memtest option: \"%s\"!", val.c_str());
-                }
-                if (secs < 15)
-                {
-                    gmx_fatal(FARGS, "Incorrect value for memtest option (%d). "
-                            "Memtest needs to run for at least 15s!", secs);
-                }
-            }
-            setOption(opt, val);
-            continue;
-        }
-
-        if (isStringEqNCase(opt, "deviceid"))
-        {
-            int id;
-            if (!from_string<int>(id, val, std::dec) )
-            {
-                gmx_fatal(FARGS, "Invalid device id: \"%s\"!", val.c_str());
-            }
-            setOption(opt, val);
-            continue;
-        }
-
-        if (isStringEqNCase(opt, "force-device"))
-        {
-            /* */
-            if (!isStringEqNCase(val, "yes") && !isStringEqNCase(val, "no"))
-            {
-                gmx_fatal(FARGS, "Invalid OpenMM force option: \"%s\"!", val.c_str());
-            }
-            setOption(opt, val);
-            continue;
-        }
-
-        if (isStringEqNCase(opt, "check-combrule"))
-        {
-            /* */
-            if (!isStringEqNCase(val, "yes") && !isStringEqNCase(val, "no"))
-            {
-                gmx_fatal(FARGS, "Invalid OpenMM force option: \"%s\"!", val.c_str());
-            }
-            setOption(opt, val);
-            continue;
-        }
-
-
-        // if we got till here something went wrong
-        gmx_fatal(FARGS, "Invalid OpenMM platform option: \"%s\"!", (*it).c_str());
-    }
-}
-
-
-/*!
- * \brief Getter function.
- * \param[in] opt   Name of the option.
- * \returns         Returns the value associated to an option. 
- */
-string GmxOpenMMPlatformOptions::getOptionValue(const string &opt)
-{
-       map<string, string> :: const_iterator it = options.find(toUpper(opt));
-       if (it != options.end())
-    {
-               return it->second;
-    }
-    else
-    {
-        return NULL;
-    }
-}
-
-/*!
- * \brief Setter function - private, only used from contructor.
- * \param[in] opt   Name of the option.
- * \param[in] val   Value for the option. 
- */
-void GmxOpenMMPlatformOptions::setOption(const string &opt, const string &val)
-{
-    options[toUpper(opt)] = val;
-}
-
-/*!
- * \brief Removes an option with its value from the map structure. If the option 
- * does not exist, returns without any action.
- * \param[in] opt   Name of the option.
- */
-void GmxOpenMMPlatformOptions::remOption(const string &opt) 
-{ 
-    options.erase(toUpper(opt)); 
-}
-
-/*!
- * \brief Print option-value pairs to a file (debugging function). 
- */
-void GmxOpenMMPlatformOptions::print()
-{
-    cout << ">> Platform options: " << endl 
-         << ">> platform     = " << getOptionValue("platform") << endl
-         << ">> deviceID     = " << getOptionValue("deviceid") << endl
-         << ">> memtest      = " << getOptionValue("memtest") << endl
-         << ">> force-device = " << getOptionValue("force-device") << endl;
-}
-
-/*!
- * \brief Container for OpenMM related data structures that represent the bridge 
- *        between the Gromacs data-structures and the OpenMM library and is but it's 
- *        only passed through the API functions as void to disable direct access. 
- */
-class OpenMMData
-{
-public:
-    System* system;      /*! The system to simulate. */
-    Context* context;   /*! The OpenMM context in which the simulation is carried out. */
-    Integrator* integrator; /*! The integrator used in the simulation. */
-    gmx_bool removeCM;          /*! If \true remove venter of motion, false otherwise. */
-    GmxOpenMMPlatformOptions *platformOpt; /*! Platform options. */
-};
-
-/*!
- *  \brief Runs memtest on the GPU that has alreaby been initialized by OpenMM.
- *  \param[in] fplog    Pointer to gromacs log file.
- *  \param[in] devId    Device id of the GPU to run the test on. 
-                        Note: as OpenMM previously creates the context,for now this is always -1.
- *  \param[in] pre_post Contains either "Pre" or "Post" just to be able to differentiate in 
- *                      stdout messages/log between memtest carried out before and after simulation.
- *  \param[in] opt      Pointer to platform options object.
- */
-static void runMemtest(FILE* fplog, int devId, const char* pre_post, GmxOpenMMPlatformOptions *opt)
-{
-    char        strout_buf[STRLEN];
-    int         which_test;
-    int         res = 0;
-    string      s = opt->getOptionValue("memtest");
-    const char  *test_type = s.c_str();
-
-    if (!gmx_strcasecmp(test_type, "off"))
-    {
-        which_test = 0;
-    }
-    else
-    {
-        if (!gmx_strcasecmp(test_type, "full"))
-        {
-            which_test = 2;
-        }
-        else
-        {
-            from_string<int>(which_test, test_type, std::dec);
-        }
-    }
-
-    if (which_test < 0) 
-    {
-        gmx_fatal(FARGS, "Amount of seconds for memetest is negative (%d). ", which_test);
-    }
-
-    switch (which_test)
-    {
-        case 0: /* no memtest */
-            sprintf(strout_buf, "%s-simulation GPU memtest skipped. Note, that faulty memory can cause "
-                "incorrect results!", pre_post);
-            fprintf(fplog, "%s\n", strout_buf);
-            gmx_warning(strout_buf);
-            break; /* case 0 */
-
-        case 1: /* quick memtest */
-            fprintf(fplog,  "%s-simulation %s GPU memtest in progress...\n", pre_post, test_type);
-            fprintf(stdout, "\n%s-simulation %s GPU memtest in progress...", pre_post, test_type);
-            fflush(fplog);
-            fflush(stdout);
-            res = do_quick_memtest(devId);
-            break; /* case 1 */
-
-        case 2: /* full memtest */
-            fprintf(fplog,  "%s-simulation %s memtest in progress...\n", pre_post, test_type);
-            fprintf(stdout, "\n%s-simulation %s memtest in progress...", pre_post, test_type);
-            fflush(fplog);
-            fflush(stdout);
-            res = do_full_memtest(devId);
-            break; /* case 2 */
-
-        default: /* timed memtest */
-            fprintf(fplog,  "%s-simulation ~%ds memtest in progress...\n", pre_post, which_test);
-            fprintf(stdout, "\n%s-simulation ~%ds memtest in progress...", pre_post, which_test);
-            fflush(fplog);
-            fflush(stdout);
-            res = do_timed_memtest(devId, which_test);
-        }
-
-        if (which_test != 0)
-        {
-            if (res != 0)
-            {
-                gmx_fatal(FARGS, MEM_ERR_MSG(pre_post));
-            }
-            else
-            {
-                fprintf(fplog,  "Memory test completed without errors.\n");
-                fflush(fplog);
-                fprintf(stdout, "done, no errors detected\n");
-                fflush(stdout);           
-            }
-        }
-}
-
-/*!
- * \brief Convert Lennard-Jones parameters c12 and c6 to sigma and epsilon.
- * 
- * \param[in] c12
- * \param[in] c6
- * \param[out] sigma 
- * \param[out] epsilon
- */
-static void convert_c_12_6(double c12, double c6, double *sigma, double *epsilon)
-{
-    if (c12 == 0 && c6 == 0)
-    {
-        *epsilon    = 0.0;        
-        *sigma      = 1.0;
-    }
-    else if (c12 > 0 && c6 > 0)
-    {
-        *epsilon    = (c6*c6)/(4.0*c12);
-        *sigma      = pow(c12/c6, 1.0/6.0);
-    }
-    else 
-    {
-        gmx_fatal(FARGS,"OpenMM only supports c6 > 0 and c12 > 0 or c6 = c12 = 0.");
-    } 
-}
-
-/*!
- * \brief Does gromacs option checking.
- *
- * Checks the gromacs mdp options for features unsupported in OpenMM, case in which 
- * interrupts the execution. It also warns the user about pecularities of OpenMM 
- * implementations.
- * \param[in] fplog         Gromacs log file pointer.
- * \param[in] ir            Gromacs input parameters, see ::t_inputrec
- * \param[in] top           Gromacs node local topology, \see gmx_localtop_t
- * \param[in] state         Gromacs state structure \see ::t_state
- * \param[in] mdatoms       Gromacs atom parameters, \see ::t_mdatoms
- * \param[in] fr            \see ::t_forcerec
- * \param[in] state         Gromacs systems state, \see ::t_state
- */
-static void checkGmxOptions(FILE* fplog, GmxOpenMMPlatformOptions *opt,
-                            t_inputrec *ir, gmx_localtop_t *top,
-                            t_forcerec *fr, t_state *state)
-{
-    char    warn_buf[STRLEN];
-    int     i, j, natoms;
-    double  c6, c12;
-    double  sigma_ij=0, sigma_ji=0, sigma_ii=0, sigma_jj=0, sigma_comb;
-    double  eps_ij=0, eps_ji=0, eps_ii=0, eps_jj=0, eps_comb;
-
-    /* Abort if unsupported critical options are present */
-
-    /* Integrator */
-    if (ir->eI ==  eiMD)
-    {
-        gmx_warning( "OpenMM does not support leap-frog, will use velocity-verlet integrator.");
-    }
-
-    if (    (ir->eI !=  eiMD)   &&
-            (ir->eI !=  eiVV)   &&
-            (ir->eI !=  eiVVAK) &&
-            (ir->eI !=  eiSD1)  &&
-            (ir->eI !=  eiSD2)  &&
-            (ir->eI !=  eiBD) )
-    {
-        gmx_fatal(FARGS, "OpenMM supports only the following integrators: md/md-vv/md-vv-avek, sd/sd1, and bd.");
-    }
-
-    /* Electroctstics */
-    if (   !(ir->coulombtype == eelPME   ||
-             EEL_RF(ir->coulombtype)     ||
-             ir->coulombtype == eelRF    ||
-             ir->coulombtype == eelEWALD ||
-             // no-cutoff
-             (ir->coulombtype == eelCUT && ir->rcoulomb == 0 &&  ir->rvdw == 0) ||
-             // we could have cut-off combined with GBSA (openmm will use RF)
-             ir->implicit_solvent == eisGBSA)   )
-    {
-        gmx_fatal(FARGS,"OpenMM supports only the following methods for electrostatics: "
-                "NoCutoff (i.e. rcoulomb = rvdw = 0 ),Reaction-Field, Ewald or PME.");
-    }
-
-    if (EEL_RF(ir->coulombtype) && ir->epsilon_rf != 0)
-    {
-        // openmm has epsilon_rf=inf hard-coded
-        gmx_warning("OpenMM will use a Reaction-Field epsilon of infinity instead of %g.",ir->epsilon_rf);
-    }
-
-    if (ir->etc != etcNO &&
-        ir->eI  != eiSD1 &&
-        ir->eI  != eiSD2 &&
-        ir->eI  != eiBD )
-    {
-        gmx_warning("OpenMM supports only Andersen thermostat with the md/md-vv/md-vv-avek integrators.");
-    }
-
-    if (ir->implicit_solvent == eisGBSA &&
-        ir->gb_algorithm != egbOBC  )
-    {
-        gmx_warning("OpenMM does not support the specified algorithm for Generalized Born, will use OBC instead.");
-    }
-
-    if (ir->opts.ngtc > 1)
-        gmx_fatal(FARGS,"OpenMM does not support multiple temperature coupling groups.");
-
-    if (ir->epc != epcNO)
-        gmx_warning("OpenMM supports only Monte Carlo barostat for pressure coupling.");
-
-    if (ir->opts.annealing[0])
-        gmx_fatal(FARGS,"OpenMM does not support simulated annealing.");
-    
-    if (top->idef.il[F_CONSTR].nr > 0 && ir->eConstrAlg != econtSHAKE)
-        gmx_warning("OpenMM provides contraints as a combination "
-                    "of SHAKE, SETTLE and CCMA. Accuracy is based on the SHAKE tolerance set "
-                    "by the \"shake_tol\" option.");
-
-    if (ir->nwall != 0)
-        gmx_fatal(FARGS,"OpenMM does not support walls.");
-
-    if (ir->ePull != epullNO)
-        gmx_fatal(FARGS,"OpenMM does not support pulling.");
-
-    /* check for interaction types */
-    for (i = 0; i < F_EPOT; i++)
-    {
-        if (!(i == F_CONSTR ||
-            i == F_SETTLE   ||
-            i == F_BONDS    ||            
-            i == F_HARMONIC ||
-            i == F_UREY_BRADLEY ||
-            i == F_ANGLES   ||
-            i == F_PDIHS    ||
-            i == F_RBDIHS   ||
-            i == F_PIDIHS   ||
-            i == F_IDIHS    ||
-            i == F_LJ14     ||
-            i == F_GB12     || /* The GB parameters are hardcoded both in */
-            i == F_GB13     || /* Gromacs and OpenMM */
-            i == F_GB14   ) &&
-            top->idef.il[i].nr > 0)
-        {
-            gmx_fatal(FARGS, "OpenMM does not support (some) of the provided interaction " 
-                    "type(s) (%s) ", interaction_function[i].longname);
-        }
-    }
-
-    if (ir->efep != efepNO)
-        gmx_fatal(FARGS,"OpenMM does not support free energy calculations.");
-
-    if (ir->opts.ngacc > 1)
-        gmx_fatal(FARGS,"OpenMM does not support non-equilibrium MD (accelerated groups).");
-
-    if (IR_ELEC_FIELD(*ir))
-        gmx_fatal(FARGS,"OpenMM does not support electric fields.");
-
-    if (ir->bQMMM)
-        gmx_fatal(FARGS,"OpenMM does not support QMMM calculations.");
-
-    if (ir->rcoulomb != ir->rvdw)
-        gmx_fatal(FARGS,"OpenMM uses a single cutoff for both Coulomb "
-                  "and VdW interactions. Please set rcoulomb equal to rvdw.");
-    
-    if (EEL_FULL(ir->coulombtype))
-    {
-        if (ir->ewald_geometry == eewg3DC)
-            gmx_fatal(FARGS,"OpenMM supports only Ewald 3D geometry.");
-        if (ir->epsilon_surface != 0)
-            gmx_fatal(FARGS,"OpenMM does not support dipole correction in Ewald summation.");
-    }
-
-    if (TRICLINIC(state->box))        
-    {
-        gmx_fatal(FARGS,"OpenMM does not support triclinic unit cells.");
-    }
-
-    /* XXX this is just debugging code to disable the combination rule check */
-    if ( isStringEqNCase(opt->getOptionValue("check-combrule"), "yes") )
-    {
-    /* As OpenMM by default uses hardcoded combination rules 
-       sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j)
-       we need to check whether the force field params obey this 
-       and if not, we can't use this force field so we exit 
-       grace-fatal-fully. */
-    real *nbfp = fr->nbfp;
-    natoms = fr->ntype;
-    if (debug) 
-    {   
-        fprintf(debug, ">> Atom parameters: <<\n%10s%5s %5s %5s %5s COMB\n", 
-                "", "i-j", "j-i", "i-i", "j-j");
-    }
-    /* loop over all i-j atom pairs and verify if 
-       sigma_ij = sigma_ji = sigma_comb and eps_ij = eps_ji = eps_comb */
-    for (i = 0; i < natoms; i++)
-    {
-        /* i-i */
-        /* nbfp now includes the 6.0/12.0 prefactors to save flops in kernels */
-        c12 = C12(nbfp, natoms, i, i)/12.0;
-        c6  = C6(nbfp,  natoms, i, i)/6.0;
-        convert_c_12_6(c12, c6, &sigma_ii, &eps_ii);
-
-        for (j = 0; j < i; j++)
-        {
-            /* i-j */
-            c12 = C12(nbfp, natoms, i, j)/12.0;
-            c6  = C6(nbfp,  natoms, i, j)/6.0;
-            convert_c_12_6(c12, c6, &sigma_ij, &eps_ij);
-            /* j-i */
-            c12 = C12(nbfp, natoms, j, i)/12.0;
-            c6  = C6(nbfp,  natoms, j, i)/6.0;
-            convert_c_12_6(c12, c6, &sigma_ji, &eps_ji);
-            /* j-j */
-            c12 = C12(nbfp, natoms, j, j)/12.0;
-            c6  = C6(nbfp,  natoms, j, j)/6.0;
-            convert_c_12_6(c12, c6, &sigma_jj, &eps_jj);
-            /* OpenMM hardcoded combination rules */
-            sigma_comb = COMBRULE_SIGMA(sigma_ii, sigma_jj);
-            eps_comb = COMBRULE_EPS(eps_ii, eps_jj);
-  
-            if (debug)
-            {
-                fprintf(debug, "i=%-3d j=%-3d", i, j);
-                fprintf(debug, "%-11s", "sigma");
-                fprintf(debug, "%5.3f %5.3f %5.3f %5.3f %5.3f\n",  
-                        sigma_ij, sigma_ji, sigma_ii, sigma_jj, sigma_comb);
-                fprintf(debug, "%11s%-11s", "", "epsilon");
-                fprintf(debug, "%5.3f %5.3f %5.3f %5.3f %5.3f\n", 
-                        eps_ij, eps_ji, eps_ii, eps_jj, eps_comb);
-            }
-
-            /* check the values against the rule used by omm */
-            if((fabs(eps_ij) > COMBRULE_CHK_TOL && 
-                fabs(eps_ji) > COMBRULE_CHK_TOL) &&
-               (fabs(sigma_comb - sigma_ij) > COMBRULE_CHK_TOL ||
-               fabs(sigma_comb - sigma_ji) > COMBRULE_CHK_TOL ||
-               fabs(eps_comb - eps_ij) > COMBRULE_CHK_TOL ||
-               fabs(eps_comb - eps_ji) > COMBRULE_CHK_TOL ))
-            {
-                gmx_fatal(FARGS,
-                        "The combination rules of the used force-field do not "
-                        "match the one supported by OpenMM:  "
-                        "sigma_ij = (sigma_i + sigma_j)/2, eps_ij = sqrt(eps_i * eps_j). "
-                        "Switch to a force-field that uses these rules in order to "
-                        "simulate this system using OpenMM.\n");                        
-            }
-        }
-    }
-    if (debug) { fprintf(debug, ">><<\n\n"); }
-
-    /* if we got here, log that everything is fine */
-    if (debug)
-    {
-        fprintf(debug, ">> The combination rule of the used force matches the one used by OpenMM.\n");
-    }
-    fprintf(fplog, "The combination rule of the used force field matches the one used by OpenMM.\n");   
-
-    } /* if (are we checking the combination rules) ... */
-}
-
-
-/*!
- * \brief Initialize OpenMM, run sanity/consistency checks, and return a pointer to 
- * the OpenMMData.
- * 
- * Various gromacs data structures are passed that contain the parameters, state and 
- * other porperties of the system to simulate. These serve as input for initializing 
- * OpenMM. Besides, a set of misc action are taken:
- *  - OpenMM plugins are loaded;
- *  - platform options in \p platformOptStr are parsed and checked; 
- *  - Gromacs parameters are checked for OpenMM support and consistency;
- *  - after the OpenMM is initialized memtest executed in the same GPU context.
- * 
- * \param[in] fplog             Gromacs log file handler.
- * \param[in] platformOptStr    Platform option string. 
- * \param[in] ir                The Gromacs input parameters, see ::t_inputrec
- * \param[in] top_global        Gromacs system toppology, \see ::gmx_mtop_t
- * \param[in] top               Gromacs node local topology, \see gmx_localtop_t
- * \param[in] mdatoms           Gromacs atom parameters, \see ::t_mdatoms
- * \param[in] fr                \see ::t_forcerec
- * \param[in] state             Gromacs systems state, \see ::t_state
- * \returns                     Pointer to a 
- * 
- */
-void* openmm_init(FILE *fplog, const char *platformOptStr,
-                  t_inputrec *ir,
-                  gmx_mtop_t *top_global, gmx_localtop_t *top,
-                  t_mdatoms *mdatoms, t_forcerec *fr, t_state *state)
-{
-
-    char warn_buf[STRLEN];
-    static gmx_bool hasLoadedPlugins = false;
-    string usedPluginDir;
-    int devId;
-
-    try
-    {
-        if (!hasLoadedPlugins)
-        {
-            vector<string> loadedPlugins;
-            /*  Look for OpenMM plugins at various locations (listed in order of priority):
-                - on the path in OPENMM_PLUGIN_DIR environment variable if this is specified
-                - on the path in the OPENMM_PLUGIN_DIR macro that is set by the build script
-                - at the default location assumed by OpenMM
-            */
-            /* env var */
-            char *pluginDir = getenv("OPENMM_PLUGIN_DIR");
-            trim(pluginDir);
-            /* no env var or empty */
-            if (pluginDir != NULL && *pluginDir != '\0')
-            {
-                loadedPlugins = Platform::loadPluginsFromDirectory(pluginDir);
-                if (!loadedPlugins.empty())
-                {
-                    hasLoadedPlugins = true;
-                    usedPluginDir = pluginDir;
-                }
-                else
-                {
-                    gmx_fatal(FARGS, "The directory provided in the OPENMM_PLUGIN_DIR environment variable "
-                              "(%s) does not contain valid OpenMM plugins. Check your OpenMM installation!", 
-                              pluginDir);
-                }
-            }
-
-            /* macro set at build time  */
-#ifdef OPENMM_PLUGIN_DIR
-            if (!hasLoadedPlugins)
-            {
-                loadedPlugins = Platform::loadPluginsFromDirectory(OPENMM_PLUGIN_DIR);
-                if (!loadedPlugins.empty())
-                {
-                    hasLoadedPlugins = true;
-                    usedPluginDir = OPENMM_PLUGIN_DIR;
-                }
-            }
-#endif
-            /* default loocation */
-            if (!hasLoadedPlugins)
-            {
-                loadedPlugins = Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
-                if (!loadedPlugins.empty())
-                {
-                    hasLoadedPlugins = true;
-                    usedPluginDir = Platform::getDefaultPluginsDirectory();
-                }
-            }
-
-            /* if there are still no plugins loaded there won't be any */
-            if (!hasLoadedPlugins)
-            {
-                gmx_fatal(FARGS, "No OpenMM plugins were found! You can provide the"
-                          " plugin directory in the OPENMM_PLUGIN_DIR environment variable.", pluginDir);
-            }
-
-            fprintf(fplog, "\nOpenMM plugins loaded from directory %s:\t", usedPluginDir.c_str());
-            for (int i = 0; i < (int)loadedPlugins.size(); i++)
-            {
-                fprintf(fplog, "%s, ", loadedPlugins[i].c_str());
-            }
-            fprintf(fplog, "\n");
-        }
-
-        /* parse option string */
-        GmxOpenMMPlatformOptions *opt = new GmxOpenMMPlatformOptions(platformOptStr);
-        devId = atoi(opt->getOptionValue("deviceid").c_str());
-
-        if (debug)
-        {
-            opt->print();
-        }
-
-        /* check wheter Gromacs options compatibility with OpenMM */
-        checkGmxOptions(fplog, opt, ir, top, fr, state);
-
-        /* Create the system. */
-        const t_idef& idef = top->idef;
-        const int numAtoms = top_global->natoms;
-        const int numConstraints = idef.il[F_CONSTR].nr/3;
-        const int numSettle = idef.il[F_SETTLE].nr/2;
-        const int numBonds = idef.il[F_BONDS].nr/3;
-        const int numHarmonic = idef.il[F_HARMONIC].nr/3;
-        const int numUB = idef.il[F_UREY_BRADLEY].nr/4;
-        const int numAngles = idef.il[F_ANGLES].nr/4;
-        const int numPeriodic = idef.il[F_PDIHS].nr/5;
-        const int numPeriodicImproper = idef.il[F_PIDIHS].nr/5;
-        const int numRB = idef.il[F_RBDIHS].nr/5;
-        const int numImproperDih = idef.il[F_IDIHS].nr/5;
-        const int num14 = idef.il[F_LJ14].nr/3;
-        System* sys = new System();
-        if (ir->nstcomm > 0)
-            sys->addForce(new CMMotionRemover(ir->nstcomm));
-
-        /* Set bonded force field terms. */
-
-               /* 
-                * CUDA platform currently doesn't support more than one
-                * instance of a force object, so we pack all forces that
-                * use the same form into one.
-               */
-
-        const int* bondAtoms = (int*) idef.il[F_BONDS].iatoms;
-        HarmonicBondForce* bondForce = new HarmonicBondForce();
-        sys->addForce(bondForce);
-        int offset = 0;
-        for (int i = 0; i < numBonds; ++i)
-        {
-            int type = bondAtoms[offset++];
-            int atom1 = bondAtoms[offset++];
-            int atom2 = bondAtoms[offset++];
-            bondForce->addBond(atom1, atom2,
-                               idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
-        }
-
-        const int* harmonicAtoms = (int*) idef.il[F_HARMONIC].iatoms;
-        offset = 0;
-        for (int i = 0; i < numHarmonic; ++i)
-        {
-            int type = harmonicAtoms[offset++];
-            int atom1 = harmonicAtoms[offset++];
-            int atom2 = harmonicAtoms[offset++];
-            bondForce->addBond(atom1, atom2,
-                               idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
-        }
-
-               /* Set the angle force field terms */
-        const int* angleAtoms = (int*) idef.il[F_ANGLES].iatoms;
-        HarmonicAngleForce* angleForce = new HarmonicAngleForce();
-        sys->addForce(angleForce);
-        offset = 0;
-        for (int i = 0; i < numAngles; ++i)
-        {
-            int type = angleAtoms[offset++];
-            int atom1 = angleAtoms[offset++];
-            int atom2 = angleAtoms[offset++];
-            int atom3 = angleAtoms[offset++];
-            angleForce->addAngle(atom1, atom2, atom3, 
-                    idef.iparams[type].harmonic.rA*M_PI/180.0, idef.iparams[type].harmonic.krA);
-        }
-
-        /* Urey-Bradley includes both the angle and bond potential for 1-3 interactions */
-        const int* ubAtoms = (int*) idef.il[F_UREY_BRADLEY].iatoms;
-               /* HarmonicBondForce* ubBondForce = new HarmonicBondForce(); */
-               /*  HarmonicAngleForce* ubAngleForce = new HarmonicAngleForce(); */
-        /* sys->addForce(ubBondForce); */
-        /* sys->addForce(ubAngleForce); */
-        offset = 0;
-        for (int i = 0; i < numUB; ++i)
-        {
-            int type = ubAtoms[offset++];
-            int atom1 = ubAtoms[offset++];
-            int atom2 = ubAtoms[offset++];
-            int atom3 = ubAtoms[offset++];
-            /* ubBondForce->addBond(atom1, atom3, */
-            bondForce->addBond(atom1, atom3,
-                               idef.iparams[type].u_b.r13A, idef.iparams[type].u_b.kUBA);
-            /* ubAngleForce->addAngle(atom1, atom2, atom3, */ 
-            angleForce->addAngle(atom1, atom2, atom3, 
-                    idef.iparams[type].u_b.thetaA*M_PI/180.0, idef.iparams[type].u_b.kthetaA);
-        }
-
-               /* Set proper dihedral terms */
-        const int* periodicAtoms = (int*) idef.il[F_PDIHS].iatoms;
-        PeriodicTorsionForce* periodicForce = new PeriodicTorsionForce();
-        sys->addForce(periodicForce);
-        offset = 0;
-        for (int i = 0; i < numPeriodic; ++i)
-        {
-            int type = periodicAtoms[offset++];
-            int atom1 = periodicAtoms[offset++];
-            int atom2 = periodicAtoms[offset++];
-            int atom3 = periodicAtoms[offset++];
-            int atom4 = periodicAtoms[offset++];
-            periodicForce->addTorsion(atom1, atom2, atom3, atom4,
-                                      idef.iparams[type].pdihs.mult,
-                                      idef.iparams[type].pdihs.phiA*M_PI/180.0, 
-                                      idef.iparams[type].pdihs.cpA);
-        }
-
-               /* Set improper dihedral terms that are represented by a periodic function (as in AMBER FF) */
-        const int* periodicImproperAtoms = (int*) idef.il[F_PIDIHS].iatoms;
-        /* PeriodicTorsionForce* periodicImproperForce = new PeriodicTorsionForce(); */
-        /* sys->addForce(periodicImproperForce); */
-        offset = 0;
-        for (int i = 0; i < numPeriodicImproper; ++i)
-        {
-            int type = periodicImproperAtoms[offset++];
-            int atom1 = periodicImproperAtoms[offset++];
-            int atom2 = periodicImproperAtoms[offset++];
-            int atom3 = periodicImproperAtoms[offset++];
-            int atom4 = periodicImproperAtoms[offset++];
-            /* periodicImproperForce->addTorsion(atom1, atom2, atom3, atom4, */
-            periodicForce->addTorsion(atom1, atom2, atom3, atom4,
-                                      idef.iparams[type].pdihs.mult,
-                                      idef.iparams[type].pdihs.phiA*M_PI/180.0,
-                                      idef.iparams[type].pdihs.cpA);
-        }
-
-        /* Ryckaert-Bellemans dihedrals */
-        const int* rbAtoms = (int*) idef.il[F_RBDIHS].iatoms;
-        RBTorsionForce* rbForce = new RBTorsionForce();
-        sys->addForce(rbForce);
-        offset = 0;
-        for (int i = 0; i < numRB; ++i)
-        {
-            int type = rbAtoms[offset++];
-            int atom1 = rbAtoms[offset++];
-            int atom2 = rbAtoms[offset++];
-            int atom3 = rbAtoms[offset++];
-            int atom4 = rbAtoms[offset++];
-            rbForce->addTorsion(atom1, atom2, atom3, atom4,
-                                idef.iparams[type].rbdihs.rbcA[0], idef.iparams[type].rbdihs.rbcA[1],
-                                idef.iparams[type].rbdihs.rbcA[2], idef.iparams[type].rbdihs.rbcA[3],
-                                idef.iparams[type].rbdihs.rbcA[4], idef.iparams[type].rbdihs.rbcA[5]);
-        }
-
-               /* Set improper dihedral terms (as in CHARMM FF) */
-        const int* improperDihAtoms = (int*) idef.il[F_IDIHS].iatoms;
-               CustomTorsionForce* improperDihForce = new CustomTorsionForce("2.0*k*asin(sin((theta-theta0)/2))^2");
-        sys->addForce(improperDihForce);
-               improperDihForce->addPerTorsionParameter("k");
-               improperDihForce->addPerTorsionParameter("theta0");
-               vector<double> improperDihParameters(2);
-        offset = 0;
-        for (int i = 0; i < numImproperDih; ++i)
-        {
-            int type = improperDihAtoms[offset++];
-            int atom1 = improperDihAtoms[offset++];
-            int atom2 = improperDihAtoms[offset++];
-            int atom3 = improperDihAtoms[offset++];
-            int atom4 = improperDihAtoms[offset++];
-                       improperDihParameters[0] = idef.iparams[type].harmonic.krA;
-                       improperDihParameters[1] = idef.iparams[type].harmonic.rA*M_PI/180.0;
-            improperDihForce->addTorsion(atom1, atom2, atom3, atom4,
-                                improperDihParameters);
-        }
-
-        /* Set nonbonded parameters and masses. */
-        int ntypes = fr->ntype;
-        int* types = mdatoms->typeA;
-        real* nbfp = fr->nbfp;
-        real* charges = mdatoms->chargeA;
-        real* masses = mdatoms->massT;
-        NonbondedForce* nonbondedForce = new NonbondedForce();
-        sys->addForce(nonbondedForce);
-        
-        switch (ir->ePBC)
-        {
-        case epbcNONE:
-            if (ir->rcoulomb == 0)
-            {
-                nonbondedForce->setNonbondedMethod(NonbondedForce::NoCutoff);
-            }
-            else
-            {
-                nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
-            }
-            break;
-        case epbcXYZ:
-            switch (ir->coulombtype)
-            {
-            case eelCUT:
-            case eelRF:
-            case eelGRF:
-            case eelRF_NEC:
-            case eelRF_ZERO:
-                nonbondedForce->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
-                break;
-
-            case eelEWALD:
-                nonbondedForce->setNonbondedMethod(NonbondedForce::Ewald);
-                break;
-
-            case eelPME:
-                nonbondedForce->setNonbondedMethod(NonbondedForce::PME);
-                break;
-
-            default:
-                gmx_fatal(FARGS,"Internal error: you should not see this message, it means that the"
-                          "electrosatics option check failed. Please report this error!");
-            }        
-            sys->setDefaultPeriodicBoxVectors(Vec3(state->box[0][0], 0, 0),
-                                       Vec3(0, state->box[1][1], 0), Vec3(0, 0, state->box[2][2]));                    
-            nonbondedForce->setCutoffDistance(ir->rcoulomb);
-           
-            break;
-        default:            
-            gmx_fatal(FARGS,"OpenMM supports only full periodic boundary conditions "
-                              "(pbc = xyz), or none (pbc = no).");
-        }
-
-
-        /* Fix for PME and Ewald error tolerance 
-         *
-                *  OpenMM uses approximate formulas to calculate the Ewald parameter:
-                *  alpha = (1.0/cutoff)*sqrt(-log(2.0*tolerlance));
-                *  and the grid spacing for PME:
-                *  gridX = ceil(2*alpha*box[0][0]/3*(pow(tol, 0.2)))
-                *  gridY = ceil(2*alpha*box[1][1]/3*(pow(tol, 0.2)));
-                *  gridZ = ceil(2*alpha*box[2][2]/3*(pow(tol, 0.2)));
-                *
-                *  
-                *  If the default ewald_rtol=1e-5 is used we silently adjust the value to the 
-                *  OpenMM default of 5e-4 otherwise a warning is issued about the action taken. 
-                *
-               */
-        double corr_ewald_rtol = 50.0 * ir->ewald_rtol;
-        if ((ir->ePBC == epbcXYZ) && 
-            (ir->coulombtype == eelEWALD || ir->coulombtype == eelPME))
-        {
-            if (debug)
-            {
-                fprintf(debug, ">> ewald_rtol = %e (corrected = %e) \n",
-                    ir->ewald_rtol, corr_ewald_rtol);
-            }
-
-            if (fabs(ir->ewald_rtol - 1e-5) > 1e-10)
-            {
-                gmx_warning("OpenMM uses the ewald_rtol parameter with approximate formulas "
-                        "to calculate the alpha and grid spacing parameters of the Ewald "
-                        "and PME methods. This tolerance need to be corrected in order to get "
-                        "settings close to the ones used in GROMACS. Although the internal correction "
-                        "should work for any reasonable value of ewald_rtol, using values other than "
-                        "the default 1e-5 might cause incorrect behavior.");
-
-                if (corr_ewald_rtol > 1)
-                {
-                    gmx_fatal(FARGS, "The ewald_rtol accuracy term is >1 after the "
-                            "adjustment for OpenMM (%e)", corr_ewald_rtol);
-                }
-            }
-            nonbondedForce->setEwaldErrorTolerance(corr_ewald_rtol);
-        }
-
-        for (int i = 0; i < numAtoms; ++i)
-        {
-            /* nbfp now includes the 6.0/12.0 derivative prefactors to save flops in kernels*/
-            double c12 = nbfp[types[i]*2*ntypes+types[i]*2+1]/12.0;
-            double c6 = nbfp[types[i]*2*ntypes+types[i]*2]/6.0;
-            double sigma=0.0, epsilon=0.0;
-            convert_c_12_6(c12, c6, &sigma, &epsilon);
-            nonbondedForce->addParticle(charges[i], sigma, epsilon);
-            sys->addParticle(masses[i]);
-        }
-
-        // Build a table of all exclusions.
-        vector<set<int> > exclusions(numAtoms);
-        for (int i = 0; i < numAtoms; i++)
-        {
-            int start = top->excls.index[i];
-            int end = top->excls.index[i+1];
-            for (int j = start; j < end; j++)
-                exclusions[i].insert(top->excls.a[j]);
-        }
-
-        // Record the 1-4 interactions, and remove them from the list of exclusions.
-        const int* nb14Atoms = (int*) idef.il[F_LJ14].iatoms;
-        offset = 0;
-        for (int i = 0; i < num14; ++i)
-        {
-            int type = nb14Atoms[offset++];
-            int atom1 = nb14Atoms[offset++];
-            int atom2 = nb14Atoms[offset++];
-            double sigma=0, epsilon=0;
-            convert_c_12_6(idef.iparams[type].lj14.c12A, 
-                    idef.iparams[type].lj14.c6A,
-                    &sigma, &epsilon);
-            nonbondedForce->addException(atom1, atom2,
-                                         fr->fudgeQQ*charges[atom1]*charges[atom2], sigma, epsilon);
-            exclusions[atom1].erase(atom2);
-            exclusions[atom2].erase(atom1);
-        }
-
-        // Record exclusions.
-        for (int i = 0; i < numAtoms; i++)
-        {
-            for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
-            {
-                if (i < *iter)
-                {
-                    nonbondedForce->addException(i, *iter, 0.0, 1.0, 0.0);
-                }
-            }
-        }
-
-        // Add GBSA if needed.
-        if (ir->implicit_solvent == eisGBSA)
-        {
-            gmx_warning("The OBC scale factors alpha, beta and gamma are hardcoded in OpenMM with the default Gromacs values.");
-            t_atoms atoms       = gmx_mtop_global_atoms(top_global);
-            GBSAOBCForce* gbsa  = new GBSAOBCForce();
-
-            sys->addForce(gbsa);
-            gbsa->setSoluteDielectric(ir->epsilon_r);
-            gbsa->setSolventDielectric(ir->gb_epsilon_solvent);
-            gbsa->setCutoffDistance(nonbondedForce->getCutoffDistance());
-            if (nonbondedForce->getNonbondedMethod() == NonbondedForce::NoCutoff)
-                gbsa->setNonbondedMethod(GBSAOBCForce::NoCutoff);
-            else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
-                gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
-            else if (nonbondedForce->getNonbondedMethod() == NonbondedForce::CutoffPeriodic)
-                gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
-            else
-                gmx_fatal(FARGS,"OpenMM supports only Reaction-Field electrostatics with OBC/GBSA.");
-
-            for (int i = 0; i < numAtoms; ++i)
-            {
-                gbsa->addParticle(charges[i],
-                                  top_global->atomtypes.gb_radius[atoms.atom[i].type],
-                                  top_global->atomtypes.S_hct[atoms.atom[i].type]);
-            }
-            free_t_atoms(&atoms, FALSE);
-        }
-
-        // Set constraints.
-        const int* constraintAtoms = (int*) idef.il[F_CONSTR].iatoms;
-        offset = 0;
-        for (int i = 0; i < numConstraints; ++i)
-        {
-            int type = constraintAtoms[offset++];
-            int atom1 = constraintAtoms[offset++];
-            int atom2 = constraintAtoms[offset++];
-            sys->addConstraint(atom1, atom2, idef.iparams[type].constr.dA);
-        }
-        const int* settleAtoms = (int*) idef.il[F_SETTLE].iatoms;
-        offset = 0;
-        for (int i = 0; i < numSettle; ++i)
-        {
-            int type = settleAtoms[offset++];
-            int oxygen = settleAtoms[offset++];
-            sys->addConstraint(oxygen, oxygen+1, idef.iparams[type].settle.doh);
-            sys->addConstraint(oxygen, oxygen+2, idef.iparams[type].settle.doh);
-            sys->addConstraint(oxygen+1, oxygen+2, idef.iparams[type].settle.dhh);
-        }
-
-        // Create an integrator for simulating the system.
-        double friction = (ir->opts.tau_t[0] == 0.0 ? 0.0 : 1.0/ir->opts.tau_t[0]);
-        Integrator* integ;
-        if (ir->eI == eiBD)
-        {
-            integ = new BrownianIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
-            static_cast<BrownianIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed); 
-        }
-        else if (EI_SD(ir->eI))
-        {
-            integ = new LangevinIntegrator(ir->opts.ref_t[0], friction, ir->delta_t);
-            static_cast<LangevinIntegrator*>(integ)->setRandomNumberSeed(ir->ld_seed); 
-        }
-        else 
-        {
-            integ = new VerletIntegrator(ir->delta_t);
-            if ( ir->etc != etcNO)
-            {
-                AndersenThermostat* thermostat = new AndersenThermostat(ir->opts.ref_t[0], friction); 
-                sys->addForce(thermostat);
-            }           
-        }
-
-               // Add pressure coupling
-        if (ir->epc != epcNO)
-               {
-          // convert gromacs pressure tensor to a scalar
-          double pressure = (ir->ref_p[0][0] + ir->ref_p[1][1] + ir->ref_p[2][2]) / 3.0;
-          int frequency = int(ir->tau_p / ir->delta_t); // update frequency in time steps
-          if (frequency < 1) frequency = 1;
-          double temperature = ir->opts.ref_t[0]; // in kelvin
-          sys->addForce(new MonteCarloBarostat(pressure, temperature, frequency));
-               }
-
-        integ->setConstraintTolerance(ir->shake_tol);
-
-        // Create a context and initialize it.
-        Context* context = NULL;
-
-        /*      
-        OpenMM could automatically select the "best" GPU, however we're not't 
-        going to let it do that for now, as the current algorithm is very rudimentary
-        and we anyway support only CUDA.        
-        if (platformOptStr == NULL || platformOptStr == "")
-        {
-            context = new Context(*sys, *integ);
-        }
-        else
-        */        
-        {
-            /* which platform should we use */
-            for (int i = 0; i < (int)Platform::getNumPlatforms() && context == NULL; i++)
-            {
-                if (isStringEqNCase(opt->getOptionValue("platform"), Platform::getPlatform(i).getName()))
-                {
-                    Platform& platform = Platform::getPlatform(i);
-                    // set standard properties
-                    platform.setPropertyDefaultValue("CudaDevice", opt->getOptionValue("deviceid"));
-                    // TODO add extra properties
-                    context = new Context(*sys, *integ, platform);
-                }
-            }
-            if (context == NULL)
-            {
-                gmx_fatal(FARGS, "The requested platform \"%s\" could not be found.", 
-                        opt->getOptionValue("platform").c_str());
-            }
-        }
-
-        Platform& platform = context->getPlatform();
-        fprintf(fplog, "Gromacs will use the OpenMM platform: %s\n", platform.getName().c_str());
-
-        const vector<string>& properties = platform.getPropertyNames();
-        if (debug)
-        {
-            for (int i = 0; i < (int)properties.size(); i++)
-            {
-                fprintf(debug, ">> %s: %s\n", properties[i].c_str(), 
-                        platform.getPropertyValue(*context, properties[i]).c_str());
-            }
-        }
-
-        /* only for CUDA */
-        if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
-        {
-            int tmp;
-            if (!from_string<int>(tmp, platform.getPropertyValue(*context, "CudaDevice"), std::dec))
-            {
-                gmx_fatal(FARGS, "Internal error: couldn't determine the device selected by OpenMM");
-
-            }
-
-            /* For now this is just to double-check if OpenMM selected the GPU we wanted,
-            but when we'll let OpenMM select the GPU automatically, it will query the deviceId.
-            */            
-            if (tmp != devId)
-            {
-                gmx_fatal(FARGS, "Internal error: OpenMM is using device #%d"
-                        "while initialized for device #%d", tmp, devId);
-            }        
-            
-            /* check GPU compatibility */
-            char gpuname[STRLEN];
-            devId = atoi(opt->getOptionValue("deviceid").c_str());
-            if (!is_gmx_openmm_supported_gpu(-1, gpuname))
-            {
-                if (!gmx_strcasecmp(opt->getOptionValue("force-device").c_str(), "yes"))
-                {
-                    sprintf(warn_buf, "Non-supported GPU selected (#%d, %s), forced continuing."
-                            "Note, that the simulation can be slow or it migth even crash.", 
-                            devId, gpuname);
-                    fprintf(fplog, "%s\n", warn_buf);
-                    gmx_warning(warn_buf);
-                }
-                else
-                {
-                    gmx_fatal(FARGS, "The selected GPU (#%d, %s) is not supported by Gromacs! "
-                              "Most probably you have a low-end GPU which would not perform well, " 
-                              "or new hardware that has not been tested with the current release. "
-                              "If you still want to try using the device, use the force-device=yes option.", 
-                              devId, gpuname);
-                }
-            }
-            else
-            {
-                fprintf(fplog, "Gromacs will run on the GPU #%d (%s).\n", devId, gpuname);
-            }
-        }
-        
-        /* only for CUDA */
-        if (isStringEqNCase(opt->getOptionValue("platform"), "CUDA"))
-        {
-            /* pre-simulation memtest */
-            runMemtest(fplog, -1, "Pre", opt);
-        }
-
-        vector<Vec3> pos(numAtoms);
-        vector<Vec3> vel(numAtoms);
-        for (int i = 0; i < numAtoms; ++i)
-        {
-            pos[i] = Vec3(state->x[i][0], state->x[i][1], state->x[i][2]);
-            vel[i] = Vec3(state->v[i][0], state->v[i][1], state->v[i][2]);
-        }
-        context->setPositions(pos);
-        context->setVelocities(vel);
-
-        // Return a structure containing the system, integrator, and context.
-        OpenMMData* data = new OpenMMData();
-        data->system = sys;
-        data->integrator = integ;
-        data->context = context;
-        data->removeCM = (ir->nstcomm > 0);
-        data->platformOpt = opt;
-        return data;
-    }
-    catch (std::exception& e)
-    {
-        gmx_fatal(FARGS, "OpenMM exception caught while initializating: %s", e.what());
-    } 
-    return NULL; /* just to avoid warnings */
-}
-
-/*!
- * \brief Integrate one step.
- *
- * \param[in] data  OpenMMData object created by openmm_init().
- */
-void openmm_take_one_step(void* data)
-{
-    // static int step = 0; printf("----> taking step #%d\n", step++);
-    try
-    {
-        static_cast<OpenMMData*>(data)->integrator->step(1);
-    }
-    catch (std::exception& e)
-    {
-        gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
-    }
-}
-
-/*!
- * \brief Integrate n steps.
- *
- * \param[in] data  OpenMMData object created by openmm_init().
- */
-void openmm_take_steps(void* data, int nstep)
-{
-    try
-    {
-        static_cast<OpenMMData*>(data)->integrator->step(nstep);
-    }
-    catch (std::exception& e)
-    {
-        gmx_fatal(FARGS, "OpenMM exception caught while taking a step: %s", e.what());
-    }
-}
-
-/*!
- * \brief Clean up the data structures cretead for OpenMM.
- *
- * \param[in] log   Log file pointer.
- * \param[in] data  OpenMMData object created by openmm_init().
- */
-void openmm_cleanup(FILE* fplog, void* data)
-{
-    OpenMMData* d = static_cast<OpenMMData*>(data);
-    /* only for CUDA */
-    if (isStringEqNCase(d->platformOpt->getOptionValue("platform"), "CUDA"))
-    {
-        /* post-simulation memtest */
-        runMemtest(fplog, -1, "Post", d->platformOpt);
-    }
-    delete d->system;
-    delete d->integrator;
-    delete d->context;
-    delete d->platformOpt;
-    delete d;
-}
-
-/*!
- * \brief Copy the current state information from OpenMM into the Gromacs data structures.
- * 
- * This function results in the requested proprties to be copied from the 
- * GPU to host. As this represents a bottleneck, the frequency of pulling data
- * should be minimized. 
- *
- * \param[in]   data        OpenMMData object created by openmm_init().
- * \param[out]  time        Simulation time for which the state was created.
- * \param[out]  state       State of the system: coordinates and velocities.
- * \param[out]  f           Forces.
- * \param[out]  enerd       Energies.
- * \param[in]   includePos  True if coordinates are requested.
- * \param[in]   includeVel  True if velocities are requested. 
- * \param[in]   includeForce True if forces are requested. 
- * \param[in]   includeEnergy True if energies are requested. 
- */
-void openmm_copy_state(void *data,
-                       t_state *state, double *time,
-                       rvec f[], gmx_enerdata_t *enerd,
-                       gmx_bool includePos, gmx_bool includeVel, gmx_bool includeForce, gmx_bool includeEnergy)
-{
-    int types = 0;
-    if (includePos)
-        types += State::Positions;
-    if (includeVel)
-        types += State::Velocities;
-    if (includeForce)
-        types += State::Forces;
-    if (includeEnergy)
-        types += State::Energy;
-    if (types == 0)
-        return;
-    try
-    {
-        State currentState = static_cast<OpenMMData*>(data)->context->getState(types);
-        int numAtoms =  static_cast<OpenMMData*>(data)->system->getNumParticles();
-        if (includePos)
-        {
-            for (int i = 0; i < numAtoms; i++)
-            {
-                Vec3 x = currentState.getPositions()[i];
-                state->x[i][0] = x[0];
-                state->x[i][1] = x[1];
-                state->x[i][2] = x[2];
-            }
-        }
-        if (includeVel)
-        {
-            for (int i = 0; i < numAtoms; i++)
-            {
-                Vec3 v = currentState.getVelocities()[i];
-                state->v[i][0] = v[0];
-                state->v[i][1] = v[1];
-                state->v[i][2] = v[2];
-            }
-        }
-        if (includeForce)
-        {
-            for (int i = 0; i < numAtoms; i++)
-            {
-                Vec3 force = currentState.getForces()[i];
-                f[i][0] = force[0];
-                f[i][1] = force[1];
-                f[i][2] = force[2];
-            }
-        }
-        if (includeEnergy)
-        {
-            int numConstraints = static_cast<OpenMMData*>(data)->system->getNumConstraints();
-            int dof = 3*numAtoms-numConstraints;
-            if (static_cast<OpenMMData*>(data)->removeCM)
-                dof -= 3;
-            enerd->term[F_EPOT] = currentState.getPotentialEnergy();
-            enerd->term[F_EKIN] = currentState.getKineticEnergy();
-            enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
-            enerd->term[F_TEMP] = 2.0*enerd->term[F_EKIN]/dof/BOLTZ;
-        }
-        *time = currentState.getTime();
-    }
-    catch (std::exception& e)
-    {
-        gmx_fatal(FARGS, "OpenMM exception caught while retrieving state information: %s", e.what());
-    }
-}
diff --git a/src/contrib/openmm_wrapper.h b/src/contrib/openmm_wrapper.h
deleted file mode 100644 (file)
index 7013426..0000000
+++ /dev/null
@@ -1,68 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2010, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
- * Copyright (c) 2012, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-
-#ifndef _OPENMM_WRAPPER_H_
-#define _OPENMM_WRAPPER_H_
-
-#ifdef __cplusplus
-extern "C"
-{
-#endif
-
-void* openmm_init(FILE *fplog, const char *platformOptStr,
-                    t_inputrec *ir,
-                    gmx_mtop_t *top_global, gmx_localtop_t *top,
-                    t_mdatoms *mdatoms, t_forcerec *fr, t_state *state);
-
-void openmm_take_one_step(void* data);
-
-void openmm_take_steps(void* data, int nsteps);
-
-void openmm_copy_state(void *data,
-                        t_state *state, double *time,
-                        rvec f[], gmx_enerdata_t *enerd,
-                        gmx_bool includePos, gmx_bool includeVel, gmx_bool includeForce, gmx_bool includeEnergy);
-
-void openmm_cleanup(FILE *fplog, void* data);
-
-#ifdef __cplusplus
-} // extern "C"
-#endif
-
-#endif /* _OPENMM_WRAPPER_H_ */
-
diff --git a/src/contrib/runner_openmm.c b/src/contrib/runner_openmm.c
deleted file mode 100644 (file)
index 1d15fd0..0000000
+++ /dev/null
@@ -1,1993 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
- * Copyright (c) 2012, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-#if defined(HAVE_SCHED_H) && defined(HAVE_SCHED_GETAFFINITY)
-#define _GNU_SOURCE
-#include <sched.h>
-#include <sys/syscall.h>
-#endif
-#include <signal.h>
-#include <stdlib.h>
-#include <string.h>
-#include <assert.h>
-
-#include "typedefs.h"
-#include "smalloc.h"
-#include "sysstuff.h"
-#include "statutil.h"
-#include "mdrun.h"
-#include "md_logging.h"
-#include "md_support.h"
-#include "network.h"
-#include "pull.h"
-#include "names.h"
-#include "disre.h"
-#include "orires.h"
-#include "pme.h"
-#include "mdatoms.h"
-#include "repl_ex.h"
-#include "qmmm.h"
-#include "mpelogging.h"
-#include "domdec.h"
-#include "partdec.h"
-#include "coulomb.h"
-#include "constr.h"
-#include "mvdata.h"
-#include "checkpoint.h"
-#include "mtop_util.h"
-#include "sighandler.h"
-#include "tpxio.h"
-#include "txtdump.h"
-#include "gmx_detect_hardware.h"
-#include "gmx_omp_nthreads.h"
-#include "pull_rotation.h"
-#include "calc_verletbuf.h"
-#include "../mdlib/nbnxn_search.h"
-#include "../mdlib/nbnxn_consts.h"
-#include "gmx_fatal_collective.h"
-#include "membed.h"
-#include "md_openmm.h"
-#include "gmx_omp.h"
-
-#include "thread_mpi/threads.h"
-
-#include "gromacs/utility/gmxmpi.h"
-
-#ifdef GMX_FAHCORE
-#include "corewrap.h"
-#endif
-
-#include "gpu_utils.h"
-#include "nbnxn_cuda_data_mgmt.h"
-
-typedef struct { 
-    gmx_integrator_t *func;
-} gmx_intp_t;
-
-/* The array should match the eI array in include/types/enums.h */
-const gmx_intp_t integrator[eiNR] = { {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm},{do_md_openmm}};
-
-gmx_large_int_t     deform_init_init_step_tpx;
-matrix              deform_init_box_tpx;
-#ifdef GMX_THREAD_MPI
-tMPI_Thread_mutex_t deform_init_box_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
-#endif
-
-
-#ifdef GMX_THREAD_MPI
-struct mdrunner_arglist
-{
-    gmx_hw_opt_t *hw_opt;
-    FILE *fplog;
-    t_commrec *cr;
-    int nfile;
-    const t_filenm *fnm;
-    output_env_t oenv;
-    gmx_bool bVerbose;
-    gmx_bool bCompact;
-    int nstglobalcomm;
-    ivec ddxyz;
-    int dd_node_order;
-    real rdd;
-    real rconstr;
-    const char *dddlb_opt;
-    real dlb_scale;
-    const char *ddcsx;
-    const char *ddcsy;
-    const char *ddcsz;
-    const char *nbpu_opt;
-    int nsteps_cmdline;
-    int nstepout;
-    int resetstep;
-    int nmultisim;
-    int repl_ex_nst;
-    int repl_ex_nex;
-    int repl_ex_seed;
-    real pforce;
-    real cpt_period;
-    real max_hours;
-    const char *deviceOptions;
-    unsigned long Flags;
-    int ret; /* return value */
-};
-
-
-/* The function used for spawning threads. Extracts the mdrunner() 
-   arguments from its one argument and calls mdrunner(), after making
-   a commrec. */
-static void mdrunner_start_fn(void *arg)
-{
-    struct mdrunner_arglist *mda=(struct mdrunner_arglist*)arg;
-    struct mdrunner_arglist mc=*mda; /* copy the arg list to make sure 
-                                        that it's thread-local. This doesn't
-                                        copy pointed-to items, of course,
-                                        but those are all const. */
-    t_commrec *cr;  /* we need a local version of this */
-    FILE *fplog=NULL;
-    t_filenm *fnm;
-
-    fnm = dup_tfn(mc.nfile, mc.fnm);
-
-    cr = init_par_threads(mc.cr);
-
-    if (MASTER(cr))
-    {
-        fplog=mc.fplog;
-    }
-
-    mda->ret=mdrunner(mc.hw_opt, fplog, cr, mc.nfile, fnm, mc.oenv, 
-                      mc.bVerbose, mc.bCompact, mc.nstglobalcomm, 
-                      mc.ddxyz, mc.dd_node_order, mc.rdd,
-                      mc.rconstr, mc.dddlb_opt, mc.dlb_scale, 
-                      mc.ddcsx, mc.ddcsy, mc.ddcsz,
-                      mc.nbpu_opt,
-                      mc.nsteps_cmdline, mc.nstepout, mc.resetstep,
-                      mc.nmultisim, mc.repl_ex_nst, mc.repl_ex_nex, mc.repl_ex_seed, mc.pforce, 
-                      mc.cpt_period, mc.max_hours, mc.deviceOptions, mc.Flags);
-}
-
-/* called by mdrunner() to start a specific number of threads (including 
-   the main thread) for thread-parallel runs. This in turn calls mdrunner()
-   for each thread. 
-   All options besides nthreads are the same as for mdrunner(). */
-static t_commrec *mdrunner_start_threads(gmx_hw_opt_t *hw_opt, 
-              FILE *fplog,t_commrec *cr,int nfile, 
-              const t_filenm fnm[], const output_env_t oenv, gmx_bool bVerbose,
-              gmx_bool bCompact, int nstglobalcomm,
-              ivec ddxyz,int dd_node_order,real rdd,real rconstr,
-              const char *dddlb_opt,real dlb_scale,
-              const char *ddcsx,const char *ddcsy,const char *ddcsz,
-              const char *nbpu_opt,
-              int nsteps_cmdline, int nstepout,int resetstep,
-              int nmultisim,int repl_ex_nst,int repl_ex_nex, int repl_ex_seed,
-              real pforce,real cpt_period, real max_hours, 
-              const char *deviceOptions, unsigned long Flags)
-{
-    int ret;
-    struct mdrunner_arglist *mda;
-    t_commrec *crn; /* the new commrec */
-    t_filenm *fnmn;
-
-    /* first check whether we even need to start tMPI */
-    if (hw_opt->nthreads_tmpi < 2)
-    {
-        return cr;
-    }
-
-    /* a few small, one-time, almost unavoidable memory leaks: */
-    snew(mda,1);
-    fnmn=dup_tfn(nfile, fnm);
-
-    /* fill the data structure to pass as void pointer to thread start fn */
-    mda->hw_opt=hw_opt;
-    mda->fplog=fplog;
-    mda->cr=cr;
-    mda->nfile=nfile;
-    mda->fnm=fnmn;
-    mda->oenv=oenv;
-    mda->bVerbose=bVerbose;
-    mda->bCompact=bCompact;
-    mda->nstglobalcomm=nstglobalcomm;
-    mda->ddxyz[XX]=ddxyz[XX];
-    mda->ddxyz[YY]=ddxyz[YY];
-    mda->ddxyz[ZZ]=ddxyz[ZZ];
-    mda->dd_node_order=dd_node_order;
-    mda->rdd=rdd;
-    mda->rconstr=rconstr;
-    mda->dddlb_opt=dddlb_opt;
-    mda->dlb_scale=dlb_scale;
-    mda->ddcsx=ddcsx;
-    mda->ddcsy=ddcsy;
-    mda->ddcsz=ddcsz;
-    mda->nbpu_opt=nbpu_opt;
-    mda->nsteps_cmdline=nsteps_cmdline;
-    mda->nstepout=nstepout;
-    mda->resetstep=resetstep;
-    mda->nmultisim=nmultisim;
-    mda->repl_ex_nst=repl_ex_nst;
-    mda->repl_ex_nex=repl_ex_nex;
-    mda->repl_ex_seed=repl_ex_seed;
-    mda->pforce=pforce;
-    mda->cpt_period=cpt_period;
-    mda->max_hours=max_hours;
-    mda->deviceOptions=deviceOptions;
-    mda->Flags=Flags;
-
-    /* now spawn new threads that start mdrunner_start_fn(), while 
-       the main thread returns */
-    ret=tMPI_Init_fn(TRUE, hw_opt->nthreads_tmpi,
-                     (hw_opt->bThreadPinning ? TMPI_AFFINITY_ALL_CORES : TMPI_AFFINITY_NONE),
-                     mdrunner_start_fn, (void*)(mda) );
-    if (ret!=TMPI_SUCCESS)
-        return NULL;
-
-    /* make a new comm_rec to reflect the new situation */
-    crn=init_par_threads(cr);
-    return crn;
-}
-
-
-static int get_tmpi_omp_thread_division(const gmx_hw_info_t *hwinfo,
-                                        const gmx_hw_opt_t *hw_opt,
-                                        int nthreads_tot,
-                                        int ngpu)
-{
-    int nthreads_tmpi;
-
-    /* There are no separate PME nodes here, as we ensured in
-     * check_and_update_hw_opt that nthreads_tmpi>0 with PME nodes
-     * and a conditional ensures we would not have ended up here.
-     * Note that separate PME nodes might be switched on later.
-     */
-    if (ngpu > 0)
-    {
-        nthreads_tmpi = ngpu;
-        if (nthreads_tot > 0 && nthreads_tot < nthreads_tmpi)
-        {
-            nthreads_tmpi = nthreads_tot;
-        }
-    }
-    else if (hw_opt->nthreads_omp > 0)
-    {
-        /* Here we could oversubscribe, when we do, we issue a warning later */
-        nthreads_tmpi = max(1,nthreads_tot/hw_opt->nthreads_omp);
-    }
-    else
-    {
-        /* TODO choose nthreads_omp based on hardware topology
-           when we have a hardware topology detection library */
-        /* In general, when running up to 4 threads, OpenMP should be faster.
-         * Note: on AMD Bulldozer we should avoid running OpenMP over two dies.
-         * On Intel>=Nehalem running OpenMP on a single CPU is always faster,
-         * even on two CPUs it's usually faster (but with many OpenMP threads
-         * it could be faster not to use HT, currently we always use HT).
-         * On Nehalem/Westmere we want to avoid running 16 threads over
-         * two CPUs with HT, so we need a limit<16; thus we use 12.
-         * A reasonable limit for Intel Sandy and Ivy bridge,
-         * not knowing the topology, is 16 threads.
-         */
-        const int nthreads_omp_always_faster             =  4;
-        const int nthreads_omp_always_faster_Nehalem     = 12;
-        const int nthreads_omp_always_faster_SandyBridge = 16;
-        const int first_model_Nehalem     = 0x1A;
-        const int first_model_SandyBridge = 0x2A;
-        gmx_bool bIntel_Family6;
-
-        bIntel_Family6 =
-            (gmx_cpuid_vendor(hwinfo->cpuid_info) == GMX_CPUID_VENDOR_INTEL &&
-             gmx_cpuid_family(hwinfo->cpuid_info) == 6);
-
-        if (nthreads_tot <= nthreads_omp_always_faster ||
-            (bIntel_Family6 &&
-             ((gmx_cpuid_model(hwinfo->cpuid_info) >= nthreads_omp_always_faster_Nehalem && nthreads_tot <= nthreads_omp_always_faster_Nehalem) ||
-              (gmx_cpuid_model(hwinfo->cpuid_info) >= nthreads_omp_always_faster_SandyBridge && nthreads_tot <= nthreads_omp_always_faster_SandyBridge))))
-        {
-            /* Use pure OpenMP parallelization */
-            nthreads_tmpi = 1;
-        }
-        else
-        {
-            /* Don't use OpenMP parallelization */
-            nthreads_tmpi = nthreads_tot;
-        }
-    }
-
-    return nthreads_tmpi;
-}
-
-
-/* Get the number of threads to use for thread-MPI based on how many
- * were requested, which algorithms we're using,
- * and how many particles there are.
- * At the point we have already called check_and_update_hw_opt.
- * Thus all options should be internally consistent and consistent
- * with the hardware, except that ntmpi could be larger than #GPU.
- */
-static int get_nthreads_mpi(gmx_hw_info_t *hwinfo,
-                            gmx_hw_opt_t *hw_opt,
-                            t_inputrec *inputrec, gmx_mtop_t *mtop,
-                            const t_commrec *cr,
-                            FILE *fplog)
-{
-    int nthreads_hw,nthreads_tot_max,nthreads_tmpi,nthreads_new,ngpu;
-    int min_atoms_per_mpi_thread;
-    char *env;
-    char sbuf[STRLEN];
-    gmx_bool bCanUseGPU;
-
-    if (hw_opt->nthreads_tmpi > 0)
-    {
-        /* Trivial, return right away */
-        return hw_opt->nthreads_tmpi;
-    }
-
-    nthreads_hw = hwinfo->nthreads_hw_avail;
-
-    /* How many total (#tMPI*#OpenMP) threads can we start? */ 
-    if (hw_opt->nthreads_tot > 0)
-    {
-        nthreads_tot_max = hw_opt->nthreads_tot;
-    }
-    else
-    {
-        nthreads_tot_max = nthreads_hw;
-    }
-
-    bCanUseGPU = (inputrec->cutoff_scheme == ecutsVERLET && hwinfo->bCanUseGPU);
-    if (bCanUseGPU)
-    {
-        ngpu = hwinfo->gpu_info.ncuda_dev_use;
-    }
-    else
-    {
-        ngpu = 0;
-    }
-
-    nthreads_tmpi =
-        get_tmpi_omp_thread_division(hwinfo,hw_opt,nthreads_tot_max,ngpu);
-
-    if (inputrec->eI == eiNM || EI_TPI(inputrec->eI))
-    {
-        /* Steps are divided over the nodes iso splitting the atoms */
-        min_atoms_per_mpi_thread = 0;
-    }
-    else
-    {
-        if (bCanUseGPU)
-        {
-            min_atoms_per_mpi_thread = MIN_ATOMS_PER_GPU;
-        }
-        else
-        {
-            min_atoms_per_mpi_thread = MIN_ATOMS_PER_MPI_THREAD;
-        }
-    }
-
-    /* Check if an algorithm does not support parallel simulation.  */
-    if (nthreads_tmpi != 1 &&
-        ( inputrec->eI == eiLBFGS ||
-          inputrec->coulombtype == eelEWALD ) )
-    {
-        nthreads_tmpi = 1;
-
-        md_print_warn(cr,fplog,"The integration or electrostatics algorithm doesn't support parallel runs. Using a single thread-MPI thread.\n");
-        if (hw_opt->nthreads_tmpi > nthreads_tmpi)
-        {
-            gmx_fatal(FARGS,"You asked for more than 1 thread-MPI thread, but an algorithm doesn't support that");
-        }
-    }
-    else if (mtop->natoms/nthreads_tmpi < min_atoms_per_mpi_thread)
-    {
-        /* the thread number was chosen automatically, but there are too many
-           threads (too few atoms per thread) */
-        nthreads_new = max(1,mtop->natoms/min_atoms_per_mpi_thread);
-
-        /* Avoid partial use of Hyper-Threading */
-        if (gmx_cpuid_x86_smt(hwinfo->cpuid_info) == GMX_CPUID_X86_SMT_ENABLED &&
-            nthreads_new > nthreads_hw/2 && nthreads_new < nthreads_hw)
-        {
-            nthreads_new = nthreads_hw/2;
-        }
-
-        /* Avoid large prime numbers in the thread count */
-        if (nthreads_new >= 6)
-        {
-            /* Use only 6,8,10 with additional factors of 2 */
-            int fac;
-
-            fac = 2;
-            while (3*fac*2 <= nthreads_new)
-            {
-                fac *= 2;
-            }
-
-            nthreads_new = (nthreads_new/fac)*fac;
-        }
-        else
-        {
-            /* Avoid 5 */
-            if (nthreads_new == 5)
-            {
-                nthreads_new = 4;
-            }
-        }
-
-        nthreads_tmpi = nthreads_new;
-
-        fprintf(stderr,"\n");
-        fprintf(stderr,"NOTE: Parallelization is limited by the small number of atoms,\n");
-        fprintf(stderr,"      only starting %d thread-MPI threads.\n",nthreads_tmpi);
-        fprintf(stderr,"      You can use the -nt and/or -ntmpi option to optimize the number of threads.\n\n");
-    }
-
-    return nthreads_tmpi;
-}
-#endif /* GMX_THREAD_MPI */
-
-
-/* Environment variable for setting nstlist */
-static const char*  NSTLIST_ENVVAR          =  "GMX_NSTLIST";
-/* Try to increase nstlist when using a GPU with nstlist less than this */
-static const int    NSTLIST_GPU_ENOUGH      = 20;
-/* Increase nstlist until the non-bonded cost increases more than this factor */
-static const float  NBNXN_GPU_LIST_OK_FAC   = 1.25;
-/* Don't increase nstlist beyond a non-bonded cost increases of this factor */
-static const float  NBNXN_GPU_LIST_MAX_FAC  = 1.40;
-
-/* Try to increase nstlist when running on a GPU */
-static void increase_nstlist(FILE *fp,t_commrec *cr,
-                             t_inputrec *ir,const gmx_mtop_t *mtop,matrix box)
-{
-    char *env;
-    int  nstlist_orig,nstlist_prev;
-    verletbuf_list_setup_t ls;
-    real rlist_inc,rlist_ok,rlist_max,rlist_new,rlist_prev;
-    int  i;
-    t_state state_tmp;
-    gmx_bool bBox,bDD,bCont;
-    const char *nstl_fmt="\nFor optimal performance with a GPU nstlist (now %d) should be larger.\nThe optimum depends on your CPU and GPU resources.\nYou might want to try several nstlist values.\n";
-    const char *vbd_err="Can not increase nstlist for GPU run because verlet-buffer-drift is not set or used";
-    const char *box_err="Can not increase nstlist for GPU run because the box is too small";
-    const char *dd_err ="Can not increase nstlist for GPU run because of domain decomposition limitations";
-    char buf[STRLEN];
-
-    /* Number of + nstlist alternative values to try when switching  */
-    const int nstl[]={ 20, 25, 40, 50 };
-#define NNSTL  sizeof(nstl)/sizeof(nstl[0])
-
-    env = getenv(NSTLIST_ENVVAR);
-    if (env == NULL)
-    {
-        if (fp != NULL)
-        {
-            fprintf(fp,nstl_fmt,ir->nstlist);
-        }
-    }
-
-    if (ir->verletbuf_drift == 0)
-    {
-        gmx_fatal(FARGS,"You are using an old tpr file with a GPU, please generate a new tpr file with an up to date version of grompp");
-    }
-
-    if (ir->verletbuf_drift < 0)
-    {
-        if (MASTER(cr))
-        {
-            fprintf(stderr,"%s\n",vbd_err);
-        }
-        if (fp != NULL)
-        {
-            fprintf(fp,"%s\n",vbd_err);
-        }
-
-        return;
-    }
-
-    nstlist_orig = ir->nstlist;
-    if (env != NULL)
-    {
-        sprintf(buf,"Getting nstlist from environment variable GMX_NSTLIST=%s",env);
-        if (MASTER(cr))
-        {
-            fprintf(stderr,"%s\n",buf);
-        }
-        if (fp != NULL)
-        {
-            fprintf(fp,"%s\n",buf);
-        }
-        sscanf(env,"%d",&ir->nstlist);
-    }
-
-    verletbuf_get_list_setup(TRUE,&ls);
-
-    /* Allow rlist to make the list double the size of the cut-off sphere */
-    rlist_inc = nbnxn_get_rlist_effective_inc(NBNXN_GPU_CLUSTER_SIZE,mtop->natoms/det(box));
-    rlist_ok  = (max(ir->rvdw,ir->rcoulomb) + rlist_inc)*pow(NBNXN_GPU_LIST_OK_FAC,1.0/3.0) - rlist_inc;
-    rlist_max = (max(ir->rvdw,ir->rcoulomb) + rlist_inc)*pow(NBNXN_GPU_LIST_MAX_FAC,1.0/3.0) - rlist_inc;
-    if (debug)
-    {
-        fprintf(debug,"GPU nstlist tuning: rlist_inc %.3f rlist_max %.3f\n",
-                rlist_inc,rlist_max);
-    }
-
-    i = 0;
-    nstlist_prev = nstlist_orig;
-    rlist_prev   = ir->rlist;
-    do
-    {
-        if (env == NULL)
-        {
-            ir->nstlist = nstl[i];
-        }
-
-        /* Set the pair-list buffer size in ir */
-        calc_verlet_buffer_size(mtop,det(box),ir,ir->verletbuf_drift,&ls,
-                                NULL,&rlist_new);
-
-        /* Does rlist fit in the box? */
-        bBox = (sqr(rlist_new) < max_cutoff2(ir->ePBC,box));
-        bDD  = TRUE;
-        if (bBox && DOMAINDECOMP(cr))
-        {
-            /* Check if rlist fits in the domain decomposition */
-            if (inputrec2nboundeddim(ir) < DIM)
-            {
-                gmx_incons("Changing nstlist with domain decomposition and unbounded dimensions is not implemented yet");
-            }
-            copy_mat(box,state_tmp.box);
-            bDD = change_dd_cutoff(cr,&state_tmp,ir,rlist_new);
-        }
-
-        bCont = FALSE;
-
-        if (env == NULL)
-        {
-            if (bBox && bDD && rlist_new <= rlist_max)
-            {
-                /* Increase nstlist */
-                nstlist_prev = ir->nstlist;
-                rlist_prev   = rlist_new;
-                bCont = (i+1 < NNSTL && rlist_new < rlist_ok);
-            }
-            else
-            {
-                /* Stick with the previous nstlist */
-                ir->nstlist = nstlist_prev;
-                rlist_new   = rlist_prev;
-                bBox = TRUE;
-                bDD  = TRUE;
-            }
-        }
-
-        i++;
-    }
-    while (bCont);
-
-    if (!bBox || !bDD)
-    {
-        gmx_warning(!bBox ? box_err : dd_err);
-        if (fp != NULL)
-        {
-            fprintf(fp,"\n%s\n",bBox ? box_err : dd_err);
-        }
-        ir->nstlist = nstlist_orig;
-    }
-    else if (ir->nstlist != nstlist_orig || rlist_new != ir->rlist)
-    {
-        sprintf(buf,"Changing nstlist from %d to %d, rlist from %g to %g",
-                nstlist_orig,ir->nstlist,
-                ir->rlist,rlist_new);
-        if (MASTER(cr))
-        {
-            fprintf(stderr,"%s\n\n",buf);
-        }
-        if (fp != NULL)
-        {
-            fprintf(fp,"%s\n\n",buf);
-        }
-        ir->rlist     = rlist_new;
-        ir->rlistlong = rlist_new;
-    }
-}
-
-static void prepare_verlet_scheme(FILE *fplog,
-                                  gmx_hw_info_t *hwinfo,
-                                  t_commrec *cr,
-                                  gmx_hw_opt_t *hw_opt,
-                                  const char *nbpu_opt,
-                                  t_inputrec *ir,
-                                  const gmx_mtop_t *mtop,
-                                  matrix box,
-                                  gmx_bool *bUseGPU)
-{
-    /* Here we only check for GPU usage on the MPI master process,
-     * as here we don't know how many GPUs we will use yet.
-     * We check for a GPU on all processes later.
-     */
-    *bUseGPU = hwinfo->bCanUseGPU || (getenv("GMX_EMULATE_GPU") != NULL);
-
-    if (ir->verletbuf_drift > 0)
-    {
-        /* Update the Verlet buffer size for the current run setup */
-        verletbuf_list_setup_t ls;
-        real rlist_new;
-
-        /* Here we assume CPU acceleration is on. But as currently
-         * calc_verlet_buffer_size gives the same results for 4x8 and 4x4
-         * and 4x2 gives a larger buffer than 4x4, this is ok.
-         */
-        verletbuf_get_list_setup(*bUseGPU,&ls);
-
-        calc_verlet_buffer_size(mtop,det(box),ir,
-                                ir->verletbuf_drift,&ls,
-                                NULL,&rlist_new);
-        if (rlist_new != ir->rlist)
-        {
-            if (fplog != NULL)
-            {
-                fprintf(fplog,"\nChanging rlist from %g to %g for non-bonded %dx%d atom kernels\n\n",
-                        ir->rlist,rlist_new,
-                        ls.cluster_size_i,ls.cluster_size_j);
-            }
-            ir->rlist     = rlist_new;
-            ir->rlistlong = rlist_new;
-        }
-    }
-
-    /* With GPU or emulation we should check nstlist for performance */
-    if ((EI_DYNAMICS(ir->eI) &&
-         *bUseGPU &&
-         ir->nstlist < NSTLIST_GPU_ENOUGH) ||
-        getenv(NSTLIST_ENVVAR) != NULL)
-    {
-        /* Choose a better nstlist */
-        increase_nstlist(fplog,cr,ir,mtop,box);
-    }
-}
-
-static void convert_to_verlet_scheme(FILE *fplog,
-                                     t_inputrec *ir,
-                                     gmx_mtop_t *mtop,real box_vol)
-{
-    char *conv_mesg="Converting input file with group cut-off scheme to the Verlet cut-off scheme";
-
-    md_print_warn(NULL,fplog,"%s\n",conv_mesg);
-
-    ir->cutoff_scheme   = ecutsVERLET;
-    ir->verletbuf_drift = 0.005;
-
-    if (ir->rcoulomb != ir->rvdw)
-    {
-        gmx_fatal(FARGS,"The VdW and Coulomb cut-offs are different, whereas the Verlet scheme only supports equal cut-offs");
-    }
-
-    if (ir->vdwtype == evdwUSER || EEL_USER(ir->coulombtype))
-    {
-        gmx_fatal(FARGS,"User non-bonded potentials are not (yet) supported with the Verlet scheme");
-    }
-    else if (EVDW_SWITCHED(ir->vdwtype) || EEL_SWITCHED(ir->coulombtype))
-    {
-        md_print_warn(NULL,fplog,"Converting switched or shifted interactions to a shifted potential (without force shift), this will lead to slightly different interaction potentials");
-
-        if (EVDW_SWITCHED(ir->vdwtype))
-        {
-            ir->vdwtype = evdwCUT;
-        }
-        if (EEL_SWITCHED(ir->coulombtype))
-        {
-            if (EEL_FULL(ir->coulombtype))
-            {
-                /* With full electrostatic only PME can be switched */
-                ir->coulombtype = eelPME;
-            }
-            else
-            {
-                md_print_warn(NULL,fplog,"NOTE: Replacing %s electrostatics with reaction-field with epsilon-rf=inf\n",eel_names[ir->coulombtype]);
-                ir->coulombtype = eelRF;
-                ir->epsilon_rf  = 0.0;
-            }
-        }
-
-        /* We set the target energy drift to a small number.
-         * Note that this is only for testing. For production the user
-         * should think about this and set the mdp options.
-         */
-        ir->verletbuf_drift = 1e-4;
-    }
-
-    if (inputrec2nboundeddim(ir) != 3)
-    {
-        gmx_fatal(FARGS,"Can only convert old tpr files to the Verlet cut-off scheme with 3D pbc");
-    }
-
-    if (ir->efep != efepNO || ir->implicit_solvent != eisNO)
-    {
-        gmx_fatal(FARGS,"Will not convert old tpr files to the Verlet cut-off scheme with free-energy calculations or implicit solvent");
-    }
-
-    if (EI_DYNAMICS(ir->eI) && !(EI_MD(ir->eI) && ir->etc == etcNO))
-    {
-        verletbuf_list_setup_t ls;
-
-        verletbuf_get_list_setup(FALSE,&ls);
-        calc_verlet_buffer_size(mtop,box_vol,ir,ir->verletbuf_drift,&ls,
-                                NULL,&ir->rlist);
-    }
-    else
-    {
-        ir->verletbuf_drift = -1;
-        ir->rlist           = 1.05*max(ir->rvdw,ir->rcoulomb);
-    }
-
-    gmx_mtop_remove_chargegroups(mtop);
-}
-
-/* Check the process affinity mask. If it is non-zero, something
- * else has set the affinity, and mdrun should honor that and
- * not attempt to do its own thread pinning.
- *
- * This function should be called twice. Once before the OpenMP
- * library gets initialized with bAfterOpenMPInit=FALSE (which will
- * detect affinity set by external tools like taskset), and again
- * later, after the OpenMP initialization, with bAfterOpenMPInit=TRUE
- * (which will detect affinity changes made by the OpenMP library).
- *
- * Note that this will only work on Linux, because we use a GNU
- * feature. */
-static void check_cpu_affinity_set(FILE *fplog, const t_commrec *cr,
-                                   gmx_hw_opt_t *hw_opt, int ncpus,
-                                   gmx_bool bAfterOpenmpInit)
-{
-#ifdef HAVE_SCHED_GETAFFINITY
-    cpu_set_t mask_current;
-    int       i, ret, cpu_count, cpu_set;
-    gmx_bool  bAllSet;
-
-    assert(hw_opt);
-    if (!hw_opt->bThreadPinning)
-    {
-        /* internal affinity setting is off, don't bother checking process affinity */
-        return;
-    }
-
-    CPU_ZERO(&mask_current);
-    if ((ret = sched_getaffinity(0, sizeof(cpu_set_t), &mask_current)) != 0)
-    {
-        /* failed to query affinity mask, will just return */
-        if (debug)
-        {
-            fprintf(debug, "Failed to query affinity mask (error %d)", ret);
-        }
-        return;
-    }
-
-    /* Before proceeding with the actual check, make sure that the number of
-     * detected CPUs is >= the CPUs in the current set.
-     * We need to check for CPU_COUNT as it was added only in glibc 2.6. */
-#ifdef CPU_COUNT
-    if (ncpus < CPU_COUNT(&mask_current))
-    {
-        if (debug)
-        {
-            fprintf(debug, "%d CPUs detected, but %d was returned by CPU_COUNT",
-                    ncpus, CPU_COUNT(&mask_current));
-        }
-        return;
-    }
-#endif /* CPU_COUNT */
-
-    bAllSet = TRUE;
-    for (i = 0; (i < ncpus && i < CPU_SETSIZE); i++)
-    {
-        bAllSet = bAllSet && (CPU_ISSET(i, &mask_current) != 0);
-    }
-
-    if (!bAllSet)
-    {
-        if (!bAfterOpenmpInit)
-        {
-            md_print_warn(cr, fplog,
-                          "%s detected a non-default process affinity, "
-                          "so it will not attempt to pin its threads", ShortProgram());
-        }
-        else
-        {
-            md_print_warn(cr, fplog,
-                          "%s detected a non-default process affinity, "
-                          "probably set by the OpenMP library, "
-                          "so it will not attempt to pin its threads", ShortProgram());
-        }
-        hw_opt->bThreadPinning = FALSE;
-
-        if (debug)
-        {
-            fprintf(debug, "Non-default affinity mask found, mdrun will not pin threads\n");
-        }
-    }
-    else
-    {
-        if (debug)
-        {
-            fprintf(debug, "Default affinity mask found\n");
-        }
-    }
-#endif /* HAVE_SCHED_GETAFFINITY */
-}
-
-/* Set CPU affinity. Can be important for performance.
-   On some systems (e.g. Cray) CPU Affinity is set by default.
-   But default assigning doesn't work (well) with only some ranks
-   having threads. This causes very low performance.
-   External tools have cumbersome syntax for setting affinity
-   in the case that only some ranks have threads.
-   Thus it is important that GROMACS sets the affinity internally
-   if only PME is using threads.
-*/
-static void set_cpu_affinity(FILE *fplog,
-                             const t_commrec *cr,
-                             gmx_hw_opt_t *hw_opt,
-                             int nthreads_pme,
-                             const gmx_hw_info_t *hwinfo,
-                             const t_inputrec *inputrec)
-{
-#if defined GMX_THREAD_MPI
-    /* With the number of TMPI threads equal to the number of cores
-     * we already pinned in thread-MPI, so don't pin again here.
-     */
-    if (hw_opt->nthreads_tmpi == tMPI_Thread_get_hw_number())
-    {
-        return;
-    }
-#endif
-
-#ifndef __APPLE__
-    /* If the tMPI thread affinity setting is not supported encourage the user
-     * to report it as it's either a bug or an exotic platform which we might
-     * want to support. */
-    if (tMPI_Thread_setaffinity_support() != TMPI_SETAFFINITY_SUPPORT_YES)
-    {
-        md_print_warn(NULL, fplog,
-                      "Can not set thread affinities on the current plarform. On NUMA systems this\n"
-                      "can cause performance degradation. If you think your platform should support\n"
-                      "setting affinities, contact the GROMACS developers.");
-        return;
-    }
-#endif /* __APPLE__ */
-
-    if (hw_opt->bThreadPinning)
-    {
-        int nth_affinity_set, thread_id_node, thread_id,
-            nthread_local, nthread_node, nthread_hw_max, nphyscore;
-        int offset;
-        char *env;
-
-        /* threads on this MPI process or TMPI thread */
-        if (cr->duty & DUTY_PP)
-        {
-            nthread_local = gmx_omp_nthreads_get(emntNonbonded);
-        }
-        else
-        {
-            nthread_local = gmx_omp_nthreads_get(emntPME);
-        }
-
-        /* map the current process to cores */
-        thread_id_node = 0;
-        nthread_node = nthread_local;
-#ifdef GMX_MPI
-        if (PAR(cr) || MULTISIM(cr))
-        {
-            /* We need to determine a scan of the thread counts in this
-             * compute node.
-             */
-            MPI_Comm comm_intra;
-
-            MPI_Comm_split(MPI_COMM_WORLD,gmx_hostname_num(),cr->rank_intranode,
-                           &comm_intra);
-            MPI_Scan(&nthread_local,&thread_id_node,1,MPI_INT,MPI_SUM,comm_intra);
-            /* MPI_Scan is inclusive, but here we need exclusive */
-            thread_id_node -= nthread_local;
-            /* Get the total number of threads on this physical node */
-            MPI_Allreduce(&nthread_local,&nthread_node,1,MPI_INT,MPI_SUM,comm_intra);
-            MPI_Comm_free(&comm_intra);
-        }
-#endif
-
-        offset = 0;
-        if (hw_opt->core_pinning_offset > 0)
-        {
-            offset = hw_opt->core_pinning_offset;
-            if (SIMMASTER(cr))
-            {
-                fprintf(stderr, "Applying core pinning offset %d\n", offset);
-            }
-            if (fplog)
-            {
-                fprintf(fplog, "Applying core pinning offset %d\n", offset);
-            }
-        }
-
-        /* With Intel Hyper-Threading enabled, we want to pin consecutive
-         * threads to physical cores when using more threads than physical
-         * cores or when the user requests so.
-         */
-        nthread_hw_max = hwinfo->nthreads_hw_avail;
-        nphyscore = -1;
-        if (hw_opt->bPinHyperthreading ||
-            (gmx_cpuid_x86_smt(hwinfo->cpuid_info) == GMX_CPUID_X86_SMT_ENABLED &&
-             nthread_node > nthread_hw_max/2 && getenv("GMX_DISABLE_PINHT") == NULL))
-        {
-            if (gmx_cpuid_x86_smt(hwinfo->cpuid_info) != GMX_CPUID_X86_SMT_ENABLED)
-            {
-                /* We print to stderr on all processes, as we might have
-                 * different settings on different physical nodes.
-                 */
-                if (gmx_cpuid_vendor(hwinfo->cpuid_info) != GMX_CPUID_VENDOR_INTEL)
-                {
-                    md_print_warn(NULL, fplog, "Pinning for Hyper-Threading layout requested, "
-                                  "but non-Intel CPU detected (vendor: %s)\n",
-                                  gmx_cpuid_vendor_string[gmx_cpuid_vendor(hwinfo->cpuid_info)]);
-                }
-                else
-                {
-                    md_print_warn(NULL, fplog, "Pinning for Hyper-Threading layout requested, "
-                                  "but the CPU detected does not have Intel Hyper-Threading support "
-                                  "(or it is turned off)\n");
-                }
-            }
-            nphyscore = nthread_hw_max/2;
-
-            if (SIMMASTER(cr))
-            {
-                fprintf(stderr, "Pinning to Hyper-Threading cores with %d physical cores in a compute node\n",
-                        nphyscore);
-            }
-            if (fplog)
-            {
-                fprintf(fplog, "Pinning to Hyper-Threading cores with %d physical cores in a compute node\n",
-                        nphyscore);
-            }
-        }
-
-        /* Set the per-thread affinity. In order to be able to check the success
-         * of affinity settings, we will set nth_affinity_set to 1 on threads
-         * where the affinity setting succeded and to 0 where it failed.
-         * Reducing these 0/1 values over the threads will give the total number
-         * of threads on which we succeeded.
-         */
-         nth_affinity_set = 0;
-#pragma omp parallel firstprivate(thread_id_node) num_threads(nthread_local) \
-                     reduction(+:nth_affinity_set)
-        {
-            int      core;
-            gmx_bool setaffinity_ret;
-
-            thread_id       = gmx_omp_get_thread_num();
-            thread_id_node += thread_id;
-            if (nphyscore <= 0)
-            {
-                core = offset + thread_id_node;
-            }
-            else
-            {
-                /* Lock pairs of threads to the same hyperthreaded core */
-                core = offset + thread_id_node/2 + (thread_id_node % 2)*nphyscore;
-            }
-
-            setaffinity_ret = tMPI_Thread_setaffinity_single(tMPI_Thread_self(), core);
-
-            /* store the per-thread success-values of the setaffinity */
-            nth_affinity_set = (setaffinity_ret == 0);
-
-            if (debug)
-            {
-                fprintf(debug, "On rank %2d, thread %2d, core %2d the affinity setting returned %d\n",
-                        cr->nodeid, gmx_omp_get_thread_num(), core, setaffinity_ret);
-            }
-        }
-
-        if (nth_affinity_set > nthread_local)
-        {
-            char msg[STRLEN];
-
-            sprintf(msg, "Looks like we have set affinity for more threads than "
-                    "we have (%d > %d)!\n", nth_affinity_set, nthread_local);
-            gmx_incons(msg);
-        }
-        else
-        {
-            /* check & warn if some threads failed to set their affinities */
-            if (nth_affinity_set != nthread_local)
-            {
-                char sbuf1[STRLEN], sbuf2[STRLEN];
-
-                /* sbuf1 contains rank info, while sbuf2 OpenMP thread info */
-                sbuf1[0] = sbuf2[0] = '\0';
-#ifdef GMX_MPI
-#ifdef GMX_THREAD_MPI
-                sprintf(sbuf1, "In thread-MPI thread #%d: ", cr->nodeid);
-#else /* GMX_LIB_MPI */
-                sprintf(sbuf1, "In MPI process #%d: ", cr->nodeid);
-#endif
-#endif /* GMX_MPI */
-
-                if (nthread_local > 1)
-                {
-                    sprintf(sbuf2, "of %d/%d thread%s ",
-                            nthread_local - nth_affinity_set, nthread_local,
-                            (nthread_local - nth_affinity_set) > 1 ? "s" : "");
-                }
-
-                md_print_warn(NULL, fplog,
-                              "NOTE: %sAffinity setting %sfailed.\n"
-                              "      This can cause performance degradation!",
-                              sbuf1, sbuf2);
-            }
-        }
-    }
-}
-
-
-static void check_and_update_hw_opt(gmx_hw_opt_t *hw_opt,
-                                    int cutoff_scheme)
-{
-    gmx_omp_nthreads_read_env(&hw_opt->nthreads_omp);
-
-#ifndef GMX_THREAD_MPI
-    if (hw_opt->nthreads_tot > 0)
-    {
-        gmx_fatal(FARGS,"Setting the total number of threads is only supported with thread-MPI and Gromacs was compiled without thread-MPI");
-    }
-    if (hw_opt->nthreads_tmpi > 0)
-    {
-        gmx_fatal(FARGS,"Setting the number of thread-MPI threads is only supported with thread-MPI and Gromacs was compiled without thread-MPI");
-    }
-#endif
-
-    if (hw_opt->nthreads_tot > 0 && hw_opt->nthreads_omp_pme <= 0)
-    {
-        /* We have the same number of OpenMP threads for PP and PME processes,
-         * thus we can perform several consistency checks.
-         */
-        if (hw_opt->nthreads_tmpi > 0 &&
-            hw_opt->nthreads_omp > 0 &&
-            hw_opt->nthreads_tot != hw_opt->nthreads_tmpi*hw_opt->nthreads_omp)
-        {
-            gmx_fatal(FARGS,"The total number of threads requested (%d) does not match the thread-MPI threads (%d) times the OpenMP threads (%d) requested",
-                      hw_opt->nthreads_tot,hw_opt->nthreads_tmpi,hw_opt->nthreads_omp);
-        }
-
-        if (hw_opt->nthreads_tmpi > 0 &&
-            hw_opt->nthreads_tot % hw_opt->nthreads_tmpi != 0)
-        {
-            gmx_fatal(FARGS,"The total number of threads requested (%d) is not divisible by the number of thread-MPI threads requested (%d)",
-                      hw_opt->nthreads_tot,hw_opt->nthreads_tmpi);
-        }
-
-        if (hw_opt->nthreads_omp > 0 &&
-            hw_opt->nthreads_tot % hw_opt->nthreads_omp != 0)
-        {
-            gmx_fatal(FARGS,"The total number of threads requested (%d) is not divisible by the number of OpenMP threads requested (%d)",
-                      hw_opt->nthreads_tot,hw_opt->nthreads_omp);
-        }
-
-        if (hw_opt->nthreads_tmpi > 0 &&
-            hw_opt->nthreads_omp <= 0)
-        {
-            hw_opt->nthreads_omp = hw_opt->nthreads_tot/hw_opt->nthreads_tmpi;
-        }
-    }
-
-#ifndef GMX_OPENMP
-    if (hw_opt->nthreads_omp > 1)
-    {
-        gmx_fatal(FARGS,"OpenMP threads are requested, but Gromacs was compiled without OpenMP support");
-    }
-#endif
-
-    if (cutoff_scheme == ecutsGROUP)
-    {
-        /* We only have OpenMP support for PME only nodes */
-        if (hw_opt->nthreads_omp > 1)
-        {
-            gmx_fatal(FARGS,"OpenMP threads have been requested with cut-off scheme %s, but these are only supported with cut-off scheme %s",
-                      ecutscheme_names[cutoff_scheme],
-                      ecutscheme_names[ecutsVERLET]);
-        }
-        hw_opt->nthreads_omp = 1;
-    }
-
-    if (hw_opt->nthreads_omp_pme > 0 && hw_opt->nthreads_omp <= 0)
-    {
-        gmx_fatal(FARGS,"You need to specify -ntomp in addition to -ntomp_pme");
-    }
-
-    if (hw_opt->nthreads_tot == 1)
-    {
-        hw_opt->nthreads_tmpi = 1;
-
-        if (hw_opt->nthreads_omp > 1)
-        {
-            gmx_fatal(FARGS,"You requested %d OpenMP threads with %d total threads",
-                      hw_opt->nthreads_tmpi,hw_opt->nthreads_tot);
-        }
-        hw_opt->nthreads_omp = 1;
-    }
-
-    if (hw_opt->nthreads_omp_pme <= 0 && hw_opt->nthreads_omp > 0)
-    {
-        hw_opt->nthreads_omp_pme = hw_opt->nthreads_omp;
-    }
-
-    if (debug)
-    {
-        fprintf(debug,"hw_opt: nt %d ntmpi %d ntomp %d ntomp_pme %d gpu_id '%s'\n",
-                hw_opt->nthreads_tot,
-                hw_opt->nthreads_tmpi,
-                hw_opt->nthreads_omp,
-                hw_opt->nthreads_omp_pme,
-                hw_opt->gpu_id!=NULL ? hw_opt->gpu_id : "");
-                
-    }
-}
-
-
-/* Override the value in inputrec with value passed on the command line (if any) */
-static void override_nsteps_cmdline(FILE *fplog,
-                                    int nsteps_cmdline,
-                                    t_inputrec *ir,
-                                    const t_commrec *cr)
-{
-    assert(ir);
-    assert(cr);
-
-    /* override with anything else than the default -2 */
-    if (nsteps_cmdline > -2)
-    {
-        char stmp[STRLEN];
-
-        ir->nsteps = nsteps_cmdline;
-        if (EI_DYNAMICS(ir->eI))
-        {
-            sprintf(stmp, "Overriding nsteps with value passed on the command line: %d steps, %.3f ps",
-                    nsteps_cmdline, nsteps_cmdline*ir->delta_t);
-        }
-        else
-        {
-            sprintf(stmp, "Overriding nsteps with value passed on the command line: %d steps",
-                    nsteps_cmdline);
-        }
-
-        md_print_warn(cr, fplog, "%s\n", stmp);
-    }
-}
-
-/* Data structure set by SIMMASTER which needs to be passed to all nodes
- * before the other nodes have read the tpx file and called gmx_detect_hardware.
- */
-typedef struct {
-    int cutoff_scheme; /* The cutoff scheme from inputrec_t */
-    gmx_bool bUseGPU;       /* Use GPU or GPU emulation          */
-} master_inf_t;
-
-int mdrunner(gmx_hw_opt_t *hw_opt,
-             FILE *fplog,t_commrec *cr,int nfile,
-             const t_filenm fnm[], const output_env_t oenv, gmx_bool bVerbose,
-             gmx_bool bCompact, int nstglobalcomm,
-             ivec ddxyz,int dd_node_order,real rdd,real rconstr,
-             const char *dddlb_opt,real dlb_scale,
-             const char *ddcsx,const char *ddcsy,const char *ddcsz,
-             const char *nbpu_opt,
-             int nsteps_cmdline, int nstepout,int resetstep,
-             int nmultisim,int repl_ex_nst,int repl_ex_nex,
-             int repl_ex_seed, real pforce,real cpt_period,real max_hours,
-             const char *deviceOptions, unsigned long Flags)
-{
-    gmx_bool   bForceUseGPU,bTryUseGPU;
-    double     nodetime=0,realtime;
-    t_inputrec *inputrec;
-    t_state    *state=NULL;
-    matrix     box;
-    gmx_ddbox_t ddbox={0};
-    int        npme_major,npme_minor;
-    real       tmpr1,tmpr2;
-    t_nrnb     *nrnb;
-    gmx_mtop_t *mtop=NULL;
-    t_mdatoms  *mdatoms=NULL;
-    t_forcerec *fr=NULL;
-    t_fcdata   *fcd=NULL;
-    real       ewaldcoeff=0;
-    gmx_pme_t  *pmedata=NULL;
-    gmx_vsite_t *vsite=NULL;
-    gmx_constr_t constr;
-    int        i,m,nChargePerturbed=-1,status,nalloc;
-    char       *gro;
-    gmx_wallcycle_t wcycle;
-    gmx_bool       bReadRNG,bReadEkin;
-    int        list;
-    gmx_runtime_t runtime;
-    int        rc;
-    gmx_large_int_t reset_counters;
-    gmx_edsam_t ed=NULL;
-    t_commrec   *cr_old=cr; 
-    int         nthreads_pme=1;
-    int         nthreads_pp=1;
-    gmx_membed_t membed=NULL;
-    gmx_hw_info_t *hwinfo=NULL;
-    master_inf_t minf={-1,FALSE};
-
-    /* CAUTION: threads may be started later on in this function, so
-       cr doesn't reflect the final parallel state right now */
-    snew(inputrec,1);
-    snew(mtop,1);
-    
-    if (Flags & MD_APPENDFILES) 
-    {
-        fplog = NULL;
-    }
-
-    bForceUseGPU = (strncmp(nbpu_opt, "gpu", 3) == 0);
-    bTryUseGPU   = (strncmp(nbpu_opt, "auto", 4) == 0) || bForceUseGPU;
-
-    snew(state,1);
-    if (SIMMASTER(cr)) 
-    {
-        /* Read (nearly) all data required for the simulation */
-        read_tpx_state(ftp2fn(efTPX,nfile,fnm),inputrec,state,NULL,mtop);
-
-        if (inputrec->cutoff_scheme != ecutsVERLET &&
-            ((Flags & MD_TESTVERLET) || getenv("GMX_VERLET_SCHEME") != NULL))
-        {
-            convert_to_verlet_scheme(fplog,inputrec,mtop,det(state->box));
-        }
-
-        /* Detect hardware, gather information. With tMPI only thread 0 does it
-         * and after threads are started broadcasts hwinfo around. */
-        snew(hwinfo, 1);
-        gmx_detect_hardware(fplog, hwinfo, cr,
-                            bForceUseGPU, bTryUseGPU, hw_opt->gpu_id);
-
-        minf.cutoff_scheme = inputrec->cutoff_scheme;
-        minf.bUseGPU       = FALSE;
-
-        if (inputrec->cutoff_scheme == ecutsVERLET)
-        {
-            prepare_verlet_scheme(fplog,hwinfo,cr,hw_opt,nbpu_opt,
-                                  inputrec,mtop,state->box,
-                                  &minf.bUseGPU);
-        }
-        else if (hwinfo->bCanUseGPU)
-        {
-            md_print_warn(cr,fplog,
-                          "NOTE: GPU(s) found, but the current simulation can not use GPUs\n"
-                          "      To use a GPU, set the mdp option: cutoff-scheme = Verlet\n"
-                          "      (for quick performance testing you can use the -testverlet option)\n");
-
-            if (bForceUseGPU)
-            {
-                gmx_fatal(FARGS,"GPU requested, but can't be used without cutoff-scheme=Verlet");
-            }
-        }
-    }
-#ifndef GMX_THREAD_MPI
-    if (PAR(cr))
-    {
-        gmx_bcast_sim(sizeof(minf),&minf,cr);
-    }
-#endif
-    if (minf.bUseGPU && cr->npmenodes == -1)
-    {
-        /* Don't automatically use PME-only nodes with GPUs */
-        cr->npmenodes = 0;
-    }
-
-    /* Check for externally set OpenMP affinity and turn off internal
-     * pinning if any is found. We need to do this check early to tell
-     * thread-MPI whether it should do pinning when spawning threads.
-     */
-    gmx_omp_check_thread_affinity(fplog, cr, hw_opt);
-
-#ifdef GMX_THREAD_MPI
-    /* With thread-MPI inputrec is only set here on the master thread */
-    if (SIMMASTER(cr))
-#endif
-    {
-        check_and_update_hw_opt(hw_opt,minf.cutoff_scheme);
-
-#ifdef GMX_THREAD_MPI
-        /* Early check for externally set process affinity. Can't do over all
-         * MPI processes because hwinfo is not available everywhere, but with
-         * thread-MPI it's needed as pinning might get turned off which needs
-         * to be known before starting thread-MPI. */
-        check_cpu_affinity_set(fplog,
-                               NULL,
-                               hw_opt, hwinfo->nthreads_hw_avail, FALSE);
-#endif
-
-#ifdef GMX_THREAD_MPI
-        if (cr->npmenodes > 0 && hw_opt->nthreads_tmpi <= 0)
-        {
-            gmx_fatal(FARGS,"You need to explicitly specify the number of MPI threads (-ntmpi) when using separate PME nodes");
-        }
-#endif
-
-        if (hw_opt->nthreads_omp_pme != hw_opt->nthreads_omp &&
-            cr->npmenodes <= 0)
-        {
-            gmx_fatal(FARGS,"You need to explicitly specify the number of PME nodes (-npme) when using different number of OpenMP threads for PP and PME nodes");
-        }
-    }
-
-#ifdef GMX_THREAD_MPI
-    if (SIMMASTER(cr))
-    {
-        /* NOW the threads will be started: */
-        hw_opt->nthreads_tmpi = get_nthreads_mpi(hwinfo,
-                                                 hw_opt,
-                                                 inputrec, mtop,
-                                                 cr, fplog);
-        if (hw_opt->nthreads_tot > 0 && hw_opt->nthreads_omp <= 0)
-        {
-            hw_opt->nthreads_omp = hw_opt->nthreads_tot/hw_opt->nthreads_tmpi;
-        }
-
-        if (hw_opt->nthreads_tmpi > 1)
-        {
-            /* now start the threads. */
-            cr=mdrunner_start_threads(hw_opt, fplog, cr_old, nfile, fnm, 
-                                      oenv, bVerbose, bCompact, nstglobalcomm, 
-                                      ddxyz, dd_node_order, rdd, rconstr, 
-                                      dddlb_opt, dlb_scale, ddcsx, ddcsy, ddcsz,
-                                      nbpu_opt,
-                                      nsteps_cmdline, nstepout, resetstep, nmultisim, 
-                                      repl_ex_nst, repl_ex_nex, repl_ex_seed, pforce,
-                                      cpt_period, max_hours, deviceOptions, 
-                                      Flags);
-            /* the main thread continues here with a new cr. We don't deallocate
-               the old cr because other threads may still be reading it. */
-            if (cr == NULL)
-            {
-                gmx_comm("Failed to spawn threads");
-            }
-        }
-    }
-#endif
-    /* END OF CAUTION: cr is now reliable */
-
-    /* g_membed initialisation *
-     * Because we change the mtop, init_membed is called before the init_parallel *
-     * (in case we ever want to make it run in parallel) */
-    if (opt2bSet("-membed",nfile,fnm))
-    {
-        if (MASTER(cr))
-        {
-            fprintf(stderr,"Initializing membed");
-        }
-        membed = init_membed(fplog,nfile,fnm,mtop,inputrec,state,cr,&cpt_period);
-    }
-
-    if (PAR(cr))
-    {
-        /* now broadcast everything to the non-master nodes/threads: */
-        init_parallel(fplog, cr, inputrec, mtop);
-
-        /* This check needs to happen after get_nthreads_mpi() */
-        if (inputrec->cutoff_scheme == ecutsVERLET && (Flags & MD_PARTDEC))
-        {
-            gmx_fatal_collective(FARGS,cr,NULL,
-                                 "The Verlet cut-off scheme is not supported with particle decomposition.\n"
-                                 "You can achieve the same effect as particle decomposition by running in parallel using only OpenMP threads.");
-        }
-    }
-    if (fplog != NULL)
-    {
-        pr_inputrec(fplog,0,"Input Parameters",inputrec,FALSE);
-    }
-
-#if defined GMX_THREAD_MPI
-    /* With tMPI we detected on thread 0 and we'll just pass the hwinfo pointer
-     * to the other threads  -- slightly uncool, but works fine, just need to
-     * make sure that the data doesn't get freed twice. */
-    if (cr->nnodes > 1)
-    {
-        if (!SIMMASTER(cr))
-        {
-            snew(hwinfo, 1);
-        }
-        gmx_bcast(sizeof(&hwinfo), &hwinfo, cr);
-    }
-#else
-    if (PAR(cr) && !SIMMASTER(cr))
-    {
-        /* now we have inputrec on all nodes, can run the detection */
-        /* TODO: perhaps it's better to propagate within a node instead? */
-        snew(hwinfo, 1);
-        gmx_detect_hardware(fplog, hwinfo, cr,
-                                 bForceUseGPU, bTryUseGPU, hw_opt->gpu_id);
-    }
-
-    /* Now do the affinity check with MPI/no-MPI (done earlier with thread-MPI). */
-    check_cpu_affinity_set(fplog, cr,
-                           hw_opt, hwinfo->nthreads_hw_avail, FALSE);
-#endif
-
-    /* now make sure the state is initialized and propagated */
-    set_state_entries(state,inputrec,cr->nnodes);
-
-    /* remove when vv and rerun works correctly! */
-    if (PAR(cr) && EI_VV(inputrec->eI) && ((Flags & MD_RERUN) || (Flags & MD_RERUN_VSITE)))
-    {
-        gmx_fatal(FARGS,
-                  "Currently can't do velocity verlet with rerun in parallel.");
-    }
-
-    /* A parallel command line option consistency check that we can
-       only do after any threads have started. */
-    if (!PAR(cr) &&
-        (ddxyz[XX] > 1 || ddxyz[YY] > 1 || ddxyz[ZZ] > 1 || cr->npmenodes > 0))
-    {
-        gmx_fatal(FARGS,
-                  "The -dd or -npme option request a parallel simulation, "
-#ifndef GMX_MPI
-                  "but %s was compiled without threads or MPI enabled"
-#else
-#ifdef GMX_THREAD_MPI
-                  "but the number of threads (option -nt) is 1"
-#else
-                  "but %s was not started through mpirun/mpiexec or only one process was requested through mpirun/mpiexec"
-#endif
-#endif
-                  , ShortProgram()
-            );
-    }
-
-    if ((Flags & MD_RERUN) &&
-        (EI_ENERGY_MINIMIZATION(inputrec->eI) || eiNM == inputrec->eI))
-    {
-        gmx_fatal(FARGS, "The .mdp file specified an energy mininization or normal mode algorithm, and these are not compatible with mdrun -rerun");
-    }
-
-    if (can_use_allvsall(inputrec,mtop,TRUE,cr,fplog) && PAR(cr))
-    {
-        /* All-vs-all loops do not work with domain decomposition */
-        Flags |= MD_PARTDEC;
-    }
-
-    if (!EEL_PME(inputrec->coulombtype) || (Flags & MD_PARTDEC))
-    {
-        if (cr->npmenodes > 0)
-        {
-            if (!EEL_PME(inputrec->coulombtype))
-            {
-                gmx_fatal_collective(FARGS,cr,NULL,
-                                     "PME nodes are requested, but the system does not use PME electrostatics");
-            }
-            if (Flags & MD_PARTDEC)
-            {
-                gmx_fatal_collective(FARGS,cr,NULL,
-                                     "PME nodes are requested, but particle decomposition does not support separate PME nodes");
-            }
-        }
-
-        cr->npmenodes = 0;
-    }
-
-#ifdef GMX_FAHCORE
-    fcRegisterSteps(inputrec->nsteps,inputrec->init_step);
-#endif
-
-    /* NMR restraints must be initialized before load_checkpoint,
-     * since with time averaging the history is added to t_state.
-     * For proper consistency check we therefore need to extend
-     * t_state here.
-     * So the PME-only nodes (if present) will also initialize
-     * the distance restraints.
-     */
-    snew(fcd,1);
-
-    /* This needs to be called before read_checkpoint to extend the state */
-    init_disres(fplog,mtop,inputrec,cr,Flags & MD_PARTDEC,fcd,state);
-
-    if (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0)
-    {
-        if (PAR(cr) && !(Flags & MD_PARTDEC))
-        {
-            gmx_fatal(FARGS,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
-        }
-        /* Orientation restraints */
-        if (MASTER(cr))
-        {
-            init_orires(fplog,mtop,state->x,inputrec,cr->ms,&(fcd->orires),
-                        state);
-        }
-    }
-
-    if (DEFORM(*inputrec))
-    {
-        /* Store the deform reference box before reading the checkpoint */
-        if (SIMMASTER(cr))
-        {
-            copy_mat(state->box,box);
-        }
-        if (PAR(cr))
-        {
-            gmx_bcast(sizeof(box),box,cr);
-        }
-        /* Because we do not have the update struct available yet
-         * in which the reference values should be stored,
-         * we store them temporarily in static variables.
-         * This should be thread safe, since they are only written once
-         * and with identical values.
-         */
-#ifdef GMX_THREAD_MPI
-        tMPI_Thread_mutex_lock(&deform_init_box_mutex);
-#endif
-        deform_init_init_step_tpx = inputrec->init_step;
-        copy_mat(box,deform_init_box_tpx);
-#ifdef GMX_THREAD_MPI
-        tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
-#endif
-    }
-
-    if (opt2bSet("-cpi",nfile,fnm)) 
-    {
-        /* Check if checkpoint file exists before doing continuation.
-         * This way we can use identical input options for the first and subsequent runs...
-         */
-        if( gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr) )
-        {
-            load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
-                            cr,Flags & MD_PARTDEC,ddxyz,
-                            inputrec,state,&bReadRNG,&bReadEkin,
-                            (Flags & MD_APPENDFILES),
-                            (Flags & MD_APPENDFILESSET));
-            
-            if (bReadRNG)
-            {
-                Flags |= MD_READ_RNG;
-            }
-            if (bReadEkin)
-            {
-                Flags |= MD_READ_EKIN;
-            }
-        }
-    }
-
-    if (((MASTER(cr) || (Flags & MD_SEPPOT)) && (Flags & MD_APPENDFILES))
-#ifdef GMX_THREAD_MPI
-        /* With thread MPI only the master node/thread exists in mdrun.c,
-         * therefore non-master nodes need to open the "seppot" log file here.
-         */
-        || (!MASTER(cr) && (Flags & MD_SEPPOT))
-#endif
-        )
-    {
-        gmx_log_open(ftp2fn(efLOG,nfile,fnm),cr,!(Flags & MD_SEPPOT),
-                             Flags,&fplog);
-    }
-
-    /* override nsteps with value from cmdline */
-    override_nsteps_cmdline(fplog, nsteps_cmdline, inputrec, cr);
-
-    if (SIMMASTER(cr)) 
-    {
-        copy_mat(state->box,box);
-    }
-
-    if (PAR(cr)) 
-    {
-        gmx_bcast(sizeof(box),box,cr);
-    }
-
-    /* Essential dynamics */
-    if (opt2bSet("-ei",nfile,fnm))
-    {
-        /* Open input and output files, allocate space for ED data structure */
-        ed = ed_open(nfile,fnm,Flags,cr);
-    }
-
-    if (PAR(cr) && !((Flags & MD_PARTDEC) ||
-                     EI_TPI(inputrec->eI) ||
-                     inputrec->eI == eiNM))
-    {
-        cr->dd = init_domain_decomposition(fplog,cr,Flags,ddxyz,rdd,rconstr,
-                                           dddlb_opt,dlb_scale,
-                                           ddcsx,ddcsy,ddcsz,
-                                           mtop,inputrec,
-                                           box,state->x,
-                                           &ddbox,&npme_major,&npme_minor);
-
-        make_dd_communicators(fplog,cr,dd_node_order);
-
-        /* Set overallocation to avoid frequent reallocation of arrays */
-        set_over_alloc_dd(TRUE);
-    }
-    else
-    {
-        /* PME, if used, is done on all nodes with 1D decomposition */
-        cr->npmenodes = 0;
-        cr->duty = (DUTY_PP | DUTY_PME);
-        npme_major = 1;
-        npme_minor = 1;
-        if (!EI_TPI(inputrec->eI))
-        {
-            npme_major = cr->nnodes;
-        }
-        
-        if (inputrec->ePBC == epbcSCREW)
-        {
-            gmx_fatal(FARGS,
-                      "pbc=%s is only implemented with domain decomposition",
-                      epbc_names[inputrec->ePBC]);
-        }
-    }
-
-    if (PAR(cr))
-    {
-        /* After possible communicator splitting in make_dd_communicators.
-         * we can set up the intra/inter node communication.
-         */
-        gmx_setup_nodecomm(fplog,cr);
-    }
-
-    /* Initialize per-physical-node MPI process/thread ID and counters. */
-    gmx_init_intranode_counters(cr);
-
-#ifdef GMX_MPI
-    md_print_info(cr,fplog,"Using %d MPI %s\n",
-                  cr->nnodes,
-#ifdef GMX_THREAD_MPI
-                  cr->nnodes==1 ? "thread" : "threads"
-#else
-                  cr->nnodes==1 ? "process" : "processes"
-#endif
-                  );
-    fflush(stderr);
-#endif
-
-    gmx_omp_nthreads_init(fplog, cr,
-                          hwinfo->nthreads_hw_avail,
-                          hw_opt->nthreads_omp,
-                          hw_opt->nthreads_omp_pme,
-                          (cr->duty & DUTY_PP) == 0,
-                          inputrec->cutoff_scheme == ecutsVERLET);
-
-    gmx_check_hw_runconf_consistency(fplog, hwinfo, cr, hw_opt->nthreads_tmpi, minf.bUseGPU);
-
-    /* getting number of PP/PME threads
-       PME: env variable should be read only on one node to make sure it is 
-       identical everywhere;
-     */
-    /* TODO nthreads_pp is only used for pinning threads.
-     * This is a temporary solution until we have a hw topology library.
-     */
-    nthreads_pp  = gmx_omp_nthreads_get(emntNonbonded);
-    nthreads_pme = gmx_omp_nthreads_get(emntPME);
-
-    wcycle = wallcycle_init(fplog,resetstep,cr,nthreads_pp,nthreads_pme);
-
-    if (PAR(cr))
-    {
-        /* Master synchronizes its value of reset_counters with all nodes 
-         * including PME only nodes */
-        reset_counters = wcycle_get_reset_counters(wcycle);
-        gmx_bcast_sim(sizeof(reset_counters),&reset_counters,cr);
-        wcycle_set_reset_counters(wcycle, reset_counters);
-    }
-
-    snew(nrnb,1);
-    if (cr->duty & DUTY_PP)
-    {
-        /* For domain decomposition we allocate dynamically
-         * in dd_partition_system.
-         */
-        if (DOMAINDECOMP(cr))
-        {
-            bcast_state_setup(cr,state);
-        }
-        else
-        {
-            if (PAR(cr))
-            {
-                bcast_state(cr,state,TRUE);
-            }
-        }
-
-        /* Initiate forcerecord */
-        fr = mk_forcerec();
-        fr->hwinfo = hwinfo;
-        init_forcerec(fplog,oenv,fr,fcd,inputrec,mtop,cr,box,FALSE,
-                      opt2fn("-table",nfile,fnm),
-                      opt2fn("-tabletf",nfile,fnm),
-                      opt2fn("-tablep",nfile,fnm),
-                      opt2fn("-tableb",nfile,fnm),
-                      nbpu_opt,
-                      FALSE,pforce);
-
-        /* version for PCA_NOT_READ_NODE (see md.c) */
-        /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
-          "nofile","nofile","nofile","nofile",FALSE,pforce);
-          */        
-        fr->bSepDVDL = ((Flags & MD_SEPPOT) == MD_SEPPOT);
-
-        /* Initialize QM-MM */
-        if(fr->bQMMM)
-        {
-            init_QMMMrec(cr,box,mtop,inputrec,fr);
-        }
-
-        /* Initialize the mdatoms structure.
-         * mdatoms is not filled with atom data,
-         * as this can not be done now with domain decomposition.
-         */
-        mdatoms = init_mdatoms(fplog,mtop,inputrec->efep!=efepNO);
-
-        /* Initialize the virtual site communication */
-        vsite = init_vsite(mtop,cr,FALSE);
-
-        calc_shifts(box,fr->shift_vec);
-
-        /* With periodic molecules the charge groups should be whole at start up
-         * and the virtual sites should not be far from their proper positions.
-         */
-        if (!inputrec->bContinuation && MASTER(cr) &&
-            !(inputrec->ePBC != epbcNONE && inputrec->bPeriodicMols))
-        {
-            /* Make molecules whole at start of run */
-            if (fr->ePBC != epbcNONE)
-            {
-                do_pbc_first_mtop(fplog,inputrec->ePBC,box,mtop,state->x);
-            }
-            if (vsite)
-            {
-                /* Correct initial vsite positions are required
-                 * for the initial distribution in the domain decomposition
-                 * and for the initial shell prediction.
-                 */
-                construct_vsites_mtop(fplog,vsite,mtop,state->x);
-            }
-        }
-
-        if (EEL_PME(fr->eeltype))
-        {
-            ewaldcoeff = fr->ewaldcoeff;
-            pmedata = &fr->pmedata;
-        }
-        else
-        {
-            pmedata = NULL;
-        }
-    }
-    else
-    {
-        /* This is a PME only node */
-
-        /* We don't need the state */
-        done_state(state);
-
-        ewaldcoeff = calc_ewaldcoeff(inputrec->rcoulomb, inputrec->ewald_rtol);
-        snew(pmedata,1);
-    }
-
-    /* Before setting affinity, check whether the affinity has changed
-     * - which indicates that probably the OpenMP library has changed it since
-     * we first checked). */
-    check_cpu_affinity_set(fplog, cr, hw_opt, hwinfo->nthreads_hw_avail, TRUE);
-
-    /* Set the CPU affinity */
-    set_cpu_affinity(fplog,cr,hw_opt,nthreads_pme,hwinfo,inputrec);
-
-    /* Initiate PME if necessary,
-     * either on all nodes or on dedicated PME nodes only. */
-    if (EEL_PME(inputrec->coulombtype))
-    {
-        if (mdatoms)
-        {
-            nChargePerturbed = mdatoms->nChargePerturbed;
-        }
-        if (cr->npmenodes > 0)
-        {
-            /* The PME only nodes need to know nChargePerturbed */
-            gmx_bcast_sim(sizeof(nChargePerturbed),&nChargePerturbed,cr);
-        }
-
-        if (cr->duty & DUTY_PME)
-        {
-            status = gmx_pme_init(pmedata,cr,npme_major,npme_minor,inputrec,
-                                  mtop ? mtop->natoms : 0,nChargePerturbed,
-                                  (Flags & MD_REPRODUCIBLE),nthreads_pme);
-            if (status != 0) 
-            {
-                gmx_fatal(FARGS,"Error %d initializing PME",status);
-            }
-        }
-    }
-
-
-    if (integrator[inputrec->eI].func == do_md
-        ||
-        integrator[inputrec->eI].func == do_md_openmm
-        )
-    {
-        /* Turn on signal handling on all nodes */
-        /*
-         * (A user signal from the PME nodes (if any)
-         * is communicated to the PP nodes.
-         */
-        signal_handler_install();
-    }
-
-    if (cr->duty & DUTY_PP)
-    {
-        if (inputrec->ePull != epullNO)
-        {
-            /* Initialize pull code */
-            init_pull(fplog,inputrec,nfile,fnm,mtop,cr,oenv, inputrec->fepvals->init_lambda,
-                      EI_DYNAMICS(inputrec->eI) && MASTER(cr),Flags);
-        }
-        
-        if (inputrec->bRot)
-        {
-           /* Initialize enforced rotation code */
-           init_rot(fplog,inputrec,nfile,fnm,cr,state->x,box,mtop,oenv,
-                    bVerbose,Flags);
-        }
-
-        constr = init_constraints(fplog,mtop,inputrec,ed,state,cr);
-
-        if (DOMAINDECOMP(cr))
-        {
-            dd_init_bondeds(fplog,cr->dd,mtop,vsite,constr,inputrec,
-                            Flags & MD_DDBONDCHECK,fr->cginfo_mb);
-
-            set_dd_parameters(fplog,cr->dd,dlb_scale,inputrec,fr,&ddbox);
-
-            setup_dd_grid(fplog,cr->dd);
-        }
-
-        /* Now do whatever the user wants us to do (how flexible...) */
-        integrator[inputrec->eI].func(fplog,cr,nfile,fnm,
-                                      oenv,bVerbose,bCompact,
-                                      nstglobalcomm,
-                                      vsite,constr,
-                                      nstepout,inputrec,mtop,
-                                      fcd,state,
-                                      mdatoms,nrnb,wcycle,ed,fr,
-                                      repl_ex_nst,repl_ex_nex,repl_ex_seed,
-                                      membed,
-                                      cpt_period,max_hours,
-                                      deviceOptions,
-                                      Flags,
-                                      &runtime);
-
-        if (inputrec->ePull != epullNO)
-        {
-            finish_pull(fplog,inputrec->pull);
-        }
-        
-        if (inputrec->bRot)
-        {
-            finish_rot(inputrec->rot);
-        }
-
-    } 
-    else 
-    {
-        /* do PME only */
-        gmx_pmeonly(*pmedata,cr,nrnb,wcycle,ewaldcoeff,FALSE,inputrec);
-    }
-
-    if (EI_DYNAMICS(inputrec->eI) || EI_TPI(inputrec->eI))
-    {
-        /* Some timing stats */  
-        if (SIMMASTER(cr))
-        {
-            if (runtime.proc == 0)
-            {
-                runtime.proc = runtime.real;
-            }
-        }
-        else
-        {
-            runtime.real = 0;
-        }
-    }
-
-    wallcycle_stop(wcycle,ewcRUN);
-
-    /* Finish up, write some stuff
-     * if rerunMD, don't write last frame again 
-     */
-    finish_run(fplog,cr,ftp2fn(efSTO,nfile,fnm),
-               inputrec,nrnb,wcycle,&runtime,
-               fr != NULL && fr->nbv != NULL && fr->nbv->bUseGPU ?
-                 nbnxn_cuda_get_timings(fr->nbv->cu_nbv) : NULL,
-               nthreads_pp, 
-               EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
-
-    if ((cr->duty & DUTY_PP) && fr->nbv != NULL && fr->nbv->bUseGPU)
-    {
-        char gpu_err_str[STRLEN];
-
-        /* free GPU memory and uninitialize GPU (by destroying the context) */
-        nbnxn_cuda_free(fplog, fr->nbv->cu_nbv);
-
-        if (!free_gpu(gpu_err_str))
-        {
-            gmx_warning("On node %d failed to free GPU #%d: %s",
-                        cr->nodeid, get_current_gpu_device_id(), gpu_err_str);
-        }
-    }
-
-    if (opt2bSet("-membed",nfile,fnm))
-    {
-        sfree(membed);
-    }
-
-#ifdef GMX_THREAD_MPI
-    if (PAR(cr) && SIMMASTER(cr))
-#endif
-    {
-        gmx_hardware_info_free(hwinfo);
-    }
-
-    /* Does what it says */  
-    print_date_and_time(fplog,cr->nodeid,"Finished mdrun",&runtime);
-
-    /* Close logfile already here if we were appending to it */
-    if (MASTER(cr) && (Flags & MD_APPENDFILES))
-    {
-        gmx_log_close(fplog);
-    }  
-
-    rc=(int)gmx_get_stop_condition();
-
-#ifdef GMX_THREAD_MPI
-    /* we need to join all threads. The sub-threads join when they
-       exit this function, but the master thread needs to be told to 
-       wait for that. */
-    if (PAR(cr) && MASTER(cr))
-    {
-        tMPI_Finalize();
-    }
-#endif
-
-    return rc;
-}
index 37fa948faf5a085c357dab4b1cd482b7f6f80415..e9491b838b574c00d2db62e80427602dbb5a04aa 100644 (file)
@@ -188,10 +188,6 @@ t_mdebin *init_mdebin(ener_file_t       fp_ene,
         md->bEInd[i] = FALSE;
     }
 
-    /* Even though the OpenMM build has moved to contrib, it's not
-     * practical to move/remove this code fragment, because of the
-     * fundamental mess that is the GROMACS library structure. */
-#ifndef GMX_OPENMM
     for (i = 0; i < F_NRE; i++)
     {
         md->bEner[i] = FALSE;
@@ -303,13 +299,6 @@ t_mdebin *init_mdebin(ener_file_t       fp_ene,
             md->bEner[i] = (gmx_mtop_ftype_count(mtop, i) > 0);
         }
     }
-#else
-    /* OpenMM always produces only the following 4 energy terms */
-    md->bEner[F_EPOT] = TRUE;
-    md->bEner[F_EKIN] = TRUE;
-    md->bEner[F_ETOT] = TRUE;
-    md->bEner[F_TEMP] = TRUE;
-#endif
 
     /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
     if (ir->bAdress && !debug)
index 6c2b258e47db58ad8d7ad5b3fe58cb684fbba268..de0acc0efc313987134ad95fc6e3bb52af2e905e 100644 (file)
@@ -52,7 +52,7 @@
 const unsigned supercl_interaction_mask = ((1U << NCL_PER_SUPERCL) - 1U);
 
 /*! Interpolate Ewald coulomb force using the table through the tex_nbfp texture.
- *  Original idea: OpenMM
+ *  Original idea: from the OpenMM project
  */
 static inline __device__
 float interpolate_coulomb_force_r(float r, float scale)
index 090710ff61783b9bb163ba0eae14f9321c1dfac5..ec64f9ad314e6900619e2b4e1cf27a16565873e0 100644 (file)
@@ -5,16 +5,6 @@ set(MDRUN_SOURCES
     md.c          mdrun.cpp     membed.c
     pme_loadbal.c repl_ex.c     runner.c    xutils.c)
 
-if(GMX_OPENMM)
-    # Even though the OpenMM build has "moved to contrib", many things
-    # have be be done from within the scope of the CMakeLists.txt that
-    # builds its mdrun, and that is here
-    list(APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/src/contrib)
-    find_package(OpenMM)
-    include_directories(${CMAKE_CURRENT_SOURCE_DIR})
-    include(${CMAKE_SOURCE_DIR}/src/contrib/BuildMdrunOpenMM.cmake)
-endif(GMX_OPENMM)
-
 if(GMX_FAHCORE)
     add_library(fahcore ${MDRUN_SOURCES})
 else(GMX_FAHCORE)
@@ -24,9 +14,6 @@ else(GMX_FAHCORE)
         ${GMX_EXE_LINKER_FLAGS})
     set_target_properties(mdrun PROPERTIES OUTPUT_NAME "mdrun${GMX_BINARY_SUFFIX}"
         COMPILE_FLAGS "${OpenMP_C_FLAGS}")
-    if(GMX_OPENMM)
-        target_link_libraries(mdrun openmm_api_wrapper)
-    endif()
     install(TARGETS mdrun DESTINATION ${BIN_INSTALL_DIR} COMPONENT mdrun)
 
     # Create the custom install-mdrun target