Commit 7f1da383 authored by Peter Eastman's avatar Peter Eastman
Browse files

Fixed bugs in OpenCL RPMD

parent a713aabf
...@@ -664,6 +664,10 @@ int OpenCLIntegrationUtilities::prepareRandomNumbers(int numValues) { ...@@ -664,6 +664,10 @@ int OpenCLIntegrationUtilities::prepareRandomNumbers(int numValues) {
randomPos += numValues; randomPos += numValues;
return oldPos; return oldPos;
} }
if (numValues > random->getSize()) {
delete random;
random = new OpenCLArray<mm_float4>(context, numValues, "random");
}
randomKernel.setArg<cl_int>(0, random->getSize()); randomKernel.setArg<cl_int>(0, random->getSize());
randomKernel.setArg<cl::Buffer>(1, random->getDeviceBuffer()); randomKernel.setArg<cl::Buffer>(1, random->getDeviceBuffer());
randomKernel.setArg<cl::Buffer>(2, randomSeed->getDeviceBuffer()); randomKernel.setArg<cl::Buffer>(2, randomSeed->getDeviceBuffer());
......
...@@ -130,7 +130,7 @@ void OpenCLIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDInte ...@@ -130,7 +130,7 @@ void OpenCLIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDInte
// Apply the PILE-L thermostat. // Apply the PILE-L thermostat.
const double dt = integrator.getStepSize(); const double dt = integrator.getStepSize();
pileKernel.setArg<cl_uint>(5, integration.prepareRandomNumbers(cl.getPaddedNumAtoms())); pileKernel.setArg<cl_uint>(5, integration.prepareRandomNumbers(numParticles*numCopies));
pileKernel.setArg<cl_float>(6, dt); pileKernel.setArg<cl_float>(6, dt);
pileKernel.setArg<cl_float>(7, integrator.getTemperature()*BOLTZ); pileKernel.setArg<cl_float>(7, integrator.getTemperature()*BOLTZ);
pileKernel.setArg<cl_float>(8, integrator.getFriction()); pileKernel.setArg<cl_float>(8, integrator.getFriction());
...@@ -158,6 +158,7 @@ void OpenCLIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDInte ...@@ -158,6 +158,7 @@ void OpenCLIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDInte
// Apply the PILE-L thermostat again. // Apply the PILE-L thermostat again.
pileKernel.setArg<cl_uint>(5, integration.prepareRandomNumbers(numParticles*numCopies));
cl.executeKernel(pileKernel, numParticles*numCopies, workgroupSize); cl.executeKernel(pileKernel, numParticles*numCopies, workgroupSize);
// Update the time and step count. // Update the time and step count.
...@@ -203,6 +204,8 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable ...@@ -203,6 +204,8 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable
int L = size; int L = size;
int m = 1; int m = 1;
string sign = (forward ? "1.0f" : "-1.0f"); string sign = (forward ? "1.0f" : "-1.0f");
string multReal = (forward ? "multiplyComplexRealPart" : "multiplyComplexRealPartConj");
string multImag = (forward ? "multiplyComplexImagPart" : "multiplyComplexImagPartConj");
source<<"{\n"; source<<"{\n";
source<<"__local float4* real0 = "<<variable<<"real;\n"; source<<"__local float4* real0 = "<<variable<<"real;\n";
...@@ -257,14 +260,14 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable ...@@ -257,14 +260,14 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable
source<<"float4 d10i = "<<sign<<"*(d3r-"<<coeff<<"*d2r);\n"; source<<"float4 d10i = "<<sign<<"*(d3r-"<<coeff<<"*d2r);\n";
source<<"real"<<output<<"[i+4*j*"<<m<<"] = c0r+d4r;\n"; source<<"real"<<output<<"[i+4*j*"<<m<<"] = c0r+d4r;\n";
source<<"imag"<<output<<"[i+4*j*"<<m<<"] = c0i+d4i;\n"; source<<"imag"<<output<<"[i+4*j*"<<m<<"] = c0i+d4i;\n";
source<<"real"<<output<<"[i+(4*j+1)*"<<m<<"] = multiplyComplexRealPart(w[j*"<<size<<"/"<<(5*L)<<"], d7r+d9r, d7i+d9i);\n"; source<<"real"<<output<<"[i+(4*j+1)*"<<m<<"] = "<<multReal<<"(w[j*"<<size<<"/"<<(5*L)<<"], d7r+d9r, d7i+d9i);\n";
source<<"imag"<<output<<"[i+(4*j+1)*"<<m<<"] = multiplyComplexImagPart(w[j*"<<size<<"/"<<(5*L)<<"], d7r+d9r, d7i+d9i);\n"; source<<"imag"<<output<<"[i+(4*j+1)*"<<m<<"] = "<<multImag<<"(w[j*"<<size<<"/"<<(5*L)<<"], d7r+d9r, d7i+d9i);\n";
source<<"real"<<output<<"[i+(4*j+2)*"<<m<<"] = multiplyComplexRealPart(w[j*"<<(2*size)<<"/"<<(5*L)<<"], d8r+d10r, d8i+d10i);\n"; source<<"real"<<output<<"[i+(4*j+2)*"<<m<<"] = "<<multReal<<"(w[j*"<<(2*size)<<"/"<<(5*L)<<"], d8r+d10r, d8i+d10i);\n";
source<<"imag"<<output<<"[i+(4*j+2)*"<<m<<"] = multiplyComplexImagPart(w[j*"<<(2*size)<<"/"<<(5*L)<<"], d8r+d10r, d8i+d10i);\n"; source<<"imag"<<output<<"[i+(4*j+2)*"<<m<<"] = "<<multImag<<"(w[j*"<<(2*size)<<"/"<<(5*L)<<"], d8r+d10r, d8i+d10i);\n";
source<<"real"<<output<<"[i+(4*j+3)*"<<m<<"] = multiplyComplexRealPart(w[j*"<<(3*size)<<"/"<<(5*L)<<"], d8r-d10r, d8i-d10i);\n"; source<<"real"<<output<<"[i+(4*j+3)*"<<m<<"] = "<<multReal<<"(w[j*"<<(3*size)<<"/"<<(5*L)<<"], d8r-d10r, d8i-d10i);\n";
source<<"imag"<<output<<"[i+(4*j+3)*"<<m<<"] = multiplyComplexImagPart(w[j*"<<(3*size)<<"/"<<(5*L)<<"], d8r-d10r, d8i-d10i);\n"; source<<"imag"<<output<<"[i+(4*j+3)*"<<m<<"] = "<<multImag<<"(w[j*"<<(3*size)<<"/"<<(5*L)<<"], d8r-d10r, d8i-d10i);\n";
source<<"real"<<output<<"[i+(4*j+4)*"<<m<<"] = multiplyComplexRealPart(w[j*"<<(4*size)<<"/"<<(5*L)<<"], d7r-d9r, d7i-d9i);\n"; source<<"real"<<output<<"[i+(4*j+4)*"<<m<<"] = "<<multReal<<"(w[j*"<<(4*size)<<"/"<<(5*L)<<"], d7r-d9r, d7i-d9i);\n";
source<<"imag"<<output<<"[i+(4*j+4)*"<<m<<"] = multiplyComplexImagPart(w[j*"<<(4*size)<<"/"<<(5*L)<<"], d7r-d9r, d7i-d9i);\n"; source<<"imag"<<output<<"[i+(4*j+4)*"<<m<<"] = "<<multImag<<"(w[j*"<<(4*size)<<"/"<<(5*L)<<"], d7r-d9r, d7i-d9i);\n";
source<<"}\n"; source<<"}\n";
m = m*5; m = m*5;
unfactored /= 5; unfactored /= 5;
...@@ -293,12 +296,12 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable ...@@ -293,12 +296,12 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable
source<<"float4 d3i = "<<sign<<"*(c3r-c1r);\n"; source<<"float4 d3i = "<<sign<<"*(c3r-c1r);\n";
source<<"real"<<output<<"[i+3*j*"<<m<<"] = d0r+d2r;\n"; source<<"real"<<output<<"[i+3*j*"<<m<<"] = d0r+d2r;\n";
source<<"imag"<<output<<"[i+3*j*"<<m<<"] = d0i+d2i;\n"; source<<"imag"<<output<<"[i+3*j*"<<m<<"] = d0i+d2i;\n";
source<<"real"<<output<<"[i+(3*j+1)*"<<m<<"] = multiplyComplexRealPart(w[j*"<<size<<"/"<<(4*L)<<"], d1r+d3r, d1i+d3i);\n"; source<<"real"<<output<<"[i+(3*j+1)*"<<m<<"] = "<<multReal<<"(w[j*"<<size<<"/"<<(4*L)<<"], d1r+d3r, d1i+d3i);\n";
source<<"imag"<<output<<"[i+(3*j+1)*"<<m<<"] = multiplyComplexImagPart(w[j*"<<size<<"/"<<(4*L)<<"], d1r+d3r, d1i+d3i);\n"; source<<"imag"<<output<<"[i+(3*j+1)*"<<m<<"] = "<<multImag<<"(w[j*"<<size<<"/"<<(4*L)<<"], d1r+d3r, d1i+d3i);\n";
source<<"real"<<output<<"[i+(3*j+2)*"<<m<<"] = multiplyComplexRealPart(w[j*"<<(2*size)<<"/"<<(4*L)<<"], d0r-d2r, d0i-d2i);\n"; source<<"real"<<output<<"[i+(3*j+2)*"<<m<<"] = "<<multReal<<"(w[j*"<<(2*size)<<"/"<<(4*L)<<"], d0r-d2r, d0i-d2i);\n";
source<<"imag"<<output<<"[i+(3*j+2)*"<<m<<"] = multiplyComplexImagPart(w[j*"<<(2*size)<<"/"<<(4*L)<<"], d0r-d2r, d0i-d2i);\n"; source<<"imag"<<output<<"[i+(3*j+2)*"<<m<<"] = "<<multImag<<"(w[j*"<<(2*size)<<"/"<<(4*L)<<"], d0r-d2r, d0i-d2i);\n";
source<<"real"<<output<<"[i+(3*j+3)*"<<m<<"] = multiplyComplexRealPart(w[j*"<<(3*size)<<"/"<<(4*L)<<"], d1r-d3r, d1i-d3i);\n"; source<<"real"<<output<<"[i+(3*j+3)*"<<m<<"] = "<<multReal<<"(w[j*"<<(3*size)<<"/"<<(4*L)<<"], d1r-d3r, d1i-d3i);\n";
source<<"imag"<<output<<"[i+(3*j+3)*"<<m<<"] = multiplyComplexImagPart(w[j*"<<(3*size)<<"/"<<(4*L)<<"], d1r-d3r, d1i-d3i);\n"; source<<"imag"<<output<<"[i+(3*j+3)*"<<m<<"] = "<<multImag<<"(w[j*"<<(3*size)<<"/"<<(4*L)<<"], d1r-d3r, d1i-d3i);\n";
source<<"}\n"; source<<"}\n";
m = m*4; m = m*4;
unfactored /= 4; unfactored /= 4;
...@@ -323,10 +326,10 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable ...@@ -323,10 +326,10 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable
source<<"float4 d2i = "<<sign<<"*"<<OpenCLExpressionUtilities::doubleToString(sin(M_PI/3.0))<<"*(c2r-c1r);\n"; source<<"float4 d2i = "<<sign<<"*"<<OpenCLExpressionUtilities::doubleToString(sin(M_PI/3.0))<<"*(c2r-c1r);\n";
source<<"real"<<output<<"[i+2*j*"<<m<<"] = c0r+d0r;\n"; source<<"real"<<output<<"[i+2*j*"<<m<<"] = c0r+d0r;\n";
source<<"imag"<<output<<"[i+2*j*"<<m<<"] = c0i+d0i;\n"; source<<"imag"<<output<<"[i+2*j*"<<m<<"] = c0i+d0i;\n";
source<<"real"<<output<<"[i+(2*j+1)*"<<m<<"] = multiplyComplexRealPart(w[j*"<<size<<"/"<<(3*L)<<"], d1r+d2r, d1i+d2i);\n"; source<<"real"<<output<<"[i+(2*j+1)*"<<m<<"] = "<<multReal<<"(w[j*"<<size<<"/"<<(3*L)<<"], d1r+d2r, d1i+d2i);\n";
source<<"imag"<<output<<"[i+(2*j+1)*"<<m<<"] = multiplyComplexImagPart(w[j*"<<size<<"/"<<(3*L)<<"], d1r+d2r, d1i+d2i);\n"; source<<"imag"<<output<<"[i+(2*j+1)*"<<m<<"] = "<<multImag<<"(w[j*"<<size<<"/"<<(3*L)<<"], d1r+d2r, d1i+d2i);\n";
source<<"real"<<output<<"[i+(2*j+2)*"<<m<<"] = multiplyComplexRealPart(w[j*"<<(2*size)<<"/"<<(3*L)<<"], d1r-d2r, d1i-d2i);\n"; source<<"real"<<output<<"[i+(2*j+2)*"<<m<<"] = "<<multReal<<"(w[j*"<<(2*size)<<"/"<<(3*L)<<"], d1r-d2r, d1i-d2i);\n";
source<<"imag"<<output<<"[i+(2*j+2)*"<<m<<"] = multiplyComplexImagPart(w[j*"<<(2*size)<<"/"<<(3*L)<<"], d1r-d2r, d1i-d2i);\n"; source<<"imag"<<output<<"[i+(2*j+2)*"<<m<<"] = "<<multImag<<"(w[j*"<<(2*size)<<"/"<<(3*L)<<"], d1r-d2r, d1i-d2i);\n";
source<<"}\n"; source<<"}\n";
m = m*3; m = m*3;
unfactored /= 3; unfactored /= 3;
...@@ -343,8 +346,8 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable ...@@ -343,8 +346,8 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable
source<<"float4 c1i = imag"<<input<<"[i+"<<(L*m)<<"];\n"; source<<"float4 c1i = imag"<<input<<"[i+"<<(L*m)<<"];\n";
source<<"real"<<output<<"[i+j*"<<m<<"] = c0r+c1r;\n"; source<<"real"<<output<<"[i+j*"<<m<<"] = c0r+c1r;\n";
source<<"imag"<<output<<"[i+j*"<<m<<"] = c0i+c1i;\n"; source<<"imag"<<output<<"[i+j*"<<m<<"] = c0i+c1i;\n";
source<<"real"<<output<<"[i+(j+1)*"<<m<<"] = multiplyComplexRealPart(w[j*"<<size<<"/"<<(2*L)<<"], c0r-c1r, c0i-c1i);\n"; source<<"real"<<output<<"[i+(j+1)*"<<m<<"] = "<<multReal<<"(w[j*"<<size<<"/"<<(2*L)<<"], c0r-c1r, c0i-c1i);\n";
source<<"imag"<<output<<"[i+(j+1)*"<<m<<"] = multiplyComplexImagPart(w[j*"<<size<<"/"<<(2*L)<<"], c0r-c1r, c0i-c1i);\n"; source<<"imag"<<output<<"[i+(j+1)*"<<m<<"] = "<<multImag<<"(w[j*"<<size<<"/"<<(2*L)<<"], c0r-c1r, c0i-c1i);\n";
source<<"}\n"; source<<"}\n";
m = m*2; m = m*2;
unfactored /= 2; unfactored /= 2;
...@@ -359,8 +362,8 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable ...@@ -359,8 +362,8 @@ string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable
// Create the kernel. // Create the kernel.
if (stage%2 == 1) { if (stage%2 == 1) {
source<<variable<<"real[indexInBlock] = real1[indexInBlock];\n"; source<<"real0[indexInBlock] = real1[indexInBlock];\n";
source<<variable<<"imag[indexInBlock] = imag1[indexInBlock];\n"; source<<"imag0[indexInBlock] = imag1[indexInBlock];\n";
} }
source<<"}\n"; source<<"}\n";
return source.str(); return source.str();
......
...@@ -10,6 +10,14 @@ float4 multiplyComplexImagPart(float2 c1, float4 c2r, float4 c2i) { ...@@ -10,6 +10,14 @@ float4 multiplyComplexImagPart(float2 c1, float4 c2r, float4 c2i) {
return c1.x*c2i+c1.y*c2r; return c1.x*c2i+c1.y*c2r;
} }
float4 multiplyComplexRealPartConj(float2 c1, float4 c2r, float4 c2i) {
return c1.x*c2r+c1.y*c2i;
}
float4 multiplyComplexImagPartConj(float2 c1, float4 c2r, float4 c2i) {
return c1.x*c2i-c1.y*c2r;
}
/** /**
* Apply the PILE-L thermostat. * Apply the PILE-L thermostat.
*/ */
...@@ -25,8 +33,9 @@ __kernel void applyPileThermostat(__global float4* velm, __local float4* v, __lo ...@@ -25,8 +33,9 @@ __kernel void applyPileThermostat(__global float4* velm, __local float4* v, __lo
__local float4* vreal = &v[blockStart]; __local float4* vreal = &v[blockStart];
__local float4* vimag = &v[blockStart+get_local_size(0)]; __local float4* vimag = &v[blockStart+get_local_size(0)];
if (get_local_id(0) < NUM_COPIES) if (get_local_id(0) < NUM_COPIES)
w[get_local_id(0)] = (float2) (cos(get_local_id(0)*2*M_PI/NUM_COPIES), sin(-get_local_id(0)*2*M_PI/NUM_COPIES)); w[indexInBlock] = (float2) (cos(-indexInBlock*2*M_PI/NUM_COPIES), sin(-indexInBlock*2*M_PI/NUM_COPIES));
barrier(CLK_LOCAL_MEM_FENCE); barrier(CLK_LOCAL_MEM_FENCE);
randomIndex += blockStart;
for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) { for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
float4 particleVelm = velm[particle+indexInBlock*PADDED_NUM_ATOMS]; float4 particleVelm = velm[particle+indexInBlock*PADDED_NUM_ATOMS];
float invMass = particleVelm.w; float invMass = particleVelm.w;
...@@ -66,7 +75,7 @@ __kernel void applyPileThermostat(__global float4* velm, __local float4* v, __lo ...@@ -66,7 +75,7 @@ __kernel void applyPileThermostat(__global float4* velm, __local float4* v, __lo
FFT_V_BACKWARD FFT_V_BACKWARD
velm[particle+indexInBlock*PADDED_NUM_ATOMS].xyz = SCALE*vreal[indexInBlock].xyz; velm[particle+indexInBlock*PADDED_NUM_ATOMS].xyz = SCALE*vreal[indexInBlock].xyz;
randomIndex += numBlocks*NUM_COPIES; randomIndex += get_global_size(0);
} }
} }
...@@ -97,7 +106,7 @@ __kernel void integrateStep(__global float4* posq, __global float4* velm, __glob ...@@ -97,7 +106,7 @@ __kernel void integrateStep(__global float4* posq, __global float4* velm, __glob
__local float4* vreal = &v[blockStart]; __local float4* vreal = &v[blockStart];
__local float4* vimag = &v[blockStart+get_local_size(0)]; __local float4* vimag = &v[blockStart+get_local_size(0)];
if (get_local_id(0) < NUM_COPIES) if (get_local_id(0) < NUM_COPIES)
w[get_local_id(0)] = (float2) (cos(get_local_id(0)*2*M_PI/NUM_COPIES), sin(-get_local_id(0)*2*M_PI/NUM_COPIES)); w[indexInBlock] = (float2) (cos(-indexInBlock*2*M_PI/NUM_COPIES), sin(-indexInBlock*2*M_PI/NUM_COPIES));
barrier(CLK_LOCAL_MEM_FENCE); barrier(CLK_LOCAL_MEM_FENCE);
for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) { for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
float4 particlePosq = posq[particle+indexInBlock*PADDED_NUM_ATOMS]; float4 particlePosq = posq[particle+indexInBlock*PADDED_NUM_ATOMS];
......
...@@ -6,6 +6,8 @@ ENABLE_TESTING() ...@@ -6,6 +6,8 @@ ENABLE_TESTING()
INCLUDE_DIRECTORIES(${OPENCL_INCLUDE_DIR}) INCLUDE_DIRECTORIES(${OPENCL_INCLUDE_DIR})
SET(SHARED_OPENMM_RPMD_TARGET OpenMMRPMD)
# Automatically create tests using files named "Test*.cpp" # Automatically create tests using files named "Test*.cpp"
FILE(GLOB TEST_PROGS "*Test*.cpp") FILE(GLOB TEST_PROGS "*Test*.cpp")
FOREACH(TEST_PROG ${TEST_PROGS}) FOREACH(TEST_PROG ${TEST_PROGS})
...@@ -13,7 +15,7 @@ FOREACH(TEST_PROG ${TEST_PROGS}) ...@@ -13,7 +15,7 @@ FOREACH(TEST_PROG ${TEST_PROGS})
# Link with shared library # Link with shared library
ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG}) ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET}) TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_RPMD_TARGET})
ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT}) ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
ENDFOREACH(TEST_PROG ${TEST_PROGS}) ENDFOREACH(TEST_PROG ${TEST_PROGS})
/* -------------------------------------------------------------------------- *
* 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) 2011 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 OpenCL implementation of RPMDIntegrator.
*/
#include "../../../tests/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/Platform.h"
#include "openmm/System.h"
#include "openmm/RPMDIntegrator.h"
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
void testIntegration() {
const int numParticles = 1;
const int numCopies = 25;
const double temperature = 300.0;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(i+1);
if (numParticles > 1) {
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
for (int i = 0; i < numParticles-1; i++)
bonds->addBond(i, i+1, 1.0, 100.0);
}
RPMDIntegrator integ(numCopies, temperature, 1.0, 0.001);
Platform& platform = Platform::getPlatformByName("OpenCL");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numCopies; i++)
{
for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(j+0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt));
integ.setPositions(i, positions);
}
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
vector<double> rg(numParticles, 0.0);
const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
const double wn = numCopies*BOLTZ*temperature/hbar;
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Positions | State::Velocities | State::Energy);
for (int j = 0; j < numCopies; j++)
ke[j] += state[j].getKineticEnergy();
double totalEnergy = 0.0;
for (int j = 0; j < numCopies; j++) {
totalEnergy += state[j].getKineticEnergy()+state[j].getPotentialEnergy();
for (int k = 0; k < numParticles; k++) {
Vec3 delta = state[j].getPositions()[k]-state[j == 0 ? numCopies-1 : j-1].getPositions()[k];
totalEnergy += 0.5*system.getParticleMass(k)*wn*wn*delta.dot(delta);
}
}
for (int j = 0; j < numParticles; j++) {
double rg2 = 0.0;
for (int k = 0; k < numCopies; k++)
for (int m = 0; m < numCopies; m++) {
Vec3 delta = state[k].getPositions()[j]-state[m].getPositions()[j];
rg2 += delta.dot(delta);
}
rg[j] += sqrt(rg2/(2*numCopies*numCopies));
}
}
for (int i = 0; i < numCopies; i++) {
double value = ke[i]/numSteps;
double expected = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
printf("%d: %g %g %g\n", i, value, expected, value/expected);
}
for (int i = 0; i < numParticles; i++) {
double value = rg[i]/numSteps;
double expected = hbar/(2*sqrt(system.getParticleMass(i)*BOLTZ*temperature));
printf("%d: %g %g %g\n", i, value, expected, value/expected);
}
}
int main() {
try {
Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
testIntegration();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "Done" << std::endl;
return 0;
}
...@@ -8,8 +8,7 @@ INCLUDE_DIRECTORIES(${OPENMM_DIR}/openmmapi/include/openmm) ...@@ -8,8 +8,7 @@ INCLUDE_DIRECTORIES(${OPENMM_DIR}/openmmapi/include/openmm)
INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/src) INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/src)
#INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/src/kernels) #INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/src/kernels)
Set( SHARED_OPENMM_RPMD_TARGET OpenMMRPMD) SET(SHARED_OPENMM_RPMD_TARGET OpenMMRPMD)
Set( SHARED_CUDA_TARGET OpenMMRPMDCuda )
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug) IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(SHARED_CUDA_TARGET ${SHARED_CUDA_TARGET}_d) SET(SHARED_CUDA_TARGET ${SHARED_CUDA_TARGET}_d)
...@@ -26,6 +25,6 @@ FOREACH(TEST_PROG ${TEST_PROGS}) ...@@ -26,6 +25,6 @@ FOREACH(TEST_PROG ${TEST_PROGS})
# Link with shared library # Link with shared library
ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG}) ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_OPENMM_RPMD_TARGET} ) TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_OPENMM_RPMD_TARGET})
ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT}) ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
ENDFOREACH(TEST_PROG ${TEST_PROGS}) ENDFOREACH(TEST_PROG ${TEST_PROGS})
...@@ -59,14 +59,15 @@ void testIntegration() { ...@@ -59,14 +59,15 @@ void testIntegration() {
for (int i = 0; i < numParticles-1; i++) for (int i = 0; i < numParticles-1; i++)
bonds->addBond(i, i+1, 1.0, 100.0); bonds->addBond(i, i+1, 1.0, 100.0);
RPMDIntegrator integ(numCopies, temperature, 1.0, 0.001); RPMDIntegrator integ(numCopies, temperature, 1.0, 0.001);
Context context(system, integ); Platform& platform = Platform::getPlatformByName("Reference");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt; OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt); init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles); vector<Vec3> positions(numParticles);
for (int i = 0; i < numCopies; i++) for (int i = 0; i < numCopies; i++)
{ {
for (int j = 0; j < numParticles; j++) for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(i+0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt)); positions[j] = Vec3(j+0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt));
integ.setPositions(i, positions); integ.setPositions(i, positions);
} }
const int numSteps = 100000; const int numSteps = 100000;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment