Early beginnings of creating an optimized CPU platform

CMakeLists.txt

@@ -374,6 +374,13 @@ IF(OPENMM_BUILD_OPENCL_LIB)
     ADD_SUBDIRECTORY(platforms/opencl)
 ENDIF(OPENMM_BUILD_OPENCL_LIB)
 
+# Optimized CPU platform
+
+SET(OPENMM_BUILD_CPU_LIB ON CACHE BOOL "Build optimized CPU platform")
+IF(OPENMM_BUILD_CPU_LIB)
+   ADD_SUBDIRECTORY(platforms/cpu)
+ENDIF(OPENMM_BUILD_CPU_LIB)
+
 # Amoeba plugin
 
 SET(OPENMM_BUILD_AMOEBA_PLUGIN ON CACHE BOOL "Build Amoeba plugin")

platforms/cpu/CMakeLists.txt

+#---------------------------------------------------
+# OpenMM CPU Platform
+#
+# Creates OpenMM library, base name=OpenMMCPU.
+# Default libraries are shared & optimized. Variants
+# are created for static (_static) and debug (_d).
+#
+# Windows:
+#   OpenMMCPU[_d].dll
+#   OpenMMCPU[_d].lib
+#   OpenMMCPU_static[_d].lib
+# Unix:
+#   libOpenMMCPU[_d].so
+#   libOpenMMCPU_static[_d].a
+#----------------------------------------------------
+
+IF (APPLE)
+    SET (CMAKE_OSX_DEPLOYMENT_TARGET "10.6")
+ENDIF (APPLE)
+
+SUBDIRS (tests)
+
+# The source is organized into subdirectories, but we handle them all from
+# this CMakeLists file rather than letting CMake visit them as SUBDIRS.
+SET(OPENMM_SOURCE_SUBDIRS .)
+
+
+# Collect up information about the version of the OpenMM library we're building
+# and make it available to the code so it can be built into the binaries.
+
+SET(OPENMMCPU_LIBRARY_NAME OpenMMCPU)
+
+SET(SHARED_TARGET ${OPENMMCPU_LIBRARY_NAME})
+SET(STATIC_TARGET ${OPENMMCPU_LIBRARY_NAME}_static)
+
+
+# Ensure that debug libraries have "_d" appended to their names.
+# CMake gets this right on Windows automatically with this definition.
+IF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
+    SET(CMAKE_DEBUG_POSTFIX "_d" CACHE INTERNAL "" FORCE)
+ENDIF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
+
+# But on Unix or Cygwin we have to add the suffix manually
+IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+    SET(SHARED_TARGET ${SHARED_TARGET}_d)
+    SET(STATIC_TARGET ${STATIC_TARGET}_d)
+ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+
+
+# These are all the places to search for header files which are
+# to be part of the API.
+SET(API_INCLUDE_DIRS) # start empty
+FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
+    # append
+    SET(API_INCLUDE_DIRS ${API_INCLUDE_DIRS}
+                         ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include
+                         ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include/internal)
+ENDFOREACH(subdir)
+
+# We'll need both *relative* path names, starting with their API_INCLUDE_DIRS,
+# and absolute pathnames.
+SET(API_REL_INCLUDE_FILES)   # start these out empty
+SET(API_ABS_INCLUDE_FILES)
+
+FOREACH(dir ${API_INCLUDE_DIRS})
+    FILE(GLOB fullpaths ${dir}/*.h)	# returns full pathnames
+    SET(API_ABS_INCLUDE_FILES ${API_ABS_INCLUDE_FILES} ${fullpaths})
+
+    FOREACH(pathname ${fullpaths})
+        GET_FILENAME_COMPONENT(filename ${pathname} NAME)
+        SET(API_REL_INCLUDE_FILES ${API_REL_INCLUDE_FILES} ${dir}/${filename})
+    ENDFOREACH(pathname)
+ENDFOREACH(dir)
+
+# collect up source files
+SET(SOURCE_FILES) # empty
+SET(SOURCE_INCLUDE_FILES)
+
+FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
+    FILE(GLOB_RECURSE src_files  ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.c)
+    FILE(GLOB incl_files ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.h)
+    SET(SOURCE_FILES         ${SOURCE_FILES}         ${src_files})   #append
+    IF(MSVC)
+         FILE(GLOB_RECURSE kernel_files ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/kernels/*.cu)
+         SET(SOURCE_FILES ${SOURCE_FILES} ${kernel_files})
+    ENDIF(MSVC)
+    SET(SOURCE_INCLUDE_FILES ${SOURCE_INCLUDE_FILES} ${incl_files})
+    INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include)
+ENDFOREACH(subdir)
+
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src)
+
+# Install headers
+
+FILE(GLOB CORE_HEADERS include/*.h)
+INSTALL_FILES(/include/openmm/cpu FILES ${CORE_HEADERS})
+
+SUBDIRS (sharedTarget)

platforms/cpu/include/CpuKernelFactory.h

+#ifndef OPENMM_CPUKERNELFACTORY_H_
+#define OPENMM_CPUKERNELFACTORY_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: Peter Eastman                                                     *
+ * 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 "openmm/KernelFactory.h"
+
+namespace OpenMM {
+
+/**
+ * This KernelFactory creates all kernels for CpuPlatform.
+ */
+
+class CpuKernelFactory : public KernelFactory {
+public:
+    KernelImpl* createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CPUKERNELFACTORY_H_*/

platforms/cpu/include/CpuKernels.h

+#ifndef OPENMM_CPUKERNELS_H_
+#define OPENMM_CPUKERNELS_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: Peter Eastman                                                     *
+ * 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 "CpuPlatform.h"
+#include "CpuNeighborList.h"
+#include "openmm/kernels.h"
+#include "openmm/System.h"
+
+namespace OpenMM {
+
+/**
+ * This kernel is invoked by NonbondedForce to calculate the forces acting on the system.
+ */
+class CpuCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
+public:
+    CpuCalcNonbondedForceKernel(std::string name, const Platform& platform) : CalcNonbondedForceKernel(name, platform) {
+    }
+    ~CpuCalcNonbondedForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the NonbondedForce this kernel will be used for
+     */
+    void initialize(const System& system, const NonbondedForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @param includeDirect  true if direct space interactions should be included
+     * @param includeReciprocal  true if reciprocal space interactions should be included
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the NonbondedForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const NonbondedForce& force);
+private:
+    int numParticles, num14;
+    int **bonded14IndexArray;
+    double **particleParamArray, **bonded14ParamArray;
+    double nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, dispersionCoefficient;
+    int kmax[3], gridSize[3];
+    bool useSwitchingFunction;
+    std::vector<std::set<int> > exclusions;
+    std::vector<float> posq;
+    NonbondedMethod nonbondedMethod;
+    CpuNeighborList neighborList;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CPUKERNELS_H_*/
+

platforms/cpu/include/CpuNeighborList.h

+#ifndef OPENMM_CPU_NEIGHBORLIST_H_
+#define OPENMM_CPU_NEIGHBORLIST_H_
+
+#include "windowsExportCpu.h"
+#include <set>
+#include <utility>
+#include <vector>
+
+namespace OpenMM {
+    
+class OPENMM_EXPORT_CPU CpuNeighborList {
+public:
+    void computeNeighborList(int nAtoms,
+                             const std::vector<float>& atomLocations, 
+                             const std::vector<std::set<int> >& exclusions,
+                             const float* periodicBoxSize,
+                             bool usePeriodic,
+                             float maxDistance,
+                             float minDistance = 0.0f,
+                             bool reportSymmetricPairs = false);
+    const std::vector<std::pair<int, int> >& getNeighbors();
+private:
+    std::vector<std::pair<int, int> > neighbors;
+};
+
+} // namespace OpenMM
+
+#endif // OPENMM_REFERENCE_NEIGHBORLIST_H_

platforms/cpu/include/CpuPlatform.h

+#ifndef OPENMM_CPUPLATFORM_H_
+#define OPENMM_CPUPLATFORM_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: Peter Eastman                                                     *
+ * 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 "ReferencePlatform.h"
+#include "windowsExportCpu.h"
+
+namespace OpenMM {
+    
+/**
+ * This Platform subclass uses CPU implementations of the OpenMM kernels.
+ */
+
+class OPENMM_EXPORT_CPU CpuPlatform : public ReferencePlatform {
+public:
+    CpuPlatform();
+    const std::string& getName() const {
+        static const std::string name = "CPU";
+        return name;
+    }
+    double getSpeed() const;
+    bool supportsDoublePrecision() const;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CPUPLATFORM_H_*/

platforms/cpu/include/windowsExportCpu.h

+#ifndef OPENMM_WINDOWSEXPORTCPU_H_
+#define OPENMM_WINDOWSEXPORTCPU_H_
+
+/*
+ * Shared libraries are messy in Visual Studio. We have to distinguish three
+ * cases:
+ *   (1) this header is being used to build the OpenMM shared library
+ *       (dllexport)
+ *   (2) this header is being used by a *client* of the OpenMM shared
+ *       library (dllimport)
+ *   (3) we are building the OpenMM static library, or the client is
+ *       being compiled with the expectation of linking with the
+ *       OpenMM static library (nothing special needed)
+ * In the CMake script for building this library, we define one of the symbols
+ *     OPENMM_CPU_BUILDING_{SHARED|STATIC}_LIBRARY
+ * Client code normally has no special symbol defined, in which case we'll
+ * assume it wants to use the shared library. However, if the client defines
+ * the symbol OPENMM_USE_STATIC_LIBRARIES we'll suppress the dllimport so
+ * that the client code can be linked with static libraries. Note that
+ * the client symbol is not library dependent, while the library symbols
+ * affect only the OpenMM library, meaning that other libraries can
+ * be clients of this one. However, we are assuming all-static or all-shared.
+ */
+
+#ifdef _MSC_VER
+    // We don't want to hear about how sprintf is "unsafe".
+    #pragma warning(disable:4996)
+    // Keep MS VC++ quiet about lack of dll export of private members.
+    #pragma warning(disable:4251)
+    #if defined(OPENMM_CPU_BUILDING_SHARED_LIBRARY)
+        #define OPENMM_EXPORT_CPU __declspec(dllexport)
+    #elif defined(OPENMM_CPU_BUILDING_STATIC_LIBRARY) || defined(OPENMM_CPU_USE_STATIC_LIBRARIES)
+        #define OPENMM_EXPORT_CPU
+    #else
+        #define OPENMM_EXPORT_CPU __declspec(dllimport)   // i.e., a client of a shared library
+    #endif
+#else
+    #define OPENMM_EXPORT_CPU // Linux, Mac
+#endif
+
+#endif // OPENMM_WINDOWSEXPORTCPU_H_

platforms/cpu/sharedTarget/CMakeLists.txt

+ADD_LIBRARY(${SHARED_TARGET} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
+
+IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+    SET(MAIN_OPENMM_LIB ${OPENMM_LIBRARY_NAME}_d)
+ELSE (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+    SET(MAIN_OPENMM_LIB ${OPENMM_LIBRARY_NAME})
+ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+TARGET_LINK_LIBRARIES(${SHARED_TARGET} ${MAIN_OPENMM_LIB} ${PTHREADS_LIB})
+SET_TARGET_PROPERTIES(${SHARED_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_CPU_BUILDING_SHARED_LIBRARY")
+
+INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${SHARED_TARGET})

platforms/cpu/src/CpuKernelFactory.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: Peter Eastman                                                     *
+ * 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 "CpuKernelFactory.h"
+#include "CpuKernels.h"
+#include "CpuPlatform.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/OpenMMException.h"
+
+using namespace OpenMM;
+
+KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
+    ReferencePlatform::PlatformData& data = *static_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    if (name == CalcNonbondedForceKernel::Name())
+        return new CpuCalcNonbondedForceKernel(name, platform);
+    throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '") + name + "'").c_str());
+}

platforms/cpu/src/CpuKernels.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: Peter Eastman                                                     *
+ * 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 "CpuKernels.h"
+#include "openmm/Context.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/NonbondedForceImpl.h"
+#include "RealVec.h"
+
+using namespace OpenMM;
+using namespace std;
+
+static vector<RealVec>& extractPositions(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return *((vector<RealVec>*) data->positions);
+}
+
+static vector<RealVec>& extractVelocities(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return *((vector<RealVec>*) data->velocities);
+}
+
+static vector<RealVec>& extractForces(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return *((vector<RealVec>*) data->forces);
+}
+
+static RealVec& extractBoxSize(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return *(RealVec*) data->periodicBoxSize;
+}
+
+CpuCalcNonbondedForceKernel::~CpuCalcNonbondedForceKernel() {
+}
+
+void CpuCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
+
+    // Identify which exceptions are 1-4 interactions.
+
+    numParticles = force.getNumParticles();
+    exclusions.resize(numParticles);
+    vector<int> nb14s;
+    for (int i = 0; i < force.getNumExceptions(); i++) {
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
+        exclusions[particle1].insert(particle2);
+        exclusions[particle2].insert(particle1);
+        if (chargeProd != 0.0 || epsilon != 0.0)
+            nb14s.push_back(i);
+    }
+
+    // Record the particle parameters.
+
+    num14 = nb14s.size();
+//    bonded14IndexArray = allocateIntArray(num14, 2);
+//    bonded14ParamArray = allocateRealArray(num14, 3);
+//    particleParamArray = allocateRealArray(numParticles, 3);
+    posq.resize(4*numParticles, 0);
+    for (int i = 0; i < numParticles; ++i) {
+        double charge, radius, depth;
+        force.getParticleParameters(i, charge, radius, depth);
+        posq[4*i+3] = (float) charge;
+//        particleParamArray[i][0] = static_cast<RealOpenMM>(0.5*radius);
+//        particleParamArray[i][1] = static_cast<RealOpenMM>(2.0*sqrt(depth));
+//        particleParamArray[i][2] = static_cast<RealOpenMM>(charge);
+    }
+//    this->exclusions = exclusions;
+//    for (int i = 0; i < num14; ++i) {
+//        int particle1, particle2;
+//        double charge, radius, depth;
+//        force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth);
+//        bonded14IndexArray[i][0] = particle1;
+//        bonded14IndexArray[i][1] = particle2;
+//        bonded14ParamArray[i][0] = static_cast<RealOpenMM>(radius);
+//        bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
+//        bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge);
+//    }
+    nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
+    nonbondedCutoff = force.getCutoffDistance();
+    if (nonbondedMethod == NoCutoff)
+        useSwitchingFunction = false;
+    else {
+        useSwitchingFunction = force.getUseSwitchingFunction();
+        switchingDistance = force.getSwitchingDistance();
+    }
+    if (nonbondedMethod == Ewald) {
+        double alpha;
+        NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmax[0], kmax[1], kmax[2]);
+        ewaldAlpha = alpha;
+    }
+    else if (nonbondedMethod == PME) {
+        double alpha;
+        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2]);
+        ewaldAlpha = alpha;
+    }
+    rfDielectric = force.getReactionFieldDielectric();
+    if (force.getUseDispersionCorrection())
+        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
+    else
+        dispersionCoefficient = 0.0;
+}
+
+double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
+    vector<RealVec>& posData = extractPositions(context);
+    vector<RealVec>& forceData = extractForces(context);
+    RealVec boxSize = extractBoxSize(context);
+    float floatBoxSize[3] = {(float) boxSize[0], (float) boxSize[1], (float) boxSize[2]};
+    double energy = 0;
+//    CpuLJCoulombIxn clj;
+    bool periodic = (nonbondedMethod == CutoffPeriodic);
+    bool ewald  = (nonbondedMethod == Ewald);
+    bool pme  = (nonbondedMethod == PME);
+    
+    // Convert the positions to single precision.
+    
+    if (periodic)
+        for (int i = 0; i < numParticles; i++)
+            for (int j = 0; j < 3; j++) {
+                RealOpenMM x = posData[i][j];
+                double base = floor(x/boxSize[j]+0.5)*boxSize[j];
+                posq[4*i+j] = (float) (x-base);
+            }
+    else
+        for (int i = 0; i < numParticles; i++) {
+            posq[4*i] = (float) posData[i][0];
+            posq[4*i+1] = (float) posData[i][1];
+            posq[4*i+2] = (float) posData[i][2];
+        }
+    neighborList.computeNeighborList(numParticles, posq, exclusions, floatBoxSize, periodic || ewald || pme, nonbondedCutoff, 0.0);
+//    if (nonbondedMethod != NoCutoff) {
+//        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic || ewald || pme, nonbondedCutoff, 0.0);
+//        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
+//    }
+//    if (periodic || ewald || pme) {
+//        RealVec& box = extractBoxSize(context);
+//        double minAllowedSize = 1.999999*nonbondedCutoff;
+//        if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
+//            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
+//        clj.setPeriodic(box);
+//    }
+//    if (ewald)
+//        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
+//    if (pme)
+//        clj.setUsePME(ewaldAlpha, gridSize);
+//    if (useSwitchingFunction)
+//        clj.setUseSwitchingFunction(switchingDistance);
+//    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusions, 0, forceData, 0, includeEnergy ? &energy : NULL, includeDirect, includeReciprocal);
+//    if (includeDirect) {
+//        CpuBondForce refBondForce;
+//        CpuLJCoulomb14 nonbonded14;
+//        refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
+//        if (periodic || ewald || pme) {
+//            RealVec& boxSize = extractBoxSize(context);
+//            energy += dispersionCoefficient/(boxSize[0]*boxSize[1]*boxSize[2]);
+//        }
+//    }
+//    return energy;
+    return 0.0;
+}
+
+void CpuCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const NonbondedForce& force) {
+//    if (force.getNumParticles() != numParticles)
+//        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+//    vector<int> nb14s;
+//    for (int i = 0; i < force.getNumExceptions(); i++) {
+//        int particle1, particle2;
+//        double chargeProd, sigma, epsilon;
+//        force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
+//        if (chargeProd != 0.0 || epsilon != 0.0)
+//            nb14s.push_back(i);
+//    }
+//    if (nb14s.size() != num14)
+//        throw OpenMMException("updateParametersInContext: The number of non-excluded exceptions has changed");
+//
+//    // Record the values.
+//
+//    for (int i = 0; i < numParticles; ++i) {
+//        double charge, radius, depth;
+//        force.getParticleParameters(i, charge, radius, depth);
+//        particleParamArray[i][0] = static_cast<RealOpenMM>(0.5*radius);
+//        particleParamArray[i][1] = static_cast<RealOpenMM>(2.0*sqrt(depth));
+//        particleParamArray[i][2] = static_cast<RealOpenMM>(charge);
+//    }
+//    for (int i = 0; i < num14; ++i) {
+//        int particle1, particle2;
+//        double charge, radius, depth;
+//        force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth);
+//        bonded14IndexArray[i][0] = particle1;
+//        bonded14IndexArray[i][1] = particle2;
+//        bonded14ParamArray[i][0] = static_cast<RealOpenMM>(radius);
+//        bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
+//        bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge);
+//    }
+//    
+//    // Recompute the coefficient for the dispersion correction.
+//
+//    NonbondedForce::NonbondedMethod method = force.getNonbondedMethod();
+//    if (force.getUseDispersionCorrection() && (method == NonbondedForce::CutoffPeriodic || method == NonbondedForce::Ewald || method == NonbondedForce::PME))
+//        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
+}

platforms/cpu/src/CpuNeighborList.cpp

+#include "CpuNeighborList.h"
+#include <set>
+#include <map>
+#include <cmath>
+
+using namespace std;
+
+namespace OpenMM {
+
+static float periodicDifference(float val1, float val2, float period) {
+    float diff = val1-val2;
+    float base = floorf(diff/period+0.5f)*period;
+    return diff-base;
+}
+
+// squared distance between two points
+static float compPairDistanceSquared(const float* pos1, const float* pos2, const float* periodicBoxSize, bool usePeriodic) {
+    float dx, dy, dz;
+    if (!usePeriodic) {
+        dx = pos2[0] - pos1[0];
+        dy = pos2[1] - pos1[1];
+        dz = pos2[2] - pos1[2];
+    }
+    else {
+        dx = periodicDifference(pos2[0], pos1[0], periodicBoxSize[0]);
+        dy = periodicDifference(pos2[1], pos1[1], periodicBoxSize[1]);
+        dz = periodicDifference(pos2[2], pos1[2], periodicBoxSize[2]);
+    }
+    return dx*dx + dy*dy + dz*dz;
+}
+
+class VoxelIndex 
+{
+public:
+    VoxelIndex(int xx, int yy, int zz) : x(xx), y(yy), z(zz) {}
+
+    // operator<() needed for map
+    bool operator<(const VoxelIndex& other) const {
+        if      (x < other.x) return true;
+        else if (x > other.x) return false;
+        else if (y < other.y) return true;
+        else if (y > other.y) return false;
+        else if (z < other.z) return true;
+        else return false;
+    }
+
+    int x;
+    int y;
+    int z;
+};
+
+typedef std::pair<const float*, int> VoxelItem;
+typedef std::vector< VoxelItem > Voxel;
+
+class VoxelHash
+{
+public:
+    VoxelHash(float vsx, float vsy, float vsz, const float* periodicBoxSize, bool usePeriodic) :
+            voxelSizeX(vsx), voxelSizeY(vsy), voxelSizeZ(vsz), periodicBoxSize(periodicBoxSize), usePeriodic(usePeriodic) {
+        if (usePeriodic) {
+            nx = (int) floorf(periodicBoxSize[0]/voxelSizeX+0.5f);
+            ny = (int) floorf(periodicBoxSize[1]/voxelSizeY+0.5f);
+            nz = (int) floorf(periodicBoxSize[2]/voxelSizeZ+0.5f);
+        }
+    }
+
+    void insert(const int& item, const float* location)
+    {
+        VoxelIndex voxelIndex = getVoxelIndex(location);
+        if (voxelMap.find(voxelIndex) == voxelMap.end()) voxelMap[voxelIndex] = Voxel(); 
+        Voxel& voxel = voxelMap.find(voxelIndex)->second;
+        voxel.push_back(VoxelItem(location, item));
+    }
+
+
+    VoxelIndex getVoxelIndex(const float* location) const {
+        float xperiodic, yperiodic, zperiodic;
+        if (!usePeriodic) {
+            xperiodic = location[0];
+            yperiodic = location[1];
+            zperiodic = location[2];
+        }
+        else {
+            xperiodic = location[0]-periodicBoxSize[0]*floorf(location[0]/periodicBoxSize[0]);
+            yperiodic = location[1]-periodicBoxSize[1]*floorf(location[1]/periodicBoxSize[1]);
+            zperiodic = location[2]-periodicBoxSize[2]*floorf(location[2]/periodicBoxSize[2]);
+        }
+        int x = int(floorf(xperiodic / voxelSizeX));
+        int y = int(floorf(yperiodic / voxelSizeY));
+        int z = int(floorf(zperiodic / voxelSizeZ));
+        
+        return VoxelIndex(x, y, z);
+    }
+
+    void getNeighbors(
+            vector<pair<int, int> >& neighbors, 
+            const VoxelItem& referencePoint, 
+            const vector<set<int> >& exclusions,
+            bool reportSymmetricPairs,
+            float maxDistance, 
+            float minDistance) const 
+    {
+
+        // Loop over neighboring voxels
+        // TODO use more clever selection of neighboring voxels
+
+        const int atomI = referencePoint.second;
+        const float* locationI = referencePoint.first;
+        
+        float maxDistanceSquared = maxDistance * maxDistance;
+        float minDistanceSquared = minDistance * minDistance;
+
+        int dIndexX = int(maxDistance / voxelSizeX) + 1; // How may voxels away do we have to look?
+        int dIndexY = int(maxDistance / voxelSizeY) + 1;
+        int dIndexZ = int(maxDistance / voxelSizeZ) + 1;
+        VoxelIndex centerVoxelIndex = getVoxelIndex(locationI);
+        int lastx = centerVoxelIndex.x+dIndexX;
+        int lasty = centerVoxelIndex.y+dIndexY;
+        int lastz = centerVoxelIndex.z+dIndexZ;
+        if (usePeriodic) {
+            lastx = min(lastx, centerVoxelIndex.x-dIndexX+nx-1);
+            lasty = min(lasty, centerVoxelIndex.y-dIndexY+ny-1);
+            lastz = min(lastz, centerVoxelIndex.z-dIndexZ+nz-1);
+        }
+        for (int x = centerVoxelIndex.x - dIndexX; x <= lastx; ++x)
+        {
+            for (int y = centerVoxelIndex.y - dIndexY; y <= lasty; ++y)
+            {
+                for (int z = centerVoxelIndex.z - dIndexZ; z <= lastz; ++z)
+                {
+                    VoxelIndex voxelIndex(x, y, z);
+                    if (usePeriodic) {
+                        voxelIndex.x = (x+nx)%nx;
+                        voxelIndex.y = (y+ny)%ny;
+                        voxelIndex.z = (z+nz)%nz;
+                    }
+                    if (voxelMap.find(voxelIndex) == voxelMap.end()) continue; // no such voxel; skip
+                    const Voxel& voxel = voxelMap.find(voxelIndex)->second;
+                    for (Voxel::const_iterator itemIter = voxel.begin(); itemIter != voxel.end(); ++itemIter)
+                    {
+                        const int atomJ = itemIter->second;
+                        const float* locationJ = itemIter->first;
+                        
+                        // Ignore self hits
+                        if (atomI == atomJ) continue;
+                        
+                        // Ignore exclusions.
+                        if (exclusions[atomI].find(atomJ) != exclusions[atomI].end()) continue;
+                        
+                        float dSquared = compPairDistanceSquared(locationI, locationJ, periodicBoxSize, usePeriodic);
+                        if (dSquared > maxDistanceSquared) continue;
+                        if (dSquared < minDistanceSquared) continue;
+                        
+                        neighbors.push_back(make_pair(atomI, atomJ));
+                        if (reportSymmetricPairs)
+                            neighbors.push_back(make_pair(atomJ, atomI));
+                    }
+                }
+            }
+        }
+    }
+
+private:
+    float voxelSizeX, voxelSizeY, voxelSizeZ;
+    int nx, ny, nz;
+    const float* periodicBoxSize;
+    const bool usePeriodic;
+    std::map<VoxelIndex, Voxel> voxelMap;
+};
+
+
+// O(n) neighbor list method using voxel hash data structure
+void CpuNeighborList::computeNeighborList(
+                              int nAtoms,
+                              const vector<float>& atomLocations, 
+                              const vector<set<int> >& exclusions,
+                              const float* periodicBoxSize,
+                              bool usePeriodic,
+                              float maxDistance,
+                              float minDistance,
+                              bool reportSymmetricPairs)
+{
+    neighbors.clear();
+
+    float edgeSizeX, edgeSizeY, edgeSizeZ;
+    if (!usePeriodic)
+        edgeSizeX = edgeSizeY = edgeSizeZ = maxDistance; // TODO - adjust this as needed
+    else {
+        edgeSizeX = periodicBoxSize[0]/floorf(periodicBoxSize[0]/maxDistance);
+        edgeSizeY = periodicBoxSize[1]/floorf(periodicBoxSize[1]/maxDistance);
+        edgeSizeZ = periodicBoxSize[2]/floorf(periodicBoxSize[2]/maxDistance);
+    }
+    VoxelHash voxelHash(edgeSizeX, edgeSizeY, edgeSizeZ, periodicBoxSize, usePeriodic);
+    for (int atomJ = 0; atomJ < (int) nAtoms; ++atomJ) // use "j", because j > i for pairs
+    {
+        // 1) Find other atoms that are close to this one
+        const float location[3] = {atomLocations[4*atomJ], atomLocations[4*atomJ+1], atomLocations[4*atomJ+2]};
+        voxelHash.getNeighbors(
+            neighbors, 
+            VoxelItem(location, atomJ),
+            exclusions,
+            reportSymmetricPairs, 
+            maxDistance, 
+            minDistance);
+            
+        // 2) Add this atom to the voxelHash
+        voxelHash.insert(atomJ, location);
+    }
+}
+
+const vector<pair<int, int> >& CpuNeighborList::getNeighbors() {
+    return neighbors;
+}
+
+} // namespace OpenMM

platforms/cpu/src/CpuPlatform.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: Peter Eastman                                                     *
+ * 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 "CpuPlatform.h"
+#include "CpuKernelFactory.h"
+#include "CpuKernels.h"
+
+using namespace OpenMM;
+
+extern "C" OPENMM_EXPORT_CPU void registerPlatforms() {
+    Platform::registerPlatform(new CpuPlatform());
+}
+
+CpuPlatform::CpuPlatform() {
+    CpuKernelFactory* factory = new CpuKernelFactory();
+    registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
+}
+
+double CpuPlatform::getSpeed() const {
+    return 10;
+}
+
+bool CpuPlatform::supportsDoublePrecision() const {
+    return false;
+}

platforms/cpu/tests/CMakeLists.txt

+#
+# Testing
+#
+
+ENABLE_TESTING()
+
+SET( INCLUDE_SERIALIZATION FALSE )
+#SET( INCLUDE_SERIALIZATION TRUE )
+
+IF( INCLUDE_SERIALIZATION )
+    INCLUDE_DIRECTORIES(${OPENMM_DIR}/serialization/include)
+    SET( SHARED_OPENMM_SERIALIZATION "OpenMMSerialization" )
+ENDIF( INCLUDE_SERIALIZATION )
+
+# Automatically create tests using files named "Test*.cpp"
+FILE(GLOB TEST_PROGS "*Test*.cpp")
+FOREACH(TEST_PROG ${TEST_PROGS})
+    GET_FILENAME_COMPONENT(TEST_ROOT ${TEST_PROG} NAME_WE)
+
+    # Link with shared library
+    ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
+    TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET})
+    ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} single)
+
+ENDFOREACH(TEST_PROG ${TEST_PROGS})

platforms/cpu/tests/TestCpuEwald.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) 2008-2013 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Ewald summation method CPU implementation of NonbondedForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CpuPlatform.h"
+#include "ReferencePlatform.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include "openmm/VerletIntegrator.h"
+#include "openmm/internal/ContextImpl.h"
+#include "SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+CpuPlatform platform;
+
+const double TOL = 1e-5;
+
+void testEwaldPME(bool includeExceptions) {
+
+//      Use amorphous NaCl system for the tests
+
+    const int numParticles = 894;
+    const double cutoff = 1.2;
+    const double boxSize = 3.00646;
+    double tol = 1e-5;
+
+    ReferencePlatform reference;
+    System system;
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setEwaldErrorTolerance(tol);
+
+    for (int i = 0; i < numParticles/2; i++)
+        system.addParticle(22.99);
+    for (int i = 0; i < numParticles/2; i++)
+        system.addParticle(35.45);
+    for (int i = 0; i < numParticles/2; i++)
+        nonbonded->addParticle(1.0, 1.0,0.0);
+    for (int i = 0; i < numParticles/2; i++)
+        nonbonded->addParticle(-1.0, 1.0,0.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(nonbonded);
+
+    vector<Vec3> positions(numParticles);
+    #include "nacl_amorph.dat"
+    if (includeExceptions) {
+        // Add some exclusions.
+
+        for (int i = 0; i < numParticles-1; i++) {
+            Vec3 delta = positions[i]-positions[i+1];
+            if (sqrt(delta.dot(delta)) < 0.5*cutoff)
+                nonbonded->addException(i, i+1, i%2 == 0 ? 0.0 : 0.5, 1.0, 0.0);
+        }
+    }
+
+//    (1)  Check whether the Reference and CUDA platforms agree when using Ewald Method
+
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context cuContext(system, integrator1, platform);
+    Context referenceContext(system, integrator2, reference);
+    cuContext.setPositions(positions);
+    referenceContext.setPositions(positions);
+    State cuState = cuContext.getState(State::Forces | State::Energy);
+    State referenceState = referenceContext.getState(State::Forces | State::Energy);
+    tol = 1e-2;
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(referenceState.getForces()[i], cuState.getForces()[i], tol);
+    }
+    tol = 1e-5;
+    ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), cuState.getPotentialEnergy(), tol);
+
+//    (2) Check whether Ewald method in CUDA is self-consistent
+
+    double norm = 0.0;
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 f = cuState.getForces()[i];
+        norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
+    }
+
+    norm = std::sqrt(norm);
+    const double delta = 5e-3;
+    double step = delta/norm;
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = cuState.getForces()[i];
+        positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+    }
+    VerletIntegrator integrator3(0.01);
+    Context cuContext2(system, integrator3, platform);
+    cuContext2.setPositions(positions);
+
+    tol = 1e-2;
+    State cuState2 = cuContext2.getState(State::Energy);
+    ASSERT_EQUAL_TOL(norm, (cuState2.getPotentialEnergy()-cuState.getPotentialEnergy())/delta, tol)
+
+//    (3)  Check whether the Reference and CUDA platforms agree when using PME
+
+    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    cuContext.reinitialize();
+    referenceContext.reinitialize();
+    cuContext.setPositions(positions);
+    referenceContext.setPositions(positions);
+    cuState = cuContext.getState(State::Forces | State::Energy);
+    referenceState = referenceContext.getState(State::Forces | State::Energy);
+    tol = 1e-2;
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(referenceState.getForces()[i], cuState.getForces()[i], tol);
+    }
+    tol = 1e-5;
+    ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), cuState.getPotentialEnergy(), tol);
+
+//    (4) Check whether PME method in CUDA is self-consistent
+
+    norm = 0.0;
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 f = cuState.getForces()[i];
+        norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
+    }
+
+    norm = std::sqrt(norm);
+    step = delta/norm;
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = cuState.getForces()[i];
+        positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+    }
+    VerletIntegrator integrator4(0.01);
+    Context cuContext3(system, integrator4, platform);
+    cuContext3.setPositions(positions);
+
+    tol = 1e-2;
+    State cuState3 = cuContext3.getState(State::Energy);
+    ASSERT_EQUAL_TOL(norm, (cuState3.getPotentialEnergy()-cuState.getPotentialEnergy())/delta, tol)
+}
+
+void testEwald2Ions() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->addParticle(-1.0, 1, 0);
+    nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
+    const double cutoff = 2.0;
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setEwaldErrorTolerance(TOL);
+    system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(3.048000,2.764000,3.156000);
+    positions[1] = Vec3(2.809000,2.888000,2.571000);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+
+    ASSERT_EQUAL_VEC(Vec3(-123.711,  64.1877, -302.716), forces[0], 10*TOL);
+    ASSERT_EQUAL_VEC(Vec3( 123.711, -64.1877,  302.716), forces[1], 10*TOL);
+    ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 0.01/*10*TOL*/);
+}
+
+void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
+    // Create a cloud of random point charges.
+
+    const int numParticles = 51;
+    const double boxWidth = 5.0;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
+    NonbondedForce* force = new NonbondedForce();
+    system.addForce(force);
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
+        positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
+    }
+    force->setNonbondedMethod(method);
+
+    // For various values of the cutoff and error tolerance, see if the actual error is reasonable.
+
+    for (double cutoff = 1.0; cutoff < boxWidth/2; cutoff *= 1.2) {
+        force->setCutoffDistance(cutoff);
+        vector<Vec3> refForces;
+        double norm = 0.0;
+        for (double tol = 5e-5; tol < 1e-3; tol *= 2.0) {
+            force->setEwaldErrorTolerance(tol);
+            VerletIntegrator integrator(0.01);
+            Context context(system, integrator, platform);
+            context.setPositions(positions);
+            State state = context.getState(State::Forces);
+            if (refForces.size() == 0) {
+                refForces = state.getForces();
+                for (int i = 0; i < numParticles; i++)
+                    norm += refForces[i].dot(refForces[i]);
+                norm = sqrt(norm);
+            }
+            else {
+                double diff = 0.0;
+                for (int i = 0; i < numParticles; i++) {
+                    Vec3 delta = refForces[i]-state.getForces()[i];
+                    diff += delta.dot(delta);
+                }
+                diff = sqrt(diff)/norm;
+                ASSERT(diff < 2*tol);
+            }
+        }
+    }
+}
+
+int main(int argc, char* argv[]) {
+    try {
+        testEwaldPME(false);
+        testEwaldPME(true);
+//        testEwald2Ions();
+        testErrorTolerance(NonbondedForce::Ewald);
+        testErrorTolerance(NonbondedForce::PME);
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cpu/tests/TestCpuNonbondedForce.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) 2008-2013 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests all the different force terms in the CUDA implementation of NonbondedForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CpuPlatform.h"
+#include "ReferencePlatform.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include "openmm/VerletIntegrator.h"
+#include "openmm/internal/ContextImpl.h"
+#include "SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+CpuPlatform platform;
+
+const double TOL = 1e-5;
+
+void testCoulomb() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* forceField = new NonbondedForce();
+    forceField->addParticle(0.5, 1, 0);
+    forceField->addParticle(-1.5, 1, 0);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double force = ONE_4PI_EPS0*(-0.75)/4.0;
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(-0.75)/2.0, state.getPotentialEnergy(), TOL);
+}
+
+void testLJ() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* forceField = new NonbondedForce();
+    forceField->addParticle(0, 1.2, 1);
+    forceField->addParticle(0, 1.4, 2);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double x = 1.3/2.0;
+    double eps = SQRT_TWO;
+    double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/2.0;
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0)), state.getPotentialEnergy(), TOL);
+}
+
+void testExclusionsAnd14() {
+    System system;
+    NonbondedForce* nonbonded = new NonbondedForce();
+    for (int i = 0; i < 5; ++i) {
+        system.addParticle(1.0);
+        nonbonded->addParticle(0, 1.5, 0);
+    }
+    vector<pair<int, int> > bonds;
+    bonds.push_back(pair<int, int>(0, 1));
+    bonds.push_back(pair<int, int>(1, 2));
+    bonds.push_back(pair<int, int>(2, 3));
+    bonds.push_back(pair<int, int>(3, 4));
+    nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
+    int first14, second14;
+    for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
+        if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
+            first14 = i;
+        if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
+            second14 = i;
+    }
+    system.addForce(nonbonded);
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    for (int i = 1; i < 5; ++i) {
+
+        // Test LJ forces
+
+        vector<Vec3> positions(5);
+        const double r = 1.0;
+        for (int j = 0; j < 5; ++j) {
+            nonbonded->setParticleParameters(j, 0, 1.5, 0);
+            positions[j] = Vec3(0, j, 0);
+        }
+        nonbonded->setParticleParameters(0, 0, 1.5, 1);
+        nonbonded->setParticleParameters(i, 0, 1.5, 1);
+        nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
+        nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
+        positions[i] = Vec3(r, 0, 0);
+        context.reinitialize();
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = state.getForces();
+        double x = 1.5/r;
+        double eps = 1.0;
+        double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
+        double energy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
+        if (i == 3) {
+            force *= 0.5;
+            energy *= 0.5;
+        }
+        if (i < 3) {
+            force = 0;
+            energy = 0;
+        }
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+
+        // Test Coulomb forces
+
+        nonbonded->setParticleParameters(0, 2, 1.5, 0);
+        nonbonded->setParticleParameters(i, 2, 1.5, 0);
+        nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
+        nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
+        context.reinitialize();
+        context.setPositions(positions);
+        state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces2 = state.getForces();
+        force = ONE_4PI_EPS0*4/(r*r);
+        energy = ONE_4PI_EPS0*4/r;
+        if (i == 3) {
+            force /= 1.2;
+            energy /= 1.2;
+        }
+        if (i < 3) {
+            force = 0;
+            energy = 0;
+        }
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+    }
+}
+
+void testCutoff() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* forceField = new NonbondedForce();
+    forceField->addParticle(1.0, 1, 0);
+    forceField->addParticle(1.0, 1, 0);
+    forceField->addParticle(1.0, 1, 0);
+    forceField->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    const double cutoff = 2.9;
+    forceField->setCutoffDistance(cutoff);
+    const double eps = 50.0;
+    forceField->setReactionFieldDielectric(eps);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 2, 0);
+    positions[2] = Vec3(0, 3, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
+    const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
+    const double force1 = ONE_4PI_EPS0*(1.0)*(0.25-2.0*krf*2.0);
+    const double force2 = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
+    ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL);
+    const double energy1 = ONE_4PI_EPS0*(1.0)*(0.5+krf*4.0-crf);
+    const double energy2 = ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf);
+    ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL);
+}
+
+void testCutoff14() {
+    System system;
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    for (int i = 0; i < 5; ++i) {
+        system.addParticle(1.0);
+        nonbonded->addParticle(0, 1.5, 0);
+    }
+    const double cutoff = 3.5;
+    nonbonded->setCutoffDistance(cutoff);
+    const double eps = 30.0;
+    nonbonded->setReactionFieldDielectric(eps);
+    vector<pair<int, int> > bonds;
+    bonds.push_back(pair<int, int>(0, 1));
+    bonds.push_back(pair<int, int>(1, 2));
+    bonds.push_back(pair<int, int>(2, 3));
+    bonds.push_back(pair<int, int>(3, 4));
+    nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
+    int first14, second14;
+    for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
+        if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
+            first14 = i;
+        if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
+            second14 = i;
+    }
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(5);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    positions[2] = Vec3(2, 0, 0);
+    positions[3] = Vec3(3, 0, 0);
+    positions[4] = Vec3(4, 0, 0);
+    for (int i = 1; i < 5; ++i) {
+
+        // Test LJ forces
+
+        nonbonded->setParticleParameters(0, 0, 1.5, 1);
+        for (int j = 1; j < 5; ++j)
+            nonbonded->setParticleParameters(j, 0, 1.5, 0);
+        nonbonded->setParticleParameters(i, 0, 1.5, 1);
+        nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
+        nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
+        context.reinitialize();
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = state.getForces();
+        double r = positions[i][0];
+        double x = 1.5/r;
+        double e = 1.0;
+        double force = 4.0*e*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
+        double energy = 4.0*e*(std::pow(x, 12.0)-std::pow(x, 6.0));
+        if (i == 3) {
+            force *= 0.5;
+            energy *= 0.5;
+        }
+        if (i < 3 || r > cutoff) {
+            force = 0;
+            energy = 0;
+        }
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+
+        // Test Coulomb forces
+
+        const double q = 0.7;
+        nonbonded->setParticleParameters(0, q, 1.5, 0);
+        nonbonded->setParticleParameters(i, q, 1.5, 0);
+        nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
+        nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
+        context.reinitialize();
+        context.setPositions(positions);
+        state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces2 = state.getForces();
+        force = ONE_4PI_EPS0*q*q/(r*r);
+        energy = ONE_4PI_EPS0*q*q/r;
+        if (i == 3) {
+            force /= 1.2;
+            energy /= 1.2;
+        }
+        if (i < 3 || r > cutoff) {
+            force = 0;
+            energy = 0;
+        }
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+    }
+}
+
+void testPeriodic() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->addException(0, 1, 0.0, 1.0, 0.0);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    const double cutoff = 2.0;
+    nonbonded->setCutoffDistance(cutoff);
+    system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    positions[2] = Vec3(3, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    const double eps = 78.3;
+    const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
+    const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
+    const double force = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
+}
+
+
+void testLargeSystem() {
+    const int numMolecules = 600;
+    const int numParticles = numMolecules*2;
+    const double cutoff = 2.0;
+    const double boxSize = 20.0;
+    const double tol = 2e-3;
+    ReferencePlatform reference;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < numMolecules; i++) {
+        if (i < numMolecules/2) {
+            nonbonded->addParticle(-1.0, 0.2, 0.1);
+            nonbonded->addParticle(1.0, 0.1, 0.1);
+        }
+        else {
+            nonbonded->addParticle(-1.0, 0.2, 0.2);
+            nonbonded->addParticle(1.0, 0.1, 0.2);
+        }
+        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        velocities[2*i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
+        velocities[2*i+1] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
+        bonds->addBond(2*i, 2*i+1, 1.0, 0.1);
+        nonbonded->addException(2*i, 2*i+1, 0.0, 0.15, 0.0);
+    }
+
+    // Try with cutoffs but not periodic boundary conditions, and make sure the cl and Reference
+    // platforms agree.
+
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    nonbonded->setCutoffDistance(cutoff);
+    system.addForce(nonbonded);
+    system.addForce(bonds);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context cuContext(system, integrator1, platform);
+    Context referenceContext(system, integrator2, reference);
+    cuContext.setPositions(positions);
+    cuContext.setVelocities(velocities);
+    referenceContext.setPositions(positions);
+    referenceContext.setVelocities(velocities);
+    State cuState = cuContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
+    State referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(cuState.getPositions()[i], referenceState.getPositions()[i], tol);
+        ASSERT_EQUAL_VEC(cuState.getVelocities()[i], referenceState.getVelocities()[i], tol);
+        ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
+    }
+    ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
+
+    // Now do the same thing with periodic boundary conditions.
+
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    cuContext.reinitialize();
+    referenceContext.reinitialize();
+    cuContext.setPositions(positions);
+    cuContext.setVelocities(velocities);
+    referenceContext.setPositions(positions);
+    referenceContext.setVelocities(velocities);
+    cuState = cuContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
+    referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
+    for (int i = 0; i < numParticles; i++) {
+        double dx = cuState.getPositions()[i][0]-referenceState.getPositions()[i][0];
+        double dy = cuState.getPositions()[i][1]-referenceState.getPositions()[i][1];
+        double dz = cuState.getPositions()[i][2]-referenceState.getPositions()[i][2];
+        ASSERT_EQUAL_TOL(fmod(cuState.getPositions()[i][0]-referenceState.getPositions()[i][0], boxSize), 0, tol);
+        ASSERT_EQUAL_TOL(fmod(cuState.getPositions()[i][1]-referenceState.getPositions()[i][1], boxSize), 0, tol);
+        ASSERT_EQUAL_TOL(fmod(cuState.getPositions()[i][2]-referenceState.getPositions()[i][2], boxSize), 0, tol);
+        ASSERT_EQUAL_VEC(cuState.getVelocities()[i], referenceState.getVelocities()[i], tol);
+        ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
+    }
+    ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
+}
+
+void testDispersionCorrection() {
+    // Create a box full of identical particles.
+
+    int gridSize = 5;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = gridSize*0.7;
+    double cutoff = boxSize/3;
+    System system;
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    vector<Vec3> positions(numParticles);
+    int index = 0;
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                system.addParticle(1.0);
+                nonbonded->addParticle(0, 1.1, 0.5);
+                positions[index] = Vec3(i*boxSize/gridSize, j*boxSize/gridSize, k*boxSize/gridSize);
+                index++;
+            }
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nonbonded->setCutoffDistance(cutoff);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(nonbonded);
+
+    // See if the correction has the correct value.
+
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    double energy1 = context.getState(State::Energy).getPotentialEnergy();
+    nonbonded->setUseDispersionCorrection(false);
+    context.reinitialize();
+    context.setPositions(positions);
+    double energy2 = context.getState(State::Energy).getPotentialEnergy();
+    double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
+    double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
+    double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
+    ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
+
+    // Now modify half the particles to be different, and see if it is still correct.
+
+    int numType2 = 0;
+    for (int i = 0; i < numParticles; i += 2) {
+        nonbonded->setParticleParameters(i, 0, 1, 1);
+        numType2++;
+    }
+    int numType1 = numParticles-numType2;
+    nonbonded->updateParametersInContext(context);
+    energy2 = context.getState(State::Energy).getPotentialEnergy();
+    nonbonded->setUseDispersionCorrection(true);
+    context.reinitialize();
+    context.setPositions(positions);
+    energy1 = context.getState(State::Energy).getPotentialEnergy();
+    term1 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
+    term2 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
+    term1 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 12)/pow(cutoff, 9))/9;
+    term2 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 6)/pow(cutoff, 3))/3;
+    double combinedSigma = 0.5*(1+1.1);
+    double combinedEpsilon = sqrt(1*0.5);
+    term1 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 12)/pow(cutoff, 9))/9;
+    term2 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 6)/pow(cutoff, 3))/3;
+    term1 /= (numParticles*(numParticles+1))/2;
+    term2 /= (numParticles*(numParticles+1))/2;
+    expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
+    ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
+}
+
+void testChangingParameters() {
+    const int numMolecules = 600;
+    const int numParticles = numMolecules*2;
+    const double cutoff = 2.0;
+    const double boxSize = 20.0;
+    const double tol = 2e-3;
+    ReferencePlatform reference;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < numMolecules; i++) {
+        if (i < numMolecules/2) {
+            nonbonded->addParticle(-1.0, 0.2, 0.1);
+            nonbonded->addParticle(1.0, 0.1, 0.1);
+        }
+        else {
+            nonbonded->addParticle(-1.0, 0.2, 0.2);
+            nonbonded->addParticle(1.0, 0.1, 0.2);
+        }
+        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        system.addConstraint(2*i, 2*i+1, 1.0);
+        nonbonded->addException(2*i, 2*i+1, 0.0, 0.15, 0.0);
+    }
+    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    nonbonded->setCutoffDistance(cutoff);
+    system.addForce(nonbonded);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    
+    // See if Reference and CPU give the same forces and energies.
+    
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context cuContext(system, integrator1, platform);
+    Context referenceContext(system, integrator2, reference);
+    cuContext.setPositions(positions);
+    referenceContext.setPositions(positions);
+    State cuState = cuContext.getState(State::Forces | State::Energy);
+    State referenceState = referenceContext.getState(State::Forces | State::Energy);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
+    ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
+    
+    // Now modify parameters and see if they still agree.
+
+    for (int i = 0; i < numParticles; i += 5) {
+        double charge, sigma, epsilon;
+        nonbonded->getParticleParameters(i, charge, sigma, epsilon);
+        nonbonded->setParticleParameters(i, 1.5*charge, 1.1*sigma, 1.7*epsilon);
+    }
+    nonbonded->updateParametersInContext(cuContext);
+    nonbonded->updateParametersInContext(referenceContext);
+    cuState = cuContext.getState(State::Forces | State::Energy);
+    referenceState = referenceContext.getState(State::Forces | State::Energy);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
+    ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
+}
+
+void testSwitchingFunction(NonbondedForce::NonbondedMethod method) {
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->addParticle(0, 1.2, 1);
+    nonbonded->addParticle(0, 1.4, 2);
+    nonbonded->setNonbondedMethod(method);
+    nonbonded->setCutoffDistance(2.0);
+    nonbonded->setUseSwitchingFunction(true);
+    nonbonded->setSwitchingDistance(1.5);
+    nonbonded->setUseDispersionCorrection(false);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    double eps = SQRT_TWO;
+    
+    // Compute the interaction at various distances.
+    
+    for (double r = 1.0; r < 2.5; r += 0.1) {
+        positions[1] = Vec3(r, 0, 0);
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        
+        // See if the energy is correct.
+        
+        double x = 1.3/r;
+        double expectedEnergy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
+        double switchValue;
+        if (r <= 1.5)
+            switchValue = 1;
+        else if (r >= 2.0)
+            switchValue = 0;
+        else {
+            double t = (r-1.5)/0.5;
+            switchValue = 1+t*t*t*(-10+t*(15-t*6));
+        }
+        ASSERT_EQUAL_TOL(switchValue*expectedEnergy, state.getPotentialEnergy(), TOL);
+        
+        // See if the force is the gradient of the energy.
+        
+        double delta = 1e-3;
+        positions[1] = Vec3(r-delta, 0, 0);
+        context.setPositions(positions);
+        double e1 = context.getState(State::Energy).getPotentialEnergy();
+        positions[1] = Vec3(r+delta, 0, 0);
+        context.setPositions(positions);
+        double e2 = context.getState(State::Energy).getPotentialEnergy();
+        ASSERT_EQUAL_TOL((e2-e1)/(2*delta), state.getForces()[0][0], 1e-3);
+    }
+}
+
+int main(int argc, char* argv[]) {
+    try {
+        testCoulomb();
+        testLJ();
+        testExclusionsAnd14();
+        testCutoff();
+        testCutoff14();
+        testPeriodic();
+        testLargeSystem();
+        testDispersionCorrection();
+        testChangingParameters();
+        testSwitchingFunction(NonbondedForce::CutoffNonPeriodic);
+        testSwitchingFunction(NonbondedForce::PME);
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cpu/tests/nacl_amorph.dat

+positions[0] = Vec3(1.066000,1.628000,0.835000);
+positions[1] = Vec3(1.072000,0.428000,0.190000);
+positions[2] = Vec3(0.524000,1.442000,1.160000);
+positions[3] = Vec3(2.383000,1.524000,1.119000);
+positions[4] = Vec3(0.390000,1.441000,0.575000);
+positions[5] = Vec3(0.618000,0.399000,0.819000);
+positions[6] = Vec3(1.003000,1.257000,1.543000);
+positions[7] = Vec3(2.933000,1.569000,0.642000);
+positions[8] = Vec3(0.849000,0.739000,0.089000);
+positions[9] = Vec3(0.060000,0.794000,0.766000);
+positions[10] = Vec3(1.652000,1.405000,1.010000);
+positions[11] = Vec3(2.843000,1.533000,1.781000);
+positions[12] = Vec3(0.952000,1.309000,0.996000);
+positions[13] = Vec3(1.847000,1.402000,0.313000);
+positions[14] = Vec3(2.674000,0.083000,1.691000);
+positions[15] = Vec3(1.763000,2.104000,0.728000);
+positions[16] = Vec3(0.914000,0.574000,0.982000);
+positions[17] = Vec3(0.514000,0.078000,0.891000);
+positions[18] = Vec3(0.538000,0.766000,1.110000);
+positions[19] = Vec3(0.808000,0.676000,0.570000);
+positions[20] = Vec3(0.178000,0.014000,0.628000);
+positions[21] = Vec3(1.329000,1.333000,0.339000);
+positions[22] = Vec3(1.029000,1.678000,0.503000);
+positions[23] = Vec3(1.423000,1.767000,1.104000);
+positions[24] = Vec3(1.966000,1.051000,0.282000);
+positions[25] = Vec3(1.596000,1.971000,0.194000);
+positions[26] = Vec3(1.025000,1.043000,2.809000);
+positions[27] = Vec3(1.628000,2.614000,0.088000);
+positions[28] = Vec3(0.440000,0.606000,0.141000);
+positions[29] = Vec3(1.050000,2.821000,2.517000);
+positions[30] = Vec3(0.644000,1.604000,0.770000);
+positions[31] = Vec3(0.637000,0.917000,0.392000);
+positions[32] = Vec3(0.611000,2.768000,0.013000);
+positions[33] = Vec3(1.892000,0.660000,0.473000);
+positions[34] = Vec3(1.052000,2.081000,0.982000);
+positions[35] = Vec3(1.508000,2.300000,0.439000);
+positions[36] = Vec3(2.617000,0.328000,1.099000);
+positions[37] = Vec3(0.910000,0.040000,0.259000);
+positions[38] = Vec3(1.195000,1.494000,1.202000);
+positions[39] = Vec3(2.657000,0.997000,0.564000);
+positions[40] = Vec3(1.465000,1.580000,0.648000);
+positions[41] = Vec3(0.154000,2.538000,1.331000);
+positions[42] = Vec3(0.849000,1.476000,1.365000);
+positions[43] = Vec3(0.898000,0.987000,1.178000);
+positions[44] = Vec3(0.958000,0.656000,1.358000);
+positions[45] = Vec3(1.067000,0.934000,0.211000);
+positions[46] = Vec3(1.030000,0.319000,1.281000);
+positions[47] = Vec3(2.709000,0.807000,0.240000);
+positions[48] = Vec3(0.837000,1.362000,0.588000);
+positions[49] = Vec3(2.080000,0.791000,2.947000);
+positions[50] = Vec3(0.200000,0.266000,1.474000);
+positions[51] = Vec3(0.848000,0.379000,1.625000);
+positions[52] = Vec3(0.637000,1.071000,0.821000);
+positions[53] = Vec3(1.324000,0.757000,2.951000);
+positions[54] = Vec3(2.666000,0.935000,1.373000);
+positions[55] = Vec3(1.584000,1.025000,1.703000);
+positions[56] = Vec3(1.699000,0.636000,0.038000);
+positions[57] = Vec3(1.099000,1.644000,1.879000);
+positions[58] = Vec3(2.897000,1.302000,1.522000);
+positions[59] = Vec3(1.753000,0.949000,2.885000);
+positions[60] = Vec3(2.502000,1.321000,0.752000);
+positions[61] = Vec3(0.545000,0.193000,1.959000);
+positions[62] = Vec3(1.098000,2.646000,1.706000);
+positions[63] = Vec3(0.001000,1.205000,0.670000);
+positions[64] = Vec3(2.997000,0.061000,1.040000);
+positions[65] = Vec3(0.662000,0.828000,1.535000);
+positions[66] = Vec3(1.252000,1.246000,0.780000);
+positions[67] = Vec3(1.173000,0.472000,0.810000);
+positions[68] = Vec3(0.124000,0.622000,2.992000);
+positions[69] = Vec3(1.036000,0.883000,0.848000);
+positions[70] = Vec3(1.423000,2.146000,1.340000);
+positions[71] = Vec3(2.391000,1.136000,1.165000);
+positions[72] = Vec3(1.189000,2.961000,0.425000);
+positions[73] = Vec3(1.584000,2.500000,0.782000);
+positions[74] = Vec3(0.565000,1.122000,1.240000);
+positions[75] = Vec3(1.733000,1.716000,1.763000);
+positions[76] = Vec3(1.548000,1.522000,0.041000);
+positions[77] = Vec3(1.485000,0.561000,0.369000);
+positions[78] = Vec3(0.350000,1.661000,0.928000);
+positions[79] = Vec3(1.653000,1.223000,0.578000);
+positions[80] = Vec3(0.648000,1.349000,0.253000);
+positions[81] = Vec3(0.340000,1.820000,0.483000);
+positions[82] = Vec3(2.926000,0.119000,1.421000);
+positions[83] = Vec3(1.512000,1.084000,0.156000);
+positions[84] = Vec3(1.600000,2.115000,1.792000);
+positions[85] = Vec3(1.089000,0.934000,1.584000);
+positions[86] = Vec3(1.276000,1.104000,1.230000);
+positions[87] = Vec3(0.485000,0.305000,0.428000);
+positions[88] = Vec3(1.317000,1.261000,1.795000);
+positions[89] = Vec3(0.039000,1.413000,1.085000);
+positions[90] = Vec3(0.453000,0.701000,0.605000);
+positions[91] = Vec3(1.283000,1.937000,0.752000);
+positions[92] = Vec3(0.212000,1.416000,1.447000);
+positions[93] = Vec3(0.203000,0.358000,0.723000);
+positions[94] = Vec3(0.556000,0.445000,1.364000);
+positions[95] = Vec3(1.436000,0.861000,0.911000);
+positions[96] = Vec3(0.358000,0.966000,0.176000);
+positions[97] = Vec3(1.478000,2.715000,0.427000);
+positions[98] = Vec3(1.581000,0.575000,0.809000);
+positions[99] = Vec3(1.007000,2.153000,2.887000);
+positions[100] = Vec3(2.343000,0.663000,2.513000);
+positions[101] = Vec3(2.105000,0.649000,1.635000);
+positions[102] = Vec3(0.875000,0.743000,2.459000);
+positions[103] = Vec3(0.229000,1.315000,1.879000);
+positions[104] = Vec3(0.285000,0.935000,1.700000);
+positions[105] = Vec3(2.269000,1.284000,2.234000);
+positions[106] = Vec3(1.406000,1.149000,2.767000);
+positions[107] = Vec3(1.076000,0.220000,1.849000);
+positions[108] = Vec3(2.001000,1.532000,2.881000);
+positions[109] = Vec3(2.893000,0.485000,1.860000);
+positions[110] = Vec3(1.621000,1.786000,2.624000);
+positions[111] = Vec3(0.500000,0.616000,1.818000);
+positions[112] = Vec3(0.938000,2.978000,2.104000);
+positions[113] = Vec3(0.550000,2.081000,0.454000);
+positions[114] = Vec3(1.121000,0.685000,2.196000);
+positions[115] = Vec3(1.088000,1.385000,2.184000);
+positions[116] = Vec3(1.122000,2.705000,2.080000);
+positions[117] = Vec3(0.918000,1.767000,2.861000);
+positions[118] = Vec3(2.748000,0.232000,2.126000);
+positions[119] = Vec3(1.238000,2.766000,0.109000);
+positions[120] = Vec3(1.380000,0.785000,1.961000);
+positions[121] = Vec3(1.236000,1.757000,0.150000);
+positions[122] = Vec3(1.339000,2.187000,2.592000);
+positions[123] = Vec3(1.414000,0.342000,2.714000);
+positions[124] = Vec3(1.310000,0.770000,2.589000);
+positions[125] = Vec3(1.686000,0.765000,2.321000);
+positions[126] = Vec3(1.659000,1.367000,2.780000);
+positions[127] = Vec3(0.141000,0.095000,1.903000);
+positions[128] = Vec3(2.084000,1.002000,2.520000);
+positions[129] = Vec3(2.819000,1.286000,2.626000);
+positions[130] = Vec3(1.257000,1.044000,2.401000);
+positions[131] = Vec3(1.064000,0.546000,2.839000);
+positions[132] = Vec3(0.078000,1.246000,0.010000);
+positions[133] = Vec3(1.506000,0.420000,2.223000);
+positions[134] = Vec3(1.778000,0.699000,1.920000);
+positions[135] = Vec3(1.315000,1.721000,2.733000);
+positions[136] = Vec3(0.114000,0.281000,0.279000);
+positions[137] = Vec3(1.082000,1.421000,2.596000);
+positions[138] = Vec3(3.001000,0.592000,2.276000);
+positions[139] = Vec3(1.384000,2.227000,2.992000);
+positions[140] = Vec3(1.353000,0.044000,1.985000);
+positions[141] = Vec3(1.367000,1.832000,2.383000);
+positions[142] = Vec3(0.853000,1.119000,2.230000);
+positions[143] = Vec3(1.675000,1.482000,2.295000);
+positions[144] = Vec3(1.334000,1.890000,1.904000);
+positions[145] = Vec3(1.630000,0.140000,2.939000);
+positions[146] = Vec3(0.195000,1.290000,2.300000);
+positions[147] = Vec3(2.178000,1.173000,3.001000);
+positions[148] = Vec3(0.637000,0.655000,2.126000);
+positions[149] = Vec3(0.993000,1.796000,2.529000);
+positions[150] = Vec3(0.910000,0.701000,1.845000);
+positions[151] = Vec3(0.191000,2.128000,0.355000);
+positions[152] = Vec3(0.692000,1.163000,2.578000);
+positions[153] = Vec3(1.154000,1.052000,1.974000);
+positions[154] = Vec3(1.682000,0.335000,2.509000);
+positions[155] = Vec3(0.569000,1.032000,1.895000);
+positions[156] = Vec3(1.800000,2.796000,1.295000);
+positions[157] = Vec3(2.517000,2.347000,2.878000);
+positions[158] = Vec3(0.639000,2.470000,1.678000);
+positions[159] = Vec3(0.634000,2.006000,1.829000);
+positions[160] = Vec3(0.892000,0.215000,0.566000);
+positions[161] = Vec3(1.800000,2.143000,1.491000);
+positions[162] = Vec3(1.898000,0.226000,2.765000);
+positions[163] = Vec3(0.791000,1.738000,0.260000);
+positions[164] = Vec3(0.437000,1.740000,2.048000);
+positions[165] = Vec3(1.687000,2.438000,1.166000);
+positions[166] = Vec3(1.337000,2.304000,1.643000);
+positions[167] = Vec3(1.270000,2.397000,1.033000);
+positions[168] = Vec3(0.702000,2.429000,0.591000);
+positions[169] = Vec3(0.842000,1.976000,0.724000);
+positions[170] = Vec3(1.965000,0.095000,1.206000);
+positions[171] = Vec3(0.355000,2.710000,0.618000);
+positions[172] = Vec3(0.745000,1.434000,2.781000);
+positions[173] = Vec3(0.707000,2.171000,1.502000);
+positions[174] = Vec3(1.294000,2.696000,0.847000);
+positions[175] = Vec3(1.143000,2.075000,0.276000);
+positions[176] = Vec3(1.111000,2.474000,0.312000);
+positions[177] = Vec3(1.144000,2.316000,0.633000);
+positions[178] = Vec3(1.335000,0.292000,0.515000);
+positions[179] = Vec3(1.926000,2.813000,2.703000);
+positions[180] = Vec3(0.559000,2.314000,2.904000);
+positions[181] = Vec3(1.308000,1.605000,1.534000);
+positions[182] = Vec3(0.773000,2.913000,1.217000);
+positions[183] = Vec3(1.612000,0.082000,1.027000);
+positions[184] = Vec3(1.510000,0.287000,1.787000);
+positions[185] = Vec3(0.716000,1.424000,1.843000);
+positions[186] = Vec3(1.267000,0.685000,1.160000);
+positions[187] = Vec3(0.306000,1.653000,1.717000);
+positions[188] = Vec3(0.349000,0.020000,1.275000);
+positions[189] = Vec3(0.166000,1.979000,0.804000);
+positions[190] = Vec3(1.523000,2.992000,0.711000);
+positions[191] = Vec3(1.998000,2.146000,0.088000);
+positions[192] = Vec3(0.047000,2.513000,0.642000);
+positions[193] = Vec3(0.501000,1.793000,1.438000);
+positions[194] = Vec3(1.099000,2.010000,1.626000);
+positions[195] = Vec3(2.580000,2.854000,1.328000);
+positions[196] = Vec3(1.080000,2.779000,1.190000);
+positions[197] = Vec3(0.901000,2.561000,0.948000);
+positions[198] = Vec3(0.920000,2.990000,0.844000);
+positions[199] = Vec3(0.819000,2.924000,1.711000);
+positions[200] = Vec3(0.434000,1.516000,0.063000);
+positions[201] = Vec3(1.470000,0.058000,0.231000);
+positions[202] = Vec3(0.530000,3.005000,1.550000);
+positions[203] = Vec3(0.447000,2.330000,1.277000);
+positions[204] = Vec3(1.632000,2.683000,1.593000);
+positions[205] = Vec3(0.885000,1.835000,2.072000);
+positions[206] = Vec3(0.868000,2.601000,1.425000);
+positions[207] = Vec3(0.720000,2.242000,0.907000);
+positions[208] = Vec3(1.194000,0.144000,1.065000);
+positions[209] = Vec3(0.448000,2.485000,0.959000);
+positions[210] = Vec3(1.377000,2.694000,1.352000);
+positions[211] = Vec3(1.305000,2.928000,2.713000);
+positions[212] = Vec3(1.784000,2.456000,1.981000);
+positions[213] = Vec3(0.354000,2.136000,1.563000);
+positions[214] = Vec3(0.489000,2.000000,1.108000);
+positions[215] = Vec3(1.884000,2.221000,0.461000);
+positions[216] = Vec3(1.860000,2.540000,0.306000);
+positions[217] = Vec3(1.753000,2.335000,2.768000);
+positions[218] = Vec3(1.536000,2.441000,2.344000);
+positions[219] = Vec3(0.531000,0.025000,2.235000);
+positions[220] = Vec3(0.809000,0.011000,2.834000);
+positions[221] = Vec3(0.289000,2.614000,2.879000);
+positions[222] = Vec3(0.613000,1.891000,2.337000);
+positions[223] = Vec3(0.507000,0.037000,2.694000);
+positions[224] = Vec3(0.882000,2.185000,2.583000);
+positions[225] = Vec3(0.503000,2.051000,2.615000);
+positions[226] = Vec3(1.907000,1.956000,2.831000);
+positions[227] = Vec3(2.833000,2.769000,1.644000);
+positions[228] = Vec3(1.141000,0.113000,2.945000);
+positions[229] = Vec3(0.600000,1.338000,2.200000);
+positions[230] = Vec3(0.904000,2.360000,1.952000);
+positions[231] = Vec3(0.738000,1.568000,2.437000);
+positions[232] = Vec3(1.136000,2.535000,2.805000);
+positions[233] = Vec3(1.430000,2.767000,2.321000);
+positions[234] = Vec3(1.078000,2.470000,2.385000);
+positions[235] = Vec3(0.296000,2.376000,2.560000);
+positions[236] = Vec3(0.719000,0.300000,0.075000);
+positions[237] = Vec3(0.518000,1.911000,0.080000);
+positions[238] = Vec3(0.381000,1.570000,2.450000);
+positions[239] = Vec3(0.716000,2.581000,2.697000);
+positions[240] = Vec3(1.473000,2.617000,1.936000);
+positions[241] = Vec3(0.421000,2.449000,0.229000);
+positions[242] = Vec3(0.425000,2.817000,1.910000);
+positions[243] = Vec3(1.312000,2.294000,2.057000);
+positions[244] = Vec3(1.239000,0.007000,1.539000);
+positions[245] = Vec3(0.822000,0.379000,2.086000);
+positions[246] = Vec3(0.560000,2.562000,2.227000);
+positions[247] = Vec3(0.863000,2.417000,0.050000);
+positions[248] = Vec3(1.263000,0.151000,2.332000);
+positions[249] = Vec3(0.237000,0.208000,2.336000);
+positions[250] = Vec3(0.437000,2.370000,1.910000);
+positions[251] = Vec3(1.119000,2.058000,2.207000);
+positions[252] = Vec3(1.960000,1.749000,0.118000);
+positions[253] = Vec3(2.415000,0.870000,2.757000);
+positions[254] = Vec3(1.781000,0.342000,0.366000);
+positions[255] = Vec3(2.172000,1.279000,1.421000);
+positions[256] = Vec3(1.986000,0.715000,1.301000);
+positions[257] = Vec3(1.657000,1.804000,0.810000);
+positions[258] = Vec3(2.405000,1.202000,0.416000);
+positions[259] = Vec3(1.932000,1.457000,0.786000);
+positions[260] = Vec3(1.901000,1.271000,1.207000);
+positions[261] = Vec3(1.864000,0.301000,0.810000);
+positions[262] = Vec3(1.658000,0.673000,1.558000);
+positions[263] = Vec3(2.637000,2.247000,0.396000);
+positions[264] = Vec3(1.353000,0.369000,1.438000);
+positions[265] = Vec3(0.530000,2.688000,1.346000);
+positions[266] = Vec3(0.237000,0.485000,1.047000);
+positions[267] = Vec3(2.806000,0.601000,0.822000);
+positions[268] = Vec3(1.617000,2.018000,2.136000);
+positions[269] = Vec3(2.000000,2.898000,0.022000);
+positions[270] = Vec3(2.049000,1.883000,1.001000);
+positions[271] = Vec3(2.477000,0.355000,1.786000);
+positions[272] = Vec3(1.646000,0.983000,1.266000);
+positions[273] = Vec3(1.683000,2.097000,1.114000);
+positions[274] = Vec3(2.161000,0.921000,1.065000);
+positions[275] = Vec3(2.099000,0.463000,1.942000);
+positions[276] = Vec3(2.561000,1.638000,0.572000);
+positions[277] = Vec3(2.205000,0.395000,1.005000);
+positions[278] = Vec3(2.836000,0.203000,0.698000);
+positions[279] = Vec3(2.662000,0.909000,0.966000);
+positions[280] = Vec3(0.334000,0.350000,2.767000);
+positions[281] = Vec3(2.241000,2.934000,1.248000);
+positions[282] = Vec3(2.599000,2.953000,0.921000);
+positions[283] = Vec3(2.219000,0.262000,0.058000);
+positions[284] = Vec3(0.274000,0.656000,1.456000);
+positions[285] = Vec3(1.814000,1.008000,0.882000);
+positions[286] = Vec3(2.793000,1.395000,0.316000);
+positions[287] = Vec3(0.773000,1.753000,1.639000);
+positions[288] = Vec3(2.321000,0.994000,1.591000);
+positions[289] = Vec3(2.243000,2.255000,1.690000);
+positions[290] = Vec3(0.178000,1.342000,0.327000);
+positions[291] = Vec3(1.623000,1.756000,1.426000);
+positions[292] = Vec3(2.252000,0.109000,0.375000);
+positions[293] = Vec3(3.003000,1.895000,1.895000);
+positions[294] = Vec3(0.407000,0.831000,2.756000);
+positions[295] = Vec3(2.193000,0.956000,0.632000);
+positions[296] = Vec3(2.405000,0.641000,1.107000);
+positions[297] = Vec3(2.361000,0.958000,0.162000);
+positions[298] = Vec3(2.173000,1.544000,0.528000);
+positions[299] = Vec3(1.565000,1.380000,1.428000);
+positions[300] = Vec3(2.342000,0.538000,0.253000);
+positions[301] = Vec3(1.910000,0.701000,0.954000);
+positions[302] = Vec3(2.910000,0.288000,2.938000);
+positions[303] = Vec3(0.257000,1.189000,0.958000);
+positions[304] = Vec3(0.134000,1.775000,1.243000);
+positions[305] = Vec3(2.476000,1.583000,1.956000);
+positions[306] = Vec3(1.838000,1.791000,2.354000);
+positions[307] = Vec3(1.906000,1.338000,1.696000);
+positions[308] = Vec3(2.413000,2.869000,0.166000);
+positions[309] = Vec3(3.006000,1.038000,1.322000);
+positions[310] = Vec3(1.961000,0.962000,1.605000);
+positions[311] = Vec3(0.082000,2.857000,0.020000);
+positions[312] = Vec3(2.408000,1.499000,0.062000);
+positions[313] = Vec3(2.349000,0.267000,1.415000);
+positions[314] = Vec3(2.327000,1.717000,2.350000);
+positions[315] = Vec3(2.928000,0.810000,1.582000);
+positions[316] = Vec3(2.150000,0.336000,0.576000);
+positions[317] = Vec3(2.664000,1.085000,2.962000);
+positions[318] = Vec3(2.851000,0.670000,1.174000);
+positions[319] = Vec3(1.954000,1.013000,1.975000);
+positions[320] = Vec3(2.474000,1.542000,1.545000);
+positions[321] = Vec3(2.826000,0.455000,1.490000);
+positions[322] = Vec3(2.140000,2.826000,0.558000);
+positions[323] = Vec3(2.151000,1.684000,1.780000);
+positions[324] = Vec3(0.174000,0.673000,0.397000);
+positions[325] = Vec3(0.066000,1.708000,0.160000);
+positions[326] = Vec3(2.158000,0.303000,2.582000);
+positions[327] = Vec3(2.602000,1.611000,2.632000);
+positions[328] = Vec3(2.566000,1.138000,2.465000);
+positions[329] = Vec3(2.026000,1.443000,2.477000);
+positions[330] = Vec3(2.365000,0.309000,2.255000);
+positions[331] = Vec3(1.636000,1.107000,2.058000);
+positions[332] = Vec3(2.522000,2.584000,2.399000);
+positions[333] = Vec3(2.537000,2.900000,2.158000);
+positions[334] = Vec3(2.660000,0.537000,2.577000);
+positions[335] = Vec3(2.679000,1.158000,1.724000);
+positions[336] = Vec3(0.220000,1.894000,2.498000);
+positions[337] = Vec3(2.266000,1.248000,1.837000);
+positions[338] = Vec3(0.055000,1.656000,2.128000);
+positions[339] = Vec3(2.899000,1.902000,2.823000);
+positions[340] = Vec3(0.085000,2.994000,2.720000);
+positions[341] = Vec3(0.013000,0.889000,2.468000);
+positions[342] = Vec3(1.804000,0.372000,1.636000);
+positions[343] = Vec3(0.201000,1.616000,2.824000);
+positions[344] = Vec3(0.369000,1.273000,2.699000);
+positions[345] = Vec3(2.996000,0.355000,2.596000);
+positions[346] = Vec3(2.867000,1.314000,2.107000);
+positions[347] = Vec3(2.611000,0.563000,2.140000);
+positions[348] = Vec3(2.676000,2.954000,2.955000);
+positions[349] = Vec3(0.256000,0.848000,2.062000);
+positions[350] = Vec3(2.530000,0.028000,2.528000);
+positions[351] = Vec3(0.537000,1.273000,1.596000);
+positions[352] = Vec3(0.004000,1.004000,0.401000);
+positions[353] = Vec3(1.676000,1.060000,2.463000);
+positions[354] = Vec3(2.622000,1.473000,2.257000);
+positions[355] = Vec3(2.917000,2.991000,2.316000);
+positions[356] = Vec3(0.672000,1.123000,2.984000);
+positions[357] = Vec3(2.229000,1.806000,2.673000);
+positions[358] = Vec3(0.463000,0.951000,2.383000);
+positions[359] = Vec3(2.126000,0.049000,2.037000);
+positions[360] = Vec3(2.868000,0.876000,2.015000);
+positions[361] = Vec3(2.720000,2.582000,0.079000);
+positions[362] = Vec3(1.966000,0.693000,2.624000);
+positions[363] = Vec3(1.971000,0.398000,2.318000);
+positions[364] = Vec3(0.337000,0.630000,2.458000);
+positions[365] = Vec3(2.562000,1.044000,2.040000);
+positions[366] = Vec3(2.817000,1.485000,2.963000);
+positions[367] = Vec3(2.514000,0.621000,2.992000);
+positions[368] = Vec3(3.000000,1.551000,2.496000);
+positions[369] = Vec3(0.698000,2.167000,2.180000);
+positions[370] = Vec3(2.693000,0.849000,2.389000);
+positions[371] = Vec3(2.092000,2.565000,2.986000);
+positions[372] = Vec3(2.010000,3.001000,0.819000);
+positions[373] = Vec3(2.392000,2.622000,1.636000);
+positions[374] = Vec3(2.086000,2.325000,1.340000);
+positions[375] = Vec3(2.578000,2.971000,0.502000);
+positions[376] = Vec3(1.871000,2.789000,2.225000);
+positions[377] = Vec3(2.230000,2.985000,1.594000);
+positions[378] = Vec3(2.860000,2.788000,0.729000);
+positions[379] = Vec3(2.051000,1.928000,1.472000);
+positions[380] = Vec3(2.307000,2.219000,1.067000);
+positions[381] = Vec3(2.369000,2.572000,1.289000);
+positions[382] = Vec3(2.206000,1.924000,0.693000);
+positions[383] = Vec3(1.984000,2.058000,2.005000);
+positions[384] = Vec3(2.287000,1.854000,0.317000);
+positions[385] = Vec3(2.525000,0.345000,0.686000);
+positions[386] = Vec3(2.933000,1.920000,1.053000);
+positions[387] = Vec3(0.324000,2.324000,0.601000);
+positions[388] = Vec3(2.042000,1.576000,1.277000);
+positions[389] = Vec3(0.031000,2.376000,0.949000);
+positions[390] = Vec3(2.519000,2.250000,1.465000);
+positions[391] = Vec3(0.221000,2.722000,1.652000);
+positions[392] = Vec3(2.409000,2.361000,2.051000);
+positions[393] = Vec3(2.472000,1.917000,1.673000);
+positions[394] = Vec3(0.999000,2.715000,0.562000);
+positions[395] = Vec3(1.669000,0.017000,1.508000);
+positions[396] = Vec3(1.924000,1.777000,0.542000);
+positions[397] = Vec3(2.635000,2.634000,1.905000);
+positions[398] = Vec3(2.042000,2.628000,1.025000);
+positions[399] = Vec3(2.694000,1.974000,2.009000);
+positions[400] = Vec3(2.988000,2.221000,1.333000);
+positions[401] = Vec3(1.772000,0.196000,1.978000);
+positions[402] = Vec3(2.893000,2.961000,0.283000);
+positions[403] = Vec3(2.615000,0.261000,0.245000);
+positions[404] = Vec3(2.797000,2.521000,1.412000);
+positions[405] = Vec3(0.013000,2.497000,0.246000);
+positions[406] = Vec3(1.875000,2.861000,1.801000);
+positions[407] = Vec3(2.800000,2.617000,1.049000);
+positions[408] = Vec3(2.824000,1.858000,1.487000);
+positions[409] = Vec3(2.434000,1.868000,1.275000);
+positions[410] = Vec3(2.814000,0.526000,0.384000);
+positions[411] = Vec3(2.844000,2.545000,2.246000);
+positions[412] = Vec3(1.896000,2.587000,0.719000);
+positions[413] = Vec3(0.350000,0.055000,0.076000);
+positions[414] = Vec3(2.686000,1.784000,0.222000);
+positions[415] = Vec3(2.724000,1.604000,0.989000);
+positions[416] = Vec3(0.807000,1.761000,1.122000);
+positions[417] = Vec3(2.120000,2.382000,2.226000);
+positions[418] = Vec3(2.058000,1.587000,2.067000);
+positions[419] = Vec3(2.904000,2.571000,2.686000);
+positions[420] = Vec3(2.228000,2.910000,2.410000);
+positions[421] = Vec3(2.797000,2.142000,0.114000);
+positions[422] = Vec3(2.905000,1.875000,0.480000);
+positions[423] = Vec3(1.881000,2.565000,2.469000);
+positions[424] = Vec3(2.404000,1.929000,2.999000);
+positions[425] = Vec3(2.389000,2.814000,2.782000);
+positions[426] = Vec3(2.520000,0.301000,2.815000);
+positions[427] = Vec3(2.726000,1.907000,2.339000);
+positions[428] = Vec3(2.880000,2.273000,2.500000);
+positions[429] = Vec3(2.574000,2.045000,2.716000);
+positions[430] = Vec3(2.988000,2.288000,2.001000);
+positions[431] = Vec3(0.011000,2.341000,2.904000);
+positions[432] = Vec3(0.215000,2.265000,2.257000);
+positions[433] = Vec3(2.268000,2.311000,0.234000);
+positions[434] = Vec3(2.462000,2.621000,0.550000);
+positions[435] = Vec3(1.530000,2.540000,2.728000);
+positions[436] = Vec3(2.162000,2.306000,2.687000);
+positions[437] = Vec3(2.748000,2.301000,1.734000);
+positions[438] = Vec3(2.334000,1.976000,2.041000);
+positions[439] = Vec3(1.981000,2.076000,2.443000);
+positions[440] = Vec3(2.301000,1.367000,2.665000);
+positions[441] = Vec3(2.399000,2.164000,2.403000);
+positions[442] = Vec3(0.244000,2.713000,2.257000);
+positions[443] = Vec3(0.683000,0.488000,2.781000);
+positions[444] = Vec3(2.194000,2.711000,1.993000);
+positions[445] = Vec3(2.947000,2.848000,2.001000);
+positions[446] = Vec3(0.223000,1.981000,2.913000);
+positions[447] = Vec3(0.010000,1.226000,0.917000);
+positions[448] = Vec3(1.911000,0.426000,0.582000);
+positions[449] = Vec3(2.204000,0.015000,0.136000);
+positions[450] = Vec3(0.927000,0.138000,1.645000);
+positions[451] = Vec3(0.155000,0.885000,1.479000);
+positions[452] = Vec3(1.550000,1.933000,1.261000);
+positions[453] = Vec3(1.304000,0.407000,0.287000);
+positions[454] = Vec3(0.270000,1.384000,2.910000);
+positions[455] = Vec3(0.516000,1.817000,1.695000);
+positions[456] = Vec3(1.458000,2.879000,1.523000);
+positions[457] = Vec3(1.702000,1.670000,0.593000);
+positions[458] = Vec3(1.974000,1.380000,0.534000);
+positions[459] = Vec3(2.835000,1.185000,0.479000);
+positions[460] = Vec3(0.548000,2.979000,1.126000);
+positions[461] = Vec3(1.202000,2.174000,1.466000);
+positions[462] = Vec3(1.237000,1.701000,0.653000);
+positions[463] = Vec3(2.939000,0.761000,0.349000);
+positions[464] = Vec3(1.667000,2.119000,0.377000);
+positions[465] = Vec3(1.196000,0.552000,1.372000);
+positions[466] = Vec3(1.416000,0.901000,1.178000);
+positions[467] = Vec3(0.519000,1.577000,2.227000);
+positions[468] = Vec3(1.214000,1.281000,1.063000);
+positions[469] = Vec3(0.822000,0.433000,1.375000);
+positions[470] = Vec3(0.095000,2.760000,0.604000);
+positions[471] = Vec3(1.325000,2.144000,1.848000);
+positions[472] = Vec3(0.681000,0.896000,1.285000);
+positions[473] = Vec3(0.406000,2.936000,0.717000);
+positions[474] = Vec3(0.565000,1.852000,0.349000);
+positions[475] = Vec3(0.597000,1.651000,1.020000);
+positions[476] = Vec3(1.236000,0.170000,1.335000);
+positions[477] = Vec3(0.586000,0.441000,1.980000);
+positions[478] = Vec3(1.443000,1.208000,1.575000);
+positions[479] = Vec3(0.247000,0.243000,0.502000);
+positions[480] = Vec3(1.386000,1.564000,0.236000);
+positions[481] = Vec3(0.871000,1.063000,0.930000);
+positions[482] = Vec3(0.136000,0.992000,0.621000);
+positions[483] = Vec3(0.889000,0.986000,0.010000);
+positions[484] = Vec3(1.107000,2.731000,1.452000);
+positions[485] = Vec3(0.942000,2.471000,0.517000);
+positions[486] = Vec3(0.989000,0.652000,0.747000);
+positions[487] = Vec3(0.899000,1.235000,2.707000);
+positions[488] = Vec3(1.105000,0.684000,0.068000);
+positions[489] = Vec3(1.660000,1.235000,2.276000);
+positions[490] = Vec3(1.593000,1.883000,1.915000);
+positions[491] = Vec3(1.528000,1.613000,0.920000);
+positions[492] = Vec3(0.459000,1.046000,1.011000);
+positions[493] = Vec3(0.213000,0.612000,0.644000);
+positions[494] = Vec3(0.078000,1.392000,1.676000);
+positions[495] = Vec3(0.605000,0.491000,0.574000);
+positions[496] = Vec3(0.990000,1.586000,1.076000);
+positions[497] = Vec3(0.297000,1.434000,1.028000);
+positions[498] = Vec3(1.101000,1.471000,1.443000);
+positions[499] = Vec3(0.072000,0.139000,1.653000);
+positions[500] = Vec3(0.633000,0.884000,0.645000);
+positions[501] = Vec3(0.352000,2.841000,1.463000);
+positions[502] = Vec3(0.418000,0.774000,0.350000);
+positions[503] = Vec3(2.641000,0.198000,0.869000);
+positions[504] = Vec3(0.608000,1.341000,0.695000);
+positions[505] = Vec3(1.778000,1.151000,1.830000);
+positions[506] = Vec3(1.669000,0.342000,2.768000);
+positions[507] = Vec3(1.256000,0.994000,0.798000);
+positions[508] = Vec3(1.068000,0.375000,1.036000);
+positions[509] = Vec3(0.910000,0.758000,1.589000);
+positions[510] = Vec3(0.243000,2.452000,0.805000);
+positions[511] = Vec3(1.018000,0.764000,1.122000);
+positions[512] = Vec3(2.464000,1.089000,1.404000);
+positions[513] = Vec3(0.670000,0.564000,0.034000);
+positions[514] = Vec3(0.030000,1.296000,1.310000);
+positions[515] = Vec3(1.210000,1.785000,1.691000);
+positions[516] = Vec3(0.022000,0.620000,0.974000);
+positions[517] = Vec3(1.499000,1.277000,2.986000);
+positions[518] = Vec3(1.227000,1.896000,1.006000);
+positions[519] = Vec3(0.528000,1.022000,1.635000);
+positions[520] = Vec3(1.887000,2.670000,0.089000);
+positions[521] = Vec3(1.661000,0.825000,0.793000);
+positions[522] = Vec3(0.831000,1.494000,0.374000);
+positions[523] = Vec3(1.356000,0.613000,0.930000);
+positions[524] = Vec3(0.667000,0.600000,0.968000);
+positions[525] = Vec3(1.154000,1.702000,2.925000);
+positions[526] = Vec3(1.420000,1.581000,1.289000);
+positions[527] = Vec3(1.383000,0.041000,0.932000);
+positions[528] = Vec3(1.727000,0.140000,1.725000);
+positions[529] = Vec3(0.711000,1.215000,2.004000);
+positions[530] = Vec3(1.061000,1.067000,1.366000);
+positions[531] = Vec3(0.377000,0.597000,1.224000);
+positions[532] = Vec3(0.274000,0.719000,1.842000);
+positions[533] = Vec3(0.840000,1.658000,1.874000);
+positions[534] = Vec3(0.877000,0.290000,0.311000);
+positions[535] = Vec3(2.130000,1.153000,1.196000);
+positions[536] = Vec3(1.028000,1.379000,0.747000);
+positions[537] = Vec3(1.107000,2.450000,2.079000);
+positions[538] = Vec3(1.419000,1.333000,0.585000);
+positions[539] = Vec3(0.430000,1.305000,1.369000);
+positions[540] = Vec3(0.775000,1.363000,1.596000);
+positions[541] = Vec3(1.522000,2.009000,0.736000);
+positions[542] = Vec3(0.857000,1.722000,0.696000);
+positions[543] = Vec3(0.722000,2.831000,1.478000);
+positions[544] = Vec3(0.411000,1.673000,0.681000);
+positions[545] = Vec3(1.511000,0.456000,0.597000);
+positions[546] = Vec3(2.684000,0.820000,2.996000);
+positions[547] = Vec3(1.593000,1.713000,2.369000);
+positions[548] = Vec3(1.113000,0.803000,1.958000);
+positions[549] = Vec3(1.267000,1.095000,0.254000);
+positions[550] = Vec3(2.120000,0.540000,2.477000);
+positions[551] = Vec3(0.566000,1.409000,2.588000);
+positions[552] = Vec3(0.261000,0.872000,2.546000);
+positions[553] = Vec3(1.878000,1.446000,2.680000);
+positions[554] = Vec3(0.878000,1.606000,2.658000);
+positions[555] = Vec3(1.564000,0.749000,1.786000);
+positions[556] = Vec3(1.412000,1.942000,2.625000);
+positions[557] = Vec3(1.660000,1.114000,2.710000);
+positions[558] = Vec3(1.118000,0.813000,2.424000);
+positions[559] = Vec3(1.482000,0.893000,2.434000);
+positions[560] = Vec3(1.093000,1.129000,1.740000);
+positions[561] = Vec3(2.163000,0.849000,2.709000);
+positions[562] = Vec3(1.201000,1.429000,1.957000);
+positions[563] = Vec3(0.235000,2.258000,2.002000);
+positions[564] = Vec3(0.413000,1.444000,0.314000);
+positions[565] = Vec3(0.164000,0.450000,2.408000);
+positions[566] = Vec3(1.551000,0.851000,0.033000);
+positions[567] = Vec3(0.659000,0.228000,2.807000);
+positions[568] = Vec3(0.741000,0.131000,2.124000);
+positions[569] = Vec3(0.455000,0.567000,2.682000);
+positions[570] = Vec3(0.729000,0.971000,2.408000);
+positions[571] = Vec3(1.487000,2.820000,0.162000);
+positions[572] = Vec3(1.855000,0.700000,2.858000);
+positions[573] = Vec3(0.305000,1.074000,1.926000);
+positions[574] = Vec3(1.300000,0.153000,1.747000);
+positions[575] = Vec3(1.272000,1.249000,2.568000);
+positions[576] = Vec3(0.431000,0.743000,2.238000);
+positions[577] = Vec3(0.493000,0.240000,0.184000);
+positions[578] = Vec3(1.734000,0.506000,2.317000);
+positions[579] = Vec3(0.874000,0.631000,2.692000);
+positions[580] = Vec3(0.473000,2.790000,2.161000);
+positions[581] = Vec3(1.310000,0.571000,2.759000);
+positions[582] = Vec3(0.677000,0.798000,1.916000);
+positions[583] = Vec3(0.944000,0.442000,1.858000);
+positions[584] = Vec3(3.006000,2.098000,2.976000);
+positions[585] = Vec3(0.864000,0.592000,2.231000);
+positions[586] = Vec3(1.366000,0.611000,2.147000);
+positions[587] = Vec3(2.871000,1.217000,2.880000);
+positions[588] = Vec3(1.674000,2.664000,2.336000);
+positions[589] = Vec3(1.757000,0.879000,2.101000);
+positions[590] = Vec3(1.293000,2.939000,2.457000);
+positions[591] = Vec3(1.108000,1.131000,2.206000);
+positions[592] = Vec3(1.207000,1.658000,2.498000);
+positions[593] = Vec3(0.850000,1.373000,2.312000);
+positions[594] = Vec3(1.413000,1.060000,1.939000);
+positions[595] = Vec3(1.138000,0.140000,2.102000);
+positions[596] = Vec3(0.752000,1.307000,1.190000);
+positions[597] = Vec3(1.254000,0.942000,2.790000);
+positions[598] = Vec3(1.544000,1.614000,2.800000);
+positions[599] = Vec3(2.128000,0.302000,2.833000);
+positions[600] = Vec3(0.300000,1.744000,0.027000);
+positions[601] = Vec3(1.878000,2.986000,2.060000);
+positions[602] = Vec3(1.528000,0.233000,2.045000);
+positions[603] = Vec3(1.146000,1.817000,2.067000);
+positions[604] = Vec3(1.037000,2.746000,0.813000);
+positions[605] = Vec3(0.524000,0.610000,1.566000);
+positions[606] = Vec3(0.945000,2.964000,0.503000);
+positions[607] = Vec3(1.788000,2.565000,0.965000);
+positions[608] = Vec3(0.471000,2.510000,0.491000);
+positions[609] = Vec3(0.512000,2.043000,1.371000);
+positions[610] = Vec3(2.316000,2.423000,1.494000);
+positions[611] = Vec3(1.575000,2.394000,2.953000);
+positions[612] = Vec3(2.845000,2.869000,0.985000);
+positions[613] = Vec3(1.016000,2.335000,1.003000);
+positions[614] = Vec3(0.998000,2.830000,1.879000);
+positions[615] = Vec3(0.624000,2.508000,0.075000);
+positions[616] = Vec3(1.362000,2.808000,2.069000);
+positions[617] = Vec3(1.747000,0.068000,0.810000);
+positions[618] = Vec3(1.768000,2.355000,0.661000);
+positions[619] = Vec3(1.535000,2.410000,2.085000);
+positions[620] = Vec3(0.844000,2.004000,1.646000);
+positions[621] = Vec3(1.124000,0.280000,0.649000);
+positions[622] = Vec3(0.689000,2.170000,0.648000);
+positions[623] = Vec3(0.849000,2.666000,1.175000);
+positions[624] = Vec3(2.975000,1.963000,1.308000);
+positions[625] = Vec3(1.074000,2.082000,0.714000);
+positions[626] = Vec3(1.284000,2.651000,1.110000);
+positions[627] = Vec3(1.669000,0.205000,0.180000);
+positions[628] = Vec3(1.716000,0.047000,1.253000);
+positions[629] = Vec3(0.501000,2.241000,1.043000);
+positions[630] = Vec3(1.038000,1.833000,0.305000);
+positions[631] = Vec3(0.646000,2.431000,1.424000);
+positions[632] = Vec3(1.383000,2.059000,2.230000);
+positions[633] = Vec3(0.370000,2.566000,1.192000);
+positions[634] = Vec3(1.355000,2.006000,0.120000);
+positions[635] = Vec3(2.113000,0.075000,0.589000);
+positions[636] = Vec3(1.850000,0.448000,1.890000);
+positions[637] = Vec3(1.215000,2.704000,0.405000);
+positions[638] = Vec3(0.575000,2.997000,1.798000);
+positions[639] = Vec3(0.967000,2.586000,2.603000);
+positions[640] = Vec3(0.276000,1.669000,1.430000);
+positions[641] = Vec3(1.483000,2.284000,1.128000);
+positions[642] = Vec3(0.983000,3.003000,1.099000);
+positions[643] = Vec3(0.539000,2.222000,1.720000);
+positions[644] = Vec3(0.648000,2.826000,2.751000);
+positions[645] = Vec3(0.803000,1.994000,0.993000);
+positions[646] = Vec3(0.451000,0.216000,1.438000);
+positions[647] = Vec3(1.604000,2.512000,0.334000);
+positions[648] = Vec3(1.980000,2.022000,0.588000);
+positions[649] = Vec3(1.843000,2.834000,1.544000);
+positions[650] = Vec3(1.835000,3.005000,2.858000);
+positions[651] = Vec3(0.679000,2.499000,0.838000);
+positions[652] = Vec3(0.012000,2.637000,1.524000);
+positions[653] = Vec3(0.314000,2.065000,0.602000);
+positions[654] = Vec3(1.157000,0.004000,0.173000);
+positions[655] = Vec3(0.736000,1.705000,1.382000);
+positions[656] = Vec3(1.511000,2.736000,0.690000);
+positions[657] = Vec3(1.330000,2.541000,1.735000);
+positions[658] = Vec3(0.744000,0.170000,0.785000);
+positions[659] = Vec3(2.593000,2.794000,0.703000);
+positions[660] = Vec3(0.275000,1.872000,1.043000);
+positions[661] = Vec3(1.624000,2.608000,1.341000);
+positions[662] = Vec3(1.382000,0.122000,2.855000);
+positions[663] = Vec3(1.326000,2.434000,0.783000);
+positions[664] = Vec3(0.117000,0.116000,1.254000);
+positions[665] = Vec3(1.045000,2.970000,2.748000);
+positions[666] = Vec3(1.341000,2.692000,2.799000);
+positions[667] = Vec3(1.797000,2.586000,2.709000);
+positions[668] = Vec3(0.890000,2.484000,1.716000);
+positions[669] = Vec3(2.373000,2.558000,1.889000);
+positions[670] = Vec3(1.566000,2.323000,2.574000);
+positions[671] = Vec3(1.257000,2.280000,0.399000);
+positions[672] = Vec3(0.679000,2.130000,2.434000);
+positions[673] = Vec3(2.016000,2.334000,2.462000);
+positions[674] = Vec3(1.077000,2.213000,2.416000);
+positions[675] = Vec3(0.581000,1.950000,2.081000);
+positions[676] = Vec3(0.805000,2.315000,2.810000);
+positions[677] = Vec3(0.844000,1.787000,2.322000);
+positions[678] = Vec3(0.980000,2.205000,0.129000);
+positions[679] = Vec3(2.468000,0.603000,2.740000);
+positions[680] = Vec3(2.366000,2.403000,2.299000);
+positions[681] = Vec3(0.337000,2.487000,2.329000);
+positions[682] = Vec3(2.007000,2.793000,2.452000);
+positions[683] = Vec3(1.072000,2.571000,0.063000);
+positions[684] = Vec3(1.217000,2.283000,2.806000);
+positions[685] = Vec3(0.459000,2.477000,2.728000);
+positions[686] = Vec3(0.958000,1.975000,2.710000);
+positions[687] = Vec3(0.914000,2.111000,2.052000);
+positions[688] = Vec3(0.768000,2.958000,0.075000);
+positions[689] = Vec3(0.474000,1.805000,2.533000);
+positions[690] = Vec3(1.313000,2.552000,2.395000);
+positions[691] = Vec3(1.853000,2.014000,2.229000);
+positions[692] = Vec3(2.405000,2.230000,2.658000);
+positions[693] = Vec3(0.727000,1.781000,0.016000);
+positions[694] = Vec3(0.974000,2.791000,2.271000);
+positions[695] = Vec3(0.438000,0.096000,2.457000);
+positions[696] = Vec3(0.652000,2.392000,2.064000);
+positions[697] = Vec3(1.972000,2.209000,2.834000);
+positions[698] = Vec3(0.333000,0.141000,2.088000);
+positions[699] = Vec3(1.813000,1.952000,0.063000);
+positions[700] = Vec3(0.166000,2.838000,1.877000);
+positions[701] = Vec3(1.772000,0.487000,0.951000);
+positions[702] = Vec3(1.924000,1.404000,1.434000);
+positions[703] = Vec3(2.734000,0.348000,1.712000);
+positions[704] = Vec3(2.874000,0.729000,1.811000);
+positions[705] = Vec3(1.841000,0.877000,1.137000);
+positions[706] = Vec3(2.327000,1.491000,1.768000);
+positions[707] = Vec3(1.916000,1.483000,1.057000);
+positions[708] = Vec3(2.783000,0.850000,0.745000);
+positions[709] = Vec3(1.829000,1.526000,0.085000);
+positions[710] = Vec3(2.426000,1.082000,0.652000);
+positions[711] = Vec3(1.645000,1.241000,1.217000);
+positions[712] = Vec3(2.286000,0.725000,0.084000);
+positions[713] = Vec3(2.755000,0.691000,1.421000);
+positions[714] = Vec3(2.651000,0.591000,1.023000);
+positions[715] = Vec3(2.040000,0.863000,0.442000);
+positions[716] = Vec3(0.035000,0.109000,2.497000);
+positions[717] = Vec3(0.127000,1.410000,0.572000);
+positions[718] = Vec3(2.174000,0.357000,0.307000);
+positions[719] = Vec3(2.705000,1.508000,0.758000);
+positions[720] = Vec3(2.223000,1.407000,2.913000);
+positions[721] = Vec3(2.528000,1.722000,1.088000);
+positions[722] = Vec3(2.860000,0.345000,0.198000);
+positions[723] = Vec3(2.580000,1.789000,1.479000);
+positions[724] = Vec3(2.779000,0.295000,1.295000);
+positions[725] = Vec3(0.097000,0.434000,2.826000);
+positions[726] = Vec3(2.952000,1.654000,1.091000);
+positions[727] = Vec3(0.119000,1.878000,0.343000);
+positions[728] = Vec3(1.718000,1.173000,0.327000);
+positions[729] = Vec3(2.833000,0.016000,0.527000);
+positions[730] = Vec3(2.085000,1.779000,2.888000);
+positions[731] = Vec3(2.754000,2.952000,1.485000);
+positions[732] = Vec3(2.826000,0.935000,1.162000);
+positions[733] = Vec3(1.564000,1.585000,1.615000);
+positions[734] = Vec3(2.132000,0.645000,1.093000);
+positions[735] = Vec3(2.294000,1.490000,1.350000);
+positions[736] = Vec3(0.081000,0.490000,1.479000);
+positions[737] = Vec3(2.118000,1.165000,1.642000);
+positions[738] = Vec3(2.141000,0.121000,1.390000);
+positions[739] = Vec3(2.385000,0.389000,1.196000);
+positions[740] = Vec3(0.049000,0.166000,0.817000);
+positions[741] = Vec3(1.993000,0.806000,1.814000);
+positions[742] = Vec3(0.006000,1.450000,0.171000);
+positions[743] = Vec3(2.297000,0.428000,0.764000);
+positions[744] = Vec3(2.851000,0.469000,2.114000);
+positions[745] = Vec3(1.814000,1.957000,0.945000);
+positions[746] = Vec3(0.386000,0.327000,0.902000);
+positions[747] = Vec3(2.452000,1.070000,1.807000);
+positions[748] = Vec3(2.309000,1.537000,2.159000);
+positions[749] = Vec3(2.712000,1.497000,2.007000);
+positions[750] = Vec3(1.727000,0.924000,1.503000);
+positions[751] = Vec3(0.861000,0.801000,0.344000);
+positions[752] = Vec3(1.740000,1.245000,0.819000);
+positions[753] = Vec3(0.117000,0.042000,0.197000);
+positions[754] = Vec3(2.557000,0.996000,0.317000);
+positions[755] = Vec3(2.228000,1.588000,2.548000);
+positions[756] = Vec3(2.849000,1.557000,2.708000);
+positions[757] = Vec3(0.152000,1.107000,0.219000);
+positions[758] = Vec3(2.460000,1.318000,1.002000);
+positions[759] = Vec3(2.405000,1.436000,0.528000);
+positions[760] = Vec3(2.135000,1.179000,2.046000);
+positions[761] = Vec3(1.726000,0.588000,0.286000);
+positions[762] = Vec3(2.831000,1.053000,1.538000);
+positions[763] = Vec3(1.932000,1.556000,1.833000);
+positions[764] = Vec3(2.423000,0.900000,1.064000);
+positions[765] = Vec3(3.001000,1.807000,0.709000);
+positions[766] = Vec3(0.578000,1.095000,0.223000);
+positions[767] = Vec3(2.215000,1.160000,0.252000);
+positions[768] = Vec3(2.050000,0.921000,0.835000);
+positions[769] = Vec3(2.080000,1.682000,0.738000);
+positions[770] = Vec3(2.851000,1.753000,0.027000);
+positions[771] = Vec3(0.203000,0.509000,0.202000);
+positions[772] = Vec3(1.967000,1.018000,0.018000);
+positions[773] = Vec3(1.869000,0.878000,2.472000);
+positions[774] = Vec3(1.917000,0.228000,2.507000);
+positions[775] = Vec3(0.316000,0.795000,2.991000);
+positions[776] = Vec3(2.175000,1.229000,2.472000);
+positions[777] = Vec3(2.405000,1.062000,2.931000);
+positions[778] = Vec3(2.501000,0.511000,2.369000);
+positions[779] = Vec3(2.641000,0.819000,2.141000);
+positions[780] = Vec3(0.649000,1.384000,3.006000);
+positions[781] = Vec3(1.012000,0.329000,2.963000);
+positions[782] = Vec3(2.755000,0.350000,2.718000);
+positions[783] = Vec3(2.315000,0.153000,2.454000);
+positions[784] = Vec3(2.583000,1.696000,2.389000);
+positions[785] = Vec3(0.439000,2.593000,1.776000);
+positions[786] = Vec3(2.630000,1.390000,0.116000);
+positions[787] = Vec3(2.854000,0.669000,2.478000);
+positions[788] = Vec3(2.551000,1.342000,2.621000);
+positions[789] = Vec3(2.533000,2.734000,2.987000);
+positions[790] = Vec3(2.772000,2.446000,2.875000);
+positions[791] = Vec3(2.817000,1.051000,2.498000);
+positions[792] = Vec3(2.688000,1.404000,1.621000);
+positions[793] = Vec3(0.083000,2.737000,2.775000);
+positions[794] = Vec3(2.514000,0.322000,2.041000);
+positions[795] = Vec3(2.470000,0.900000,2.504000);
+positions[796] = Vec3(2.790000,0.444000,0.624000);
+positions[797] = Vec3(0.040000,0.840000,2.202000);
+positions[798] = Vec3(0.530000,1.067000,2.764000);
+positions[799] = Vec3(0.191000,1.385000,2.541000);
+positions[800] = Vec3(2.465000,0.363000,0.051000);
+positions[801] = Vec3(1.850000,1.902000,2.592000);
+positions[802] = Vec3(1.432000,0.306000,2.449000);
+positions[803] = Vec3(2.259000,0.489000,1.753000);
+positions[804] = Vec3(2.803000,1.118000,1.956000);
+positions[805] = Vec3(2.426000,0.147000,1.636000);
+positions[806] = Vec3(2.880000,1.846000,2.133000);
+positions[807] = Vec3(2.862000,2.110000,1.867000);
+positions[808] = Vec3(0.424000,1.184000,2.299000);
+positions[809] = Vec3(2.518000,1.218000,2.228000);
+positions[810] = Vec3(2.153000,0.834000,1.468000);
+positions[811] = Vec3(0.105000,1.397000,2.088000);
+positions[812] = Vec3(2.579000,0.601000,0.316000);
+positions[813] = Vec3(2.594000,2.106000,2.968000);
+positions[814] = Vec3(0.448000,1.435000,1.783000);
+positions[815] = Vec3(2.125000,0.299000,2.132000);
+positions[816] = Vec3(2.849000,1.402000,2.356000);
+positions[817] = Vec3(2.956000,0.091000,2.078000);
+positions[818] = Vec3(0.156000,1.696000,2.357000);
+positions[819] = Vec3(1.566000,2.211000,1.557000);
+positions[820] = Vec3(2.047000,0.194000,0.985000);
+positions[821] = Vec3(1.947000,2.680000,0.488000);
+positions[822] = Vec3(2.343000,2.796000,1.447000);
+positions[823] = Vec3(2.006000,2.332000,0.265000);
+positions[824] = Vec3(2.396000,1.834000,0.546000);
+positions[825] = Vec3(2.538000,2.059000,2.207000);
+positions[826] = Vec3(0.110000,2.360000,0.447000);
+positions[827] = Vec3(2.198000,2.448000,1.136000);
+positions[828] = Vec3(2.420000,2.121000,1.271000);
+positions[829] = Vec3(0.422000,2.192000,0.260000);
+positions[830] = Vec3(2.145000,2.767000,2.839000);
+positions[831] = Vec3(2.434000,2.398000,0.421000);
+positions[832] = Vec3(2.489000,2.175000,1.718000);
+positions[833] = Vec3(2.870000,2.527000,0.814000);
+positions[834] = Vec3(2.741000,2.016000,0.337000);
+positions[835] = Vec3(1.997000,2.574000,2.107000);
+positions[836] = Vec3(0.002000,2.128000,0.932000);
+positions[837] = Vec3(2.787000,2.375000,0.234000);
+positions[838] = Vec3(2.235000,1.852000,1.620000);
+positions[839] = Vec3(2.782000,1.642000,0.422000);
+positions[840] = Vec3(2.915000,1.760000,1.699000);
+positions[841] = Vec3(2.047000,2.178000,1.549000);
+positions[842] = Vec3(1.808000,1.878000,1.556000);
+positions[843] = Vec3(2.224000,2.043000,0.913000);
+positions[844] = Vec3(2.619000,2.611000,1.237000);
+positions[845] = Vec3(2.916000,2.726000,0.168000);
+positions[846] = Vec3(2.021000,2.833000,1.176000);
+positions[847] = Vec3(2.967000,2.308000,2.258000);
+positions[848] = Vec3(2.778000,2.270000,1.477000);
+positions[849] = Vec3(2.121000,1.834000,2.002000);
+positions[850] = Vec3(2.097000,2.752000,0.808000);
+positions[851] = Vec3(1.897000,0.566000,1.501000);
+positions[852] = Vec3(0.359000,2.802000,0.036000);
+positions[853] = Vec3(2.966000,2.454000,1.186000);
+positions[854] = Vec3(2.461000,2.964000,1.132000);
+positions[855] = Vec3(2.093000,1.821000,1.243000);
+positions[856] = Vec3(1.706000,2.659000,1.841000);
+positions[857] = Vec3(2.074000,1.709000,0.342000);
+positions[858] = Vec3(2.137000,2.894000,1.813000);
+positions[859] = Vec3(0.223000,2.293000,1.417000);
+positions[860] = Vec3(2.637000,0.007000,0.197000);
+positions[861] = Vec3(1.416000,0.050000,0.483000);
+positions[862] = Vec3(1.845000,2.250000,1.251000);
+positions[863] = Vec3(2.906000,0.034000,2.896000);
+positions[864] = Vec3(2.481000,0.204000,0.474000);
+positions[865] = Vec3(2.234000,2.051000,0.158000);
+positions[866] = Vec3(0.185000,2.453000,0.055000);
+positions[867] = Vec3(2.509000,0.048000,2.786000);
+positions[868] = Vec3(2.202000,2.206000,2.027000);
+positions[869] = Vec3(0.061000,2.367000,2.656000);
+positions[870] = Vec3(3.003000,2.755000,2.241000);
+positions[871] = Vec3(0.297000,2.131000,2.463000);
+positions[872] = Vec3(1.553000,0.429000,1.573000);
+positions[873] = Vec3(2.506000,1.832000,1.911000);
+positions[874] = Vec3(2.472000,1.814000,2.759000);
+positions[875] = Vec3(1.922000,1.563000,2.278000);
+positions[876] = Vec3(2.623000,2.666000,2.169000);
+positions[877] = Vec3(0.120000,1.834000,2.723000);
+positions[878] = Vec3(0.294000,0.103000,2.826000);
+positions[879] = Vec3(2.364000,2.821000,0.417000);
+positions[880] = Vec3(2.446000,1.734000,0.153000);
+positions[881] = Vec3(2.777000,2.037000,2.565000);
+positions[882] = Vec3(2.837000,2.477000,1.924000);
+positions[883] = Vec3(2.221000,1.961000,2.443000);
+positions[884] = Vec3(2.284000,2.895000,2.157000);
+positions[885] = Vec3(2.728000,2.880000,1.861000);
+positions[886] = Vec3(0.454000,2.080000,2.868000);
+positions[887] = Vec3(2.430000,2.790000,2.524000);
+positions[888] = Vec3(1.808000,2.213000,1.899000);
+positions[889] = Vec3(2.666000,0.053000,2.309000);
+positions[890] = Vec3(2.290000,2.408000,2.995000);
+positions[891] = Vec3(2.646000,2.592000,1.625000);
+positions[892] = Vec3(2.750000,2.508000,2.489000);
+positions[893] = Vec3(0.211000,1.753000,1.939000);