Eliminated Stream and related classes, replacing their functionality with UpdateStateDataKernel.

openmmapi/include/openmm/internal/RBTorsionForceImpl.h

@@ -56,7 +56,7 @@ public:
     void updateContextState(ContextImpl& context) {
         // This force field doesn't update the state directly.
     }
-    void calcForces(ContextImpl& context, Stream& forces);
+    void calcForces(ContextImpl& context);
     double calcEnergy(ContextImpl& context);
     std::map<std::string, double> getDefaultParameters() {
         return std::map<std::string, double>(); // This force field doesn't define any parameters.

openmmapi/src/Context.cpp

@@ -79,16 +79,16 @@ State Context::getState(int types) const {
         state.setEnergy(impl->calcKineticEnergy(), impl->calcPotentialEnergy());
     if (types&State::Forces) {
         impl->calcForces();
-        impl->getForces().saveToArray(&state.updForces()[0]);
+        impl->getForces(state.updForces());
     }
     if (types&State::Parameters) {
         for (map<string, double>::const_iterator iter = impl->parameters.begin(); iter != impl->parameters.end(); iter++)
             state.updParameters()[iter->first] = iter->second;
     }
     if (types&State::Positions)
-        impl->getPositions().saveToArray(&state.updPositions()[0]);
+        impl->getPositions(state.updPositions());
     if (types&State::Velocities)
-        impl->getVelocities().saveToArray(&state.updVelocities()[0]);
+        impl->getVelocities(state.updVelocities());
     return state;
 }
 
@@ -99,13 +99,13 @@ void Context::setTime(double time) {
 void Context::setPositions(const vector<Vec3>& positions) {
     if ((int) positions.size() != impl->getSystem().getNumParticles())
         throw OpenMMException("Called setPositions() on a Context with the wrong number of positions");
-    impl->getPositions().loadFromArray(&positions[0]);
+    impl->setPositions(positions);
 }
 
 void Context::setVelocities(const vector<Vec3>& velocities) {
     if ((int) velocities.size() != impl->getSystem().getNumParticles())
         throw OpenMMException("Called setVelocities() on a Context with the wrong number of velocities");
-    impl->getVelocities().loadFromArray(&velocities[0]);
+    impl->setVelocities(velocities);
 }
 
 double Context::getParameter(const string& name) {

openmmapi/src/ContextImpl.cpp

@@ -45,14 +45,11 @@ using std::vector;
 using std::string;
 
 ContextImpl::ContextImpl(Context& owner, System& system, Integrator& integrator, Platform* platform) :
-				owner(owner), system(system), 
-				integrator(integrator), platform(platform), 
-				platformData(NULL)
-{
+         owner(owner), system(system), integrator(integrator), platform(platform), platformData(NULL) {
     vector<string> kernelNames;
     kernelNames.push_back(CalcKineticEnergyKernel::Name());
     kernelNames.push_back(CalcForcesAndEnergyKernel::Name());
-    kernelNames.push_back(UpdateTimeKernel::Name());
+    kernelNames.push_back(UpdateStateDataKernel::Name());
     for (int i = 0; i < system.getNumForces(); ++i) {
         forceImpls.push_back(system.getForce(i).createImpl());
         map<string, double> forceParameters = forceImpls[forceImpls.size()-1]->getDefaultParameters();
@@ -71,16 +68,12 @@ ContextImpl::ContextImpl(Context& owner, System& system, Integrator& integrator,
     dynamic_cast<CalcForcesAndEnergyKernel&>(initializeForcesKernel.getImpl()).initialize(system);
     kineticEnergyKernel = platform->createKernel(CalcKineticEnergyKernel::Name(), *this);
     dynamic_cast<CalcKineticEnergyKernel&>(kineticEnergyKernel.getImpl()).initialize(system);
-    updateTimeKernel = platform->createKernel(UpdateTimeKernel::Name(), *this);
-    dynamic_cast<UpdateTimeKernel&>(updateTimeKernel.getImpl()).initialize(system);
+    updateStateDataKernel = platform->createKernel(UpdateStateDataKernel::Name(), *this);
+    dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).initialize(system);
     for (size_t i = 0; i < forceImpls.size(); ++i)
         forceImpls[i]->initialize(*this);
     integrator.initialize(*this);
-    positions = platform->createStream("particlePositions", system.getNumParticles(), Stream::Double3, *this);
-    velocities = platform->createStream("particleVelocities", system.getNumParticles(), Stream::Double3, *this);
-    forces = platform->createStream("particleForces", system.getNumParticles(), Stream::Double3, *this);
-    double zero[] = {0.0, 0.0, 0.0};
-    velocities.fillWithValue(&zero);
+    dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).setVelocities(*this, vector<Vec3>(system.getNumParticles()));
 }
 
 ContextImpl::~ContextImpl() {
@@ -90,11 +83,31 @@ ContextImpl::~ContextImpl() {
 }
 
 double ContextImpl::getTime() const {
-    return dynamic_cast<const UpdateTimeKernel&>(updateTimeKernel.getImpl()).getTime(*this);
+    return dynamic_cast<const UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).getTime(*this);
 }
 
 void ContextImpl::setTime(double t) {
-    dynamic_cast<UpdateTimeKernel&>(updateTimeKernel.getImpl()).setTime(*this, t);
+    dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).setTime(*this, t);
+}
+
+void ContextImpl::getPositions(std::vector<Vec3>& positions) {
+    dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).getPositions(*this, positions);
+}
+
+void ContextImpl::setPositions(const std::vector<Vec3>& positions) {
+    dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).setPositions(*this, positions);
+}
+
+void ContextImpl::getVelocities(std::vector<Vec3>& velocities) {
+    dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).getVelocities(*this, velocities);
+}
+
+void ContextImpl::setVelocities(const std::vector<Vec3>& velocities) {
+    dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).setVelocities(*this, velocities);
+}
+
+void ContextImpl::getForces(std::vector<Vec3>& forces) {
+    dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).getForces(*this, forces);
 }
 
 double ContextImpl::getParameter(std::string name) {
@@ -113,7 +126,7 @@ void ContextImpl::calcForces() {
     CalcForcesAndEnergyKernel& kernel = dynamic_cast<CalcForcesAndEnergyKernel&>(initializeForcesKernel.getImpl());
     kernel.beginForceComputation(*this);
     for (int i = 0; i < (int) forceImpls.size(); ++i)
-        forceImpls[i]->calcForces(*this, forces);
+        forceImpls[i]->calcForces(*this);
     kernel.finishForceComputation(*this);
 }
 

openmmapi/src/CustomNonbondedForceImpl.cpp

@@ -104,7 +104,7 @@ void CustomNonbondedForceImpl::initialize(ContextImpl& context) {
     dynamic_cast<CalcCustomNonbondedForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }
 
-void CustomNonbondedForceImpl::calcForces(ContextImpl& context, Stream& forces) {
+void CustomNonbondedForceImpl::calcForces(ContextImpl& context) {
     dynamic_cast<CalcCustomNonbondedForceKernel&>(kernel.getImpl()).executeForces(context);
 }
 

openmmapi/src/GBSAOBCForceImpl.cpp

@@ -48,7 +48,7 @@ void GBSAOBCForceImpl::initialize(ContextImpl& context) {
     dynamic_cast<CalcGBSAOBCForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }
 
-void GBSAOBCForceImpl::calcForces(ContextImpl& context, Stream& forces) {
+void GBSAOBCForceImpl::calcForces(ContextImpl& context) {
     dynamic_cast<CalcGBSAOBCForceKernel&>(kernel.getImpl()).executeForces(context);
 }
 

openmmapi/src/GBVIForceImpl.cpp

@@ -188,7 +188,7 @@ void GBVIForceImpl::findScaledRadii( int numberOfParticles, const std::vector<st
 
 }
 
-void GBVIForceImpl::calcForces(ContextImpl& context, Stream& forces) {
+void GBVIForceImpl::calcForces(ContextImpl& context) {
     dynamic_cast<CalcGBVIForceKernel&>(kernel.getImpl()).executeForces(context);
 }
 

openmmapi/src/HarmonicAngleForceImpl.cpp

@@ -49,7 +49,7 @@ void HarmonicAngleForceImpl::initialize(ContextImpl& context) {
     dynamic_cast<CalcHarmonicAngleForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }
 
-void HarmonicAngleForceImpl::calcForces(ContextImpl& context, Stream& forces) {
+void HarmonicAngleForceImpl::calcForces(ContextImpl& context) {
     dynamic_cast<CalcHarmonicAngleForceKernel&>(kernel.getImpl()).executeForces(context);
 }
 

openmmapi/src/HarmonicBondForceImpl.cpp

@@ -49,7 +49,7 @@ void HarmonicBondForceImpl::initialize(ContextImpl& context) {
     dynamic_cast<CalcHarmonicBondForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }
 
-void HarmonicBondForceImpl::calcForces(ContextImpl& context, Stream& forces) {
+void HarmonicBondForceImpl::calcForces(ContextImpl& context) {
     dynamic_cast<CalcHarmonicBondForceKernel&>(kernel.getImpl()).executeForces(context);
 }
 

openmmapi/src/NonbondedForceImpl.cpp

@@ -86,7 +86,7 @@ void NonbondedForceImpl::initialize(ContextImpl& context) {
     dynamic_cast<CalcNonbondedForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }
 
-void NonbondedForceImpl::calcForces(ContextImpl& context, Stream& forces) {
+void NonbondedForceImpl::calcForces(ContextImpl& context) {
     dynamic_cast<CalcNonbondedForceKernel&>(kernel.getImpl()).executeForces(context);
 }
 

openmmapi/src/PeriodicTorsionForceImpl.cpp

@@ -49,7 +49,7 @@ void PeriodicTorsionForceImpl::initialize(ContextImpl& context) {
     dynamic_cast<CalcPeriodicTorsionForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }
 
-void PeriodicTorsionForceImpl::calcForces(ContextImpl& context, Stream& forces) {
+void PeriodicTorsionForceImpl::calcForces(ContextImpl& context) {
     dynamic_cast<CalcPeriodicTorsionForceKernel&>(kernel.getImpl()).executeForces(context);
 }
 

openmmapi/src/RBTorsionForceImpl.cpp

@@ -49,7 +49,7 @@ void RBTorsionForceImpl::initialize(ContextImpl& context) {
     dynamic_cast<CalcRBTorsionForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }
 
-void RBTorsionForceImpl::calcForces(ContextImpl& context, Stream& forces) {
+void RBTorsionForceImpl::calcForces(ContextImpl& context) {
     dynamic_cast<CalcRBTorsionForceKernel&>(kernel.getImpl()).executeForces(context);
 }
 

platforms/brook/src/BrookKernelFactory.cpp

@@ -60,7 +60,7 @@ KernelImpl* BrookKernelFactory::createKernelImpl( std::string name, const Platfo
 
    // update time
 
-	} else if( name == UpdateTimeKernel::Name() ){
+	} else if( name == UpdateStateDataKernel::Name() ){
 
       return new BrookUpdateTimeKernel( name, platform, openMMBrookInterface );
 

platforms/brook/src/BrookPlatform.cpp

@@ -288,7 +288,7 @@ void BrookPlatform::_initializeKernelFactory( void ){
    registerKernelFactory( CalcGBSAOBCForceKernel::Name(),         factory );
    registerKernelFactory( IntegrateVerletStepKernel::Name(),      factory );
    registerKernelFactory( IntegrateLangevinStepKernel::Name(),    factory );
-   registerKernelFactory( UpdateTimeKernel::Name(),               factory );
+   registerKernelFactory( UpdateStateDataKernel::Name(),               factory );
    // registerKernelFactory( IntegrateBrownianStepKernel::Name(),    factory );
    //registerKernelFactory( ApplyAndersenThermostatKernel::Name(),  factory );
    registerKernelFactory( CalcKineticEnergyKernel::Name(),        factory );

platforms/brook/src/BrookUpdateTimeKernel.cpp

@@ -42,7 +42,7 @@ using namespace std;
  */
 
 BrookUpdateTimeKernel::BrookUpdateTimeKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface ) :
-                     UpdateTimeKernel( name, platform ), _openMMBrookInterface(openMMBrookInterface){
+                     UpdateStateDataKernel( name, platform ), _openMMBrookInterface(openMMBrookInterface){
 
 // ---------------------------------------------------------------------------------------
 

platforms/brook/src/BrookUpdateTimeKernel.h

@@ -35,7 +35,7 @@ namespace OpenMM {
 /**
  * This kernel initializes the forces
  */
-class BrookUpdateTimeKernel : public UpdateTimeKernel {
+class BrookUpdateTimeKernel : public UpdateStateDataKernel {
 
    public:
 

platforms/cuda/include/CudaPlatform.h

@@ -28,7 +28,6 @@
  * -------------------------------------------------------------------------- */
 
 #include "openmm/Platform.h"
-#include "CudaStreamFactory.h"
 
 struct _gpuContext;
 
@@ -52,7 +51,6 @@ public:
     bool supportsDoublePrecision() const;
     const std::string& getPropertyValue(const Context& context, const std::string& property) const;
     void setPropertyValue(Context& context, const std::string& property, const std::string& value) const;
-    const StreamFactory& getDefaultStreamFactory() const;
     void contextCreated(ContextImpl& context) const;
     void contextDestroyed(ContextImpl& context) const;
     /**
@@ -69,8 +67,6 @@ public:
         static const std::string key = "CudaUseBlockingSync";
         return key;
     }
-private:
-    CudaStreamFactory defaultStreamFactory;
 };
 
 class CudaPlatform::PlatformData {

platforms/cuda/include/CudaStreamFactory.h

-#ifndef OPENMM_CUDASTREAMFACTORY_H_
-#define OPENMM_CUDASTREAMFACTORY_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) 2008 Stanford University and the Authors.           *
- * Authors: Peter Eastman                                                     *
- * Contributors:                                                              *
- *                                                                            *
- * 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 3 of the License, or       *
- * (at your option) any later version.                                        *
- *                                                                            *
- * This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.      *
- * -------------------------------------------------------------------------- */
-
-#include "openmm/StreamFactory.h"
-
-namespace OpenMM {
-
-/**
- * This StreamFactory creates all streams for CudaPlatform.
- */
-
-class OPENMM_EXPORT CudaStreamFactory : public StreamFactory {
-public:
-    StreamImpl* createStreamImpl(std::string name, int size, Stream::DataType type, const Platform& platform, ContextImpl& context) const;
-};
-
-} // namespace OpenMM
-
-#endif /*OPENMM_CUDASTREAMFACTORY_H_*/

platforms/cuda/src/CudaKernelFactory.cpp

@@ -35,8 +35,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
     CudaPlatform::PlatformData& data = *static_cast<CudaPlatform::PlatformData*>(context.getPlatformData());
     if (name == CalcForcesAndEnergyKernel::Name())
         return new CudaCalcForcesAndEnergyKernel(name, platform, data);
-    if (name == UpdateTimeKernel::Name())
-        return new CudaUpdateTimeKernel(name, platform, data);
+    if (name == UpdateStateDataKernel::Name())
+        return new CudaUpdateStateDataKernel(name, platform, data);
     if (name == CalcHarmonicBondForceKernel::Name())
         return new CudaCalcHarmonicBondForceKernel(name, platform, data, context.getSystem());
     if (name == CalcHarmonicAngleForceKernel::Name())
@@ -64,7 +64,7 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
     if (name == ApplyAndersenThermostatKernel::Name())
         return new CudaApplyAndersenThermostatKernel(name, platform, data);
     if (name == CalcKineticEnergyKernel::Name())
-        return new CudaCalcKineticEnergyKernel(name, platform);
+        return new CudaCalcKineticEnergyKernel(name, platform, data);
     if (name == RemoveCMMotionKernel::Name())
         return new CudaRemoveCMMotionKernel(name, platform, data);
     throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());

platforms/cuda/src/CudaKernels.cpp

@@ -25,7 +25,6 @@
  * -------------------------------------------------------------------------- */
 
 #include "CudaKernels.h"
-#include "CudaStreamImpl.h"
 #include "openmm/LangevinIntegrator.h"
 #include "openmm/Context.h"
 #include "openmm/internal/ContextImpl.h"
@@ -91,17 +90,84 @@ double CudaCalcForcesAndEnergyKernel::finishEnergyComputation(ContextImpl& conte
     return kReduceEnergy(gpu)+data.ewaldSelfEnergy;
 }
 
-void CudaUpdateTimeKernel::initialize(const System& system) {
+void CudaUpdateStateDataKernel::initialize(const System& system) {
 }
 
-double CudaUpdateTimeKernel::getTime(const ContextImpl& context) const {
+double CudaUpdateStateDataKernel::getTime(const ContextImpl& context) const {
     return data.time;
 }
 
-void CudaUpdateTimeKernel::setTime(ContextImpl& context, double time) {
+void CudaUpdateStateDataKernel::setTime(ContextImpl& context, double time) {
     data.time = time;
 }
 
+void CudaUpdateStateDataKernel::getPositions(ContextImpl& context, std::vector<Vec3>& positions) {
+    _gpuContext* gpu = data.gpu;
+    gpu->psPosq4->Download();
+    int* order = gpu->psAtomIndex->_pSysData;
+    int numParticles = context.getSystem().getNumParticles();
+    positions.resize(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        float4 pos = (*gpu->psPosq4)[i];
+        int3 offset = gpu->posCellOffsets[i];
+        positions[order[i]] = Vec3(pos.x-offset.x*gpu->sim.periodicBoxSizeX, pos.y-offset.y*gpu->sim.periodicBoxSizeY, pos.z-offset.z*gpu->sim.periodicBoxSizeZ);
+    }
+}
+
+void CudaUpdateStateDataKernel::setPositions(ContextImpl& context, const std::vector<Vec3>& positions) {
+    _gpuContext* gpu = data.gpu;
+    int* order = gpu->psAtomIndex->_pSysData;
+    int numParticles = context.getSystem().getNumParticles();
+    for (int i = 0; i < numParticles; ++i) {
+        float4& pos = (*gpu->psPosq4)[i];
+        const Vec3& p = positions[order[i]];
+        pos.x = p[0];
+        pos.y = p[1];
+        pos.z = p[2];
+    }
+    gpu->psPosq4->Upload();
+    for (int i = 0; i < gpu->posCellOffsets.size(); i++)
+        gpu->posCellOffsets[i] = make_int3(0, 0, 0);
+}
+
+void CudaUpdateStateDataKernel::getVelocities(ContextImpl& context, std::vector<Vec3>& velocities) {
+    _gpuContext* gpu = data.gpu;
+    gpu->psVelm4->Download();
+    int* order = gpu->psAtomIndex->_pSysData;
+    int numParticles = context.getSystem().getNumParticles();
+    velocities.resize(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        float4 vel = (*gpu->psVelm4)[i];
+        velocities[order[i]] = Vec3(vel.x, vel.y, vel.z);
+    }
+}
+
+void CudaUpdateStateDataKernel::setVelocities(ContextImpl& context, const std::vector<Vec3>& velocities) {
+    _gpuContext* gpu = data.gpu;
+    int* order = gpu->psAtomIndex->_pSysData;
+    int numParticles = context.getSystem().getNumParticles();
+    for (int i = 0; i < numParticles; ++i) {
+        float4& vel = (*gpu->psVelm4)[i];
+        const Vec3& v = velocities[order[i]];
+        vel.x = v[0];
+        vel.y = v[1];
+        vel.z = v[2];
+    }
+    gpu->psVelm4->Upload();
+}
+
+void CudaUpdateStateDataKernel::getForces(ContextImpl& context, std::vector<Vec3>& forces) {
+    _gpuContext* gpu = data.gpu;
+    int* order = gpu->psAtomIndex->_pSysData;
+    gpu->psForce4->Download();
+    int numParticles = context.getSystem().getNumParticles();
+    forces.resize(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        float4 force = (*gpu->psForce4)[i];
+        forces[order[i]] = Vec3(force.x, force.y, force.z);
+    }
+}
+
 CudaCalcHarmonicBondForceKernel::~CudaCalcHarmonicBondForceKernel() {
 }
 
@@ -772,13 +838,13 @@ double CudaCalcKineticEnergyKernel::execute(ContextImpl& context) {
     // We don't currently have a GPU kernel to do this, so we retrieve the velocities and calculate the energy
     // on the CPU.
     
-    const Stream& velocities = context.getVelocities();
-    double* v = new double[velocities.getSize()*3];
-    velocities.saveToArray(v);
+    _gpuContext* gpu = data.gpu;
+    gpu->psVelm4->Download();
     double energy = 0.0;
-    for (size_t i = 0; i < masses.size(); ++i)
-        energy += masses[i]*(v[i*3]*v[i*3]+v[i*3+1]*v[i*3+1]+v[i*3+2]*v[i*3+2]);
-    delete v;
+    for (int i = 0; i < (int) masses.size(); ++i) {
+        float4 v = (*gpu->psVelm4)[i];
+        energy += masses[i]*(v.x*v.x+v.y*v.y+v.z*v.z);
+    }
     return 0.5*energy;
 }
 

platforms/cuda/src/CudaKernels.h

@@ -91,11 +91,12 @@ private:
 };
 
 /**
- * This kernel is invoked to get or set the current time.
+ * This kernel provides methods for setting and retrieving various state data: time, positions,
+ * velocities, and forces.
  */
-class CudaUpdateTimeKernel : public UpdateTimeKernel {
+class CudaUpdateStateDataKernel : public UpdateStateDataKernel {
 public:
-    CudaUpdateTimeKernel(std::string name, const Platform& platform, CudaPlatform::PlatformData& data) : UpdateTimeKernel(name, platform), data(data) {
+    CudaUpdateStateDataKernel(std::string name, const Platform& platform, CudaPlatform::PlatformData& data) : UpdateStateDataKernel(name, platform), data(data) {
     }
     /**
      * Initialize the kernel.
@@ -115,6 +116,36 @@ public:
      * @param context    the context in which to execute this kernel
      */
     void setTime(ContextImpl& context, double time);
+    /**
+     * Get the positions of all particles.
+     *
+     * @param positions  on exit, this contains the particle positions
+     */
+    void getPositions(ContextImpl& context, std::vector<Vec3>& positions);
+    /**
+     * Set the positions of all particles.
+     *
+     * @param positions  a vector containg the particle positions
+     */
+    void setPositions(ContextImpl& context, const std::vector<Vec3>& positions);
+    /**
+     * Get the velocities of all particles.
+     *
+     * @param velocities  on exit, this contains the particle velocities
+     */
+    void getVelocities(ContextImpl& context, std::vector<Vec3>& velocities);
+    /**
+     * Set the velocities of all particles.
+     *
+     * @param velocities  a vector containg the particle velocities
+     */
+    void setVelocities(ContextImpl& context, const std::vector<Vec3>& velocities);
+    /**
+     * Get the current forces on all particles.
+     *
+     * @param forces  on exit, this contains the forces
+     */
+    void getForces(ContextImpl& context, std::vector<Vec3>& forces);
 private:
     CudaPlatform::PlatformData& data;
 };
@@ -526,7 +557,7 @@ private:
  */
 class CudaCalcKineticEnergyKernel : public CalcKineticEnergyKernel {
 public:
-    CudaCalcKineticEnergyKernel(std::string name, const Platform& platform) : CalcKineticEnergyKernel(name, platform) {
+    CudaCalcKineticEnergyKernel(std::string name, const Platform& platform, CudaPlatform::PlatformData& data) : CalcKineticEnergyKernel(name, platform), data(data) {
     }
     /**
      * Initialize the kernel.
@@ -541,6 +572,7 @@ public:
      */
     double execute(ContextImpl& context);
 private:
+    CudaPlatform::PlatformData& data;
     std::vector<double> masses;
 };