Created CUDA implementation of MonteCarloBarostat

platforms/cuda/src/CudaKernelFactory.cpp

@@ -73,6 +73,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaIntegrateVariableLangevinStepKernel(name, platform, data);
     if (name == ApplyAndersenThermostatKernel::Name())
         return new CudaApplyAndersenThermostatKernel(name, platform, data);
+    if (name == ApplyMonteCarloBarostatKernel::Name())
+        return new CudaApplyMonteCarloBarostatKernel(name, platform, data);
     if (name == CalcKineticEnergyKernel::Name())
         return new CudaCalcKineticEnergyKernel(name, platform, data);
     if (name == RemoveCMMotionKernel::Name())

platforms/cuda/src/CudaKernels.cpp

@@ -1042,6 +1042,46 @@ void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
     kCalculateAndersenThermostat(gpu);
 }
 
+CudaApplyMonteCarloBarostatKernel::~CudaApplyMonteCarloBarostatKernel() {
+    if (moleculeAtoms != NULL)
+        delete moleculeAtoms;
+    if (moleculeStartIndex != NULL)
+        delete moleculeStartIndex;
+}
+
+void CudaApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& thermostat) {
+}
+
+void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
+    if (!hasInitializedMolecules) {
+        hasInitializedMolecules = true;
+
+        // Create the arrays with the molecule definitions.
+
+        vector<vector<int> > molecules = context.getMolecules();
+        numMolecules = molecules.size();
+        moleculeAtoms = new CUDAStream<int>(context.getSystem().getNumParticles(), 1, "moleculeAtoms");
+        moleculeStartIndex = new CUDAStream<int>(numMolecules+1, 1, "moleculeStartIndex");
+        int index = 0;
+        for (int i = 0; i < numMolecules; i++) {
+            (*moleculeStartIndex)[i] = index;
+            for (int j = 0; j < (int) molecules[i].size(); j++)
+                (*moleculeAtoms)[index++] = molecules[i][j];
+        }
+        (*moleculeStartIndex)[numMolecules] = index;
+        moleculeAtoms->Upload();
+        moleculeStartIndex->Upload();
+    }
+    _gpuContext* gpu = data.gpu;
+    gpu->psPosqP4->CopyFrom(*gpu->psPosq4);
+    kScaleAtomCoordinates(gpu, scale, *moleculeAtoms, *moleculeStartIndex);
+}
+
+void CudaApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
+    _gpuContext* gpu = data.gpu;
+    gpu->psPosq4->CopyFrom(*gpu->psPosqP4);
+}
+
 void CudaCalcKineticEnergyKernel::initialize(const System& system) {
     int numParticles = system.getNumParticles();
     masses.resize(numParticles);

platforms/cuda/src/CudaKernels.h

@@ -740,6 +740,47 @@ private:
     double prevTemp, prevFrequency, prevStepSize;
 };
 
+/**
+ * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
+ */
+class CudaApplyMonteCarloBarostatKernel : public ApplyMonteCarloBarostatKernel {
+public:
+    CudaApplyMonteCarloBarostatKernel(std::string name, const Platform& platform, CudaPlatform::PlatformData& data) : ApplyMonteCarloBarostatKernel(name, platform), data(data),
+            moleculeAtoms(NULL), moleculeStartIndex(NULL) {
+    }
+    ~CudaApplyMonteCarloBarostatKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param barostat   the MonteCarloBarostat this kernel will be used for
+     */
+    void initialize(const System& system, const MonteCarloBarostat& barostat);
+    /**
+     * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
+     * This is called BEFORE the periodic box size is modified.  It should begin by translating each particle
+     * or cluster into the first periodic box, so that coordinates will still be correct after the box size
+     * is changed.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param scale      the scale factor by which to multiply particle positions
+     */
+    void scaleCoordinates(ContextImpl& context, double scale);
+    /**
+     * Reject the most recent Monte Carlo step, restoring the particle positions to where they were before
+     * scaleCoordinates() was last called.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    void restoreCoordinates(ContextImpl& context);
+private:
+    CudaPlatform::PlatformData& data;
+    bool hasInitializedMolecules;
+    int numMolecules;
+    CUDAStream<int>* moleculeAtoms;
+    CUDAStream<int>* moleculeStartIndex;
+};
+
 /**
  * This kernel is invoked to calculate the kinetic energy of the system.
  */

platforms/cuda/src/CudaPlatform.cpp

@@ -65,6 +65,7 @@ CudaPlatform::CudaPlatform() {
     registerKernelFactory(IntegrateVariableVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);
     registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
+    registerKernelFactory(ApplyMonteCarloBarostatKernel::Name(), factory);
     registerKernelFactory(CalcKineticEnergyKernel::Name(), factory);
     registerKernelFactory(RemoveCMMotionKernel::Name(), factory);
     platformProperties.push_back(CudaDevice());

platforms/cuda/src/kernels/cudaKernels.h

@@ -68,6 +68,7 @@ extern void kVerletUpdatePart2(gpuContext gpu);
 extern void kSelectVerletStepSize(gpuContext gpu, float maxTimeStep);
 extern void kBrownianUpdatePart1(gpuContext gpu);
 extern void kBrownianUpdatePart2(gpuContext gpu);
+extern void kScaleAtomCoordinates(gpuContext gpu, float scale, CUDAStream<int>& moleculeAtoms, CUDAStream<int>& moleculeStartIndex);
 
 // Extras
 extern void kReduceForces(gpuContext gpu);

platforms/cuda/src/kernels/cudatypes.h

@@ -82,6 +82,7 @@ struct CUDAStream : public SoADeviceObject
     void Deallocate();
     void Upload();
     void Download();
+    void CopyFrom(const CUDAStream<T>& src);
     void Collapse(unsigned int newstreams = 1, unsigned int interleave = 1);
     T& operator[](int index);
 };
@@ -166,6 +167,14 @@ void CUDAStream<T>::Download()
     RTERROR(status, (_name+": cudaMemcpy in CUDAStream::Download failed").c_str());
 }
 
+template <typename T>
+void CUDAStream<T>::CopyFrom(const CUDAStream<T>& src)
+{
+    cudaError_t status;
+    status = cudaMemcpy(_pDevData, src._pDevData, _stride * _subStreams * sizeof(T), cudaMemcpyDeviceToDevice);
+    RTERROR(status, (_name+": cudaMemcpy in CUDAStream::Copy failed").c_str());
+}
+
 template <typename T>
 void CUDAStream<T>::Collapse(unsigned int newstreams, unsigned int interleave)
 {

platforms/cuda/src/kernels/kMonteCarloBarostat.cu

+/* -------------------------------------------------------------------------- *
+ *                                   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) 2010 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include <stdio.h>
+#include <cuda.h>
+#include <vector_functions.h>
+#include <cstdlib>
+#include <string>
+#include <iostream>
+//#include <fstream>
+using namespace std;
+
+#include "gputypes.h"
+
+__global__ void kScaleAtomCoordinates_kernel(float scale, int numMolecules, float3 periodicBoxSize, float4* posq, int* moleculeAtoms, int* moleculeStartIndex) {
+    float3 invPeriodicBoxSize = make_float3(1.0f/periodicBoxSize.x, 1.0f/periodicBoxSize.y, 1.0f/periodicBoxSize.z);
+    for (int index = threadIdx.x+blockIdx.x*blockDim.x; index < numMolecules; index += blockDim.x*gridDim.x) {
+        int first = moleculeStartIndex[index];
+        int last = moleculeStartIndex[index+1];
+        int numAtoms = last-first;
+
+        // Find the center of each molecule.
+
+        float3 center = make_float3(0, 0, 0);
+        for (int atom = first; atom < last; atom++) {
+            float4 pos = posq[moleculeAtoms[atom]];
+            center.x += pos.x;
+            center.y += pos.y;
+            center.z += pos.z;
+        }
+        center.x /= (float) numAtoms;
+        center.y /= (float) numAtoms;
+        center.z /= (float) numAtoms;
+
+        // Move it into the first periodic box.
+
+        int xcell = (int) floor(center.x*invPeriodicBoxSize.x);
+        int ycell = (int) floor(center.y*invPeriodicBoxSize.y);
+        int zcell = (int) floor(center.z*invPeriodicBoxSize.z);
+        float3 delta = make_float3(xcell*periodicBoxSize.x, ycell*periodicBoxSize.y, zcell*periodicBoxSize.z);
+        center.x -= delta.x;
+        center.y -= delta.y;
+        center.z -= delta.z;
+
+        // Now scale the position of the molecule center.
+
+        delta.x = center.x*(scale-1)-delta.x;
+        delta.y = center.y*(scale-1)-delta.y;
+        delta.z = center.z*(scale-1)-delta.z;
+        for (int atom = first; atom < last; atom++) {
+            float4 pos = posq[moleculeAtoms[atom]];
+            pos.x += delta.x;
+            pos.y += delta.y;
+            pos.z += delta.z;
+            posq[moleculeAtoms[atom]] = pos;
+        }
+    }
+}
+
+void kScaleAtomCoordinates(gpuContext gpu, float scale, CUDAStream<int>& moleculeAtoms, CUDAStream<int>& moleculeStartIndex)
+{
+//    printf("kScaleAtomCoordinates\n");
+    kScaleAtomCoordinates_kernel<<<gpu->sim.blocks, gpu->sim.update_threads_per_block>>>(scale, moleculeStartIndex._length-1,
+            make_float3(gpu->sim.periodicBoxSizeX, gpu->sim.periodicBoxSizeY, gpu->sim.periodicBoxSizeZ), gpu->sim.pPosq,
+            moleculeAtoms._pDevData, moleculeStartIndex._pDevData);
+    LAUNCHERROR("kScaleAtomCoordinates");
+}
+

platforms/cuda/tests/TestCudaMonteCarloBarostat.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-2010 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 Cuda implementation of MonteCarloBarostat.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/MonteCarloBarostat.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testChangingBoxSize() {
+    CudaPlatform platform;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 5, 0), Vec3(0, 0, 6));
+    system.addParticle(1.0);
+    LangevinIntegrator integrator(300.0, 1.0, 0.01);
+    Context context(system, integrator, platform);
+    Vec3 x, y, z;
+    context.getPeriodicBoxVectors(x, y, z);
+    ASSERT_EQUAL_VEC(Vec3(4, 0, 0), x, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 5, 0), y, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 6), z, 0);
+    context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 9));
+    context.getPeriodicBoxVectors(x, y, z);
+    ASSERT_EQUAL_VEC(Vec3(7, 0, 0), x, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 8, 0), y, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 9), z, 0);
+}
+
+void testIdealGas() {
+    const int numParticles = 64;
+    const int frequency = 10;
+    const int steps = 1000;
+    const double pressure = 1.5;
+    const double pressureInMD = pressure/(AVOGADRO*1e-25);
+    const double temp[] = {300.0, 600.0, 1000.0};
+    const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
+    const double initialLength = std::pow(initialVolume, 1.0/3.0);
+
+    // Create a gas of noninteracting particles.
+
+    CudaPlatform platform;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(1.0);
+        positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
+    }
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency);
+    system.addForce(barostat);
+
+    // Test it for three different temperatures.
+
+    for (int i = 0; i < 3; i++) {
+        barostat->setTemperature(temp[i]);
+        LangevinIntegrator integrator(temp[i], 0.1, 0.01);
+        Context context(system, integrator, platform);
+        context.setPositions(positions);
+
+        // Let it equilibrate.
+
+        integrator.step(10000);
+
+        // Now run it for a while and see if the volume is correct.
+
+        double volume = 0.0;
+        for (int j = 0; j < steps; ++j) {
+            Vec3 box[3];
+            context.getPeriodicBoxVectors(box[0], box[1], box[2]);
+            volume += box[0][0]*box[1][1]*box[2][2];
+            ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
+            ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
+            integrator.step(frequency);
+        }
+        volume /= steps;
+        double expected = numParticles*BOLTZ*temp[i]/pressureInMD;//+numParticles*(4*M_PI/3)*sigma*sigma*sigma;
+        ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt(steps));
+    }
+}
+
+void testRandomSeed() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    const double pressure = 1.5;
+    CudaPlatform platform;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(8, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 8));
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* forceField = new NonbondedForce();
+    forceField->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(2.0);
+        forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
+    }
+    system.addForce(forceField);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, 1);
+    system.addForce(barostat);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
+        velocities[i] = Vec3(0, 0, 0);
+    }
+
+    // Try twice with the same random seed.
+
+    barostat->setRandomNumberSeed(5);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state1 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state2 = context.getState(State::Positions);
+
+    // Try twice with a different random seed.
+
+    barostat->setRandomNumberSeed(10);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state3 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state4 = context.getState(State::Positions);
+
+    // Compare the results.
+
+    for (int i = 0; i < numParticles; i++) {
+        for (int j = 0; j < 3; j++) {
+            ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
+            ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
+            ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
+        }
+    }
+}
+
+void testWater() {
+    const int gridSize = 8;
+    const int numMolecules = gridSize*gridSize*gridSize;
+    const int frequency = 10;
+    const int steps = 400;
+    const double temp = 273.15;
+    const double pressure = 3;
+    const double spacing = 0.32;
+    const double angle = 109.47*M_PI/180;
+    const double dOH = 0.1;
+    const double dHH = dOH*2*std::sin(0.5*angle);
+
+    // Create a box of SPC water molecules.
+
+    CudaPlatform platform;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(gridSize*spacing, 0, 0), Vec3(0, gridSize*spacing, 0), Vec3(0, 0, gridSize*spacing));
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    vector<Vec3> positions;
+    Vec3 offset1(dOH, 0, 0);
+    Vec3 offset2(dOH*std::cos(angle), dOH*std::sin(angle), 0);
+    for (int i = 0; i < gridSize; ++i) {
+        for (int j = 0; j < gridSize; ++j) {
+            for (int k = 0; k < gridSize; ++k) {
+                int firstParticle = system.getNumParticles();
+                system.addParticle(16.0);
+                system.addParticle(1.0);
+                system.addParticle(1.0);
+                nonbonded->addParticle(-0.82, 0.316557, 0.650194);
+                nonbonded->addParticle(0.41, 1, 0);
+                nonbonded->addParticle(0.41, 1, 0);
+                Vec3 pos = Vec3(spacing*i, spacing*j, spacing*k);
+                positions.push_back(pos);
+                positions.push_back(pos+offset1);
+                positions.push_back(pos+offset2);
+                system.addConstraint(firstParticle, firstParticle+1, dOH);
+                system.addConstraint(firstParticle, firstParticle+2, dOH);
+                system.addConstraint(firstParticle+1, firstParticle+2, dHH);
+                nonbonded->addException(firstParticle, firstParticle+1, 0, 1, 0);
+                nonbonded->addException(firstParticle, firstParticle+2, 0, 1, 0);
+                nonbonded->addException(firstParticle+1, firstParticle+2, 0, 1, 0);
+            }
+        }
+    }
+    system.addForce(nonbonded);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency);
+    system.addForce(barostat);
+
+    // Simulate it and see if the density matches the expected value (1 g/mL).
+
+    LangevinIntegrator integrator(temp, 1.0, 0.002);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    integrator.step(2000);
+    double volume = 0.0;
+    for (int j = 0; j < steps; ++j) {
+        Vec3 box[3];
+        context.getPeriodicBoxVectors(box[0], box[1], box[2]);
+        volume += box[0][0]*box[1][1]*box[2][2];
+        integrator.step(frequency);
+    }
+    volume /= steps;
+    double density = numMolecules*18/(AVOGADRO*volume*1e-21);
+    ASSERT_USUALLY_EQUAL_TOL(1, density, 0.1);
+}
+
+int main() {
+    try {
+        testChangingBoxSize();
+        testIdealGas();
+        testRandomSeed();
+        testWater();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+

platforms/opencl/src/kernels/monteCarloBarostat.cl

@@ -21,7 +21,7 @@ __kernel void scalePositions(float scale, int numMolecules, float4 periodicBoxSi
         int xcell = (int) floor(center.x*invPeriodicBoxSize.x);
         int ycell = (int) floor(center.y*invPeriodicBoxSize.y);
         int zcell = (int) floor(center.z*invPeriodicBoxSize.z);
-        float4 delta = (float) (xcell*periodicBoxSize.x, ycell*periodicBoxSize.y, zcell*periodicBoxSize.z, 0);
+        float4 delta = (float4) (xcell*periodicBoxSize.x, ycell*periodicBoxSize.y, zcell*periodicBoxSize.z, 0);
         center -= delta;
 
         // Now scale the position of the molecule center.