Merge remote-tracking branch 'origin/master' into qc

CMakeLists.txt

@@ -289,7 +289,8 @@ ENDIF(WIN32)
 SET(OPENMM_BUILD_STATIC_LIB OFF CACHE BOOL "Whether to build static OpenMM libraries")
 IF(OPENMM_BUILD_STATIC_LIB)
     ADD_LIBRARY(${STATIC_TARGET} STATIC ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
-    SET_TARGET_PROPERTIES(${STATIC_TARGET} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_USE_STATIC_LIBRARIES -DOPENMM_BUILDING_STATIC_LIBRARY -DLEPTON_USE_STATIC_LIBRARIES -DLEPTON_BUILDING_STATIC_LIBRARY -DOPENMMM_VALIDATE_BUILDING_STATIC_LIBRARY -DOPENMM_VALIDATE_BUILDING_STATIC_LIBRARY")
+    SET(EXTRA_COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_USE_STATIC_LIBRARIES -DLEPTON_USE_STATIC_LIBRARIES -DPTW32_STATIC_LIB")
+    SET_TARGET_PROPERTIES(${STATIC_TARGET} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_BUILDING_STATIC_LIBRARY -DLEPTON_BUILDING_STATIC_LIBRARY -DOPENMMM_VALIDATE_BUILDING_STATIC_LIBRARY -DOPENMM_VALIDATE_BUILDING_STATIC_LIBRARY")
 ENDIF(OPENMM_BUILD_STATIC_LIB)
 
 IF(OPENMM_BUILD_C_AND_FORTRAN_WRAPPERS)

docs/usersguide/application.rst

@@ -660,7 +660,9 @@ amber99sbildn.xml  AMBER99SB plus improved side chain torsions\ :cite:`Lindorff-
 amber99sbnmr.xml   AMBER99SB with modifications to fit NMR data\ :cite:`Li2010`
 amber03.xml        AMBER03\ :cite:`Duan2003`
 amber10.xml        AMBER10
-amoeba2009.xml     AMOEBA\ :cite:`Ren2002`                                                       
+amoeba2009.xml     AMOEBA 2009\ :cite:`Ren2002`.  This force field is deprecated.  It is 
+                   recommended to use AMOEBA 2013 instead.
+amoeba2013.xml     AMOEBA 2013\ :cite:`Shi2013`
 =================  ================================================================================
 
 
@@ -675,7 +677,9 @@ files:
 File         Water Model                                 
 ===========  ============================================
 tip3p.xml    TIP3P water model\ :cite:`Jorgensen1983`  
+tip3pfb.xml  TIP3P-FB water model\ :cite:`Wang2014`    
 tip4pew.xml  TIP4P-Ew water model\ :cite:`Horn2004`    
+tip4pfb.xml  TIP4P-FB water model\ :cite:`Wang2014`    
 tip5p.xml    TIP5P water model\ :cite:`Mahoney2000`    
 spce.xml     SPC/E water model\ :cite:`Berendsen1987`  
 swm4ndp.xml  SWM4-NDP water model\ :cite:`Lamoureux2006`
@@ -692,15 +696,16 @@ the following files:
 
 .. tabularcolumns:: |l|L|
 
-=================  ==============================================================================================
+=================  =================================================================================================
 File               Implicit Solvation Model                                                                      
-=================  ==============================================================================================
+=================  =================================================================================================
 amber96_obc.xml    GBSA-OBC solvation model\ :cite:`Onufriev2004` for use with AMBER96 force field
 amber99_obc.xml    GBSA-OBC solvation model for use with AMBER99 force fields
 amber03_obc.xml    GBSA-OBC solvation model for use with AMBER03 force field
 amber10_obc.xml    GBSA-OBC solvation model for use with AMBER10 force field
-amoeba2009_gk.xml  Generalized Kirkwood solvation model\ :cite:`Schnieders2007` for use with AMOEBA force field
-=================  ==============================================================================================
+amoeba2009_gk.xml  Generalized Kirkwood solvation model\ :cite:`Schnieders2007` for use with AMOEBA 2009 force field
+amoeba2013_gk.xml  Generalized Kirkwood solvation model for use with AMOEBA 2013 force field
+=================  =================================================================================================
 
 
 For example, to use the GBSA-OBC solvation model with the Amber99SB force field,
@@ -1297,8 +1302,9 @@ water models:
 
     modeller.addSolvent(forcefield, model='tip5p')
 
-Allowed values for the :code:`model` option are 'tip3p', 'spce', 'tip4pew',
-and 'tip5p'.  Be sure to include the single quotes around the value.
+Allowed values for the :code:`model` option are 'tip3p', 'tip3pfb', 'spce', 
+'tip4pew', 'tip4pfb', and 'tip5p'.  Be sure to include the single quotes 
+around the value.
 
 Another option is to add extra ion pairs to give a desired total ionic strength.
 For example:
@@ -2180,6 +2186,49 @@ second atom has class OS and the third has class P:
 
     <Proper class1="" class2="OS" class3="P" class4="" per="3" phase="0.0" k="0.66944"/>
 
+<CustomNonbondedForce>
+===============
+
+To add a CustomNonbondedForce to the System, include a tag that looks like this:
+
+.. code-block:: xml
+
+    <CustomNonbondedForce energy="scale*epsilon1*epsilon2*((sigma1+sigma2)/r)^12" bondCutoff="3">
+     <GlobalParameter name="scale" defaultValue="1"/>
+     <PerParticleParameter name="sigma"/>
+     <PerParticleParameter name="epsilon"/>
+     <Atom type="0" sigma="0.3249" epsilon="0.7112"/>
+     <Atom type="1" sigma="0.1069" epsilon="0.0656"/>
+     <Atom type="2" sigma="0.3399" epsilon="0.4577"/>
+     ...
+    </CustomNonbondedForce>
+
+The energy expression for the CustomNonbondedForce is specified by the
+:code:`energy` attribute.  This is a mathematical expression that gives the
+energy of each pairwise interaction as a function of the distance *r*\ .  It
+also may depend on an arbitrary list of global or per-particle parameters.  Use
+a :code:`<GlobalParameter>` tag to define a global parameter, and a
+:code:`<PerParticleParameter>` tag to define a per-particle parameter.
+
+Exclusions are created automatically based on the :code:`bondCutoff` attribute.
+After setting the nonbonded parameters for all atoms, the force field calls
+:code:`createExclusionsFromBonds()` on the CustomNonbondedForce, passing in this
+value as its argument.  To avoid creating exclusions, set :code:`bondCutoff` to 0.
+
+Each :code:`<Atom>` tag specifies the parameters for one atom type
+(specified with the :code:`type` attribute) or atom class (specified with
+the :code:`class` attribute).  It is fine to mix these two methods, having
+some tags specify a type and others specify a class.  However you do it, you
+must make sure that a unique set of parameters is defined for every atom type.
+The remaining attributes are the values to use for the per-atom parameters. All
+per-atom parameters must be specified for every :code:`<Atom>` tag, and the
+attribute name must match the name of the parameter.  For instance, if there is
+a per-atom parameter with the name “radius”, then every :code:`<Atom>` tag
+must include an attribute called :code:`radius`\ .
+
+CustomNonbondedForce also allows you to define tabulated functions.  See section
+:ref:`tabulated-functions` for details.
+
 <CustomGBForce>
 ===============
 
@@ -2237,30 +2286,8 @@ attribute name must match the name of the parameter.  For instance, if there is
 a per-atom parameter with the name “radius”, then every :code:`<Atom>` tag
 must include an attribute called :code:`radius`\ .
 
-CustomGBForce also allows you to define tabulated functions.  To define a
-function, include a :code:`<Function>` tag inside the
-:code:`<CustomGBForce>` tag:
-
-.. code-block:: xml
-
-    <Function name="myfn" min="-5" max="5">
-    0.983674857694 -0.980096396266 -0.975743130031 -0.970451936613 -0.964027580076
-    -0.956237458128 -0.946806012846 -0.935409070603 -0.921668554406 -0.905148253645
-    -0.885351648202 -0.861723159313 -0.833654607012 -0.800499021761 -0.761594155956
-    -0.716297870199 -0.664036770268 -0.604367777117 -0.537049566998 -0.46211715726
-    -0.379948962255 -0.291312612452 -0.197375320225 -0.099667994625 0.0
-    0.099667994625 0.197375320225 0.291312612452 0.379948962255 0.46211715726
-    0.537049566998 0.604367777117 0.664036770268 0.716297870199 0.761594155956
-    0.800499021761 0.833654607012 0.861723159313 0.885351648202 0.905148253645
-    0.921668554406 0.935409070603 0.946806012846 0.956237458128 0.964027580076
-    0.970451936613 0.975743130031 0.980096396266 0.983674857694 0.986614298151
-    0.989027402201
-    </Function>
-
-The tag’s attributes define the name of the function and the range of values for
-which it is defined.  The tabulated values are listed inside the body of the
-tag, with successive values separated by white space.  Again, see the API
-documentation for more details.
+CustomGBForce also allows you to define tabulated functions.  See section
+:ref:`tabulated-functions` for details.
 
 Writing Custom Expressions
 ==========================
@@ -2297,6 +2324,55 @@ is exactly equivalent to
 The definition of an intermediate value may itself involve other intermediate
 values.  All uses of a value must appear *before* that value’s definition.
 
+.. _tabulated-functions:
+
+TabulatedFunctions
+==================
+
+Some forces, such as CustomNonbondedForce and CustomGBForce, allow you to define
+tabulated functions.  To define a function, include a :code:`<Function>` tag inside the
+:code:`<CustomNonbondedForce>` or :code:`<CustomGBForce>` tag:
+
+.. code-block:: xml
+
+    <Function name="myfn" type="Continuous1D" min="-5" max="5">
+    0.983674857694 -0.980096396266 -0.975743130031 -0.970451936613 -0.964027580076
+    -0.956237458128 -0.946806012846 -0.935409070603 -0.921668554406 -0.905148253645
+    -0.885351648202 -0.861723159313 -0.833654607012 -0.800499021761 -0.761594155956
+    -0.716297870199 -0.664036770268 -0.604367777117 -0.537049566998 -0.46211715726
+    -0.379948962255 -0.291312612452 -0.197375320225 -0.099667994625 0.0
+    0.099667994625 0.197375320225 0.291312612452 0.379948962255 0.46211715726
+    0.537049566998 0.604367777117 0.664036770268 0.716297870199 0.761594155956
+    0.800499021761 0.833654607012 0.861723159313 0.885351648202 0.905148253645
+    0.921668554406 0.935409070603 0.946806012846 0.956237458128 0.964027580076
+    0.970451936613 0.975743130031 0.980096396266 0.983674857694 0.986614298151
+    0.989027402201
+    </Function>
+
+The tag’s attributes define the name of the function, the type of function, and
+the range of values for which it is defined.  The required set of attributed
+depends on the function type:
+
+.. tabularcolumns:: |l|L|
+
+============  =======================================================
+Type          Required Attributes
+============  =======================================================
+Continuous1D  min, max
+Continuous2D  xmin, ymin, xmax, ymax, xsize, ysize
+Continuous3D  xmin, ymin, zmin, xmax, ymax, zmax, xsize, ysize, zsize
+Discrete1D
+Discrete2D    xsize, ysize
+Discrete3D    xsize, ysize, zsize
+============  =======================================================
+
+
+The "min" and "max" attributes define the range of the independent variables for
+a continuous function.  The "size" attributes define the size of the table along
+each axis.  The tabulated values are listed inside the body of the tag, with
+successive values separated by white space.  See the API documentation for more
+details.
+
 
 Using Multiple Files
 ********************

docs/usersguide/references.bib

@@ -336,10 +336,20 @@
    type = {Journal Article}
 }
 
+@article{Shi2013
+   author = {Shi, Yue and Xia, Zhen and Zhang, Jiajing and Best, Robert and Wu, Chuanjie and Ponder, Jay W. and Ren, Pengyu},
+   title = {Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins},
+   journal = {Journal of Chemical Theory and Computation},
+   volume = {9},
+   number = {9},
+   pages = {4046-4063},
+   year = {2013},
+   type = {Journal Article}
+}
+
 @article{Shirts2007
    author = {Shirts, Michael R. and Mobley, David L. and Chodera, John D. and Pande, Vijay S.},
-   title = {Accurate and Efficient Corrections for Missing Dispersion Interactions in Molecular
-Simulations},
+   title = {Accurate and Efficient Corrections for Missing Dispersion Interactions in Molecular Simulations},
    journal = {Journal of Physical Chemistry B},
    volume = {111},
    pages = {13052-13063},
@@ -399,3 +409,12 @@ Simulations},
    type = {Journal Article}
 }
 
+@article{Wang2014
+   author = {Wang, Lee-Ping and Martinez, Todd J. and Pande, Vijay S.},
+   title = {Building force fields: an automatic, systematic, and reproducible approach},
+   journal = {Journal of Physical Chemistry Letters},
+   volume = {5},
+   pages = {1885-1891},
+   year = {2014},
+   type = {Journal Article}
+}

openmmapi/include/openmm/Context.h

@@ -35,12 +35,12 @@
 #include "Integrator.h"
 #include "State.h"
 #include "System.h"
-#include <ctime>
 #include <iosfwd>
 #include <map>
 #include <string>
 #include <vector>
 #include "internal/windowsExport.h"
+#include "internal/OSRngSeed.h"
 
 namespace OpenMM {
 
@@ -164,7 +164,7 @@ public:
      * @param temperature    the temperature for which to select the velocities (measured in Kelvin)
      * @param randomSeed     the random number seed to use when selecting velocities
      */
-    void setVelocitiesToTemperature(double temperature, int randomSeed=time(NULL));
+    void setVelocitiesToTemperature(double temperature, int randomSeed=osrngseed());
     /**
      * Get the value of an adjustable parameter defined by a Force object in the System.
      * 

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -344,6 +344,8 @@ public:
     /**
      * Add a particle pair to the list of interactions that should be excluded.
      * 
+     * In many cases, you can use createExclusionsFromBonds() rather than adding each exclusion explicitly.
+     *
      * @param particle1  the index of the first particle in the pair
      * @param particle2  the index of the second particle in the pair
      * @return the index of the exclusion that was added
@@ -365,6 +367,15 @@ public:
      * @param particle2  the index of the second particle in the pair
      */
     void setExclusionParticles(int index, int particle1, int particle2);
+    /**
+     * Identify exclusions based on the molecular topology.  Particles which are separated by up to a specified number of
+     * bonds are added as exclusions.
+     *
+     * @param bonds       the set of bonds based on which to construct exclusions.  Each element specifies the indices of
+     *                    two particles that are bonded to each other.
+     * @param bondCutoff  pairs of particles that are separated by this many bonds or fewer are added to the list of exclusions
+     */
+    void createExclusionsFromBonds(const std::vector<std::pair<int, int> >& bonds, int bondCutoff);
     /**
      * Add a tabulated function that may appear in the energy expression.
      *

openmmapi/include/openmm/internal/OSRngSeed.h

+#ifndef OPENMM_OSRNGSEED_H_
+#define OPENMM_OSRNGSEED_H_
+
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2013 Stanford University and the Authors.           *
+ * Authors: Robert T. McGibbon                                                *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+
+
+/**
+ * Return an integer int for use as a seed for a random number generator.
+ * 
+ * The behavior of this function is platform dependent. On Windows machines,
+ * this uses CryptGenRandom from the CryptoAPI to get a single int. On other
+ * platforms (*nix, apple), we read from /dev/urandom
+ */
+int osrngseed(void);
+
+#endif /* OPENMM_OSRNGSEED_H_ */
\ No newline at end of file

openmmapi/src/AndersenThermostat.cpp

@@ -31,13 +31,13 @@
 
 #include "openmm/AndersenThermostat.h"
 #include "openmm/internal/AndersenThermostatImpl.h"
-#include <ctime>
+#include "openmm/internal/OSRngSeed.h"
 
 using namespace OpenMM;
 
 AndersenThermostat::AndersenThermostat(double defaultTemperature, double defaultCollisionFrequency) :
         defaultTemp(defaultTemperature), defaultFreq(defaultCollisionFrequency) {
-    setRandomNumberSeed((int) time(NULL));
+    setRandomNumberSeed(osrngseed());
 }
 
 ForceImpl* AndersenThermostat::createImpl() const {

openmmapi/src/BrownianIntegrator.cpp

@@ -33,8 +33,8 @@
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/OSRngSeed.h"
 #include "openmm/kernels.h"
-#include <ctime>
 #include <string>
 
 using namespace OpenMM;
@@ -46,7 +46,7 @@ BrownianIntegrator::BrownianIntegrator(double temperature, double frictionCoeff,
     setFriction(frictionCoeff);
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
-    setRandomNumberSeed((int) time(NULL));
+    setRandomNumberSeed(osrngseed());
 }
 
 void BrownianIntegrator::initialize(ContextImpl& contextRef) {

openmmapi/src/CustomIntegrator.cpp

@@ -34,8 +34,8 @@
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/AssertionUtilities.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/OSRngSeed.h"
 #include "openmm/kernels.h"
-#include <ctime>
 #include <string>
 
 using namespace OpenMM;
@@ -45,7 +45,7 @@ using std::vector;
 CustomIntegrator::CustomIntegrator(double stepSize) : globalsAreCurrent(true), forcesAreValid(false) {
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
-    setRandomNumberSeed((int) time(NULL));
+    setRandomNumberSeed(osrngseed());
     kineticEnergy = "m*v*v/2";
 }
 

openmmapi/src/CustomNonbondedForce.cpp

@@ -174,6 +174,33 @@ void CustomNonbondedForce::setExclusionParticles(int index, int particle1, int p
     exclusions[index].particle1 = particle1;
     exclusions[index].particle2 = particle2;
 }
+
+void CustomNonbondedForce::createExclusionsFromBonds(const vector<pair<int, int> >& bonds, int bondCutoff) {
+    if (bondCutoff < 1)
+        return;
+    vector<set<int> > exclusions(particles.size());
+    vector<set<int> > bonded12(exclusions.size());
+    for (int i = 0; i < (int) bonds.size(); ++i) {
+        int p1 = bonds[i].first;
+        int p2 = bonds[i].second;
+        exclusions[p1].insert(p2);
+        exclusions[p2].insert(p1);
+        bonded12[p1].insert(p2);
+        bonded12[p2].insert(p1);
+    }
+    for (int level = 0; level < bondCutoff-1; level++) {
+        vector<set<int> > currentExclusions = exclusions;
+        for (int i = 0; i < (int) particles.size(); i++) {
+            for (set<int>::const_iterator iter = currentExclusions[i].begin(); iter != currentExclusions[i].end(); ++iter)
+                exclusions[*iter].insert(bonded12[i].begin(), bonded12[i].end());
+        }
+    }
+    for (int i = 0; i < (int) exclusions.size(); ++i)
+        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
+            if (*iter < i)
+                addExclusion(*iter, i);
+}
+
 int CustomNonbondedForce::addTabulatedFunction(const std::string& name, TabulatedFunction* function) {
     functions.push_back(FunctionInfo(name, function));
     return functions.size()-1;

openmmapi/src/LangevinIntegrator.cpp

@@ -33,8 +33,8 @@
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/OSRngSeed.h"
 #include "openmm/kernels.h"
-#include <ctime>
 #include <string>
 
 using namespace OpenMM;
@@ -46,7 +46,7 @@ LangevinIntegrator::LangevinIntegrator(double temperature, double frictionCoeff,
     setFriction(frictionCoeff);
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
-    setRandomNumberSeed((int) time(NULL));
+    setRandomNumberSeed(osrngseed());
 }
 
 void LangevinIntegrator::initialize(ContextImpl& contextRef) {

openmmapi/src/MonteCarloAnisotropicBarostat.cpp

@@ -31,13 +31,13 @@
 
 #include "openmm/MonteCarloAnisotropicBarostat.h"
 #include "openmm/internal/MonteCarloAnisotropicBarostatImpl.h"
-#include <ctime>
+#include "openmm/internal/OSRngSeed.h"
 
 using namespace OpenMM;
 
 MonteCarloAnisotropicBarostat::MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double temperature, bool scaleX, bool scaleY, bool scaleZ, int frequency) :
         defaultPressure(defaultPressure), temperature(temperature), scaleX(scaleX), scaleY(scaleY), scaleZ(scaleZ), frequency(frequency) {
-    setRandomNumberSeed((int) time(NULL));
+    setRandomNumberSeed(osrngseed());
 }
 
 ForceImpl* MonteCarloAnisotropicBarostat::createImpl() const {

openmmapi/src/MonteCarloBarostat.cpp

@@ -31,13 +31,13 @@
 
 #include "openmm/MonteCarloBarostat.h"
 #include "openmm/internal/MonteCarloBarostatImpl.h"
-#include <ctime>
+#include "openmm/internal/OSRngSeed.h"
 
 using namespace OpenMM;
 
 MonteCarloBarostat::MonteCarloBarostat(double defaultPressure, double temperature, int frequency) :
         defaultPressure(defaultPressure), temperature(temperature), frequency(frequency) {
-    setRandomNumberSeed((int) time(NULL));
+    setRandomNumberSeed(osrngseed());
 }
 
 ForceImpl* MonteCarloBarostat::createImpl() const {

openmmapi/src/OSRngSeed.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2013 Stanford University and the Authors.           *
+ * Authors: Robert T. McGibbon                                                *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include <stdexcept>
+#if defined(_WIN32) || defined(__CYGWIN__)
+#include <windows.h>
+static HCRYPTPROV hCryptProv = 0;
+#pragma comment(lib, "advapi32.lib")
+#else
+#include <fcntl.h>
+#include <unistd.h>
+#endif
+#include "openmm/internal/OSRngSeed.h"
+
+int osrngseed(void) {
+    int value;
+#if defined(_WIN32) || defined(__CYGWIN__)
+    if (!::CryptAcquireContextW(&hCryptProv, 0, 0, PROV_RSA_FULL, CRYPT_VERIFYCONTEXT | CRYPT_SILENT)) {
+        throw std::runtime_error("Failed to initialize Windows random API (CryptoGen)");
+    }
+    if (!CryptGenRandom(hCryptProv, sizeof(int), (BYTE*) &value)) {
+        ::CryptReleaseContext(hCryptProv, 0);
+        throw std::runtime_error("Failed to get random numbers");
+    }
+    if (!::CryptReleaseContext(hCryptProv, 0)) {
+        throw std::runtime_error("Failed to release Windows random API context");
+    }
+#else
+    int m_fd = open("/dev/urandom", O_RDONLY);
+    if (m_fd == -1) {
+        throw std::runtime_error("Failed to open /dev/urandom");
+    }
+    if (read(m_fd, &value, sizeof(int)) != sizeof(int)) {
+        throw std::runtime_error("Failed to read bytes from /dev/urandom");
+    }
+    close(m_fd);
+#endif
+    return value;
+}
\ No newline at end of file

openmmapi/src/VariableLangevinIntegrator.cpp

@@ -33,10 +33,11 @@
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/OSRngSeed.h"
 #include "openmm/kernels.h"
 #include <limits>
 #include <string>
-#include <ctime>
+
 
 using namespace OpenMM;
 using std::string;
@@ -47,7 +48,7 @@ VariableLangevinIntegrator::VariableLangevinIntegrator(double temperature, doubl
     setFriction(frictionCoeff);
     setErrorTolerance(errorTol);
     setConstraintTolerance(1e-5);
-    setRandomNumberSeed((int) time(NULL));
+    setRandomNumberSeed(osrngseed());
 }
 
 void VariableLangevinIntegrator::initialize(ContextImpl& contextRef) {

platforms/cuda/tests/TestCudaVirtualSites.cpp

@@ -524,6 +524,47 @@ void testReordering() {
     }
 }
 
+/**
+ * Test a System where multiple virtual sites are all calculated from the same particles.
+ */
+void testOverlappingSites() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    system.addForce(nonbonded);
+    nonbonded->addParticle(1.0, 0.0, 0.0);
+    nonbonded->addParticle(-0.5, 0.0, 0.0);
+    nonbonded->addParticle(-0.5, 0.0, 0.0);
+    vector<Vec3> positions;
+    positions.push_back(Vec3(0, 0, 0));
+    positions.push_back(Vec3(10, 0, 0));
+    positions.push_back(Vec3(0, 10, 0));
+    for (int i = 0; i < 20; i++) {
+        system.addParticle(0.0);
+        double u = 0.1*((i+1)%4);
+        double v = 0.05*i;
+        system.setVirtualSite(3+i, new ThreeParticleAverageSite(0, 1, 2, u, v, 1-u-v));
+        nonbonded->addParticle(i%2 == 0 ? -1.0 : 1.0, 0.0, 0.0);
+        positions.push_back(Vec3());
+    }
+    VerletIntegrator i1(0.002);
+    VerletIntegrator i2(0.002);
+    Context c1(system, i1, Platform::getPlatformByName("Reference"));
+    Context c2(system, i2, platform);
+    c1.setPositions(positions);
+    c2.setPositions(positions);
+    c1.applyConstraints(0.0001);
+    c2.applyConstraints(0.0001);
+    State s1 = c1.getState(State::Positions | State::Forces);
+    State s2 = c2.getState(State::Positions | State::Forces);
+    for (int i = 0; i < system.getNumParticles(); i++)
+        ASSERT_EQUAL_VEC(s1.getPositions()[i], s2.getPositions()[i], 1e-5);
+    for (int i = 0; i < 3; i++)
+        ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -535,6 +576,7 @@ int main(int argc, char* argv[]) {
         testLocalCoordinates();
         testConservationLaws();
         testReordering();
+        testOverlappingSites();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/opencl/include/OpenCLIntegrationUtilities.h

@@ -121,7 +121,7 @@ private:
     cl::Kernel ccmaPosForceKernel, ccmaVelForceKernel;
     cl::Kernel ccmaMultiplyKernel;
     cl::Kernel ccmaPosUpdateKernel, ccmaVelUpdateKernel;
-    cl::Kernel vsitePositionKernel, vsiteForceKernel;
+    cl::Kernel vsitePositionKernel, vsiteForceKernel, vsiteAddForcesKernel;
     cl::Kernel randomKernel, timeShiftKernel;
     OpenCLArray* posDelta;
     OpenCLArray* settleAtoms;
@@ -152,7 +152,7 @@ private:
     OpenCLArray* vsiteLocalCoordsParams;
     int randomPos;
     int lastSeed, numVsites;
-    bool hasInitializedPosConstraintKernels, hasInitializedVelConstraintKernels, ccmaUseDirectBuffer;
+    bool hasInitializedPosConstraintKernels, hasInitializedVelConstraintKernels, ccmaUseDirectBuffer, hasOverlappingVsites;
     struct ShakeCluster;
     struct ConstraintOrderer;
 };

platforms/opencl/src/OpenCLIntegrationUtilities.cpp

@@ -101,7 +101,7 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
         ccmaConstraintMatrixValue(NULL), ccmaDelta1(NULL), ccmaDelta2(NULL), ccmaConverged(NULL), ccmaConvergedHostBuffer(NULL),
         vsite2AvgAtoms(NULL), vsite2AvgWeights(NULL), vsite3AvgAtoms(NULL), vsite3AvgWeights(NULL),
         vsiteOutOfPlaneAtoms(NULL), vsiteOutOfPlaneWeights(NULL), vsiteLocalCoordsAtoms(NULL), vsiteLocalCoordsParams(NULL),
-        hasInitializedPosConstraintKernels(false), hasInitializedVelConstraintKernels(false) {
+        hasInitializedPosConstraintKernels(false), hasInitializedVelConstraintKernels(false), hasOverlappingVsites(false) {
     // Create workspace arrays.
 
     if (context.getUseDoublePrecision() || context.getUseMixedPrecision()) {
@@ -649,6 +649,7 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
     int num3Avg = vsite3AvgAtomVec.size();
     int numOutOfPlane = vsiteOutOfPlaneAtomVec.size();
     int numLocalCoords = vsiteLocalCoordsAtomVec.size();
+    numVsites = num2Avg+num3Avg+numOutOfPlane+numLocalCoords;
     vsite2AvgAtoms = OpenCLArray::create<mm_int4>(context, max(1, num2Avg), "vsite2AvgAtoms");
     vsite3AvgAtoms = OpenCLArray::create<mm_int4>(context, max(1, num3Avg), "vsite3AvgAtoms");
     vsiteOutOfPlaneAtoms = OpenCLArray::create<mm_int4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneAtoms");
@@ -706,6 +707,20 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
         }
     }
     
+    // If multiple virtual sites depend on the same particle, make sure the force distribution
+    // can be done safely.
+    
+    vector<int> atomCounts(numAtoms, 0);
+    for (int i = 0; i < numAtoms; i++)
+        if (system.isVirtualSite(i))
+            for (int j = 0; j < system.getVirtualSite(i).getNumParticles(); j++)
+                atomCounts[system.getVirtualSite(i).getParticle(j)]++;
+    for (int i = 0; i < numAtoms; i++)
+        if (atomCounts[i] > 1)
+            hasOverlappingVsites = true;
+    if (hasOverlappingVsites && context.getUseDoublePrecision() && !context.getSupports64BitGlobalAtomics())
+        throw OpenMMException("This device does not support 64 bit atomics.  Cannot use double precision when multiple virtual sites depend on the same atom.");
+    
     // Create the kernels for virtual sites.
 
     map<string, string> defines;
@@ -713,6 +728,10 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
     defines["NUM_3_AVERAGE"] = context.intToString(num3Avg);
     defines["NUM_OUT_OF_PLANE"] = context.intToString(numOutOfPlane);
     defines["NUM_LOCAL_COORDS"] = context.intToString(numLocalCoords);
+    defines["NUM_ATOMS"] = context.intToString(numAtoms);
+    defines["PADDED_NUM_ATOMS"] = context.intToString(context.getPaddedNumAtoms());
+    if (hasOverlappingVsites)
+        defines["HAS_OVERLAPPING_VSITES"] = "1";
     cl::Program vsiteProgram = context.createProgram(OpenCLKernelSources::virtualSites, defines);
     vsitePositionKernel = cl::Kernel(vsiteProgram, "computeVirtualSites");
     int index = 0;
@@ -731,6 +750,8 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
     index = 0;
     vsiteForceKernel.setArg<cl::Buffer>(index++, context.getPosq().getDeviceBuffer());
     index++; // Skip argument 1: the force array hasn't been created yet.
+    if (context.getSupports64BitGlobalAtomics())
+        index++; // Skip argument 2: the force array hasn't been created yet.
     if (context.getUseMixedPrecision())
         vsiteForceKernel.setArg<cl::Buffer>(index++, context.getPosqCorrection().getDeviceBuffer());
     vsiteForceKernel.setArg<cl::Buffer>(index++, vsite2AvgAtoms->getDeviceBuffer());
@@ -741,7 +762,8 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
     vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteOutOfPlaneWeights->getDeviceBuffer());
     vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsAtoms->getDeviceBuffer());
     vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsParams->getDeviceBuffer());
-    numVsites = num2Avg+num3Avg+numOutOfPlane+numLocalCoords;
+    if (hasOverlappingVsites && context.getSupports64BitGlobalAtomics())
+        vsiteAddForcesKernel = cl::Kernel(vsiteProgram, "addDistributedForces");
 }
 
 OpenCLIntegrationUtilities::~OpenCLIntegrationUtilities() {
@@ -941,8 +963,25 @@ void OpenCLIntegrationUtilities::computeVirtualSites() {
 
 void OpenCLIntegrationUtilities::distributeForcesFromVirtualSites() {
     if (numVsites > 0) {
+        // Set arguments that didn't exist yet in the constructor.
+        
         vsiteForceKernel.setArg<cl::Buffer>(1, context.getForce().getDeviceBuffer());
+        if (context.getSupports64BitGlobalAtomics()) {
+            vsiteForceKernel.setArg<cl::Buffer>(2, context.getLongForceBuffer().getDeviceBuffer());
+            if (hasOverlappingVsites) {
+                // We'll be using 64 bit atomics for the force redistribution, so clear the buffer.
+                
+                context.clearBuffer(context.getLongForceBuffer());
+            }
+        }
         context.executeKernel(vsiteForceKernel, numVsites);
+        if (context.getSupports64BitGlobalAtomics() && hasOverlappingVsites) {
+            // Add the redistributed forces from the virtual sites to the main force array.
+            
+            vsiteAddForcesKernel.setArg<cl::Buffer>(0, context.getLongForceBuffer().getDeviceBuffer());
+            vsiteAddForcesKernel.setArg<cl::Buffer>(1, context.getForce().getDeviceBuffer());
+            context.executeKernel(vsiteAddForcesKernel, context.getNumAtoms());
+        }
     }
 }
 

platforms/opencl/src/kernels/virtualSites.cl

@@ -108,10 +108,66 @@ __kernel void computeVirtualSites(__global real4* restrict posq,
     }
 }
 
+#ifdef HAS_OVERLAPPING_VSITES
+    #ifdef SUPPORTS_64_BIT_ATOMICS
+        // We will use 64 bit atomics for force redistribution.
+
+        #define ADD_FORCE(index, f) addForce(index, f, longForce);
+
+        #pragma OPENCL EXTENSION cl_khr_int64_base_atomics : enable
+
+        void addForce(int index, real4 f,  __global long* longForce) {
+            atom_add(&longForce[index], (long) (f.x*0x100000000));
+            atom_add(&longForce[index+PADDED_NUM_ATOMS], (long) (f.y*0x100000000));
+            atom_add(&longForce[index+2*PADDED_NUM_ATOMS], (long) (f.z*0x100000000));
+        }
+
+        __kernel void addDistributedForces(__global const long* restrict longForces, __global real4* restrict forces) {
+            real scale = 1/(real) 0x100000000;
+            for (int index = get_global_id(0); index < NUM_ATOMS; index += get_global_size(0)) {
+                real4 f = (real4) (scale*longForces[index], scale*longForces[index+PADDED_NUM_ATOMS], scale*longForces[index+2*PADDED_NUM_ATOMS], 0);
+                forces[index] += f;
+            }
+        }
+    #else
+        // 64 bit atomics aren't supported, so we have to use atomic_cmpxchg() for force redistribution.
+        
+        #define ADD_FORCE(index, f) addForce(index, f, force);
+
+        void atomicAddFloat(__global float* p, float v) {
+            __global int* ip = (__global int*) p;
+            while (true) {
+                union {
+                    float f;
+                    int i;
+                } oldval, newval;
+                oldval.f = *p;
+                newval.f = oldval.f+v;
+                if (atomic_cmpxchg(ip, oldval.i, newval.i) == oldval.i)
+                    return;
+            }
+        }
+
+        void addForce(int index, float4 f, __global float4* force) {
+            __global float* components = (__global float*) force;
+            atomicAddFloat(&components[4*index], f.x);
+            atomicAddFloat(&components[4*index+1], f.y);
+            atomicAddFloat(&components[4*index+2], f.z);
+        }
+    #endif
+#else
+    // There are no overlapping virtual sites, so we can just store forces directly.
+
+    #define ADD_FORCE(index, f) force[index].xyz += (f).xyz;
+#endif
+
 /**
  * Distribute forces from virtual sites to the atoms they are based on.
  */
 __kernel void distributeForces(__global const real4* restrict posq, __global real4* restrict force,
+#ifdef SUPPORTS_64_BIT_ATOMICS
+        __global long* restrict longForce,
+#endif
 #ifdef USE_MIXED_PRECISION
         __global real4* restrict posqCorrection,
 #endif
@@ -129,12 +185,8 @@ __kernel void distributeForces(__global const real4* restrict posq, __global rea
         int4 atoms = avg2Atoms[index];
         real2 weights = avg2Weights[index];
         real4 f = force[atoms.x];
-        real4 f1 = force[atoms.y];
-        real4 f2 = force[atoms.z];
-        f1.xyz += f.xyz*weights.x;
-        f2.xyz += f.xyz*weights.y;
-        force[atoms.y] = f1;
-        force[atoms.z] = f2;
+        ADD_FORCE(atoms.y, f*weights.x);
+        ADD_FORCE(atoms.z, f*weights.y);
     }
     
     // Three particle average sites.
@@ -143,15 +195,9 @@ __kernel void distributeForces(__global const real4* restrict posq, __global rea
         int4 atoms = avg3Atoms[index];
         real4 weights = avg3Weights[index];
         real4 f = force[atoms.x];
-        real4 f1 = force[atoms.y];
-        real4 f2 = force[atoms.z];
-        real4 f3 = force[atoms.w];
-        f1.xyz += f.xyz*weights.x;
-        f2.xyz += f.xyz*weights.y;
-        f3.xyz += f.xyz*weights.z;
-        force[atoms.y] = f1;
-        force[atoms.z] = f2;
-        force[atoms.w] = f3;
+        ADD_FORCE(atoms.y, f*weights.x);
+        ADD_FORCE(atoms.z, f*weights.y);
+        ADD_FORCE(atoms.w, f*weights.z);
     }
     
     // Out of plane sites.
@@ -165,21 +211,15 @@ __kernel void distributeForces(__global const real4* restrict posq, __global rea
         mixed4 v12 = pos2-pos1;
         mixed4 v13 = pos3-pos1;
         real4 f = force[atoms.x];
-        real4 f1 = force[atoms.y];
-        real4 f2 = force[atoms.z];
-        real4 f3 = force[atoms.w];
         real4 fp2 = (real4) (weights.x*f.x - weights.z*v13.z*f.y + weights.z*v13.y*f.z,
                    weights.z*v13.z*f.x + weights.x*f.y - weights.z*v13.x*f.z,
                   -weights.z*v13.y*f.x + weights.z*v13.x*f.y + weights.x*f.z, 0.0f);
         real4 fp3 = (real4) (weights.y*f.x + weights.z*v12.z*f.y - weights.z*v12.y*f.z,
                   -weights.z*v12.z*f.x + weights.y*f.y + weights.z*v12.x*f.z,
                    weights.z*v12.y*f.x - weights.z*v12.x*f.y + weights.y*f.z, 0.0f);
-        f1.xyz += f.xyz-fp2.xyz-fp3.xyz;
-        f2.xyz += fp2.xyz;
-        f3.xyz += fp3.xyz;
-        force[atoms.y] = f1;
-        force[atoms.z] = f2;
-        force[atoms.w] = f3;
+        ADD_FORCE(atoms.y, f-fp2-fp3);
+        ADD_FORCE(atoms.z, fp2);
+        ADD_FORCE(atoms.w, fp3);
     }
     
     // Local coordinates sites.
@@ -230,9 +270,9 @@ __kernel void distributeForces(__global const real4* restrict posq, __global rea
         mixed sz3 = t32*dz.x-t31*dz.y;
         mixed4 wxScaled = wx*invNormXdir;
         real4 f = force[atoms.x];
-        real4 f1 = force[atoms.y];
-        real4 f2 = force[atoms.z];
-        real4 f3 = force[atoms.w];
+        real4 f1 = 0;
+        real4 f2 = 0;
+        real4 f3 = 0;
         mixed4 fp1 = localPosition*f.x;
         mixed4 fp2 = localPosition*f.y;
         mixed4 fp3 = localPosition*f.z;
@@ -263,8 +303,8 @@ __kernel void distributeForces(__global const real4* restrict posq, __global rea
         f3.x += fp3.x*wxScaled.z*( -dx.z*dx.x) + fp3.z*(dz.z*sx3+t32) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.z + dy.z*sx3 + dx.x*t33);
         f3.y += fp3.x*wxScaled.z*( -dx.z*dx.y) + fp3.z*(dz.z*sy3-t31) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.z + dy.z*sy3 + dx.y*t33);
         f3.z += fp3.x*wxScaled.z*(1-dx.z*dx.z) + fp3.z*(dz.z*sz3    ) + fp3.y*((-dx.z*dy.z     )*wxScaled.z + dy.z*sz3 - dx.x*t31 - dx.y*t32) + f.z*originWeights.z;
-        force[atoms.y] = f1;
-        force[atoms.z] = f2;
-        force[atoms.w] = f3;
+        ADD_FORCE(atoms.y, f1);
+        ADD_FORCE(atoms.z, f2);
+        ADD_FORCE(atoms.w, f3);
     }
 }