Skip to content

GitLab

Commit fe1e6ffa authored by Rossen Apostolov's avatar Rossen Apostolov
Browse files

Completed Cuda implementation of Ewald method

parent 57606845
Hide whitespace changes
Inline Side-by-side
Showing with 1457 additions and 409 deletions
platforms/cuda/src/kernels/kCalculateCDLJEwaldReciprocal.h 0 → 100644
View file @ fe1e6ffa
/* -------------------------------------------------------------------------- *
* 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) 2009 Stanford University and the Authors. *
* Authors: Rossen P. Apostolov, 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/>. *
* -------------------------------------------------------------------------- */
/**
* This file contains the kernel for evaluating nonbonded forces using the
* Ewald summation method (Reciprocal space summation).
*/
__global__ void kCalculateCDLJEwaldReciprocalForces_kernel()
{
float alphaEwald = 3.123413;
float PI = 3.14159265358979323846f;
float TWO_PI = 2.0 * PI;
float SQRT_PI = sqrt(PI);
float eps0 = 5.72765E-4;
int numRx = 20+1;
int numRy = 20+1;
int numRz = 20+1;
int lowry, lowrz;
float kx, ky, kz, k2, ek;
float Qi, Qj, SinI, SinJ, CosI, CosJ;
float factorEwald = -1 / (4*alphaEwald*alphaEwald);
float recipBoxSizeX = TWO_PI / cSim.periodicBoxSizeX;
float recipBoxSizeY = TWO_PI / cSim.periodicBoxSizeY;
float recipBoxSizeZ = TWO_PI / cSim.periodicBoxSizeZ;
float V = cSim.periodicBoxSizeX * cSim.periodicBoxSizeY * cSim.periodicBoxSizeZ;
float4 apos1, apos2 ;
unsigned int atomID1 = threadIdx.x + blockIdx.x * blockDim.x;
while (atomID1 < cSim.atoms)
{
apos1 = cSim.pPosq[atomID1];
unsigned int atomID2 = 0;
while (atomID2 < cSim.atoms)
{
apos2 = cSim.pPosq[atomID2];
lowry = 0;
lowrz = 1;
for(int rx = 0; rx < numRx; rx++) {
kx = rx * recipBoxSizeX;
for(int ry = lowry; ry < numRy; ry++) {
ky = ry * recipBoxSizeY;
for (int rz = lowrz; rz < numRz; rz++) {
kz = rz * recipBoxSizeZ;
k2 = kx * kx + ky * ky + kz * kz;
ek = exp ( k2 * factorEwald);
Qi = apos1.w ;
Qj = apos2.w ;
SinI = sin ( kx * apos1.x + ky * apos1.y + kz * apos1.z );
SinJ = sin ( kx * apos2.x + ky * apos2.y + kz * apos2.z );
CosI = cos ( kx * apos1.x + ky * apos1.y + kz * apos1.z );
CosJ = cos ( kx * apos2.x + ky * apos2.y + kz * apos2.z );
cSim.pForce4[atomID1].x -= (2.0 / (V * eps0 )) * Qi * ( kx/k2) * ek * ( - SinI * Qj * CosJ + CosI * Qj * SinJ);
cSim.pForce4[atomID1].y -= (2.0 / (V * eps0 )) * Qi * ( ky/k2) * ek * ( - SinI * Qj * CosJ + CosI * Qj * SinJ);
cSim.pForce4[atomID1].z -= (2.0 / (V * eps0 )) * Qi * ( kz/k2) * ek * ( - SinI * Qj * CosJ + CosI * Qj * SinJ);
lowrz = 1 - numRz;
}
lowry = 1 - numRy;
}
}
atomID2++;
}
atomID1 += blockDim.x * gridDim.x;
}
}
platforms/cuda/src/kernels/kCalculateCDLJForces.cu
View file @ fe1e6ffa
...@@ -104,8 +104,7 @@ void GetCalculateCDLJForcesSim(gpuContext gpu) ...@@ -104,8 +104,7 @@ void GetCalculateCDLJForcesSim(gpuContext gpu)
// Include version of the kernel with Ewald method // Include version of the kernel with Ewald method
// Real Space Ewald uses almost the same kernels as Periodic // Real Space Ewald summation utilizes almost the same kernel as Periodic
#undef METHOD_NAME #undef METHOD_NAME
#undef USE_OUTPUT_BUFFER_PER_WARP #undef USE_OUTPUT_BUFFER_PER_WARP
#define USE_PERIODIC #define USE_PERIODIC
...@@ -118,8 +117,9 @@ void GetCalculateCDLJForcesSim(gpuContext gpu) ...@@ -118,8 +117,9 @@ void GetCalculateCDLJForcesSim(gpuContext gpu)
#define METHOD_NAME(a, b) a##EwaldDirectByWarp##b #define METHOD_NAME(a, b) a##EwaldDirectByWarp##b
#include "kCalculateCDLJForces.h" #include "kCalculateCDLJForces.h"
// Reciprocal Space Ewald summation is in a separate kernel // Reciprocal Space Ewald summation is in a separate kernel
//#include "kCalculateEwaldReciprocal.h" #include "kCalculateCDLJEwaldReciprocal.h"
__global__ extern void kCalculateCDLJCutoffForces_12_kernel(); __global__ extern void kCalculateCDLJCutoffForces_12_kernel();
...@@ -197,12 +197,20 @@ void kCalculateCDLJForces(gpuContext gpu) ...@@ -197,12 +197,20 @@ void kCalculateCDLJForces(gpuContext gpu)
kFindInteractionsWithinBlocksEwaldDirect_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block, kFindInteractionsWithinBlocksEwaldDirect_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
sizeof(unsigned int)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit); sizeof(unsigned int)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
if (gpu->bOutputBufferPerWarp) if (gpu->bOutputBufferPerWarp)
{
kCalculateCDLJEwaldDirectByWarpForces_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block, kCalculateCDLJEwaldDirectByWarpForces_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
(sizeof(Atom)+sizeof(float3))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit); (sizeof(Atom)+sizeof(float3))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
kCalculateCDLJEwaldReciprocalForces_kernel<<<gpu->sim.blocks, gpu->sim.update_threads_per_block>>>();
LAUNCHERROR("kCalculateCDLJEwaldReciprocalForces");
}
else else
{
kCalculateCDLJEwaldDirectForces_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block, kCalculateCDLJEwaldDirectForces_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
(sizeof(Atom)+sizeof(float3))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit); (sizeof(Atom)+sizeof(float3))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
LAUNCHERROR("kCalculateCDLJEwaldDirectForces"); LAUNCHERROR("kCalculateCDLJEwaldDirectForces");
kCalculateCDLJEwaldReciprocalForces_kernel<<<gpu->sim.blocks, gpu->sim.update_threads_per_block>>>();
LAUNCHERROR("kCalculateCDLJEwaldReciprocalForces");
}
} }
} }
platforms/cuda/src/kernels/kCalculateCDLJForces.h
View file @ fe1e6ffa
...@@ -43,9 +43,9 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni ...@@ -43,9 +43,9 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni
#endif #endif
#ifdef USE_EWALD #ifdef USE_EWALD
float alphaEwald = 3.123413; float alphaEwald = 3.123413;
float PI = 3.14159265358979323846f; float PI = 3.14159265358979323846f;
float SQRT_PI = sqrt(PI); float SQRT_PI = sqrt(PI);
#endif #endif
unsigned int lasty = 0xFFFFFFFF; unsigned int lasty = 0xFFFFFFFF;
...@@ -113,7 +113,7 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni ...@@ -113,7 +113,7 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni
#ifdef USE_EWALD #ifdef USE_EWALD
float r = sqrt(r2); float r = sqrt(r2);
float alphaR = alphaEwald * r; float alphaR = alphaEwald * r;
dEdR += apos.w * psA[tj].q * invR * invR * invR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ); dEdR += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
#else #else
dEdR += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2); dEdR += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2);
#endif #endif
...@@ -162,7 +162,7 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni ...@@ -162,7 +162,7 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni
#ifdef USE_EWALD #ifdef USE_EWALD
float r = sqrt(r2); float r = sqrt(r2);
float alphaR = alphaEwald * r; float alphaR = alphaEwald * r;
dEdR += apos.w * psA[tj].q * invR * invR * invR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ); dEdR += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
#else #else
dEdR += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2); dEdR += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2);
#endif #endif
...@@ -260,7 +260,7 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni ...@@ -260,7 +260,7 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni
#ifdef USE_EWALD #ifdef USE_EWALD
float r = sqrt(r2); float r = sqrt(r2);
float alphaR = alphaEwald * r; float alphaR = alphaEwald * r;
dEdR += apos.w * psA[tj].q * invR * invR * invR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ); dEdR += apos.w * psA[tj].q * invR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
#else #else
dEdR += apos.w * psA[tj].q * (invR - 2.0f * cSim.reactionFieldK * r2); dEdR += apos.w * psA[tj].q * (invR - 2.0f * cSim.reactionFieldK * r2);
#endif #endif
...@@ -315,7 +315,7 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni ...@@ -315,7 +315,7 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni
#ifdef USE_EWALD #ifdef USE_EWALD
float r = sqrt(r2); float r = sqrt(r2);
float alphaR = alphaEwald * r; float alphaR = alphaEwald * r;
dEdR += apos.w * psA[tj].q * invR * invR * invR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ); dEdR += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
#else #else
dEdR += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2); dEdR += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2);
#endif #endif
...@@ -406,7 +406,7 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni ...@@ -406,7 +406,7 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni
#ifdef USE_EWALD #ifdef USE_EWALD
float r = sqrt(r2); float r = sqrt(r2);
float alphaR = alphaEwald * r; float alphaR = alphaEwald * r;
dEdR += apos.w * psA[tj].q * invR * invR * invR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ); dEdR += apos.w * psA[tj].q * invR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
#else #else
dEdR += apos.w * psA[tj].q * (invR - 2.0f * cSim.reactionFieldK * r2); dEdR += apos.w * psA[tj].q * (invR - 2.0f * cSim.reactionFieldK * r2);
#endif #endif
......
platforms/cuda/src/kernels/kCalculateEwaldReciprocal.h deleted 100755 → 0
View file @ 57606845
/* -------------------------------------------------------------------------- *
* 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) 2009 Stanford University and the Authors. *
* Authors: Scott Le Grand, 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/>. *
* -------------------------------------------------------------------------- */
/**
* This file contains the kernel for evalauating nonbonded forces
* using the Ewald summation method.
*/
#include <cuComplex.h>
/* Define complex multiply operations */
__device__ cuComplex ComplexMul(cuComplex a, cuComplex b)
{
cuComplex c;
c.x = a.x * b.x - a.y * b.y;
c.y = a.x * b.y + a.y * b.x;
return c;
}
__device__ cuComplex ComplexConjMul(cuComplex a, cuComplex b)
{
cuComplex c;
c.x = a.x*b.x + a.y*b.y;
c.y = a.y*b.x - a.x*b.y;
return c;
}
__device__ cuComplex FloatComplexMul(float r, cuComplex a)
{
cuComplex b;
b.x = r*a.x;
b.y = r*a.y;
return b;
}
/* This kernel is under development */
__global__ void kCalculateCDLJEwaldForces_kernel(unsigned int* workUnit, int numWorkUnits)
{
// *******************************************************************
float alphaEwald = 3.123413f;
float factorEwald = -1.0 / (4*alphaEwald*alphaEwald);
float PI = 3.14159265358979323846f;
float SQRT_PI = sqrt(PI);
float TWO_PI = 2.0f * PI;
float epsilon = 1.0f;
/////##############################################################################
float recipCoeff = 4.0*PI/(cSim.periodicBoxSizeX * cSim.periodicBoxSizeY * cSim.periodicBoxSizeZ) /epsilon;
// setup reciprocal box
float recipBoxSizeX = TWO_PI / cSim.periodicBoxSizeX;
float recipBoxSizeY = TWO_PI / cSim.periodicBoxSizeY;
float recipBoxSizeZ = TWO_PI / cSim.periodicBoxSizeZ;
// setup K-vectors
unsigned int numRx = 60+1;
unsigned int numRy = 60+1;
unsigned int numRz = 60+1;
const int kmax = 61;
unsigned int pos = threadIdx.x + blockIdx.x * blockDim.x;
cuComplex eir[kmax][cSim.atoms][3];
cuComplex tab_xy[cSim.atoms];
cuComplex tab_qxyz[cSim.atoms];
while (pos < cSim.atoms)
{
float4 apos = cSim.pPosq[pos];
for(unsigned int m = 0; (m < 3); m++) {
eir[0][pos][m].x = 1;
eir[0][pos][m].y = 0;
}
eir[1][pos][0].x = cos(apos.x*recipBoxSizeX);
eir[1][pos][0].y = sin(apos.x*recipBoxSizeX);
eir[1][pos][1].x = cos(apos.y*recipBoxSizeY);
eir[1][pos][1].y = sin(apos.y*recipBoxSizeY);
eir[1][pos][2].x = cos(apos.z*recipBoxSizeZ);
eir[1][pos][2].y = sin(apos.z*recipBoxSizeZ);
for(unsigned int j=2; (j<kmax); j++) {
for(unsigned int m=0; (m<3); m++) {
eir[j][pos][m] = ComplexMul (eir[j-1][pos][m] , eir[1][pos][m]);
}
}
pos += blockDim.x * gridDim.x;
}
// *******************************************************************
int lowry = 0;
int lowrz = 1;
for(int rx = 0; rx < numRx; rx++) {
float kx = rx * recipBoxSizeX;
for(int ry = lowry; ry < numRy; ry++) {
float ky = ry * recipBoxSizeY;
if(ry >= 0) {
while (pos < cSim.atoms)
{
tab_xy[pos] = ComplexMul (eir[rx][pos][0] , eir[ry][pos][1]);
pos += blockDim.x * gridDim.x;
}
}
else {
while (pos < cSim.atoms)
{
tab_xy[pos]= ComplexConjMul (eir[rx][pos][0] , (eir[-ry][pos][1]));
pos += blockDim.x * gridDim.x;
}
}
for (int rz = lowrz; rz < numRz; rz++) {
float kz = rz * recipBoxSizeZ;
float k2 = kx * kx + ky * ky + kz * kz;
float ak = exp(k2*factorEwald) / k2;
float akv = 2.0 * ak * (1.0/k2 - factorEwald);
if( rz >= 0) {
while (pos < cSim.atoms)
{
float4 apos = cSim.pPosq[pos];
apos.w *= cSim.epsfac;
tab_qxyz[pos] = FloatComplexMul ( apos.w * ComplexMul (tab_xy[pos] , eir[rz][pos][2]));
pos += blockDim.x * gridDim.x;
}
}
else {
while (pos < cSim.atoms)
{
float4 apos = cSim.pPosq[pos];
apos.w *= cSim.epsfac;
tab_qxyz[pos] = FloatComplexMul( apos.w * ComplexConjMul (tab_xy[pos] , conj(eir[-rz][pos][2])) );
pos += blockDim.x * gridDim.x;
}
}
float cs = 0.0f;
float ss = 0.0f;
while (pos < cSim.atoms)
{
cs += tab_qxyz[pos].x;
ss += tab_qxyz[pos].y;
pos += blockDim.x * gridDim.x;
}
recipEnergy += ak * ( cs * cs + ss * ss);
float vir = akv * ( cs * cs + ss * ss);
while (pos < cSim.atoms)
{
float4 force = cSim.pForce4[pos];
float dEdR = ak * (cs * tab_qxyz[pos].y - ss * tab_qxyz[pos].x);
force.x += 2.0 * recipCoeff * dEdR * kx ;
force.y += 2.0 * recipCoeff * dEdR * ky ;
force.z += 2.0 * recipCoeff * dEdR * kz ;
}
lowrz = 1 - numRz;
}
lowry = 1 - numRy;
}
}
//###########################################################################
//###########################################################################
// END EWALD RECIP SPACE
//###########################################################################
}
platforms/cuda/tests/TestCudaEwald.cpp deleted 100755 → 0
View file @ 57606845
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* 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 reference implementation of NonbondedForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "openmm/OpenMMContext.h"
#include "CudaPlatform.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/OpenMMContextImpl.h"
#include "kernels/gputypes.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "../src/sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testEwald() {
CudaPlatform platform;
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->setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
system.addForce(nonbonded);
OpenMMContext 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();
cout << " force 0: " << forces[0] << endl;
cout << " force 1: " << forces[1] << endl;
cout << " energyPoten: " << state.getPotentialEnergy() << endl;
ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(123.711, -64.1877, 302.716), forces[1], TOL);
//ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
}
void testPeriodic() {
CudaPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
HarmonicBondForce* bonds = new HarmonicBondForce();
bonds->addBond(0, 1, 1, 0);
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
nonbonded->setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
system.addForce(nonbonded);
OpenMMContext 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();
cout << "" << endl;
cout << " force 0: " << forces[0] << endl;
cout << " force 1: " << forces[1] << endl;
cout << " force 2: " << forces[2] << endl;
cout << "energyPoten: " << state.getPotentialEnergy() << endl;
cout << "" << endl;
cout << " no cutoff force: 138.935485" << endl;
cout << "" << endl;
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 = 138.935485*(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*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
}
int main() {
try {
cout << "" << endl;
cout << "Executing Periodic..." << endl;
testPeriodic();
cout << "" << endl;
cout << "Executing Ewald..." << endl;
testEwald();
cout << "" << endl;
cout << "Done" << endl;
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cuda/tests/TestCudaNonbondedForce.cpp
View file @ fe1e6ffa
...@@ -357,6 +357,34 @@ void testPeriodic() { ...@@ -357,6 +357,34 @@ void testPeriodic() {
ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
} }
void testEwald() {
CudaPlatform platform;
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->setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
system.addForce(nonbonded);
OpenMMContext 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], TOL);
ASSERT_EQUAL_VEC(Vec3( 123.711, -64.1877, 302.716), forces[1], TOL);
ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), TOL);
}
void testLargeSystem() { void testLargeSystem() {
const int numMolecules = 600; const int numMolecules = 600;
const int numParticles = numMolecules*2; const int numParticles = numMolecules*2;
...@@ -602,6 +630,7 @@ int main() { ...@@ -602,6 +630,7 @@ int main() {
testCutoff(); testCutoff();
testCutoff14(); testCutoff14();
testPeriodic(); testPeriodic();
testEwald();
testLargeSystem(); testLargeSystem();
testBlockInteractions(false); testBlockInteractions(false);
testBlockInteractions(true); testBlockInteractions(true);
......
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
View file @ fe1e6ffa
...@@ -234,14 +234,13 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at ...@@ -234,14 +234,13 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
#include "../SimTKUtilities/RealTypeSimTk.h" #include "../SimTKUtilities/RealTypeSimTk.h"
typedef std::complex<RealOpenMM> d_complex; typedef std::complex<RealOpenMM> d_complex;
typedef std::complex<int> int_complex;
// Number of R-vectors (real space vectors) // Number of R-vectors (real space vectors)
// to be calculated automatically eventually from alphaEwald and desired precision // to be calculated automatically eventually from alphaEwald and desired precision
int numRx = 60+1; int numRx = 20+1;
int numRy = 60+1; int numRy = 20+1;
int numRz = 60+1; int numRz = 20+1;
int kmax = std::max(numRx, std::max(numRy,numRz)); int kmax = std::max(numRx, std::max(numRy,numRz));
static const RealOpenMM alphaEwald = (RealOpenMM) 3.123413; static const RealOpenMM alphaEwald = (RealOpenMM) 3.123413;
...@@ -253,7 +252,7 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at ...@@ -253,7 +252,7 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
static const RealOpenMM epsilon = 1.0; static const RealOpenMM epsilon = 1.0;
static const RealOpenMM one = 1.0; static const RealOpenMM one = 1.0;
RealOpenMM recipCoeff = (RealOpenMM)(4*M_PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon); RealOpenMM recipCoeff = (RealOpenMM)(4*PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
RealOpenMM selfEwaldEnergy = 0.0; RealOpenMM selfEwaldEnergy = 0.0;
RealOpenMM realSpaceEwaldEnergy = 0.0; RealOpenMM realSpaceEwaldEnergy = 0.0;
...@@ -283,9 +282,6 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at ...@@ -283,9 +282,6 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
d_complex* eir = new d_complex[kmax*numberOfAtoms*3]; d_complex* eir = new d_complex[kmax*numberOfAtoms*3];
d_complex* tab_xy = new d_complex[numberOfAtoms]; d_complex* tab_xy = new d_complex[numberOfAtoms];
d_complex* tab_qxyz = new d_complex[numberOfAtoms]; d_complex* tab_qxyz = new d_complex[numberOfAtoms];
d_complex a,b,c;
int i,j,m;
if (kmax < 1) { if (kmax < 1) {
std::stringstream message; std::stringstream message;
...@@ -293,16 +289,16 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at ...@@ -293,16 +289,16 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
SimTKOpenMMLog::printError( message ); SimTKOpenMMLog::printError( message );
} }
for(i = 0; (i < numberOfAtoms); i++) { for(int i = 0; (i < numberOfAtoms); i++) {
for(m = 0; (m < 3); m++) for(int m = 0; (m < 3); m++)
EIR(0, i, m) = d_complex(1,0); EIR(0, i, m) = d_complex(1,0);
for(m=0; (m<3); m++) for(int m=0; (m<3); m++)
EIR(1, i, m) = d_complex(cos(atomCoordinates[i][m]*recipBoxSize[m]), EIR(1, i, m) = d_complex(cos(atomCoordinates[i][m]*recipBoxSize[m]),
sin(atomCoordinates[i][m]*recipBoxSize[m])); sin(atomCoordinates[i][m]*recipBoxSize[m]));
for(j=2; (j<kmax); j++) for(int j=2; (j<kmax); j++)
for(m=0; (m<3); m++) for(int m=0; (m<3); m++)
EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m); EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
} }
...@@ -353,10 +349,8 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at ...@@ -353,10 +349,8 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
RealOpenMM kz = rz * recipBoxSize[2]; RealOpenMM kz = rz * recipBoxSize[2];
RealOpenMM k2 = kx * kx + ky * ky + kz * kz; RealOpenMM k2 = kx * kx + ky * ky + kz * kz;
RealOpenMM ak = exp(k2*factorEwald) / k2; RealOpenMM ak = exp(k2*factorEwald) / k2;
RealOpenMM akv = 2 * ak * (1/k2 - factorEwald);
recipEnergy += ak * ( cs * cs + ss * ss); recipEnergy += ak * ( cs * cs + ss * ss);
RealOpenMM vir = akv * ( cs * cs + ss * ss);
for(int n = 0; n < numberOfAtoms; n++) { for(int n = 0; n < numberOfAtoms; n++) {
RealOpenMM force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real()); RealOpenMM force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real());
......
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp.Ewald-optimized 0 → 100644
View file @ fe1e6ffa
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* 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 <string.h>
#include <sstream>
#include <complex>
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceLJCoulombIxn.h"
#include "ReferenceForce.h"
// In case we're using some primitive version of Visual Studio this will
// make sure that erf() and erfc() are defined.
#include "MSVC_erfc.h"
/**---------------------------------------------------------------------------------------
ReferenceLJCoulombIxn constructor
--------------------------------------------------------------------------------------- */
ReferenceLJCoulombIxn::ReferenceLJCoulombIxn( ) : cutoff(false), periodic(false) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceLJCoulombIxn::ReferenceLJCoulombIxn";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
ReferenceLJCoulombIxn destructor
--------------------------------------------------------------------------------------- */
ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceLJCoulombIxn::~ReferenceLJCoulombIxn";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
@param distance the cutoff distance
@param neighbors the neighbor list to use
@param solventDielectric the dielectric constant of the bulk solvent
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric ) {
cutoff = true;
cutoffDistance = distance;
neighborList = &neighbors;
krf = pow(cutoffDistance, -3.0f)*(solventDielectric-1.0f)/(2.0f*solventDielectric+1.0f);
crf = (1.0f/cutoffDistance)*(3.0f*solventDielectric)/(2.0f*solventDielectric+1.0f);
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions. This requires that a cutoff has
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::setPeriodic( RealOpenMM* boxSize ) {
assert(cutoff);
assert(boxSize[0] >= 2.0*cutoffDistance);
assert(boxSize[1] >= 2.0*cutoffDistance);
assert(boxSize[2] >= 2.0*cutoffDistance);
periodic = true;
periodicBoxSize[0] = boxSize[0];
periodicBoxSize[1] = boxSize[1];
periodicBoxSize[2] = boxSize[2];
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate parameters for LJ Coulomb ixn
@param c6 c6
@param c12 c12
@param q1 q1 charge atom 1
@param epsfac epsfacSqrt ????????????
@param parameters output parameters:
parameter[SigIndex] = 0.5*( (c12/c6)**1/6 ) (sigma/2)
parameter[EpsIndex] = sqrt(c6*c6/c12) (2*sqrt(epsilon))
parameter[QIndex] = epsfactorSqrt*q1
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM epsfacSqrt,
RealOpenMM* parameters ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceLJCoulombIxn::getDerivedParameters";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM six = 6.0;
static const RealOpenMM half = 0.5;
static const RealOpenMM oneSixth = one/six;
static const RealOpenMM oneTweleth = half*oneSixth;
// ---------------------------------------------------------------------------------------
if( c12 <= 0.0 ){
parameters[EpsIndex] = zero;
parameters[SigIndex] = half;
} else {
parameters[EpsIndex] = c6*SQRT( one/c12 );
parameters[SigIndex] = POW( (c12/c6), oneSixth );
parameters[SigIndex] *= half;
}
parameters[QIndex] = epsfacSqrt*q1;
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Set the reciprocal space vectors to use with Ewald
// Currently a dumb routine, vectors are set in calculateEwaldIxn
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::setRecipVectors() {
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate the Ewald parameter based on the cutoff and the desired tolerance using
erfc( alpha*cutoff )/cutoff < tolerance
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
/* int ReferenceLJCoulombIxn::setAlphaEwald(RealOpenMM cutoff, RealOpenMM tolerance,
RealOpenMM alphaEwald) {
alphaEwald = 1.0f;
while ( erfc(alphaEwald*cutoff) >= tolerance*cutoff) {
alphaEwald= 1.2 * alphaEwald;
}
return ReferenceForce::DefaultReturn;
}
*/
/**---------------------------------------------------------------------------------------
Calculate Ewald ixn
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, RealOpenMM** forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
#include "../SimTKUtilities/RealTypeSimTk.h"
typedef std::complex<RealOpenMM> d_complex;
// Number of R-vectors (real space vectors)
// to be calculated automatically eventually from alphaEwald and desired precision
int numRx = 20+1;
int numRy = 20+1;
int numRz = 20+1;
int kmax = std::max(numRx, std::max(numRy,numRz));
static const RealOpenMM alphaEwald = (RealOpenMM) 3.123413;
RealOpenMM factorEwald = -1 / (4*alphaEwald*alphaEwald);
RealOpenMM SQRT_PI = sqrt(PI);
RealOpenMM TWO_PI = 2.0 * PI;
static const RealOpenMM epsilon = 1.0;
static const RealOpenMM one = 1.0;
RealOpenMM recipCoeff = (RealOpenMM)(4*PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
RealOpenMM selfEwaldEnergy = 0.0;
RealOpenMM realSpaceEwaldEnergy = 0.0;
RealOpenMM recipEnergy = 0.0;
// **************************************************************************************
// SELF ENERGY
// **************************************************************************************
for( int atomID = 0; atomID < numberOfAtoms; atomID++ ){
selfEwaldEnergy = selfEwaldEnergy + atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex];
}
selfEwaldEnergy = selfEwaldEnergy * alphaEwald/SQRT_PI ;
// **************************************************************************************
// RECIPROCAL SPACE EWALD ENERGY AND FORCES
// **************************************************************************************
// setup reciprocal box
RealOpenMM recipBoxSize[3] = { TWO_PI / periodicBoxSize[0], TWO_PI / periodicBoxSize[1], TWO_PI / periodicBoxSize[2]};
// setup K-vectors
#define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
d_complex* eir = new d_complex[kmax*numberOfAtoms*3];
d_complex* tab_xy = new d_complex[numberOfAtoms];
d_complex* tab_qxyz = new d_complex[numberOfAtoms];
if (kmax < 1) {
std::stringstream message;
message << " kmax < 1 , Aborting" << std::endl;
SimTKOpenMMLog::printError( message );
}
for(int i = 0; (i < numberOfAtoms); i++) {
for(int m = 0; (m < 3); m++)
EIR(0, i, m) = d_complex(1,0);
for(int m=0; (m<3); m++)
EIR(1, i, m) = d_complex(cos(atomCoordinates[i][m]*recipBoxSize[m]),
sin(atomCoordinates[i][m]*recipBoxSize[m]));
for(int j=2; (j<kmax); j++)
for(int m=0; (m<3); m++)
EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
}
// calculate reciprocal space energy and forces
int lowry = 0;
int lowrz = 1;
for(int rx = 0; rx < numRx; rx++) {
RealOpenMM kx = rx * recipBoxSize[0];
for(int ry = lowry; ry < numRy; ry++) {
RealOpenMM ky = ry * recipBoxSize[1];
if(ry >= 0) {
for(int n = 0; n < numberOfAtoms; n++)
tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
}
else {
for(int n = 0; n < numberOfAtoms; n++)
tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
}
for (int rz = lowrz; rz < numRz; rz++) {
if( rz >= 0) {
for( int n = 0; n < numberOfAtoms; n++)
tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2));
}
else {
for( int n = 0; n < numberOfAtoms; n++)
tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2)));
}
RealOpenMM cs = 0.0f;
RealOpenMM ss = 0.0f;
for( int n = 0; n < numberOfAtoms; n++) {
cs += tab_qxyz[n].real();
ss += tab_qxyz[n].imag();
}
RealOpenMM kz = rz * recipBoxSize[2];
RealOpenMM k2 = kx * kx + ky * ky + kz * kz;
RealOpenMM ak = exp(k2*factorEwald) / k2;
recipEnergy += ak * ( cs * cs + ss * ss);
for(int n = 0; n < numberOfAtoms; n++) {
RealOpenMM force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real());
forces[n][0] += 2 * recipCoeff * force * kx ;
forces[n][1] += 2 * recipCoeff * force * ky ;
forces[n][2] += 2 * recipCoeff * force * kz ;
}
lowrz = 1 - numRz;
}
lowry = 1 - numRy;
}
}
recipEnergy *= recipCoeff;
// **************************************************************************************
// SHORT-RANGE ENERGY AND FORCES
// **************************************************************************************
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
for( int atomID1 = 0; atomID1 < numberOfAtoms; atomID1++ ){
for( int atomID2 = atomID1 + 1; atomID2 < numberOfAtoms; atomID2++ ){
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomID2], atomCoordinates[atomID1], periodicBoxSize, deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
realSpaceEwaldEnergy =
(RealOpenMM)(realSpaceEwaldEnergy + atomParameters[atomID1][QIndex]*atomParameters[atomID2][QIndex]*inverseR*erfc(alphaEwald*r));
}
}
// allocate and initialize exclusion array
int* exclusionIndices = new int[numberOfAtoms];
for( int ii = 0; ii < numberOfAtoms; ii++ ){
exclusionIndices[ii] = -1;
}
for( int ii = 0; ii < numberOfAtoms; ii++ ){
// set exclusions
for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
exclusionIndices[exclusions[ii][jj]] = ii;
}
// loop over atom pairs
for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
if( exclusionIndices[jj] != ii ){
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r;
realSpaceEwaldEnergy =
(RealOpenMM)(realSpaceEwaldEnergy + atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR));
RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
dEdR = (RealOpenMM)(dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
for( int kk = 0; kk < 3; kk++ ){
RealOpenMM force = dEdR*deltaR[0][kk];
forces[ii][kk] += force;
forces[jj][kk] -= force;
}
}
}
}
delete[] exclusionIndices;
delete[] eir;
delete[] tab_xy;
delete[] tab_qxyz;
// ***********************************************************************
if( totalEnergy ) {
*totalEnergy += recipEnergy + realSpaceEwaldEnergy - selfEwaldEnergy;
}
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate LJ Coulomb pair ixn
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, RealOpenMM** forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
if (cutoff) {
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
}
}
else {
// allocate and initialize exclusion array
int* exclusionIndices = new int[numberOfAtoms];
for( int ii = 0; ii < numberOfAtoms; ii++ ){
exclusionIndices[ii] = -1;
}
for( int ii = 0; ii < numberOfAtoms; ii++ ){
// set exclusions
for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
exclusionIndices[exclusions[ii][jj]] = ii;
}
// loop over atom pairs
for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
if( exclusionIndices[jj] != ii ){
calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
}
}
}
delete[] exclusionIndices;
}
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate LJ Coulomb pair ixn between two atoms
@param ii the index of the first atom
@param jj the index of the second atom
@param atomCoordinates atom coordinates
@param atomParameters atom parameters (charges, c6, c12, ...) atomParameters[atomIndex][paramterIndex]
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoordinates,
RealOpenMM** atomParameters, RealOpenMM** forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\nReferenceLJCoulombIxn::calculateOneIxn";
// ---------------------------------------------------------------------------------------
// constants -- reduce Visual Studio warnings regarding conversions between float & double
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM six = 6.0;
static const RealOpenMM twelve = 12.0;
static const RealOpenMM oneM = -1.0;
static const int threeI = 3;
// debug flag
static const int debug = -1;
static const int LastAtomIndex = 2;
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
// get deltaR, R2, and R between 2 atoms
if (periodic)
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
else
ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] );
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM sig = atomParameters[ii][SigIndex] + atomParameters[jj][SigIndex];
RealOpenMM sig2 = inverseR*sig;
sig2 *= sig2;
RealOpenMM sig6 = sig2*sig2*sig2;
RealOpenMM eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
RealOpenMM dEdR = eps*( twelve*sig6 - six )*sig6;
if (cutoff)
dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
else
dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
dEdR *= inverseR*inverseR;
// accumulate forces
for( int kk = 0; kk < 3; kk++ ){
RealOpenMM force = dEdR*deltaR[0][kk];
forces[ii][kk] += force;
forces[jj][kk] -= force;
}
RealOpenMM energy = 0.0;
// accumulate energies
if( totalEnergy || energyByAtom ) {
if (cutoff)
energy = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf);
else
energy = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
energy += eps*(sig6-one)*sig6;
if( totalEnergy )
*totalEnergy += energy;
if( energyByAtom ){
energyByAtom[ii] += energy;
energyByAtom[jj] += energy;
}
}
// debug
if( debug == ii ){
static bool printHeader = false;
std::stringstream message;
message << methodName;
message << std::endl;
int pairArray[2] = { ii, jj };
if( !printHeader ){
printHeader = true;
message << std::endl;
message << methodName.c_str() << " a0 k [c q p s] r1 r2 angle dt rp p[] dot cosine angle dEdR*r F[]" << std::endl;
}
message << std::endl;
for( int kk = 0; kk < 2; kk++ ){
message << " Atm " << pairArray[kk] << " [" << atomCoordinates[pairArray[kk]][0] << " " << atomCoordinates[pairArray[kk]][1] << " " << atomCoordinates[pairArray[kk]][2] << "] ";
}
message << std::endl << " Delta:";
for( int kk = 0; kk < (LastAtomIndex - 1); kk++ ){
message << " [";
for( int jj = 0; jj < ReferenceForce::LastDeltaRIndex; jj++ ){
message << deltaR[kk][jj] << " ";
}
message << "]";
}
message << std::endl;
for( int kk = 0; kk < 2; kk++ ){
message << " p" << pairArray[kk] << " [";
message << atomParameters[pairArray[kk]][0] << " " << atomParameters[pairArray[kk]][1] << " " << atomParameters[pairArray[kk]][2];
message << "]";
}
message << std::endl;
message << " dEdR=" << dEdR;
message << " E=" << energy << " force factors: ";
message << "F=compute force; f=cumulative force";
message << std::endl << " ";
message << " f" << ii << "[";
SimTKOpenMMUtilities::formatRealStringStream( message, deltaR[0], threeI, dEdR );
message << "]";
for( int kk = 0; kk < 2; kk++ ){
message << " F" << pairArray[kk] << " [";
SimTKOpenMMUtilities::formatRealStringStream( message, forces[pairArray[kk]], threeI );
message << "]";
}
SimTKOpenMMLog::printMessage( message );
}
return ReferenceForce::DefaultReturn;
}
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp.Ewald-vanilla 0 → 100644
View file @ fe1e6ffa
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* 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 <string.h>
#include <sstream>
#include <complex>
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceLJCoulombIxn.h"
#include "ReferenceForce.h"
// In case we're using some primitive version of Visual Studio this will
// make sure that erf() and erfc() are defined.
#include "MSVC_erfc.h"
/**---------------------------------------------------------------------------------------
ReferenceLJCoulombIxn constructor
--------------------------------------------------------------------------------------- */
ReferenceLJCoulombIxn::ReferenceLJCoulombIxn( ) : cutoff(false), periodic(false) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceLJCoulombIxn::ReferenceLJCoulombIxn";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
ReferenceLJCoulombIxn destructor
--------------------------------------------------------------------------------------- */
ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceLJCoulombIxn::~ReferenceLJCoulombIxn";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
@param distance the cutoff distance
@param neighbors the neighbor list to use
@param solventDielectric the dielectric constant of the bulk solvent
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric ) {
cutoff = true;
cutoffDistance = distance;
neighborList = &neighbors;
krf = pow(cutoffDistance, -3.0f)*(solventDielectric-1.0f)/(2.0f*solventDielectric+1.0f);
crf = (1.0f/cutoffDistance)*(3.0f*solventDielectric)/(2.0f*solventDielectric+1.0f);
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions. This requires that a cutoff has
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::setPeriodic( RealOpenMM* boxSize ) {
assert(cutoff);
assert(boxSize[0] >= 2.0*cutoffDistance);
assert(boxSize[1] >= 2.0*cutoffDistance);
assert(boxSize[2] >= 2.0*cutoffDistance);
periodic = true;
periodicBoxSize[0] = boxSize[0];
periodicBoxSize[1] = boxSize[1];
periodicBoxSize[2] = boxSize[2];
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate parameters for LJ Coulomb ixn
@param c6 c6
@param c12 c12
@param q1 q1 charge atom 1
@param epsfac epsfacSqrt ????????????
@param parameters output parameters:
parameter[SigIndex] = 0.5*( (c12/c6)**1/6 ) (sigma/2)
parameter[EpsIndex] = sqrt(c6*c6/c12) (2*sqrt(epsilon))
parameter[QIndex] = epsfactorSqrt*q1
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM epsfacSqrt,
RealOpenMM* parameters ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceLJCoulombIxn::getDerivedParameters";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM six = 6.0;
static const RealOpenMM half = 0.5;
static const RealOpenMM oneSixth = one/six;
static const RealOpenMM oneTweleth = half*oneSixth;
// ---------------------------------------------------------------------------------------
if( c12 <= 0.0 ){
parameters[EpsIndex] = zero;
parameters[SigIndex] = half;
} else {
parameters[EpsIndex] = c6*SQRT( one/c12 );
parameters[SigIndex] = POW( (c12/c6), oneSixth );
parameters[SigIndex] *= half;
}
parameters[QIndex] = epsfacSqrt*q1;
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Set the reciprocal space vectors to use with Ewald
// Currently a dumb routine, vectors are set in calculateEwaldIxn
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::setRecipVectors() {
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate the Ewald parameter based on the cutoff and the desired tolerance using
erfc( alpha*cutoff )/cutoff < tolerance
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
/* int ReferenceLJCoulombIxn::setAlphaEwald(RealOpenMM cutoff, RealOpenMM tolerance,
RealOpenMM alphaEwald) {
alphaEwald = 1.0f;
while ( erfc(alphaEwald*cutoff) >= tolerance*cutoff) {
alphaEwald= 1.2 * alphaEwald;
}
return ReferenceForce::DefaultReturn;
}
*/
/**---------------------------------------------------------------------------------------
Calculate Ewald ixn
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, RealOpenMM** forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
#include "../SimTKUtilities/RealTypeSimTk.h"
typedef std::complex<RealOpenMM> d_complex;
// Number of R-vectors (real space vectors)
// to be calculated automatically eventually from alphaEwald and desired precision
int numRx = 20+1;
int numRy = 20+1;
int numRz = 20+1;
int kmax = std::max(numRx, std::max(numRy,numRz));
static const RealOpenMM alphaEwald = (RealOpenMM) 3.123413;
RealOpenMM factorEwald = -1 / (4*alphaEwald*alphaEwald);
RealOpenMM SQRT_PI = sqrt(PI);
RealOpenMM TWO_PI = 2.0 * PI;
static const RealOpenMM epsilon = 1.0;
static const RealOpenMM one = 1.0;
RealOpenMM recipCoeff = (RealOpenMM)(4*PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
RealOpenMM selfEwaldEnergy = 0.0;
RealOpenMM realSpaceEwaldEnergy = 0.0;
RealOpenMM recipEnergy = 0.0;
// **************************************************************************************
// SELF ENERGY
// **************************************************************************************
for( int atomID = 0; atomID < numberOfAtoms; atomID++ ){
selfEwaldEnergy = selfEwaldEnergy + atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex];
}
selfEwaldEnergy = selfEwaldEnergy * alphaEwald/SQRT_PI ;
// **************************************************************************************
// RECIPROCAL SPACE EWALD ENERGY AND FORCES
// **************************************************************************************
// setup reciprocal box
RealOpenMM recipBoxSize[3] = { TWO_PI / periodicBoxSize[0], TWO_PI / periodicBoxSize[1], TWO_PI / periodicBoxSize[2]};
// calculate reciprocal space energy and forces
RealOpenMM pi4eps = 11.787089;
RealOpenMM eps0 = 5.72765E-4;
RealOpenMM V = periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2];
for( int atomID1 = 0; atomID1 < numberOfAtoms; atomID1++ ){
for( int atomID2 = 0; atomID2 < numberOfAtoms; atomID2++ ){
int lowry = 0;
int lowrz = 1;
for(int rx = 0; rx < numRx; rx++) {
RealOpenMM kx = rx * recipBoxSize[0];
for(int ry = lowry; ry < numRy; ry++) {
RealOpenMM ky = ry * recipBoxSize[1];
for (int rz = lowrz; rz < numRz; rz++) {
RealOpenMM kz = rz * recipBoxSize[2];
RealOpenMM k2 = kx * kx + ky * ky + kz * kz;
RealOpenMM ek = exp ( k2 * factorEwald);
RealOpenMM Qi = atomParameters[atomID1][QIndex] / pi4eps;
RealOpenMM Qj = atomParameters[atomID2][QIndex] / pi4eps;
RealOpenMM Xi = atomCoordinates[atomID1][0];
RealOpenMM Yi = atomCoordinates[atomID1][1];
RealOpenMM Zi = atomCoordinates[atomID1][2];
RealOpenMM Xj = atomCoordinates[atomID2][0];
RealOpenMM Yj = atomCoordinates[atomID2][1];
RealOpenMM Zj = atomCoordinates[atomID2][2];
RealOpenMM SinI = sin ( kx * Xi + ky * Yi + kz * Zi );
RealOpenMM SinJ = sin ( kx * Xj + ky * Yj + kz * Zj );
RealOpenMM CosI = cos ( kx * Xi + ky * Yi + kz * Zi );
RealOpenMM CosJ = cos ( kx * Xj + ky * Yj + kz * Zj );
RealOpenMM ax = (2.0 / (V * eps0 )) * Qi * ( kx/k2) * ek * ( - SinI * Qj * CosJ + CosI * Qj * SinJ);
RealOpenMM ay = (2.0 / (V * eps0 )) * Qi * ( ky/k2) * ek * ( - SinI * Qj * CosJ + CosI * Qj * SinJ);
RealOpenMM az = (2.0 / (V * eps0 )) * Qi * ( kz/k2) * ek * ( - SinI * Qj * CosJ + CosI * Qj * SinJ);
forces[atomID1][0] -= ax;
forces[atomID1][1] -= ay;
forces[atomID1][2] -= az;
lowrz = 1 - numRz;
}
lowry = 1 - numRy;
}
}
}
}
recipEnergy *= recipCoeff;
// **************************************************************************************
// SHORT-RANGE ENERGY AND FORCES
// **************************************************************************************
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
for( int atomID1 = 0; atomID1 < numberOfAtoms; atomID1++ ){
for( int atomID2 = atomID1 + 1; atomID2 < numberOfAtoms; atomID2++ ){
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomID2], atomCoordinates[atomID1], periodicBoxSize, deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
realSpaceEwaldEnergy =
(RealOpenMM)(realSpaceEwaldEnergy + atomParameters[atomID1][QIndex]*atomParameters[atomID2][QIndex]*inverseR*erfc(alphaEwald*r));
}
}
// allocate and initialize exclusion array
int* exclusionIndices = new int[numberOfAtoms];
for( int ii = 0; ii < numberOfAtoms; ii++ ){
exclusionIndices[ii] = -1;
}
for( int ii = 0; ii < numberOfAtoms; ii++ ){
// set exclusions
for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
exclusionIndices[exclusions[ii][jj]] = ii;
}
// loop over atom pairs
for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
if( exclusionIndices[jj] != ii ){
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r;
realSpaceEwaldEnergy =
(RealOpenMM)(realSpaceEwaldEnergy + atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR));
RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
dEdR = (RealOpenMM)(dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
for( int kk = 0; kk < 3; kk++ ){
RealOpenMM force = dEdR*deltaR[0][kk];
forces[ii][kk] += force;
forces[jj][kk] -= force;
}
}
}
}
delete[] exclusionIndices;
// ***********************************************************************
if( totalEnergy ) {
*totalEnergy += recipEnergy + realSpaceEwaldEnergy - selfEwaldEnergy;
}
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate LJ Coulomb pair ixn
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, RealOpenMM** forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
if (cutoff) {
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
}
}
else {
// allocate and initialize exclusion array
int* exclusionIndices = new int[numberOfAtoms];
for( int ii = 0; ii < numberOfAtoms; ii++ ){
exclusionIndices[ii] = -1;
}
for( int ii = 0; ii < numberOfAtoms; ii++ ){
// set exclusions
for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
exclusionIndices[exclusions[ii][jj]] = ii;
}
// loop over atom pairs
for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
if( exclusionIndices[jj] != ii ){
calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
}
}
}
delete[] exclusionIndices;
}
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate LJ Coulomb pair ixn between two atoms
@param ii the index of the first atom
@param jj the index of the second atom
@param atomCoordinates atom coordinates
@param atomParameters atom parameters (charges, c6, c12, ...) atomParameters[atomIndex][paramterIndex]
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoordinates,
RealOpenMM** atomParameters, RealOpenMM** forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\nReferenceLJCoulombIxn::calculateOneIxn";
// ---------------------------------------------------------------------------------------
// constants -- reduce Visual Studio warnings regarding conversions between float & double
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM six = 6.0;
static const RealOpenMM twelve = 12.0;
static const RealOpenMM oneM = -1.0;
static const int threeI = 3;
// debug flag
static const int debug = -1;
static const int LastAtomIndex = 2;
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
// get deltaR, R2, and R between 2 atoms
if (periodic)
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
else
ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] );
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM sig = atomParameters[ii][SigIndex] + atomParameters[jj][SigIndex];
RealOpenMM sig2 = inverseR*sig;
sig2 *= sig2;
RealOpenMM sig6 = sig2*sig2*sig2;
RealOpenMM eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
RealOpenMM dEdR = eps*( twelve*sig6 - six )*sig6;
if (cutoff)
dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
else
dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
dEdR *= inverseR*inverseR;
// accumulate forces
for( int kk = 0; kk < 3; kk++ ){
RealOpenMM force = dEdR*deltaR[0][kk];
forces[ii][kk] += force;
forces[jj][kk] -= force;
}
RealOpenMM energy = 0.0;
// accumulate energies
if( totalEnergy || energyByAtom ) {
if (cutoff)
energy = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf);
else
energy = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
energy += eps*(sig6-one)*sig6;
if( totalEnergy )
*totalEnergy += energy;
if( energyByAtom ){
energyByAtom[ii] += energy;
energyByAtom[jj] += energy;
}
}
// debug
if( debug == ii ){
static bool printHeader = false;
std::stringstream message;
message << methodName;
message << std::endl;
int pairArray[2] = { ii, jj };
if( !printHeader ){
printHeader = true;
message << std::endl;
message << methodName.c_str() << " a0 k [c q p s] r1 r2 angle dt rp p[] dot cosine angle dEdR*r F[]" << std::endl;
}
message << std::endl;
for( int kk = 0; kk < 2; kk++ ){
message << " Atm " << pairArray[kk] << " [" << atomCoordinates[pairArray[kk]][0] << " " << atomCoordinates[pairArray[kk]][1] << " " << atomCoordinates[pairArray[kk]][2] << "] ";
}
message << std::endl << " Delta:";
for( int kk = 0; kk < (LastAtomIndex - 1); kk++ ){
message << " [";
for( int jj = 0; jj < ReferenceForce::LastDeltaRIndex; jj++ ){
message << deltaR[kk][jj] << " ";
}
message << "]";
}
message << std::endl;
for( int kk = 0; kk < 2; kk++ ){
message << " p" << pairArray[kk] << " [";
message << atomParameters[pairArray[kk]][0] << " " << atomParameters[pairArray[kk]][1] << " " << atomParameters[pairArray[kk]][2];
message << "]";
}
message << std::endl;
message << " dEdR=" << dEdR;
message << " E=" << energy << " force factors: ";
message << "F=compute force; f=cumulative force";
message << std::endl << " ";
message << " f" << ii << "[";
SimTKOpenMMUtilities::formatRealStringStream( message, deltaR[0], threeI, dEdR );
message << "]";
for( int kk = 0; kk < 2; kk++ ){
message << " F" << pairArray[kk] << " [";
SimTKOpenMMUtilities::formatRealStringStream( message, forces[pairArray[kk]], threeI );
message << "]";
}
SimTKOpenMMLog::printMessage( message );
}
return ReferenceForce::DefaultReturn;
}
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